madgraph.various.process

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 """Several different checks for processes (and hence models): 16 permutation tests, gauge invariance tests, lorentz invariance 17 tests. Also class for evaluation of Python matrix elements, 18 MatrixElementEvaluator.""" 19 20 from __future__ import division 21 22 import array 23 import copy 24 import fractions 25 import itertools 26 import logging 27 import math 28 import os 29 import re 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.helas_call_writers as helas_call_writers 36 import models.import_ufo as import_ufo 37 import madgraph.iolibs.save_load_object as save_load_object 38 39 import madgraph.core.base_objects as base_objects 40 import madgraph.core.color_algebra as color 41 import madgraph.core.color_amp as color_amp 42 import madgraph.core.helas_objects as helas_objects 43 import madgraph.core.diagram_generation as diagram_generation 44 45 import madgraph.various.rambo as rambo 46 import madgraph.various.misc as misc 47 48 49 from madgraph import MG5DIR, InvalidCmd 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, txxxxx 55 56 #=============================================================================== 57 # Logger for process_checks 58 #=============================================================================== 59 60 logger = logging.getLogger('madgraph.various.process_checks') 61 62 #=============================================================================== 63 # Helper class MatrixElementEvaluator 64 #===============================================================================

65 -class MatrixElementEvaluator(object):

66 """Class taking care of matrix element evaluation, storing 67 relevant quantities for speedup.""" 68

69 - def __init__(self, model, param_card = None, 70 auth_skipping = False, reuse = True):

71 """Initialize object with stored_quantities, helas_writer, 72 model, etc. 73 auth_skipping = True means that any identical matrix element will be 74 evaluated only once 75 reuse = True means that the matrix element corresponding to a 76 given process can be reused (turn off if you are using 77 different models for the same process)""" 78 79 # Writer for the Python matrix elements 80 self.helas_writer = helas_call_writers.PythonUFOHelasCallWriter(model) 81 82 # Read a param_card and calculate couplings 83 self.full_model = model_reader.ModelReader(model) 84 self.full_model.set_parameters_and_couplings(param_card) 85 86 self.auth_skipping = auth_skipping 87 self.reuse = reuse 88 self.store_aloha = [] 89 self.stored_quantities = {}

90 91 #=============================================================================== 92 # Helper function evaluate_matrix_element 93 #===============================================================================

94 - def evaluate_matrix_element(self, matrix_element, p=None, full_model=None, 95 gauge_check=False, auth_skipping=None, output='m2'):

96 """Calculate the matrix element and evaluate it for a phase space point 97 output is either m2, amp, jamp 98 """ 99 if full_model: 100 self.full_model = full_model 101 102 process = matrix_element.get('processes')[0] 103 model = process.get('model') 104 105 if "matrix_elements" not in self.stored_quantities: 106 self.stored_quantities['matrix_elements'] = [] 107 matrix_methods = {} 108 109 if self.reuse and "Matrix_%s" % process.shell_string() in globals() and p: 110 # Evaluate the matrix element for the momenta p 111 matrix = eval("Matrix_%s()" % process.shell_string()) 112 me_value = matrix.smatrix(p, self.full_model) 113 if output == "m2": 114 return matrix.smatrix(p, self.full_model), matrix.amp2 115 else: 116 m2 = matrix.smatrix(p, self.full_model) 117 return {'m2': m2, output:getattr(matrix, output)} 118 119 if (auth_skipping or self.auth_skipping) and matrix_element in \ 120 self.stored_quantities['matrix_elements']: 121 # Exactly the same matrix element has been tested 122 logger.info("Skipping %s, " % process.nice_string() + \ 123 "identical matrix element already tested" \ 124 ) 125 return None 126 127 self.stored_quantities['matrix_elements'].append(matrix_element) 128 129 # Create an empty color basis, and the list of raw 130 # colorize objects (before simplification) associated 131 # with amplitude 132 if "list_colorize" not in self.stored_quantities: 133 self.stored_quantities["list_colorize"] = [] 134 if "list_color_basis" not in self.stored_quantities: 135 self.stored_quantities["list_color_basis"] = [] 136 if "list_color_matrices" not in self.stored_quantities: 137 self.stored_quantities["list_color_matrices"] = [] 138 139 col_basis = color_amp.ColorBasis() 140 new_amp = matrix_element.get_base_amplitude() 141 matrix_element.set('base_amplitude', new_amp) 142 colorize_obj = col_basis.create_color_dict_list(new_amp) 143 144 try: 145 # If the color configuration of the ME has 146 # already been considered before, recycle 147 # the information 148 col_index = self.stored_quantities["list_colorize"].index(colorize_obj) 149 except ValueError: 150 # If not, create color basis and color 151 # matrix accordingly 152 self.stored_quantities['list_colorize'].append(colorize_obj) 153 col_basis.build() 154 self.stored_quantities['list_color_basis'].append(col_basis) 155 col_matrix = color_amp.ColorMatrix(col_basis) 156 self.stored_quantities['list_color_matrices'].append(col_matrix) 157 col_index = -1 158 159 # Set the color for the matrix element 160 matrix_element.set('color_basis', 161 self.stored_quantities['list_color_basis'][col_index]) 162 matrix_element.set('color_matrix', 163 self.stored_quantities['list_color_matrices'][col_index]) 164 165 # Create the needed aloha routines 166 if "used_lorentz" not in self.stored_quantities: 167 self.stored_quantities["used_lorentz"] = [] 168 169 me_used_lorentz = set(matrix_element.get_used_lorentz()) 170 me_used_lorentz = [lorentz for lorentz in me_used_lorentz \ 171 if lorentz not in self.store_aloha] 172 173 aloha_model = create_aloha.AbstractALOHAModel(model.get('name')) 174 aloha_model.compute_subset(me_used_lorentz) 175 176 # Write out the routines in Python 177 aloha_routines = [] 178 for routine in aloha_model.values(): 179 aloha_routines.append(routine.write(output_dir = None, 180 mode='mg5', 181 language = 'Python')) 182 for routine in aloha_model.external_routines: 183 aloha_routines.append( 184 open(aloha_model.locate_external(routine, 'Python')).read()) 185 186 187 # Define the routines to be available globally 188 for routine in aloha_routines: 189 exec(routine, globals()) 190 191 # Add the defined Aloha routines to used_lorentz 192 self.store_aloha.extend(me_used_lorentz) 193 194 # Export the matrix element to Python calls 195 exporter = export_python.ProcessExporterPython(matrix_element, 196 self.helas_writer) 197 try: 198 matrix_methods = exporter.get_python_matrix_methods(\ 199 gauge_check=gauge_check) 200 except helas_call_writers.HelasWriterError, error: 201 logger.info(error) 202 return None 203 204 if self.reuse: 205 # Define the routines (globally) 206 exec(matrix_methods[process.shell_string()], globals()) 207 else: 208 # Define the routines (locally is enough) 209 exec(matrix_methods[process.shell_string()]) 210 211 # Generate phase space point to use 212 if not p: 213 p, w_rambo = self.get_momenta(process) 214 215 # Evaluate the matrix element for the momenta p 216 exec("data = Matrix_%s()" % process.shell_string()) 217 if output == "m2": 218 return data.smatrix(p, self.full_model), data.amp2 219 else: 220 m2 = data.smatrix(p, self.full_model) 221 return {'m2': m2, output:getattr(data, output)}

222 223 #=============================================================================== 224 # Helper function get_momenta 225 #===============================================================================

226 - def get_momenta(self, process, energy = 1000.):

227 """Get a point in phase space for the external states in the given 228 process, with the CM energy given. The incoming particles are 229 assumed to be oriented along the z axis, with particle 1 along the 230 positive z axis.""" 231 232 if not (isinstance(process, base_objects.Process) and \ 233 isinstance(energy, float)): 234 raise rambo.RAMBOError, "Not correct type for arguments to get_momenta" 235 236 sorted_legs = sorted(process.get('legs'), lambda l1, l2:\ 237 l1.get('number') - l2.get('number')) 238 239 nincoming = len([leg for leg in sorted_legs if leg.get('state') == False]) 240 nfinal = len(sorted_legs) - nincoming 241 242 # Find masses of particles 243 mass_strings = [self.full_model.get_particle(l.get('id')).get('mass') \ 244 for l in sorted_legs] 245 mass = [abs(self.full_model.get('parameter_dict')[m]) for m in mass_strings] 246 247 # Make sure energy is large enough for incoming and outgoing particles 248 energy = max(energy, sum(mass[:nincoming]) + 200., 249 sum(mass[nincoming:]) + 200.) 250 251 e2 = energy**2 252 m1 = mass[0] 253 254 p = [] 255 256 masses = rambo.FortranList(nfinal) 257 for i in range(nfinal): 258 masses[i+1] = mass[nincoming + i] 259 260 if nincoming == 1: 261 262 # Momenta for the incoming particle 263 p.append([m1, 0., 0., 0.]) 264 265 p_rambo, w_rambo = rambo.RAMBO(nfinal, m1, masses) 266 267 # Reorder momenta from px,py,pz,E to E,px,py,pz scheme 268 for i in range(1, nfinal+1): 269 momi = [p_rambo[(4,i)], p_rambo[(1,i)], 270 p_rambo[(2,i)], p_rambo[(3,i)]] 271 p.append(momi) 272 273 return p, w_rambo 274 275 if nincoming != 2: 276 raise rambo.RAMBOError('Need 1 or 2 incoming particles') 277 278 if nfinal == 1: 279 energy = masses[0] 280 281 m2 = mass[1] 282 283 mom = math.sqrt((e2**2 - 2*e2*m1**2 + m1**4 - 2*e2*m2**2 - \ 284 2*m1**2*m2**2 + m2**4) / (4*e2)) 285 e1 = math.sqrt(mom**2+m1**2) 286 e2 = math.sqrt(mom**2+m2**2) 287 # Set momenta for incoming particles 288 p.append([e1, 0., 0., mom]) 289 p.append([e2, 0., 0., -mom]) 290 291 if nfinal == 1: 292 p.append([energy, 0., 0., 0.]) 293 return p, 1. 294 295 p_rambo, w_rambo = rambo.RAMBO(nfinal, energy, masses) 296 297 # Reorder momenta from px,py,pz,E to E,px,py,pz scheme 298 for i in range(1, nfinal+1): 299 momi = [p_rambo[(4,i)], p_rambo[(1,i)], 300 p_rambo[(2,i)], p_rambo[(3,i)]] 301 p.append(momi) 302 303 return p, w_rambo

304 305 #=============================================================================== 306 # Global helper function run_multiprocs 307 #===============================================================================

308 -def run_multiprocs_no_crossings(function, multiprocess, stored_quantities, 309 *args):

310 """A wrapper function for running an iteration of a function over 311 a multiprocess, without having to first create a process list 312 (which makes a big difference for very large multiprocesses. 313 stored_quantities is a dictionary for any quantities that we want 314 to reuse between runs.""" 315 316 model = multiprocess.get('model') 317 isids = [leg.get('ids') for leg in multiprocess.get('legs') \ 318 if not leg.get('state')] 319 fsids = [leg.get('ids') for leg in multiprocess.get('legs') \ 320 if leg.get('state')] 321 # Create dictionary between isids and antiids, to speed up lookup 322 id_anti_id_dict = {} 323 for id in set(tuple(sum(isids+fsids, []))): 324 id_anti_id_dict[id] = model.get_particle(id).get_anti_pdg_code() 325 id_anti_id_dict[model.get_particle(id).get_anti_pdg_code()] = id 326 327 sorted_ids = [] 328 results = [] 329 for is_prod in apply(itertools.product, isids): 330 for fs_prod in apply(itertools.product, fsids): 331 332 # Check if we have already checked the process 333 if check_already_checked(is_prod, fs_prod, sorted_ids, 334 multiprocess, model, id_anti_id_dict): 335 continue 336 337 # Generate process based on the selected ids 338 process = base_objects.Process({\ 339 'legs': base_objects.LegList(\ 340 [base_objects.Leg({'id': id, 'state':False}) for \ 341 id in is_prod] + \ 342 [base_objects.Leg({'id': id, 'state':True}) for \ 343 id in fs_prod]), 344 'model':multiprocess.get('model'), 345 'id': multiprocess.get('id'), 346 'orders': multiprocess.get('orders'), 347 'required_s_channels': \ 348 multiprocess.get('required_s_channels'), 349 'forbidden_s_channels': \ 350 multiprocess.get('forbidden_s_channels'), 351 'forbidden_particles': \ 352 multiprocess.get('forbidden_particles'), 353 'is_decay_chain': \ 354 multiprocess.get('is_decay_chain'), 355 'overall_orders': \ 356 multiprocess.get('overall_orders')}) 357 358 result = function(process, stored_quantities, *args) 359 360 if result: 361 results.append(result) 362 363 return results

364 365 #=============================================================================== 366 # Helper function check_already_checked 367 #===============================================================================

368 -def check_already_checked(is_ids, fs_ids, sorted_ids, process, model, 369 id_anti_id_dict = {}):

370 """Check if process already checked, if so return True, otherwise add 371 process and antiprocess to sorted_ids.""" 372 373 # Check if process is already checked 374 if id_anti_id_dict: 375 is_ids = [id_anti_id_dict[id] for id in \ 376 is_ids] 377 else: 378 is_ids = [model.get_particle(id).get_anti_pdg_code() for id in \ 379 is_ids] 380 381 ids = array.array('i', sorted(is_ids + list(fs_ids)) + \ 382 [process.get('id')]) 383 384 if ids in sorted_ids: 385 # We have already checked (a crossing of) this process 386 return True 387 388 # Add this process to tested_processes 389 sorted_ids.append(ids) 390 391 # Skip adding antiprocess below, since might be relevant too 392 return False

393 394 #=============================================================================== 395 # check_processes 396 #===============================================================================

397 -def check_processes(processes, param_card = None, quick = []):

398 """Check processes by generating them with all possible orderings 399 of particles (which means different diagram building and Helas 400 calls), and comparing the resulting matrix element values.""" 401 402 if isinstance(processes, base_objects.ProcessDefinition): 403 # Generate a list of unique processes 404 # Extract IS and FS ids 405 multiprocess = processes 406 model = multiprocess.get('model') 407 408 # Initialize matrix element evaluation 409 evaluator = MatrixElementEvaluator(model, 410 auth_skipping = True, reuse = False) 411 results = run_multiprocs_no_crossings(check_process, 412 multiprocess, 413 evaluator, 414 quick) 415 416 if "used_lorentz" not in evaluator.stored_quantities: 417 evaluator.stored_quantities["used_lorentz"] = [] 418 return results, evaluator.stored_quantities["used_lorentz"] 419 420 elif isinstance(processes, base_objects.Process): 421 processes = base_objects.ProcessList([processes]) 422 elif isinstance(processes, base_objects.ProcessList): 423 pass 424 else: 425 raise InvalidCmd("processes is of non-supported format") 426 427 if not processes: 428 raise InvalidCmd("No processes given") 429 430 model = processes[0].get('model') 431 432 # Initialize matrix element evaluation 433 evaluator = MatrixElementEvaluator(model, param_card, 434 auth_skipping = True, reuse = False) 435 436 # Keep track of tested processes, matrix elements, color and already 437 # initiated Lorentz routines, to reuse as much as possible 438 sorted_ids = [] 439 comparison_results = [] 440 441 # Check process by process 442 for process in processes: 443 444 # Check if we already checked process 445 if check_already_checked([l.get('id') for l in process.get('legs') if \ 446 not l.get('state')], 447 [l.get('id') for l in process.get('legs') if \ 448 l.get('state')], 449 sorted_ids, process, model): 450 continue 451 # Get process result 452 res = check_process(process, evaluator, quick) 453 if res: 454 comparison_results.append(res) 455 456 if "used_lorentz" not in evaluator.stored_quantities: 457 evaluator.stored_quantities["used_lorentz"] = [] 458 return comparison_results, evaluator.stored_quantities["used_lorentz"]

459 460

461 -def check_process(process, evaluator, quick):

462 """Check the helas calls for a process by generating the process 463 using all different permutations of the process legs (or, if 464 quick, use a subset of permutations), and check that the matrix 465 element is invariant under this.""" 466 467 model = process.get('model') 468 469 # Ensure that leg numbers are set 470 for i, leg in enumerate(process.get('legs')): 471 leg.set('number', i+1) 472 473 logger.info("Checking %s" % \ 474 process.nice_string().replace('Process', 'process')) 475 476 process_matrix_elements = [] 477 478 # For quick checks, only test twp permutations with leg "1" in 479 # each position 480 if quick: 481 leg_positions = [[] for leg in process.get('legs')] 482 quick = range(1,len(process.get('legs')) + 1) 483 484 values = [] 485 486 # Now, generate all possible permutations of the legs 487 for legs in itertools.permutations(process.get('legs')): 488 489 order = [l.get('number') for l in legs] 490 491 if quick: 492 found_leg = True 493 for num in quick: 494 # Only test one permutation for each position of the 495 # specified legs 496 leg_position = legs.index([l for l in legs if \ 497 l.get('number') == num][0]) 498 499 if not leg_position in leg_positions[num-1]: 500 found_leg = False 501 leg_positions[num-1].append(leg_position) 502 503 if found_leg: 504 continue 505 506 legs = base_objects.LegList(legs) 507 508 if order != range(1,len(legs) + 1): 509 logger.info("Testing permutation: %s" % \ 510 order) 511 512 # Generate a process with these legs 513 newproc = copy.copy(process) 514 newproc.set('legs', legs) 515 516 # Generate the amplitude for this process 517 try: 518 amplitude = diagram_generation.Amplitude(newproc) 519 except InvalidCmd: 520 result=False 521 else: 522 result = amplitude.get('diagrams') 523 524 if not result: 525 # This process has no diagrams; go to next process 526 logging.info("No diagrams for %s" % \ 527 process.nice_string().replace('Process', 'process')) 528 break 529 530 if order == range(1,len(legs) + 1): 531 # Generate phase space point to use 532 p, w_rambo = evaluator.get_momenta(process) 533 534 # Generate the HelasMatrixElement for the process 535 matrix_element = helas_objects.HelasMatrixElement(amplitude, 536 gen_color=False) 537 538 if matrix_element in process_matrix_elements: 539 # Exactly the same matrix element has been tested 540 # for other permutation of same process 541 continue 542 543 process_matrix_elements.append(matrix_element) 544 545 res = evaluator.evaluate_matrix_element(matrix_element, p = p) 546 if res == None: 547 break 548 549 values.append(res[0]) 550 551 # Check if we failed badly (1% is already bad) - in that 552 # case done for this process 553 if abs(max(values)) + abs(min(values)) > 0 and \ 554 2 * abs(max(values) - min(values)) / \ 555 (abs(max(values)) + abs(min(values))) > 0.01: 556 break 557 558 # Check if process was interrupted 559 if not values: 560 return None 561 562 # Done with this process. Collect values, and store 563 # process and momenta 564 diff = 0 565 if abs(max(values)) + abs(min(values)) > 0: 566 diff = 2* abs(max(values) - min(values)) / \ 567 (abs(max(values)) + abs(min(values))) 568 569 passed = diff < 1.e-8 570 571 return {"process": process, 572 "momenta": p, 573 "values": values, 574 "difference": diff, 575 "passed": passed}

576 577

578 -def output_comparisons(comparison_results):

579 """Present the results of a comparison in a nice list format 580 mode short: return the number of fail process 581 """ 582 583 proc_col_size = 17 584 585 for proc in comparison_results: 586 if len(proc['process'].base_string()) + 1 > proc_col_size: 587 proc_col_size = len(proc['process'].base_string()) + 1 588 589 col_size = 17 590 591 pass_proc = 0 592 fail_proc = 0 593 no_check_proc = 0 594 595 failed_proc_list = [] 596 no_check_proc_list = [] 597 598 res_str = fixed_string_length("Process", proc_col_size) + \ 599 fixed_string_length("Min element", col_size) + \ 600 fixed_string_length("Max element", col_size) + \ 601 fixed_string_length("Relative diff.", col_size) + \ 602 "Result" 603 604 for result in comparison_results: 605 proc = result['process'].base_string() 606 values = result['values'] 607 608 if len(values) <= 1: 609 res_str += '\n' + fixed_string_length(proc, proc_col_size) + \ 610 " * No permutations, process not checked *" 611 no_check_proc += 1 612 no_check_proc_list.append(result['process'].nice_string()) 613 continue 614 615 passed = result['passed'] 616 617 res_str += '\n' + fixed_string_length(proc, proc_col_size) + \ 618 fixed_string_length("%1.10e" % min(values), col_size) + \ 619 fixed_string_length("%1.10e" % max(values), col_size) + \ 620 fixed_string_length("%1.10e" % result['difference'], 621 col_size) 622 if passed: 623 pass_proc += 1 624 res_str += "Passed" 625 else: 626 fail_proc += 1 627 failed_proc_list.append(result['process'].nice_string()) 628 res_str += "Failed" 629 630 res_str += "\nSummary: %i/%i passed, %i/%i failed" % \ 631 (pass_proc, pass_proc + fail_proc, 632 fail_proc, pass_proc + fail_proc) 633 634 if fail_proc != 0: 635 res_str += "\nFailed processes: %s" % ', '.join(failed_proc_list) 636 if no_check_proc != 0: 637 res_str += "\nNot checked processes: %s" % ', '.join(no_check_proc_list) 638 639 return res_str

640

641 -def fixed_string_length(mystr, length):

642 """Helper function to fix the length of a string by cutting it 643 or adding extra space.""" 644 645 if len(mystr) > length: 646 return mystr[0:length] 647 else: 648 return mystr + " " * (length - len(mystr))

649 650 651 #=============================================================================== 652 # check_gauge 653 #===============================================================================

654 -def check_gauge(processes, param_card = None):

655 """Check gauge invariance of the processes by using the BRS check. 656 For one of the massless external bosons (e.g. gluon or photon), 657 replace the polarization vector (epsilon_mu) with its momentum (p_mu) 658 """ 659 660 if isinstance(processes, base_objects.ProcessDefinition): 661 # Generate a list of unique processes 662 # Extract IS and FS ids 663 multiprocess = processes 664 665 model = multiprocess.get('model') 666 667 # Initialize matrix element evaluation 668 evaluator = MatrixElementEvaluator(model, param_card, 669 auth_skipping = True, reuse = False) 670 671 # Set all widths to zero for gauge check 672 for particle in evaluator.full_model.get('particles'): 673 if particle.get('width') != 'ZERO': 674 evaluator.full_model.get('parameter_dict')[particle.get('width')] = 0. 675 676 return run_multiprocs_no_crossings(check_gauge_process, 677 multiprocess, 678 evaluator) 679 680 elif isinstance(processes, base_objects.Process): 681 processes = base_objects.ProcessList([processes]) 682 elif isinstance(processes, base_objects.ProcessList): 683 pass 684 else: 685 raise InvalidCmd("processes is of non-supported format") 686 687 assert processes, "No processes given" 688 689 model = processes[0].get('model') 690 691 # Initialize matrix element evaluation 692 evaluator = MatrixElementEvaluator(model, param_card, 693 auth_skipping = True, reuse = False) 694 695 comparison_results = [] 696 697 # For each process, make sure we have set up leg numbers: 698 for process in processes: 699 # Check if we already checked process 700 #if check_already_checked([l.get('id') for l in process.get('legs') if \ 701 # not l.get('state')], 702 ## [l.get('id') for l in process.get('legs') if \ 703 # l.get('state')], 704 # sorted_ids, process, model): 705 # continue 706 707 # Get process result 708 result = check_gauge_process(process, evaluator) 709 if result: 710 comparison_results.append(result) 711 712 return comparison_results

713 714

715 -def check_gauge_process(process, evaluator):

716 """Check gauge invariance for the process, unless it is already done.""" 717 718 model = process.get('model') 719 720 # Check that there are massless vector bosons in the process 721 found_gauge = False 722 for i, leg in enumerate(process.get('legs')): 723 part = model.get_particle(leg.get('id')) 724 if part.get('spin') == 3 and part.get('mass').lower() == 'zero': 725 found_gauge = True 726 break 727 728 if not found_gauge: 729 # This process can't be checked 730 return None 731 732 for i, leg in enumerate(process.get('legs')): 733 leg.set('number', i+1) 734 735 logger.info("Checking gauge %s" % \ 736 process.nice_string().replace('Process', 'process')) 737 738 legs = process.get('legs') 739 # Generate a process with these legs 740 # Generate the amplitude for this process 741 try: 742 amplitude = diagram_generation.Amplitude(process) 743 except InvalidCmd: 744 logging.info("No diagrams for %s" % \ 745 process.nice_string().replace('Process', 'process')) 746 return None 747 748 if not amplitude.get('diagrams'): 749 # This process has no diagrams; go to next process 750 logging.info("No diagrams for %s" % \ 751 process.nice_string().replace('Process', 'process')) 752 return None 753 754 # Generate the HelasMatrixElement for the process 755 matrix_element = helas_objects.HelasMatrixElement(amplitude, 756 gen_color = False) 757 758 brsvalue = evaluator.evaluate_matrix_element(matrix_element, gauge_check = True, 759 output='jamp') 760 761 762 matrix_element = helas_objects.HelasMatrixElement(amplitude, 763 gen_color = False) 764 mvalue = evaluator.evaluate_matrix_element(matrix_element, gauge_check = False, 765 output='jamp') 766 767 if mvalue and mvalue['m2']: 768 return {'process':process.base_string(),'value':mvalue,'brs':brsvalue}

769

770 -def output_gauge(comparison_results, output='text'):

771 """Present the results of a comparison in a nice list format""" 772 773 proc_col_size = 17 774 for one_comp in comparison_results: 775 proc = one_comp['process'] 776 mvalue = one_comp['value'] 777 brsvalue = one_comp['brs'] 778 if len(proc) + 1 > proc_col_size: 779 proc_col_size = len(proc) + 1 780 781 col_size = 17 782 783 pass_proc = 0 784 fail_proc = 0 785 786 failed_proc_list = [] 787 no_check_proc_list = [] 788 789 res_str = fixed_string_length("Process", proc_col_size) + \ 790 fixed_string_length("matrix", col_size) + \ 791 fixed_string_length("BRS", col_size) + \ 792 fixed_string_length("ratio", col_size) + \ 793 "Result" 794 795 for one_comp in comparison_results: 796 proc = one_comp['process'] 797 mvalue = one_comp['value'] 798 brsvalue = one_comp['brs'] 799 ratio = (brsvalue['m2']/mvalue['m2']) 800 res_str += '\n' + fixed_string_length(proc, proc_col_size) + \ 801 fixed_string_length("%1.10e" % mvalue['m2'], col_size)+ \ 802 fixed_string_length("%1.10e" % brsvalue['m2'], col_size)+ \ 803 fixed_string_length("%1.10e" % ratio, col_size) 804 805 if ratio > 1e-12: 806 fail_proc += 1 807 proc_succeed = False 808 failed_proc_list.append(proc) 809 res_str += "Failed" 810 else: 811 pass_proc += 1 812 proc_succeed = True 813 res_str += "Passed" 814 815 #check all the JAMP 816 # loop over jamp 817 for k in range(len(mvalue['jamp'][0])): 818 m_sum = 0 819 brs_sum = 0 820 # loop over helicity 821 for j in range(len(mvalue['jamp'])): 822 #values for the different lorentz boost 823 m_sum += abs(mvalue['jamp'][j][k])**2 824 brs_sum += abs(brsvalue['jamp'][j][k])**2 825 826 # Compare the different helicity 827 if not m_sum: 828 continue 829 ratio = brs_sum / m_sum 830 831 tmp_str = '\n' + fixed_string_length(' JAMP %s'%k , proc_col_size) + \ 832 fixed_string_length("%1.10e" % m_sum, col_size) + \ 833 fixed_string_length("%1.10e" % brs_sum, col_size) + \ 834 fixed_string_length("%1.10e" % ratio, col_size) 835 836 if ratio > 1e-15: 837 if not len(failed_proc_list) or failed_proc_list[-1] != proc: 838 fail_proc += 1 839 pass_proc -= 1 840 failed_proc_list.append(proc) 841 res_str += tmp_str + "Failed" 842 elif not proc_succeed: 843 res_str += tmp_str + "Passed" 844 845 846 res_str += "\nSummary: %i/%i passed, %i/%i failed" % \ 847 (pass_proc, pass_proc + fail_proc, 848 fail_proc, pass_proc + fail_proc) 849 850 if fail_proc != 0: 851 res_str += "\nFailed processes: %s" % ', '.join(failed_proc_list) 852 853 if output=='text': 854 return res_str 855 else: 856 return fail_proc

857 #=============================================================================== 858 # check_lorentz 859 #===============================================================================

860 -def check_lorentz(processes, param_card = None):

861 """ Check if the square matrix element (sum over helicity) is lorentz 862 invariant by boosting the momenta with different value.""" 863 864 if isinstance(processes, base_objects.ProcessDefinition): 865 # Generate a list of unique processes 866 # Extract IS and FS ids 867 multiprocess = processes 868 869 model = multiprocess.get('model') 870 871 # Initialize matrix element evaluation 872 evaluator = MatrixElementEvaluator(model, 873 auth_skipping = False, reuse = True) 874 875 # Set all widths to zero for lorentz check 876 for particle in evaluator.full_model.get('particles'): 877 if particle.get('width') != 'ZERO': 878 evaluator.full_model.get('parameter_dict')[\ 879 particle.get('width')] = 0. 880 return run_multiprocs_no_crossings(check_lorentz_process, 881 multiprocess, 882 evaluator) 883 elif isinstance(processes, base_objects.Process): 884 processes = base_objects.ProcessList([processes]) 885 elif isinstance(processes, base_objects.ProcessList): 886 pass 887 else: 888 raise InvalidCmd("processes is of non-supported format") 889 890 assert processes, "No processes given" 891 892 model = processes[0].get('model') 893 894 # Initialize matrix element evaluation 895 evaluator = MatrixElementEvaluator(model, param_card, 896 auth_skipping = False, reuse = True) 897 898 comparison_results = [] 899 900 # For each process, make sure we have set up leg numbers: 901 for process in processes: 902 # Check if we already checked process 903 #if check_already_checked([l.get('id') for l in process.get('legs') if \ 904 # not l.get('state')], 905 # [l.get('id') for l in process.get('legs') if \ 906 # l.get('state')], 907 # sorted_ids, process, model): 908 # continue 909 910 # Get process result 911 result = check_lorentz_process(process, evaluator) 912 if result: 913 comparison_results.append(result) 914 915 return comparison_results

916 917

918 -def check_lorentz_process(process, evaluator):

919 """Check gauge invariance for the process, unless it is already done.""" 920 921 amp_results = [] 922 model = process.get('model') 923 924 for i, leg in enumerate(process.get('legs')): 925 leg.set('number', i+1) 926 927 logger.info("Checking lorentz %s" % \ 928 process.nice_string().replace('Process', 'process')) 929 930 legs = process.get('legs') 931 # Generate a process with these legs 932 # Generate the amplitude for this process 933 try: 934 amplitude = diagram_generation.Amplitude(process) 935 except InvalidCmd: 936 logging.info("No diagrams for %s" % \ 937 process.nice_string().replace('Process', 'process')) 938 return None 939 940 if not amplitude.get('diagrams'): 941 # This process has no diagrams; go to next process 942 logging.info("No diagrams for %s" % \ 943 process.nice_string().replace('Process', 'process')) 944 return None 945 946 # Generate phase space point to use 947 p, w_rambo = evaluator.get_momenta(process) 948 949 # Generate the HelasMatrixElement for the process 950 matrix_element = helas_objects.HelasMatrixElement(amplitude, 951 gen_color = True) 952 953 data = evaluator.evaluate_matrix_element(matrix_element, p=p, output='jamp', 954 auth_skipping = True) 955 956 if data and data['m2']: 957 results = [data] 958 else: 959 return {'process':process.base_string(), 'results':'pass'} 960 961 for boost in range(1,4): 962 boost_p = boost_momenta(p, boost) 963 results.append(evaluator.evaluate_matrix_element(matrix_element, 964 p=boost_p, 965 output='jamp')) 966 967 968 return {'process': process.base_string(), 'results': results}

969 970

971 -def boost_momenta(p, boost_direction=1, beta=0.5):

972 """boost the set momenta in the 'boost direction' by the 'beta' 973 factor""" 974 boost_p = [] 975 gamma = 1/ math.sqrt(1 - beta**2) 976 for imp in p: 977 bosst_p = imp[boost_direction] 978 E, px, py, pz = imp 979 boost_imp = [] 980 # Energy: 981 boost_imp.append(gamma * E - gamma * beta * bosst_p) 982 # PX 983 if boost_direction == 1: 984 boost_imp.append(-gamma * beta * E + gamma * px) 985 else: 986 boost_imp.append(px) 987 # PY 988 if boost_direction == 2: 989 boost_imp.append(-gamma * beta * E + gamma * py) 990 else: 991 boost_imp.append(py) 992 # PZ 993 if boost_direction == 3: 994 boost_imp.append(-gamma * beta * E + gamma * pz) 995 else: 996 boost_imp.append(pz) 997 #Add the momenta to the list 998 boost_p.append(boost_imp) 999 1000 return boost_p

1001

1002 -def output_lorentz_inv(comparison_results, output='text'):

1003 """Present the results of a comparison in a nice list format 1004 if output='fail' return the number of failed process -- for test-- 1005 """ 1006 1007 proc_col_size = 17 1008 for proc, values in comparison_results: 1009 if len(proc) + 1 > proc_col_size: 1010 proc_col_size = len(proc) + 1 1011 1012 col_size = 17 1013 1014 pass_proc = 0 1015 fail_proc = 0 1016 no_check_proc = 0 1017 1018 failed_proc_list = [] 1019 no_check_proc_list = [] 1020 1021 res_str = fixed_string_length("Process", proc_col_size) + \ 1022 fixed_string_length("Min element", col_size) + \ 1023 fixed_string_length("Max element", col_size) + \ 1024 fixed_string_length("Relative diff.", col_size) + \ 1025 "Result" 1026 1027 for one_comp in comparison_results: 1028 proc = one_comp['process'] 1029 data = one_comp['results'] 1030 1031 if data == 'pass': 1032 no_check_proc += 1 1033 no_check_proc_list.append(proc) 1034 continue 1035 1036 values = [data[i]['m2'] for i in range(len(data))] 1037 1038 min_val = min(values) 1039 max_val = max(values) 1040 diff = (max_val - min_val) / max_val 1041 1042 res_str += '\n' + fixed_string_length(proc, proc_col_size) + \ 1043 fixed_string_length("%1.10e" % min_val, col_size) + \ 1044 fixed_string_length("%1.10e" % max_val, col_size) + \ 1045 fixed_string_length("%1.10e" % diff, col_size) 1046 1047 if diff < 1e-8: 1048 pass_proc += 1 1049 proc_succeed = True 1050 res_str += "Passed" 1051 else: 1052 fail_proc += 1 1053 proc_succeed = False 1054 failed_proc_list.append(proc) 1055 res_str += "Failed" 1056 1057 #check all the JAMP 1058 # loop over jamp 1059 for k in range(len(data[0]['jamp'][0])): 1060 sum = [0] * len(data) 1061 # loop over helicity 1062 for j in range(len(data[0]['jamp'])): 1063 #values for the different lorentz boost 1064 values = [abs(data[i]['jamp'][j][k])**2 for i in range(len(data))] 1065 sum = [sum[i] + values[i] for i in range(len(values))] 1066 1067 # Compare the different lorentz boost 1068 min_val = min(sum) 1069 max_val = max(sum) 1070 if not max_val: 1071 continue 1072 diff = (max_val - min_val) / max_val 1073 1074 tmp_str = '\n' + fixed_string_length(' JAMP %s'%k , proc_col_size) + \ 1075 fixed_string_length("%1.10e" % min_val, col_size) + \ 1076 fixed_string_length("%1.10e" % max_val, col_size) + \ 1077 fixed_string_length("%1.10e" % diff, col_size) 1078 1079 if diff > 1e-10: 1080 if not len(failed_proc_list) or failed_proc_list[-1] != proc: 1081 fail_proc += 1 1082 pass_proc -= 1 1083 failed_proc_list.append(proc) 1084 res_str += tmp_str + "Failed" 1085 elif not proc_succeed: 1086 res_str += tmp_str + "Passed" 1087 1088 1089 1090 res_str += "\nSummary: %i/%i passed, %i/%i failed" % \ 1091 (pass_proc, pass_proc + fail_proc, 1092 fail_proc, pass_proc + fail_proc) 1093 1094 if fail_proc != 0: 1095 res_str += "\nFailed processes: %s" % ', '.join(failed_proc_list) 1096 if no_check_proc: 1097 res_str += "\nNot checked processes: %s" % ', '.join(no_check_proc_list) 1098 1099 if output == 'text': 1100 return res_str 1101 else: 1102 return fail_proc

1103

Source Code for Module madgraph.various.process_checks