~maddevelopers/mg5amcnlo/new_clustering

« back to all changes in this revision

Viewing changes to madgraph/fks/fks_base.py

  • Committer: Rikkert Frederix
  • Date: 2021-09-09 15:51:40 UTC
  • mfrom: (78.75.502 3.2.1)
  • Revision ID: frederix@physik.uzh.ch-20210909155140-rg6umfq68h6h47cf
merge with 3.2.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
15
15
 
16
16
"""Definitions of the objects needed for the implementation of MadFKS"""
17
17
 
 
18
from __future__ import absolute_import
 
19
from __future__ import print_function
18
20
import madgraph.core.base_objects as MG
19
21
import madgraph.core.helas_objects as helas_objects
20
22
import madgraph.core.diagram_generation as diagram_generation
27
29
import array
28
30
import madgraph.various.misc as misc
29
31
from madgraph import InvalidCmd
 
32
from six.moves import range
30
33
 
31
34
logger = logging.getLogger('madgraph.fks_base')
32
35
 
36
39
#===============================================================================
37
40
# FKS Process
38
41
#===============================================================================
 
42
 
 
43
 
39
44
class FKSMultiProcess(diagram_generation.MultiProcess): #test written
40
45
    """A multi process class that contains informations on the born processes 
41
46
    and the reals.
42
47
    """
43
 
    
 
48
 
44
49
    def default_setup(self):
45
50
        """Default values for all properties"""
46
51
        super(FKSMultiProcess, self).default_setup()
47
52
        self['real_amplitudes'] = diagram_generation.AmplitudeList()
48
53
        self['pdgs'] = []
49
54
        self['born_processes'] = FKSProcessList()
50
 
        if not 'OLP' in self.keys():
 
55
 
 
56
        if not 'OLP' in list(self.keys()):
51
57
            self['OLP'] = 'MadLoop'
52
58
            self['ncores_for_proc_gen'] = 0
53
59
    
55
61
        """Return particle property names as a nicely sorted list."""
56
62
        keys = super(FKSMultiProcess, self).get_sorted_keys()
57
63
        keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr', 
58
 
                 'has_fsr', 'OLP', 'ncores_for_proc_gen']
 
64
                 'has_fsr', 'spltting_types', 'OLP', 'ncores_for_proc_gen']
59
65
        return keys
60
66
 
61
67
    def filter(self, name, value):
63
69
 
64
70
        if name == 'born_processes':
65
71
            if not isinstance(value, FKSProcessList):
66
 
                raise self.PhysicsObjectError, \
67
 
                        "%s is not a valid list for born_processes " % str(value)                             
 
72
                raise self.PhysicsObjectError("%s is not a valid list for born_processes " % str(value))                             
68
73
 
69
74
        if name == 'real_amplitudes':
70
75
            if not isinstance(value, diagram_generation.AmplitudeList):
71
 
                raise self.PhysicsObjectError, \
72
 
                        "%s is not a valid list for real_amplitudes " % str(value)
 
76
                raise self.PhysicsObjectError("%s is not a valid list for real_amplitudes " % str(value))
73
77
                                                  
74
78
        if name == 'real_pdgs':
75
79
            if not isinstance(value, list):
76
 
                raise self.PhysicsObjectError, \
77
 
                        "%s is not a valid list for real_amplitudes " % str(value)
 
80
                raise self.PhysicsObjectError("%s is not a valid list for real_amplitudes " % str(value))
78
81
        
79
82
        if name == 'OLP':
80
83
            if not isinstance(value,str):
81
 
                raise self.PhysicsObjectError, \
82
 
                    "%s is not a valid string for OLP " % str(value)
 
84
                raise self.PhysicsObjectError("%s is not a valid string for OLP " % str(value))
83
85
 
84
86
        if name == 'ncores_for_proc_gen':
85
87
            if not isinstance(value,int):
86
 
                raise self.PhysicsObjectError, \
87
 
                    "%s is not a valid value for ncores_for_proc_gen " % str(value)
 
88
                raise self.PhysicsObjectError("%s is not a valid value for ncores_for_proc_gen " % str(value))
88
89
                                                     
89
90
        return super(FKSMultiProcess,self).filter(name, value)
 
91
 
 
92
 
 
93
    def check_ij_confs(self):
 
94
        """check that there is no duplicate FKS ij configuration"""
 
95
        ijconfs_dict = {}
 
96
        for born in self['born_processes']:
 
97
            # the copy.copy is needed as duplicate configurations will be removed on the fly
 
98
            for real in copy.copy(born.real_amps):
 
99
                pdgs = ' '.join([ '%d' % pdg for pdg in real.pdgs])
 
100
                for info in copy.copy(real.fks_infos):
 
101
                    ij = [info['i'], info['j']]
 
102
                    try:
 
103
                        if ij in ijconfs_dict[pdgs]:
 
104
                            logger.debug('Duplicate FKS configuration found for %s : ij = %s' %
 
105
                                    (real.process.nice_string(), str(ij)))
 
106
                            #remove the configuration
 
107
                            born.real_amps[born.real_amps.index(real)].fks_infos.remove(info)
 
108
                        else:
 
109
                            ijconfs_dict[pdgs].append(ij)
 
110
                    except KeyError:
 
111
                        ijconfs_dict[pdgs] = [ij]
 
112
                # check if any FKS configuration remains for the real emission, otherwise
 
113
                # remove it
 
114
                if not born.real_amps[born.real_amps.index(real)].fks_infos:
 
115
                    logger.debug('Removing real %s from born %s' % \
 
116
                            (real.process.nice_string(), born.born_amp['process'].nice_string()))
 
117
                    born.real_amps.remove(real)
 
118
 
90
119
    
91
120
    def __init__(self, procdef=None, options={}):
92
121
        """Initializes the original multiprocess, then generates the amps for the 
95
124
        legs (stored in pdgs, so that they need to be generated only once and then reicycled
96
125
        """
97
126
 
 
127
 
 
128
        if 'nlo_mixed_expansion' in options:
 
129
            self['nlo_mixed_expansion'] = options['nlo_mixed_expansion']
 
130
            del options['nlo_mixed_expansion']
 
131
        else:
 
132
            self['nlo_mixed_expansion'] = True
 
133
 
 
134
 
98
135
        #swhich the other loggers off
99
136
        loggers_off = [logging.getLogger('madgraph.diagram_generation'), 
100
137
                       logging.getLogger('madgraph.loop_diagram_generation')]
102
139
        for logg in loggers_off:
103
140
            logg.setLevel(logging.WARNING)
104
141
        
 
142
        self['real_amplitudes'] = diagram_generation.AmplitudeList()
 
143
        self['pdgs'] = []
 
144
 
105
145
        # OLP option
106
146
        olp='MadLoop'
107
 
        if 'OLP' in options.keys():
 
147
        if 'OLP' in list(options.keys()):
108
148
            olp = options['OLP']
109
149
            del options['OLP']
110
150
 
 
151
        self['init_lep_split']=False
 
152
        if 'init_lep_split' in list(options.keys()):
 
153
            self['init_lep_split']=options['init_lep_split']
 
154
            del options['init_lep_split']
 
155
 
111
156
        ncores_for_proc_gen = 0
112
157
        # ncores_for_proc_gen has the following meaning
113
158
        #   0 : do things the old way
114
159
        #   > 0 use ncores_for_proc_gen
115
160
        #   -1 : use all cores
116
 
        if 'ncores_for_proc_gen' in options.keys():
 
161
        if 'ncores_for_proc_gen' in list(options.keys()):
117
162
            ncores_for_proc_gen = options['ncores_for_proc_gen']
118
163
            del options['ncores_for_proc_gen']
119
164
 
123
168
 
124
169
        except diagram_generation.NoDiagramException as error:
125
170
            # If no born, then this process most likely does not have any.
126
 
            raise NoBornException, "Born diagrams could not be generated for the "+\
 
171
            raise NoBornException("Born diagrams could not be generated for the "+\
127
172
               self['process_definitions'][0].nice_string().replace('Process',\
128
173
               'process')+". Notice that aMC@NLO does not handle loop-induced"+\
129
174
               " processes yet, but you can still use MadLoop if you want to "+\
130
175
               "only generate them."+\
131
 
               " For this, use the 'virt=' mode, without multiparticle labels."
 
176
               " For this, use the 'virt=' mode, without multiparticle labels.")
132
177
 
133
178
        self['OLP'] = olp
134
179
        self['ncores_for_proc_gen'] = ncores_for_proc_gen
135
 
            
 
180
 
136
181
        #check process definition(s):
137
182
        # a process such as g g > g g will lead to real emissions 
138
183
        #   (e.g: u g > u g g ) which will miss some corresponding born,
156
201
            for leg in procdef['legs']:
157
202
                if any([id in soft_particles for id in leg['ids']]) \
158
203
                        and sorted(leg['ids']) != soft_particles:
159
 
                    logger.warning(('%s can have real emission processes ' + \
160
 
            'which are not finite.\nTo avoid this, please use multiparticles ' + \
161
 
            'when generating the process and be sure to include all the following ' + \
162
 
            'particles in the multiparticle definition:\n %s' ) \
163
 
               % (procdef.nice_string(), soft_particles_string) )
 
204
                    logger.warning('Use of multiparticles is non-trivial for NLO '+ \
 
205
                         'process generation and depends on the orders included, '+ \
 
206
                         'the process considered, as well as the PDF set chosen. '+ \
 
207
                         'See appendix D of arXiv:1804.10017 [hep-ph] for some '+ \
 
208
                         'guidance.')
164
209
                    break
165
 
        for procdef in self['process_definitions']:
166
 
            procdef.set('orders', diagram_generation.MultiProcess.find_optimal_process_orders(procdef))
167
210
 
168
211
        amps = self.get('amplitudes')
 
212
 
 
213
        # get the list of leptons from the model, in order to discard 
 
214
        # lepton-initiated processes unless the init_lep_split flag is specified
 
215
        if self['process_definitions']:
 
216
            leptons = self['process_definitions'][0]['model'].get_lepton_pdgs()
 
217
        else:
 
218
            leptons = []
 
219
 
 
220
        #generate reals, but combine them after having combined the borns
169
221
        for i, amp in enumerate(amps):
 
222
            # skip amplitudes with two initial leptons unless the init_lep_split option is True
 
223
            if not self['init_lep_split'] and \
 
224
                    all([l['id'] in leptons \
 
225
                    for l in [ll for ll in amp.get('process').get('legs') if not ll['state']]]):
 
226
                logger.info(('Discarding process%s.\n  If you want to include it, set the \n' + \
 
227
                             '  \'include_lepton_initiated_processes\' option to True') % \
 
228
                             amp.get('process').nice_string().replace('Process', ''))
 
229
                continue
 
230
 
170
231
            logger.info("Generating FKS-subtracted matrix elements for born process%s (%d / %d)" \
171
 
                % (amp['process'].nice_string(print_weighted=False).replace(\
172
 
                                                                 'Process', ''),
173
 
                 i + 1, len(amps)))
 
232
                % (amp['process'].nice_string(print_weighted=False, print_perturbated=False).replace('Process', ''),
 
233
                   i + 1, len(amps)))
174
234
 
175
 
            born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'])
 
235
            born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'], \
 
236
                                   init_lep_split=self['init_lep_split'])
176
237
            self['born_processes'].append(born)
177
 
            born.generate_reals(self['pdgs'], self['real_amplitudes'])
 
238
 
 
239
            born.generate_reals(self['pdgs'], self['real_amplitudes'], combine = False)
 
240
 
 
241
            # finally combine the real amplitudes
 
242
            born.combine_real_amplitudes()
178
243
 
179
244
        if not self['ncores_for_proc_gen']:
180
245
            # old generation mode 
181
246
 
182
 
            born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
 
247
            born_pdg_list = [[l['id'] for l in born.get_leglist()] \
183
248
                    for born in self['born_processes'] ]
184
249
 
185
250
            for born in self['born_processes']:
186
251
                for real in born.real_amps:
187
252
                    real.find_fks_j_from_i(born_pdg_list)
188
253
            if amps:
189
 
                if self['process_definitions'][0].get('NLO_mode') == 'all':
 
254
                if self['process_definitions'][0].get('NLO_mode') in ['all']:
190
255
                    self.generate_virtuals()
191
256
                
192
257
                elif not self['process_definitions'][0].get('NLO_mode') in ['all', 'real','LOonly']:
230
295
        self['OLP'] = other['OLP']
231
296
        self['ncores_for_proc_gen'] = other['ncores_for_proc_gen']
232
297
 
 
298
 
233
299
    def get_born_amplitudes(self):
234
300
        """return an amplitudelist with the born amplitudes"""
235
 
        return diagram_generation.AmplitudeList([born.born_amp \
236
 
                for born in self['born_processes']])
 
301
        return diagram_generation.AmplitudeList([
 
302
                born.born_amp for \
 
303
                born in self['born_processes']])
 
304
 
237
305
 
238
306
    def get_virt_amplitudes(self):
239
307
        """return an amplitudelist with the virt amplitudes"""
240
308
        return diagram_generation.AmplitudeList([born.virt_amp \
241
309
                for born in self['born_processes'] if born.virt_amp])
242
310
 
 
311
 
243
312
    def get_real_amplitudes(self):
244
313
        """return an amplitudelist with the real amplitudes"""
245
314
        return self.get('real_amplitudes')
257
326
                                     '%s at the output stage only.'%self['OLP'])
258
327
            return
259
328
 
260
 
        # determine the orders to be used to generate the loop
261
 
        loop_orders = {}
262
 
        for  born in self['born_processes']:
263
 
            for coup, val in fks_common.find_orders(born.born_amp).items():
264
 
                try:
265
 
                    loop_orders[coup] = max([loop_orders[coup], val])
266
 
                except KeyError:
267
 
                    loop_orders[coup] = val
 
329
        if not self['nlo_mixed_expansion']:
 
330
            # determine the orders to be used to generate the loop
 
331
            loop_orders = {}
 
332
            for  born in self['born_processes']:
 
333
                for coup, val in fks_common.find_orders(born.born_amp).items():
 
334
                    try:
 
335
                        loop_orders[coup] = max([loop_orders[coup], val])
 
336
                    except KeyError:
 
337
                        loop_orders[coup] = val
 
338
 
268
339
 
269
340
        for i, born in enumerate(self['born_processes']):
270
 
            logger.info('Generating virtual matrix elements using MadLoop:')
271
 
            myproc = copy.copy(born.born_proc)
 
341
            myproc = copy.copy(born.born_amp['process'])
 
342
            #misc.sprint(born.born_proc)
 
343
            #misc.sprint(myproc.input_string())
 
344
            #misc.sprint(myproc['orders'])
 
345
            # if [orders] are not specified, then
 
346
            # include all particles in the loops
 
347
            # i.e. allow all orders to be perturbed
 
348
            # (this is the case for EW corrections, where only squared oders 
 
349
            # are imposed)
 
350
            if not self['nlo_mixed_expansion']:
 
351
                myproc['orders'] = loop_orders
 
352
            elif not myproc['orders']:
 
353
                    myproc['perturbation_couplings'] = myproc['model']['coupling_orders']
272
354
            # take the orders that are actually used bu the matrix element
273
 
            myproc['orders'] = loop_orders
274
355
            myproc['legs'] = fks_common.to_legs(copy.copy(myproc['legs']))
275
356
            logger.info('Generating virtual matrix element with MadLoop for process%s (%d / %d)' \
276
 
                    % (myproc.nice_string(print_weighted = False).replace(\
 
357
                    % (myproc.nice_string(print_weighted= False, print_perturbated= False).replace(\
277
358
                                                             'Process', ''),
278
359
                        i + 1, len(self['born_processes'])))
279
 
            myamp = loop_diagram_generation.LoopAmplitude(myproc)
280
 
            if myamp.get('diagrams'):
 
360
            try:
 
361
                myamp = loop_diagram_generation.LoopAmplitude(myproc)
281
362
                born.virt_amp = myamp
 
363
            except InvalidCmd:
 
364
                logger.debug('invalid command for loop')
 
365
                pass
282
366
 
283
367
 
284
368
class FKSRealProcess(object): 
286
370
    -- fks_infos (list containing the possible fks configs for a given process
287
371
    -- amplitude 
288
372
    -- is_to_integrate
289
 
    -- leg permutation<<REMOVED!.
290
373
    """
291
374
    
292
 
    def __init__(self, born_proc, leglist, ij, ijglu,
 
375
    def __init__(self, born_proc, leglist, ij, ij_id, born_pdgs, splitting_type,
293
376
                 perturbed_orders = ['QCD']): #test written
294
377
        """Initializes the real process based on born_proc and leglist.
295
378
        Stores the fks informations into the list of dictionaries fks_infos
296
379
        """      
 
380
        #safety check
 
381
        assert type(splitting_type) == list and not type(splitting_type) == str 
297
382
        self.fks_infos = []
298
383
        for leg in leglist:
299
384
            if leg.get('fks') == 'i':
301
386
                # i is a gluon or a photon
302
387
                need_color_links = leg.get('massless') \
303
388
                        and leg.get('spin') == 3 \
304
 
                        and leg.get('self_antipart')
 
389
                        and leg.get('self_antipart') \
 
390
                        and leg.get('color') == 8
 
391
                need_charge_links = leg.get('massless') \
 
392
                        and leg.get('spin') == 3 \
 
393
                        and leg.get('self_antipart') \
 
394
                        and leg.get('color') == 1
305
395
            if leg.get('fks') == 'j':
306
396
                j_fks = leg.get('number')
307
 
        self.fks_infos.append({'i' : i_fks, 
308
 
                               'j' : j_fks, 
309
 
                               'ij' : ij, 
310
 
                               'ij_glu': ijglu, 
311
 
                               'need_color_links' : need_color_links})
 
397
        self.fks_infos.append({'i': i_fks, 
 
398
                               'j': j_fks, 
 
399
                               'ij': ij, 
 
400
                               'ij_id': ij_id, 
 
401
                               'underlying_born': born_pdgs,
 
402
                               'splitting_type': splitting_type,
 
403
                               'need_color_links': need_color_links,
 
404
                               'need_charge_links': need_charge_links,
 
405
                               'extra_cnt_index': -1})
312
406
 
313
407
        self.process = copy.copy(born_proc)
314
 
        orders = copy.copy(born_proc.get('orders'))
315
 
        # compute the weighted order if not present
316
 
        if not 'WEIGHTED' in orders:
317
 
            orders['WEIGHTED'] = sum([v * born_proc.get('model').get('order_hierarchy')[o] \
318
 
                        for o, v in orders.items()])
319
 
 
320
 
        for order in perturbed_orders:
321
 
            try:
322
 
                orders[order] +=1
323
 
            except KeyError:
324
 
                pass
325
 
            orders['WEIGHTED'] += born_proc.get('model').get('order_hierarchy')[order]
326
 
 
327
 
        self.process.set('orders', orders)
 
408
        self.process['perturbation_couplings'] = \
 
409
                copy.copy(born_proc['perturbation_couplings'])
 
410
        for o in splitting_type:
 
411
            if o not in self.process['perturbation_couplings']:
 
412
                self.process['perturbation_couplings'].append(o)
 
413
        # set the orders to empty, to force the use of the squared_orders
 
414
        self.process['orders'] = copy.copy(born_proc['orders'])
 
415
 
328
416
        legs = [(leg.get('id'), leg) for leg in leglist]
329
417
        self.pdgs = array.array('i',[s[0] for s in legs])
330
 
        if 'QCD' in perturbed_orders:
331
 
            self.colors = [leg['color'] for leg in leglist]
332
 
            # in QCD, charges are irrelevant !
333
 
            self.charges = [0. for leg in leglist]
334
 
            self.perturbation = 'QCD'
335
 
        else:
336
 
            self.colors = [leg['color'] for leg in leglist]
 
418
        self.colors = [leg['color'] for leg in leglist]
 
419
        if not self.process['perturbation_couplings'] == ['QCD']:
337
420
            self.charges = [leg['charge'] for leg in leglist]
338
 
            self.perturbation = 'QED'
 
421
        else:
 
422
            self.charges = [0.] * len(leglist)
 
423
        self.perturbation = 'QCD'
339
424
        self.process.set('legs', MG.LegList(leglist))
340
425
        self.process.set('legs_with_decays', MG.LegList())
341
426
        self.amplitude = diagram_generation.Amplitude()
342
427
        self.is_to_integrate = True
343
428
        self.is_nbody_only = False
344
429
        self.fks_j_from_i = {}
 
430
        self.missing_borns = []
345
431
 
346
432
 
347
433
    def generate_real_amplitude(self):
354
440
        """Returns a dictionary with the entries i : [j_from_i], if the born pdgs are in 
355
441
        born_pdg_list"""
356
442
        fks_j_from_i = {}
357
 
        dict = {}
358
443
        for i in self.process.get('legs'):
359
444
            fks_j_from_i[i.get('number')] = []
360
445
            if i.get('state'):
361
446
                for j in [l for l in self.process.get('legs') if \
362
447
                        l.get('number') != i.get('number')]:
363
 
                    ijlist = fks_common.combine_ij(i, j, self.process.get('model'), dict,\
364
 
                                                   pert=self.perturbation)
365
 
                    for ij in ijlist:
366
 
                        born_leglist = fks_common.to_fks_legs(
367
 
                                      copy.deepcopy(self.process.get('legs')), 
368
 
                                      self.process.get('model'))
369
 
                        born_leglist.remove(i)
370
 
                        born_leglist.remove(j)
371
 
                        born_leglist.insert(ij.get('number') - 1, ij)
372
 
                        born_leglist.sort(pert = self.perturbation)
373
 
                        if [l['id'] for l in born_leglist] in born_pdg_list:
374
 
                            fks_j_from_i[i.get('number')].append(\
375
 
                                                    j.get('number'))                                
 
448
                    for pert_order in self.process.get('perturbation_couplings'):
 
449
                        ijlist = fks_common.combine_ij(i, j, self.process.get('model'), {},\
 
450
                                                       pert=pert_order)
 
451
                        for ij in ijlist:
 
452
                            born_leglist = fks_common.to_fks_legs(
 
453
                                          copy.deepcopy(self.process.get('legs')), 
 
454
                                          self.process.get('model'))
 
455
                            born_leglist.remove(i)
 
456
                            born_leglist.remove(j)
 
457
                            born_leglist.insert(ij.get('number') - 1, ij)
 
458
                            born_leglist.sort(pert = self.perturbation)
 
459
                            if [leg['id'] for leg in born_leglist] in born_pdg_list \
 
460
                               and not j.get('number') in fks_j_from_i[i.get('number')]:
 
461
                                fks_j_from_i[i.get('number')].append(\
 
462
                                                        j.get('number'))                                
376
463
 
377
464
        self.fks_j_from_i = fks_j_from_i
378
465
        return fks_j_from_i
406
493
class FKSProcess(object):
407
494
    """The class for a FKS process. Starts from the born process and finds
408
495
    all the possible splittings."""  
 
496
 
 
497
 
 
498
#helper functions
 
499
 
 
500
    def get_colors(self):
 
501
        """return the list of color representations 
 
502
        for each leg in born_amp"""
 
503
        return [leg.get('color') for \
 
504
                    leg in self.born_amp['process']['legs']]                    
 
505
 
 
506
 
 
507
    def get_charges(self):
 
508
        """return the list of charges
 
509
        for each leg in born_amp"""
 
510
        return [leg.get('charge') for \
 
511
                    leg in self.born_amp['process']['legs']]                    
 
512
 
 
513
 
 
514
    def get_nlegs(self):
 
515
        """return the number of born legs"""
 
516
        return len(self.born_amp['process']['legs'])
 
517
 
 
518
 
 
519
    def get_born_nice_string(self):
 
520
        """Return the nice string for the born process.
 
521
        """
 
522
        return self.born_amp['process'].nice_string()
 
523
 
 
524
 
 
525
    def get_pdg_codes(self):
 
526
        """return the list of the pdg codes
 
527
        of each leg in born_amp"""
 
528
        return [leg.get('id') for \
 
529
                    leg in self.born_amp['process']['legs']]                    
 
530
 
 
531
 
 
532
    def get_leglist(self):
 
533
        """return the leg list
 
534
        for the born amp"""
 
535
        return fks_common.to_fks_legs( \
 
536
                self.born_amp['process']['legs'], \
 
537
                self.born_amp['process']['model'])
 
538
 
 
539
 
 
540
###############################################################################
409
541
    
410
 
    def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0):
 
542
    def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0, init_lep_split = False):
411
543
        """initialization: starts either from an amplitude or a process,
412
544
        then init the needed variables.
413
545
        remove_borns tells if the borns not needed for integration will be removed
418
550
           -1 : use all cores
419
551
        """
420
552
                
421
 
        self.splittings = {}
422
553
        self.reals = []
423
 
        self.fks_dirs = []
424
 
        self.leglist = []
425
554
        self.myorders = {}
426
 
        self.pdg_codes = []
427
 
        self.colors = [] # color
428
 
        self.charges = [] # charge
429
 
        self.nlegs = 0
430
 
        self.fks_ipos = []
431
 
        self.fks_j_from_i = {}
432
555
        self.real_amps = []
433
556
        self.remove_reals = remove_reals
 
557
        self.init_lep_split = init_lep_split
434
558
        self.nincoming = 0
435
559
        self.virt_amp = None
436
560
        self.perturbation = 'QCD'
 
561
        self.born_amp = diagram_generation.Amplitude()
 
562
        self.extra_cnt_amp_list = diagram_generation.AmplitudeList()
437
563
        self.ncores_for_proc_gen = ncores_for_proc_gen
438
564
 
439
565
        if not remove_reals in [True, False]:
441
567
                    'Not valid type for remove_reals in FKSProcess')
442
568
 
443
569
        if start_proc:
 
570
            #initilaize with process definition (for test purporses)
444
571
            if isinstance(start_proc, MG.Process):
445
572
                pertur = start_proc['perturbation_couplings']
446
573
                if pertur:
447
574
                    self.perturbation = sorted(pertur)[0]
448
 
                self.born_proc = fks_common.sort_proc(start_proc,pert = self.perturbation)
449
 
                # filter in Amplitude will legs sorted in bornproc
450
 
                bornproc = copy.copy(self.born_proc) # deepcopy might let T -> array
451
 
                assert bornproc==self.born_proc
452
 
                self.born_amp = diagram_generation.Amplitude(bornproc)
 
575
                self.born_amp = diagram_generation.Amplitude(\
 
576
                                copy.copy(fks_common.sort_proc(\
 
577
                                        start_proc, pert = self.perturbation)))
 
578
            #initialize with an amplitude
453
579
            elif isinstance(start_proc, diagram_generation.Amplitude):
454
580
                pertur = start_proc.get('process')['perturbation_couplings']
455
 
                if pertur:
456
 
                    self.perturbation = sorted(pertur)[0]
457
 
                self.born_proc = fks_common.sort_proc(start_proc.get('process'),\
458
 
                                                      pert = self.perturbation)
459
 
                # filter in Amplitude will legs sorted in bornproc
460
 
                bornproc = copy.copy(self.born_proc)
461
 
                assert bornproc == self.born_proc
462
 
                self.born_amp = diagram_generation.Amplitude(bornproc)
 
581
                self.born_amp = diagram_generation.Amplitude(\
 
582
                                copy.copy(fks_common.sort_proc(\
 
583
                                    start_proc['process'], 
 
584
                                    pert = self.perturbation)))
463
585
            else:
464
586
                raise fks_common.FKSProcessError(\
465
587
                    'Not valid start_proc in FKSProcess')
466
 
            self.born_proc.set('legs_with_decays', MG.LegList())
 
588
            self.born_amp['process'].set('legs_with_decays', MG.LegList())
467
589
 
468
 
            self.leglist = fks_common.to_fks_legs(
469
 
                                    self.born_proc['legs'], self.born_proc['model'])
470
 
            self.nlegs = len(self.leglist)
471
 
            self.pdg_codes = [leg.get('id') for leg in self.leglist]
472
 
            if self.perturbation == 'QCD':
473
 
                self.colors = [leg.get('color') for leg in self.leglist]
474
 
                # in QCD, charge is irrelevant but impact the combine process !
475
 
                self.charges = [0. for leg in self.leglist]
476
 
                color = 'color'
477
 
                zero = 1
478
 
            elif self.perturbation == 'QED':
479
 
                self.colors = [leg.get('color') for leg in self.leglist]
480
 
                self.charges = [leg.get('charge') for leg in self.leglist]
481
 
                color = 'charge'
482
 
                zero = 0.
483
590
            # special treatment of photon is needed !
484
 
            self.isr = set([leg.get(color) for leg in self.leglist if not leg.get('state')]) != set([zero])
485
 
            self.fsr = set([leg.get(color) for leg in self.leglist if leg.get('state')]) != set([zero])
486
 
            for leg in self.leglist:
487
 
                if not leg['state']:
488
 
                    self.nincoming += 1
489
 
            self.orders = self.born_amp['process']['orders']
490
 
            # this is for cases in which the user specifies e.g. QED=0
491
 
            if sum(self.orders.values()) == 0:
492
 
                self.orders = fks_common.find_orders(self.born_amp)
 
591
            #MZ to be fixed
 
592
            ###self.isr = set([leg.get(color) for leg in self.leglist if not leg.get('state')]) != set([zero])
 
593
            ###self.fsr = set([leg.get(color) for leg in self.leglist if leg.get('state')]) != set([zero])
 
594
            self.isr = False
 
595
            self.fsr = False
 
596
            #######
 
597
            self.nincoming = len([l for l in self.born_amp['process']['legs'] \
 
598
                                  if not l['state']])
493
599
                
494
600
            self.ndirs = 0
495
601
            # generate reals, when the mode is not LOonly
496
602
            # when is LOonly it is supposed to be a 'fake' NLO process
497
603
            # e.g. to be used in merged sampels at high multiplicities
498
 
            if self.born_proc['NLO_mode'] != 'LOonly':
499
 
                for order in self.born_proc.get('perturbation_couplings'):
500
 
                    self.find_reals(order)
 
604
            if self.born_amp['process']['NLO_mode'] != 'LOonly':
 
605
                self.find_reals()
501
606
 
502
607
 
503
608
    def generate_real_amplitudes(self, pdg_list, real_amp_list):
504
609
        """generates the real amplitudes for all the real emission processes, using pdgs and real_amps
505
 
        to avoid multiple generation of the same amplitude"""
 
610
        to avoid multiple generation of the same amplitude.
 
611
        Amplitude without diagrams are discarded at this stage"""
506
612
 
 
613
        no_diags_amps = []
507
614
        for amp in self.real_amps:
508
615
            try:
509
616
                amp.amplitude = real_amp_list[pdg_list.index(amp.pdgs)]
510
617
            except ValueError:
511
 
                pdg_list.append(amp.pdgs)
512
 
                real_amp_list.append(amp.generate_real_amplitude())
 
618
                amplitude = amp.generate_real_amplitude()
 
619
                if amplitude['diagrams']:
 
620
                    pdg_list.append(amp.pdgs)
 
621
                    real_amp_list.append(amplitude)
 
622
                else:
 
623
                    no_diags_amps.append(amp)
 
624
        
 
625
        for amp in no_diags_amps:
 
626
            self.real_amps.remove(amp)
 
627
 
513
628
 
514
629
 
515
630
    def combine_real_amplitudes(self):
535
650
        if combine is true, FKS_real_processes having the same pdgs (i.e. real amplitude)
536
651
        are combined together
537
652
        """
538
 
 
539
 
        born_proc = copy.copy(self.born_proc)
540
 
        born_proc['orders'] = self.orders
541
 
        for i, list in enumerate(self.reals):
542
 
            # i is a gluon or photon,i.e. g(a) > ffx or gg
543
 
            if self.leglist[i]['massless'] and self.leglist[i]['spin'] == 3:
544
 
                ijglu = i + 1
545
 
            else:
546
 
                ijglu = 0
547
 
            for l in list:
548
 
                ij = self.leglist[i].get('number')
 
653
        #copy the born process
 
654
        born_proc = copy.copy(self.born_amp['process'])
 
655
        born_pdgs = self.get_pdg_codes()
 
656
        leglist = self.get_leglist()
 
657
        extra_cnt_pdgs = []
 
658
        for i, real_list in enumerate(self.reals):
 
659
            # i is the born leg which splits
 
660
            # keep track of the id of the mother (will be used to constrct the
 
661
            # spin-correlated borns)
 
662
            ij_id = leglist[i].get('id')
 
663
            ij = leglist[i].get('number')
 
664
            for real_dict in real_list:
 
665
                nmom = 0
 
666
 
 
667
                # check first if other counterterms need to be generated
 
668
                # (e.g. g/a > q qbar)
 
669
                # this is quite a tricky business, as double counting
 
670
                # the singular configuration must be avoided. 
 
671
                # Let real, born, cnt be the real emission, the born process
 
672
                # (that will give the name to the P0_** dir) and the 
 
673
                # extra counterterm (obtained by the born process replacing
 
674
                # ij with the extra mother). 
 
675
                # If there are extra mothers, first check that
 
676
                # 1) born, at order born[squared_orders] - 
 
677
                #    2 * (the perturbation type of the real emission) has diagrams
 
678
                # 2) cnt at order born[squared_orders] - 
 
679
                #    2 * (the perturbation type of the extra mom) has diagrams
 
680
 
 
681
                cnt_amp = diagram_generation.Amplitude()
 
682
                born_cnt_amp = diagram_generation.Amplitude()
 
683
                mom_cnt = 0
 
684
                cnt_ord = None
 
685
 
 
686
                # check condition 1) above (has_coll_sing_born)
 
687
                born_proc_coll_sing = copy.copy(born_proc)
 
688
                born_proc_coll_sing['squared_orders'] = copy.copy(born_proc['squared_orders'])
 
689
                if born_proc_coll_sing['squared_orders'][real_dict['perturbation'][0]] < 2:
 
690
                    has_coll_sing_born = False
 
691
                else:
 
692
                    born_proc_coll_sing['squared_orders'][real_dict['perturbation'][0]] += -2
 
693
                    has_coll_sing_born = bool(diagram_generation.Amplitude(born_proc_coll_sing)['diagrams'])
 
694
 
 
695
                # check that there is at most one extra mother
 
696
                allmothers = []
 
697
                for order, mothers in real_dict['extra_mothers'].items():
 
698
                    allmothers += mothers
 
699
                    if mothers:
 
700
                        cnt_ord = order
 
701
 
 
702
                if len(allmothers) > 1:
 
703
                    raise fks_common.FKSProcessError(\
 
704
                            'Error, more than one extra mother has been found: %d', len(allmothers))
 
705
                # here we are sure to have just one extra mother
 
706
                    
 
707
                has_coll_sing_cnt = False
 
708
                if allmothers:
 
709
                    mom_cnt = allmothers[0]
 
710
 
 
711
                    # generate a new process with the mother particle 
 
712
                    # replaced by the new mother and with the
 
713
                    # squared orders changed accordingly
 
714
 
 
715
                    cnt_process = copy.copy(born_proc)
 
716
                    cnt_process['legs'] = copy.deepcopy(born_proc['legs'])
 
717
                    cnt_process['legs'][i]['id'] = mom_cnt
 
718
                    cnt_process['legs'] = fks_common.to_fks_legs(
 
719
                            cnt_process['legs'], cnt_process['model'])
 
720
                    cnt_process['squared_orders'] = \
 
721
                            copy.copy(born_proc['squared_orders'])
 
722
 
 
723
                    # check if the cnt amplitude exists with the current 
 
724
                    # squared orders (i.e. if it will appear as a P0 dir)
 
725
                    # if it does not exist, then no need to worry about anything, as all further
 
726
                    # checks will have stricter orders than here
 
727
                    
 
728
                    cnt_process_for_amp = copy.copy(cnt_process)
 
729
                    cnt_process_for_amp['squared_orders'] = copy.copy(cnt_process['squared_orders'])
 
730
                    cnt_amp = diagram_generation.Amplitude(cnt_process_for_amp)
 
731
 
 
732
                    if bool(cnt_amp['diagrams']) and \
 
733
                       cnt_process['squared_orders'][cnt_ord] >= 2:
 
734
 
 
735
                        # check condition 2) above (has_coll_sing_cnt)
 
736
                        # MZMZ17062014 beware that the Amplitude reorders the legs
 
737
                        cnt_process['squared_orders'][cnt_ord] += -2
 
738
                        born_cnt_amp = diagram_generation.Amplitude(cnt_process)
 
739
                        has_coll_sing_cnt = bool(born_cnt_amp['diagrams'])
 
740
 
 
741
                # remember there is at most one mom
 
742
                # now, one of these cases can occur
 
743
                #  a) no real collinear singularity exists (e.g. just the interference
 
744
                #     has to be integrated for the real emission). Add the real process to
 
745
                #     this born process if ij_id < mom_cnt. No extra_cnt is needed
 
746
                #  b) a collinear singularity exists, with the underlying born being the born
 
747
                #     process of this dir, while no singularity is there for the extra_cnt.
 
748
                #     In this case keep the real emission in this directory, without any 
 
749
                #     extra_cnt
 
750
                #  c) a collinear singularity exists, with the underlying born being the 
 
751
                #     extra_cnt, while no singularity is there for the born of the dir.
 
752
                #     In this case skip the real emission, it will be included in the
 
753
                #     directory of the extra cnt 
 
754
                #  d) a collinear singularity exists, and both the process-dir born and
 
755
                #     the extra cnt are needed to subtract it. Add the real process to
 
756
                #     this born process if ij_id < mom_cnt and keeping the extra_cnt 
 
757
                #
 
758
                # in all cases, remember that mom_cnt is set to 0 if no extra mother is there
 
759
 
 
760
                # the real emission has to be skipped if mom_cnt(!=0) < ij_id and either a) or d)
 
761
                if mom_cnt and mom_cnt < ij_id:
 
762
                    if ((not has_coll_sing_born and not has_coll_sing_cnt) or \
 
763
                        (has_coll_sing_born and has_coll_sing_cnt)):
 
764
                        continue
 
765
 
 
766
                # the real emission has also to be skipped in case of c)
 
767
                if has_coll_sing_cnt and not has_coll_sing_born:
 
768
                    continue
 
769
 
 
770
                # if we arrive here, we need to keep this real mession
 
771
                ij = leglist[i].get('number')
549
772
                self.real_amps.append(FKSRealProcess( \
550
 
                        born_proc, l, ij, ijglu,\
551
 
                        perturbed_orders = [self.perturbation]))
 
773
                        born_proc, real_dict['leglist'], ij, ij_id, \
 
774
                        [born_pdgs], 
 
775
                        real_dict['perturbation'], \
 
776
                        perturbed_orders = born_proc['perturbation_couplings']))
 
777
 
 
778
                # keep the extra_cnt if needed
 
779
                if has_coll_sing_cnt:
 
780
                    # for the moment we just check the pdgs, regardless of any
 
781
                    # permutation in the final state
 
782
                    try:
 
783
                        indx = extra_cnt_pdgs.index([l['id'] for l in cnt_process['legs']])
 
784
                    except ValueError:
 
785
                        extra_cnt_pdgs.append([l['id'] for l in cnt_process['legs']])
 
786
                        assert cnt_amp != None
 
787
                        self.extra_cnt_amp_list.append(cnt_amp)
 
788
                        indx = len(self.extra_cnt_amp_list) - 1
 
789
                      
 
790
                    # update the fks infos
 
791
                    self.real_amps[-1].fks_infos[-1]['extra_cnt_index'] = indx
 
792
                    self.real_amps[-1].fks_infos[-1]['underlying_born'].append(\
 
793
                                            [l['id'] for l in cnt_process['legs']])
 
794
                    self.real_amps[-1].fks_infos[-1]['splitting_type'].append(cnt_ord)
 
795
 
552
796
        self.find_reals_to_integrate()
553
797
        if combine:
554
798
            self.combine_real_amplitudes()
561
805
        """create the rb_links in the real matrix element to find 
562
806
        which configuration in the real correspond to which in the born
563
807
        """
 
808
 
 
809
        # check that all splitting types are of ['QCD'] type, otherwise return
 
810
        for real in self.real_amps:
 
811
            for info in real.fks_infos:
 
812
                if info['splitting_type'] != ['QCD']:
 
813
                    logger.info('link_born_real: skipping because not all splittings are QCD')
 
814
                    return
 
815
 
 
816
 
564
817
        for real in self.real_amps:
565
818
            for info in real.fks_infos:
566
819
                info['rb_links'] = fks_common.link_rb_configs(\
568
821
                        info['i'], info['j'], info['ij'])
569
822
 
570
823
 
571
 
    def find_reals(self, pert_order):
572
 
        """finds the FKS real configurations for a given process"""
573
 
        if range(len(self.leglist)) != [l['number']-1 for l in self.leglist]:
 
824
 
 
825
    def find_reals(self, pert_orders = []):
 
826
        """finds the FKS real configurations for a given process.
 
827
        self.reals[i] is a list of dictionaries corresponding to the real 
 
828
        emissions obtained splitting leg i.
 
829
        The dictionaries contain the leglist, the type (order) of the
 
830
        splitting and extra born particles which can give the same
 
831
        splitting (e.g. gluon/photon -> qqbar).
 
832
        If pert orders is empty, all the orders of the model will be used
 
833
        """
 
834
        
 
835
        model = self.born_amp['process']['model']
 
836
        # if [orders] are not specified, then
 
837
        # include all kind of splittings
 
838
        # i.e. allow all orders to be perturbed
 
839
        # (this is the case for EW corrections, where only squared oders 
 
840
        # are imposed)
 
841
        if not pert_orders:
 
842
            if not self.born_amp['process']['orders']:
 
843
                pert_orders = model['coupling_orders']
 
844
            else:
 
845
                pert_orders = self.born_amp['process']['perturbation_couplings']
 
846
 
 
847
        leglist = self.get_leglist()
 
848
        if list(range(len(leglist))) != [l['number']-1 for l in leglist]:
574
849
            raise fks_common.FKSProcessError('Disordered numbers of leglist')
575
850
 
576
 
        if [ i['state'] for i in self.leglist].count(False) == 1:
 
851
        if [ i['state'] for i in leglist].count(False) == 1:
577
852
            decay_process=True
578
853
        else:
579
854
            decay_process=False
580
855
 
581
 
        for i in self.leglist:
 
856
        # count the number of initial-state leptons
 
857
        ninit_lep = [l['id'] in model.get_lepton_pdgs() and not l['state'] for l in leglist].count(True)
 
858
 
 
859
        for i in leglist:
582
860
            i_i = i['number'] - 1
583
861
            self.reals.append([])
584
 
            if decay_process and not i['state']:
585
 
                self.splittings[i_i]=[]
586
 
            else:
587
 
                self.splittings[i_i] = fks_common.find_splittings(i, self.born_proc['model'], {}, pert_order)
588
 
            for split in self.splittings[i_i]:
589
 
                self.reals[i_i].append(
590
 
                            fks_common.insert_legs(self.leglist, i, split,pert=pert_order))
591
 
                
 
862
 
 
863
            # for 2->1 processes, map only the initial-state singularities
 
864
            # this is because final-state mapping preserves shat, which
 
865
            # is not possible in 2->1
 
866
            if len(leglist) == 3 and not decay_process and i['state']: 
 
867
                continue
 
868
            for pert_order in pert_orders:
 
869
                # no splittings for initial states in decay processes
 
870
                if decay_process and not i['state']:
 
871
                    splittings=[]
 
872
                # if there are leptons in the initial state and init_lep_split is False,
 
873
                # only split initial state leptons; do nothing for any other particle
 
874
                elif not self.init_lep_split and ninit_lep >= 1 and \
 
875
                  (i['state'] or i['id'] not in model.get_lepton_pdgs()):
 
876
                    splittings=[]
 
877
                else:
 
878
                    splittings = fks_common.find_splittings( \
 
879
                        i, model, {}, pert_order, \
 
880
                        include_init_leptons=self.init_lep_split)
 
881
                for split in splittings:
 
882
                    # find other 'mother' particles which can end up in the same splitting
 
883
                    extra_mothers = {}
 
884
                    for pert in pert_orders:
 
885
                        extra_mothers[pert] = fks_common.find_mothers(split[0], split[1], model, pert=pert,
 
886
                                        mom_mass=model.get('particle_dict')[i['id']]['mass'].lower())
 
887
 
 
888
                    #remove the current mother from the extra mothers    
 
889
                    if i['state']:
 
890
                        extra_mothers[pert_order].remove(i['id'])
 
891
                    else:
 
892
                        extra_mothers[pert_order].remove(model.get('particle_dict')[i['id']].get_anti_pdg_code())
 
893
 
 
894
                    self.reals[i_i].append({
 
895
                        'leglist': fks_common.insert_legs(leglist, i, split ,pert=pert_order),
 
896
                        'perturbation': [pert_order], 
 
897
                        'extra_mothers': extra_mothers})
592
898
 
593
899
 
594
900
    def find_reals_to_integrate(self): #test written
610
916
 
611
917
                i_m = real_m.fks_infos[0]['i']
612
918
                j_m = real_m.fks_infos[0]['j']
613
 
                i_n = real_n.fks_infos[0]['i']
 
919
                i_n = real_n.fks_infos[-1]['i']
614
920
                j_n = real_n.fks_infos[0]['j']
 
921
                ij_id_m = real_m.fks_infos[0]['ij_id']
 
922
                ij_id_n = real_n.fks_infos[0]['ij_id']
615
923
                if j_m > self.nincoming and j_n > self.nincoming:
 
924
                    # make sure i and j in the two real emissions have the same mother 
 
925
                    if (ij_id_m != ij_id_n):
 
926
                        continue
616
927
                    if (real_m.get_leg_i()['id'] == real_n.get_leg_i()['id'] \
617
928
                        and \
618
929
                        real_m.get_leg_j()['id'] == real_n.get_leg_j()['id']) \
621
932
                        and \
622
933
                        real_m.get_leg_j()['id'] == real_n.get_leg_i()['id']):
623
934
                        if i_m > i_n:
624
 
                            print real_m.get_leg_i()['id'], real_m.get_leg_j()['id']
625
935
                            if real_m.get_leg_i()['id'] == -real_m.get_leg_j()['id']:
626
936
                                self.real_amps[m].is_to_integrate = False
627
937
                            else:
628
938
                                self.real_amps[n].is_to_integrate = False
629
939
                        elif i_m == i_n and j_m > j_n:
630
 
                            print real_m.get_leg_i()['id'], real_m.get_leg_j()['id']
631
940
                            if real_m.get_leg_i()['id'] == -real_m.get_leg_j()['id']:
632
941
                                self.real_amps[m].is_to_integrate = False
633
942
                            else: