models.import

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 from compiler.ast import Continue 16 """ How to import a UFO model to the MG5 format """ 17 18 19 import fractions 20 import logging 21 import os 22 import re 23 import sys 24 25 from madgraph import MadGraph5Error, MG5DIR 26 import madgraph.core.base_objects as base_objects 27 import madgraph.core.color_algebra as color 28 import madgraph.iolibs.files as files 29 import madgraph.iolibs.save_load_object as save_load_object 30 from madgraph.core.color_algebra import * 31 32 import madgraph.various.misc as misc 33 34 import aloha.create_aloha as create_aloha 35 import aloha.aloha_fct as aloha_fct 36 37 import models as ufomodels 38 import models.model_reader as model_reader 39 logger = logging.getLogger('madgraph.model') 40 logger_mod = logging.getLogger('madgraph.model') 41 42 root_path = os.path.dirname(os.path.realpath( __file__ )) 43 sys.path.append(root_path) 44 45 sys.path.append(os.path.join(root_path, os.path.pardir, 'Template', 'bin', 'internal')) 46 import check_param_card 47 48 49

50 -class UFOImportError(MadGraph5Error):

51 """ a error class for wrong import of UFO model"""

52

53 -class InvalidModel(MadGraph5Error):

54 """ a class for invalid Model """

55

56 -def find_ufo_path(model_name):

57 """ find the path to a model """ 58 59 # Check for a valid directory 60 if model_name.startswith('./') and os.path.isdir(model_name): 61 model_path = model_name 62 elif os.path.isdir(os.path.join(MG5DIR, 'models', model_name)): 63 model_path = os.path.join(MG5DIR, 'models', model_name) 64 elif os.path.isdir(model_name): 65 model_path = model_name 66 else: 67 raise UFOImportError("Path %s is not a valid pathname" % model_name) 68 69 return model_path

70

71 -def import_model(model_name):

72 """ a practical and efficient way to import a model""" 73 74 # check if this is a valid path or if this include restriction file 75 try: 76 model_path = find_ufo_path(model_name) 77 except UFOImportError: 78 if '-' not in model_name: 79 raise 80 split = model_name.split('-') 81 model_name = '-'.join([text for text in split[:-1]]) 82 model_path = find_ufo_path(model_name) 83 restrict_name = split[-1] 84 85 restrict_file = os.path.join(model_path, 'restrict_%s.dat'% restrict_name) 86 87 #if restriction is full, then we by pass restriction (avoid default) 88 if split[-1] == 'full': 89 restrict_file = None 90 else: 91 # Check if by default we need some restrictions 92 restrict_name = "" 93 if os.path.exists(os.path.join(model_path,'restrict_default.dat')): 94 restrict_file = os.path.join(model_path,'restrict_default.dat') 95 else: 96 restrict_file = None 97 98 #import the FULL model 99 model = import_full_model(model_path) 100 # restore the model name 101 if restrict_name: 102 model["name"] += '-' + restrict_name 103 104 #restrict it if needed 105 if restrict_file: 106 try: 107 logger.info('Restrict model %s with file %s .' % (model_name, os.path.relpath(restrict_file))) 108 except OSError: 109 # sometimes has trouble with relative path 110 logger.info('Restrict model %s with file %s .' % (model_name, restrict_file)) 111 112 if logger_mod.getEffectiveLevel() > 10: 113 logger.info('Run \"set stdout_level DEBUG\" before import for more information.') 114 115 # Modify the mother class of the object in order to allow restriction 116 model = RestrictModel(model) 117 model.restrict_model(restrict_file) 118 119 return model

120 121 _import_once = []

122 -def import_full_model(model_path):

123 """ a practical and efficient way to import one of those models 124 (no restriction file use)""" 125 126 assert model_path == find_ufo_path(model_path) 127 128 # Check the validity of the model 129 files_list_prov = ['couplings.py','lorentz.py','parameters.py', 130 'particles.py', 'vertices.py'] 131 files_list = [] 132 for filename in files_list_prov: 133 filepath = os.path.join(model_path, filename) 134 if not os.path.isfile(filepath): 135 raise UFOImportError, "%s directory is not a valid UFO model: \n %s is missing" % \ 136 (model_path, filename) 137 files_list.append(filepath) 138 139 # use pickle files if defined and up-to-date 140 if files.is_uptodate(os.path.join(model_path, 'model.pkl'), files_list): 141 try: 142 model = save_load_object.load_from_file( \ 143 os.path.join(model_path, 'model.pkl')) 144 except Exception, error: 145 logger.info('failed to load model from pickle file. Try importing UFO from File') 146 else: 147 # check path is correct 148 if model.has_key('version_tag') and model.get('version_tag') == os.path.realpath(model_path) + str(misc.get_pkg_info()): 149 _import_once.append(model_path) 150 return model 151 152 if model_path in _import_once: 153 raise MadGraph5Error, 'This model is modified on disk. To reload it you need to quit/relaunch mg5' 154 155 # Load basic information 156 ufo_model = ufomodels.load_model(model_path) 157 ufo2mg5_converter = UFOMG5Converter(ufo_model) 158 model = ufo2mg5_converter.load_model() 159 160 if model_path[-1] == '/': model_path = model_path[:-1] #avoid empty name 161 model.set('name', os.path.split(model_path)[-1]) 162 model.set('version_tag', os.path.realpath(model_path) + str(misc.get_pkg_info())) 163 164 # Load the Parameter/Coupling in a convinient format. 165 parameters, couplings = OrganizeModelExpression(ufo_model).main() 166 model.set('parameters', parameters) 167 model.set('couplings', couplings) 168 model.set('functions', ufo_model.all_functions) 169 170 # save in a pickle files to fasten future usage 171 save_load_object.save_to_file(os.path.join(model_path, 'model.pkl'), model) 172 173 #if default and os.path.exists(os.path.join(model_path, 'restrict_default.dat')): 174 # restrict_file = os.path.join(model_path, 'restrict_default.dat') 175 # model = import_ufo.RestrictModel(model) 176 # model.restrict_model(restrict_file) 177 return model

178 179

180 -class UFOMG5Converter(object):

181 """Convert a UFO model to the MG5 format""" 182 183 use_lower_part_names = False 184

185 - def __init__(self, model, auto=False):

186 """ initialize empty list for particles/interactions """ 187 188 self.particles = base_objects.ParticleList() 189 self.interactions = base_objects.InteractionList() 190 self.model = base_objects.Model() 191 self.model.set('particles', self.particles) 192 self.model.set('interactions', self.interactions) 193 self.conservecharge = set(['charge']) 194 195 self.ufomodel = model 196 self.checked_lor = set() 197 198 if auto: 199 self.load_model()

200

201 - def load_model(self):

202 """load the different of the model first particles then interactions""" 203 204 # Check the validity of the model 205 # 1) check that all lhablock are single word. 206 for param in self.ufomodel.all_parameters: 207 if param.nature == "external": 208 if len(param.lhablock.split())>1: 209 raise InvalidModel, '''LHABlock should be single word which is not the case for 210 \'%s\' parameter with lhablock \'%s\'''' % (param.name, param.lhablock) 211 212 213 logger.info('load particles') 214 # Check if multiple particles have the same name but different case. 215 # Otherwise, we can use lowercase particle names. 216 if len(set([p.name for p in self.ufomodel.all_particles] + \ 217 [p.antiname for p in self.ufomodel.all_particles])) == \ 218 len(set([p.name.lower() for p in self.ufomodel.all_particles] + \ 219 [p.antiname.lower() for p in self.ufomodel.all_particles])): 220 self.use_lower_part_names = True 221 222 # check which of the fermion/anti-fermion should be set as incoming 223 self.detect_incoming_fermion() 224 225 for particle_info in self.ufomodel.all_particles: 226 self.add_particle(particle_info) 227 228 # Find which particles is in the 3/3bar color states (retrun {id: 3/-3}) 229 color_info = self.find_color_anti_color_rep() 230 231 232 logger.info('load vertices') 233 for interaction_info in self.ufomodel.all_vertices: 234 self.add_interaction(interaction_info, color_info) 235 236 self.model.set('conserved_charge', self.conservecharge) 237 238 # If we deal with a Loop model here, the order hierarchy MUST be 239 # defined in the file coupling_orders.py and we import it from 240 # there. 241 242 hierarchy={} 243 try: 244 all_orders = self.ufomodel.all_orders 245 for order in all_orders: 246 hierarchy[order.name]=order.hierarchy 247 except AttributeError: 248 pass 249 else: 250 self.model.set('order_hierarchy', hierarchy) 251 252 # Also set expansion_order, i.e., maximum coupling order per process 253 254 expansion_order={} 255 try: 256 all_orders = self.ufomodel.all_orders 257 for order in all_orders: 258 expansion_order[order.name]=order.expansion_order 259 except AttributeError: 260 pass 261 else: 262 self.model.set('expansion_order', expansion_order) 263 264 #clean memory 265 del self.checked_lor 266 267 return self.model

268 269

270 - def add_particle(self, particle_info):

271 """ convert and add a particle in the particle list """ 272 273 # MG5 have only one entry for particle and anti particles. 274 #UFO has two. use the color to avoid duplictions 275 pdg = particle_info.pdg_code 276 if pdg in self.incoming or (pdg not in self.outcoming and pdg <0): 277 return 278 279 # MG5 doesn't use ghost (use unitary gauges) 280 if particle_info.spin < 0: 281 return 282 283 # MG5 doesn't use goldstone boson 284 if hasattr(particle_info, 'GoldstoneBoson'): 285 if particle_info.GoldstoneBoson: 286 return 287 288 # Initialize a particles 289 particle = base_objects.Particle() 290 291 nb_property = 0 #basic check that the UFO information is complete 292 # Loop over the element defining the UFO particles 293 for key,value in particle_info.__dict__.items(): 294 # Check if we use it in the MG5 definition of a particles 295 if key in base_objects.Particle.sorted_keys: 296 nb_property +=1 297 if key in ['name', 'antiname']: 298 if self.use_lower_part_names: 299 particle.set(key, value.lower()) 300 else: 301 particle.set(key, value) 302 elif key == 'charge': 303 particle.set(key, float(value)) 304 elif key in ['mass','width']: 305 particle.set(key, str(value)) 306 else: 307 particle.set(key, value) 308 elif key.lower() not in ('ghostnumber','selfconjugate','goldstoneboson'): 309 # add charge -we will check later if those are conserve 310 self.conservecharge.add(key) 311 particle.set(key,value, force=True) 312 313 assert(12 == nb_property) #basic check that all the information is there 314 315 # Identify self conjugate particles 316 if particle_info.name == particle_info.antiname: 317 particle.set('self_antipart', True) 318 319 # Add the particles to the list 320 self.particles.append(particle)

321 322

323 - def find_color_anti_color_rep(self):

324 """find which color are in the 3/3bar states""" 325 # method look at the 3 3bar 8 configuration. 326 # If the color is T(3,2,1) and the interaction F1 F2 V 327 # Then set F1 to anticolor (and F2 to color) 328 # if this is T(3,1,2) set the opposite 329 output = {} 330 331 for interaction_info in self.ufomodel.all_vertices: 332 if len(interaction_info.particles) != 3: 333 continue 334 colors = [abs(p.color) for p in interaction_info.particles] 335 if colors[:2] == [3,3]: 336 if 'T(3,2,1)' in interaction_info.color: 337 color, anticolor, other = interaction_info.particles 338 elif 'T(3,1,2)' in interaction_info.color: 339 anticolor, color, other = interaction_info.particles 340 else: 341 continue 342 elif colors[1:] == [3,3]: 343 if 'T(1,2,3)' in interaction_info.color: 344 other, anticolor, color = interaction_info.particles 345 elif 'T(1,3,2)' in interaction_info.color: 346 other, color, anticolor = interaction_info.particles 347 else: 348 continue 349 350 elif colors.count(3) == 2: 351 if 'T(2,3,1)' in interaction_info.color: 352 color, other, anticolor = interaction_info.particles 353 elif 'T(2,1,3)' in interaction_info.color: 354 anticolor, other, color = interaction_info.particles 355 else: 356 continue 357 else: 358 continue 359 360 # Check/assign for the color particle 361 if color.pdg_code in output: 362 if output[color.pdg_code] == -3: 363 raise InvalidModel, 'Particles %s is sometimes in the 3 and sometimes in the 3bar representations' \ 364 % color.name 365 else: 366 output[color.pdg_code] = 3 367 368 # Check/assign for the anticolor particle 369 if anticolor.pdg_code in output: 370 if output[anticolor.pdg_code] == 3: 371 raise InvalidModel, 'Particles %s is sometimes set as in the 3 and sometimes in the 3bar representations' \ 372 % anticolor.name 373 else: 374 output[anticolor.pdg_code] = -3 375 376 return output

377

378 - def detect_incoming_fermion(self):

379 """define which fermion should be incoming 380 for that we look at F F~ X interactions 381 """ 382 self.incoming = [] 383 self.outcoming = [] 384 for interaction_info in self.ufomodel.all_vertices: 385 # check if the interaction meet requirements: 386 pdg = [p.pdg_code for p in interaction_info.particles if p.spin in [2,4]] 387 if len(pdg) % 2: 388 raise InvalidModel, 'Odd number of fermion in vertex: %s' % [p.pdg_code for p in interaction_info.particles] 389 for i in range(0, len(pdg),2): 390 if pdg[i] == - pdg[i+1]: 391 if pdg[i] in self.outcoming: 392 raise InvalidModel, '%s has not coherent incoming/outcoming status between interactions' %\ 393 [p for p in interaction_info.particles if p.spin in [2,4]][i].name 394 395 elif not pdg[i] in self.incoming: 396 self.incoming.append(pdg[i]) 397 self.outcoming.append(pdg[i+1])

398 399

400 - def add_interaction(self, interaction_info, color_info):

401 """add an interaction in the MG5 model. interaction_info is the 402 UFO vertices information.""" 403 404 # Import particles content: 405 particles = [self.model.get_particle(particle.pdg_code) \ 406 for particle in interaction_info.particles] 407 408 if None in particles: 409 # Interaction with a ghost/goldstone 410 return 411 412 particles = base_objects.ParticleList(particles) 413 414 # Check the coherence of the Fermion Flow 415 nb_fermion = sum([p.is_fermion() and 1 or 0 for p in particles]) 416 try: 417 if nb_fermion: 418 [aloha_fct.check_flow_validity(helas.structure, nb_fermion) \ 419 for helas in interaction_info.lorentz 420 if helas.name not in self.checked_lor] 421 self.checked_lor.update(set([helas.name for helas in interaction_info.lorentz])) 422 except aloha_fct.WrongFermionFlow, error: 423 text = 'Fermion Flow error for interactions %s: %s: %s\n %s' % \ 424 (', '.join([p.name for p in interaction_info.particles]), 425 helas.name, helas.structure, error) 426 raise InvalidModel, text 427 428 # Import Lorentz content: 429 lorentz = [helas.name for helas in interaction_info.lorentz] 430 431 # Import color information: 432 colors = [self.treat_color(color_obj, interaction_info, color_info) for color_obj in \ 433 interaction_info.color] 434 435 order_to_int={} 436 437 for key, couplings in interaction_info.couplings.items(): 438 if not isinstance(couplings, list): 439 couplings = [couplings] 440 for coupling in couplings: 441 order = tuple(coupling.order.items()) 442 if '1' in order: 443 raise InvalidModel, '''Some couplings have \'1\' order. 444 This is not allowed in MG. 445 Please defines an additional coupling to your model''' 446 if order in order_to_int: 447 order_to_int[order].get('couplings')[key] = coupling.name 448 else: 449 # Initialize a new interaction with a new id tag 450 interaction = base_objects.Interaction({'id':len(self.interactions)+1}) 451 interaction.set('particles', particles) 452 interaction.set('lorentz', lorentz) 453 interaction.set('couplings', {key: coupling.name}) 454 interaction.set('orders', coupling.order) 455 interaction.set('color', colors) 456 order_to_int[order] = interaction 457 # add to the interactions 458 self.interactions.append(interaction) 459 460 # check if this interaction conserve the charge defined 461 for charge in list(self.conservecharge): #duplicate to allow modification 462 total = 0 463 for part in interaction_info.particles: 464 try: 465 total += getattr(part, charge) 466 except AttributeError: 467 pass 468 if abs(total) > 1e-12: 469 logger.info('The model has interaction violating the charge: %s' % charge) 470 self.conservecharge.discard(charge)

471 472 _pat_T = re.compile(r'T\((?P<first>\d*),(?P<second>\d*)\)') 473 _pat_id = re.compile(r'Identity\((?P<first>\d*),(?P<second>\d*)\)') 474

475 - def treat_color(self, data_string, interaction_info, color_info):

476 """ convert the string to ColorString""" 477 478 #original = copy.copy(data_string) 479 #data_string = p.sub('color.T(\g<first>,\g<second>)', data_string) 480 481 482 output = [] 483 factor = 1 484 for term in data_string.split('*'): 485 pattern = self._pat_id.search(term) 486 if pattern: 487 particle = interaction_info.particles[int(pattern.group('first'))-1] 488 particle2 = interaction_info.particles[int(pattern.group('second'))-1] 489 if particle.color == particle2.color and particle.color in [-6, 6]: 490 error_msg = 'UFO model have inconsistency in the format:\n' 491 error_msg += 'interactions for particles %s has color information %s\n' 492 error_msg += ' but both fermion are in the same representation %s' 493 raise UFOFormatError, error_msg % (', '.join([p.name for p in interaction_info.particles]),data_string, particle.color) 494 if particle.color == particle2.color and particle.color in [-3, 3]: 495 if particle.pdg_code in color_info and particle2.pdg_code in color_info: 496 if color_info[particle.pdg_code] == color_info[particle2.pdg_code]: 497 error_msg = 'UFO model have inconsistency in the format:\n' 498 error_msg += 'interactions for particles %s has color information %s\n' 499 error_msg += ' but both fermion are in the same representation %s' 500 raise UFOFormatError, error_msg % (', '.join([p.name for p in interaction_info.particles]),data_string, particle.color) 501 elif particle.pdg_code in color_info: 502 color_info[particle2.pdg_code] = -particle.pdg_code 503 elif particle2.pdg_code in color_info: 504 color_info[particle.pdg_code] = -particle2.pdg_code 505 else: 506 error_msg = 'UFO model have inconsistency in the format:\n' 507 error_msg += 'interactions for particles %s has color information %s\n' 508 error_msg += ' but both fermion are in the same representation %s' 509 raise UFOFormatError, error_msg % (', '.join([p.name for p in interaction_info.particles]),data_string, particle.color) 510 511 512 if particle.color == 6: 513 output.append(self._pat_id.sub('color.T6(\g<first>,\g<second>)', term)) 514 elif particle.color == -6 : 515 output.append(self._pat_id.sub('color.T6(\g<second>,\g<first>)', term)) 516 elif particle.color == 8: 517 output.append(self._pat_id.sub('color.Tr(\g<first>,\g<second>)', term)) 518 factor *= 2 519 elif particle.color in [-3,3]: 520 if particle.pdg_code not in color_info: 521 logger.debug('Not able to find the 3/3bar rep from the interactions for particle %s' % particle.name) 522 color_info[particle.pdg_code] = particle.color 523 if particle2.pdg_code not in color_info: 524 logger.debug('Not able to find the 3/3bar rep from the interactions for particle %s' % particle2.name) 525 color_info[particle2.pdg_code] = particle2.color 526 527 528 if color_info[particle.pdg_code] == 3 : 529 output.append(self._pat_id.sub('color.T(\g<second>,\g<first>)', term)) 530 elif color_info[particle.pdg_code] == -3: 531 output.append(self._pat_id.sub('color.T(\g<first>,\g<second>)', term)) 532 else: 533 raise MadGraph5Error, \ 534 "Unknown use of Identity for particle with color %d" \ 535 % particle.color 536 else: 537 output.append(term) 538 data_string = '*'.join(output) 539 540 # Change convention for summed indices 541 p = re.compile(r'''\'\w(?P<number>\d+)\'''') 542 data_string = p.sub('-\g<number>', data_string) 543 544 # Shift indices by -1 545 new_indices = {} 546 new_indices = dict([(j,i) for (i,j) in \ 547 enumerate(range(1, 548 len(interaction_info.particles)+1))]) 549 550 551 output = data_string.split('*') 552 output = color.ColorString([eval(data) \ 553 for data in output if data !='1']) 554 output.coeff = fractions.Fraction(factor) 555 for col_obj in output: 556 col_obj.replace_indices(new_indices) 557 558 return output

559

560 -class OrganizeModelExpression:

561 """Organize the couplings/parameters of a model""" 562 563 track_dependant = ['aS','aEWM1'] # list of variable from which we track 564 #dependencies those variables should be define 565 #as external parameters 566 567 # regular expression to shorten the expressions 568 complex_number = re.compile(r'''complex\((?P<real>[^,\(\)]+),(?P<imag>[^,\(\)]+)\)''') 569 expo_expr = re.compile(r'''(?P<expr>[\w.]+)\s*\*\*\s*(?P<expo>[\d.+-]+)''') 570 cmath_expr = re.compile(r'''cmath.(?P<operation>\w+)\((?P<expr>\w+)\)''') 571 #operation is usualy sqrt / sin / cos / tan 572 conj_expr = re.compile(r'''complexconjugate\((?P<expr>\w+)\)''') 573 574 #RE expression for is_event_dependent 575 separator = re.compile(r'''[+,\-*/()]''') 576

577 - def __init__(self, model):

578 579 self.model = model # UFOMODEL 580 self.params = {} # depend on -> ModelVariable 581 self.couplings = {} # depend on -> ModelVariable 582 self.all_expr = {} # variable_name -> ModelVariable

583

584 - def main(self):

585 """Launch the actual computation and return the associate 586 params/couplings.""" 587 588 self.analyze_parameters() 589 self.analyze_couplings() 590 return self.params, self.couplings

591 592

593 - def analyze_parameters(self):

594 """ separate the parameters needed to be recomputed events by events and 595 the others""" 596 597 for param in self.model.all_parameters: 598 if param.nature == 'external': 599 parameter = base_objects.ParamCardVariable(param.name, param.value, \ 600 param.lhablock, param.lhacode) 601 602 else: 603 expr = self.shorten_expr(param.value) 604 depend_on = self.find_dependencies(expr) 605 parameter = base_objects.ModelVariable(param.name, expr, param.type, depend_on) 606 607 self.add_parameter(parameter)

608 609

610 - def add_parameter(self, parameter):

611 """ add consistently the parameter in params and all_expr. 612 avoid duplication """ 613 614 assert isinstance(parameter, base_objects.ModelVariable) 615 616 if parameter.name in self.all_expr.keys(): 617 return 618 619 self.all_expr[parameter.name] = parameter 620 try: 621 self.params[parameter.depend].append(parameter) 622 except: 623 self.params[parameter.depend] = [parameter]

624

625 - def add_coupling(self, coupling):

626 """ add consistently the coupling in couplings and all_expr. 627 avoid duplication """ 628 629 assert isinstance(coupling, base_objects.ModelVariable) 630 631 if coupling.name in self.all_expr.keys(): 632 return 633 self.all_expr[coupling.value] = coupling 634 try: 635 self.coupling[coupling.depend].append(coupling) 636 except: 637 self.coupling[coupling.depend] = [coupling]

638 639 640

641 - def analyze_couplings(self):

642 """creates the shortcut for all special function/parameter 643 separate the couplings dependent of track variables of the others""" 644 645 for coupling in self.model.all_couplings: 646 647 # shorten expression, find dependencies, create short object 648 expr = self.shorten_expr(coupling.value) 649 depend_on = self.find_dependencies(expr) 650 parameter = base_objects.ModelVariable(coupling.name, expr, 'complex', depend_on) 651 652 # Add consistently in the couplings/all_expr 653 try: 654 self.couplings[depend_on].append(parameter) 655 except KeyError: 656 self.couplings[depend_on] = [parameter] 657 self.all_expr[coupling.value] = parameter

658 659

660 - def find_dependencies(self, expr):

661 """check if an expression should be evaluated points by points or not 662 """ 663 depend_on = set() 664 665 # Treat predefined result 666 #if name in self.track_dependant: 667 # return tuple() 668 669 # Split the different part of the expression in order to say if a 670 #subexpression is dependent of one of tracked variable 671 expr = self.separator.sub(' ',expr) 672 673 # look for each subexpression 674 for subexpr in expr.split(): 675 if subexpr in self.track_dependant: 676 depend_on.add(subexpr) 677 678 elif subexpr in self.all_expr.keys() and self.all_expr[subexpr].depend: 679 [depend_on.add(value) for value in self.all_expr[subexpr].depend 680 if self.all_expr[subexpr].depend != ('external',)] 681 682 if depend_on: 683 return tuple(depend_on) 684 else: 685 return tuple()

686 687

688 - def shorten_expr(self, expr):

689 """ apply the rules of contraction and fullfill 690 self.params with dependent part""" 691 692 expr = self.complex_number.sub(self.shorten_complex, expr) 693 expr = self.expo_expr.sub(self.shorten_expo, expr) 694 expr = self.cmath_expr.sub(self.shorten_cmath, expr) 695 expr = self.conj_expr.sub(self.shorten_conjugate, expr) 696 return expr

697 698

699 - def shorten_complex(self, matchobj):

700 """add the short expression, and return the nice string associate""" 701 702 real = float(matchobj.group('real')) 703 imag = float(matchobj.group('imag')) 704 if real == 0 and imag ==1: 705 new_param = base_objects.ModelVariable('complexi', 'complex(0,1)', 'complex') 706 self.add_parameter(new_param) 707 return 'complexi' 708 else: 709 return 'complex(%s, %s)' % (real, imag)

710 711

712 - def shorten_expo(self, matchobj):

713 """add the short expression, and return the nice string associate""" 714 715 expr = matchobj.group('expr') 716 exponent = matchobj.group('expo') 717 new_exponent = exponent.replace('.','_').replace('+','').replace('-','_m_') 718 output = '%s__exp__%s' % (expr, new_exponent) 719 old_expr = '%s**%s' % (expr,exponent) 720 721 if expr.startswith('cmath'): 722 return old_expr 723 724 if expr.isdigit(): 725 output = 'nb__' + output #prevent to start with a number 726 new_param = base_objects.ModelVariable(output, old_expr,'real') 727 else: 728 depend_on = self.find_dependencies(expr) 729 type = self.search_type(expr) 730 new_param = base_objects.ModelVariable(output, old_expr, type, depend_on) 731 self.add_parameter(new_param) 732 return output

733

734 - def shorten_cmath(self, matchobj):

735 """add the short expression, and return the nice string associate""" 736 737 expr = matchobj.group('expr') 738 operation = matchobj.group('operation') 739 output = '%s__%s' % (operation, expr) 740 old_expr = ' cmath.%s(%s) ' % (operation, expr) 741 if expr.isdigit(): 742 new_param = base_objects.ModelVariable(output, old_expr , 'real') 743 else: 744 depend_on = self.find_dependencies(expr) 745 type = self.search_type(expr) 746 new_param = base_objects.ModelVariable(output, old_expr, type, depend_on) 747 self.add_parameter(new_param) 748 749 return output

750

751 - def shorten_conjugate(self, matchobj):

752 """add the short expression, and retrun the nice string associate""" 753 754 expr = matchobj.group('expr') 755 output = 'conjg__%s' % (expr) 756 old_expr = ' complexconjugate(%s) ' % expr 757 depend_on = self.find_dependencies(expr) 758 type = 'complex' 759 new_param = base_objects.ModelVariable(output, old_expr, type, depend_on) 760 self.add_parameter(new_param) 761 762 return output

763 764 765

766 - def search_type(self, expr, dep=''):

767 """return the type associate to the expression if define""" 768 769 try: 770 return self.all_expr[expr].type 771 except: 772 return 'complex'

773

774 -class RestrictModel(model_reader.ModelReader):

775 """ A class for restricting a model for a given param_card. 776 rules applied: 777 - Vertex with zero couplings are throw away 778 - external parameter with zero/one input are changed into internal parameter. 779 - identical coupling/mass/width are replace in the model by a unique one 780 """ 781

782 - def default_setup(self):

783 """define default value""" 784 self.del_coup = [] 785 super(RestrictModel, self).default_setup() 786 self.rule_card = check_param_card.ParamCardRule()

787

788 - def restrict_model(self, param_card):

789 """apply the model restriction following param_card""" 790 791 # Reset particle dict to ensure synchronized particles and interactions 792 self.set('particles', self.get('particles')) 793 794 # compute the value of all parameters 795 self.set_parameters_and_couplings(param_card) 796 # associte to each couplings the associated vertex: def self.coupling_pos 797 self.locate_coupling() 798 # deal with couplings 799 zero_couplings, iden_couplings = self.detect_identical_couplings() 800 801 # remove the out-dated interactions 802 self.remove_interactions(zero_couplings) 803 804 # replace in interactions identical couplings 805 for iden_coups in iden_couplings: 806 self.merge_iden_couplings(iden_coups) 807 808 # remove zero couplings and other pointless couplings 809 self.del_coup += zero_couplings 810 self.remove_couplings(self.del_coup) 811 812 # deal with parameters 813 parameters = self.detect_special_parameters() 814 self.fix_parameter_values(*parameters) 815 816 # deal with identical parameters 817 iden_parameters = self.detect_identical_parameters() 818 for iden_param in iden_parameters: 819 self.merge_iden_parameters(iden_param) 820 821 # change value of default parameter if they have special value: 822 # 9.999999e-1 -> 1.0 823 # 0.000001e-99 -> 0 Those value are used to avoid restriction 824 for name, value in self['parameter_dict'].items(): 825 if value == 9.999999e-1: 826 self['parameter_dict'][name] = 1 827 elif value == 0.000001e-99: 828 self['parameter_dict'][name] = 0

829

830 - def locate_coupling(self):

831 """ create a dict couplings_name -> vertex """ 832 833 self.coupling_pos = {} 834 for vertex in self['interactions']: 835 for key, coupling in vertex['couplings'].items(): 836 if coupling in self.coupling_pos: 837 if vertex not in self.coupling_pos[coupling]: 838 self.coupling_pos[coupling].append(vertex) 839 else: 840 self.coupling_pos[coupling] = [vertex] 841 842 return self.coupling_pos

843

844 - def detect_identical_couplings(self, strict_zero=False):

845 """return a list with the name of all vanishing couplings""" 846 847 dict_value_coupling = {} 848 iden_key = set() 849 zero_coupling = [] 850 iden_coupling = [] 851 852 for name, value in self['coupling_dict'].items(): 853 if value == 0: 854 zero_coupling.append(name) 855 continue 856 elif not strict_zero and abs(value) < 1e-13: 857 logger.debug('coupling with small value %s: %s treated as zero' % 858 (name, value)) 859 zero_coupling.append(name) 860 elif not strict_zero and abs(value) < 1e-10: 861 return self.detect_identical_couplings(strict_zero=True) 862 863 864 if value in dict_value_coupling: 865 iden_key.add(value) 866 dict_value_coupling[value].append(name) 867 else: 868 dict_value_coupling[value] = [name] 869 870 for key in iden_key: 871 iden_coupling.append(dict_value_coupling[key]) 872 873 return zero_coupling, iden_coupling

874 875

876 - def detect_special_parameters(self):

877 """ return the list of (name of) parameter which are zero """ 878 879 null_parameters = [] 880 one_parameters = [] 881 for name, value in self['parameter_dict'].items(): 882 if value == 0 and name != 'ZERO': 883 null_parameters.append(name) 884 elif value == 1: 885 one_parameters.append(name) 886 887 return null_parameters, one_parameters

888

889 - def detect_identical_parameters(self):

890 """ return the list of tuple of name of parameter with the same 891 input value """ 892 893 # Extract external parameters 894 external_parameters = self['parameters'][('external',)] 895 896 # define usefull variable to detect identical input 897 block_value_to_var={} #(lhablok, value): list_of_var 898 mult_param = set([]) # key of the previous dict with more than one 899 #parameter. 900 901 #detect identical parameter and remove the duplicate parameter 902 for param in external_parameters[:]: 903 value = self['parameter_dict'][param.name] 904 if value == 0: 905 continue 906 key = (param.lhablock, value) 907 mkey = (param.lhablock, -value) 908 if key in block_value_to_var: 909 block_value_to_var[key].append((param,1)) 910 mult_param.add(key) 911 elif mkey in block_value_to_var: 912 block_value_to_var[mkey].append((param,-1)) 913 mult_param.add(mkey) 914 else: 915 block_value_to_var[key] = [(param,1)] 916 917 output=[] 918 for key in mult_param: 919 output.append(block_value_to_var[key]) 920 921 return output

922 923

924 - def merge_iden_couplings(self, couplings):

925 """merge the identical couplings in the interactions""" 926 927 928 logger_mod.debug(' Fuse the Following coupling (they have the same value): %s '% \ 929 ', '.join([obj for obj in couplings])) 930 931 main = couplings[0] 932 self.del_coup += couplings[1:] # add the other coupl to the suppress list 933 934 for coupling in couplings[1:]: 935 # check if param is linked to an interaction 936 if coupling not in self.coupling_pos: 937 continue 938 # replace the coupling, by checking all coupling of the interaction 939 vertices = self.coupling_pos[coupling] 940 for vertex in vertices: 941 for key, value in vertex['couplings'].items(): 942 if value == coupling: 943 vertex['couplings'][key] = main

944 945

946 - def merge_iden_parameters(self, parameters):

947 """ merge the identical parameters given in argument """ 948 949 logger_mod.debug('Parameters set to identical values: %s '% \ 950 ', '.join(['%s*%s' % (f, obj.name) for (obj,f) in parameters])) 951 952 # Extract external parameters 953 external_parameters = self['parameters'][('external',)] 954 for i, (obj, factor) in enumerate(parameters): 955 # Keeped intact the first one and store information 956 if i == 0: 957 obj.info = 'set of param :' + \ 958 ', '.join([str(f)+'*'+param.name for (param, f) in parameters]) 959 expr = obj.name 960 continue 961 # Add a Rule linked to the param_card 962 if factor ==1: 963 self.rule_card.add_identical(obj.lhablock.lower(), obj.lhacode, 964 parameters[0][0].lhacode ) 965 else: 966 self.rule_card.add_opposite(obj.lhablock.lower(), obj.lhacode, 967 parameters[0][0].lhacode ) 968 # delete the old parameters 969 external_parameters.remove(obj) 970 # replace by the new one pointing of the first obj of the class 971 new_param = base_objects.ModelVariable(obj.name, '%s*%s' %(factor, expr), 'real') 972 self['parameters'][()].insert(0, new_param) 973 974 # For Mass-Width, we need also to replace the mass-width in the particles 975 #This allows some optimization for multi-process. 976 if parameters[0][0].lhablock in ['MASS','DECAY']: 977 new_name = parameters[0][0].name 978 if parameters[0][0].lhablock == 'MASS': 979 arg = 'mass' 980 else: 981 arg = 'width' 982 change_name = [p.name for (p,f) in parameters[1:]] 983 [p.set(arg, new_name) for p in self['particle_dict'].values() 984 if p[arg] in change_name]

985

986 - def remove_interactions(self, zero_couplings):

987 """ remove the interactions associated to couplings""" 988 989 mod = [] 990 for coup in zero_couplings: 991 # some coupling might be not related to any interactions 992 if coup not in self.coupling_pos: 993 coup, self.coupling_pos.keys() 994 continue 995 for vertex in self.coupling_pos[coup]: 996 modify = False 997 for key, coupling in vertex['couplings'].items(): 998 if coupling in zero_couplings: 999 modify=True 1000 del vertex['couplings'][key] 1001 if modify: 1002 mod.append(vertex) 1003 1004 # print usefull log and clean the empty interaction 1005 for vertex in mod: 1006 part_name = [part['name'] for part in vertex['particles']] 1007 orders = ['%s=%s' % (order,value) for order,value in vertex['orders'].items()] 1008 1009 if not vertex['couplings']: 1010 logger_mod.debug('remove interactions: %s at order: %s' % \ 1011 (' '.join(part_name),', '.join(orders))) 1012 self['interactions'].remove(vertex) 1013 else: 1014 logger_mod.debug('modify interactions: %s at order: %s' % \ 1015 (' '.join(part_name),', '.join(orders))) 1016 1017 return

1018

1019 - def remove_couplings(self, couplings):

1020 #clean the coupling list: 1021 for name, data in self['couplings'].items(): 1022 for coupling in data[:]: 1023 if coupling.name in couplings: 1024 data.remove(coupling)

1025 1026

1027 - def fix_parameter_values(self, zero_parameters, one_parameters):

1028 """ Remove all instance of the parameters in the model and replace it by 1029 zero when needed.""" 1030 1031 # Add a rule for zero/one parameter 1032 external_parameters = self['parameters'][('external',)] 1033 for param in external_parameters[:]: 1034 value = self['parameter_dict'][param.name] 1035 block = param.lhablock.lower() 1036 if value == 0: 1037 self.rule_card.add_zero(block, param.lhacode) 1038 elif value == 1: 1039 self.rule_card.add_one(block, param.lhacode) 1040 1041 1042 1043 1044 special_parameters = zero_parameters + one_parameters 1045 1046 # treat specific cases for masses and width 1047 for particle in self['particles']: 1048 if particle['mass'] in zero_parameters: 1049 particle['mass'] = 'ZERO' 1050 if particle['width'] in zero_parameters: 1051 particle['width'] = 'ZERO' 1052 for pdg, particle in self['particle_dict'].items(): 1053 if particle['mass'] in zero_parameters: 1054 particle['mass'] = 'ZERO' 1055 if particle['width'] in zero_parameters: 1056 particle['width'] = 'ZERO' 1057 1058 # check if the parameters is still usefull: 1059 re_str = '|'.join(special_parameters) 1060 re_pat = re.compile(r'''\b(%s)\b''' % re_str) 1061 used = set() 1062 # check in coupling 1063 for name, coupling_list in self['couplings'].items(): 1064 for coupling in coupling_list: 1065 for use in re_pat.findall(coupling.expr): 1066 used.add(use) 1067 1068 # simplify the regular expression 1069 re_str = '|'.join([param for param in special_parameters 1070 if param not in used]) 1071 re_pat = re.compile(r'''\b(%s)\b''' % re_str) 1072 1073 param_info = {} 1074 # check in parameters 1075 for dep, param_list in self['parameters'].items(): 1076 for tag, parameter in enumerate(param_list): 1077 # update information concerning zero/one parameters 1078 if parameter.name in special_parameters: 1079 param_info[parameter.name]= {'dep': dep, 'tag': tag, 1080 'obj': parameter} 1081 continue 1082 1083 # Bypass all external parameter 1084 if isinstance(parameter, base_objects.ParamCardVariable): 1085 continue 1086 1087 # check the presence of zero/one parameter 1088 for use in re_pat.findall(parameter.expr): 1089 used.add(use) 1090 1091 # modify the object for those which are still used 1092 for param in used: 1093 data = self['parameters'][param_info[param]['dep']] 1094 data.remove(param_info[param]['obj']) 1095 tag = param_info[param]['tag'] 1096 data = self['parameters'][()] 1097 if param in zero_parameters: 1098 data.insert(0, base_objects.ModelVariable(param, '0.0', 'real')) 1099 else: 1100 data.insert(0, base_objects.ModelVariable(param, '1.0', 'real')) 1101 1102 # remove completely useless parameters 1103 for param in special_parameters: 1104 #by pass parameter still in use 1105 if param in used: 1106 logger_mod.debug('fix parameter value: %s' % param) 1107 continue 1108 logger_mod.debug('remove parameters: %s' % param) 1109 data = self['parameters'][param_info[param]['dep']] 1110 data.remove(param_info[param]['obj'])

1111

Source Code for Module models.import_ufo