1
################################################################################
3
# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors
5
# This file is a part of the MadGraph5_aMC@NLO project, an application which
6
# automatically generates Feynman diagrams and matrix elements for arbitrary
7
# high-energy processes in the Standard Model and beyond.
9
# It is subject to the MadGraph5_aMC@NLO license which should accompany this
12
# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
14
################################################################################
15
"""Definitions of all basic objects used in the core code: particle,
16
interaction, model, leg, vertex, process, ..."""
18
from __future__ import absolute_import
28
import madgraph.core.color_algebra as color
30
from madgraph import MadGraph5Error, MG5DIR, InvalidCmd
31
import madgraph.various.misc as misc
32
from six.moves import range
33
from six.moves import zip
34
from functools import reduce
37
logger = logging.getLogger('madgraph.base_objects')
40
#===============================================================================
42
#===============================================================================
43
class PhysicsObject(dict):
44
"""A parent class for all physics objects."""
46
class PhysicsObjectError(Exception):
47
"""Exception raised if an error occurs in the definition
48
or the execution of a physics object."""
51
def __init__(self, init_dict={}):
52
"""Creates a new particle object. If a dictionary is given, tries to
53
use it to give values to properties."""
58
assert isinstance(init_dict, dict), \
59
"Argument %s is not a dictionary" % repr(init_dict)
62
for item in init_dict.keys():
63
self.set(item, init_dict[item])
66
def __getitem__(self, name):
67
""" force the check that the property exist before returning the
68
value associated to value. This ensure that the correct error
73
return dict.__getitem__(self, name)
75
self.is_valid_prop(name) #raise the correct error
78
def default_setup(self):
79
"""Function called to create and setup default values for all object
83
def is_valid_prop(self, name):
84
"""Check if a given property name is valid"""
86
assert isinstance(name, str), \
87
"Property name %s is not a string" % repr(name)
89
if name not in list(self.keys()):
90
raise self.PhysicsObjectError("""%s is not a valid property for this object: %s\n
91
Valid property are %s""" % (name,self.__class__.__name__, list(self.keys())))
95
"""Get the value of the property name."""
99
def set(self, name, value, force=False):
100
"""Set the value of the property name. First check if value
101
is a valid value for the considered property. Return True if the
102
value has been correctly set, False otherwise."""
103
if not __debug__ or force:
107
if self.is_valid_prop(name):
109
self.filter(name, value)
112
except self.PhysicsObjectError as why:
113
logger.warning("Property " + name + " cannot be changed:" + \
117
def filter(self, name, value):
118
"""Checks if the proposed value is valid for a given property
119
name. Returns True if OK. Raises an error otherwise."""
123
def get_sorted_keys(self):
124
"""Returns the object keys sorted in a certain way. By default,
127
return list(self.keys()).sort()
130
"""String representation of the object. Outputs valid Python
131
with improved format."""
134
for prop in self.get_sorted_keys():
135
if isinstance(self[prop], str):
136
mystr = mystr + ' \'' + prop + '\': \'' + \
138
elif isinstance(self[prop], float):
139
mystr = mystr + ' \'' + prop + '\': %.2f,\n' % self[prop]
141
mystr = mystr + ' \'' + prop + '\': ' + \
142
repr(self[prop]) + ',\n'
143
mystr = mystr.rstrip(',\n')
144
mystr = mystr + '\n}'
151
#===============================================================================
153
#===============================================================================
154
class PhysicsObjectList(list):
155
"""A class to store lists of physics object."""
157
class PhysicsObjectListError(Exception):
158
"""Exception raised if an error occurs in the definition
159
or execution of a physics object list."""
162
def __init__(self, init_list=None):
163
"""Creates a new particle list object. If a list of physics
164
object is given, add them."""
168
if init_list is not None:
169
for object in init_list:
172
def append(self, object):
173
"""Appends an element, but test if valid before."""
175
assert self.is_valid_element(object), \
176
"Object %s is not a valid object for the current list" % repr(object)
178
list.append(self, object)
181
def is_valid_element(self, obj):
182
"""Test if object obj is a valid element for the list."""
186
"""String representation of the physics object list object.
187
Outputs valid Python with improved format."""
192
mystr = mystr + str(obj) + ',\n'
194
mystr = mystr.rstrip(',\n')
198
#===============================================================================
200
#===============================================================================
201
class Particle(PhysicsObject):
202
"""The particle object containing the whole set of information required to
203
univocally characterize a given type of physical particle: name, spin,
204
color, mass, width, charge,... The is_part flag tells if the considered
205
particle object is a particle or an antiparticle. The self_antipart flag
206
tells if the particle is its own antiparticle."""
208
sorted_keys = ['name', 'antiname', 'spin', 'color',
209
'charge', 'mass', 'width', 'pdg_code',
210
'line', 'propagator',
211
'is_part', 'self_antipart', 'type', 'counterterm']
213
def default_setup(self):
214
"""Default values for all properties"""
216
self['name'] = 'none'
217
self['antiname'] = 'none'
221
self['mass'] = 'ZERO'
222
self['width'] = 'ZERO'
224
#self['texname'] = 'none'
225
#self['antitexname'] = 'none'
226
self['line'] = 'dashed'
227
#self['propagating'] = True -> removed in favor or 'line' = None
228
self['propagator'] = ''
229
self['is_part'] = True
230
self['self_antipart'] = False
231
# True if ghost, False otherwise
232
#self['ghost'] = False
233
self['type'] = '' # empty means normal can also be ghost or goldstone
234
# Counterterm defined as a dictionary with format:
235
# ('ORDER_OF_COUNTERTERM',((Particle_list_PDG))):{laurent_order:CTCouplingName}
236
self['counterterm'] = {}
241
return self['type'] == 'ghost'
242
elif name == 'goldstone':
243
return self['type'] == 'goldstone'
244
elif name == 'propagating':
245
return self['line'] not in ['None',None]
247
return super(Particle, self).get(name)
249
def set(self, name, value, force=False):
251
if name in ['texname', 'antitexname']:
253
elif name == 'propagating':
255
return self.set('line', None, force=force)
256
elif not self.get('line'):
257
return self.set('line', 'dashed',force=force)
259
elif name in ['ghost', 'goldstone']:
260
if self.get('type') == name:
264
return self.set('type', '', force=force)
267
return self.set('type', name, force=force)
270
return super(Particle, self).set(name, value,force=force)
273
def filter(self, name, value):
274
"""Filter for valid particle property values."""
276
if name in ['name', 'antiname']:
277
# Forbid special character but +-~_
278
p=re.compile('''^[\w\-\+~_]+$''')
279
if not p.match(value):
280
raise self.PhysicsObjectError("%s is not a valid particle name" % value)
283
if not isinstance(value,bool):
284
raise self.PhysicsObjectError("%s is not a valid bool for the 'ghost' attribute" % str(value))
286
if name == 'counterterm':
287
if not isinstance(value,dict):
288
raise self.PhysicsObjectError("counterterm %s is not a valid dictionary" % repr(value))
289
for key, val in value.items():
290
if not isinstance(key,tuple):
291
raise self.PhysicsObjectError("key %s is not a valid tuple for counterterm key" % repr(key))
292
if not isinstance(key[0],str):
293
raise self.PhysicsObjectError("%s is not a valid string" % repr(key[0]))
294
if not isinstance(key[1],tuple):
295
raise self.PhysicsObjectError("%s is not a valid list" % repr(key[1]))
297
if not isinstance(elem,tuple):
298
raise self.PhysicsObjectError("%s is not a valid list" % repr(elem))
300
if not isinstance(partPDG,int):
301
raise self.PhysicsObjectError("%s is not a valid integer for PDG" % repr(partPDG))
303
raise self.PhysicsObjectError("%s is not a valid positive PDG" % repr(partPDG))
304
if not isinstance(val,dict):
305
raise self.PhysicsObjectError("value %s is not a valid dictionary for counterterm value" % repr(val))
306
for vkey, vvalue in val.items():
307
if vkey not in [0,-1,-2]:
308
raise self.PhysicsObjectError("Key %s is not a valid laurent serie order" % repr(vkey))
309
if not isinstance(vvalue,str):
310
raise self.PhysicsObjectError("Coupling %s is not a valid string" % repr(vvalue))
312
if not isinstance(value, int):
313
raise self.PhysicsObjectError("Spin %s is not an integer" % repr(value))
314
if (value < 1 or value > 5) and value != 99:
315
raise self.PhysicsObjectError("Spin %i not valid" % value)
318
if not isinstance(value, int):
319
raise self.PhysicsObjectError("Color %s is not an integer" % repr(value))
320
if value not in [1, 3, 6, 8]:
321
raise self.PhysicsObjectError("Color %i is not valid" % value)
323
if name in ['mass', 'width']:
324
# Must start with a letter, followed by letters, digits or _
325
p = re.compile('\A[a-zA-Z]+[\w\_]*\Z')
326
if not p.match(value):
327
raise self.PhysicsObjectError("%s is not a valid name for mass/width variable" % \
330
if name == 'pdg_code':
331
if not isinstance(value, int):
332
raise self.PhysicsObjectError("PDG code %s is not an integer" % repr(value))
335
if not isinstance(value, str):
336
raise self.PhysicsObjectError("Line type %s is not a string" % repr(value))
337
if value not in ['None','dashed', 'straight', 'wavy', 'curly', 'double','swavy','scurly','dotted']:
338
raise self.PhysicsObjectError("Line type %s is unknown" % value)
341
if not isinstance(value, float):
342
raise self.PhysicsObjectError("Charge %s is not a float" % repr(value))
344
if name == 'propagating':
345
if not isinstance(value, bool):
346
raise self.PhysicsObjectError("Propagating tag %s is not a boolean" % repr(value))
348
if name in ['is_part', 'self_antipart']:
349
if not isinstance(value, bool):
350
raise self.PhysicsObjectError("%s tag %s is not a boolean" % (name, repr(value)))
354
def get_sorted_keys(self):
355
"""Return particle property names as a nicely sorted list."""
357
return self.sorted_keys
361
def is_perturbating(self,order,model):
362
"""Returns wether this particle contributes in perturbation of the order passed
363
in argument given the model specified. It is very fast for usual models"""
365
for int in model['interactions'].get_type('base'):
366
# We discard the interactions with more than one type of orders
367
# contributing because it then doesn't necessarly mean that this
368
# particle (self) is charged under the group corresponding to the
369
# coupling order 'order'. The typical example is in SUSY which
370
# features a ' photon-gluon-squark-antisquark ' interaction which
371
# has coupling orders QED=1, QCD=1 and would induce the photon
372
# to be considered as a valid particle to circulate in a loop of
374
if len(int.get('orders'))>1:
376
if order in list(int.get('orders').keys()) and self.get('pdg_code') in \
377
[part.get('pdg_code') for part in int.get('particles')]:
382
def get_pdg_code(self):
383
"""Return the PDG code with a correct minus sign if the particle is its
386
if not self['is_part'] and not self['self_antipart']:
387
return - self['pdg_code']
389
return self['pdg_code']
391
def get_anti_pdg_code(self):
392
"""Return the PDG code of the antiparticle with a correct minus sign
393
if the particle is its own antiparticle"""
395
if not self['self_antipart']:
396
return - self.get_pdg_code()
398
return self['pdg_code']
401
"""Return the color code with a correct minus sign"""
403
if not self['is_part'] and abs(self['color']) in [3, 6]:
404
return - self['color']
408
def get_anti_color(self):
409
"""Return the color code of the antiparticle with a correct minus sign
412
if self['is_part'] and self['color'] not in [1, 8]:
413
return - self['color']
417
def get_charge(self):
418
"""Return the charge code with a correct minus sign"""
420
if not self['is_part']:
421
return - self['charge']
423
return self['charge']
425
def get_anti_charge(self):
426
"""Return the charge code of the antiparticle with a correct minus sign
430
return - self['charge']
432
return self['charge']
435
"""Return the name if particle, antiname if antiparticle"""
437
if not self['is_part'] and not self['self_antipart']:
438
return self['antiname']
442
def get_helicity_states(self, allow_reverse=True):
443
"""Return a list of the helicity states for the onshell particle"""
445
spin = self.get('spin')
452
elif spin == 3 and self.get('mass').lower() == 'zero':
458
elif spin == 4 and self.get('mass').lower() == 'zero':
464
elif spin == 5 and self.get('mass').lower() == 'zero':
467
elif spin in [5, 99]:
469
res = [-2, -1, 0, 1, 2]
471
raise self.PhysicsObjectError("No helicity state assignment for spin %d particles" % spin)
473
if allow_reverse and not self.get('is_part'):
479
def is_fermion(self):
480
"""Returns True if this is a fermion, False if boson"""
482
return self['spin'] % 2 == 0
485
"""Returns True if this is a boson, False if fermion"""
487
return self['spin'] % 2 == 1
489
#===============================================================================
491
#===============================================================================
492
class ParticleList(PhysicsObjectList):
493
"""A class to store lists of particles."""
495
def is_valid_element(self, obj):
496
"""Test if object obj is a valid Particle for the list."""
497
return isinstance(obj, Particle)
499
def get_copy(self, name):
500
"""Try to find a particle with the given name. Check both name
501
and antiname. If a match is found, return the a copy of the
502
corresponding particle (first one in the list), with the
503
is_part flag set accordingly. None otherwise."""
505
assert isinstance(name, str)
507
part = self.find_name(name)
509
# Then try to look for a particle with that PDG
516
if p.get_pdg_code()==pdg:
518
part.set('is_part', True)
520
elif p.get_anti_pdg_code()==pdg:
522
part.set('is_part', False)
526
part = copy.copy(part)
528
if part.get('name') == name:
529
part.set('is_part', True)
531
elif part.get('antiname') == name:
532
part.set('is_part', False)
536
def find_name(self, name):
537
"""Try to find a particle with the given name. Check both name
538
and antiname. If a match is found, return the a copy of the
539
corresponding particle (first one in the list), with the
540
is_part flag set accordingly. None otherwise."""
542
assert isinstance(name, str), "%s is not a valid string" % str(name)
545
if part.get('name') == name:
547
elif part.get('antiname') == name:
552
def generate_ref_dict(self):
553
"""Generate a dictionary of part/antipart pairs (as keys) and
559
ref_dict_to0[(part.get_pdg_code(), part.get_anti_pdg_code())] = [0]
560
ref_dict_to0[(part.get_anti_pdg_code(), part.get_pdg_code())] = [0]
564
def generate_dict(self):
565
"""Generate a dictionary from particle id to particle.
566
Include antiparticles.
571
for particle in self:
572
particle_dict[particle.get('pdg_code')] = particle
573
if not particle.get('self_antipart'):
574
antipart = copy.deepcopy(particle)
575
antipart.set('is_part', False)
576
particle_dict[antipart.get_pdg_code()] = antipart
581
#===============================================================================
583
#===============================================================================
584
class Interaction(PhysicsObject):
585
"""The interaction object containing the whole set of information
586
required to univocally characterize a given type of physical interaction:
588
particles: a list of particle ids
589
color: a list of string describing all the color structures involved
590
lorentz: a list of variable names describing all the Lorentz structure
592
couplings: dictionary listing coupling variable names. The key is a
593
2-tuple of integers referring to color and Lorentz structures
594
orders: dictionary listing order names (as keys) with their value
597
sorted_keys = ['id', 'particles', 'color', 'lorentz', 'couplings',
598
'orders','loop_particles','type','perturbation_type']
600
def default_setup(self):
601
"""Default values for all properties"""
604
self['particles'] = []
607
self['couplings'] = { (0, 0):'none'}
609
# The type of interactions can be 'base', 'UV' or 'R2'.
610
# For 'UV' or 'R2', one can always specify the loop it corresponds
611
# to by a tag in the second element of the list. If the tag is an
612
# empty list, then the R2/UV interaction will be recognized only
613
# based on the nature of the identity of the particles branching
614
# off the loop and the loop orders.
615
# Otherwise, the tag can be specified and it will be used when
616
# identifying the R2/UV interaction corresponding to a given loop
618
# The format is [(lp1ID,int1ID),(lp1ID,int1ID),(lp1ID,int1ID),etc...]
619
# Example of a tag for the following loop
621
# ___34_____ The ';' line is a gluon with ID 21
622
# 45/ ; The '|' line is a d-quark with ID 1
623
# ------< ; The numbers are the interactions ID
624
# \___;______ The tag for this loop would be:
625
# 12 ((21,34),(1,45),(1,12))
627
# This tag is equivalent to all its cyclic permutations. This is why
628
# it must be specified in the canonical order which is defined with
629
# by putting in front of the tag the lowest 2-tuple it contains.
630
# (the order relation is defined by comparing the particle ID first
631
# and the interaction ID after in case the particle ID are the same).
632
# In case there are two identical lowest 2-tuple in the tag, the
633
# tag chosen is such that it has the lowest second 2-tuple. The procedure
634
# is repeated again with the subsequent 2-tuple until there is only
635
# one cyclic permutation remaining and the ambiguity is resolved.
636
# This insures to have one unique unambiguous canonical tag chosen.
637
# In the example above, it would be:
638
# ((1,12),(21,34),(1,45))
639
# PS: Notice that in the UFO model, the tag-information is limited to
640
# the minimally relevant one which are the loop particles specified in
641
# in the attribute below. In this case, 'loop_particles' is the list of
642
# all the loops giving this same counterterm contribution.
643
# Each loop being represented by a set of the PDG of the particles
644
# (not repeated) constituting it. In the example above, it would simply
645
# be (1,21). In the UFO, if the loop particles are not specified then
646
# MG5 will account for this counterterm only once per concerned vertex.
647
# Taking the example of the three gluon vertex counterterm, one can
648
# possibly have in the ufo:
649
# VertexB = blabla, loop_particles = (b)
650
# VertexT = blabla, loop_particles = (t)
652
# VertexALL = blabla, loop_particles = ()
653
# In the first case UFO specifies the specific counterterm to the three-
654
# gluon loop with the bottom running in (VertexB) and with the top running
655
# in (VertexT). So MG5 will associate these counterterm vertices once to
656
# each of the two loop.
657
# In the case where UFO defined VertexALL, then whenever MG5 encounters
658
# a triangle three-gluon loop (say the bottom one), it will associate to
659
# it the vertex VertexALL but will not do so again when encountering the
660
# same loop with the top quark running in. This, because it assumes that
661
# the UFO vertexALL comprises all contributions already.
663
self['loop_particles']=[[]]
664
self['type'] = 'base'
665
self['perturbation_type'] = None
667
def filter(self, name, value):
668
"""Filter for valid interaction property values."""
671
#Should be an integer
672
if not isinstance(value, int):
673
raise self.PhysicsObjectError("%s is not a valid integer" % str(value))
675
if name == 'particles':
676
#Should be a list of valid particle names
677
if not isinstance(value, ParticleList):
678
raise self.PhysicsObjectError("%s is not a valid list of particles" % str(value))
680
if name == 'perturbation_type':
681
if value!=None and not isinstance(value, str):
682
raise self.PhysicsObjectError("%s is not a valid string" % str(value))
686
if not isinstance(value, str):
687
raise self.PhysicsObjectError("%s is not a valid string" % str(value))
688
if name == 'loop_particles':
689
if isinstance(value,list):
691
if isinstance(l,list):
693
if not isinstance(part,int):
694
raise self.PhysicsObjectError("%s is not a valid integer" % str(part))
696
raise self.PhysicsObjectError("%s is not a valid positive integer" % str(part))
699
#Should be a dict with valid order names ask keys and int as values
700
if not isinstance(value, dict):
701
raise self.PhysicsObjectError("%s is not a valid dict for coupling orders" % \
703
for order in value.keys():
704
if not isinstance(order, str):
705
raise self.PhysicsObjectError("%s is not a valid string" % str(order))
706
if not isinstance(value[order], int):
707
raise self.PhysicsObjectError("%s is not a valid integer" % str(value[order]))
709
if name in ['color']:
710
#Should be a list of list strings
711
if not isinstance(value, list):
712
raise self.PhysicsObjectError("%s is not a valid list of Color Strings" % str(value))
713
for mycolstring in value:
714
if not isinstance(mycolstring, color.ColorString):
715
raise self.PhysicsObjectError("%s is not a valid list of Color Strings" % str(value))
717
if name in ['lorentz']:
718
#Should be a list of list strings
719
if not isinstance(value, list):
720
raise self.PhysicsObjectError("%s is not a valid list of strings" % str(value))
722
if not isinstance(mystr, str):
723
raise self.PhysicsObjectError("%s is not a valid string" % str(mystr))
725
if name == 'couplings':
726
#Should be a dictionary of strings with (i,j) keys
727
if not isinstance(value, dict):
728
raise self.PhysicsObjectError("%s is not a valid dictionary for couplings" % \
731
for key in value.keys():
732
if not isinstance(key, tuple):
733
raise self.PhysicsObjectError("%s is not a valid tuple" % str(key))
735
raise self.PhysicsObjectError("%s is not a valid tuple with 2 elements" % str(key))
736
if not isinstance(key[0], int) or not isinstance(key[1], int):
737
raise self.PhysicsObjectError("%s is not a valid tuple of integer" % str(key))
738
if not isinstance(value[key], str):
739
raise self.PhysicsObjectError("%s is not a valid string" % value[key])
743
def get_sorted_keys(self):
744
"""Return particle property names as a nicely sorted list."""
746
return self.sorted_keys
748
def is_perturbating(self, orders_considered):
749
""" Returns if this interaction comes from the perturbation of one of
750
the order listed in the argument """
752
if self['perturbation_type']==None:
755
return (self['perturbation_type'] in orders_considered)
758
""" Returns if the interaction is of R2 type."""
760
# Precaution only useful because some tests have a predefined model
761
# bypassing the default_setup and for which type was not defined.
762
if 'type' in list(self.keys()):
763
return (len(self['type'])>=2 and self['type'][:2]=='R2')
768
""" Returns if the interaction is of UV type."""
770
# Precaution only useful because some tests have a predefined model
771
# bypassing the default_setup and for which type was not defined.
772
if 'type' in list(self.keys()):
773
return (len(self['type'])>=2 and self['type'][:2]=='UV')
778
""" Returns if the interaction is of UVmass type."""
780
# Precaution only useful because some tests have a predefined model
781
# bypassing the default_setup and for which type was not defined.
782
if 'type' in list(self.keys()):
783
return (len(self['type'])>=6 and self['type'][:6]=='UVmass')
788
""" Returns if the interaction is of UVmass type."""
790
# Precaution only useful because some tests have a predefined model
791
# bypassing the default_setup and for which type was not defined.
792
if 'type' in list(self.keys()):
793
return (len(self['type'])>=6 and self['type'][:6]=='UVloop')
798
""" Returns if the interaction is of UVmass type."""
800
# Precaution only useful because some tests have a predefined model
801
# bypassing the default_setup and for which type was not defined.
802
if 'type' in list(self.keys()):
803
return (len(self['type'])>=6 and self['type'][:6]=='UVtree')
808
""" Returns if the interaction is of the UVCT type which means that
809
it has been selected as a possible UV counterterm interaction for this
810
process. Such interactions are marked by having the 'UVCT_SPECIAL' order
811
key in their orders."""
813
# Precaution only useful because some tests have a predefined model
814
# bypassing the default_setup and for which type was not defined.
815
if 'UVCT_SPECIAL' in list(self['orders'].keys()):
820
def get_epsilon_order(self):
821
""" Returns 0 if this interaction contributes to the finite part of the
822
amplitude and 1 (2) is it contributes to its single (double) pole """
824
if 'type' in list(self.keys()):
825
if '1eps' in self['type']:
827
elif '2eps' in self['type']:
834
def generate_dict_entries(self, ref_dict_to0, ref_dict_to1):
835
"""Add entries corresponding to the current interactions to
836
the reference dictionaries (for n>0 and n-1>1)"""
838
# Create n>0 entries. Format is (p1,p2,p3,...):interaction_id.
839
# We are interested in the unordered list, so use sorted()
841
pdg_tuple = tuple(sorted([p.get_pdg_code() for p in self['particles']]))
842
if pdg_tuple not in list(ref_dict_to0.keys()):
843
ref_dict_to0[pdg_tuple] = [self['id']]
845
ref_dict_to0[pdg_tuple].append(self['id'])
847
# Create n-1>1 entries. Note that, in the n-1 > 1 dictionary,
848
# the n-1 entries should have opposite sign as compared to
849
# interaction, since the interaction has outgoing particles,
850
# while in the dictionary we treat the n-1 particles as
853
for part in self['particles']:
855
# We are interested in the unordered list, so use sorted()
856
pdg_tuple = tuple(sorted([p.get_pdg_code() for (i, p) in \
857
enumerate(self['particles']) if \
858
i != self['particles'].index(part)]))
859
pdg_part = part.get_anti_pdg_code()
860
if pdg_tuple in list(ref_dict_to1.keys()):
861
if (pdg_part, self['id']) not in ref_dict_to1[pdg_tuple]:
862
ref_dict_to1[pdg_tuple].append((pdg_part, self['id']))
864
ref_dict_to1[pdg_tuple] = [(pdg_part, self['id'])]
866
def get_WEIGHTED_order(self, model):
867
"""Get the WEIGHTED order for this interaction, for equivalent
868
3-particle vertex. Note that it can be fractional."""
870
return float(sum([model.get('order_hierarchy')[key]*self.get('orders')[key]\
871
for key in self.get('orders')]))/ \
872
max((len(self.get('particles'))-2), 1)
875
"""String representation of an interaction. Outputs valid Python
876
with improved format. Overrides the PhysicsObject __str__ to only
877
display PDG code of involved particles."""
881
for prop in self.get_sorted_keys():
882
if isinstance(self[prop], str):
883
mystr = mystr + ' \'' + prop + '\': \'' + \
885
elif isinstance(self[prop], float):
886
mystr = mystr + ' \'' + prop + '\': %.2f,\n' % self[prop]
887
elif isinstance(self[prop], ParticleList):
888
mystr = mystr + ' \'' + prop + '\': [%s],\n' % \
889
','.join([str(part.get_pdg_code()) for part in self[prop]])
891
mystr = mystr + ' \'' + prop + '\': ' + \
892
repr(self[prop]) + ',\n'
893
mystr = mystr.rstrip(',\n')
894
mystr = mystr + '\n}'
898
def canonical_repr(self):
899
""" Returns a string representation that allows to order CT vertices in
900
a diagram so that identical HelasMatrix element (i.e. typically different
901
flavors with identical masses and couplings) can be matched even though
902
the order of the couplings specified in the UFO is different. """
904
return '%s|%s'%(self['type'],
905
'&'.join( sorted('%s_%s_%s'%(self['color'][k[0]],self['lorentz'][k[1]],v)
906
for k,v in self['couplings'].items() ) ) )
908
def get_canonical_couplings_keys_order(self):
909
""" Returns a list of the keys to the 'couplings' dictionary, canonically
910
ordered so so that identical HelasMatrix element (i.e. typically different
911
flavors with identical masses and couplings) can be matched even though
912
the order of the couplings specified in the UFO is different. """
914
return sorted(list(self['couplings'].keys()), key=lambda k:
915
'%s_%s_%s'%(self['color'][k[0]],self['lorentz'][k[1]],self['couplings'][k]))
918
#===============================================================================
920
#===============================================================================
921
class InteractionList(PhysicsObjectList):
922
"""A class to store lists of interactionss."""
924
def is_valid_element(self, obj):
925
"""Test if object obj is a valid Interaction for the list."""
927
return isinstance(obj, Interaction)
929
def generate_ref_dict(self,useR2UV=False, useUVCT=False):
930
"""Generate the reference dictionaries from interaction list.
931
Return a list where the first element is the n>0 dictionary and
932
the second one is n-1>1."""
939
if useR2UV or (not inter.is_UV() and not inter.is_R2() and \
940
not inter.is_UVCT()):
941
inter.generate_dict_entries(ref_dict_to0, ref_dict_to1)
942
if useUVCT and inter.is_UVCT():
943
inter.generate_dict_entries(ref_dict_to0, ref_dict_to1)
945
return [ref_dict_to0, ref_dict_to1]
947
def generate_dict(self):
948
"""Generate a dictionary from interaction id to interaction.
951
interaction_dict = {}
954
interaction_dict[inter.get('id')] = inter
956
return interaction_dict
958
def synchronize_interactions_with_particles(self, particle_dict):
959
"""Make sure that the particles in the interactions are those
960
in the particle_dict, and that there are no interactions
961
refering to particles that don't exist. To be called when the
962
particle_dict is updated in a model.
966
while iint < len(self):
968
particles = inter.get('particles')
970
for ipart, part in enumerate(particles):
971
particles[ipart] = particle_dict[part.get_pdg_code()]
974
# This interaction has particles that no longer exist
977
def get_type(self, type):
978
""" return all interactions in the list of type 'type' """
979
return InteractionList([int for int in self if int.get('type')==type])
982
""" return all interactions in the list of type R2 """
983
return InteractionList([int for int in self if int.is_R2()])
986
""" return all interactions in the list of type UV """
987
return InteractionList([int for int in self if int.is_UV()])
989
def get_UVmass(self):
990
""" return all interactions in the list of type UVmass """
991
return InteractionList([int for int in self if int.is_UVmass()])
993
def get_UVtree(self):
994
""" return all interactions in the list of type UVtree """
995
return InteractionList([int for int in self if int.is_UVtree()])
997
def get_UVloop(self):
998
""" return all interactions in the list of type UVloop """
999
return InteractionList([int for int in self if int.is_UVloop()])
1001
#===============================================================================
1003
#===============================================================================
1004
class Model(PhysicsObject):
1005
"""A class to store all the model information."""
1007
mg5_name = False #store if particle name follow mg5 convention
1009
def __init__(self, init_dict={}):
1010
"""Creates a new particle object. If a dictionary is given, tries to
1011
use it to give values to properties."""
1014
self.default_setup()
1016
assert isinstance(init_dict, dict), \
1017
"Argument %s is not a dictionary" % repr(init_dict)
1020
for item in init_dict.keys():
1021
self[item] = init_dict[item]
1023
def default_setup(self):
1026
self['particles'] = ParticleList()
1027
self['interactions'] = InteractionList()
1028
self['parameters'] = None
1029
self['functions'] = None
1030
self['couplings'] = None
1031
self['lorentz'] = None
1032
self['particle_dict'] = {}
1033
self['interaction_dict'] = {}
1034
self['ref_dict_to0'] = {}
1035
self['ref_dict_to1'] = {}
1036
self['got_majoranas'] = None
1037
self['order_hierarchy'] = {}
1038
self['conserved_charge'] = set()
1039
self['coupling_orders'] = None
1040
self['expansion_order'] = None
1041
self['version_tag'] = None # position of the directory (for security)
1042
self['gauge'] = [0, 1]
1043
self['case_sensitive'] = True
1044
self['allow_pickle'] = True
1045
self['limitations'] = [] # MLM means that the model can sometimes have issue with MLM/default scale.
1046
# attribute which might be define if needed
1047
#self['name2pdg'] = {'name': pdg}
1051
def filter(self, name, value):
1052
"""Filter for model property values"""
1054
if name in ['name']:
1055
if not isinstance(value, str):
1056
raise self.PhysicsObjectError("Object of type %s is not a string" %type(value))
1058
elif name == 'particles':
1059
if not isinstance(value, ParticleList):
1060
raise self.PhysicsObjectError("Object of type %s is not a ParticleList object" % \
1062
elif name == 'interactions':
1063
if not isinstance(value, InteractionList):
1064
raise self.PhysicsObjectError("Object of type %s is not a InteractionList object" % \
1066
elif name == 'particle_dict':
1067
if not isinstance(value, dict):
1068
raise self.PhysicsObjectError("Object of type %s is not a dictionary" % \
1070
elif name == 'interaction_dict':
1071
if not isinstance(value, dict):
1072
raise self.PhysicsObjectError("Object of type %s is not a dictionary" % type(value))
1074
elif name == 'ref_dict_to0':
1075
if not isinstance(value, dict):
1076
raise self.PhysicsObjectError("Object of type %s is not a dictionary" % type(value))
1078
elif name == 'ref_dict_to1':
1079
if not isinstance(value, dict):
1080
raise self.PhysicsObjectError("Object of type %s is not a dictionary" % type(value))
1082
elif name == 'got_majoranas':
1083
if not (isinstance(value, bool) or value == None):
1084
raise self.PhysicsObjectError("Object of type %s is not a boolean" % type(value))
1086
elif name == 'conserved_charge':
1087
if not (isinstance(value, set)):
1088
raise self.PhysicsObjectError("Object of type %s is not a set" % type(value))
1090
elif name == 'version_tag':
1091
if not (isinstance(value, str)):
1092
raise self.PhysicsObjectError("Object of type %s is not a string" % type(value))
1094
elif name == 'order_hierarchy':
1095
if not isinstance(value, dict):
1096
raise self.PhysicsObjectError("Object of type %s is not a dictionary" % \
1098
for key in value.keys():
1099
if not isinstance(value[key],int):
1100
raise self.PhysicsObjectError("Object of type %s is not an integer" % \
1102
elif name == 'gauge':
1103
if not (isinstance(value, list)):
1104
raise self.PhysicsObjectError("Object of type %s is not a list" % type(value))
1106
elif name == 'case_sensitive':
1107
if not value in [True ,False]:
1108
raise self.PhysicsObjectError("Object of type %s is not a boolean" % type(value))
1113
def get(self, name):
1114
"""Get the value of the property name."""
1116
if (name == 'ref_dict_to0' or name == 'ref_dict_to1') and \
1118
if self['interactions']:
1119
[self['ref_dict_to0'], self['ref_dict_to1']] = \
1120
self['interactions'].generate_ref_dict()
1121
self['ref_dict_to0'].update(
1122
self['particles'].generate_ref_dict())
1124
if (name == 'particle_dict') and not self[name]:
1125
if self['particles']:
1126
self['particle_dict'] = self['particles'].generate_dict()
1127
if self['interactions']:
1128
self['interactions'].synchronize_interactions_with_particles(\
1129
self['particle_dict'])
1130
if name == 'modelpath':
1131
modeldir = self.get('version_tag').rsplit('##',1)[0]
1132
if os.path.exists(modeldir):
1133
modeldir = os.path.expanduser(modeldir)
1136
raise Exception("path %s not valid anymore." % modeldir)
1137
#modeldir = os.path.join(os.path.dirname(modeldir),
1138
# os.path.basename(modeldir).rsplit("-",1)[0])
1139
#if os.path.exists(modeldir):
1141
#raise Exception, 'Invalid Path information: %s' % self.get('version_tag')
1142
elif name == 'modelpath+restriction':
1143
modeldir = self.get('version_tag').rsplit('##',1)[0]
1144
modelname = self['name']
1145
if not os.path.exists(modeldir):
1146
raise Exception("path %s not valid anymore" % modeldir)
1147
modeldir = os.path.dirname(modeldir)
1148
modeldir = pjoin(modeldir, modelname)
1149
modeldir = os.path.expanduser(modeldir)
1151
elif name == 'restrict_name':
1152
modeldir = self.get('version_tag').rsplit('##',1)[0]
1153
modelname = self['name']
1154
basename = os.path.basename(modeldir)
1155
restriction = modelname[len(basename)+1:]
1158
if (name == 'interaction_dict') and not self[name]:
1159
if self['interactions']:
1160
self['interaction_dict'] = self['interactions'].generate_dict()
1162
if (name == 'got_majoranas') and self[name] == None:
1163
if self['particles']:
1164
self['got_majoranas'] = self.check_majoranas()
1166
if (name == 'coupling_orders') and self[name] == None:
1167
if self['interactions']:
1168
self['coupling_orders'] = self.get_coupling_orders()
1170
if (name == 'order_hierarchy') and not self[name]:
1171
if self['interactions']:
1172
self['order_hierarchy'] = self.get_order_hierarchy()
1174
if (name == 'expansion_order') and self[name] == None:
1175
if self['interactions']:
1176
self['expansion_order'] = \
1177
dict([(order, -1) for order in self.get('coupling_orders')])
1179
if (name == 'name2pdg') and 'name2pdg' not in self:
1180
self['name2pdg'] = {}
1181
for p in self.get('particles'):
1182
self['name2pdg'][p.get('antiname')] = -1*p.get('pdg_code')
1183
self['name2pdg'][p.get('name')] = p.get('pdg_code')
1185
return Model.__bases__[0].get(self, name) # call the mother routine
1187
def set(self, name, value, force = False):
1188
"""Special set for particles and interactions - need to
1189
regenerate dictionaries."""
1191
if name == 'particles':
1192
# Ensure no doublets in particle list
1194
# Reset dictionaries
1195
self['particle_dict'] = {}
1196
self['ref_dict_to0'] = {}
1197
self['got_majoranas'] = None
1199
if name == 'interactions':
1200
# Ensure no doublets in interaction list
1202
# Reset dictionaries
1203
self['interaction_dict'] = {}
1204
self['ref_dict_to1'] = {}
1205
self['ref_dict_to0'] = {}
1206
self['got_majoranas'] = None
1207
self['coupling_orders'] = None
1208
self['order_hierarchy'] = {}
1209
self['expansion_order'] = None
1211
if name == 'name2pdg':
1212
self['name2pgg'] = value
1215
result = Model.__bases__[0].set(self, name, value, force) # call the mother routine
1217
if name == 'particles':
1218
# Recreate particle_dict
1219
self.get('particle_dict')
1223
def actualize_dictionaries(self):
1224
"""This function actualizes the dictionaries"""
1226
[self['ref_dict_to0'], self['ref_dict_to1']] = \
1227
self['interactions'].generate_ref_dict()
1228
self['ref_dict_to0'].update(
1229
self['particles'].generate_ref_dict())
1231
def get_sorted_keys(self):
1232
"""Return process property names as a nicely sorted list."""
1234
return ['name', 'particles', 'parameters', 'interactions',
1235
'couplings','lorentz', 'gauge']
1237
def get_particle(self, id):
1238
"""Return the particle corresponding to the id / name"""
1241
return self["particle_dict"][id]
1243
if isinstance(id, int):
1245
return self.get("particle_dict")[id]
1246
except Exception as error:
1249
if not hasattr(self, 'name2part'):
1250
self.create_name2part()
1252
return self.name2part[id]
1256
def create_name2part(self):
1257
"""create a dictionary name 2 part"""
1260
for part in self.get("particle_dict").values():
1261
self.name2part[part.get('name')] = part
1262
self.name2part[part.get('antiname')] = part
1264
def get_lorentz(self, name):
1265
"""return the lorentz object from the associate name"""
1266
if hasattr(self, 'lorentz_name2obj'):
1267
return self.lorentz_name2obj[name]
1269
self.create_lorentz_dict()
1270
return self.lorentz_name2obj[name]
1272
def create_lorentz_dict(self):
1273
"""create the dictionary linked to the lorentz structure"""
1274
self.lorentz_name2obj = {}
1275
self.lorentz_expr2name = {}
1276
if not self.get('lorentz'):
1278
for lor in self.get('lorentz'):
1279
self.lorentz_name2obj[lor.name] = lor
1280
self.lorentz_expr2name[lor.structure] = lor.name
1282
def get_interaction(self, id):
1283
"""Return the interaction corresponding to the id"""
1286
return self.get("interaction_dict")[id]
1290
def get_parameter(self, name):
1291
"""Return the parameter associated to the name NAME"""
1293
# If information is saved
1294
if hasattr(self, 'parameters_dict') and self.parameters_dict:
1296
return self.parameters_dict[name]
1298
# try to reload it before crashing
1301
# Else first build the dictionary
1302
self.parameters_dict = {}
1303
for data in self['parameters'].values():
1304
[self.parameters_dict.__setitem__(p.name,p) for p in data]
1306
return self.parameters_dict[name]
1308
def get_coupling_orders(self):
1309
"""Determine the coupling orders of the model"""
1310
return set(sum([list(i.get('orders').keys()) for i in \
1311
self.get('interactions')], []))
1313
def get_order_hierarchy(self):
1314
"""Set a default order hierarchy for the model if not set by the UFO."""
1315
# Set coupling hierachy
1316
hierarchy = dict([(order, 1) for order in self.get('coupling_orders')])
1317
# Special case for only QCD and QED couplings, unless already set
1318
if self.get('coupling_orders') == set(['QCD', 'QED']):
1319
hierarchy['QED'] = 2
1323
def get_nflav(self):
1324
"""returns the number of light quark flavours in the model."""
1325
return len([p for p in self.get('particles') \
1326
if p['spin'] == 2 and p['is_part'] and \
1327
p ['color'] != 1 and p['mass'].lower() == 'zero'])
1330
def get_quark_pdgs(self):
1331
"""returns the PDG codes of the light quarks and antiquarks"""
1332
pdg_list = [p['pdg_code'] for p in self.get('particles') \
1333
if p['spin'] == 2 and \
1334
p['color'] == 3 and \
1335
p['charge'] != 0. and p['mass'].lower() == 'zero']
1337
for p in pdg_list[:]:
1338
if not self.get('particle_dict')[p]['self_antipart']:
1339
pdg_list.append(self.get('particle_dict')[p].get_anti_pdg_code())
1341
return sorted(pdg_list)
1344
def get_nleps(self):
1345
"""returns the number of light lepton flavours in the model."""
1346
return len([p for p in self.get('particles') \
1347
if p['spin'] == 2 and p['is_part'] and \
1348
p['color'] == 1 and \
1349
p['charge'] != 0. and p['mass'].lower() == 'zero'])
1352
def get_lepton_pdgs(self):
1353
"""returns the PDG codes of the light leptons and antileptons"""
1354
pdg_list = [p['pdg_code'] for p in self.get('particles') \
1355
if p['spin'] == 2 and \
1356
p['color'] == 1 and \
1357
p['charge'] != 0. and p['mass'].lower() == 'zero']
1359
for p in pdg_list[:]:
1360
if not self.get('particle_dict')[p]['self_antipart']:
1361
pdg_list.append(self.get('particle_dict')[p].get_anti_pdg_code())
1363
return sorted(pdg_list)
1366
def get_particles_hierarchy(self):
1367
"""Returns the order hierarchies of the model and the
1368
particles which have interactions in at least this hierarchy
1369
(used in find_optimal_process_orders in MultiProcess diagram
1372
Check the coupling hierarchy of the model. Assign all
1373
particles to the different coupling hierarchies so that a
1374
particle is considered to be in the highest hierarchy (i.e.,
1375
with lowest value) where it has an interaction.
1378
# Find coupling orders in model
1379
coupling_orders = self.get('coupling_orders')
1380
# Loop through the different coupling hierarchy values, so we
1381
# start with the most dominant and proceed to the least dominant
1382
hierarchy = sorted(list(set([self.get('order_hierarchy')[k] for \
1383
k in coupling_orders])))
1385
# orders is a rising list of the lists of orders with a given hierarchy
1387
for value in hierarchy:
1388
orders.append([ k for (k, v) in \
1389
self.get('order_hierarchy').items() if \
1392
# Extract the interaction that correspond to the different
1393
# coupling hierarchies, and the corresponding particles
1396
for iorder, order in enumerate(orders):
1397
sum_orders = sum(orders[:iorder+1], [])
1398
sum_interactions = sum(interactions[:iorder], [])
1399
sum_particles = sum([list(p) for p in particles[:iorder]], [])
1400
# Append all interactions that have only orders with at least
1402
interactions.append([i for i in self.get('interactions') if \
1403
not i in sum_interactions and \
1404
not any([k not in sum_orders for k in \
1405
i.get('orders').keys()])])
1406
# Append the corresponding particles, excluding the
1407
# particles that have already been added
1408
particles.append(set(sum([[p.get_pdg_code() for p in \
1409
inter.get('particles') if \
1410
p.get_pdg_code() not in sum_particles] \
1411
for inter in interactions[-1]], [])))
1413
return particles, hierarchy
1415
def get_max_WEIGHTED(self):
1416
"""Return the maximum WEIGHTED order for any interaction in the model,
1417
for equivalent 3-particle vertices. Note that it can be fractional."""
1419
return max([inter.get_WEIGHTED_order(self) for inter in \
1420
self.get('interactions')])
1423
def check_majoranas(self):
1424
"""Return True if there is fermion flow violation, False otherwise"""
1426
if any([part.is_fermion() and part.get('self_antipart') \
1427
for part in self.get('particles')]):
1430
# No Majorana particles, but may still be fermion flow
1431
# violating interactions
1432
for inter in self.get('interactions'):
1433
# Do not look at UV Wfct renormalization counterterms
1434
if len(inter.get('particles'))==1:
1436
fermions = [p for p in inter.get('particles') if p.is_fermion()]
1437
for i in range(0, len(fermions), 2):
1438
if fermions[i].get('is_part') == \
1439
fermions[i+1].get('is_part'):
1440
# This is a fermion flow violating interaction
1442
# No fermion flow violations
1445
def reset_dictionaries(self):
1446
"""Reset all dictionaries and got_majoranas. This is necessary
1447
whenever the particle or interaction content has changed. If
1448
particles or interactions are set using the set routine, this
1449
is done automatically."""
1451
self['particle_dict'] = {}
1452
self['ref_dict_to0'] = {}
1453
self['got_majoranas'] = None
1454
self['interaction_dict'] = {}
1455
self['ref_dict_to1'] = {}
1456
self['ref_dict_to0'] = {}
1458
def pass_particles_name_in_mg_default(self):
1459
"""Change the name of the particles such that all SM and MSSM particles
1460
follows the MG convention"""
1462
self.mg5_name = True
1464
# Check that default name/antiname is not already use
1465
def check_name_free(self, name):
1466
""" check if name is not use for a particle in the model if it is
1467
raise an MadGraph5error"""
1468
part = self['particles'].find_name(name)
1471
'%s particles with pdg code %s is in conflict with MG ' + \
1472
'convention name for particle %s.\n Use -modelname in order ' + \
1473
'to use the particles name defined in the model and not the ' + \
1474
'MadGraph5_aMC@NLO convention'
1476
raise MadGraph5Error(error_text % \
1477
(part.get_name(), part.get_pdg_code(), pdg))
1479
default = self.load_default_name()
1481
for pdg in default.keys():
1482
part = self.get_particle(pdg)
1485
antipart = self.get_particle(-pdg)
1486
name = part.get_name()
1487
if name != default[pdg]:
1488
check_name_free(self, default[pdg])
1489
if part.get('is_part'):
1490
part.set('name', default[pdg])
1492
antipart.set('name', default[pdg])
1494
part.set('antiname', default[pdg])
1496
part.set('antiname', default[pdg])
1498
antipart.set('antiname', default[pdg])
1500
#additional check for the Higgs in the mssm
1501
if self.get('name') == 'mssm' or self.get('name').startswith('mssm-'):
1502
part = self.get_particle(25)
1503
part.set('name', 'h1')
1504
part.set('antiname', 'h1')
1508
def change_parameter_name_with_prefix(self, prefix='mdl_'):
1509
""" Change all model parameter by a given prefix.
1510
Modify the parameter if some of them are identical up to the case"""
1514
keys = list(self.get('parameters').keys())
1516
for param in self['parameters'][key]:
1517
lower_name = param.name.lower()
1521
lower_dict[lower_name].append(param)
1523
lower_dict[lower_name] = [param]
1525
duplicate.add(lower_name)
1526
logger.debug('%s is defined both as lower case and upper case.'
1529
if prefix == '' and not duplicate:
1532
re_expr = r'''\b(%s)\b'''
1535
# recast all parameter in prefix_XX
1537
for param in self['parameters'][key]:
1538
value = param.name.lower()
1539
if value in ['as','mu_r', 'zero','aewm1','g']:
1541
elif value.startswith(prefix):
1543
elif value in duplicate:
1544
continue # handle later
1546
change[param.name] = '%s%s' % (prefix,param.name)
1547
to_change.append(param.name)
1548
param.name = change[param.name]
1550
for value in duplicate:
1551
for i, var in enumerate(lower_dict[value]):
1552
to_change.append(var.name)
1553
new_name = '%s%s%s' % (prefix, var.name.lower(),
1554
('__%d'%(i+1) if i>0 else ''))
1555
change[var.name] = new_name
1557
to_change.append(var.name)
1558
assert 'zero' not in to_change
1559
replace = lambda match_pattern: change[match_pattern.groups()[0]]
1564
if 'parameter_dict' in self:
1565
new_dict = dict( (change[name] if (name in change) else name, value) for
1566
name, value in self['parameter_dict'].items())
1567
self['parameter_dict'] = new_dict
1569
if hasattr(self,'map_CTcoup_CTparam'):
1570
# If the map for the dependence of couplings to CTParameters has
1571
# been defined, we must apply the renaming there as well.
1573
self.map_CTcoup_CTparam = dict( (coup_name,
1574
[change[name] if (name in change) else name for name in params])
1575
for coup_name, params in self.map_CTcoup_CTparam.items() )
1578
while i*1000 <= len(to_change):
1579
one_change = to_change[i*1000: min((i+1)*1000,len(to_change))]
1581
rep_pattern = re.compile('\\b%s\\b'% (re_expr % ('\\b|\\b'.join(one_change))))
1585
if key == ('external',):
1587
for param in self['parameters'][key]:
1588
param.expr = rep_pattern.sub(replace, param.expr)
1590
for key in self['couplings'].keys():
1591
for coup in self['couplings'][key]:
1592
coup.expr = rep_pattern.sub(replace, coup.expr)
1594
# change form-factor
1595
ff = [l.formfactors for l in self['lorentz'] if hasattr(l, 'formfactors')]
1596
ff = set(sum(ff,[])) # here we have the list of ff used in the model
1598
f.value = rep_pattern.sub(replace, f.value)
1601
for part in self['particles']:
1602
if str(part.get('mass')) in one_change:
1603
part.set('mass', rep_pattern.sub(replace, str(part.get('mass'))))
1604
if str(part.get('width')) in one_change:
1605
part.set('width', rep_pattern.sub(replace, str(part.get('width'))))
1606
if hasattr(part, 'partial_widths'):
1607
for key, value in part.partial_widths.items():
1608
part.partial_widths[key] = rep_pattern.sub(replace, value)
1610
#ensure that the particle_dict is up-to-date
1611
self['particle_dict'] =''
1612
self.get('particle_dict')
1616
def get_first_non_pdg(self):
1617
"""Return the first positive number that is not a valid PDG code"""
1618
return [c for c in range(1, len(self.get('particles')) + 1) if \
1619
c not in list(self.get('particle_dict').keys())][0]
1622
def write_param_card(self, filepath=None):
1623
"""Write out the param_card, and return as string."""
1625
import models.write_param_card as writer
1627
out = StringIO.StringIO() # it's suppose to be written in a file
1630
param = writer.ParamCardWriter(self, filepath=out)
1632
return out.getvalue()
1637
def load_default_name():
1638
""" load the default for name convention """
1640
logger.info('Change particles name to pass to MG5 convention')
1642
for line in open(os.path.join(MG5DIR, 'input', \
1643
'particles_name_default.txt')):
1644
line = line.lstrip()
1645
if line.startswith('#'):
1650
logger.warning('Invalid syntax in interface/default_name:\n %s' % line)
1652
default[int(args[0])] = args[1].lower()
1656
def change_electroweak_mode(self, mode):
1657
"""Change the electroweak mode. The only valid mode now is external.
1658
Where in top of the default MW and sw2 are external parameters."""
1660
if isinstance(mode, str) and "_" in mode:
1661
mode = set([s.lower() for s in mode.split('_')])
1663
assert mode in ["external",set(['mz','mw','alpha'])]
1666
W = self.get('particle_dict')[24]
1668
raise InvalidCmd('No W particle in the model impossible to '+
1669
'change the EW scheme!')
1671
if mode=='external':
1672
MW = self.get_parameter(W.get('mass'))
1673
if not isinstance(MW, ParamCardVariable):
1674
newMW = ParamCardVariable(MW.name, MW.value, 'MASS', [24])
1676
newMW.value = 80.385
1677
#remove the old definition
1678
self.get('parameters')[MW.depend].remove(MW)
1680
self.add_param(newMW, ['external'])
1682
# Now check for sw2. if not define bypass this
1684
sw2 = self.get_parameter('sw2')
1687
sw2 = self.get_parameter('mdl_sw2')
1692
newsw2 = ParamCardVariable(sw2.name,sw2.value, 'SMINPUTS', [4])
1693
if not newsw2.value:
1694
newsw2.value = 0.222246485786
1695
#remove the old definition
1696
self.get('parameters')[sw2.depend].remove(sw2)
1698
self.add_param(newsw2, ['external'])
1699
# Force a refresh of the parameter dictionary
1700
self.parameters_dict = None
1703
elif mode==set(['mz','mw','alpha']):
1704
# For now, all we support is to go from mz, Gf, alpha to mz, mw, alpha
1705
W = self.get('particle_dict')[24]
1706
mass = self.get_parameter(W.get('mass'))
1707
mass_expr = 'cmath.sqrt(%(prefix)sMZ__exp__2/2. + cmath.sqrt('+\
1708
'%(prefix)sMZ__exp__4/4. - (%(prefix)saEW*cmath.pi*%(prefix)s'+\
1709
'MZ__exp__2)/(%(prefix)sGf*%(prefix)ssqrt__2)))'
1710
if 'external' in mass.depend:
1711
# Nothing to be done
1714
if mass.expr == mass_expr%{'prefix':''}:
1717
elif mass.expr == mass_expr%{'prefix':'mdl_'}:
1721
MW = ParamCardVariable(mass.name, mass.value, 'MASS', [24])
1724
self.get('parameters')[('external',)].append(MW)
1725
self.get('parameters')[mass.depend].remove(mass)
1726
# Make Gf an internal parameter
1727
new_param = ModelVariable('Gf',
1728
'-%(prefix)saEW*%(prefix)sMZ**2*cmath.pi/(cmath.sqrt(2)*%(MW)s**2*(%(MW)s**2 - %(prefix)sMZ**2))' %\
1729
{'MW': mass.name,'prefix':prefix}, 'complex', mass.depend)
1730
Gf = self.get_parameter('%sGf'%prefix)
1731
self.get('parameters')[('external',)].remove(Gf)
1732
self.add_param(new_param, ['%saEW'%prefix])
1733
# Force a refresh of the parameter dictionary
1734
self.parameters_dict = None
1739
def change_mass_to_complex_scheme(self, toCMS=True):
1740
"""modify the expression changing the mass to complex mass scheme"""
1742
# 1) Change the 'CMSParam' of loop_qcd_qed model to 1.0 so as to remove
1743
# the 'real' prefix fromall UVCT wf renormalization expressions.
1744
# If toCMS is False, then it makes sure CMSParam is 0.0 and returns
1746
# 2) Find All input parameter mass and width associated
1747
# Add a internal parameter and replace mass with that param
1748
# 3) Find All mass fixed by the model and width associated
1749
# -> Both need to be fixed with a real() /Imag()
1750
# 4) Find All width set by the model
1751
# -> Need to be set with a real()
1752
# 5) Fix the Yukawa mass to the value of the complex mass/ real mass
1753
# 6) Loop through all expression and modify those accordingly
1754
# Including all parameter expression as complex
1757
CMSParam = self.get_parameter('CMSParam')
1760
CMSParam = self.get_parameter('mdl_CMSParam')
1764
# Handle the case where we want to make sure the CMS is turned off
1767
CMSParam.expr = '0.0'
1770
# Now handle the case where we want to turn to CMS.
1772
CMSParam.expr = '1.0'
1775
mass_widths = [] # parameter which should stay real
1776
for particle in self.get('particles'):
1777
m = particle.get('width')
1778
if m in mass_widths:
1780
mass_widths.append(particle.get('width'))
1781
mass_widths.append(particle.get('mass'))
1782
width = self.get_parameter(particle.get('width'))
1783
if (isinstance(width.value, (complex,float)) and abs(width.value)==0.0) or \
1784
width.name.lower() =='zero':
1785
#everything is fine since the width is zero
1787
if not isinstance(width, ParamCardVariable):
1788
width.expr = 're(%s)' % width.expr
1789
mass = self.get_parameter(particle.get('mass'))
1790
if (isinstance(width.value, (complex,float)) and abs(width.value)!=0.0) or \
1791
mass.name.lower() != 'zero':
1792
# special SM treatment to change the gauge scheme automatically.
1793
if particle.get('pdg_code') == 24 and isinstance(mass,
1795
status = self.change_electroweak_mode(
1796
set(['mz','mw','alpha']))
1797
# Use the newly defined parameter for the W mass
1798
mass = self.get_parameter(particle.get('mass'))
1800
logger.warning('The W mass is not an external '+
1801
'parameter in this model and the automatic change of'+
1802
' electroweak scheme changed. This is not advised for '+
1803
'applying the complex mass scheme.')
1805
# Add A new parameter CMASS
1806
#first compute the dependencies (as,...)
1807
depend = list(set(mass.depend + width.depend))
1808
if len(depend)>1 and 'external' in depend:
1809
depend.remove('external')
1810
depend = tuple(depend)
1811
if depend == ('external',):
1814
# Create the new parameter
1815
if isinstance(mass, ParamCardVariable):
1816
New_param = ModelVariable('CMASS_'+mass.name,
1817
'cmath.sqrt(%(mass)s**2 - complex(0,1) * %(mass)s * %(width)s)' \
1818
% {'mass': mass.name, 'width': width.name},
1821
New_param = ModelVariable('CMASS_'+mass.name,
1822
mass.expr, 'complex', depend)
1823
# Modify the treatment of the width in this case
1824
if not isinstance(width, ParamCardVariable):
1825
width.expr = '- im(%s**2) / cmath.sqrt(re(%s**2))' % (mass.expr, mass.expr)
1827
# Remove external parameter from the param_card
1828
New_width = ModelVariable(width.name,
1829
'-1 * im(CMASS_%s**2) / %s' % (mass.name, mass.name), 'real', mass.depend)
1830
self.get('parameters')[('external',)].remove(width)
1831
self.add_param(New_param, (mass,))
1832
self.add_param(New_width, (New_param,))
1833
mass.expr = 'cmath.sqrt(re(%s**2))' % mass.expr
1834
to_change[mass.name] = New_param.name
1837
mass.expr = 're(%s)' % mass.expr
1838
self.add_param(New_param, (mass, width))
1839
to_change[mass.name] = New_param.name
1841
# Remove the Yukawa and fix those accordingly to the mass/complex mass
1842
yukawas = [p for p in self.get('parameters')[('external',)]
1843
if p.lhablock.lower() == 'yukawa']
1844
for yukawa in yukawas:
1845
# clean the pevious parameter
1846
self.get('parameters')[('external',)].remove(yukawa)
1848
particle = self.get_particle(yukawa.lhacode[0])
1849
mass = self.get_parameter(particle.get('mass'))
1851
# add the new parameter in the correct category
1852
if mass.depend == ('external',):
1855
depend = mass.depend
1857
New_param = ModelVariable(yukawa.name, mass.name, 'real', depend)
1859
# Add it in the model at the correct place (for the dependences)
1860
if mass.name in to_change:
1861
expr = 'CMASS_%s' % mass.name
1864
param_depend = self.get_parameter(expr)
1865
self.add_param(New_param, [param_depend])
1871
# So at this stage we still need to modify all parameters depending of
1872
# particle's mass. In addition all parameter (but mass/width/external
1873
# parameter) should be pass in complex mode.
1874
pat = '|'.join(list(to_change.keys()))
1875
pat = r'(%s)\b' % pat
1876
pat = re.compile(pat)
1878
return to_change[match.group()]
1880
# Modify the parameters
1881
for dep, list_param in self['parameters'].items():
1882
for param in list_param:
1883
if param.name.startswith('CMASS_') or param.name in mass_widths or\
1884
isinstance(param, ParamCardVariable):
1886
param.type = 'complex'
1887
# print param.expr, to_change
1889
param.expr = pat.sub(replace, param.expr)
1891
# Modify the couplings
1892
for dep, list_coup in self['couplings'].items():
1893
for coup in list_coup:
1894
coup.expr = pat.sub(replace, coup.expr)
1896
def add_param(self, new_param, depend_param):
1897
"""add the parameter in the list of parameter in a correct position"""
1900
for i,param in enumerate(self.get('parameters')[new_param.depend]):
1901
if param.name in depend_param:
1903
self.get('parameters')[new_param.depend].insert(pos, new_param)
1906
#def __repr__(self):
1909
# return "Model(%s)" % self.get_name()
1911
################################################################################
1912
# Class for Parameter / Coupling
1913
################################################################################
1914
class ModelVariable(object):
1915
"""A Class for storing the information about coupling/ parameter"""
1917
def __init__(self, name, expression, type, depend=()):
1918
"""Initialize a new parameter/coupling"""
1921
self.expr = expression # python expression
1922
self.type = type # real/complex
1923
self.depend = depend # depend on some other parameter -tuple-
1926
def __eq__(self, other):
1927
"""Object with same name are identical, If the object is a string we check
1928
if the attribute name is equal to this string"""
1931
return other.name == self.name
1933
return other == self.name
1935
class ParamCardVariable(ModelVariable):
1936
""" A class for storing the information linked to all the parameter
1937
which should be define in the param_card.dat"""
1939
depend = ('external',)
1942
def __init__(self, name, value, lhablock, lhacode):
1943
"""Initialize a new ParamCardVariable
1944
name: name of the variable
1945
value: default numerical value
1946
lhablock: name of the block in the param_card.dat
1947
lhacode: code associate to the variable
1951
self.lhablock = lhablock
1952
self.lhacode = lhacode
1955
#===============================================================================
1956
# Classes used in diagram generation and process definition:
1957
# Leg, Vertex, Diagram, Process
1958
#===============================================================================
1960
#===============================================================================
1962
#===============================================================================
1963
class Leg(PhysicsObject):
1964
"""Leg object: id (Particle), number, I/F state, flag from_group
1967
def default_setup(self):
1968
"""Default values for all properties"""
1972
# state: True = final, False = initial (boolean to save memory)
1973
self['state'] = True
1974
#self['loop_line'] = False
1975
self['loop_line'] = False
1976
# from_group: Used in diagram generation
1977
self['from_group'] = True
1978
# onshell: decaying leg (True), forbidden s-channel (False), none (None)
1979
self['onshell'] = None
1980
# filter on the helicty
1981
self['polarization'] = []
1983
def filter(self, name, value):
1984
"""Filter for valid leg property values."""
1986
if name in ['id', 'number']:
1987
if not isinstance(value, int):
1988
raise self.PhysicsObjectError("%s is not a valid integer for leg id" % str(value))
1990
elif name == 'state':
1991
if not isinstance(value, bool):
1992
raise self.PhysicsObjectError("%s is not a valid leg state (True|False)" % \
1995
elif name == 'from_group':
1996
if not isinstance(value, bool) and value != None:
1997
raise self.PhysicsObjectError("%s is not a valid boolean for leg flag from_group" % \
2000
elif name == 'loop_line':
2001
if not isinstance(value, bool) and value != None:
2002
raise self.PhysicsObjectError("%s is not a valid boolean for leg flag loop_line" % \
2005
elif name == 'onshell':
2006
if not isinstance(value, bool) and value != None:
2007
raise self.PhysicsObjectError("%s is not a valid boolean for leg flag onshell" % \
2010
elif name == 'polarization':
2011
if not isinstance(value, list):
2012
raise self.PhysicsObjectError( \
2013
"%s is not a valid list" % str(value))
2015
if i not in [-1, 1, 2,-2, 3,-3, 0, 99]:
2016
raise self.PhysicsObjectError( \
2017
"%s is not a valid polarization" % str(value))
2021
def get_sorted_keys(self):
2022
"""Return particle property names as a nicely sorted list."""
2024
return ['id', 'number', 'state', 'from_group', 'loop_line', 'onshell', 'polarization']
2026
def is_fermion(self, model):
2027
"""Returns True if the particle corresponding to the leg is a
2030
assert isinstance(model, Model), "%s is not a model" % str(model)
2032
return model.get('particle_dict')[self['id']].is_fermion()
2034
def is_incoming_fermion(self, model):
2035
"""Returns True if leg is an incoming fermion, i.e., initial
2036
particle or final antiparticle"""
2038
assert isinstance(model, Model), "%s is not a model" % str(model)
2040
part = model.get('particle_dict')[self['id']]
2041
return part.is_fermion() and \
2042
(self.get('state') == False and part.get('is_part') or \
2043
self.get('state') == True and not part.get('is_part'))
2045
def is_outgoing_fermion(self, model):
2046
"""Returns True if leg is an outgoing fermion, i.e., initial
2047
antiparticle or final particle"""
2049
assert isinstance(model, Model), "%s is not a model" % str(model)
2051
part = model.get('particle_dict')[self['id']]
2052
return part.is_fermion() and \
2053
(self.get('state') == True and part.get('is_part') or \
2054
self.get('state') == False and not part.get('is_part'))
2056
# Helper function. We don't overload the == operator because it might be useful
2057
# to define it differently than that later.
2059
def same(self, leg):
2060
""" Returns true if the leg in argument has the same ID and the same numer """
2062
# In case we want to check this leg with an integer in the tagging procedure,
2063
# then it only has to match the leg number.
2064
if isinstance(leg,int):
2065
if self['number']==leg:
2070
# If using a Leg object instead, we also want to compare the other relevant
2072
elif isinstance(leg, Leg):
2073
if self['id']==leg.get('id') and \
2074
self['number']==leg.get('number') and \
2075
self['loop_line']==leg.get('loop_line') :
2083
# Make sure sort() sorts lists of legs according to 'number'
2084
def __lt__(self, other):
2085
return self['number'] < other['number']
2087
#===============================================================================
2089
#===============================================================================
2090
class LegList(PhysicsObjectList):
2091
"""List of Leg objects
2094
def is_valid_element(self, obj):
2095
"""Test if object obj is a valid Leg for the list."""
2097
return isinstance(obj, Leg)
2099
# Helper methods for diagram generation
2101
def from_group_elements(self):
2102
"""Return all elements which have 'from_group' True"""
2104
return [leg for leg in self if leg.get('from_group')]
2106
def minimum_one_from_group(self):
2107
"""Return True if at least one element has 'from_group' True"""
2109
return len(self.from_group_elements()) > 0
2111
def minimum_two_from_group(self):
2112
"""Return True if at least two elements have 'from_group' True"""
2114
return len(self.from_group_elements()) > 1
2116
def can_combine_to_1(self, ref_dict_to1):
2117
"""If has at least one 'from_group' True and in ref_dict_to1,
2118
return the return list from ref_dict_to1, otherwise return False"""
2119
if self.minimum_one_from_group():
2120
return tuple(sorted([leg.get('id') for leg in self])) in ref_dict_to1
2124
def can_combine_to_0(self, ref_dict_to0, is_decay_chain=False):
2125
"""If has at least two 'from_group' True and in ref_dict_to0,
2127
return the vertex (with id from ref_dict_to0), otherwise return None
2129
If is_decay_chain = True, we only allow clustering of the
2130
initial leg, since we want this to be the last wavefunction to
2134
# Special treatment - here we only allow combination to 0
2135
# if the initial leg (marked by from_group = None) is
2136
# unclustered, since we want this to stay until the very
2138
return any(leg.get('from_group') == None for leg in self) and \
2139
tuple(sorted([leg.get('id') \
2140
for leg in self])) in ref_dict_to0
2142
if self.minimum_two_from_group():
2143
return tuple(sorted([leg.get('id') for leg in self])) in ref_dict_to0
2147
def get_outgoing_id_list(self, model):
2148
"""Returns the list of ids corresponding to the leglist with
2149
all particles outgoing"""
2153
assert isinstance(model, Model), "Error! model not model"
2157
if leg.get('state') == False:
2158
res.append(model.get('particle_dict')[leg.get('id')].get_anti_pdg_code())
2160
res.append(leg.get('id'))
2164
def sort(self,*args, **opts):
2165
"""Match with FKSLegList"""
2166
Opts=copy.copy(opts)
2167
if 'pert' in list(Opts.keys()):
2169
return super(LegList,self).sort(*args, **Opts)
2172
#===============================================================================
2174
#===============================================================================
2175
class MultiLeg(PhysicsObject):
2176
"""MultiLeg object: ids (Particle or particles), I/F state
2179
def default_setup(self):
2180
"""Default values for all properties"""
2183
self['state'] = True
2184
self['polarization'] = []
2186
def filter(self, name, value):
2187
"""Filter for valid multileg property values."""
2190
if not isinstance(value, list):
2191
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2193
if not isinstance(i, int):
2194
raise self.PhysicsObjectError("%s is not a valid list of integers" % str(value))
2196
if name == 'polarization':
2197
if not isinstance(value, list):
2198
raise self.PhysicsObjectError( \
2199
"%s is not a valid list" % str(value))
2201
if i not in [-1, 1, 2, -2, 3, -3, 0, 99]:
2202
raise self.PhysicsObjectError( \
2203
"%s is not a valid polarization" % str(value))
2206
if not isinstance(value, bool):
2207
raise self.PhysicsObjectError("%s is not a valid leg state (initial|final)" % \
2212
def get_sorted_keys(self):
2213
"""Return particle property names as a nicely sorted list."""
2215
return ['ids', 'state','polarization']
2217
#===============================================================================
2219
#===============================================================================
2220
class MultiLegList(PhysicsObjectList):
2221
"""List of MultiLeg objects
2224
def is_valid_element(self, obj):
2225
"""Test if object obj is a valid MultiLeg for the list."""
2227
return isinstance(obj, MultiLeg)
2229
#===============================================================================
2231
#===============================================================================
2232
class Vertex(PhysicsObject):
2233
"""Vertex: list of legs (ordered), id (Interaction)
2236
sorted_keys = ['id', 'legs']
2238
# This sets what are the ID's of the vertices that must be ignored for the
2239
# purpose of the multi-channeling. 0 and -1 are ID's of various technical
2240
# vertices which have no relevance from the perspective of the diagram
2241
# topology, while -2 is the ID of a vertex that results from a shrunk loop
2242
# (for loop-induced integration with MadEvent) and one may or may not want
2243
# to consider these higher point loops for the purpose of the multi-channeling.
2244
# So, adding -2 to the list below makes sur that all loops are considered
2245
# for multichanneling.
2246
ID_to_veto_for_multichanneling = [0,-1,-2]
2248
# For loop-induced integration, considering channels from up to box loops
2249
# typically leads to better efficiencies. Beyond that, it is detrimental
2250
# because the phase-space generation is not suited to map contact interactions
2251
# This parameter controls up to how many legs should loop-induced diagrams
2252
# be considered for multichanneling.
2253
# Notice that, in the grouped subprocess case mode, if -2 is not added to
2254
# the list ID_to_veto_for_multichanneling then all loop are considered by
2255
# default and the constraint below is not applied.
2256
max_n_loop_for_multichanneling = 4
2259
def default_setup(self):
2260
"""Default values for all properties"""
2262
# The 'id' of the vertex corresponds to the interaction ID it is made of.
2263
# Notice that this 'id' can take the special values :
2264
# -1 : A two-point vertex which either 'sews' the two L-cut particles
2265
# together or simply merges two wavefunctions to create an amplitude
2266
# (in the case of tree-level diagrams).
2267
# -2 : The id given to the ContractedVertices (i.e. a shrunk loop) so
2268
# that it can be easily identified when constructing the DiagramChainLinks.
2270
self['legs'] = LegList()
2272
def filter(self, name, value):
2273
"""Filter for valid vertex property values."""
2276
if not isinstance(value, int):
2277
raise self.PhysicsObjectError("%s is not a valid integer for vertex id" % str(value))
2280
if not isinstance(value, LegList):
2281
raise self.PhysicsObjectError("%s is not a valid LegList object" % str(value))
2285
def get_sorted_keys(self):
2286
"""Return particle property names as a nicely sorted list."""
2288
return self.sorted_keys #['id', 'legs']
2290
def nice_string(self):
2291
"""return a nice string"""
2294
for leg in self['legs']:
2295
mystr.append( str(leg['number']) + '(%s)' % str(leg['id']))
2296
mystr = '(%s,id=%s ,obj_id:%s)' % (', '.join(mystr), self['id'], id(self))
2301
def get_s_channel_id(self, model, ninitial):
2302
"""Returns the id for the last leg as an outgoing
2303
s-channel. Returns 0 if leg is t-channel, or if identity
2304
vertex. Used to check for required and forbidden s-channel
2307
leg = self.get('legs')[-1]
2310
# For one initial particle, all legs are s-channel
2311
# Only need to flip particle id if state is False
2312
if leg.get('state') == True:
2313
return leg.get('id')
2315
return model.get('particle_dict')[leg.get('id')].\
2318
# Number of initial particles is at least 2
2319
if self.get('id') == 0 or \
2320
leg.get('state') == False:
2321
# identity vertex or t-channel particle
2324
if leg.get('loop_line'):
2325
# Loop lines never count as s-channel
2328
# Check if the particle number is <= ninitial
2329
# In that case it comes from initial and we should switch direction
2330
if leg.get('number') > ninitial:
2331
return leg.get('id')
2333
return model.get('particle_dict')[leg.get('id')].\
2336
## Check if the other legs are initial or final.
2337
## If the latter, return leg id, if the former, return -leg id
2338
#if self.get('legs')[0].get('state') == True:
2339
# return leg.get('id')
2341
# return model.get('particle_dict')[leg.get('id')].\
2342
# get_anti_pdg_code()
2344
#===============================================================================
2346
#===============================================================================
2347
class VertexList(PhysicsObjectList):
2348
"""List of Vertex objects
2353
def is_valid_element(self, obj):
2354
"""Test if object obj is a valid Vertex for the list."""
2356
return isinstance(obj, Vertex)
2358
def __init__(self, init_list=None, orders=None):
2359
"""Creates a new list object, with an optional dictionary of
2364
if init_list is not None:
2365
for object in init_list:
2368
if isinstance(orders, dict):
2369
self.orders = orders
2371
#===============================================================================
2373
#===============================================================================
2374
class ContractedVertex(Vertex):
2375
"""ContractedVertex: When contracting a loop to a given vertex, the created
2376
vertex object is then a ContractedVertex object which has additional
2377
information with respect to a regular vertex object. For example, it contains
2378
the PDG of the particles attached to it. (necessary because the contracted
2379
vertex doesn't have an interaction ID which would allow to retrieve such
2383
def default_setup(self):
2384
"""Default values for all properties"""
2387
self['loop_tag'] = tuple()
2388
self['loop_orders'] = {}
2389
super(ContractedVertex, self).default_setup()
2391
def filter(self, name, value):
2392
"""Filter for valid vertex property values."""
2395
if isinstance(value, list):
2397
if not isinstance(elem,int):
2398
raise self.PhysicsObjectError("%s is not a valid integer for leg PDG" % str(elem))
2400
raise self.PhysicsObjectError("%s is not a valid list for contracted vertex PDGs"%str(value))
2401
if name == 'loop_tag':
2402
if isinstance(value, tuple):
2404
if not (isinstance(elem,int) or isinstance(elem,tuple)):
2405
raise self.PhysicsObjectError("%s is not a valid int or tuple for loop tag element"%str(elem))
2407
raise self.PhysicsObjectError("%s is not a valid tuple for a contracted vertex loop_tag."%str(value))
2408
if name == 'loop_orders':
2409
Interaction.filter(Interaction(), 'orders', value)
2411
return super(ContractedVertex, self).filter(name, value)
2415
def get_sorted_keys(self):
2416
"""Return particle property names as a nicely sorted list."""
2418
return super(ContractedVertex, self).get_sorted_keys()+['PDGs']
2420
#===============================================================================
2422
#===============================================================================
2423
class Diagram(PhysicsObject):
2424
"""Diagram: list of vertices (ordered)
2427
def default_setup(self):
2428
"""Default values for all properties"""
2430
self['vertices'] = VertexList()
2433
def filter(self, name, value):
2434
"""Filter for valid diagram property values."""
2436
if name == 'vertices':
2437
if not isinstance(value, VertexList):
2438
raise self.PhysicsObjectError("%s is not a valid VertexList object" % str(value))
2440
if name == 'orders':
2441
Interaction.filter(Interaction(), 'orders', value)
2445
def get_sorted_keys(self):
2446
"""Return particle property names as a nicely sorted list."""
2448
return ['vertices', 'orders']
2450
def nice_string(self):
2451
"""Returns a nicely formatted string of the diagram content."""
2454
removed_index = set()
2456
if self['vertices']:
2458
for vert in self['vertices']:
2461
for leg in vert['legs'][:-1]:
2462
if leg.get('polarization'):
2463
mystr = mystr + str(leg['number']) + '(%s{%s})' % (str(leg['id']),leg['polarization']) + ','
2465
mystr = mystr + str(leg['number']) + '(%s)' % str(leg['id']) + ','
2467
used_leg.append(leg['number'])
2468
if __debug__ and len(used_leg) != len(set(used_leg)):
2470
responsible = id(vert)
2471
if __debug__ and any(l['number'] in removed_index for l in vert['legs']):
2473
responsible = id(vert)
2475
if self['vertices'].index(vert) < len(self['vertices']) - 1:
2476
# Do not want ">" in the last vertex
2477
mystr = mystr[:-1] + '>'
2479
if vert['legs'][-1]['number'] != min([l['number'] for l in vert['legs'][:-1]]):
2481
responsible = id(vert)
2482
for l in vert['legs']:
2483
removed_index.add(l['number'])
2484
removed_index.remove(vert['legs'][-1]['number'])
2486
lastleg = vert['legs'][-1]
2487
if lastleg['polarization']:
2488
mystr = mystr + str(lastleg['number']) + '(%s{%s})' % (str(lastleg['id']), lastleg['polarization']) + ','
2490
mystr = mystr + str(lastleg['number']) + '(%s)' % str(lastleg['id']) + ','
2491
mystr = mystr + 'id:' + str(vert['id']) + '),'
2493
mystr = mystr[:-1] + ')'
2494
mystr += " (%s)" % (",".join(["%s=%d" % (key, self['orders'][key]) \
2495
for key in sorted(self['orders'].keys())]))
2498
raise Exception("invalid diagram: %s. vert_id: %s" % (mystr, responsible))
2504
def calculate_orders(self, model):
2505
"""Calculate the actual coupling orders of this diagram. Note
2506
that the special order WEIGTHED corresponds to the sum of
2507
hierarchys for the couplings."""
2509
coupling_orders = dict([(c, 0) for c in model.get('coupling_orders')])
2511
for vertex in self['vertices']:
2512
if vertex.get('id') in [0,-1]: continue
2513
if vertex.get('id') == -2:
2514
couplings = vertex.get('loop_orders')
2516
couplings = model.get('interaction_dict')[vertex.get('id')].\
2518
for coupling in couplings:
2519
coupling_orders[coupling] += couplings[coupling]
2520
weight += sum([model.get('order_hierarchy')[c]*n for \
2521
(c,n) in couplings.items()])
2522
coupling_orders['WEIGHTED'] = weight
2523
self.set('orders', coupling_orders)
2525
def pass_squared_order_constraints(self, diag_multiplier, squared_orders,
2527
""" Returns wether the contributiong consisting in the current diagram
2528
multiplied by diag_multiplier passes the *positive* squared_orders
2529
specified ( a dictionary ) of types sq_order_types (a dictionary whose
2530
values are the relational operator used to define the constraint of the
2533
for order, value in squared_orders.items():
2536
combined_order = self.get_order(order) + \
2537
diag_multiplier.get_order(order)
2538
if ( sq_orders_types[order]=='==' and combined_order != value ) or \
2539
( sq_orders_types[order] in ['=', '<='] and combined_order > value) or \
2540
( sq_orders_types[order]=='>' and combined_order <= value) :
2544
def get_order(self, order):
2545
"""Return the order of this diagram. It returns 0 if it is not present."""
2548
return self['orders'][order]
2552
def get_contracted_loop_diagram(self, struct_rep=None):
2553
""" Returns a Diagram which correspond to the loop diagram with the
2554
loop shrunk to a point. Of course for a instance of base_objects.Diagram
2555
one must simply return self."""
2559
def get_external_legs(self):
2560
""" Return the list of external legs of this diagram """
2562
external_legs = LegList([])
2563
for leg in sum([vert.get('legs') for vert in self.get('vertices')],[]):
2564
if not leg.get('number') in [l.get('number') for l in external_legs]:
2565
external_legs.append(leg)
2567
return external_legs
2569
def renumber_legs(self, perm_map, leg_list):
2570
"""Renumber legs in all vertices according to perm_map"""
2572
vertices = VertexList()
2573
min_dict = copy.copy(perm_map)
2574
# Dictionary from leg number to state
2575
state_dict = dict([(l.get('number'), l.get('state')) for l in leg_list])
2576
# First renumber all legs in the n-1->1 vertices
2577
for vertex in self.get('vertices')[:-1]:
2578
vertex = copy.copy(vertex)
2579
leg_list = LegList([copy.copy(l) for l in vertex.get('legs')])
2580
for leg in leg_list[:-1]:
2581
leg.set('number', min_dict[leg.get('number')])
2582
leg.set('state', state_dict[leg.get('number')])
2583
min_number = min([leg.get('number') for leg in leg_list[:-1]])
2585
min_dict[leg.get('number')] = min_number
2586
# resulting leg is initial state if there is exactly one
2587
# initial state leg among the incoming legs
2588
state_dict[min_number] = len([l for l in leg_list[:-1] if \
2589
not l.get('state')]) != 1
2590
leg.set('number', min_number)
2591
leg.set('state', state_dict[min_number])
2592
vertex.set('legs', leg_list)
2593
vertices.append(vertex)
2594
# Now renumber the legs in final vertex
2595
vertex = copy.copy(self.get('vertices')[-1])
2596
leg_list = LegList([copy.copy(l) for l in vertex.get('legs')])
2597
for leg in leg_list:
2598
leg.set('number', min_dict[leg.get('number')])
2599
leg.set('state', state_dict[leg.get('number')])
2600
vertex.set('legs', leg_list)
2601
vertices.append(vertex)
2602
# Finally create new diagram
2603
new_diag = copy.copy(self)
2604
new_diag.set('vertices', vertices)
2605
state_dict = {True:'T',False:'F'}
2608
def get_nb_t_channel(self):
2609
"""return number of t-channel propagator in this diagram
2610
This is used to filter multi-channel.
2613
for v in self['vertices'][:-1]:
2614
l = v.get('legs')[-1]
2615
if not l.get('state'):
2622
def get_vertex_leg_numbers(self,
2623
veto_inter_id=Vertex.ID_to_veto_for_multichanneling,
2625
"""Return a list of the number of legs in the vertices for
2627
This function is only used for establishing the multi-channeling, so that
2628
we exclude from it all the fake vertices and the vertices resulting from
2629
shrunk loops (id=-2)"""
2633
max_n_loop = int(Vertex.max_n_loop_for_multichanneling)
2635
res = [len(v.get('legs')) for v in self.get('vertices') if (v.get('id') \
2636
not in veto_inter_id) or (v.get('id')==-2 and
2637
len(v.get('legs'))>max_n_loop)]
2641
def get_num_configs(self, model, ninitial):
2642
"""Return the maximum number of configs from this diagram,
2643
given by 2^(number of non-zero width s-channel propagators)"""
2645
s_channels = [v.get_s_channel_id(model,ninitial) for v in \
2646
self.get('vertices')[:-1]]
2647
num_props = len([i for i in s_channels if i != 0 and \
2648
model.get_particle(i).get('width').lower() != 'zero'])
2655
def get_flow_charge_diff(self, model):
2656
"""return the difference of total diff of charge occuring on the
2657
lofw of the initial parton. return [None,None] if the two initial parton
2658
are connected and the (partial) value if None if the initial parton is
2661
import madgraph.core.drawing as drawing
2662
drawdiag = drawing.FeynmanDiagram(self, model)
2663
drawdiag.load_diagram()
2666
for v in drawdiag.initial_vertex:
2667
init_part = v.lines[0]
2668
if not init_part.is_fermion():
2672
init_charge = model.get_particle(init_part.id).get('charge')
2676
vcurrent = l_last.end
2678
vcurrent = l_last.begin
2680
while not vcurrent.is_external():
2682
raise Exception('wrong diagram')
2683
next_l = [l for l in vcurrent.lines if l is not l_last and l.is_fermion()][0]
2685
if next_v == vcurrent:
2686
next_v = next_l.begin
2687
l_last, vcurrent = next_l, next_v
2688
if vcurrent in drawdiag.initial_vertex:
2691
out.append(model.get_particle(l_last.id).get('charge') - init_charge)
2695
#===============================================================================
2697
#===============================================================================
2698
class DiagramList(PhysicsObjectList):
2699
"""List of Diagram objects
2702
def is_valid_element(self, obj):
2703
"""Test if object obj is a valid Diagram for the list."""
2705
return isinstance(obj, Diagram)
2707
def nice_string(self, indent=0):
2708
"""Returns a nicely formatted string"""
2709
mystr = " " * indent + str(len(self)) + ' diagrams:\n'
2710
for i, diag in enumerate(self):
2711
mystr = mystr + " " * indent + str(i+1) + " " + \
2712
diag.nice_string() + '\n'
2717
def get_max_order(self,order):
2718
""" Return the order of the diagram in the list with the maximum coupling
2719
order for the coupling specified """
2723
if order in list(diag['orders'].keys()):
2724
if max_order==-1 or diag['orders'][order] > max_order:
2725
max_order = diag['orders'][order]
2729
def apply_negative_sq_order(self, ref_diag_list, order, value, order_type):
2730
""" This function returns a fitlered version of the diagram list self
2731
which satisfy the negative squared_order constraint 'order' with negative
2732
value 'value' and of type 'order_type', assuming that the diagram_list
2733
it must be squared against is 'reg_diag_list'. It also returns the
2734
new postive target squared order which correspond to this negative order
2735
constraint. Example: u u~ > d d~ QED^2<=-2 means that one wants to
2736
pick terms only up to the the next-to-leading order contributiong in QED,
2737
which is QED=2 in this case, so that target_order=4 is returned."""
2739
# First we must compute all contributions to that order
2740
target_order = min(ref_diag_list.get_order_values(order))+\
2741
min(self.get_order_values(order))+2*(-value-1)
2743
new_list = self.apply_positive_sq_orders(ref_diag_list,
2744
{order:target_order}, {order:order_type})
2746
return new_list, target_order
2748
def apply_positive_sq_orders(self, ref_diag_list, sq_orders, sq_order_types):
2749
""" This function returns a filtered version of self which contain
2750
only the diagram which satisfy the positive squared order constraints
2751
sq_orders of type sq_order_types and assuming that the diagrams are
2752
multiplied with those of the reference diagram list ref_diag_list."""
2754
new_diag_list = DiagramList()
2755
for tested_diag in self:
2756
for ref_diag in ref_diag_list:
2757
if tested_diag.pass_squared_order_constraints(ref_diag,
2758
sq_orders,sq_order_types):
2759
new_diag_list.append(tested_diag)
2761
return new_diag_list
2763
def filter_constrained_orders(self, order, value, operator):
2764
""" This function modifies the current object and remove the diagram
2765
which do not obey the condition """
2768
for tested_diag in self:
2769
if operator == '==':
2770
if tested_diag['orders'][order] == value:
2771
new.append(tested_diag)
2772
elif operator == '>':
2773
if tested_diag['orders'][order] > value:
2774
new.append(tested_diag)
2779
def get_min_order(self,order):
2780
""" Return the order of the diagram in the list with the mimimum coupling
2781
order for the coupling specified """
2784
if order in list(diag['orders'].keys()):
2785
if min_order==-1 or diag['orders'][order] < min_order:
2786
min_order = diag['orders'][order]
2792
def get_order_values(self, order):
2793
""" Return the list of possible values appearing in the diagrams of this
2794
list for the order given in argument """
2798
if order in list(diag['orders'].keys()):
2799
values.add(diag['orders'][order])
2805
#===============================================================================
2807
#===============================================================================
2808
class Process(PhysicsObject):
2809
"""Process: list of legs (ordered)
2810
dictionary of orders
2815
def default_setup(self):
2816
"""Default values for all properties"""
2818
self['legs'] = LegList()
2819
# These define the orders restrict the born and loop amplitudes.
2821
self['model'] = Model()
2822
# Optional number to identify the process
2824
self['uid'] = 0 # should be a uniq id number
2825
# Required s-channels are given as a list of id lists. Only
2826
# diagrams with all s-channels in any of the lists are
2827
# allowed. This enables generating e.g. Z/gamma as s-channel
2829
self['required_s_channels'] = []
2830
self['forbidden_onsh_s_channels'] = []
2831
self['forbidden_s_channels'] = []
2832
self['forbidden_particles'] = []
2833
self['is_decay_chain'] = False
2834
self['overall_orders'] = {}
2835
# Decay chain processes associated with this process
2836
self['decay_chains'] = ProcessList()
2837
# Legs with decay chains substituted in
2838
self['legs_with_decays'] = LegList()
2839
# Loop particles if the process is to be computed at NLO
2840
self['perturbation_couplings']=[]
2841
# These orders restrict the order of the squared amplitude.
2842
# This dictionary possibly contains a key "WEIGHTED" which
2843
# gives the upper bound for the total weighted order of the
2844
# squared amplitude.
2845
self['squared_orders'] = {}
2846
# The squared order (sqorders) constraints above can either be upper
2847
# bound (<=) or exact match (==) depending on how they were specified
2848
# in the user input. This choice is stored in the dictionary below.
2849
# Notice that the upper bound is the default
2850
self['sqorders_types'] = {}
2851
# other type of constraint at amplitude level
2852
self['constrained_orders'] = {} # {QED: (4,'>')}
2853
self['has_born'] = True
2854
# The NLO_mode is always None for a tree-level process and can be
2855
# 'all', 'real', 'virt' for a loop process.
2856
self['NLO_mode'] = 'tree'
2857
# in the context of QED or QED+QCD perturbation, it is useful to
2858
# keep track of the orders that have been explicitly asked by the
2859
# user, because other borns will appear used for the subtraction
2861
self['born_sq_orders'] = {}
2862
# The user might want to have the individual matrix element evaluations
2863
# for specific values of the coupling orders. The list below specifies
2864
# what are the coupling names which need be individually treated.
2865
# For example, for the process p p > j j [] QED=2 (QED=2 is
2866
# then a squared order constraint), then QED will appear in the
2867
# 'split_orders' list so that the subroutine in matrix.f return the
2868
# evaluation of the matrix element individually for the pure QCD
2869
# contribution 'QCD=4 QED=0', the pure interference 'QCD=2 QED=2' and
2870
# the pure QED contribution of order 'QCD=0 QED=4'.
2871
self['split_orders'] = []
2873
def filter(self, name, value):
2874
"""Filter for valid process property values."""
2876
if name in ['legs', 'legs_with_decays'] :
2877
if not isinstance(value, LegList):
2878
raise self.PhysicsObjectError("%s is not a valid LegList object" % str(value))
2880
if name in ['orders', 'overall_orders','squared_orders', 'born_sq_orders']:
2881
Interaction.filter(Interaction(), 'orders', value)
2883
if name == 'constrained_orders':
2884
if not isinstance(value, dict):
2885
raise self.PhysicsObjectError("%s is not a valid dictionary" % str(value))
2887
if name == 'sqorders_types':
2888
if not isinstance(value, dict):
2889
raise self.PhysicsObjectError("%s is not a valid dictionary" % str(value))
2890
for order in list(value.keys())+list(value.values()):
2891
if not isinstance(order, str):
2892
raise self.PhysicsObjectError("%s is not a valid string" % str(value))
2894
if name == 'split_orders':
2895
if not isinstance(value, list):
2896
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2898
if not isinstance(order, str):
2899
raise self.PhysicsObjectError("%s is not a valid string" % str(value))
2902
if not isinstance(value, Model):
2903
raise self.PhysicsObjectError("%s is not a valid Model object" % str(value))
2904
if name in ['id', 'uid']:
2905
if not isinstance(value, int):
2906
raise self.PhysicsObjectError("Process %s %s is not an integer" % (name, repr(value)))
2908
if name == 'required_s_channels':
2909
if not isinstance(value, list):
2910
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2912
if not isinstance(l, list):
2913
raise self.PhysicsObjectError("%s is not a valid list of lists" % str(value))
2915
if not isinstance(i, int):
2916
raise self.PhysicsObjectError("%s is not a valid list of integers" % str(l))
2918
raise self.PhysicsObjectError("Not valid PDG code %d for s-channel particle" % i)
2920
if name in ['forbidden_onsh_s_channels', 'forbidden_s_channels']:
2921
if not isinstance(value, list):
2922
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2924
if not isinstance(i, int):
2925
raise self.PhysicsObjectError("%s is not a valid list of integers" % str(value))
2927
raise self.PhysicsObjectError("Not valid PDG code %d for s-channel particle" % str(value))
2929
if name == 'forbidden_particles':
2930
if not isinstance(value, list):
2931
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2933
if not isinstance(i, int):
2934
raise self.PhysicsObjectError("%s is not a valid list of integers" % str(value))
2936
raise self.PhysicsObjectError("Forbidden particles should have a positive PDG code" % str(value))
2938
if name == 'perturbation_couplings':
2939
if not isinstance(value, list):
2940
raise self.PhysicsObjectError("%s is not a valid list" % str(value))
2942
if not isinstance(order, str):
2943
raise self.PhysicsObjectError("%s is not a valid string" % str(value))
2945
if name == 'is_decay_chain':
2946
if not isinstance(value, bool):
2947
raise self.PhysicsObjectError("%s is not a valid bool" % str(value))
2949
if name == 'has_born':
2950
if not isinstance(value, bool):
2951
raise self.PhysicsObjectError("%s is not a valid bool" % str(value))
2953
if name == 'decay_chains':
2954
if not isinstance(value, ProcessList):
2955
raise self.PhysicsObjectError("%s is not a valid ProcessList" % str(value))
2957
if name == 'NLO_mode':
2958
import madgraph.interface.madgraph_interface as mg
2959
if value not in mg.MadGraphCmd._valid_nlo_modes:
2960
raise self.PhysicsObjectError("%s is not a valid NLO_mode" % str(value))
2963
def has_multiparticle_label(self):
2964
""" A process, not being a ProcessDefinition never carries multiple
2969
def set(self, name, value):
2970
"""Special set for forbidden particles - set to abs value."""
2972
if name == 'forbidden_particles':
2974
value = [abs(i) for i in value]
2978
if name == 'required_s_channels':
2979
# Required s-channels need to be a list of lists of ids
2980
if value and isinstance(value, list) and \
2981
not isinstance(value[0], list):
2984
return super(Process, self).set(name, value) # call the mother routine
2986
def get_squared_order_type(self, order):
2987
""" Return what kind of squared order constraint was specified for the
2990
if order in list(self['sqorders_types'].keys()):
2991
return self['sqorders_types'][order]
2993
# Default behavior '=' is interpreted as upper bound '<='
2996
def get(self, name):
2997
"""Special get for legs_with_decays"""
2999
if name == 'legs_with_decays':
3000
self.get_legs_with_decays()
3002
if name == 'sqorders_types':
3003
# We must make sure that there is a type for each sqorder defined
3004
for order in self['squared_orders'].keys():
3005
if order not in self['sqorders_types']:
3006
# Then assign its type to the default '='
3007
self['sqorders_types'][order]='='
3009
return super(Process, self).get(name) # call the mother routine
3013
def get_sorted_keys(self):
3014
"""Return process property names as a nicely sorted list."""
3016
return ['legs', 'orders', 'overall_orders', 'squared_orders',
3017
'constrained_orders',
3018
'model', 'id', 'required_s_channels',
3019
'forbidden_onsh_s_channels', 'forbidden_s_channels',
3020
'forbidden_particles', 'is_decay_chain', 'decay_chains',
3021
'legs_with_decays', 'perturbation_couplings', 'has_born',
3022
'NLO_mode', 'split_orders', 'born_sq_orders']
3024
def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturbated=True):
3025
"""Returns a nicely formated string about current process
3026
content. Since the WEIGHTED order is automatically set and added to
3027
the user-defined list of orders, it can be ommitted for some info
3030
if isinstance(prefix, bool) and prefix:
3031
mystr = " " * indent + "Process: "
3032
elif isinstance(prefix, str):
3037
for leg in self['legs']:
3038
mypart = self['model'].get('particle_dict')[leg['id']]
3039
if prevleg and prevleg['state'] == False \
3040
and leg['state'] == True:
3041
# Separate initial and final legs by >
3042
mystr = mystr + '> '
3043
# Add required s-channels
3044
if self['required_s_channels'] and \
3045
self['required_s_channels'][0]:
3046
mystr += "|".join([" ".join([self['model'].\
3047
get('particle_dict')[req_id].get_name() \
3048
for req_id in id_list]) \
3049
for id_list in self['required_s_channels']])
3050
mystr = mystr + ' > '
3052
mystr = mystr + mypart.get_name()
3053
if leg.get('polarization'):
3054
if leg.get('polarization') in [[-1,1],[1,-1]]:
3055
mystr = mystr + '{T} '
3056
elif leg.get('polarization') == [-1]:
3057
mystr = mystr + '{L} '
3058
elif leg.get('polarization') == [1]:
3059
mystr = mystr + '{R} '
3061
mystr = mystr + '{%s} ' %','.join([str(p) for p in leg.get('polarization')])
3064
#mystr = mystr + '(%i) ' % leg['number']
3070
for key in sorted(self['orders'].keys()):
3071
if not print_weighted and key == 'WEIGHTED':
3073
value = int(self['orders'][key])
3074
if key in self['squared_orders']:
3075
if self.get_squared_order_type(key) in ['<=', '==', '='] and \
3076
self['squared_orders'][key] == value:
3078
if self.get_squared_order_type(key) in ['>'] and value == 99:
3080
if key in self['constrained_orders']:
3081
if value == self['constrained_orders'][key][0] and\
3082
self['constrained_orders'][key][1] in ['=', '<=', '==']:
3085
to_add.append('%s=0' % key)
3087
to_add.append('%s<=%s' % (key,value))
3090
mystr = mystr + " ".join(to_add) + ' '
3092
if self['constrained_orders']:
3093
mystr = mystr + " ".join('%s%s%d' % (key,
3094
self['constrained_orders'][key][1], self['constrained_orders'][key][0])
3095
for key in sorted(self['constrained_orders'].keys())) + ' '
3097
# Add perturbation_couplings
3098
if print_perturbated and self['perturbation_couplings']:
3099
mystr = mystr + '[ '
3100
if self['NLO_mode']!='tree':
3101
if self['NLO_mode']=='virt' and not self['has_born']:
3102
mystr = mystr + 'sqrvirt = '
3104
mystr = mystr + self['NLO_mode'] + ' = '
3105
for order in sorted(self['perturbation_couplings']):
3106
mystr = mystr + order + ' '
3107
mystr = mystr + '] '
3109
# Add squared orders
3110
if self['squared_orders']:
3112
for key in sorted(self['squared_orders'].keys()):
3113
if not print_weighted and key == 'WEIGHTED':
3115
if key in self['constrained_orders']:
3116
if self['constrained_orders'][key][0] == self['squared_orders'][key]/2 and \
3117
self['constrained_orders'][key][1] == self.get_squared_order_type(key):
3119
to_add.append(key + '^2%s%d'%\
3120
(self.get_squared_order_type(key),self['squared_orders'][key]))
3123
mystr = mystr + " ".join(to_add) + ' '
3126
# Add forbidden s-channels
3127
if self['forbidden_onsh_s_channels']:
3128
mystr = mystr + '$ '
3129
for forb_id in self['forbidden_onsh_s_channels']:
3130
forbpart = self['model'].get('particle_dict')[forb_id]
3131
mystr = mystr + forbpart.get_name() + ' '
3133
# Add double forbidden s-channels
3134
if self['forbidden_s_channels']:
3135
mystr = mystr + '$$ '
3136
for forb_id in self['forbidden_s_channels']:
3137
forbpart = self['model'].get('particle_dict')[forb_id]
3138
mystr = mystr + forbpart.get_name() + ' '
3140
# Add forbidden particles
3141
if self['forbidden_particles']:
3142
mystr = mystr + '/ '
3143
for forb_id in self['forbidden_particles']:
3144
forbpart = self['model'].get('particle_dict')[forb_id]
3145
mystr = mystr + forbpart.get_name() + ' '
3150
if self.get('id') or self.get('overall_orders'):
3151
mystr += " @%d" % self.get('id')
3152
if self.get('overall_orders'):
3153
mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
3154
for key in sorted(self['orders'])]) + ' '
3156
if not self.get('decay_chains'):
3159
for decay in self['decay_chains']:
3160
mystr = mystr + '\n' + \
3161
decay.nice_string(indent + 2).replace('Process', 'Decay')
3165
def input_string(self):
3166
"""Returns a process string corresponding to the input string
3167
in the command line interface."""
3172
for leg in self['legs']:
3173
mypart = self['model'].get('particle_dict')[leg['id']]
3174
if prevleg and prevleg['state'] == False \
3175
and leg['state'] == True:
3176
# Separate initial and final legs by ">"
3177
mystr = mystr + '> '
3178
# Add required s-channels
3179
if self['required_s_channels'] and \
3180
self['required_s_channels'][0]:
3181
mystr += "|".join([" ".join([self['model'].\
3182
get('particle_dict')[req_id].get_name() \
3183
for req_id in id_list]) \
3184
for id_list in self['required_s_channels']])
3185
mystr = mystr + '> '
3187
mystr = mystr + mypart.get_name()
3188
if leg.get('polarization'):
3189
if leg.get('polarization') in [[-1,1],[1,-1]]:
3190
mystr = mystr + '{T} '
3191
elif leg.get('polarization') == [-1]:
3192
mystr = mystr + '{L} '
3193
elif leg.get('polarization') == [1]:
3194
mystr = mystr + '{R} '
3196
mystr = mystr + '{%s} ' %','.join([str(p) for p in leg.get('polarization')])
3200
#mystr = mystr + '(%i) ' % leg['number']
3204
keys = list(self['orders'].keys())
3205
keys.sort(reverse=True)
3206
mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
3207
for key in keys]) + ' '
3209
# Add squared orders
3210
if self['squared_orders']:
3211
mystr = mystr + " ".join([key + '^2=' + repr(self['squared_orders'][key]) \
3212
for key in self['squared_orders']]) + ' '
3214
# Add perturbation orders
3215
if self['perturbation_couplings']:
3216
mystr = mystr + '[ '
3217
if self['NLO_mode']:
3218
mystr = mystr + self['NLO_mode']
3219
if not self['has_born'] and self['NLO_mode'] != 'noborn':
3220
mystr = mystr + '^2'
3221
mystr = mystr + '= '
3223
for order in sorted(self['perturbation_couplings']):
3224
mystr = mystr + order + ' '
3225
mystr = mystr + '] '
3228
# Add forbidden s-channels
3229
if self['forbidden_onsh_s_channels']:
3230
mystr = mystr + '$ '
3231
for forb_id in self['forbidden_onsh_s_channels']:
3232
forbpart = self['model'].get('particle_dict')[forb_id]
3233
mystr = mystr + forbpart.get_name() + ' '
3235
# Add double forbidden s-channels
3236
if self['forbidden_s_channels']:
3237
mystr = mystr + '$$ '
3238
for forb_id in self['forbidden_s_channels']:
3239
forbpart = self['model'].get('particle_dict')[forb_id]
3240
mystr = mystr + forbpart.get_name() + ' '
3242
# Add forbidden particles
3243
if self['forbidden_particles']:
3244
mystr = mystr + '/ '
3245
for forb_id in self['forbidden_particles']:
3246
forbpart = self['model'].get('particle_dict')[forb_id]
3247
mystr = mystr + forbpart.get_name() + ' '
3252
if self.get('overall_orders'):
3253
mystr += " @%d" % self.get('id')
3254
if self.get('overall_orders'):
3255
mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
3256
for key in sorted(self['orders'])]) + ' '
3258
if not self.get('decay_chains'):
3261
for decay in self['decay_chains']:
3264
if decay.get('decay_chains'):
3267
mystr += ', ' + paren1 + decay.input_string() + paren2
3271
def base_string(self):
3272
"""Returns a string containing only the basic process (w/o decays)."""
3276
for leg in self.get_legs_with_decays():
3277
mypart = self['model'].get('particle_dict')[leg['id']]
3278
if prevleg and prevleg['state'] == False \
3279
and leg['state'] == True:
3280
# Separate initial and final legs by ">"
3281
mystr = mystr + '> '
3282
mystr = mystr + mypart.get_name()
3283
if leg.get('polarization'):
3284
if leg.get('polarization') in [[-1,1],[1,-1]]:
3285
mystr = mystr + '{T} '
3286
elif leg.get('polarization') == [-1]:
3287
mystr = mystr + '{L} '
3288
elif leg.get('polarization') == [1]:
3289
mystr = mystr + '{R} '
3291
mystr = mystr + '{%s} ' %','.join([str(p) for p in leg.get('polarization')])
3299
def shell_string(self, schannel=True, forbid=True, main=True, pdg_order=False,
3301
"""Returns process as string with '~' -> 'x', '>' -> '_',
3302
'+' -> 'p' and '-' -> 'm', including process number,
3303
intermediate s-channels and forbidden particles,
3304
pdg_order allow to order to leg order by pid."""
3307
if not self.get('is_decay_chain') and print_id:
3308
mystr += "%d_" % self['id']
3312
legs = [l for l in self['legs'][1:]]
3313
legs.sort(key=lambda x: x.get('id'))
3314
legs.insert(0, self['legs'][0])
3320
mypart = self['model'].get('particle_dict')[leg['id']]
3321
if prevleg and prevleg['state'] == False \
3322
and leg['state'] == True:
3323
# Separate initial and final legs by ">"
3325
# Add required s-channels
3326
if self['required_s_channels'] and \
3327
self['required_s_channels'][0] and schannel:
3328
mystr += "_or_".join(["".join([self['model'].\
3329
get('particle_dict')[req_id].get_name() \
3330
for req_id in id_list]) \
3331
for id_list in self['required_s_channels']])
3333
if mypart['is_part']:
3334
mystr = mystr + mypart['name']
3336
mystr = mystr + mypart['antiname']
3337
if leg.get('polarization'):
3338
if leg.get('polarization') in [[-1,1],[1,-1]]:
3340
elif leg.get('polarization') == [-1]:
3342
elif leg.get('polarization') == [1]:
3345
mystr = mystr + '%s ' %''.join([str(p).replace('-','m') for p in leg.get('polarization')])
3349
# Check for forbidden particles
3350
if self['forbidden_particles'] and forbid:
3351
mystr = mystr + '_no_'
3352
for forb_id in self['forbidden_particles']:
3353
forbpart = self['model'].get('particle_dict')[forb_id]
3354
mystr = mystr + forbpart.get_name()
3356
# Replace '~' with 'x'
3357
mystr = mystr.replace('~', 'x')
3358
# Replace '+' with 'p'
3359
mystr = mystr.replace('+', 'p')
3360
# Replace '-' with 'm'
3361
mystr = mystr.replace('-', 'm')
3362
# Just to be safe, remove all spaces
3363
mystr = mystr.replace(' ', '')
3365
for decay in self.get('decay_chains'):
3366
mystr = mystr + "_" + decay.shell_string(schannel,forbid, main=False,
3367
pdg_order=pdg_order)
3369
# Too long name are problematic so restrict them to a maximal of 70 char
3370
if len(mystr) > 64 and main:
3371
if schannel and forbid:
3372
out = self.shell_string(True, False, True, pdg_order)
3374
out = self.shell_string(False, False, True, pdg_order)
3377
if not out.endswith('_%s' % self['uid']):
3378
out += '_%s' % self['uid']
3383
def shell_string_v4(self):
3384
"""Returns process as v4-compliant string with '~' -> 'x' and
3387
mystr = "%d_" % self['id']
3389
for leg in self.get_legs_with_decays():
3390
mypart = self['model'].get('particle_dict')[leg['id']]
3391
if prevleg and prevleg['state'] == False \
3392
and leg['state'] == True:
3393
# Separate initial and final legs by ">"
3395
if mypart['is_part']:
3396
mystr = mystr + mypart['name']
3398
mystr = mystr + mypart['antiname']
3399
if leg.get('polarization'):
3400
if leg.get('polarization') in [[-1,1],[1,-1]]:
3402
elif leg.get('polarization') == [-1]:
3404
elif leg.get('polarization') == [1]:
3407
mystr = mystr + '%s ' %''.join([str(p).replace('-','m') for p in leg.get('polarization')])
3411
# Replace '~' with 'x'
3412
mystr = mystr.replace('~', 'x')
3413
# Just to be safe, remove all spaces
3414
mystr = mystr.replace(' ', '')
3420
def are_negative_orders_present(self):
3421
""" Check iteratively that no coupling order constraint include negative
3424
if any(val<0 for val in list(self.get('orders').values())+\
3425
list(self.get('squared_orders').values())):
3428
for procdef in self['decay_chains']:
3429
if procdef.are_negative_orders_present():
3434
def are_decays_perturbed(self):
3435
""" Check iteratively that the decayed processes are not perturbed """
3437
for procdef in self['decay_chains']:
3438
if procdef['perturbation_couplings'] or procdef.are_decays_perturbed():
3442
def decays_have_squared_orders(self):
3443
""" Check iteratively that the decayed processes are not perturbed """
3445
for procdef in self['decay_chains']:
3446
if procdef['squared_orders']!={} or procdef.decays_have_squared_orders():
3450
def get_ninitial(self):
3451
"""Gives number of initial state particles"""
3453
return len([leg for leg in self.get('legs') if leg.get('state') == False])
3455
def get_initial_ids(self):
3456
"""Gives the pdg codes for initial state particles"""
3458
return [leg.get('id') for leg in \
3459
[leg for leg in self.get('legs') if leg.get('state') == False]]
3461
def get_initial_pdg(self, number):
3462
"""Return the pdg codes for initial state particles for beam number"""
3464
legs = [leg for leg in self.get('legs') if leg.get('state') == False and\
3465
leg.get('number') == number]
3469
return legs[0].get('id')
3471
def get_initial_final_ids(self):
3472
"""return a tuple of two tuple containing the id of the initial/final
3473
state particles. Each list is ordered"""
3476
final = [l.get('id') for l in self.get('legs')\
3477
if l.get('state') or initial.append(l.get('id'))]
3480
return (tuple(initial), tuple(final))
3482
def get_initial_final_ids_after_decay(self, max_depth=-1):
3483
"""return a tuple of two tuple containing the id of the initial/final
3484
state particles. Each list is ordered"""
3486
initial = [l.get('id') for l in self.get('legs')\
3487
if not l.get('state')]
3488
final = self.get_final_ids_after_decay(max_depth=max_depth)
3491
return (tuple(initial), tuple(final))
3494
def get_final_ids_after_decay(self, max_depth=-1):
3495
"""Give the pdg code of the process including decay"""
3497
finals = self.get_final_ids()
3499
for proc in self.get('decay_chains'):
3500
init = proc.get_initial_ids()[0]
3503
pos = finals.index(init)
3506
finals[pos] = proc.get_final_ids_after_decay(max_depth-1)
3509
if isinstance(d, list):
3517
def get_final_legs(self):
3518
"""Gives the final state legs"""
3520
return [leg for leg in self.get('legs') if leg.get('state') == True]
3522
def get_final_ids(self):
3523
"""Gives the pdg codes for final state particles"""
3525
return [l.get('id') for l in self.get_final_legs()]
3528
def get_legs_with_decays(self):
3529
"""Return process with all decay chains substituted in."""
3531
if self['legs_with_decays']:
3532
return self['legs_with_decays']
3534
legs = copy.deepcopy(self.get('legs'))
3535
org_decay_chains = copy.copy(self.get('decay_chains'))
3536
sorted_decay_chains = []
3537
# Sort decay chains according to leg order
3539
if not leg.get('state'): continue
3540
org_ids = [l.get('legs')[0].get('id') for l in \
3542
if leg.get('id') in org_ids:
3543
sorted_decay_chains.append(org_decay_chains.pop(\
3544
org_ids.index(leg.get('id'))))
3545
assert not org_decay_chains
3547
for decay in sorted_decay_chains:
3548
while legs[ileg].get('state') == False or \
3549
legs[ileg].get('id') != decay.get('legs')[0].get('id'):
3551
decay_legs = decay.get_legs_with_decays()
3552
legs = legs[:ileg] + decay_legs[1:] + legs[ileg+1:]
3553
ileg = ileg + len(decay_legs) - 1
3555
# Replace legs with copies
3556
legs = [copy.copy(l) for l in legs]
3558
for ileg, leg in enumerate(legs):
3559
leg.set('number', ileg + 1)
3561
self['legs_with_decays'] = LegList(legs)
3563
return self['legs_with_decays']
3566
"""return the tag for standalone call"""
3568
initial = [] #filled in the next line
3569
final = [l.get('id') for l in self.get('legs')\
3570
if l.get('state') or initial.append(l.get('id'))]
3571
decay_finals = self.get_final_ids_after_decay()
3573
tag = (tuple(initial), tuple(decay_finals))
3577
def list_for_sort(self):
3578
"""Output a list that can be compared to other processes as:
3579
[id, sorted(initial leg ids), sorted(final leg ids),
3580
sorted(decay list_for_sorts)]"""
3582
sorted_list = [self.get('id'),
3583
sorted(self.get_initial_ids()),
3584
sorted(self.get_final_ids())]
3586
if self.get('decay_chains'):
3587
sorted_list.extend(sorted([d.list_for_sort() for d in \
3588
self.get('decay_chains')]))
3592
def compare_for_sort(self, other):
3593
"""Sorting routine which allows to sort processes for
3594
comparison. Compare only process id and legs."""
3596
if self.list_for_sort() > other.list_for_sort():
3598
if self.list_for_sort() < other.list_for_sort():
3600
assert self.list_for_sort() == other.list_for_sort()
3603
def identical_particle_factor(self):
3604
"""Calculate the denominator factor for identical final state particles
3607
final_legs = [leg for leg in self.get_legs_with_decays() if leg.get('state') == True]
3609
identical_indices = collections.defaultdict(int)
3610
for leg in final_legs:
3611
key = (leg.get('id'), tuple(leg.get('polarization')))
3612
identical_indices[key] += 1
3615
return reduce(lambda x, y: x * y, [ math.factorial(val) for val in \
3616
identical_indices.values() ], 1)
3618
def check_expansion_orders(self):
3619
"""Ensure that maximum expansion orders from the model are
3620
properly taken into account in the process"""
3622
# Ensure that expansion orders are taken into account
3623
expansion_orders = self.get('model').get('expansion_order')
3624
orders = self.get('orders')
3625
sq_orders = self.get('squared_orders')
3627
tmp = [(k,v) for (k,v) in expansion_orders.items() if 0 < v < 99]
3631
if k in list(sq_orders.keys()) and \
3632
(sq_orders[k]>v or sq_orders[k]<0):
3634
'''The process with the squared coupling order (%s^2%s%s) specified can potentially
3635
recieve contributions with powers of the coupling %s larger than the maximal
3636
value allowed by the model builder (%s). Hence, MG5_aMC sets the amplitude order
3637
for that coupling to be this maximal one. '''%(k,self.get('sqorders_types')[k],
3638
self.get('squared_orders')[k],k,v))
3641
'''The coupling order (%s=%s) specified is larger than the one allowed
3642
by the model builder. The maximal value allowed is %s.
3643
We set the %s order to this value''' % (k,orders[k],v,k))
3648
def __eq__(self, other):
3649
"""Overloading the equality operator, so that only comparison
3650
of process id and legs is being done, using compare_for_sort."""
3652
if not isinstance(other, Process):
3655
#misc.sprint("can we speed up this computation? Yes we can!")
3656
return self.compare_for_sort(other) == 0
3657
return self.list_for_sort() == other.list_for_sort()
3659
def __ne__(self, other):
3660
return not self.__eq__(other)
3662
#===============================================================================
3664
#===============================================================================
3665
class ProcessList(PhysicsObjectList):
3666
"""List of Process objects
3669
def is_valid_element(self, obj):
3670
"""Test if object obj is a valid Process for the list."""
3672
return isinstance(obj, Process)
3674
def nice_string(self, indent = 0):
3675
"""Returns a nicely formatted string of the matrix element processes."""
3677
mystr = "\n".join([p.nice_string(indent) for p in self])
3681
#===============================================================================
3683
#===============================================================================
3684
class ProcessDefinition(Process):
3685
"""ProcessDefinition: list of multilegs (ordered)
3686
dictionary of orders
3691
def default_setup(self):
3692
"""Default values for all properties"""
3694
super(ProcessDefinition, self).default_setup()
3696
self['legs'] = MultiLegList()
3697
# Decay chain processes associated with this process
3698
self['decay_chains'] = ProcessDefinitionList()
3699
if 'legs_with_decays' in self: del self['legs_with_decays']
3701
def filter(self, name, value):
3702
"""Filter for valid process property values."""
3705
if not isinstance(value, MultiLegList):
3706
raise self.PhysicsObjectError("%s is not a valid MultiLegList object" % str(value))
3707
elif name == 'decay_chains':
3708
if not isinstance(value, ProcessDefinitionList):
3709
raise self.PhysicsObjectError("%s is not a valid ProcessDefinitionList" % str(value))
3712
return super(ProcessDefinition, self).filter(name, value)
3716
def has_multiparticle_label(self):
3717
""" Check that this process definition will yield a single process, as
3718
each multileg only has one leg"""
3720
for process in self['decay_chains']:
3721
if process.has_multiparticle_label():
3724
for mleg in self['legs']:
3725
if len(mleg['ids'])>1:
3730
def check_polarization(self):
3731
""" raise a critical information if someone tries something like
3733
return True if no issue and False if some issue is found
3737
for leg in self.get('legs'):
3738
if not leg.get('state'):
3740
if leg.get('polarization'):
3741
for pid in leg.get('ids'):
3743
pol[pid] = [leg.get('polarization')]
3744
elif leg.get('polarization') in pol[pid]:
3745
# already present polarization -> no issue
3748
for p in leg.get('polarization'):
3749
if any(p in o for o in pol[pid]):
3751
pol[pid].append(leg.get('polarization'))
3753
for pid in leg.get('ids'):
3755
pol[pid] = [list(range(-3,4))]
3756
elif pol[pid] == [list(range(-3,4))]:
3763
def get_sorted_keys(self):
3764
"""Return process property names as a nicely sorted list."""
3766
keys = super(ProcessDefinition, self).get_sorted_keys()
3767
keys.remove('legs_with_decays')
3771
def get_minimum_WEIGHTED(self):
3772
"""Retrieve the minimum starting guess for WEIGHTED order, to
3773
use in find_optimal_process_orders in MultiProcess diagram
3774
generation (as well as particles and hierarchy). The algorithm:
3776
1) Pick out the legs in the multiprocess according to the
3777
highest hierarchy represented (so don't mix particles from
3778
different hierarchy classes in the same multiparticles!)
3780
2) Find the starting maximum WEIGHTED order as the sum of the
3781
highest n-2 weighted orders
3783
3) Pick out required s-channel particle hierarchies, and use
3784
the highest of the maximum WEIGHTED order from the legs and
3785
the minimum WEIGHTED order extracted from 2*s-channel
3786
hierarchys plus the n-2-2*(number of s-channels) lowest
3787
leg weighted orders.
3790
model = self.get('model')
3792
# Extract hierarchy and particles corresponding to the
3793
# different hierarchy levels from the model
3794
particles, hierarchy = model.get_particles_hierarchy()
3795
# Find legs corresponding to the different orders
3796
# making sure we look at lowest hierarchy first for each leg
3798
new_legs = copy.copy(self.get('legs'))
3799
import madgraph.core.base_objects as base_objects
3800
for parts, value in zip(particles, hierarchy):
3802
while ileg < len(new_legs):
3803
if any([id in parts for id in new_legs[ileg].get('ids')]):
3804
max_order_now.append(value)
3809
# Now remove the two lowest orders to get maximum (since the
3810
# number of interactions is n-2)
3811
max_order_now = sorted(max_order_now)[2:]
3813
# Find s-channel propagators corresponding to the different orders
3815
for idlist in self.get('required_s_channels'):
3816
max_order_prop.append([0,0])
3818
for parts, value in zip(particles, hierarchy):
3820
max_order_prop[-1][0] += 2*value
3821
max_order_prop[-1][1] += 1
3825
if len(max_order_prop) >1:
3826
max_order_prop = min(*max_order_prop, key=lambda x:x[0])
3828
max_order_prop = max_order_prop[0]
3830
# Use either the max_order from the external legs or
3831
# the maximum order from the s-channel propagators, plus
3832
# the appropriate lowest orders from max_order_now
3833
max_order_now = max(sum(max_order_now),
3834
max_order_prop[0] + \
3835
sum(max_order_now[:-2 * max_order_prop[1]]))
3837
max_order_now = sum(max_order_now)
3839
return max_order_now, particles, hierarchy
3842
"""basic way to loop over all the process definition.
3843
not used by MG which used some smarter version (use by ML)"""
3845
isids = [leg['ids'] for leg in self['legs'] \
3846
if leg['state'] == False]
3847
fsids = [leg['ids'] for leg in self['legs'] \
3848
if leg['state'] == True]
3851
# Generate all combinations for the initial state
3852
for prod in itertools.product(*isids):
3853
islegs = [Leg({'id':id, 'state': False}) for id in prod]
3854
if tuple(sorted(prod)) in red_isidlist:
3856
red_isidlist.append(tuple(sorted(prod)))
3858
for prod in itertools.product(*fsids):
3859
# Remove double counting between final states
3860
if tuple(sorted(prod)) in red_fsidlist:
3862
red_fsidlist.append(tuple(sorted(prod)))
3863
leg_list = [copy.copy(leg) for leg in islegs]
3864
leg_list.extend([Leg({'id':id, 'state': True}) for id in prod])
3865
legs = LegList(leg_list)
3866
process = self.get_process_with_legs(legs)
3869
def nice_string(self, indent=0, print_weighted=False, prefix=True):
3870
"""Returns a nicely formated string about current process
3874
mystr = " " * indent + "Process: "
3878
for leg in self['legs']:
3880
"/".join([self['model'].get('particle_dict')[id].get_name() \
3881
for id in leg.get('ids')])
3882
if prevleg and prevleg['state'] == False \
3883
and leg['state'] == True:
3884
# Separate initial and final legs by ">"
3885
mystr = mystr + '> '
3886
# Add required s-channels
3887
if self['required_s_channels'] and \
3888
self['required_s_channels'][0]:
3889
mystr += "|".join([" ".join([self['model'].\
3890
get('particle_dict')[req_id].get_name() \
3891
for req_id in id_list]) \
3892
for id_list in self['required_s_channels']])
3893
mystr = mystr + '> '
3895
mystr = mystr + myparts
3896
if leg.get('polarization'):
3897
if leg.get('polarization') in [[-1,1],[1,-1]]:
3898
mystr = mystr + '{T}'
3899
elif leg.get('polarization') == [-1]:
3900
mystr = mystr + '{L}'
3901
elif leg.get('polarization') == [1]:
3902
mystr = mystr + '{R}'
3904
mystr = mystr + '{%s} ' %''.join([str(p) for p in leg.get('polarization')])
3907
#mystr = mystr + '(%i) ' % leg['number']
3910
# Add forbidden s-channels
3911
if self['forbidden_onsh_s_channels']:
3912
mystr = mystr + '$ '
3913
for forb_id in self['forbidden_onsh_s_channels']:
3914
forbpart = self['model'].get('particle_dict')[forb_id]
3915
mystr = mystr + forbpart.get_name() + ' '
3917
# Add double forbidden s-channels
3918
if self['forbidden_s_channels']:
3919
mystr = mystr + '$$ '
3920
for forb_id in self['forbidden_s_channels']:
3921
forbpart = self['model'].get('particle_dict')[forb_id]
3922
mystr = mystr + forbpart.get_name() + ' '
3924
# Add forbidden particles
3925
if self['forbidden_particles']:
3926
mystr = mystr + '/ '
3927
for forb_id in self['forbidden_particles']:
3928
forbpart = self['model'].get('particle_dict')[forb_id]
3929
mystr = mystr + forbpart.get_name() + ' '
3932
mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
3933
for key in sorted(self['orders'])]) + ' '
3935
if self['constrained_orders']:
3936
mystr = mystr + " ".join('%s%s%d' % (key, operator, value) for
3937
(key,(value, operator))
3938
in self['constrained_orders'].items()) + ' '
3940
# Add perturbation_couplings
3941
if self['perturbation_couplings']:
3942
mystr = mystr + '[ '
3943
if self['NLO_mode']!='tree':
3944
if self['NLO_mode']=='virt' and not self['has_born']:
3945
mystr = mystr + 'sqrvirt = '
3947
mystr = mystr + self['NLO_mode'] + ' = '
3948
for order in self['perturbation_couplings']:
3949
mystr = mystr + order + ' '
3950
mystr = mystr + '] '
3952
if self['squared_orders']:
3953
mystr = mystr + " ".join([key + '^2%s%d'%\
3954
(self.get_squared_order_type(key),self['squared_orders'][key]) \
3955
for key in self['squared_orders'].keys() \
3956
if print_weighted or key!='WEIGHTED']) + ' '
3961
if self.get('id') or self.get('overall_orders'):
3962
mystr += " @%d" % self.get('id')
3963
if self.get('overall_orders'):
3964
mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
3965
for key in sorted(self['orders'])]) + ' '
3967
if not self.get('decay_chains'):
3970
for decay in self['decay_chains']:
3971
mystr = mystr + '\n' + \
3972
decay.nice_string(indent + 2).replace('Process', 'Decay')
3976
def get_process_with_legs(self, LegList):
3977
""" Return a Process object which has the same properties of this
3978
ProcessDefinition but with the specified LegList as legs attribute.
3983
'model':self.get('model'),
3984
'id': self.get('id'),
3985
'orders': self.get('orders'),
3986
'sqorders_types': self.get('sqorders_types'),
3987
'squared_orders': self.get('squared_orders'),
3988
'constrained_orders': self.get('constrained_orders'),
3989
'has_born': self.get('has_born'),
3990
'required_s_channels': self.get('required_s_channels'),
3991
'forbidden_onsh_s_channels': self.get('forbidden_onsh_s_channels'),
3992
'forbidden_s_channels': self.get('forbidden_s_channels'),
3993
'forbidden_particles': self.get('forbidden_particles'),
3994
'perturbation_couplings': self.get('perturbation_couplings'),
3995
'is_decay_chain': self.get('is_decay_chain'),
3996
'overall_orders': self.get('overall_orders'),
3997
'split_orders': self.get('split_orders'),
3998
'born_sq_orders': self.get('born_sq_orders'),
3999
'NLO_mode': self.get('NLO_mode')
4002
def get_process(self, initial_state_ids, final_state_ids):
4003
""" Return a Process object which has the same properties of this
4004
ProcessDefinition but with the specified given leg ids. """
4006
# First make sure that the desired particle ids belong to those defined
4007
# in this process definition.
4009
my_isids = [leg.get('ids') for leg in self.get('legs') \
4010
if not leg.get('state')]
4011
my_fsids = [leg.get('ids') for leg in self.get('legs') \
4012
if leg.get('state')]
4013
for i, is_id in enumerate(initial_state_ids):
4014
assert is_id in my_isids[i]
4015
for i, fs_id in enumerate(final_state_ids):
4016
assert fs_id in my_fsids[i]
4018
return self.get_process_with_legs(LegList(\
4019
[Leg({'id': id, 'state':False, 'polarization':[]}) for id in initial_state_ids] + \
4020
[Leg({'id': id, 'state':True, 'polarization':[]}) for id in final_state_ids]))
4022
def __eq__(self, other):
4023
"""Overloading the equality operator, so that only comparison
4024
of process id and legs is being done, using compare_for_sort."""
4026
return super(Process, self).__eq__(other)
4028
#===============================================================================
4029
# ProcessDefinitionList
4030
#===============================================================================
4031
class ProcessDefinitionList(PhysicsObjectList):
4032
"""List of ProcessDefinition objects
4035
def is_valid_element(self, obj):
4036
"""Test if object obj is a valid ProcessDefinition for the list."""
4038
return isinstance(obj, ProcessDefinition)
4040
#===============================================================================
4041
# Global helper functions
4042
#===============================================================================
4044
def make_unique(doubletlist):
4045
"""Make sure there are no doublets in the list doubletlist.
4046
Note that this is a slow implementation, so don't use if speed
4049
assert isinstance(doubletlist, list), \
4050
"Argument to make_unique must be list"
4054
for elem in doubletlist:
4055
if elem not in uniquelist:
4056
uniquelist.append(elem)
4058
doubletlist[:] = uniquelist[:]
added
removed
madgraph/core/base_objects.py
