Package madgraph :: Package core :: Module base_objects
[hide private]
[frames] | no frames]

Source Code for Module madgraph.core.base_objects

   1  ################################################################################ 
   2  # 
   3  # Copyright (c) 2009 The MadGraph Development team and Contributors 
   4  # 
   5  # This file is a part of the MadGraph 5 project, an application which  
   6  # automatically generates Feynman diagrams and matrix elements for arbitrary 
   7  # high-energy processes in the Standard Model and beyond. 
   8  # 
   9  # It is subject to the MadGraph license which should accompany this  
  10  # distribution. 
  11  # 
  12  # For more information, please visit: http://madgraph.phys.ucl.ac.be 
  13  # 
  14  ################################################################################ 
  15   
  16  """Definitions of all basic objects used in the core code: particle,  
  17  interaction, model, leg, vertex, process, ...""" 
  18   
  19  import copy 
  20  import itertools 
  21  import logging 
  22  import math 
  23  import numbers 
  24  import os 
  25  import re 
  26  import StringIO 
  27  import madgraph.core.color_algebra as color 
  28  from madgraph import MadGraph5Error, MG5DIR 
  29   
  30  logger = logging.getLogger('madgraph.base_objects') 
31 32 #=============================================================================== 33 # PhysicsObject 34 #=============================================================================== 35 -class PhysicsObject(dict):
36 """A parent class for all physics objects.""" 37
38 - class PhysicsObjectError(Exception):
39 """Exception raised if an error occurs in the definition 40 or the execution of a physics object.""" 41 pass
42
43 - def __init__(self, init_dict={}):
44 """Creates a new particle object. If a dictionary is given, tries to 45 use it to give values to properties.""" 46 47 dict.__init__(self) 48 self.default_setup() 49 50 assert isinstance(init_dict, dict), \ 51 "Argument %s is not a dictionary" % repr(init_dict) 52 53 54 for item in init_dict.keys(): 55 self.set(item, init_dict[item])
56 57
58 - def __getitem__(self, name):
59 """ force the check that the property exist before returning the 60 value associated to value. This ensure that the correct error 61 is always raise 62 """ 63 64 try: 65 return dict.__getitem__(self, name) 66 except KeyError: 67 self.is_valid_prop(name) #raise the correct error
68 69
70 - def default_setup(self):
71 """Function called to create and setup default values for all object 72 properties""" 73 pass
74
75 - def is_valid_prop(self, name):
76 """Check if a given property name is valid""" 77 78 assert isinstance(name, str), \ 79 "Property name %s is not a string" % repr(name) 80 81 if name not in self.keys(): 82 raise self.PhysicsObjectError, \ 83 "%s is not a valid property for this object" % name 84 85 return True
86
87 - def get(self, name):
88 """Get the value of the property name.""" 89 90 return self[name]
91
92 - def set(self, name, value, force=False):
93 """Set the value of the property name. First check if value 94 is a valid value for the considered property. Return True if the 95 value has been correctly set, False otherwise.""" 96 97 if not __debug__ or force: 98 self[name] = value 99 return True 100 101 if self.is_valid_prop(name): 102 try: 103 self.filter(name, value) 104 self[name] = value 105 return True 106 except self.PhysicsObjectError, why: 107 logger.warning("Property " + name + " cannot be changed:" + \ 108 str(why)) 109 return False
110
111 - def filter(self, name, value):
112 """Checks if the proposed value is valid for a given property 113 name. Returns True if OK. Raises an error otherwise.""" 114 115 return True
116
117 - def get_sorted_keys(self):
118 """Returns the object keys sorted in a certain way. By default, 119 alphabetical.""" 120 121 return self.keys().sort()
122
123 - def __str__(self):
124 """String representation of the object. Outputs valid Python 125 with improved format.""" 126 127 mystr = '{\n' 128 for prop in self.get_sorted_keys(): 129 if isinstance(self[prop], str): 130 mystr = mystr + ' \'' + prop + '\': \'' + \ 131 self[prop] + '\',\n' 132 elif isinstance(self[prop], float): 133 mystr = mystr + ' \'' + prop + '\': %.2f,\n' % self[prop] 134 else: 135 mystr = mystr + ' \'' + prop + '\': ' + \ 136 repr(self[prop]) + ',\n' 137 mystr = mystr.rstrip(',\n') 138 mystr = mystr + '\n}' 139 140 return mystr
141 142 __repr__ = __str__
143
144 145 #=============================================================================== 146 # PhysicsObjectList 147 #=============================================================================== 148 -class PhysicsObjectList(list):
149 """A class to store lists of physics object.""" 150
151 - class PhysicsObjectListError(Exception):
152 """Exception raised if an error occurs in the definition 153 or execution of a physics object list.""" 154 pass
155
156 - def __init__(self, init_list=None):
157 """Creates a new particle list object. If a list of physics 158 object is given, add them.""" 159 160 list.__init__(self) 161 162 if init_list is not None: 163 for object in init_list: 164 self.append(object)
165
166 - def append(self, object):
167 """Appends an element, but test if valid before.""" 168 169 assert self.is_valid_element(object), \ 170 "Object %s is not a valid object for the current list" % repr(object) 171 172 list.append(self, object)
173
174 - def is_valid_element(self, obj):
175 """Test if object obj is a valid element for the list.""" 176 return True
177
178 - def __str__(self):
179 """String representation of the physics object list object. 180 Outputs valid Python with improved format.""" 181 182 mystr = '[' 183 184 for obj in self: 185 mystr = mystr + str(obj) + ',\n' 186 187 mystr = mystr.rstrip(',\n') 188 189 return mystr + ']'
190
191 #=============================================================================== 192 # Particle 193 #=============================================================================== 194 -class Particle(PhysicsObject):
195 """The particle object containing the whole set of information required to 196 univocally characterize a given type of physical particle: name, spin, 197 color, mass, width, charge,... The is_part flag tells if the considered 198 particle object is a particle or an antiparticle. The self_antipart flag 199 tells if the particle is its own antiparticle.""" 200 201 sorted_keys = ['name', 'antiname', 'spin', 'color', 202 'charge', 'mass', 'width', 'pdg_code', 203 'texname', 'antitexname', 'line', 'propagating', 204 'is_part', 'self_antipart'] 205
206 - def default_setup(self):
207 """Default values for all properties""" 208 209 self['name'] = 'none' 210 self['antiname'] = 'none' 211 self['spin'] = 1 212 self['color'] = 1 213 self['charge'] = 1. 214 self['mass'] = 'zero' 215 self['width'] = 'zero' 216 self['pdg_code'] = 0 217 self['texname'] = 'none' 218 self['antitexname'] = 'none' 219 self['line'] = 'dashed' 220 self['propagating'] = True 221 self['is_part'] = True 222 self['self_antipart'] = False
223
224 - def filter(self, name, value):
225 """Filter for valid particle property values.""" 226 227 if name in ['name', 'antiname']: 228 # Forbid special character but +-~_ 229 p=re.compile('''^[\w\-\+~_]+$''') 230 if not p.match(value): 231 raise self.PhysicsObjectError, \ 232 "%s is not a valid particle name" % value 233 234 if name is 'spin': 235 if not isinstance(value, int): 236 raise self.PhysicsObjectError, \ 237 "Spin %s is not an integer" % repr(value) 238 if value < 1 or value > 5: 239 raise self.PhysicsObjectError, \ 240 "Spin %i is smaller than one" % value 241 242 if name is 'color': 243 if not isinstance(value, int): 244 raise self.PhysicsObjectError, \ 245 "Color %s is not an integer" % repr(value) 246 if value not in [1, 3, 6, 8]: 247 raise self.PhysicsObjectError, \ 248 "Color %i is not valid" % value 249 250 if name in ['mass', 'width']: 251 # Must start with a letter, followed by letters, digits or _ 252 p = re.compile('\A[a-zA-Z]+[\w\_]*\Z') 253 if not p.match(value): 254 raise self.PhysicsObjectError, \ 255 "%s is not a valid name for mass/width variable" % \ 256 value 257 258 if name is 'pdg_code': 259 if not isinstance(value, int): 260 raise self.PhysicsObjectError, \ 261 "PDG code %s is not an integer" % repr(value) 262 263 if name is 'line': 264 if not isinstance(value, str): 265 raise self.PhysicsObjectError, \ 266 "Line type %s is not a string" % repr(value) 267 if value not in ['dashed', 'straight', 'wavy', 'curly', 'double','swavy','scurly']: 268 raise self.PhysicsObjectError, \ 269 "Line type %s is unknown" % value 270 271 if name is 'charge': 272 if not isinstance(value, float): 273 raise self.PhysicsObjectError, \ 274 "Charge %s is not a float" % repr(value) 275 276 if name is 'propagating': 277 if not isinstance(value, bool): 278 raise self.PhysicsObjectError, \ 279 "Propagating tag %s is not a boolean" % repr(value) 280 281 if name in ['is_part', 'self_antipart']: 282 if not isinstance(value, bool): 283 raise self.PhysicsObjectError, \ 284 "%s tag %s is not a boolean" % (name, repr(value)) 285 286 return True
287
288 - def get_sorted_keys(self):
289 """Return particle property names as a nicely sorted list.""" 290 291 return self.sorted_keys
292 293 # Helper functions 294
295 - def get_pdg_code(self):
296 """Return the PDG code with a correct minus sign if the particle is its 297 own antiparticle""" 298 299 if not self['is_part'] and not self['self_antipart']: 300 return - self['pdg_code'] 301 else: 302 return self['pdg_code']
303
304 - def get_anti_pdg_code(self):
305 """Return the PDG code of the antiparticle with a correct minus sign 306 if the particle is its own antiparticle""" 307 308 if not self['self_antipart']: 309 return - self.get_pdg_code() 310 else: 311 return self['pdg_code']
312
313 - def get_color(self):
314 """Return the color code with a correct minus sign""" 315 316 if not self['is_part'] and abs(self['color']) in [3, 6]: 317 return - self['color'] 318 else: 319 return self['color']
320
321 - def get_anti_color(self):
322 """Return the color code of the antiparticle with a correct minus sign 323 """ 324 325 if self['is_part'] and self['color'] not in [1, 8]: 326 return - self['color'] 327 else: 328 return self['color']
329
330 - def get_name(self):
331 """Return the name if particle, antiname if antiparticle""" 332 333 if not self['is_part'] and not self['self_antipart']: 334 return self['antiname'] 335 else: 336 return self['name']
337
338 - def get_helicity_states(self):
339 """Return a list of the helicity states for the onshell particle""" 340 341 spin = self.get('spin') 342 if spin == 1: 343 # Scalar 344 return [ 0 ] 345 elif spin == 2: 346 # Spinor 347 return [ -1, 1 ] 348 elif spin == 3 and self.get('mass').lower() == 'zero': 349 # Massless vector 350 return [ -1, 1 ] 351 elif spin == 3: 352 # Massive vector 353 return [ -1, 0, 1 ] 354 elif spin == 5 and self.get('mass').lower() == 'zero': 355 # Massless tensor 356 return [-2, -1, 1, 2] 357 elif spin == 5: 358 # Massive tensor 359 return [-2, -1, 0, 1, 2] 360 361 raise self.PhysicsObjectError, \ 362 "No helicity state assignment for spin %d particles" % spin
363
364 - def is_fermion(self):
365 """Returns True if this is a fermion, False if boson""" 366 367 return self['spin'] % 2 == 0
368
369 - def is_boson(self):
370 """Returns True if this is a boson, False if fermion""" 371 372 return self['spin'] % 2 == 1
373
374 #=============================================================================== 375 # ParticleList 376 #=============================================================================== 377 -class ParticleList(PhysicsObjectList):
378 """A class to store lists of particles.""" 379
380 - def is_valid_element(self, obj):
381 """Test if object obj is a valid Particle for the list.""" 382 return isinstance(obj, Particle)
383 384
385 - def find_name(self, name):
386 """Try to find a particle with the given name. Check both name 387 and antiname. If a match is found, return the a copy of the 388 corresponding particle (first one in the list), with the 389 is_part flag set accordingly. None otherwise.""" 390 391 assert isinstance(name, str), "%s is not a valid string" % str(name) 392 393 for part in self: 394 mypart = copy.copy(part) 395 if part.get('name') == name: 396 mypart.set('is_part', True) 397 return mypart 398 elif part.get('antiname') == name: 399 mypart.set('is_part', False) 400 return mypart 401 402 return None
403
404 - def generate_ref_dict(self):
405 """Generate a dictionary of part/antipart pairs (as keys) and 406 0 (as value)""" 407 408 ref_dict_to0 = {} 409 410 for part in self: 411 ref_dict_to0[(part.get_pdg_code(), part.get_anti_pdg_code())] = [0] 412 ref_dict_to0[(part.get_anti_pdg_code(), part.get_pdg_code())] = [0] 413 414 return ref_dict_to0
415
416 - def generate_dict(self):
417 """Generate a dictionary from particle id to particle. 418 Include antiparticles. 419 """ 420 421 particle_dict = {} 422 423 for particle in self: 424 particle_dict[particle.get('pdg_code')] = particle 425 if not particle.get('self_antipart'): 426 antipart = copy.copy(particle) 427 antipart.set('is_part', False) 428 particle_dict[antipart.get_pdg_code()] = antipart 429 430 return particle_dict
431
432 433 #=============================================================================== 434 # Interaction 435 #=============================================================================== 436 -class Interaction(PhysicsObject):
437 """The interaction object containing the whole set of information 438 required to univocally characterize a given type of physical interaction: 439 440 particles: a list of particle ids 441 color: a list of string describing all the color structures involved 442 lorentz: a list of variable names describing all the Lorentz structure 443 involved 444 couplings: dictionary listing coupling variable names. The key is a 445 2-tuple of integers referring to color and Lorentz structures 446 orders: dictionary listing order names (as keys) with their value 447 """ 448 449 sorted_keys = ['id', 'particles', 'color', 'lorentz', 'couplings', 'orders'] 450
451 - def default_setup(self):
452 """Default values for all properties""" 453 454 self['id'] = 0 455 self['particles'] = [] 456 self['color'] = [] 457 self['lorentz'] = [] 458 self['couplings'] = { (0, 0):'none'} 459 self['orders'] = {}
460
461 - def filter(self, name, value):
462 """Filter for valid interaction property values.""" 463 464 if name == 'id': 465 #Should be an integer 466 if not isinstance(value, int): 467 raise self.PhysicsObjectError, \ 468 "%s is not a valid integer" % str(value) 469 470 if name == 'particles': 471 #Should be a list of valid particle names 472 if not isinstance(value, ParticleList): 473 raise self.PhysicsObjectError, \ 474 "%s is not a valid list of particles" % str(value) 475 476 if name == 'orders': 477 #Should be a dict with valid order names ask keys and int as values 478 if not isinstance(value, dict): 479 raise self.PhysicsObjectError, \ 480 "%s is not a valid dict for coupling orders" % \ 481 str(value) 482 for order in value.keys(): 483 if not isinstance(order, str): 484 raise self.PhysicsObjectError, \ 485 "%s is not a valid string" % str(order) 486 if not isinstance(value[order], int): 487 raise self.PhysicsObjectError, \ 488 "%s is not a valid integer" % str(value[order]) 489 490 if name in ['color']: 491 #Should be a list of list strings 492 if not isinstance(value, list): 493 raise self.PhysicsObjectError, \ 494 "%s is not a valid list of Color Strings" % str(value) 495 for mycolstring in value: 496 if not isinstance(mycolstring, color.ColorString): 497 raise self.PhysicsObjectError, \ 498 "%s is not a valid list of Color Strings" % str(value) 499 500 if name in ['lorentz']: 501 #Should be a list of list strings 502 if not isinstance(value, list): 503 raise self.PhysicsObjectError, \ 504 "%s is not a valid list of strings" % str(value) 505 for mystr in value: 506 if not isinstance(mystr, str): 507 raise self.PhysicsObjectError, \ 508 "%s is not a valid string" % str(mystr) 509 510 if name == 'couplings': 511 #Should be a dictionary of strings with (i,j) keys 512 if not isinstance(value, dict): 513 raise self.PhysicsObjectError, \ 514 "%s is not a valid dictionary for couplings" % \ 515 str(value) 516 517 for key in value.keys(): 518 if not isinstance(key, tuple): 519 raise self.PhysicsObjectError, \ 520 "%s is not a valid tuple" % str(key) 521 if len(key) != 2: 522 raise self.PhysicsObjectError, \ 523 "%s is not a valid tuple with 2 elements" % str(key) 524 if not isinstance(key[0], int) or not isinstance(key[1], int): 525 raise self.PhysicsObjectError, \ 526 "%s is not a valid tuple of integer" % str(key) 527 if not isinstance(value[key], str): 528 raise self.PhysicsObjectError, \ 529 "%s is not a valid string" % value[key] 530 531 return True
532
533 - def get_sorted_keys(self):
534 """Return particle property names as a nicely sorted list.""" 535 536 return self.sorted_keys
537 538
539 - def generate_dict_entries(self, ref_dict_to0, ref_dict_to1):
540 """Add entries corresponding to the current interactions to 541 the reference dictionaries (for n>0 and n-1>1)""" 542 543 # Create n>0 entries. Format is (p1,p2,p3,...):interaction_id. 544 # We are interested in the unordered list, so use sorted() 545 546 pdg_tuple = tuple(sorted([p.get_pdg_code() for p in self['particles']])) 547 if pdg_tuple not in ref_dict_to0.keys(): 548 ref_dict_to0[pdg_tuple] = [self['id']] 549 else: 550 ref_dict_to0[pdg_tuple].append(self['id']) 551 552 # Create n-1>1 entries. Note that, in the n-1 > 1 dictionary, 553 # the n-1 entries should have opposite sign as compared to 554 # interaction, since the interaction has outgoing particles, 555 # while in the dictionary we treat the n-1 particles as 556 # incoming 557 558 for part in self['particles']: 559 560 # We are interested in the unordered list, so use sorted() 561 pdg_tuple = tuple(sorted([p.get_pdg_code() for (i, p) in \ 562 enumerate(self['particles']) if \ 563 i != self['particles'].index(part)])) 564 pdg_part = part.get_anti_pdg_code() 565 if pdg_tuple in ref_dict_to1.keys(): 566 if (pdg_part, self['id']) not in ref_dict_to1[pdg_tuple]: 567 ref_dict_to1[pdg_tuple].append((pdg_part, self['id'])) 568 else: 569 ref_dict_to1[pdg_tuple] = [(pdg_part, self['id'])]
570
571 - def get_WEIGHTED_order(self, model):
572 """Get the WEIGHTED order for this interaction, for equivalent 573 3-particle vertex. Note that it can be fractional.""" 574 575 return float(sum([model.get('order_hierarchy')[key]*self.get('orders')[key]\ 576 for key in self.get('orders')]))/ \ 577 max((len(self.get('particles'))-2), 1)
578
579 - def __str__(self):
580 """String representation of an interaction. Outputs valid Python 581 with improved format. Overrides the PhysicsObject __str__ to only 582 display PDG code of involved particles.""" 583 584 mystr = '{\n' 585 586 for prop in self.get_sorted_keys(): 587 if isinstance(self[prop], str): 588 mystr = mystr + ' \'' + prop + '\': \'' + \ 589 self[prop] + '\',\n' 590 elif isinstance(self[prop], float): 591 mystr = mystr + ' \'' + prop + '\': %.2f,\n' % self[prop] 592 elif isinstance(self[prop], ParticleList): 593 mystr = mystr + ' \'' + prop + '\': [%s],\n' % \ 594 ','.join([str(part.get_pdg_code()) for part in self[prop]]) 595 else: 596 mystr = mystr + ' \'' + prop + '\': ' + \ 597 repr(self[prop]) + ',\n' 598 mystr = mystr.rstrip(',\n') 599 mystr = mystr + '\n}' 600 601 return mystr
602
603 #=============================================================================== 604 # InteractionList 605 #=============================================================================== 606 -class InteractionList(PhysicsObjectList):
607 """A class to store lists of interactionss.""" 608
609 - def is_valid_element(self, obj):
610 """Test if object obj is a valid Interaction for the list.""" 611 612 return isinstance(obj, Interaction)
613
614 - def generate_ref_dict(self):
615 """Generate the reference dictionaries from interaction list. 616 Return a list where the first element is the n>0 dictionary and 617 the second one is n-1>1.""" 618 619 ref_dict_to0 = {} 620 ref_dict_to1 = {} 621 622 for inter in self: 623 inter.generate_dict_entries(ref_dict_to0, ref_dict_to1) 624 625 return [ref_dict_to0, ref_dict_to1]
626
627 - def generate_dict(self):
628 """Generate a dictionary from interaction id to interaction. 629 """ 630 631 interaction_dict = {} 632 633 for inter in self: 634 interaction_dict[inter.get('id')] = inter 635 636 return interaction_dict
637
638 - def synchronize_interactions_with_particles(self, particle_dict):
639 """Make sure that the particles in the interactions are those 640 in the particle_dict, and that there are no interactions 641 refering to particles that don't exist. To be called when the 642 particle_dict is updated in a model. 643 """ 644 645 iint = 0 646 while iint < len(self): 647 inter = self[iint] 648 particles = inter.get('particles') 649 try: 650 for ipart, part in enumerate(particles): 651 particles[ipart] = particle_dict[part.get_pdg_code()] 652 iint += 1 653 except KeyError: 654 # This interaction has particles that no longer exist 655 self.pop(iint)
656
657 #=============================================================================== 658 # Model 659 #=============================================================================== 660 -class Model(PhysicsObject):
661 """A class to store all the model information.""" 662
663 - def default_setup(self):
664 665 self['name'] = "" 666 self['particles'] = ParticleList() 667 self['interactions'] = InteractionList() 668 self['parameters'] = None 669 self['functions'] = None 670 self['couplings'] = None 671 self['lorentz'] = None 672 self['particle_dict'] = {} 673 self['interaction_dict'] = {} 674 self['ref_dict_to0'] = {} 675 self['ref_dict_to1'] = {} 676 self['got_majoranas'] = None 677 self['order_hierarchy'] = {} 678 self['conserved_charge'] = set() 679 self['coupling_orders'] = None 680 self['expansion_order'] = None 681 self['version_tag'] = None # position of the directory (for security)
682
683 - def filter(self, name, value):
684 """Filter for model property values""" 685 686 if name == 'name': 687 if not isinstance(value, str): 688 raise self.PhysicsObjectError, \ 689 "Object of type %s is not a string" % \ 690 type(value) 691 elif name == 'particles': 692 if not isinstance(value, ParticleList): 693 raise self.PhysicsObjectError, \ 694 "Object of type %s is not a ParticleList object" % \ 695 type(value) 696 elif name == 'interactions': 697 if not isinstance(value, InteractionList): 698 raise self.PhysicsObjectError, \ 699 "Object of type %s is not a InteractionList object" % \ 700 type(value) 701 elif name == 'particle_dict': 702 if not isinstance(value, dict): 703 raise self.PhysicsObjectError, \ 704 "Object of type %s is not a dictionary" % \ 705 type(value) 706 elif name == 'interaction_dict': 707 if not isinstance(value, dict): 708 raise self.PhysicsObjectError, \ 709 "Object of type %s is not a dictionary" % type(value) 710 711 elif name == 'ref_dict_to0': 712 if not isinstance(value, dict): 713 raise self.PhysicsObjectError, \ 714 "Object of type %s is not a dictionary" % type(value) 715 716 elif name == 'ref_dict_to1': 717 if not isinstance(value, dict): 718 raise self.PhysicsObjectError, \ 719 "Object of type %s is not a dictionary" % type(value) 720 721 elif name == 'got_majoranas': 722 if not (isinstance(value, bool) or value == None): 723 raise self.PhysicsObjectError, \ 724 "Object of type %s is not a boolean" % type(value) 725 726 elif name == 'conserved_charge': 727 if not (isinstance(value, set)): 728 raise self.PhysicsObjectError, \ 729 "Object of type %s is not a set" % type(value) 730 731 elif name == 'version_tag': 732 if not (isinstance(value, str)): 733 raise self.PhysicsObjectError, \ 734 "Object of type %s is not a string" % type(value) 735 736 737 return True
738
739 - def get(self, name):
740 """Get the value of the property name.""" 741 742 if (name == 'ref_dict_to0' or name == 'ref_dict_to1') and \ 743 not self[name]: 744 if self['interactions']: 745 [self['ref_dict_to0'], self['ref_dict_to1']] = \ 746 self['interactions'].generate_ref_dict() 747 self['ref_dict_to0'].update( 748 self['particles'].generate_ref_dict()) 749 750 if (name == 'particle_dict') and not self[name]: 751 if self['particles']: 752 self['particle_dict'] = self['particles'].generate_dict() 753 if self['interactions']: 754 self['interactions'].synchronize_interactions_with_particles(\ 755 self['particle_dict']) 756 757 758 if (name == 'interaction_dict') and not self[name]: 759 if self['interactions']: 760 self['interaction_dict'] = self['interactions'].generate_dict() 761 762 if (name == 'got_majoranas') and self[name] == None: 763 if self['particles']: 764 self['got_majoranas'] = self.check_majoranas() 765 766 if (name == 'coupling_orders') and self[name] == None: 767 if self['interactions']: 768 self['coupling_orders'] = self.get_coupling_orders() 769 770 if (name == 'order_hierarchy') and not self[name]: 771 if self['interactions']: 772 self['order_hierarchy'] = self.get_order_hierarchy() 773 774 if (name == 'expansion_order') and self[name] == None: 775 if self['interactions']: 776 self['expansion_order'] = \ 777 dict([(order, -1) for order in self.get('coupling_orders')]) 778 779 return Model.__bases__[0].get(self, name) # call the mother routine
780
781 - def set(self, name, value):
782 """Special set for particles and interactions - need to 783 regenerate dictionaries.""" 784 785 if name == 'particles': 786 # Ensure no doublets in particle list 787 make_unique(value) 788 # Reset dictionaries 789 self['particle_dict'] = {} 790 self['ref_dict_to0'] = {} 791 self['got_majoranas'] = None 792 793 if name == 'interactions': 794 # Ensure no doublets in interaction list 795 make_unique(value) 796 # Reset dictionaries 797 self['interaction_dict'] = {} 798 self['ref_dict_to1'] = {} 799 self['ref_dict_to0'] = {} 800 self['got_majoranas'] = None 801 self['coupling_orders'] = None 802 self['order_hierarchy'] = {} 803 self['expansion_order'] = None 804 805 Model.__bases__[0].set(self, name, value) # call the mother routine 806 807 if name == 'particles': 808 # Recreate particle_dict 809 self.get('particle_dict')
810
811 - def get_sorted_keys(self):
812 """Return process property names as a nicely sorted list.""" 813 814 return ['name', 'particles', 'parameters', 'interactions', 'couplings', 815 'lorentz']
816
817 - def get_particle(self, id):
818 """Return the particle corresponding to the id""" 819 820 if id in self.get("particle_dict").keys(): 821 return self["particle_dict"][id] 822 else: 823 return None
824
825 - def get_interaction(self, id):
826 """Return the interaction corresponding to the id""" 827 828 829 if id in self.get("interaction_dict").keys(): 830 return self["interaction_dict"][id] 831 else: 832 return None
833
834 - def get_coupling_orders(self):
835 """Determine the coupling orders of the model""" 836 837 return set(sum([i.get('orders').keys() for i in \ 838 self.get('interactions')], []))
839
840 - def get_order_hierarchy(self):
841 """Set a default order hierarchy for the model if not set by the UFO.""" 842 # Set coupling hierachy 843 hierarchy = dict([(order, 1) for order in self.get('coupling_orders')]) 844 # Special case for only QCD and QED couplings, unless already set 845 if self.get('coupling_orders') == set(['QCD', 'QED']): 846 hierarchy['QED'] = 2 847 return hierarchy
848
849 - def get_particles_hierarchy(self):
850 """Returns the order hierarchies of the model and the 851 particles which have interactions in at least this hierarchy 852 (used in find_optimal_process_orders in MultiProcess diagram 853 generation): 854 855 Check the coupling hierarchy of the model. Assign all 856 particles to the different coupling hierarchies so that a 857 particle is considered to be in the highest hierarchy (i.e., 858 with lowest value) where it has an interaction. 859 """ 860 861 # Find coupling orders in model 862 coupling_orders = self.get('coupling_orders') 863 # Loop through the different coupling hierarchy values, so we 864 # start with the most dominant and proceed to the least dominant 865 hierarchy = sorted(list(set([self.get('order_hierarchy')[k] for \ 866 k in coupling_orders]))) 867 868 # orders is a rising list of the lists of orders with a given hierarchy 869 orders = [] 870 for value in hierarchy: 871 orders.append([ k for (k, v) in \ 872 self.get('order_hierarchy').items() if \ 873 v == value ]) 874 875 # Extract the interaction that correspond to the different 876 # coupling hierarchies, and the corresponding particles 877 interactions = [] 878 particles = [] 879 for iorder, order in enumerate(orders): 880 sum_orders = sum(orders[:iorder+1], []) 881 sum_interactions = sum(interactions[:iorder], []) 882 sum_particles = sum([list(p) for p in particles[:iorder]], []) 883 # Append all interactions that have only orders with at least 884 # this hierarchy 885 interactions.append([i for i in self.get('interactions') if \ 886 not i in sum_interactions and \ 887 not any([k not in sum_orders for k in \ 888 i.get('orders').keys()])]) 889 # Append the corresponding particles, excluding the 890 # particles that have already been added 891 particles.append(set(sum([[p.get_pdg_code() for p in \ 892 inter.get('particles') if \ 893 p.get_pdg_code() not in sum_particles] \ 894 for inter in interactions[-1]], []))) 895 896 return particles, hierarchy
897
898 - def get_max_WEIGHTED(self):
899 """Return the maximum WEIGHTED order for any interaction in the model, 900 for equivalent 3-particle vertices. Note that it can be fractional.""" 901 902 return max([inter.get_WEIGHTED_order(self) for inter in \ 903 self.get('interactions')])
904 905
906 - def check_majoranas(self):
907 """Return True if there is fermion flow violation, False otherwise""" 908 909 if any([part.is_fermion() and part.get('self_antipart') \ 910 for part in self.get('particles')]): 911 return True 912 913 # No Majorana particles, but may still be fermion flow 914 # violating interactions 915 for inter in self.get('interactions'): 916 fermions = [p for p in inter.get('particles') if p.is_fermion()] 917 for i in range(0, len(fermions), 2): 918 if fermions[i].get('is_part') == \ 919 fermions[i+1].get('is_part'): 920 # This is a fermion flow violating interaction 921 return True 922 # No fermion flow violations 923 return False
924
925 - def reset_dictionaries(self):
926 """Reset all dictionaries and got_majoranas. This is necessary 927 whenever the particle or interaction content has changed. If 928 particles or interactions are set using the set routine, this 929 is done automatically.""" 930 931 self['particle_dict'] = {} 932 self['ref_dict_to0'] = {} 933 self['got_majoranas'] = None 934 self['interaction_dict'] = {} 935 self['ref_dict_to1'] = {} 936 self['ref_dict_to0'] = {}
937
939 """Change the name of the particles such that all SM and MSSM particles 940 follows the MG convention""" 941 942 # Check that default name/antiname is not already use 943 def check_name_free(self, name): 944 """ check if name is not use for a particle in the model if it is 945 raise an MadGraph5error""" 946 part = self['particles'].find_name(name) 947 if part: 948 error_text = \ 949 '%s particles with pdg code %s is in conflict with MG ' + \ 950 'convention name for particle %s.\n Use -modelname in order ' + \ 951 'to use the particles name defined in the model and not the ' + \ 952 'MadGraph convention' 953 954 raise MadGraph5Error, error_text % \ 955 (part.get_name(), part.get_pdg_code(), pdg)
956 957 default = self.load_default_name() 958 959 for pdg in default.keys(): 960 part = self.get_particle(pdg) 961 if not part: 962 continue 963 antipart = self.get_particle(-pdg) 964 name = part.get_name() 965 if name != default[pdg]: 966 check_name_free(self, default[pdg]) 967 if part.get('is_part'): 968 part.set('name', default[pdg]) 969 if antipart: 970 antipart.set('name', default[pdg]) 971 else: 972 part.set('antiname', default[pdg]) 973 else: 974 part.set('antiname', default[pdg]) 975 if antipart: 976 antipart.set('antiname', default[pdg])
977
978 - def get_first_non_pdg(self):
979 """Return the first positive number that is not a valid PDG code""" 980 return [c for c in range(1, len(self.get('particles')) + 1) if \ 981 c not in self.get('particle_dict').keys()][0]
982
983 - def write_param_card(self):
984 """Write out the param_card, and return as string.""" 985 986 import models.write_param_card as writter 987 out = StringIO.StringIO() # it's suppose to be written in a file 988 param = writter.ParamCardWriter(self) 989 param.define_output_file(out) 990 param.write_card() 991 return out.getvalue()
992 993 @ staticmethod
994 - def load_default_name():
995 """ load the default for name convention """ 996 997 logger.info('Change particles name to pass to MG5 convention') 998 default = {} 999 for line in open(os.path.join(MG5DIR, 'input', \ 1000 'particles_name_default.txt')): 1001 line = line.lstrip() 1002 if line.startswith('#'): 1003 continue 1004 1005 args = line.split() 1006 if len(args) != 2: 1007 logger.warning('Invalid syntax in interface/default_name:\n %s' % line) 1008 continue 1009 default[int(args[0])] = args[1].lower() 1010 1011 return default
1012
1013 ################################################################################ 1014 # Class for Parameter / Coupling 1015 ################################################################################ 1016 -class ModelVariable(object):
1017 """A Class for storing the information about coupling/ parameter""" 1018
1019 - def __init__(self, name, expression, type, depend=()):
1020 """Initialize a new parameter/coupling""" 1021 1022 self.name = name 1023 self.expr = expression # python expression 1024 self.type = type # real/complex 1025 self.depend = depend # depend on some other parameter -tuple- 1026 self.value = None
1027
1028 - def __eq__(self, other):
1029 """Object with same name are identical, If the object is a string we check 1030 if the attribute name is equal to this string""" 1031 1032 try: 1033 return other.name == self.name 1034 except: 1035 return other == self.name
1036
1037 -class ParamCardVariable(ModelVariable):
1038 """ A class for storing the information linked to all the parameter 1039 which should be define in the param_card.dat""" 1040 1041 depend = ('external',) 1042 type = 'real' 1043
1044 - def __init__(self, name, value, lhablock, lhacode):
1045 """Initialize a new ParamCardVariable 1046 name: name of the variable 1047 value: default numerical value 1048 lhablock: name of the block in the param_card.dat 1049 lhacode: code associate to the variable 1050 """ 1051 self.name = name 1052 self.value = value 1053 self.lhablock = lhablock 1054 self.lhacode = lhacode
1055
1056 1057 #=============================================================================== 1058 # Classes used in diagram generation and process definition: 1059 # Leg, Vertex, Diagram, Process 1060 #=============================================================================== 1061 1062 #=============================================================================== 1063 # Leg 1064 #=============================================================================== 1065 -class Leg(PhysicsObject):
1066 """Leg object: id (Particle), number, I/F state, flag from_group 1067 """ 1068
1069 - def default_setup(self):
1070 """Default values for all properties""" 1071 1072 self['id'] = 0 1073 self['number'] = 0 1074 # state: True = final, False = initial (boolean to save memory) 1075 self['state'] = True 1076 # from_group: Used in diagram generation 1077 self['from_group'] = True 1078 # onshell: decaying leg (True), forbidden s-channel (False), none (None) 1079 self['onshell'] = None
1080
1081 - def filter(self, name, value):
1082 """Filter for valid leg property values.""" 1083 1084 if name in ['id', 'number']: 1085 if not isinstance(value, int): 1086 raise self.PhysicsObjectError, \ 1087 "%s is not a valid integer for leg id" % str(value) 1088 1089 if name == 'state': 1090 if not isinstance(value, bool): 1091 raise self.PhysicsObjectError, \ 1092 "%s is not a valid leg state (True|False)" % \ 1093 str(value) 1094 1095 if name == 'from_group': 1096 if not isinstance(value, bool) and value != None: 1097 raise self.PhysicsObjectError, \ 1098 "%s is not a valid boolean for leg flag from_group" % \ 1099 str(value) 1100 1101 if name == 'onshell': 1102 if not isinstance(value, bool) and value != None: 1103 raise self.PhysicsObjectError, \ 1104 "%s is not a valid boolean for leg flag onshell" % \ 1105 str(value) 1106 1107 return True
1108
1109 - def get_sorted_keys(self):
1110 """Return particle property names as a nicely sorted list.""" 1111 1112 return ['id', 'number', 'state', 'from_group', 'onshell']
1113
1114 - def is_fermion(self, model):
1115 """Returns True if the particle corresponding to the leg is a 1116 fermion""" 1117 1118 assert isinstance(model, Model), "%s is not a model" % str(model) 1119 1120 return model.get('particle_dict')[self['id']].is_fermion()
1121
1122 - def is_incoming_fermion(self, model):
1123 """Returns True if leg is an incoming fermion, i.e., initial 1124 particle or final antiparticle""" 1125 1126 assert isinstance(model, Model), "%s is not a model" % str(model) 1127 1128 part = model.get('particle_dict')[self['id']] 1129 return part.is_fermion() and \ 1130 (self.get('state') == False and part.get('is_part') or \ 1131 self.get('state') == True and not part.get('is_part'))
1132
1133 - def is_outgoing_fermion(self, model):
1134 """Returns True if leg is an outgoing fermion, i.e., initial 1135 antiparticle or final particle""" 1136 1137 assert isinstance(model, Model), "%s is not a model" % str(model) 1138 1139 part = model.get('particle_dict')[self['id']] 1140 return part.is_fermion() and \ 1141 (self.get('state') == True and part.get('is_part') or \ 1142 self.get('state') == False and not part.get('is_part'))
1143 1144 # Make sure sort() sorts lists of legs according to 'number'
1145 - def __lt__(self, other):
1146 return self['number'] < other['number']
1147
1148 #=============================================================================== 1149 # LegList 1150 #=============================================================================== 1151 -class LegList(PhysicsObjectList):
1152 """List of Leg objects 1153 """ 1154
1155 - def is_valid_element(self, obj):
1156 """Test if object obj is a valid Leg for the list.""" 1157 1158 return isinstance(obj, Leg)
1159 1160 # Helper methods for diagram generation 1161
1162 - def from_group_elements(self):
1163 """Return all elements which have 'from_group' True""" 1164 1165 return filter(lambda leg: leg.get('from_group'), self)
1166
1167 - def minimum_one_from_group(self):
1168 """Return True if at least one element has 'from_group' True""" 1169 1170 return len(self.from_group_elements()) > 0
1171
1172 - def minimum_two_from_group(self):
1173 """Return True if at least two elements have 'from_group' True""" 1174 1175 return len(self.from_group_elements()) > 1
1176
1177 - def can_combine_to_1(self, ref_dict_to1):
1178 """If has at least one 'from_group' True and in ref_dict_to1, 1179 return the return list from ref_dict_to1, otherwise return False""" 1180 if self.minimum_one_from_group(): 1181 return ref_dict_to1.has_key(tuple(sorted([leg.get('id') for leg in self]))) 1182 else: 1183 return False
1184
1185 - def can_combine_to_0(self, ref_dict_to0, is_decay_chain=False):
1186 """If has at least two 'from_group' True and in ref_dict_to0, 1187 1188 return the vertex (with id from ref_dict_to0), otherwise return None 1189 1190 If is_decay_chain = True, we only allow clustering of the 1191 initial leg, since we want this to be the last wavefunction to 1192 be evaluated. 1193 """ 1194 if is_decay_chain: 1195 # Special treatment - here we only allow combination to 0 1196 # if the initial leg (marked by from_group = None) is 1197 # unclustered, since we want this to stay until the very 1198 # end. 1199 return any(leg.get('from_group') == None for leg in self) and \ 1200 ref_dict_to0.has_key(tuple(sorted([leg.get('id') \ 1201 for leg in self]))) 1202 1203 if self.minimum_two_from_group(): 1204 return ref_dict_to0.has_key(tuple(sorted([leg.get('id') for leg in self]))) 1205 else: 1206 return False
1207
1208 - def get_outgoing_id_list(self, model):
1209 """Returns the list of ids corresponding to the leglist with 1210 all particles outgoing""" 1211 1212 res = [] 1213 1214 assert isinstance(model, Model), "Error! model not model" 1215 1216 1217 for leg in self: 1218 if leg.get('state') == False: 1219 res.append(model.get('particle_dict')[leg.get('id')].get_anti_pdg_code()) 1220 else: 1221 res.append(leg.get('id')) 1222 1223 return res
1224
1225 1226 #=============================================================================== 1227 # MultiLeg 1228 #=============================================================================== 1229 -class MultiLeg(PhysicsObject):
1230 """MultiLeg object: ids (Particle or particles), I/F state 1231 """ 1232
1233 - def default_setup(self):
1234 """Default values for all properties""" 1235 1236 self['ids'] = [] 1237 self['state'] = True
1238
1239 - def filter(self, name, value):
1240 """Filter for valid multileg property values.""" 1241 1242 if name == 'ids': 1243 if not isinstance(value, list): 1244 raise self.PhysicsObjectError, \ 1245 "%s is not a valid list" % str(value) 1246 for i in value: 1247 if not isinstance(i, int): 1248 raise self.PhysicsObjectError, \ 1249 "%s is not a valid list of integers" % str(value) 1250 1251 if name == 'state': 1252 if not isinstance(value, bool): 1253 raise self.PhysicsObjectError, \ 1254 "%s is not a valid leg state (initial|final)" % \ 1255 str(value) 1256 1257 return True
1258
1259 - def get_sorted_keys(self):
1260 """Return particle property names as a nicely sorted list.""" 1261 1262 return ['ids', 'state']
1263
1264 #=============================================================================== 1265 # LegList 1266 #=============================================================================== 1267 -class MultiLegList(PhysicsObjectList):
1268 """List of MultiLeg objects 1269 """ 1270
1271 - def is_valid_element(self, obj):
1272 """Test if object obj is a valid MultiLeg for the list.""" 1273 1274 return isinstance(obj, MultiLeg)
1275
1276 #=============================================================================== 1277 # Vertex 1278 #=============================================================================== 1279 -class Vertex(PhysicsObject):
1280 """Vertex: list of legs (ordered), id (Interaction) 1281 """ 1282 1283 sorted_keys = ['id', 'legs'] 1284
1285 - def default_setup(self):
1286 """Default values for all properties""" 1287 1288 self['id'] = 0 1289 self['legs'] = LegList()
1290
1291 - def filter(self, name, value):
1292 """Filter for valid vertex property values.""" 1293 1294 if name == 'id': 1295 if not isinstance(value, int): 1296 raise self.PhysicsObjectError, \ 1297 "%s is not a valid integer for vertex id" % str(value) 1298 1299 if name == 'legs': 1300 if not isinstance(value, LegList): 1301 raise self.PhysicsObjectError, \ 1302 "%s is not a valid LegList object" % str(value) 1303 1304 return True
1305
1306 - def get_sorted_keys(self):
1307 """Return particle property names as a nicely sorted list.""" 1308 1309 return self.sorted_keys #['id', 'legs']
1310
1311 - def get_s_channel_id(self, model, ninitial):
1312 """Returns the id for the last leg as an outgoing 1313 s-channel. Returns 0 if leg is t-channel, or if identity 1314 vertex. Used to check for required and forbidden s-channel 1315 particles.""" 1316 1317 leg = self.get('legs')[-1] 1318 1319 if ninitial == 1: 1320 # For one initial particle, all legs are s-channel 1321 # Only need to flip particle id if state is False 1322 if leg.get('state') == True: 1323 return leg.get('id') 1324 else: 1325 return model.get('particle_dict')[leg.get('id')].\ 1326 get_anti_pdg_code() 1327 1328 # Number of initial particles is at least 2 1329 if self.get('id') == 0 or \ 1330 leg.get('state') == False: 1331 # identity vertex or t-channel particle 1332 return 0 1333 1334 # Check if the particle number is <= ninitial 1335 # In that case it comes from initial and we should switch direction 1336 if leg.get('number') > ninitial: 1337 return leg.get('id') 1338 else: 1339 return model.get('particle_dict')[leg.get('id')].\ 1340 get_anti_pdg_code()
1341
1342 ## Check if the other legs are initial or final. 1343 ## If the latter, return leg id, if the former, return -leg id 1344 #if self.get('legs')[0].get('state') == True: 1345 # return leg.get('id') 1346 #else: 1347 # return model.get('particle_dict')[leg.get('id')].\ 1348 # get_anti_pdg_code() 1349 1350 #=============================================================================== 1351 # VertexList 1352 #=============================================================================== 1353 -class VertexList(PhysicsObjectList):
1354 """List of Vertex objects 1355 """ 1356 1357 orders = {} 1358
1359 - def is_valid_element(self, obj):
1360 """Test if object obj is a valid Vertex for the list.""" 1361 1362 return isinstance(obj, Vertex)
1363
1364 - def __init__(self, init_list=None, orders=None):
1365 """Creates a new list object, with an optional dictionary of 1366 coupling orders.""" 1367 1368 list.__init__(self) 1369 1370 if init_list is not None: 1371 for object in init_list: 1372 self.append(object) 1373 1374 if isinstance(orders, dict): 1375 self.orders = orders
1376
1377 1378 #=============================================================================== 1379 # Diagram 1380 #=============================================================================== 1381 -class Diagram(PhysicsObject):
1382 """Diagram: list of vertices (ordered) 1383 """ 1384
1385 - def default_setup(self):
1386 """Default values for all properties""" 1387 1388 self['vertices'] = VertexList() 1389 self['orders'] = {}
1390
1391 - def filter(self, name, value):
1392 """Filter for valid diagram property values.""" 1393 1394 if name == 'vertices': 1395 if not isinstance(value, VertexList): 1396 raise self.PhysicsObjectError, \ 1397 "%s is not a valid VertexList object" % str(value) 1398 1399 if name == 'orders': 1400 Interaction.filter(Interaction(), 'orders', value) 1401 1402 return True
1403
1404 - def get_sorted_keys(self):
1405 """Return particle property names as a nicely sorted list.""" 1406 1407 return ['vertices', 'orders']
1408
1409 - def nice_string(self):
1410 """Returns a nicely formatted string of the diagram content.""" 1411 1412 if self['vertices']: 1413 mystr = '(' 1414 for vert in self['vertices']: 1415 mystr = mystr + '(' 1416 for leg in vert['legs'][:-1]: 1417 mystr = mystr + str(leg['number']) + '(%s)' % str(leg['id']) + ',' 1418 1419 if self['vertices'].index(vert) < len(self['vertices']) - 1: 1420 # Do not want ">" in the last vertex 1421 mystr = mystr[:-1] + '>' 1422 mystr = mystr + str(vert['legs'][-1]['number']) + '(%s)' % str(vert['legs'][-1]['id']) + ',' 1423 mystr = mystr + 'id:' + str(vert['id']) + '),' 1424 mystr = mystr[:-1] + ')' 1425 mystr += " (%s)" % ",".join(["%s=%d" % (key, self['orders'][key]) \ 1426 for key in sorted(self['orders'].keys())]) 1427 return mystr 1428 else: 1429 return '()'
1430
1431 - def calculate_orders(self, model):
1432 """Calculate the actual coupling orders of this diagram. Note 1433 that the special order WEIGTHED corresponds to the sum of 1434 hierarchys for the couplings.""" 1435 1436 coupling_orders = dict([(c, 0) for c in model.get('coupling_orders')]) 1437 weight = 0 1438 for vertex in self['vertices']: 1439 if vertex.get('id') == 0: continue 1440 couplings = model.get('interaction_dict')[vertex.get('id')].\ 1441 get('orders') 1442 for coupling in couplings: 1443 coupling_orders[coupling] += couplings[coupling] 1444 weight += sum([model.get('order_hierarchy')[c]*n for \ 1445 (c,n) in couplings.items()]) 1446 coupling_orders['WEIGHTED'] = weight 1447 self.set('orders', coupling_orders)
1448
1449 - def renumber_legs(self, perm_map, leg_list):
1450 """Renumber legs in all vertices according to perm_map""" 1451 vertices = VertexList() 1452 min_dict = copy.copy(perm_map) 1453 # Dictionary from leg number to state 1454 state_dict = dict([(l.get('number'), l.get('state')) for l in leg_list]) 1455 # First renumber all legs in the n-1->1 vertices 1456 for vertex in self.get('vertices')[:-1]: 1457 vertex = copy.copy(vertex) 1458 leg_list = LegList([copy.copy(l) for l in vertex.get('legs')]) 1459 for leg in leg_list[:-1]: 1460 leg.set('number', min_dict[leg.get('number')]) 1461 leg.set('state', state_dict[leg.get('number')]) 1462 min_number = min([leg.get('number') for leg in leg_list[:-1]]) 1463 leg = leg_list[-1] 1464 min_dict[leg.get('number')] = min_number 1465 # resulting leg is initial state if there is exactly one 1466 # initial state leg among the incoming legs 1467 state_dict[min_number] = len([l for l in leg_list[:-1] if \ 1468 not l.get('state')]) != 1 1469 leg.set('number', min_number) 1470 leg.set('state', state_dict[min_number]) 1471 vertex.set('legs', leg_list) 1472 vertices.append(vertex) 1473 # Now renumber the legs in final vertex 1474 vertex = copy.copy(self.get('vertices')[-1]) 1475 leg_list = LegList([copy.copy(l) for l in vertex.get('legs')]) 1476 for leg in leg_list: 1477 leg.set('number', min_dict[leg.get('number')]) 1478 leg.set('state', state_dict[leg.get('number')]) 1479 vertex.set('legs', leg_list) 1480 vertices.append(vertex) 1481 # Finally create new diagram 1482 new_diag = copy.copy(self) 1483 new_diag.set('vertices', vertices) 1484 state_dict = {True:'T',False:'F'} 1485 return new_diag
1486
1487 - def get_vertex_leg_numbers(self):
1488 """Return a list of the number of legs in the vertices for 1489 this diagram""" 1490 1491 return [len(v.get('legs')) for v in self.get('vertices')]
1492
1493 - def get_num_configs(self, model, ninitial):
1494 """Return the maximum number of configs from this diagram, 1495 given by 2^(number of non-zero width s-channel propagators)""" 1496 1497 s_channels = [v.get_s_channel_id(model,ninitial) for v in \ 1498 self.get('vertices')[:-1]] 1499 num_props = len([i for i in s_channels if i != 0 and \ 1500 model.get_particle(i).get('width').lower() != 'zero']) 1501 1502 if num_props < 1: 1503 return 1 1504 else: 1505 return 2**num_props
1506 #===============================================================================
1507 # DiagramList 1508 #=============================================================================== 1509 -class DiagramList(PhysicsObjectList):
1510 """List of Diagram objects 1511 """ 1512
1513 - def is_valid_element(self, obj):
1514 """Test if object obj is a valid Diagram for the list.""" 1515 1516 return isinstance(obj, Diagram)
1517
1518 - def nice_string(self, indent=0):
1519 """Returns a nicely formatted string""" 1520 mystr = " " * indent + str(len(self)) + ' diagrams:\n' 1521 for i, diag in enumerate(self): 1522 mystr = mystr + " " * indent + str(i+1) + " " + \ 1523 diag.nice_string() + '\n' 1524 return mystr[:-1]
1525
1526 1527 #=============================================================================== 1528 # Process 1529 #=============================================================================== 1530 -class Process(PhysicsObject):
1531 """Process: list of legs (ordered) 1532 dictionary of orders 1533 model 1534 process id 1535 """ 1536
1537 - def default_setup(self):
1538 """Default values for all properties""" 1539 1540 self['legs'] = LegList() 1541 self['orders'] = {} 1542 self['model'] = Model() 1543 # Optional number to identify the process 1544 self['id'] = 0 1545 self['uid'] = 0 # should be a uniq id number 1546 # Required s-channels are given as a list of id lists. Only 1547 # diagrams with all s-channels in any of the lists are 1548 # allowed. This enables generating e.g. Z/gamma as s-channel 1549 # propagators. 1550 self['required_s_channels'] = [] 1551 self['forbidden_onsh_s_channels'] = [] 1552 self['forbidden_s_channels'] = [] 1553 self['forbidden_particles'] = [] 1554 self['is_decay_chain'] = False 1555 self['overall_orders'] = {} 1556 # Decay chain processes associated with this process 1557 self['decay_chains'] = ProcessList()
1558
1559 - def filter(self, name, value):
1560 """Filter for valid process property values.""" 1561 1562 if name == 'legs': 1563 if not isinstance(value, LegList): 1564 raise self.PhysicsObjectError, \ 1565 "%s is not a valid LegList object" % str(value) 1566 1567 if name in ['orders', 'overall_orders']: 1568 Interaction.filter(Interaction(), 'orders', value) 1569 1570 if name == 'model': 1571 if not isinstance(value, Model): 1572 raise self.PhysicsObjectError, \ 1573 "%s is not a valid Model object" % str(value) 1574 if name in ['id', 'uid']: 1575 if not isinstance(value, int): 1576 raise self.PhysicsObjectError, \ 1577 "Process %s %s is not an integer" % (name, repr(value)) 1578 1579 if name == 'required_s_channels': 1580 if not isinstance(value, list): 1581 raise self.PhysicsObjectError, \ 1582 "%s is not a valid list" % str(value) 1583 for l in value: 1584 if not isinstance(l, list): 1585 raise self.PhysicsObjectError, \ 1586 "%s is not a valid list of lists" % str(value) 1587 for i in l: 1588 if not isinstance(i, int): 1589 raise self.PhysicsObjectError, \ 1590 "%s is not a valid list of integers" % str(l) 1591 if i == 0: 1592 raise self.PhysicsObjectError, \ 1593 "Not valid PDG code %d for s-channel particle" % i 1594 1595 if name in ['forbidden_onsh_s_channels', 'forbidden_s_channels']: 1596 if not isinstance(value, list): 1597 raise self.PhysicsObjectError, \ 1598 "%s is not a valid list" % str(value) 1599 for i in value: 1600 if not isinstance(i, int): 1601 raise self.PhysicsObjectError, \ 1602 "%s is not a valid list of integers" % str(value) 1603 if i == 0: 1604 raise self.PhysicsObjectError, \ 1605 "Not valid PDG code %d for s-channel particle" % str(value) 1606 1607 if name == 'forbidden_particles': 1608 if not isinstance(value, list): 1609 raise self.PhysicsObjectError, \ 1610 "%s is not a valid list" % str(value) 1611 for i in value: 1612 if not isinstance(i, int): 1613 raise self.PhysicsObjectError, \ 1614 "%s is not a valid list of integers" % str(value) 1615 if i <= 0: 1616 raise self.PhysicsObjectError, \ 1617 "Forbidden particles should have a positive PDG code" % str(value) 1618 1619 if name == 'is_decay_chain': 1620 if not isinstance(value, bool): 1621 raise self.PhysicsObjectError, \ 1622 "%s is not a valid bool" % str(value) 1623 1624 if name == 'decay_chains': 1625 if not isinstance(value, ProcessList): 1626 raise self.PhysicsObjectError, \ 1627 "%s is not a valid ProcessList" % str(value) 1628 1629 return True
1630
1631 - def set(self, name, value):
1632 """Special set for forbidden particles - set to abs value.""" 1633 1634 if name == 'forbidden_particles': 1635 try: 1636 value = [abs(i) for i in value] 1637 except: 1638 pass 1639 1640 if name == 'required_s_channels': 1641 # Required s-channels need to be a list of lists of ids 1642 if value and isinstance(value, list) and \ 1643 not isinstance(value[0], list): 1644 value = [value] 1645 1646 return super(Process, self).set(name, value) # call the mother routine
1647
1648 - def get_sorted_keys(self):
1649 """Return process property names as a nicely sorted list.""" 1650 1651 return ['legs', 'orders', 'overall_orders', 'model', 'id', 1652 'required_s_channels', 'forbidden_onsh_s_channels', 1653 'forbidden_s_channels', 1654 'forbidden_particles', 'is_decay_chain', 'decay_chains']
1655
1656 - def nice_string(self, indent=0):
1657 """Returns a nicely formated string about current process 1658 content""" 1659 1660 mystr = " " * indent + "Process: " 1661 prevleg = None 1662 for leg in self['legs']: 1663 mypart = self['model'].get('particle_dict')[leg['id']] 1664 if prevleg and prevleg['state'] == False \ 1665 and leg['state'] == True: 1666 # Separate initial and final legs by ">" 1667 mystr = mystr + '> ' 1668 # Add required s-channels 1669 if self['required_s_channels'] and \ 1670 self['required_s_channels'][0]: 1671 mystr += "|".join([" ".join([self['model'].\ 1672 get('particle_dict')[req_id].get_name() \ 1673 for req_id in id_list]) \ 1674 for id_list in self['required_s_channels']]) 1675 mystr = mystr + ' > ' 1676 1677 mystr = mystr + mypart.get_name() + ' ' 1678 #mystr = mystr + '(%i) ' % leg['number'] 1679 prevleg = leg 1680 1681 # Add forbidden s-channels 1682 if self['forbidden_onsh_s_channels']: 1683 mystr = mystr + '$ ' 1684 for forb_id in self['forbidden_onsh_s_channels']: 1685 forbpart = self['model'].get('particle_dict')[forb_id] 1686 mystr = mystr + forbpart.get_name() + ' ' 1687 1688 # Add double forbidden s-channels 1689 if self['forbidden_s_channels']: 1690 mystr = mystr + '$$ ' 1691 for forb_id in self['forbidden_s_channels']: 1692 forbpart = self['model'].get('particle_dict')[forb_id] 1693 mystr = mystr + forbpart.get_name() + ' ' 1694 1695 # Add forbidden particles 1696 if self['forbidden_particles']: 1697 mystr = mystr + '/ ' 1698 for forb_id in self['forbidden_particles']: 1699 forbpart = self['model'].get('particle_dict')[forb_id] 1700 mystr = mystr + forbpart.get_name() + ' ' 1701 1702 if self['orders']: 1703 mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \ 1704 for key in sorted(self['orders'])]) + ' ' 1705 1706 # Remove last space 1707 mystr = mystr[:-1] 1708 1709 if self.get('id') or self.get('overall_orders'): 1710 mystr += " @%d" % self.get('id') 1711 if self.get('overall_orders'): 1712 mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \ 1713 for key in sorted(self['orders'])]) + ' ' 1714 1715 if not self.get('decay_chains'): 1716 return mystr 1717 1718 for decay in self['decay_chains']: 1719 mystr = mystr + '\n' + \ 1720 decay.nice_string(indent + 2).replace('Process', 'Decay') 1721 1722 return mystr
1723
1724 - def input_string(self):
1725 """Returns a process string corresponding to the input string 1726 in the command line interface.""" 1727 1728 mystr = "" 1729 prevleg = None 1730 1731 for leg in self['legs']: 1732 mypart = self['model'].get('particle_dict')[leg['id']] 1733 if prevleg and prevleg['state'] == False \ 1734 and leg['state'] == True: 1735 # Separate initial and final legs by ">" 1736 mystr = mystr + '> ' 1737 # Add required s-channels 1738 if self['required_s_channels'] and \ 1739 self['required_s_channels'][0]: 1740 mystr += "|".join([" ".join([self['model'].\ 1741 get('particle_dict')[req_id].get_name() \ 1742 for req_id in id_list]) \ 1743 for id_list in self['required_s_channels']]) 1744 mystr = mystr + '> ' 1745 1746 mystr = mystr + mypart.get_name() + ' ' 1747 #mystr = mystr + '(%i) ' % leg['number'] 1748 prevleg = leg 1749 1750 # Add forbidden s-channels 1751 if self['forbidden_onsh_s_channels']: 1752 mystr = mystr + '$ ' 1753 for forb_id in self['forbidden_onsh_s_channels']: 1754 forbpart = self['model'].get('particle_dict')[forb_id] 1755 mystr = mystr + forbpart.get_name() + ' ' 1756 1757 # Add double forbidden s-channels 1758 if self['forbidden_s_channels']: 1759 mystr = mystr + '$$ ' 1760 for forb_id in self['forbidden_s_channels']: 1761 forbpart = self['model'].get('particle_dict')[forb_id] 1762 mystr = mystr + forbpart.get_name() + ' ' 1763 1764 # Add forbidden particles 1765 if self['forbidden_particles']: 1766 mystr = mystr + '/ ' 1767 for forb_id in self['forbidden_particles']: 1768 forbpart = self['model'].get('particle_dict')[forb_id] 1769 mystr = mystr + forbpart.get_name() + ' ' 1770 1771 if self['orders']: 1772 mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \ 1773 for key in sorted(self['orders'])]) + ' ' 1774 1775 # Remove last space 1776 mystr = mystr[:-1] 1777 1778 if self.get('overall_orders'): 1779 mystr += " @%d" % self.get('id') 1780 if self.get('overall_orders'): 1781 mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \ 1782 for key in sorted(self['orders'])]) + ' ' 1783 1784 if not self.get('decay_chains'): 1785 return mystr 1786 1787 for decay in self['decay_chains']: 1788 paren1 = '' 1789 paren2 = '' 1790 if decay.get('decay_chains'): 1791 paren1 = '(' 1792 paren2 = ')' 1793 mystr += ', ' + paren1 + decay.input_string() + paren2 1794 1795 return mystr
1796
1797 - def base_string(self):
1798 """Returns a string containing only the basic process (w/o decays).""" 1799 1800 mystr = "" 1801 prevleg = None 1802 for leg in self.get_legs_with_decays(): 1803 mypart = self['model'].get('particle_dict')[leg['id']] 1804 if prevleg and prevleg['state'] == False \ 1805 and leg['state'] == True: 1806 # Separate initial and final legs by ">" 1807 mystr = mystr + '> ' 1808 mystr = mystr + mypart.get_name() + ' ' 1809 prevleg = leg 1810 1811 # Remove last space 1812 return mystr[:-1]
1813
1814 - def shell_string(self, schannel=True, forbid=True, main=True):
1815 """Returns process as string with '~' -> 'x', '>' -> '_', 1816 '+' -> 'p' and '-' -> 'm', including process number, 1817 intermediate s-channels and forbidden particles""" 1818 1819 mystr = "" 1820 if not self.get('is_decay_chain'): 1821 mystr += "%d_" % self['id'] 1822 1823 prevleg = None 1824 for leg in self['legs']: 1825 mypart = self['model'].get('particle_dict')[leg['id']] 1826 if prevleg and prevleg['state'] == False \ 1827 and leg['state'] == True: 1828 # Separate initial and final legs by ">" 1829 mystr = mystr + '_' 1830 # Add required s-channels 1831 if self['required_s_channels'] and \ 1832 self['required_s_channels'][0] and schannel: 1833 mystr += "_or_".join(["".join([self['model'].\ 1834 get('particle_dict')[req_id].get_name() \ 1835 for req_id in id_list]) \ 1836 for id_list in self['required_s_channels']]) 1837 mystr = mystr + '_' 1838 if mypart['is_part']: 1839 mystr = mystr + mypart['name'] 1840 else: 1841 mystr = mystr + mypart['antiname'] 1842 prevleg = leg 1843 1844 # Check for forbidden particles 1845 if self['forbidden_particles'] and forbid: 1846 mystr = mystr + '_no_' 1847 for forb_id in self['forbidden_particles']: 1848 forbpart = self['model'].get('particle_dict')[forb_id] 1849 mystr = mystr + forbpart.get_name() 1850 1851 # Replace '~' with 'x' 1852 mystr = mystr.replace('~', 'x') 1853 # Replace '+' with 'p' 1854 mystr = mystr.replace('+', 'p') 1855 # Replace '-' with 'm' 1856 mystr = mystr.replace('-', 'm') 1857 # Just to be safe, remove all spaces 1858 mystr = mystr.replace(' ', '') 1859 1860 for decay in self.get('decay_chains'): 1861 mystr = mystr + "_" + decay.shell_string(schannel,forbid, main=False) 1862 1863 # Too long name are problematic so restrict them to a maximal of 70 char 1864 if len(mystr) > 64 and main: 1865 if schannel and forbid: 1866 return self.shell_string(True, False, False)+ '_%s' % self['uid'] 1867 elif schannel: 1868 return self.shell_string(False, False, False)+'_%s' % self['uid'] 1869 else: 1870 return mystr[:64]+'_%s' % self['uid'] 1871 1872 1873 1874 1875 return mystr
1876
1877 - def shell_string_v4(self):
1878 """Returns process as v4-compliant string with '~' -> 'x' and 1879 '>' -> '_'""" 1880 1881 mystr = "%d_" % self['id'] 1882 prevleg = None 1883 for leg in self.get_legs_with_decays(): 1884 mypart = self['model'].get('particle_dict')[leg['id']] 1885 if prevleg and prevleg['state'] == False \ 1886 and leg['state'] == True: 1887 # Separate initial and final legs by ">" 1888 mystr = mystr + '_' 1889 if mypart['is_part']: 1890 mystr = mystr + mypart['name'] 1891 else: 1892 mystr = mystr + mypart['antiname'] 1893 prevleg = leg 1894 1895 # Replace '~' with 'x' 1896 mystr = mystr.replace('~', 'x') 1897 # Just to be safe, remove all spaces 1898 mystr = mystr.replace(' ', '') 1899 1900 return mystr
1901 1902 # Helper functions 1903
1904 - def get_ninitial(self):
1905 """Gives number of initial state particles""" 1906 1907 return len(filter(lambda leg: leg.get('state') == False, 1908 self.get('legs')))
1909
1910 - def get_initial_ids(self):
1911 """Gives the pdg codes for initial state particles""" 1912 1913 return [leg.get('id') for leg in \ 1914 filter(lambda leg: leg.get('state') == False, 1915 self.get('legs'))]
1916
1917 - def get_initial_pdg(self, number):
1918 """Return the pdg codes for initial state particles for beam number""" 1919 1920 return filter(lambda leg: leg.get('state') == False and\ 1921 leg.get('number') == number, 1922 self.get('legs'))[0].get('id')
1923
1924 - def get_final_legs(self):
1925 """Gives the final state legs""" 1926 1927 return filter(lambda leg: leg.get('state') == True, 1928 self.get('legs'))
1929
1930 - def get_final_ids(self):
1931 """Gives the pdg codes for final state particles""" 1932 1933 return [l.get('id') for l in self.get_final_legs()]
1934
1935 - def get_legs_with_decays(self):
1936 """Return process with all decay chains substituted in.""" 1937 1938 legs = copy.deepcopy(self.get('legs')) 1939 if self.get('is_decay_chain'): 1940 legs.pop(0) 1941 org_decay_chains = copy.copy(self.get('decay_chains')) 1942 sorted_decay_chains = [] 1943 # Sort decay chains according to leg order 1944 for leg in legs: 1945 if not leg.get('state'): continue 1946 org_ids = [l.get('legs')[0].get('id') for l in \ 1947 org_decay_chains] 1948 if leg.get('id') in org_ids: 1949 sorted_decay_chains.append(org_decay_chains.pop(\ 1950 org_ids.index(leg.get('id')))) 1951 assert not org_decay_chains 1952 ileg = 0 1953 for decay in sorted_decay_chains: 1954 while legs[ileg].get('state') == False or \ 1955 legs[ileg].get('id') != decay.get('legs')[0].get('id'): 1956 ileg = ileg + 1 1957 decay_legs = decay.get_legs_with_decays() 1958 legs = legs[:ileg] + decay_legs + legs[ileg+1:] 1959 ileg = ileg + len(decay_legs) 1960 1961 for ileg, leg in enumerate(legs): 1962 leg.set('number', ileg + 1) 1963 1964 return LegList(legs)
1965
1966 - def list_for_sort(self):
1967 """Output a list that can be compared to other processes as: 1968 [id, sorted(initial leg ids), sorted(final leg ids), 1969 sorted(decay list_for_sorts)]""" 1970 1971 sorted_list = [self.get('id'), 1972 sorted(self.get_initial_ids()), 1973 sorted(self.get_final_ids())] 1974 1975 if self.get('decay_chains'): 1976 sorted_list.extend(sorted([d.list_for_sort() for d in \ 1977 self.get('decay_chains')])) 1978 1979 return sorted_list
1980
1981 - def compare_for_sort(self, other):
1982 """Sorting routine which allows to sort processes for 1983 comparison. Compare only process id and legs.""" 1984 1985 if self.list_for_sort() > other.list_for_sort(): 1986 return 1 1987 if self.list_for_sort() < other.list_for_sort(): 1988 return -1 1989 return 0
1990
1991 - def identical_particle_factor(self):
1992 """Calculate the denominator factor for identical final state particles 1993 """ 1994 1995 final_legs = filter(lambda leg: leg.get('state') == True, \ 1996 self.get_legs_with_decays()) 1997 1998 identical_indices = {} 1999 for leg in final_legs: 2000 if leg.get('id') in identical_indices: 2001 identical_indices[leg.get('id')] = \ 2002 identical_indices[leg.get('id')] + 1 2003 else: 2004 identical_indices[leg.get('id')] = 1 2005 return reduce(lambda x, y: x * y, [ math.factorial(val) for val in \ 2006 identical_indices.values() ], 1)
2007
2008 - def check_expansion_orders(self):
2009 """Ensure that maximum expansion orders from the model are 2010 properly taken into account in the process""" 2011 2012 # Ensure that expansion orders are taken into account 2013 expansion_orders = self.get('model').get('expansion_order') 2014 orders = self.get('orders') 2015 2016 tmp = [(k,v) for (k,v) in expansion_orders.items() if 0 < v < 99] 2017 for (k,v) in tmp: 2018 if k in orders: 2019 if v < orders[k]: 2020 logger.warning('''The coupling order (%s=%s) specified is larger than the one allowed 2021 by the model builder. The maximal value allowed is %s. 2022 We set the %s order to this value''' % (k,orders[k],v,k)) 2023 orders[k] = v 2024 else: 2025 orders[k] = v
2026
2027 - def __eq__(self, other):
2028 """Overloading the equality operator, so that only comparison 2029 of process id and legs is being done, using compare_for_sort.""" 2030 2031 if not isinstance(other, Process): 2032 return False 2033 2034 return self.compare_for_sort(other) == 0
2035
2036 - def __ne__(self, other):
2037 return not self.__eq__(other)
2038
2039 #=============================================================================== 2040 # ProcessList 2041 #=============================================================================== 2042 -class ProcessList(PhysicsObjectList):
2043 """List of Process objects 2044 """ 2045
2046 - def is_valid_element(self, obj):
2047 """Test if object obj is a valid Process for the list.""" 2048 2049 return isinstance(obj, Process)
2050
2051 - def nice_string(self, indent = 0):
2052 """Returns a nicely formatted string of the matrix element processes.""" 2053 2054 mystr = "\n".join([p.nice_string(indent) for p in self]) 2055 2056 return mystr
2057
2058 #=============================================================================== 2059 # ProcessDefinition 2060 #=============================================================================== 2061 -class ProcessDefinition(Process):
2062 """ProcessDefinition: list of multilegs (ordered) 2063 dictionary of orders 2064 model 2065 process id 2066 """ 2067
2068 - def default_setup(self):
2069 """Default values for all properties""" 2070 2071 super(ProcessDefinition, self).default_setup() 2072 2073 self['legs'] = MultiLegList() 2074 # Decay chain processes associated with this process 2075 self['decay_chains'] = ProcessDefinitionList()
2076
2077 - def filter(self, name, value):
2078 """Filter for valid process property values.""" 2079 2080 if name == 'legs': 2081 if not isinstance(value, MultiLegList): 2082 raise self.PhysicsObjectError, \ 2083 "%s is not a valid MultiLegList object" % str(value) 2084 elif name == 'decay_chains': 2085 if not isinstance(value, ProcessDefinitionList): 2086 raise self.PhysicsObjectError, \ 2087 "%s is not a valid ProcessDefinitionList" % str(value) 2088 2089 else: 2090 return super(ProcessDefinition, self).filter(name, value) 2091 2092 return True
2093
2094 - def get_sorted_keys(self):
2095 """Return process property names as a nicely sorted list.""" 2096 2097 return super(ProcessDefinition, self).get_sorted_keys()
2098
2099 - def get_minimum_WEIGHTED(self):
2100 """Retrieve the minimum starting guess for WEIGHTED order, to 2101 use in find_optimal_process_orders in MultiProcess diagram 2102 generation (as well as particles and hierarchy). The algorithm: 2103 2104 1) Pick out the legs in the multiprocess according to the 2105 highest hierarchy represented (so don't mix particles from 2106 different hierarchy classes in the same multiparticles!) 2107 2108 2) Find the starting maximum WEIGHTED order as the sum of the 2109 highest n-2 weighted orders 2110 2111 3) Pick out required s-channel particle hierarchies, and use 2112 the highest of the maximum WEIGHTED order from the legs and 2113 the minimum WEIGHTED order extracted from 2*s-channel 2114 hierarchys plus the n-2-2*(number of s-channels) lowest 2115 leg weighted orders. 2116 """ 2117 2118 model = self.get('model') 2119 2120 # Extract hierarchy and particles corresponding to the 2121 # different hierarchy levels from the model 2122 particles, hierarchy = model.get_particles_hierarchy() 2123 2124 # Find legs corresponding to the different orders 2125 # making sure we look at lowest hierarchy first for each leg 2126 max_order_now = [] 2127 new_legs = copy.copy(self.get('legs')) 2128 for parts, value in zip(particles, hierarchy): 2129 ileg = 0 2130 while ileg < len(new_legs): 2131 if any([id in parts for id in new_legs[ileg].get('ids')]): 2132 max_order_now.append(value) 2133 new_legs.pop(ileg) 2134 else: 2135 ileg += 1 2136 2137 # Now remove the two lowest orders to get maximum (since the 2138 # number of interactions is n-2) 2139 max_order_now = sorted(max_order_now)[2:] 2140 2141 # Find s-channel propagators corresponding to the different orders 2142 max_order_prop = [] 2143 for idlist in self.get('required_s_channels'): 2144 max_order_prop.append([0,0]) 2145 for id in idlist: 2146 for parts, value in zip(particles, hierarchy): 2147 if id in parts: 2148 max_order_prop[-1][0] += 2*value 2149 max_order_prop[-1][1] += 1 2150 break 2151 2152 if max_order_prop: 2153 if len(max_order_prop) >1: 2154 max_order_prop = min(*max_order_prop, key=lambda x:x[0]) 2155 else: 2156 max_order_prop = max_order_prop[0] 2157 2158 # Use either the max_order from the external legs or 2159 # the maximum order from the s-channel propagators, plus 2160 # the appropriate lowest orders from max_order_now 2161 max_order_now = max(sum(max_order_now), 2162 max_order_prop[0] + \ 2163 sum(max_order_now[:-2 * max_order_prop[1]])) 2164 else: 2165 max_order_now = sum(max_order_now) 2166 2167 return max_order_now, particles, hierarchy
2168
2169 - def nice_string(self, indent=0):
2170 """Returns a nicely formated string about current process 2171 content""" 2172 2173 mystr = " " * indent + "Process: " 2174 prevleg = None 2175 for leg in self['legs']: 2176 myparts = \ 2177 "/".join([self['model'].get('particle_dict')[id].get_name() \ 2178 for id in leg.get('ids')]) 2179 if prevleg and prevleg['state'] == False \ 2180 and leg['state'] == True: 2181 # Separate initial and final legs by ">" 2182 mystr = mystr + '> ' 2183 # Add required s-channels 2184 if self['required_s_channels'] and \ 2185 self['required_s_channels'][0]: 2186 mystr += "|".join([" ".join([self['model'].\ 2187 get('particle_dict')[req_id].get_name() \ 2188 for req_id in id_list]) \ 2189 for id_list in self['required_s_channels']]) 2190 mystr = mystr + '> ' 2191 2192 mystr = mystr + myparts + ' ' 2193 #mystr = mystr + '(%i) ' % leg['number'] 2194 prevleg = leg 2195 2196 # Add forbidden s-channels 2197 if self['forbidden_onsh_s_channels']: 2198 mystr = mystr + '$ ' 2199 for forb_id in self['forbidden_onsh_s_channels']: 2200 forbpart = self['model'].get('particle_dict')[forb_id] 2201 mystr = mystr + forbpart.get_name() + ' ' 2202 2203 # Add double forbidden s-channels 2204 if self['forbidden_s_channels']: 2205 mystr = mystr + '$$ ' 2206 for forb_id in self['forbidden_s_channels']: 2207 forbpart = self['model'].get('particle_dict')[forb_id] 2208 mystr = mystr + forbpart.get_name() + ' ' 2209 2210 # Add forbidden particles 2211 if self['forbidden_particles']: 2212 mystr = mystr + '/ ' 2213 for forb_id in self['forbidden_particles']: 2214 forbpart = self['model'].get('particle_dict')[forb_id] 2215 mystr = mystr + forbpart.get_name() + ' ' 2216 2217 if self['orders']: 2218 mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \ 2219 for key in sorted(self['orders'])]) + ' ' 2220 2221 # Remove last space 2222 mystr = mystr[:-1] 2223 2224 if self.get('id') or self.get('overall_orders'): 2225 mystr += " @%d" % self.get('id') 2226 if self.get('overall_orders'): 2227 mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \ 2228 for key in sorted(self['orders'])]) + ' ' 2229 2230 if not self.get('decay_chains'): 2231 return mystr 2232 2233 for decay in self['decay_chains']: 2234 mystr = mystr + '\n' + \ 2235 decay.nice_string(indent + 2).replace('Process', 'Decay') 2236 2237 return mystr
2238
2239 - def __eq__(self, other):
2240 """Overloading the equality operator, so that only comparison 2241 of process id and legs is being done, using compare_for_sort.""" 2242 2243 return super(Process, self).__eq__(other)
2244
2245 #=============================================================================== 2246 # ProcessDefinitionList 2247 #=============================================================================== 2248 -class ProcessDefinitionList(PhysicsObjectList):
2249 """List of ProcessDefinition objects 2250 """ 2251
2252 - def is_valid_element(self, obj):
2253 """Test if object obj is a valid ProcessDefinition for the list.""" 2254 2255 return isinstance(obj, ProcessDefinition)
2256
2257 #=============================================================================== 2258 # Global helper functions 2259 #=============================================================================== 2260 2261 -def make_unique(doubletlist):
2262 """Make sure there are no doublets in the list doubletlist. 2263 Note that this is a slow implementation, so don't use if speed 2264 is needed""" 2265 2266 assert isinstance(doubletlist, list), \ 2267 "Argument to make_unique must be list" 2268 2269 2270 uniquelist = [] 2271 for elem in doubletlist: 2272 if elem not in uniquelist: 2273 uniquelist.append(elem) 2274 2275 doubletlist[:] = uniquelist[:]
2276