madgraph.various.diagram

1 ################################################################################ 2 # 3 # Copyright (c) 2010 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 """Module for calculation of symmetries between diagrams, by 17 evaluating amp2 values for permutations of momenta.""" 18 19 from __future__ import division 20 21 import array 22 import copy 23 import fractions 24 import itertools 25 import logging 26 import math 27 import os 28 import re 29 import signal 30 31 import aloha.aloha_writers as aloha_writers 32 import aloha.create_aloha as create_aloha 33 34 import madgraph.iolibs.export_python as export_python 35 import madgraph.iolibs.group_subprocs as group_subprocs 36 import madgraph.iolibs.helas_call_writers as helas_call_writer 37 import models.import_ufo as import_ufo 38 import madgraph.iolibs.save_load_object as save_load_object 39 40 import madgraph.core.base_objects as base_objects 41 import madgraph.core.color_algebra as color 42 import madgraph.core.color_amp as color_amp 43 import madgraph.core.helas_objects as helas_objects 44 import madgraph.core.diagram_generation as diagram_generation 45 46 import madgraph.various.process_checks as process_checks 47 import madgraph.various.misc as misc 48 49 from madgraph import MG5DIR 50 51 import models.model_reader as model_reader 52 import aloha.template_files.wavefunctions as wavefunctions 53 from aloha.template_files.wavefunctions import \ 54 ixxxxx, oxxxxx, vxxxxx, sxxxxx 55 56 #=============================================================================== 57 # Logger for process_checks 58 #=============================================================================== 59 60 logger = logging.getLogger('madgraph.various.diagram_symmetry')

61 62 #=============================================================================== 63 # find_symmetry 64 #=============================================================================== 65 66 -def find_symmetry(matrix_element):

67 """Find symmetries between amplitudes by comparing diagram tags 68 for all the diagrams in the process. Identical diagram tags 69 correspond to different external particle permutations of the same 70 diagram. 71 72 Return list of positive number corresponding to number of 73 symmetric diagrams and negative numbers corresponding to the 74 equivalent diagram (for e+e->3a, get [6, -1, -1, -1, -1, -1]), 75 list of the corresponding permutations needed, and list of all 76 permutations of identical particles.""" 77 78 if isinstance(matrix_element, group_subprocs.SubProcessGroup): 79 return find_symmetry_subproc_group(matrix_element) 80 81 nexternal, ninitial = matrix_element.get_nexternal_ninitial() 82 83 # diagram_numbers is a list of all relevant diagram numbers 84 diagram_numbers = [] 85 # Prepare the symmetry vector with non-used amp2s (due to 86 # multiparticle vertices) 87 symmetry = [] 88 permutations = [] 89 ident_perms = [] 90 process = matrix_element.get('processes')[0] 91 base_model = process.get('model') 92 diagrams = matrix_element.get('diagrams') 93 base_diagrams = matrix_element.get_base_amplitude().get('diagrams') 94 min_vert = min([max(diag.get_vertex_leg_numbers()) for diag in diagrams]) 95 for diag in matrix_element.get('diagrams'): 96 diagram_numbers.append(diag.get('number')) 97 permutations.append(range(nexternal)) 98 if max(diag.get_vertex_leg_numbers()) > min_vert: 99 # Ignore any diagrams with 4-particle vertices 100 symmetry.append(0) 101 else: 102 symmetry.append(1) 103 104 # Check for matrix elements with no identical particles 105 if matrix_element.get("identical_particle_factor") == 1: 106 return symmetry, \ 107 permutations,\ 108 [range(nexternal)] 109 110 logger.info("Finding symmetric diagrams for process %s" % \ 111 matrix_element.get('processes')[0].nice_string().\ 112 replace("Process: ", "")) 113 114 # diagram_tags is a list of unique tags 115 diagram_tags = [] 116 # diagram_classes is a list of lists of diagram numbers belonging 117 # to the different classes 118 diagram_classes = [] 119 perms = [] 120 for diag, base_diagram in zip(diagrams, base_diagrams): 121 if any([vert > min_vert for vert in 122 diag.get_vertex_leg_numbers()]): 123 # Only 3-vertices allowed in configs.inc 124 continue 125 tag = diagram_generation.DiagramTag(base_diagram) 126 try: 127 ind = diagram_tags.index(tag) 128 except ValueError: 129 diagram_classes.append([diag.get('number')]) 130 perms.append([tag.get_external_numbers()]) 131 diagram_tags.append(tag) 132 else: 133 diagram_classes[ind].append(diag.get('number')) 134 perms[ind].append(tag.get_external_numbers()) 135 136 for inum, diag_number in enumerate(diagram_numbers): 137 if symmetry[inum] == 0: 138 continue 139 idx1 = [i for i, d in enumerate(diagram_classes) if \ 140 diag_number in d][0] 141 idx2 = diagram_classes[idx1].index(diag_number) 142 if idx2 == 0: 143 symmetry[inum] = len(diagram_classes[idx1]) 144 else: 145 symmetry[inum] = -diagram_classes[idx1][0] 146 # Order permutations according to how to reach the first perm 147 permutations[inum] = diagram_generation.DiagramTag.reorder_permutation(perms[idx1][idx2], 148 perms[idx1][0]) 149 # ident_perms ordered according to order of external momenta 150 perm = diagram_generation.DiagramTag.reorder_permutation(perms[idx1][0], 151 perms[idx1][idx2]) 152 if not perm in ident_perms: 153 ident_perms.append(perm) 154 155 return (symmetry, permutations, ident_perms)

156

157 -def find_symmetry_by_evaluation(matrix_element, evaluator, max_time = 600):

158 """Find symmetries between amplitudes by comparing the squared 159 amplitudes for all permutations of identical particles. 160 161 Return list of positive number corresponding to number of 162 symmetric diagrams and negative numbers corresponding to the 163 equivalent diagram (for e+e->3a, get [6, -1, -1, -1, -1, -1]), 164 list of the corresponding permutations needed, and list of all 165 permutations of identical particles. 166 max_time gives a cutoff time for finding symmetries (in s).""" 167 168 #if isinstance(matrix_element, group_subprocs.SubProcessGroup): 169 # return find_symmetry_subproc_group(matrix_element, evaluator, max_time) 170 171 assert isinstance(matrix_element, helas_objects.HelasMatrixElement) 172 173 # Exception class and routine to handle timeout 174 class TimeOutError(Exception): 175 pass

176 def handle_alarm(signum, frame): 177 raise TimeOutError 178 179 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 180 181 # Prepare the symmetry vector with non-used amp2s (due to 182 # multiparticle vertices) 183 symmetry = [] 184 for diag in matrix_element.get('diagrams'): 185 if max(diag.get_vertex_leg_numbers()) > 3: 186 # Ignore any diagrams with 4-particle vertices 187 symmetry.append(0) 188 else: 189 symmetry.append(1) 190 191 # Check for matrix elements with no identical particles 192 if matrix_element.get("identical_particle_factor") == 1: 193 return symmetry, \ 194 [range(nexternal)]*len(symmetry),\ 195 [range(nexternal)] 196 197 logger.info("Finding symmetric diagrams for process %s" % \ 198 matrix_element.get('processes')[0].nice_string().\ 199 replace("Process: ", "")) 200 201 process = matrix_element.get('processes')[0] 202 base_model = process.get('model') 203 equivalent_process = base_objects.Process({\ 204 'legs': base_objects.LegList([base_objects.Leg({ 205 'id': wf.get('pdg_code'), 206 'state': wf.get('leg_state')}) \ 207 for wf in matrix_element.get_external_wavefunctions()]), 208 'model': base_model}) 209 210 # Get phase space point 211 p, w_rambo = evaluator.get_momenta(equivalent_process) 212 213 # Check matrix element value for all permutations 214 amp2start = [] 215 final_states = [l.get('id') for l in \ 216 equivalent_process.get('legs')[ninitial:]] 217 nperm = 0 218 perms = [] 219 ident_perms = [] 220 221 # Set timeout for max_time 222 signal.signal(signal.SIGALRM, handle_alarm) 223 signal.alarm(max_time) 224 try: 225 for perm in itertools.permutations(range(ninitial, nexternal)): 226 if [equivalent_process.get('legs')[i].get('id') for i in perm] != \ 227 final_states: 228 # Non-identical particles permutated 229 continue 230 ident_perms.append([0,1]+list(perm)) 231 nperm += 1 232 new_p = p[:ninitial] + [p[i] for i in perm] 233 234 res = evaluator.evaluate_matrix_element(matrix_element, new_p) 235 if not res: 236 break 237 me_value, amp2 = res 238 # Make a list with (8-pos value, magnitude) to easily compare 239 amp2sum = sum(amp2) 240 amp2mag = [] 241 for a in amp2: 242 a = a*me_value/max(amp2sum, 1e-30) 243 if a > 0: 244 amp2mag.append(int(math.floor(math.log10(abs(a))))) 245 else: 246 amp2mag.append(0) 247 amp2 = [(int(a*10**(8-am)), am) for (a, am) in zip(amp2, amp2mag)] 248 249 if not perms: 250 # This is the first iteration - initialize lists 251 # Initiate symmetry with all 1:s 252 symmetry = [1 for i in range(len(amp2))] 253 # Store initial amplitudes 254 amp2start = amp2 255 # Initialize list of permutations 256 perms = [range(nexternal) for i in range(len(amp2))] 257 continue 258 259 for i, val in enumerate(amp2): 260 if val == (0,0): 261 # If amp2 is 0, just set symmetry to 0 262 symmetry[i] = 0 263 continue 264 # Only compare with diagrams below this one 265 if val in amp2start[:i]: 266 ind = amp2start.index(val) 267 # Replace if 1) this amp is unmatched (symmetry[i] > 0) or 268 # 2) this amp is matched but matched to an amp larger 269 # than ind 270 if symmetry[ind] > 0 and \ 271 (symmetry[i] > 0 or \ 272 symmetry[i] < 0 and -symmetry[i] > ind + 1): 273 symmetry[i] = -(ind+1) 274 perms[i] = [0, 1] + list(perm) 275 symmetry[ind] += 1 276 except TimeOutError: 277 # Symmetry canceled due to time limit 278 logger.warning("Cancel diagram symmetry - time exceeded") 279 280 # Stop the alarm since we're done with this process 281 signal.alarm(0) 282 283 return (symmetry, perms, ident_perms) 284

285 #=============================================================================== 286 # DiagramTag class to identify matrix elements 287 #=============================================================================== 288 289 -class IdentifySGConfigTag(diagram_generation.DiagramTag):

290 """DiagramTag daughter class to identify configs giving the same 291 config. Need to compare state, spin, mass, width, and color.""" 292 293 @staticmethod

294 - def link_from_leg(leg, model):

295 """Returns the end link for a leg needed to identify matrix 296 elements: ((leg numer, state, spin, self_antipart, mass, 297 width, color, decay and is_part), number).""" 298 299 part = model.get_particle(leg.get('id')) 300 301 state = 0 302 if not leg.get('state'): 303 # Distinguish identical initial state particles 304 state = leg.get('number') 305 306 return [((state, part.get('spin'), part.get('color'), 307 part.get('mass'), part.get('width')), 308 leg.get('number'))]

309 310 @staticmethod

311 - def vertex_id_from_vertex(vertex, last_vertex, model, ninitial):

312 """Returns the info needed to identify matrix elements: 313 interaction color, lorentz, coupling, and wavefunction 314 spin, self_antipart, mass, width, color, decay and 315 is_part. Note that is_part needs to be flipped if we move the 316 final vertex around.""" 317 318 inter = model.get_interaction(vertex.get('id')) 319 ret_list = 0 320 321 if last_vertex: 322 return (ret_list,) 323 else: 324 part = model.get_particle(vertex.get('legs')[-1].get('id')) 325 return ((part.get('color'), 326 part.get('mass'), part.get('width')), 327 ret_list)

328 329 @staticmethod

330 - def flip_vertex(new_vertex, old_vertex):

331 """Move the wavefunction part of vertex id appropriately""" 332 333 if len(new_vertex) == 1 and len(old_vertex) == 2: 334 # We go from a last link to next-to-last link - add propagator info 335 return (old_vertex[0],new_vertex[0]) 336 elif len(new_vertex) == 2 and len(old_vertex) == 1: 337 # We go from next-to-last link to last link - remove propagator info 338 return (new_vertex[1],) 339 # We should not get here 340 raise diagram_generation.DiagramTag.DiagramTagError, \ 341 "Error in IdentifyConfigTag, wrong setup of vertices in link."

342

343 -def find_symmetry_subproc_group(subproc_group):

344 """Find symmetric configs by directly comparing the configurations 345 using IdentifySGConfigTag.""" 346 347 assert isinstance(subproc_group, group_subprocs.SubProcessGroup),\ 348 "Argument to find_symmetry_subproc_group has to be SubProcessGroup" 349 350 # diagram_numbers is a list of all relevant diagram numbers 351 diagram_numbers = [] 352 # Prepare the symmetry vector with non-used amp2s (due to 353 # multiparticle vertices) 354 symmetry = [] 355 permutations = [] 356 diagrams = subproc_group.get('mapping_diagrams') 357 nexternal, ninitial = \ 358 subproc_group.get('matrix_elements')[0].get_nexternal_ninitial() 359 model = subproc_group.get('matrix_elements')[0].get('processes')[0].\ 360 get('model') 361 min_vert = min([max(diag.get_vertex_leg_numbers()) for diag in diagrams]) 362 363 for idiag,diag in enumerate(diagrams): 364 diagram_numbers.append(idiag+1) 365 permutations.append(range(nexternal)) 366 if max(diag.get_vertex_leg_numbers()) > min_vert: 367 # Ignore any diagrams with 4-particle vertices 368 symmetry.append(0) 369 else: 370 symmetry.append(1) 371 372 logger.info("Finding symmetric diagrams for subprocess group %s" % \ 373 subproc_group.get('name')) 374 375 # diagram_tags is a list of unique tags 376 diagram_tags = [] 377 # diagram_classes is a list of lists of diagram numbers belonging 378 # to the different classes 379 diagram_classes = [] 380 perms = [] 381 for idiag, diag in enumerate(diagrams): 382 if any([vert > min_vert for vert in 383 diag.get_vertex_leg_numbers()]): 384 # Only include vertices up to min_vert 385 continue 386 tag = IdentifySGConfigTag(diag, model) 387 try: 388 ind = diagram_tags.index(tag) 389 except ValueError: 390 diagram_classes.append([idiag + 1]) 391 perms.append([tag.get_external_numbers()]) 392 diagram_tags.append(tag) 393 else: 394 diagram_classes[ind].append(idiag + 1) 395 perms[ind].append(tag.get_external_numbers()) 396 397 for inum, diag_number in enumerate(diagram_numbers): 398 if symmetry[inum] == 0: 399 continue 400 idx1 = [i for i, d in enumerate(diagram_classes) if \ 401 diag_number in d][0] 402 idx2 = diagram_classes[idx1].index(diag_number) 403 # Note that for subproc groups, we want symfact to be 1 404 if idx2 > 0: 405 symmetry[inum] = -diagram_classes[idx1][0] 406 # Order permutations according to how to reach the first perm 407 permutations[inum] = diagram_generation.DiagramTag.reorder_permutation(perms[idx1][idx2], 408 perms[idx1][0]) 409 return (symmetry, permutations, [permutations[0]])

410

411 412 -def old_find_symmetry_subproc_group(subproc_group):

413 """Find symmetries between the configs in the subprocess group. 414 For each config, find all matrix elements with maximum identical 415 particle factor. Then take minimal set of these matrix elements, 416 and determine symmetries based on these.""" 417 418 assert isinstance(subproc_group, group_subprocs.SubProcessGroup),\ 419 "Argument to find_symmetry_subproc_group has to be SubProcessGroup" 420 421 matrix_elements = subproc_group.get('matrix_elements') 422 423 contributing_mes, me_config_dict = \ 424 find_matrix_elements_for_configs(subproc_group) 425 426 nexternal, ninitial = matrix_elements[0].get_nexternal_ninitial() 427 428 all_symmetry = {} 429 all_perms = {} 430 431 for me_number in contributing_mes: 432 diagram_config_map = dict([(i,n) for i,n in \ 433 enumerate(subproc_group.get('diagram_maps')[me_number]) \ 434 if n > 0]) 435 symmetry, perms, ident_perms = find_symmetry(matrix_elements[me_number]) 436 437 # Go through symmetries and remove those for any diagrams 438 # where this ME is not supposed to contribute 439 for isym, sym_config in enumerate(symmetry): 440 if sym_config == 0 or isym not in diagram_config_map: 441 continue 442 config = diagram_config_map[isym] 443 if config not in me_config_dict[me_number] or \ 444 sym_config < 0 and diagram_config_map[-sym_config-1] not in \ 445 me_config_dict[me_number]: 446 symmetry[isym] = 1 447 perms[isym]=range(nexternal) 448 if sym_config < 0 and diagram_config_map[-sym_config-1] in \ 449 me_config_dict[me_number]: 450 symmetry[-sym_config-1] -= 1 451 452 # Now update the maps all_symmetry and all_perms 453 for isym, (perm, sym_config) in enumerate(zip(perms, symmetry)): 454 if sym_config in [0,1] or isym not in diagram_config_map: 455 continue 456 config = diagram_config_map[isym] 457 458 all_perms[config] = perm 459 460 if sym_config > 0: 461 all_symmetry[config] = sym_config 462 else: 463 all_symmetry[config] = -diagram_config_map[-sym_config-1] 464 465 # Fill up all_symmetry and all_perms also for configs that have no symmetry 466 for iconf in range(len(subproc_group.get('mapping_diagrams'))): 467 all_symmetry.setdefault(iconf+1, 1) 468 all_perms.setdefault(iconf+1, range(nexternal)) 469 # Since we don't want to multiply by symmetry factor here, set to 1 470 if all_symmetry[iconf+1] > 1: 471 all_symmetry[iconf+1] = 1 472 473 symmetry = [all_symmetry[key] for key in sorted(all_symmetry.keys())] 474 perms = [all_perms[key] for key in sorted(all_perms.keys())] 475 476 return symmetry, perms, [perms[0]]

477

478 479 -def find_matrix_elements_for_configs(subproc_group):

480 """For each config, find all matrix elements with maximum identical 481 particle factor. Then take minimal set of these matrix elements.""" 482 483 matrix_elements = subproc_group.get('matrix_elements') 484 485 n_mes = len(matrix_elements) 486 487 me_config_dict = {} 488 489 # Find the MEs with maximum ident factor corresponding to each config. 490 # Only include MEs with identical particles (otherwise no contribution) 491 for iconf, diagram_list in \ 492 enumerate(subproc_group.get('diagrams_for_configs')): 493 # Check if any diagrams contribute to config 494 if set(diagram_list) == set([0]): 495 continue 496 # Add list of MEs with maximum ident factor contributing to this config 497 max_ident = max([matrix_elements[i].get('identical_particle_factor') \ 498 for i in range(n_mes) if diagram_list[i] > 0]) 499 max_mes = [i for i in range(n_mes) if \ 500 matrix_elements[i].get('identical_particle_factor') == \ 501 max_ident and diagram_list[i] > 0 and max_ident > 1] 502 for me in max_mes: 503 me_config_dict.setdefault(me, [iconf+1]).append(iconf + 1) 504 505 # Make set of the configs 506 for me in me_config_dict: 507 me_config_dict[me] = sorted(set(me_config_dict[me])) 508 509 # Sort MEs according to 1) ident factor, 2) number of configs they 510 # contribute to 511 def me_sort(me1, me2): 512 return (matrix_elements[me2].get('identical_particle_factor') \ 513 - matrix_elements[me1].get('identical_particle_factor'))\ 514 or (len(me_config_dict[me2]) - len(me_config_dict[me1]))

515 516 sorted_mes = sorted([me for me in me_config_dict], me_sort) 517 518 # Reduce to minimal set of matrix elements 519 latest_me = 0 520 checked_configs = [] 521 while latest_me < len(sorted_mes): 522 checked_configs.extend(me_config_dict[sorted_mes[latest_me]]) 523 for me in sorted_mes[latest_me+1:]: 524 me_config_dict[me] = [conf for conf in me_config_dict[me] if \ 525 conf not in checked_configs] 526 if me_config_dict[me] == []: 527 del me_config_dict[me] 528 # Re-sort MEs 529 sorted_mes = sorted([me for me in me_config_dict], me_sort) 530 latest_me += 1 531 532 return sorted_mes, me_config_dict 533

Source Code for Module madgraph.various.diagram_symmetry