Package madgraph :: Package iolibs :: Module export_v4
[hide private]
[frames] | no frames]

Source Code for Module madgraph.iolibs.export_v4

   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  """Methods and classes to export matrix elements to v4 format.""" 
  16   
  17  import copy 
  18  import fractions 
  19  import glob 
  20  import logging 
  21  import math 
  22  import os 
  23  import re 
  24  import shutil 
  25  import subprocess 
  26  import sys 
  27   
  28  import madgraph.core.base_objects as base_objects 
  29  import madgraph.core.color_algebra as color 
  30  import madgraph.core.helas_objects as helas_objects 
  31  import madgraph.iolibs.drawing_eps as draw 
  32  import madgraph.iolibs.files as files 
  33  import madgraph.iolibs.group_subprocs as group_subprocs 
  34  import madgraph.iolibs.file_writers as writers 
  35  import madgraph.iolibs.gen_infohtml as gen_infohtml 
  36  import madgraph.iolibs.template_files as template_files 
  37  import madgraph.iolibs.ufo_expression_parsers as parsers 
  38  import madgraph.various.diagram_symmetry as diagram_symmetry 
  39  import madgraph.various.misc as misc 
  40  import madgraph.various.process_checks as process_checks 
  41   
  42   
  43   
  44  import aloha.create_aloha as create_aloha 
  45  import models.write_param_card as param_writer 
  46  import models.check_param_card as check_param_card 
  47   
  48  from madgraph import MadGraph5Error, MG5DIR 
  49  from madgraph.iolibs.files import cp, ln, mv 
  50  _file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/' 
  51  logger = logging.getLogger('madgraph.export_v4') 
  52   
  53  #=============================================================================== 
  54  # ProcessExporterFortran 
  55  #=============================================================================== 
56 -class ProcessExporterFortran(object):
57 """Class to take care of exporting a set of matrix elements to 58 Fortran (v4) format.""" 59
60 - def __init__(self, mgme_dir = "", dir_path = "", clean = False):
61 """Initiate the ProcessExporterFortran with directory information""" 62 self.mgme_dir = mgme_dir 63 self.dir_path = dir_path 64 self.clean = clean 65 self.model = None
66 67 #=========================================================================== 68 # copy the Template in a new directory. 69 #===========================================================================
70 - def copy_v4template(self, modelname):
71 """create the directory run_name as a copy of the MadEvent 72 Template, and clean the directory 73 """ 74 75 #First copy the full template tree if dir_path doesn't exit 76 if not os.path.isdir(self.dir_path): 77 assert self.mgme_dir, \ 78 "No valid MG_ME path given for MG4 run directory creation." 79 logger.info('initialize a new directory: %s' % \ 80 os.path.basename(self.dir_path)) 81 shutil.copytree(os.path.join(self.mgme_dir, 'Template'), 82 self.dir_path, True) 83 # Duplicate run_card and plot_card 84 for card in ['run_card', 'plot_card']: 85 try: 86 shutil.copy(os.path.join(self.dir_path, 'Cards', 87 card + '.dat'), 88 os.path.join(self.dir_path, 'Cards', 89 card + '_default.dat')) 90 except IOError: 91 info.warning("Failed to copy " + card + ".dat to default") 92 93 elif not os.path.isfile(os.path.join(self.dir_path, 'TemplateVersion.txt')): 94 assert self.mgme_dir, \ 95 "No valid MG_ME path given for MG4 run directory creation." 96 try: 97 shutil.copy(os.path.join(self.mgme_dir, 'MGMEVersion.txt'), self.dir_path) 98 except IOError: 99 MG5_version = misc.get_pkg_info() 100 open(os.path.join(self.dir_path, 'MGMEVersion.txt'), 'w').write( \ 101 "5." + MG5_version['version']) 102 103 #Ensure that the Template is clean 104 if self.clean: 105 logger.info('remove old information in %s' % \ 106 os.path.basename(self.dir_path)) 107 if os.environ.has_key('MADGRAPH_BASE'): 108 subprocess.call([os.path.join('bin', 'internal', 'clean_template'), 109 '--web'], cwd=self.dir_path) 110 else: 111 try: 112 subprocess.call([os.path.join('bin', 'internal', 'clean_template')], \ 113 cwd=self.dir_path) 114 except Exception, why: 115 raise MadGraph5Error('Failed to clean correctly %s: \n %s' \ 116 % (os.path.basename(self.dir_path),why)) 117 118 #Write version info 119 MG_version = misc.get_pkg_info() 120 open(os.path.join(self.dir_path, 'SubProcesses', 'MGVersion.txt'), 'w').write( 121 MG_version['version']) 122 123 124 # add the makefile in Source directory 125 filename = os.path.join(self.dir_path,'Source','makefile') 126 self.write_source_makefile(writers.FortranWriter(filename))
127 128 #=========================================================================== 129 # write a procdef_mg5 (an equivalent of the MG4 proc_card.dat) 130 #===========================================================================
131 - def write_procdef_mg5(self, file_pos, modelname, process_str):
132 """ write an equivalent of the MG4 proc_card in order that all the Madevent 133 Perl script of MadEvent4 are still working properly for pure MG5 run.""" 134 135 proc_card_template = template_files.mg4_proc_card.mg4_template 136 process_template = template_files.mg4_proc_card.process_template 137 process_text = '' 138 coupling = '' 139 new_process_content = [] 140 141 142 # First find the coupling and suppress the coupling from process_str 143 #But first ensure that coupling are define whithout spaces: 144 process_str = process_str.replace(' =', '=') 145 process_str = process_str.replace('= ', '=') 146 process_str = process_str.replace(',',' , ') 147 #now loop on the element and treat all the coupling 148 for info in process_str.split(): 149 if '=' in info: 150 coupling += info + '\n' 151 else: 152 new_process_content.append(info) 153 # Recombine the process_str (which is the input process_str without coupling 154 #info) 155 process_str = ' '.join(new_process_content) 156 157 #format the SubProcess 158 process_text += process_template.substitute({'process': process_str, \ 159 'coupling': coupling}) 160 161 text = proc_card_template.substitute({'process': process_text, 162 'model': modelname, 163 'multiparticle':''}) 164 ff = open(file_pos, 'w') 165 ff.write(text) 166 ff.close()
167 168 #=========================================================================== 169 # Create jpeg diagrams, html pages,proc_card_mg5.dat and madevent.tar.gz 170 #===========================================================================
171 - def finalize_v4_directory(self, matrix_elements, history = "", makejpg = False, 172 online = False, compiler='g77'):
173 """Function to finalize v4 directory, for inheritance. 174 """ 175 pass
176 177 #=========================================================================== 178 # write_matrix_element_v4 179 #===========================================================================
180 - def write_matrix_element_v4(self):
181 """Function to write a matrix.f file, for inheritance. 182 """ 183 pass
184 185 #=========================================================================== 186 # export the model 187 #===========================================================================
188 - def export_model_files(self, model_path):
189 """Configure the files/link of the process according to the model""" 190 191 # Import the model 192 for file in os.listdir(model_path): 193 if os.path.isfile(os.path.join(model_path, file)): 194 shutil.copy2(os.path.join(model_path, file), \ 195 os.path.join(self.dir_path, 'Source', 'MODEL'))
196 197 214 215 #=========================================================================== 216 # export the helas routine 217 #===========================================================================
218 - def export_helas(self, helas_path):
219 """Configure the files/link of the process according to the model""" 220 221 # Import helas routine 222 for filename in os.listdir(helas_path): 223 filepos = os.path.join(helas_path, filename) 224 if os.path.isfile(filepos): 225 if filepos.endswith('Makefile.template'): 226 cp(filepos, self.dir_path + '/Source/DHELAS/Makefile') 227 elif filepos.endswith('Makefile'): 228 pass 229 else: 230 cp(filepos, self.dir_path + '/Source/DHELAS')
231 # following lines do the same but whithout symbolic link 232 # 233 #def export_helas(mgme_dir, dir_path): 234 # 235 # # Copy the HELAS directory 236 # helas_dir = os.path.join(mgme_dir, 'HELAS') 237 # for filename in os.listdir(helas_dir): 238 # if os.path.isfile(os.path.join(helas_dir, filename)): 239 # shutil.copy2(os.path.join(helas_dir, filename), 240 # os.path.join(dir_path, 'Source', 'DHELAS')) 241 # shutil.move(os.path.join(dir_path, 'Source', 'DHELAS', 'Makefile.template'), 242 # os.path.join(dir_path, 'Source', 'DHELAS', 'Makefile')) 243 # 244 245 #=========================================================================== 246 # generate_subprocess_directory_v4 247 #===========================================================================
248 - def generate_subprocess_directory_v4(self, matrix_element, 249 fortran_model, 250 me_number):
251 """Routine to generate a subprocess directory (for inheritance)""" 252 253 pass
254 255 #=========================================================================== 256 # write_source_makefile 257 #===========================================================================
258 - def write_source_makefile(self, writer):
259 """Write the nexternal.inc file for MG4""" 260 261 262 path = os.path.join(_file_path,'iolibs','template_files','madevent_makefile_source') 263 set_of_lib = '$(LIBRARIES) $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext)' 264 text = open(path).read() % {'libraries': set_of_lib} 265 writer.write(text) 266 267 return True
268 269 #=========================================================================== 270 # write_nexternal_file 271 #===========================================================================
272 - def write_nexternal_file(self, writer, nexternal, ninitial):
273 """Write the nexternal.inc file for MG4""" 274 275 replace_dict = {} 276 277 replace_dict['nexternal'] = nexternal 278 replace_dict['ninitial'] = ninitial 279 280 file = """ \ 281 integer nexternal 282 parameter (nexternal=%(nexternal)d) 283 integer nincoming 284 parameter (nincoming=%(ninitial)d)""" % replace_dict 285 286 # Write the file 287 writer.writelines(file) 288 289 return True
290 291 #=========================================================================== 292 # write_pmass_file 293 #===========================================================================
294 - def write_pmass_file(self, writer, matrix_element):
295 """Write the pmass.inc file for MG4""" 296 297 model = matrix_element.get('processes')[0].get('model') 298 299 lines = [] 300 for wf in matrix_element.get_external_wavefunctions(): 301 mass = model.get('particle_dict')[wf.get('pdg_code')].get('mass') 302 if mass.lower() != "zero": 303 mass = "abs(%s)" % mass 304 305 lines.append("pmass(%d)=%s" % \ 306 (wf.get('number_external'), mass)) 307 308 # Write the file 309 writer.writelines(lines) 310 311 return True
312 313 #=========================================================================== 314 # write_ngraphs_file 315 #===========================================================================
316 - def write_ngraphs_file(self, writer, nconfigs):
317 """Write the ngraphs.inc file for MG4. Needs input from 318 write_configs_file.""" 319 320 file = " integer n_max_cg\n" 321 file = file + "parameter (n_max_cg=%d)" % nconfigs 322 323 # Write the file 324 writer.writelines(file) 325 326 return True
327 328 #=========================================================================== 329 # Routines to output UFO models in MG4 format 330 #=========================================================================== 331
332 - def convert_model_to_mg4(self, model, wanted_lorentz = [], 333 wanted_couplings = []):
334 """ Create a full valid MG4 model from a MG5 model (coming from UFO)""" 335 336 # create the MODEL 337 write_dir=os.path.join(self.dir_path, 'Source', 'MODEL') 338 model_builder = UFO_model_to_mg4(model, write_dir) 339 model_builder.build(wanted_couplings) 340 341 # Create and write ALOHA Routine 342 aloha_model = create_aloha.AbstractALOHAModel(model.get('name')) 343 if wanted_lorentz: 344 aloha_model.compute_subset(wanted_lorentz) 345 else: 346 aloha_model.compute_all(save=False) 347 write_dir=os.path.join(self.dir_path, 'Source', 'DHELAS') 348 aloha_model.write(write_dir, 'Fortran') 349 350 #copy Helas Template 351 cp(MG5DIR + '/aloha/template_files/Makefile_F', write_dir+'/makefile') 352 for filename in os.listdir(os.path.join(MG5DIR,'aloha','template_files')): 353 if not filename.lower().endswith('.f'): 354 continue 355 cp((MG5DIR + '/aloha/template_files/' + filename), write_dir) 356 create_aloha.write_aloha_file_inc(write_dir, '.f', '.o') 357 358 # Make final link in the Process 359 self.make_model_symbolic_link()
360 361 #=========================================================================== 362 # Helper functions 363 #===========================================================================
364 - def get_mg5_info_lines(self):
365 """Return info lines for MG5, suitable to place at beginning of 366 Fortran files""" 367 368 info = misc.get_pkg_info() 369 info_lines = "" 370 if info and info.has_key('version') and info.has_key('date'): 371 info_lines = "# Generated by MadGraph 5 v. %s, %s\n" % \ 372 (info['version'], info['date']) 373 info_lines = info_lines + \ 374 "# By the MadGraph Development Team\n" + \ 375 "# Please visit us at https://launchpad.net/madgraph5" 376 else: 377 info_lines = "# Generated by MadGraph 5\n" + \ 378 "# By the MadGraph Development Team\n" + \ 379 "# Please visit us at https://launchpad.net/madgraph5" 380 381 return info_lines
382
383 - def get_process_info_lines(self, matrix_element):
384 """Return info lines describing the processes for this matrix element""" 385 386 return"\n".join([ "C " + process.nice_string().replace('\n', '\nC * ') \ 387 for process in matrix_element.get('processes')])
388 389
390 - def get_helicity_lines(self, matrix_element):
391 """Return the Helicity matrix definition lines for this matrix element""" 392 393 helicity_line_list = [] 394 i = 0 395 for helicities in matrix_element.get_helicity_matrix(): 396 i = i + 1 397 int_list = [i, len(helicities)] 398 int_list.extend(helicities) 399 helicity_line_list.append(\ 400 ("DATA (NHEL(I,%4r),I=1,%d) /" + \ 401 ",".join(['%2r'] * len(helicities)) + "/") % tuple(int_list)) 402 403 return "\n".join(helicity_line_list)
404
405 - def get_ic_line(self, matrix_element):
406 """Return the IC definition line coming after helicities, required by 407 switchmom in madevent""" 408 409 nexternal = matrix_element.get_nexternal_ninitial()[0] 410 int_list = range(1, nexternal + 1) 411 412 return "DATA (IC(I,1),I=1,%i) /%s/" % (nexternal, 413 ",".join([str(i) for \ 414 i in int_list]))
415
416 - def get_color_data_lines(self, matrix_element, n=6):
417 """Return the color matrix definition lines for this matrix element. Split 418 rows in chunks of size n.""" 419 420 if not matrix_element.get('color_matrix'): 421 return ["DATA Denom(1)/1/", "DATA (CF(i,1),i=1,1) /1/"] 422 else: 423 ret_list = [] 424 my_cs = color.ColorString() 425 for index, denominator in \ 426 enumerate(matrix_element.get('color_matrix').\ 427 get_line_denominators()): 428 # First write the common denominator for this color matrix line 429 ret_list.append("DATA Denom(%i)/%i/" % (index + 1, denominator)) 430 # Then write the numerators for the matrix elements 431 num_list = matrix_element.get('color_matrix').\ 432 get_line_numerators(index, denominator) 433 434 for k in xrange(0, len(num_list), n): 435 ret_list.append("DATA (CF(i,%3r),i=%3r,%3r) /%s/" % \ 436 (index + 1, k + 1, min(k + n, len(num_list)), 437 ','.join(["%5r" % i for i in num_list[k:k + n]]))) 438 my_cs.from_immutable(sorted(matrix_element.get('color_basis').keys())[index]) 439 ret_list.append("C %s" % repr(my_cs)) 440 return ret_list
441 442
443 - def get_den_factor_line(self, matrix_element):
444 """Return the denominator factor line for this matrix element""" 445 446 return "DATA IDEN/%2r/" % \ 447 matrix_element.get_denominator_factor()
448
449 - def get_icolamp_lines(self, mapconfigs, matrix_element, num_matrix_element):
450 """Return the ICOLAMP matrix, showing which JAMPs contribute to 451 which configs (diagrams).""" 452 453 ret_list = [] 454 455 booldict = {False: ".false.", True: ".true."} 456 457 if not matrix_element.get('color_basis'): 458 # No color, so only one color factor. Simply write a ".true." 459 # for each config (i.e., each diagram with only 3 particle 460 # vertices 461 configs = len(mapconfigs) 462 ret_list.append("DATA(icolamp(1,i,%d),i=1,%d)/%s/" % \ 463 (num_matrix_element, configs, 464 ','.join([".true." for i in range(configs)]))) 465 return ret_list 466 467 # There is a color basis - create a list showing which JAMPs have 468 # contributions to which configs 469 470 # Only want to include leading color flows, so find max_Nc 471 color_basis = matrix_element.get('color_basis') 472 max_Nc = max(sum([[v[4] for v in val] for val in color_basis.values()], 473 [])) 474 475 # Crate dictionary between diagram number and JAMP number 476 diag_jamp = {} 477 for ijamp, col_basis_elem in \ 478 enumerate(sorted(matrix_element.get('color_basis').keys())): 479 for diag_tuple in matrix_element.get('color_basis')[col_basis_elem]: 480 # Only use color flows with Nc == max_Nc 481 if diag_tuple[4] == max_Nc: 482 diag_num = diag_tuple[0] + 1 483 # Add this JAMP number to this diag_num 484 diag_jamp[diag_num] = diag_jamp.setdefault(diag_num, []) + \ 485 [ijamp+1] 486 487 colamps = ijamp + 1 488 489 for iconfig, num_diag in enumerate(mapconfigs): 490 if num_diag == 0: 491 continue 492 493 # List of True or False 494 bool_list = [(i + 1 in diag_jamp[num_diag]) for i in \ 495 range(colamps)] 496 # Add line 497 ret_list.append("DATA(icolamp(i,%d,%d),i=1,%d)/%s/" % \ 498 (iconfig+1, num_matrix_element, colamps, 499 ','.join(["%s" % booldict[b] for b in \ 500 bool_list]))) 501 502 return ret_list
503
504 - def get_amp2_lines(self, matrix_element, config_map = []):
505 """Return the amp2(i) = sum(amp for diag(i))^2 lines""" 506 507 nexternal, ninitial = matrix_element.get_nexternal_ninitial() 508 # Get minimum legs in a vertex 509 minvert = min([max(diag.get_vertex_leg_numbers()) for diag in \ 510 matrix_element.get('diagrams')]) 511 512 ret_lines = [] 513 if config_map: 514 # In this case, we need to sum up all amplitudes that have 515 # identical topologies, as given by the config_map (which 516 # gives the topology/config for each of the diagrams 517 diagrams = matrix_element.get('diagrams') 518 # Combine the diagrams with identical topologies 519 config_to_diag_dict = {} 520 for idiag, diag in enumerate(matrix_element.get('diagrams')): 521 if config_map[idiag] == 0: 522 continue 523 try: 524 config_to_diag_dict[config_map[idiag]].append(idiag) 525 except KeyError: 526 config_to_diag_dict[config_map[idiag]] = [idiag] 527 # Write out the AMP2s summing squares of amplitudes belonging 528 # to eiher the same diagram or different diagrams with 529 # identical propagator properties. Note that we need to use 530 # AMP2 number corresponding to the first diagram number used 531 # for that AMP2. 532 for config in sorted(config_to_diag_dict.keys()): 533 534 line = "AMP2(%(num)d)=AMP2(%(num)d)+" % \ 535 {"num": (config_to_diag_dict[config][0] + 1)} 536 537 line += "+".join(["AMP(%(num)d)*dconjg(AMP(%(num)d))" % \ 538 {"num": a.get('number')} for a in \ 539 sum([diagrams[idiag].get('amplitudes') for \ 540 idiag in config_to_diag_dict[config]], [])]) 541 ret_lines.append(line) 542 else: 543 for idiag, diag in enumerate(matrix_element.get('diagrams')): 544 # Ignore any diagrams with 4-particle vertices. 545 if max(diag.get_vertex_leg_numbers()) > minvert: 546 continue 547 # Now write out the expression for AMP2, meaning the sum of 548 # squared amplitudes belonging to the same diagram 549 line = "AMP2(%(num)d)=AMP2(%(num)d)+" % {"num": (idiag + 1)} 550 line += "+".join(["AMP(%(num)d)*dconjg(AMP(%(num)d))" % \ 551 {"num": a.get('number')} for a in \ 552 diag.get('amplitudes')]) 553 ret_lines.append(line) 554 555 return ret_lines
556
557 - def get_JAMP_lines(self, matrix_element):
558 """Return the JAMP = sum(fermionfactor * AMP(i)) lines""" 559 560 res_list = [] 561 562 for i, coeff_list in \ 563 enumerate(matrix_element.get_color_amplitudes()): 564 565 res = "JAMP(%i)=" % (i + 1) 566 567 # Optimization: if all contributions to that color basis element have 568 # the same coefficient (up to a sign), put it in front 569 list_fracs = [abs(coefficient[0][1]) for coefficient in coeff_list] 570 common_factor = False 571 diff_fracs = list(set(list_fracs)) 572 if len(diff_fracs) == 1 and abs(diff_fracs[0]) != 1: 573 common_factor = True 574 global_factor = diff_fracs[0] 575 res = res + '%s(' % self.coeff(1, global_factor, False, 0) 576 577 for (coefficient, amp_number) in coeff_list: 578 if common_factor: 579 res = res + "%sAMP(%d)" % (self.coeff(coefficient[0], 580 coefficient[1] / abs(coefficient[1]), 581 coefficient[2], 582 coefficient[3]), 583 amp_number) 584 else: 585 res = res + "%sAMP(%d)" % (self.coeff(coefficient[0], 586 coefficient[1], 587 coefficient[2], 588 coefficient[3]), 589 amp_number) 590 591 if common_factor: 592 res = res + ')' 593 594 res_list.append(res) 595 596 return res_list
597
598 - def get_pdf_lines(self, matrix_element, ninitial, subproc_group = False):
599 """Generate the PDF lines for the auto_dsig.f file""" 600 601 processes = matrix_element.get('processes') 602 model = processes[0].get('model') 603 604 pdf_definition_lines = "" 605 pdf_data_lines = "" 606 pdf_lines = "" 607 608 if ninitial == 1: 609 pdf_lines = "PD(0) = 0d0\nIPROC = 0\n" 610 for i, proc in enumerate(processes): 611 process_line = proc.base_string() 612 pdf_lines = pdf_lines + "IPROC=IPROC+1 ! " + process_line 613 pdf_lines = pdf_lines + "\nPD(IPROC)=PD(IPROC-1) + 1d0\n" 614 else: 615 # Pick out all initial state particles for the two beams 616 initial_states = [sorted(list(set([p.get_initial_pdg(1) for \ 617 p in processes]))), 618 sorted(list(set([p.get_initial_pdg(2) for \ 619 p in processes])))] 620 621 # Prepare all variable names 622 pdf_codes = dict([(p, model.get_particle(p).get_name()) for p in \ 623 sum(initial_states,[])]) 624 for key,val in pdf_codes.items(): 625 pdf_codes[key] = val.replace('~','x').replace('+','p').replace('-','m') 626 627 # Set conversion from PDG code to number used in PDF calls 628 pdgtopdf = {21: 0, 22: 7} 629 630 # Fill in missing entries of pdgtopdf 631 for pdg in sum(initial_states,[]): 632 if not pdg in pdgtopdf and not pdg in pdgtopdf.values(): 633 pdgtopdf[pdg] = pdg 634 elif pdg not in pdgtopdf and pdg in pdgtopdf.values(): 635 # If any particle has pdg code 7, we need to use something else 636 pdgtopdf[pdg] = 6000000 + pdg 637 638 # Get PDF variable declarations for all initial states 639 for i in [0,1]: 640 pdf_definition_lines += "DOUBLE PRECISION " + \ 641 ",".join(["%s%d" % (pdf_codes[pdg],i+1) \ 642 for pdg in \ 643 initial_states[i]]) + \ 644 "\n" 645 646 # Get PDF data lines for all initial states 647 for i in [0,1]: 648 pdf_data_lines += "DATA " + \ 649 ",".join(["%s%d" % (pdf_codes[pdg],i+1) \ 650 for pdg in initial_states[i]]) + \ 651 "/%d*1D0/" % len(initial_states[i]) + \ 652 "\n" 653 654 # Get PDF lines for all different initial states 655 for i, init_states in enumerate(initial_states): 656 if subproc_group: 657 pdf_lines = pdf_lines + \ 658 "IF (ABS(LPP(IB(%d))).GE.1) THEN\nLP=SIGN(1,LPP(IB(%d)))\n" \ 659 % (i + 1, i + 1) 660 else: 661 pdf_lines = pdf_lines + \ 662 "IF (ABS(LPP(%d)) .GE. 1) THEN\nLP=SIGN(1,LPP(%d))\n" \ 663 % (i + 1, i + 1) 664 665 for initial_state in init_states: 666 if initial_state in pdf_codes.keys(): 667 if subproc_group: 668 pdf_lines = pdf_lines + \ 669 ("%s%d=PDG2PDF(ABS(LPP(IB(%d))),%d*LP," + \ 670 "XBK(IB(%d)),DSQRT(Q2FACT(%d)))\n") % \ 671 (pdf_codes[initial_state], 672 i + 1, i + 1, pdgtopdf[initial_state], 673 i + 1, i + 1) 674 else: 675 pdf_lines = pdf_lines + \ 676 ("%s%d=PDG2PDF(ABS(LPP(%d)),%d*LP," + \ 677 "XBK(%d),DSQRT(Q2FACT(%d)))\n") % \ 678 (pdf_codes[initial_state], 679 i + 1, i + 1, pdgtopdf[initial_state], 680 i + 1, i + 1) 681 pdf_lines = pdf_lines + "ENDIF\n" 682 683 # Add up PDFs for the different initial state particles 684 pdf_lines = pdf_lines + "PD(0) = 0d0\nIPROC = 0\n" 685 for proc in processes: 686 process_line = proc.base_string() 687 pdf_lines = pdf_lines + "IPROC=IPROC+1 ! " + process_line 688 pdf_lines = pdf_lines + "\nPD(IPROC)=PD(IPROC-1) + " 689 for ibeam in [1, 2]: 690 initial_state = proc.get_initial_pdg(ibeam) 691 if initial_state in pdf_codes.keys(): 692 pdf_lines = pdf_lines + "%s%d*" % \ 693 (pdf_codes[initial_state], ibeam) 694 else: 695 pdf_lines = pdf_lines + "1d0*" 696 # Remove last "*" from pdf_lines 697 pdf_lines = pdf_lines[:-1] + "\n" 698 699 # Remove last line break from the return variables 700 return pdf_definition_lines[:-1], pdf_data_lines[:-1], pdf_lines[:-1]
701 702 703 #=========================================================================== 704 # Global helper methods 705 #=========================================================================== 706
707 - def coeff(self, ff_number, frac, is_imaginary, Nc_power, Nc_value=3):
708 """Returns a nicely formatted string for the coefficients in JAMP lines""" 709 710 total_coeff = ff_number * frac * fractions.Fraction(Nc_value) ** Nc_power 711 712 if total_coeff == 1: 713 if is_imaginary: 714 return '+imag1*' 715 else: 716 return '+' 717 elif total_coeff == -1: 718 if is_imaginary: 719 return '-imag1*' 720 else: 721 return '-' 722 723 res_str = '%+i' % total_coeff.numerator 724 725 if total_coeff.denominator != 1: 726 # Check if total_coeff is an integer 727 res_str = res_str + './%i.' % total_coeff.denominator 728 729 if is_imaginary: 730 res_str = res_str + '*imag1' 731 732 return res_str + '*'
733
734 - def set_compiler(self, default_compiler):
735 """Set compiler based on what's available on the system""" 736 737 738 # Check for compiler 739 if misc.which(default_compiler): 740 compiler = default_compiler 741 elif misc.which('gfortran'): 742 compiler = 'gfortran' 743 elif misc.which('g77'): 744 compiler = 'g77' 745 elif misc.which('f77'): 746 compiler = 'f77' 747 else: 748 logger.warning('No Fortran Compiler detected! Please install one') 749 compiler = default_compiler 750 logger.info('Use Fortran compiler ' + compiler) 751 self.replace_make_opt_compiler(compiler) 752 # Replace also for Template but not for cluster 753 if not os.environ.has_key('MADGRAPH_DATA'): 754 self.replace_make_opt_compiler(compiler, os.path.join(MG5DIR, 'Template'))
755
756 - def replace_make_opt_compiler(self, compiler, root_dir = ""):
757 """Set FC=compiler in Source/make_opts""" 758 759 if not root_dir: 760 root_dir = self.dir_path 761 make_opts = os.path.join(root_dir, 'Source', 'make_opts') 762 lines = open(make_opts).read().split('\n') 763 FC_re = re.compile('^(\s*)FC\s*=\s*.+\s*$') 764 for iline, line in enumerate(lines): 765 FC_result = FC_re.match(line) 766 if FC_result: 767 lines[iline] = FC_result.group(1) + "FC=" + compiler 768 outfile = open(make_opts, 'w') 769 outfile.write('\n'.join(lines))
770 771 #=============================================================================== 772 # ProcessExporterFortranSA 773 #===============================================================================
774 -class ProcessExporterFortranSA(ProcessExporterFortran):
775 """Class to take care of exporting a set of matrix elements to 776 MadGraph v4 StandAlone format.""" 777
778 - def copy_v4template(self, modelname):
779 """Additional actions needed for setup of Template 780 """ 781 782 783 #First copy the full template tree if dir_path doesn't exit 784 if os.path.isdir(self.dir_path): 785 return 786 787 logger.info('initialize a new standalone directory: %s' % \ 788 os.path.basename(self.dir_path)) 789 temp_dir = os.path.join(self.mgme_dir, 'Template') 790 791 792 # Create the directory structure 793 os.mkdir(self.dir_path) 794 os.mkdir(os.path.join(self.dir_path, 'Source')) 795 os.mkdir(os.path.join(self.dir_path, 'Source', 'MODEL')) 796 os.mkdir(os.path.join(self.dir_path, 'Source', 'DHELAS')) 797 os.mkdir(os.path.join(self.dir_path, 'SubProcesses')) 798 os.mkdir(os.path.join(self.dir_path, 'bin')) 799 os.mkdir(os.path.join(self.dir_path, 'bin', 'internal')) 800 os.mkdir(os.path.join(self.dir_path, 'lib')) 801 os.mkdir(os.path.join(self.dir_path, 'Cards')) 802 803 # Information at top-level 804 #Write version info 805 shutil.copy(os.path.join(temp_dir, 'TemplateVersion.txt'), self.dir_path) 806 try: 807 shutil.copy(os.path.join(self.mgme_dir, 'MGMEVersion.txt'), self.dir_path) 808 except IOError: 809 MG5_version = misc.get_pkg_info() 810 open(os.path.join(self.dir_path, 'MGMEVersion.txt'), 'w').write( \ 811 "5." + MG5_version['version']) 812 813 # Add file in bin directory 814 shutil.copy(os.path.join(temp_dir, 'bin', 'change_compiler.py'), 815 os.path.join(self.dir_path, 'bin')) 816 817 # Add file in SubProcesses 818 shutil.copy(os.path.join(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'makefile_sa_f_sp'), 819 os.path.join(self.dir_path, 'SubProcesses', 'makefile')) 820 shutil.copy(os.path.join(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'check_sa.f'), 821 os.path.join(self.dir_path, 'SubProcesses', 'check_sa.f')) 822 823 # Add file in Source 824 shutil.copy(os.path.join(temp_dir, 'Source', 'make_opts'), 825 os.path.join(self.dir_path, 'Source')) 826 # add the makefile 827 filename = os.path.join(self.dir_path,'Source','makefile') 828 self.write_source_makefile(writers.FortranWriter(filename))
829 830 #=========================================================================== 831 # export model files 832 #===========================================================================
833 - def export_model_files(self, model_path):
834 """export the model dependent files for V4 model""" 835 836 super(ProcessExporterFortranSA,self).export_model_files(model_path) 837 # Add the routine update_as_param in v4 model 838 # This is a function created in the UFO 839 840 841 text = open(os.path.join(self.dir_path,'SubProcesses','check_sa.f')).read() 842 text = text.replace('call setpara(\'param_card.dat\')', 'call setpara(\'param_card.dat\', .true.)') 843 fsock = open(os.path.join(self.dir_path,'SubProcesses','check_sa.f'), 'w') 844 fsock.write(text) 845 fsock.close() 846 847 self.make_model_symbolic_link()
848 849 #=========================================================================== 850 # Make the Helas and Model directories for Standalone directory 851 #===========================================================================
852 - def make(self):
853 """Run make in the DHELAS and MODEL directories, to set up 854 everything for running standalone 855 """ 856 857 source_dir = os.path.join(self.dir_path, "Source") 858 logger.info("Running make for Helas") 859 misc.compile(arg=['../lib/libdhelas.a'], cwd=source_dir, mode='fortran') 860 logger.info("Running make for Model") 861 misc.compile(arg=['../lib/libmodel.a'], cwd=source_dir, mode='fortran')
862 863 #=========================================================================== 864 # Create proc_card_mg5.dat for Standalone directory 865 #===========================================================================
866 - def finalize_v4_directory(self, matrix_elements, history, makejpg = False, 867 online = False, compiler='g77'):
868 """Finalize Standalone MG4 directory by generation proc_card_mg5.dat""" 869 870 self.set_compiler(compiler) 871 self.make() 872 873 # Write command history as proc_card_mg5 874 if os.path.isdir(os.path.join(self.dir_path, 'Cards')): 875 output_file = os.path.join(self.dir_path, 'Cards', 'proc_card_mg5.dat') 876 output_file = open(output_file, 'w') 877 text = ('\n'.join(history) + '\n') % misc.get_time_info() 878 output_file.write(text) 879 output_file.close()
880 881 #=========================================================================== 882 # generate_subprocess_directory_v4 883 #===========================================================================
884 - def generate_subprocess_directory_v4(self, matrix_element, 885 fortran_model):
886 """Generate the Pxxxxx directory for a subprocess in MG4 standalone, 887 including the necessary matrix.f and nexternal.inc files""" 888 889 cwd = os.getcwd() 890 891 # Create the directory PN_xx_xxxxx in the specified path 892 dirpath = os.path.join(self.dir_path, 'SubProcesses', \ 893 "P%s" % matrix_element.get('processes')[0].shell_string()) 894 895 try: 896 os.mkdir(dirpath) 897 except os.error as error: 898 logger.warning(error.strerror + " " + dirpath) 899 900 try: 901 os.chdir(dirpath) 902 except os.error: 903 logger.error('Could not cd to directory %s' % dirpath) 904 return 0 905 906 logger.info('Creating files in directory %s' % dirpath) 907 908 # Extract number of external particles 909 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 910 911 # Create the matrix.f file and the nexternal.inc file 912 filename = 'matrix.f' 913 calls = self.write_matrix_element_v4( 914 writers.FortranWriter(filename), 915 matrix_element, 916 fortran_model) 917 918 filename = 'nexternal.inc' 919 self.write_nexternal_file(writers.FortranWriter(filename), 920 nexternal, ninitial) 921 922 filename = 'pmass.inc' 923 self.write_pmass_file(writers.FortranWriter(filename), 924 matrix_element) 925 926 filename = 'ngraphs.inc' 927 self.write_ngraphs_file(writers.FortranWriter(filename), 928 len(matrix_element.get_all_amplitudes())) 929 930 # Generate diagrams 931 filename = "matrix.ps" 932 plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\ 933 get('diagrams'), 934 filename, 935 model=matrix_element.get('processes')[0].\ 936 get('model'), 937 amplitude=True) 938 logger.info("Generating Feynman diagrams for " + \ 939 matrix_element.get('processes')[0].nice_string()) 940 plot.draw() 941 942 linkfiles = ['check_sa.f', 'coupl.inc', 'makefile'] 943 944 945 for file in linkfiles: 946 ln('../%s' % file) 947 948 # Return to original PWD 949 os.chdir(cwd) 950 951 if not calls: 952 calls = 0 953 return calls
954 955 956 #=========================================================================== 957 # write_source_makefile 958 #===========================================================================
959 - def write_source_makefile(self, writer):
960 """Write the nexternal.inc file for MG4""" 961 962 963 path = os.path.join(_file_path,'iolibs','template_files','madevent_makefile_source') 964 set_of_lib = '$(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext)' 965 text = open(path).read() % {'libraries': set_of_lib} 966 writer.write(text) 967 968 return True
969 970 971 #=========================================================================== 972 # write_matrix_element_v4 973 #===========================================================================
974 - def write_matrix_element_v4(self, writer, matrix_element, fortran_model):
975 """Export a matrix element to a matrix.f file in MG4 standalone format""" 976 977 if not matrix_element.get('processes') or \ 978 not matrix_element.get('diagrams'): 979 return 0 980 981 if not isinstance(writer, writers.FortranWriter): 982 raise writers.FortranWriter.FortranWriterError(\ 983 "writer not FortranWriter") 984 985 # Set lowercase/uppercase Fortran code 986 writers.FortranWriter.downcase = False 987 988 replace_dict = {} 989 990 # Extract version number and date from VERSION file 991 info_lines = self.get_mg5_info_lines() 992 replace_dict['info_lines'] = info_lines 993 994 # Extract process info lines 995 process_lines = self.get_process_info_lines(matrix_element) 996 replace_dict['process_lines'] = process_lines 997 998 # Extract number of external particles 999 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 1000 replace_dict['nexternal'] = nexternal 1001 1002 # Extract ncomb 1003 ncomb = matrix_element.get_helicity_combinations() 1004 replace_dict['ncomb'] = ncomb 1005 1006 # Extract helicity lines 1007 helicity_lines = self.get_helicity_lines(matrix_element) 1008 replace_dict['helicity_lines'] = helicity_lines 1009 1010 # Extract overall denominator 1011 # Averaging initial state color, spin, and identical FS particles 1012 den_factor_line = self.get_den_factor_line(matrix_element) 1013 replace_dict['den_factor_line'] = den_factor_line 1014 1015 # Extract ngraphs 1016 ngraphs = matrix_element.get_number_of_amplitudes() 1017 replace_dict['ngraphs'] = ngraphs 1018 1019 # Extract nwavefuncs 1020 nwavefuncs = matrix_element.get_number_of_wavefunctions() 1021 replace_dict['nwavefuncs'] = nwavefuncs 1022 1023 # Extract ncolor 1024 ncolor = max(1, len(matrix_element.get('color_basis'))) 1025 replace_dict['ncolor'] = ncolor 1026 1027 # Extract color data lines 1028 color_data_lines = self.get_color_data_lines(matrix_element) 1029 replace_dict['color_data_lines'] = "\n".join(color_data_lines) 1030 1031 # Extract helas calls 1032 helas_calls = fortran_model.get_matrix_element_calls(\ 1033 matrix_element) 1034 replace_dict['helas_calls'] = "\n".join(helas_calls) 1035 1036 # Extract JAMP lines 1037 jamp_lines = self.get_JAMP_lines(matrix_element) 1038 replace_dict['jamp_lines'] = '\n'.join(jamp_lines) 1039 1040 file = open(os.path.join(_file_path, \ 1041 'iolibs/template_files/matrix_standalone_v4.inc')).read() 1042 file = file % replace_dict 1043 1044 # Write the file 1045 writer.writelines(file) 1046 1047 return len(filter(lambda call: call.find('#') != 0, helas_calls))
1048 1049 #=============================================================================== 1050 # ProcessExporterFortranME 1051 #===============================================================================
1052 -class ProcessExporterFortranME(ProcessExporterFortran):
1053 """Class to take care of exporting a set of matrix elements to 1054 MadEvent format.""" 1055 1056 matrix_file = "matrix_madevent_v4.inc" 1057
1058 - def copy_v4template(self, modelname):
1059 """Additional actions needed for setup of Template 1060 """ 1061 1062 super(ProcessExporterFortranME, self).copy_v4template(modelname) 1063 1064 # File created from Template (Different in some child class) 1065 filename = os.path.join(self.dir_path,'Source','run_config.inc') 1066 self.write_run_config_file(writers.FortranWriter(filename)) 1067 1068 # The next file are model dependant (due to SLAH convention) 1069 self.model_name = modelname 1070 # Add the combine_events.f 1071 filename = os.path.join(self.dir_path,'Source','combine_events.f') 1072 self.write_combine_events(writers.FortranWriter(filename)) 1073 # Add the symmetry.f 1074 filename = os.path.join(self.dir_path,'SubProcesses','symmetry.f') 1075 self.write_symmetry(writers.FortranWriter(filename)) 1076 # Add the driver.f 1077 filename = os.path.join(self.dir_path,'SubProcesses','driver.f') 1078 self.write_driver(writers.FortranWriter(filename)) 1079 # Copy the different python file in the Template 1080 self.copy_python_file()
1081 1082 1083 1084 1085 #=========================================================================== 1086 # generate_subprocess_directory_v4 1087 #===========================================================================
1088 - def copy_python_file(self):
1089 """copy the python file require for the Template""" 1090 1091 cp(_file_path+'/interface/madevent_interface.py', 1092 self.dir_path+'/bin/internal/madevent_interface.py') 1093 cp(_file_path+'/interface/extended_cmd.py', 1094 self.dir_path+'/bin/internal/extended_cmd.py') 1095 cp(_file_path+'/various/misc.py', self.dir_path+'/bin/internal/misc.py') 1096 cp(_file_path+'/iolibs/files.py', self.dir_path+'/bin/internal/files.py') 1097 cp(_file_path+'/iolibs/save_load_object.py', 1098 self.dir_path+'/bin/internal/save_load_object.py') 1099 cp(_file_path+'../models/check_param_card.py', 1100 self.dir_path+'/bin/internal/check_param_card.py') 1101 cp(_file_path+'/__init__.py', self.dir_path+'/bin/internal/__init__.py') 1102 cp(_file_path+'/various/gen_crossxhtml.py', 1103 self.dir_path+'/bin/internal/gen_crossxhtml.py') 1104 cp(_file_path+'/various/banner.py', 1105 self.dir_path+'/bin/internal/banner.py') 1106 cp(_file_path+'/various/cluster.py', 1107 self.dir_path+'/bin/internal/cluster.py') 1108 cp(_file_path+'/various/sum_html.py', 1109 self.dir_path+'/bin/internal/sum_html.py') 1110 cp(_file_path+'/interface/.mg5_logging.conf', 1111 self.dir_path+'/bin/internal/me5_logging.conf') 1112 cp(_file_path+'/interface/coloring_logging.py', 1113 self.dir_path+'/bin/internal/coloring_logging.py')
1114 1115 #=========================================================================== 1116 # export model files 1117 #===========================================================================
1118 - def export_model_files(self, model_path):
1119 """export the model dependent files""" 1120 1121 super(ProcessExporterFortranME,self).export_model_files(model_path) 1122 # Add the routine update_as_param in v4 model 1123 # This is a function created in the UFO 1124 text=""" 1125 subroutine update_as_param() 1126 call setpara('param_card.dat',.false.) 1127 return 1128 end 1129 """ 1130 ff = open(os.path.join(self.dir_path, 'Source', 'MODEL', 'couplings.f'),'a') 1131 ff.write(text) 1132 ff.close() 1133 1134 # Add the symmetry.f 1135 filename = os.path.join(self.dir_path,'SubProcesses','symmetry.f') 1136 self.write_symmetry(writers.FortranWriter(filename), v5=False) 1137 1138 # Add the driver.f 1139 filename = os.path.join(self.dir_path,'SubProcesses','driver.f') 1140 self.write_driver(writers.FortranWriter(filename), v5=False) 1141 1142 # Modify setrun.f 1143 text = open(os.path.join(self.dir_path,'Source','setrun.f')).read() 1144 text = text.replace('call setpara(param_card_name)', 'call setpara(param_card_name, .true.)') 1145 fsock = open(os.path.join(self.dir_path,'Source','setrun.f'), 'w') 1146 fsock.write(text) 1147 fsock.close() 1148 1149 self.make_model_symbolic_link()
1150 1151 1152 #=========================================================================== 1153 # generate_subprocess_directory_v4 1154 #===========================================================================
1155 - def generate_subprocess_directory_v4(self, matrix_element, 1156 fortran_model, 1157 me_number):
1158 """Generate the Pxxxxx directory for a subprocess in MG4 madevent, 1159 including the necessary matrix.f and various helper files""" 1160 1161 cwd = os.getcwd() 1162 path = os.path.join(self.dir_path, 'SubProcesses') 1163 1164 1165 if not self.model: 1166 self.model = matrix_element.get('processes')[0].get('model') 1167 1168 1169 1170 os.chdir(path) 1171 # Create the directory PN_xx_xxxxx in the specified path 1172 subprocdir = "P%s" % matrix_element.get('processes')[0].shell_string() 1173 try: 1174 os.mkdir(subprocdir) 1175 except os.error as error: 1176 logger.warning(error.strerror + " " + subprocdir) 1177 1178 try: 1179 os.chdir(subprocdir) 1180 except os.error: 1181 logger.error('Could not cd to directory %s' % subprocdir) 1182 return 0 1183 1184 logger.info('Creating files in directory %s' % subprocdir) 1185 1186 # Extract number of external particles 1187 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 1188 1189 # Create the matrix.f file, auto_dsig.f file and all inc files 1190 filename = 'matrix.f' 1191 calls, ncolor = \ 1192 self.write_matrix_element_v4(writers.FortranWriter(filename), 1193 matrix_element, 1194 fortran_model) 1195 1196 filename = 'auto_dsig.f' 1197 self.write_auto_dsig_file(writers.FortranWriter(filename), 1198 matrix_element) 1199 1200 filename = 'configs.inc' 1201 mapconfigs, s_and_t_channels = self.write_configs_file(\ 1202 writers.FortranWriter(filename), 1203 matrix_element) 1204 1205 filename = 'config_subproc_map.inc' 1206 self.write_config_subproc_map_file(writers.FortranWriter(filename), 1207 s_and_t_channels) 1208 1209 filename = 'coloramps.inc' 1210 self.write_coloramps_file(writers.FortranWriter(filename), 1211 mapconfigs, 1212 matrix_element) 1213 1214 filename = 'get_color.f' 1215 self.write_colors_file(writers.FortranWriter(filename), 1216 matrix_element) 1217 1218 filename = 'decayBW.inc' 1219 self.write_decayBW_file(writers.FortranWriter(filename), 1220 s_and_t_channels) 1221 1222 filename = 'dname.mg' 1223 self.write_dname_file(writers.FortranWriter(filename), 1224 "P"+matrix_element.get('processes')[0].shell_string()) 1225 1226 filename = 'iproc.dat' 1227 self.write_iproc_file(writers.FortranWriter(filename), 1228 me_number) 1229 1230 filename = 'leshouche.inc' 1231 self.write_leshouche_file(writers.FortranWriter(filename), 1232 matrix_element) 1233 1234 filename = 'maxamps.inc' 1235 self.write_maxamps_file(writers.FortranWriter(filename), 1236 len(matrix_element.get('diagrams')), 1237 ncolor, 1238 len(matrix_element.get('processes')), 1239 1) 1240 1241 filename = 'mg.sym' 1242 self.write_mg_sym_file(writers.FortranWriter(filename), 1243 matrix_element) 1244 1245 filename = 'ncombs.inc' 1246 self.write_ncombs_file(writers.FortranWriter(filename), 1247 nexternal) 1248 1249 filename = 'nexternal.inc' 1250 self.write_nexternal_file(writers.FortranWriter(filename), 1251 nexternal, ninitial) 1252 1253 filename = 'ngraphs.inc' 1254 self.write_ngraphs_file(writers.FortranWriter(filename), 1255 len(mapconfigs)) 1256 1257 1258 filename = 'pmass.inc' 1259 self.write_pmass_file(writers.FortranWriter(filename), 1260 matrix_element) 1261 1262 filename = 'props.inc' 1263 self.write_props_file(writers.FortranWriter(filename), 1264 matrix_element, 1265 s_and_t_channels) 1266 1267 # Find config symmetries and permutations 1268 symmetry, perms, ident_perms = \ 1269 diagram_symmetry.find_symmetry(matrix_element) 1270 1271 filename = 'symswap.inc' 1272 self.write_symswap_file(writers.FortranWriter(filename), 1273 ident_perms) 1274 1275 filename = 'symfact.dat' 1276 self.write_symfact_file(writers.FortranWriter(filename), 1277 symmetry) 1278 1279 # Generate diagrams 1280 filename = "matrix.ps" 1281 plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\ 1282 get('diagrams'), 1283 filename, 1284 model=matrix_element.get('processes')[0].\ 1285 get('model'), 1286 amplitude=True) 1287 logger.info("Generating Feynman diagrams for " + \ 1288 matrix_element.get('processes')[0].nice_string()) 1289 plot.draw() 1290 1291 1292 linkfiles = ['addmothers.f', 1293 'cluster.f', 1294 'cluster.inc', 1295 'coupl.inc', 1296 'cuts.f', 1297 'cuts.inc', 1298 'driver.f', 1299 'genps.f', 1300 'genps.inc', 1301 'initcluster.f', 1302 'makefile', 1303 'message.inc', 1304 'myamp.f', 1305 'reweight.f', 1306 'run.inc', 1307 'maxconfigs.inc', 1308 'maxparticles.inc', 1309 'setcuts.f', 1310 'setscales.f', 1311 'sudakov.inc', 1312 'symmetry.f', 1313 'unwgt.f'] 1314 1315 for file in linkfiles: 1316 ln('../' + file , '.') 1317 1318 #import nexternal/leshouche in Source 1319 ln('nexternal.inc', '../../Source', log=False) 1320 ln('leshouche.inc', '../../Source', log=False) 1321 ln('maxamps.inc', '../../Source', log=False) 1322 1323 # Return to SubProcesses dir 1324 os.chdir(os.path.pardir) 1325 1326 # Add subprocess to subproc.mg 1327 filename = 'subproc.mg' 1328 files.append_to_file(filename, 1329 self.write_subproc, 1330 subprocdir) 1331 1332 # Return to original dir 1333 os.chdir(cwd) 1334 1335 # Generate info page 1336 gen_infohtml.make_info_html(self.dir_path) 1337 1338 1339 if not calls: 1340 calls = 0 1341 return calls
1342
1343 - def finalize_v4_directory(self, matrix_elements, history, makejpg = False, 1344 online = False, compiler='g77'):
1345 """Finalize ME v4 directory by creating jpeg diagrams, html 1346 pages,proc_card_mg5.dat and madevent.tar.gz.""" 1347 1348 modelname = self.model.get('name') 1349 if modelname == 'mssm' or modelname.startswith('mssm-'): 1350 param_card = os.path.join(self.dir_path, 'Cards','param_card.dat') 1351 mg5_param = os.path.join(self.dir_path, 'Source', 'MODEL', 'MG5_param.dat') 1352 check_param_card.convert_to_mg5card(param_card, mg5_param) 1353 check_param_card.check_valid_param_card(mg5_param) 1354 1355 1356 # Write maxconfigs.inc based on max of ME's/subprocess groups 1357 filename = os.path.join(self.dir_path,'Source','maxconfigs.inc') 1358 self.write_maxconfigs_file(writers.FortranWriter(filename), 1359 matrix_elements) 1360 1361 # Write maxparticles.inc based on max of ME's/subprocess groups 1362 filename = os.path.join(self.dir_path,'Source','maxparticles.inc') 1363 self.write_maxparticles_file(writers.FortranWriter(filename), 1364 matrix_elements) 1365 1366 # Touch "done" file 1367 os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses')) 1368 1369 # Check for compiler 1370 self.set_compiler(compiler) 1371 1372 old_pos = os.getcwd() 1373 os.chdir(os.path.join(self.dir_path, 'SubProcesses')) 1374 P_dir_list = [proc for proc in os.listdir('.') if os.path.isdir(proc) and \ 1375 proc[0] == 'P'] 1376 1377 devnull = os.open(os.devnull, os.O_RDWR) 1378 # Convert the poscript in jpg files (if authorize) 1379 if makejpg: 1380 logger.info("Generate jpeg diagrams") 1381 for Pdir in P_dir_list: 1382 os.chdir(Pdir) 1383 try: 1384 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')], 1385 stdout = devnull) 1386 except: 1387 os.system('chmod +x %s ' % os.path.join(old_pos, self.dir_path, 'bin', 'internal', '*')) 1388 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')], 1389 stdout = devnull) 1390 os.chdir(os.path.pardir) 1391 1392 logger.info("Generate web pages") 1393 # Create the WebPage using perl script 1394 1395 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], \ 1396 stdout = devnull) 1397 1398 os.chdir(os.path.pardir) 1399 1400 obj = gen_infohtml.make_info_html(self.dir_path) 1401 [mv(name, './HTML/') for name in os.listdir('.') if \ 1402 (name.endswith('.html') or name.endswith('.jpg')) and \ 1403 name != 'index.html'] 1404 if online: 1405 nb_channel = obj.rep_rule['nb_gen_diag'] 1406 open(os.path.join('./Online'),'w').write(str(nb_channel)) 1407 1408 # Write command history as proc_card_mg5 1409 if os.path.isdir('Cards'): 1410 output_file = os.path.join('Cards', 'proc_card_mg5.dat') 1411 output_file = open(output_file, 'w') 1412 text = ('\n'.join(history) + '\n') % misc.get_time_info() 1413 output_file.write(text) 1414 output_file.close() 1415 1416 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], 1417 stdout = devnull) 1418 1419 # Run "make" to generate madevent.tar.gz file 1420 if os.path.exists(os.path.join('SubProcesses', 'subproc.mg')): 1421 if os.path.exists('madevent.tar.gz'): 1422 os.remove('madevent.tar.gz') 1423 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'make_madevent_tar')], 1424 stdout = devnull) 1425 1426 subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], 1427 stdout = devnull) 1428 1429 #return to the initial dir 1430 os.chdir(old_pos)
1431 1432 #=========================================================================== 1433 # write_matrix_element_v4 1434 #===========================================================================
1435 - def write_matrix_element_v4(self, writer, matrix_element, fortran_model, 1436 proc_id = "", config_map = []):
1437 """Export a matrix element to a matrix.f file in MG4 madevent format""" 1438 1439 if not matrix_element.get('processes') or \ 1440 not matrix_element.get('diagrams'): 1441 return 0 1442 1443 if not isinstance(writer, writers.FortranWriter): 1444 raise writers.FortranWriter.FortranWriterError(\ 1445 "writer not FortranWriter") 1446 1447 # Set lowercase/uppercase Fortran code 1448 writers.FortranWriter.downcase = False 1449 1450 replace_dict = {} 1451 1452 # Extract version number and date from VERSION file 1453 info_lines = self.get_mg5_info_lines() 1454 replace_dict['info_lines'] = info_lines 1455 1456 # Extract process info lines 1457 process_lines = self.get_process_info_lines(matrix_element) 1458 replace_dict['process_lines'] = process_lines 1459 1460 # Set proc_id 1461 replace_dict['proc_id'] = proc_id 1462 1463 # Extract ncomb 1464 ncomb = matrix_element.get_helicity_combinations() 1465 replace_dict['ncomb'] = ncomb 1466 1467 # Extract helicity lines 1468 helicity_lines = self.get_helicity_lines(matrix_element) 1469 replace_dict['helicity_lines'] = helicity_lines 1470 1471 # Extract IC line 1472 ic_line = self.get_ic_line(matrix_element) 1473 replace_dict['ic_line'] = ic_line 1474 1475 # Extract overall denominator 1476 # Averaging initial state color, spin, and identical FS particles 1477 den_factor_line = self.get_den_factor_line(matrix_element) 1478 replace_dict['den_factor_line'] = den_factor_line 1479 1480 # Extract ngraphs 1481 ngraphs = matrix_element.get_number_of_amplitudes() 1482 replace_dict['ngraphs'] = ngraphs 1483 1484 # Extract ndiags 1485 ndiags = len(matrix_element.get('diagrams')) 1486 replace_dict['ndiags'] = ndiags 1487 1488 # Set define_iconfigs_lines 1489 replace_dict['define_iconfigs_lines'] = \ 1490 """INTEGER MAPCONFIG(0:LMAXCONFIGS), ICONFIG 1491 COMMON/TO_MCONFIGS/MAPCONFIG, ICONFIG""" 1492 1493 if proc_id: 1494 # Set lines for subprocess group version 1495 # Set define_iconfigs_lines 1496 replace_dict['define_iconfigs_lines'] += \ 1497 """\nINTEGER SUBDIAG(MAXSPROC),IB(2) 1498 COMMON/TO_SUB_DIAG/SUBDIAG,IB""" 1499 # Set set_amp2_line 1500 replace_dict['set_amp2_line'] = "ANS=ANS*AMP2(SUBDIAG(%s))/XTOT" % \ 1501 proc_id 1502 else: 1503 # Standard running 1504 # Set set_amp2_line 1505 replace_dict['set_amp2_line'] = "ANS=ANS*AMP2(MAPCONFIG(ICONFIG))/XTOT" 1506 1507 # Extract nwavefuncs 1508 nwavefuncs = matrix_element.get_number_of_wavefunctions() 1509 replace_dict['nwavefuncs'] = nwavefuncs 1510 1511 # Extract ncolor 1512 ncolor = max(1, len(matrix_element.get('color_basis'))) 1513 replace_dict['ncolor'] = ncolor 1514 1515 # Extract color data lines 1516 color_data_lines = self.get_color_data_lines(matrix_element) 1517 replace_dict['color_data_lines'] = "\n".join(color_data_lines) 1518 1519 # Extract helas calls 1520 helas_calls = fortran_model.get_matrix_element_calls(\ 1521 matrix_element) 1522 replace_dict['helas_calls'] = "\n".join(helas_calls) 1523 1524 # Set the size of Wavefunction 1525 if not self.model or any([p.get('spin')==5 for p in self.model.get('particles') if p]): 1526 replace_dict['wavefunctionsize'] = 18 1527 else: 1528 replace_dict['wavefunctionsize'] = 6 1529 1530 # Extract amp2 lines 1531 amp2_lines = self.get_amp2_lines(matrix_element, config_map) 1532 replace_dict['amp2_lines'] = '\n'.join(amp2_lines) 1533 1534 # Extract JAMP lines 1535 jamp_lines = self.get_JAMP_lines(matrix_element) 1536 replace_dict['jamp_lines'] = '\n'.join(jamp_lines) 1537 1538 file = open(os.path.join(_file_path, \ 1539 'iolibs/template_files/%s' % self.matrix_file)).read() 1540 file = file % replace_dict 1541 1542 1543 # Write the file 1544 writer.writelines(file) 1545 1546 return len(filter(lambda call: call.find('#') != 0, helas_calls)), ncolor
1547 1548 #=========================================================================== 1549 # write_auto_dsig_file 1550 #===========================================================================
1551 - def write_auto_dsig_file(self, writer, matrix_element, proc_id = ""):
1552 """Write the auto_dsig.f file for the differential cross section 1553 calculation, includes pdf call information""" 1554 1555 if not matrix_element.get('processes') or \ 1556 not matrix_element.get('diagrams'): 1557 return 0 1558 1559 nexternal, ninitial = matrix_element.get_nexternal_ninitial() 1560 1561 if ninitial < 1 or ninitial > 2: 1562 raise writers.FortranWriter.FortranWriterError, \ 1563 """Need ninitial = 1 or 2 to write auto_dsig file""" 1564 1565 replace_dict = {} 1566 1567 # Extract version number and date from VERSION file 1568 info_lines = self.get_mg5_info_lines() 1569 replace_dict['info_lines'] = info_lines 1570 1571 # Extract process info lines 1572 process_lines = self.get_process_info_lines(matrix_element) 1573 replace_dict['process_lines'] = process_lines 1574 1575 # Set proc_id 1576 replace_dict['proc_id'] = proc_id 1577 replace_dict['numproc'] = 1 1578 1579 # Set dsig_line 1580 if ninitial == 1: 1581 # No conversion, since result of decay should be given in GeV 1582 dsig_line = "pd(IPROC)*dsiguu" 1583 else: 1584 # Convert result (in GeV) to pb 1585 dsig_line = "pd(IPROC)*conv*dsiguu" 1586 1587 replace_dict['dsig_line'] = dsig_line 1588 1589 # Extract pdf lines 1590 pdf_vars, pdf_data, pdf_lines = \ 1591 self.get_pdf_lines(matrix_element, ninitial, proc_id != "") 1592 replace_dict['pdf_vars'] = pdf_vars 1593 replace_dict['pdf_data'] = pdf_data 1594 replace_dict['pdf_lines'] = pdf_lines 1595 1596 # Lines that differ between subprocess group and regular 1597 if proc_id: 1598 replace_dict['numproc'] = int(proc_id) 1599 replace_dict['passcuts_begin'] = "" 1600 replace_dict['passcuts_end'] = "" 1601 # Set lines for subprocess group version 1602 # Set define_iconfigs_lines 1603 replace_dict['define_subdiag_lines'] = \ 1604 """\nINTEGER SUBDIAG(MAXSPROC),IB(2) 1605 COMMON/TO_SUB_DIAG/SUBDIAG,IB""" 1606 else: 1607 replace_dict['passcuts_begin'] = "IF (PASSCUTS(PP)) THEN" 1608 replace_dict['passcuts_end'] = "ENDIF" 1609 replace_dict['define_subdiag_lines'] = "" 1610 1611 file = open(os.path.join(_file_path, \ 1612 'iolibs/template_files/auto_dsig_v4.inc')).read() 1613 file = file % replace_dict 1614 1615 # Write the file 1616 writer.writelines(file)
1617 1618 #=========================================================================== 1619 # write_coloramps_file 1620 #===========================================================================
1621 - def write_coloramps_file(self, writer, mapconfigs, matrix_element):
1622 """Write the coloramps.inc file for MadEvent""" 1623 1624 lines = self.get_icolamp_lines(mapconfigs, matrix_element, 1) 1625 lines.insert(0, "logical icolamp(%d,%d,1)" % \ 1626 (max(len(matrix_element.get('color_basis').keys()), 1), 1627 len(mapconfigs))) 1628 1629 1630 # Write the file 1631 writer.writelines(lines) 1632 1633 return True
1634 1635 #=========================================================================== 1636 # write_coloramps_file 1637 #===========================================================================
1638 - def write_colors_file(self, writer, matrix_elements):
1639 """Write the get_color.f file for MadEvent, which returns color 1640 for all particles used in the matrix element.""" 1641 1642 if isinstance(matrix_elements, helas_objects.HelasMatrixElement): 1643 matrix_elements = [matrix_elements] 1644 1645 model = matrix_elements[0].get('processes')[0].get('model') 1646 1647 # We need the both particle and antiparticle wf_ids, since the identity 1648 # depends on the direction of the wf. 1649 wf_ids = set(sum([sum([sum([[wf.get_pdg_code(),wf.get_anti_pdg_code()] \ 1650 for wf in d.get('wavefunctions')],[]) \ 1651 for d in me.get('diagrams')], []) \ 1652 for me in matrix_elements], [])) 1653 1654 leg_ids = set(sum([sum([[l.get('id') for l in \ 1655 p.get_legs_with_decays()] for p in \ 1656 me.get('processes')], []) for me in 1657 matrix_elements], [])) 1658 particle_ids = sorted(list(wf_ids.union(leg_ids))) 1659 1660 lines = """function get_color(ipdg) 1661 implicit none 1662 integer get_color, ipdg 1663 1664 if(ipdg.eq.%d)then 1665 get_color=%d 1666 return 1667 """ % (particle_ids[0], model.get_particle(particle_ids[0]).get_color()) 1668 1669 for part_id in particle_ids[1:]: 1670 lines += """else if(ipdg.eq.%d)then 1671 get_color=%d 1672 return 1673 """ % (part_id, model.get_particle(part_id).get_color()) 1674 # Dummy particle for multiparticle vertices with pdg given by 1675 # first code not in the model 1676 lines += """else if(ipdg.eq.%d)then 1677 c This is dummy particle used in multiparticle vertices 1678 get_color=2 1679 return 1680 """ % model.get_first_non_pdg() 1681 lines += """else 1682 write(*,*)'Error: No color given for pdg ',ipdg 1683 get_color=0 1684 return 1685 endif 1686 end 1687 """ 1688 1689 # Write the file 1690 writer.writelines(lines) 1691 1692 return True
1693 1694 #=========================================================================== 1695 # write_maxconfigs_file 1696 #===========================================================================
1697 - def write_maxconfigs_file(self, writer, matrix_elements):
1698 """Write the maxconfigs.inc file for MadEvent""" 1699 1700 if isinstance(matrix_elements, helas_objects.HelasMultiProcess): 1701 maxconfigs = max([me.get_num_configs() for me in \ 1702 matrix_elements.get('matrix_elements')]) 1703 else: 1704 maxconfigs = max([me.get_num_configs() for me in matrix_elements]) 1705 1706 lines = "integer lmaxconfigs\n" 1707 lines += "parameter(lmaxconfigs=%d)" % maxconfigs 1708 1709 # Write the file 1710 writer.writelines(lines) 1711 1712 return True
1713 1714 #=========================================================================== 1715 # write_maxparticles_file 1716 #===========================================================================
1717 - def write_maxparticles_file(self, writer, matrix_elements):
1718 """Write the maxparticles.inc file for MadEvent""" 1719 1720 if isinstance(matrix_elements, helas_objects.HelasMultiProcess): 1721 maxparticles = max([me.get_nexternal_ninitial()[0] for me in \ 1722 matrix_elements.get('matrix_elements')]) 1723 else: 1724 maxparticles = max([me.get_nexternal_ninitial()[0] \ 1725 for me in matrix_elements]) 1726 1727 lines = "integer max_particles\n" 1728 lines += "parameter(max_particles=%d)" % maxparticles 1729 1730 # Write the file 1731 writer.writelines(lines) 1732 1733 return True
1734 1735 1736 #=========================================================================== 1737 # write_configs_file 1738 #===========================================================================
1739 - def write_configs_file(self, writer, matrix_element):
1740 """Write the configs.inc file for MadEvent""" 1741 1742 # Extract number of external particles 1743 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 1744 1745 configs = [(i+1, d) for i,d in enumerate(matrix_element.get('diagrams'))] 1746 mapconfigs = [c[0] for c in configs] 1747 model = matrix_element.get('processes')[0].get('model') 1748 return mapconfigs, self.write_configs_file_from_diagrams(writer, 1749 [[c[1]] for c in configs], 1750 mapconfigs, 1751 nexternal, ninitial, 1752 model)
1753 1754 1755 1756 #=========================================================================== 1757 # write_run_configs_file 1758 #===========================================================================
1759 - def write_run_config_file(self, writer):
1760 """Write the run_configs.inc file for MadEvent""" 1761 1762 path = os.path.join(_file_path,'iolibs','template_files','madevent_run_config.inc') 1763 text = open(path).read() % {'chanperjob':'5'} 1764 writer.write(text) 1765 return True
1766 1767 1768 #=========================================================================== 1769 # write_configs_file_from_diagrams 1770 #===========================================================================
1771 - def write_configs_file_from_diagrams(self, writer, configs, mapconfigs, 1772 nexternal, ninitial, model):
1773 """Write the actual configs.inc file. 1774 1775 configs is the diagrams corresponding to configs (each 1776 diagrams is a list of corresponding diagrams for all 1777 subprocesses, with None if there is no corresponding diagrams 1778 for a given process). 1779 mapconfigs gives the diagram number for each config. 1780 1781 For s-channels, we need to output one PDG for each subprocess in 1782 the subprocess group, in order to be able to pick the right 1783 one for multiprocesses.""" 1784 1785 lines = [] 1786 1787 s_and_t_channels = [] 1788 1789 minvert = min([max([d for d in config if d][0].get_vertex_leg_numbers()) \ 1790 for config in configs]) 1791 1792 # Number of subprocesses 1793 nsubprocs = len(configs[0]) 1794 1795 nconfigs = 0 1796 1797 new_pdg = model.get_first_non_pdg() 1798 1799 for iconfig, helas_diags in enumerate(configs): 1800 if any([vert > minvert for vert in 1801 [d for d in helas_diags if d][0].get_vertex_leg_numbers()]): 1802 # Only 3-vertices allowed in configs.inc 1803 continue 1804 nconfigs += 1 1805 1806 # Need s- and t-channels for all subprocesses, including 1807 # those that don't contribute to this config 1808 empty_verts = [] 1809 stchannels = [] 1810 for h in helas_diags: 1811 if h: 1812 # get_s_and_t_channels gives vertices starting from 1813 # final state external particles and working inwards 1814 stchannels.append(h.get('amplitudes')[0].\ 1815 get_s_and_t_channels(ninitial, new_pdg)) 1816 else: 1817 stchannels.append((empty_verts, None)) 1818 1819 # For t-channels, just need the first non-empty one 1820 tchannels = [t for s,t in stchannels if t != None][0] 1821 1822 # For s_and_t_channels (to be used later) use only first config 1823 s_and_t_channels.append([[s for s,t in stchannels if t != None][0], 1824 tchannels]) 1825 1826 # Make sure empty_verts is same length as real vertices 1827 if any([s for s,t in stchannels]): 1828 empty_verts[:] = [None]*max([len(s) for s,t in stchannels]) 1829 1830 # Reorganize s-channel vertices to get a list of all 1831 # subprocesses for each vertex 1832 schannels = zip(*[s for s,t in stchannels]) 1833 else: 1834 schannels = [] 1835 1836 allchannels = schannels 1837 if len(tchannels) > 1: 1838 # Write out tchannels only if there are any non-trivial ones 1839 allchannels = schannels + tchannels 1840 1841 # Write out propagators for s-channel and t-channel vertices 1842 1843 lines.append("# Diagram %d" % (mapconfigs[iconfig])) 1844 # Correspondance between the config and the diagram = amp2 1845 lines.append("data mapconfig(%d)/%d/" % (nconfigs, 1846 mapconfigs[iconfig])) 1847 1848 for verts in allchannels: 1849 if verts in schannels: 1850 vert = [v for v in verts if v][0] 1851 else: 1852 vert = verts 1853 daughters = [leg.get('number') for leg in vert.get('legs')[:-1]] 1854 last_leg = vert.get('legs')[-1] 1855 lines.append("data (iforest(i,%d,%d),i=1,%d)/%s/" % \ 1856 (last_leg.get('number'), nconfigs, len(daughters), 1857 ",".join([str(d) for d in daughters]))) 1858 if verts in schannels: 1859 pdgs = [] 1860 for v in verts: 1861 if v: 1862 pdgs.append(v.get('legs')[-1].get('id')) 1863 else: 1864 pdgs.append(0) 1865 lines.append("data (sprop(i,%d,%d),i=1,%d)/%s/" % \ 1866 (last_leg.get('number'), nconfigs, nsubprocs, 1867 ",".join([str(d) for d in pdgs]))) 1868 lines.append("data tprid(%d,%d)/0/" % \ 1869 (last_leg.get('number'), nconfigs)) 1870 elif verts in tchannels[:-1]: 1871 lines.append("data tprid(%d,%d)/%d/" % \ 1872 (last_leg.get('number'), nconfigs, 1873 abs(last_leg.get('id')))) 1874 lines.append("data (sprop(i,%d,%d),i=1,%d)/%s/" % \ 1875 (last_leg.get('number'), nconfigs, nsubprocs, 1876 ",".join(['0'] * nsubprocs))) 1877 1878 # Write out number of configs 1879 lines.append("# Number of configs") 1880 lines.append("data mapconfig(0)/%d/" % nconfigs) 1881 1882 # Write the file 1883 writer.writelines(lines) 1884 1885 return s_and_t_channels
1886 1887 #=========================================================================== 1888 # write_config_subproc_map_file 1889 #===========================================================================
1890 - def write_config_subproc_map_file(self, writer, s_and_t_channels):
1891 """Write a dummy config_subproc.inc file for MadEvent""" 1892 1893 lines = [] 1894 1895 for iconfig in range(len(s_and_t_channels)): 1896 lines.append("DATA CONFSUB(1,%d)/1/" % \ 1897 (iconfig + 1)) 1898 1899 # Write the file 1900 writer.writelines(lines) 1901 1902 return True
1903 1904 #=========================================================================== 1905 # write_decayBW_file 1906 #===========================================================================
1907 - def write_decayBW_file(self, writer, s_and_t_channels):
1908 """Write the decayBW.inc file for MadEvent""" 1909 1910 lines = [] 1911 1912 booldict = {None: "0", True: "1", False: "2"} 1913 1914 for iconf, config in enumerate(s_and_t_channels): 1915 schannels = config[0] 1916 for vertex in schannels: 1917 # For the resulting leg, pick out whether it comes from 1918 # decay or not, as given by the onshell flag 1919 leg = vertex.get('legs')[-1] 1920 lines.append("data gForceBW(%d,%d)/%s/" % \ 1921 (leg.get('number'), iconf + 1, 1922 booldict[leg.get('onshell')])) 1923 1924 # Write the file 1925 writer.writelines(lines) 1926 1927 return True
1928 1929 #=========================================================================== 1930 # write_dname_file 1931 #===========================================================================
1932 - def write_dname_file(self, writer, dir_name):
1933 """Write the dname.mg file for MG4""" 1934 1935 line = "DIRNAME=%s" % dir_name 1936 1937 # Write the file 1938 writer.write(line + "\n") 1939 1940 return True
1941 1942 #=========================================================================== 1943 # write_driver 1944 #===========================================================================
1945 - def write_driver(self, writer, v5=True):
1946 """Write the SubProcess/driver.f file for MG4""" 1947 1948 path = os.path.join(_file_path,'iolibs','template_files','madevent_driver.f') 1949 1950 if self.model_name == 'mssm' or self.model_name.startswith('mssm-'): 1951 card = 'Source/MODEL/MG5_param.dat' 1952 else: 1953 card = 'param_card.dat' 1954 if v5: 1955 text = open(path).read() % {'param_card_name':card, 'secondparam':''} 1956 else: 1957 text = open(path).read() % {'param_card_name':card, 1958 'secondparam': ',.true.'} 1959 writer.write(text) 1960 1961 return True
1962 1963 #=========================================================================== 1964 # write_combine_events 1965 #===========================================================================
1966 - def write_combine_events(self, writer):
1967 """Write the SubProcess/driver.f file for MG4""" 1968 1969 path = os.path.join(_file_path,'iolibs','template_files','madevent_combine_events.f') 1970 1971 if self.model_name == 'mssm' or self.model_name.startswith('mssm-'): 1972 card = 'Source/MODEL/MG5_param.dat' 1973 else: 1974 card = 'param_card.dat' 1975 text = open(path).read() % {'param_card_name':card} 1976 1977 writer.write(text) 1978 1979 return True
1980 1981 1982 #=========================================================================== 1983 # write_symmetry 1984 #===========================================================================
1985 - def write_symmetry(self, writer, v5=True):
1986 """Write the SubProcess/driver.f file for ME""" 1987 1988 path = os.path.join(_file_path,'iolibs','template_files','madevent_symmetry.f') 1989 1990 if self.model_name == 'mssm' or self.model_name.startswith('mssm-'): 1991 card = 'Source/MODEL/MG5_param.dat' 1992 else: 1993 card = 'param_card.dat' 1994 text = open(path).read() 1995 1996 if v5: 1997 text = text % {'param_card_name':card, 'setparasecondarg':''} 1998 else: 1999 text = text % {'param_card_name':card, 'setparasecondarg':',.true.'} 2000 writer.write(text) 2001 2002 return True
2003 2004 2005 2006 2007 #=========================================================================== 2008 # write_iproc_file 2009 #===========================================================================
2010 - def write_iproc_file(self, writer, me_number):
2011 """Write the iproc.dat file for MG4""" 2012 line = "%d" % (me_number + 1) 2013 2014 # Write the file 2015 for line_to_write in writer.write_line(line): 2016 writer.write(line_to_write) 2017 return True
2018 2019 #=========================================================================== 2020 # write_leshouche_file 2021 #===========================================================================
2022 - def write_leshouche_file(self, writer, matrix_element):
2023 """Write the leshouche.inc file for MG4""" 2024 2025 # Write the file 2026 writer.writelines(self.get_leshouche_lines(matrix_element, 0)) 2027 2028 return True
2029 2030 #=========================================================================== 2031 # get_leshouche_lines 2032 #===========================================================================
2033 - def get_leshouche_lines(self, matrix_element, numproc):
2034 """Write the leshouche.inc file for MG4""" 2035 2036 # Extract number of external particles 2037 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 2038 2039 lines = [] 2040 for iproc, proc in enumerate(matrix_element.get('processes')): 2041 legs = proc.get_legs_with_decays() 2042 lines.append("DATA (IDUP(i,%d,%d),i=1,%d)/%s/" % \ 2043 (iproc + 1, numproc+1, nexternal, 2044 ",".join([str(l.get('id')) for l in legs]))) 2045 if iproc == 0 and numproc == 0: 2046 for i in [1, 2]: 2047 lines.append("DATA (MOTHUP(%d,i),i=1,%2r)/%s/" % \ 2048 (i, nexternal, 2049 ",".join([ "%3r" % 0 ] * ninitial + \ 2050 [ "%3r" % i ] * (nexternal - ninitial)))) 2051 2052 # Here goes the color connections corresponding to the JAMPs 2053 # Only one output, for the first subproc! 2054 if iproc == 0: 2055 # If no color basis, just output trivial color flow 2056 if not matrix_element.get('color_basis'): 2057 for i in [1, 2]: 2058 lines.append("DATA (ICOLUP(%d,i,1,%d),i=1,%2r)/%s/" % \ 2059 (i, numproc+1,nexternal, 2060 ",".join([ "%3r" % 0 ] * nexternal))) 2061 2062 else: 2063 # First build a color representation dictionnary 2064 repr_dict = {} 2065 for l in legs: 2066 repr_dict[l.get('number')] = \ 2067 proc.get('model').get_particle(l.get('id')).get_color()\ 2068 * (-1)**(1+l.get('state')) 2069 # Get the list of color flows 2070 color_flow_list = \ 2071 matrix_element.get('color_basis').color_flow_decomposition(repr_dict, 2072 ninitial) 2073 # And output them properly 2074 for cf_i, color_flow_dict in enumerate(color_flow_list): 2075 for i in [0, 1]: 2076 lines.append("DATA (ICOLUP(%d,i,%d,%d),i=1,%2r)/%s/" % \ 2077 (i + 1, cf_i + 1, numproc+1, nexternal, 2078 ",".join(["%3r" % color_flow_dict[l.get('number')][i] \ 2079 for l in legs]))) 2080 2081 return lines
2082 2083 #=========================================================================== 2084 # write_maxamps_file 2085 #===========================================================================
2086 - def write_maxamps_file(self, writer, maxamps, maxflows, 2087 maxproc,maxsproc):
2088 """Write the maxamps.inc file for MG4.""" 2089 2090 file = " integer maxamps, maxflow, maxproc, maxsproc\n" 2091 file = file + "parameter (maxamps=%d, maxflow=%d)\n" % \ 2092 (maxamps, maxflows) 2093 file = file + "parameter (maxproc=%d, maxsproc=%d)" % \ 2094 (maxproc, maxsproc) 2095 2096 # Write the file 2097 writer.writelines(file) 2098 2099 return True
2100 2101 #=========================================================================== 2102 # write_mg_sym_file 2103 #===========================================================================
2104 - def write_mg_sym_file(self, writer, matrix_element):
2105 """Write the mg.sym file for MadEvent.""" 2106 2107 lines = [] 2108 2109 # Extract process with all decays included 2110 final_legs = filter(lambda leg: leg.get('state') == True, 2111 matrix_element.get('processes')[0].get_legs_with_decays()) 2112 2113 ninitial = len(filter(lambda leg: leg.get('state') == False, 2114 matrix_element.get('processes')[0].get('legs'))) 2115 2116 identical_indices = {} 2117 2118 # Extract identical particle info 2119 for i, leg in enumerate(final_legs): 2120 if leg.get('id') in identical_indices: 2121 identical_indices[leg.get('id')].append(\ 2122 i + ninitial + 1) 2123 else: 2124 identical_indices[leg.get('id')] = [i + ninitial + 1] 2125 2126 # Remove keys which have only one particle 2127 for key in identical_indices.keys(): 2128 if len(identical_indices[key]) < 2: 2129 del identical_indices[key] 2130 2131 # Write mg.sym file 2132 lines.append(str(len(identical_indices.keys()))) 2133 for key in identical_indices.keys(): 2134 lines.append(str(len(identical_indices[key]))) 2135 for number in identical_indices[key]: 2136 lines.append(str(number)) 2137 2138 # Write the file 2139 writer.writelines(lines) 2140 2141 return True
2142 2143 #=========================================================================== 2144 # write_mg_sym_file 2145 #===========================================================================
2146 - def write_default_mg_sym_file(self, writer):
2147 """Write the mg.sym file for MadEvent.""" 2148 2149 lines = "0" 2150 2151 # Write the file 2152 writer.writelines(lines) 2153 2154 return True
2155 2156 #=========================================================================== 2157 # write_ncombs_file 2158 #===========================================================================
2159 - def write_ncombs_file(self, writer, nexternal):
2160 """Write the ncombs.inc file for MadEvent.""" 2161 2162 # ncomb (used for clustering) is 2^nexternal 2163 file = " integer n_max_cl\n" 2164 file = file + "parameter (n_max_cl=%d)" % (2 ** nexternal) 2165 2166 # Write the file 2167 writer.writelines(file) 2168 2169 return True
2170 2171 #=========================================================================== 2172 # write_props_file 2173 #===========================================================================
2174 - def write_props_file(self, writer, matrix_element, s_and_t_channels):
2175 """Write the props.inc file for MadEvent. Needs input from 2176 write_configs_file.""" 2177 2178 lines = [] 2179 2180 particle_dict = matrix_element.get('processes')[0].get('model').\ 2181 get('particle_dict') 2182 2183 for iconf, configs in enumerate(s_and_t_channels): 2184 for vertex in configs[0] + configs[1][:-1]: 2185 leg = vertex.get('legs')[-1] 2186 if leg.get('id') not in particle_dict: 2187 # Fake propagator used in multiparticle vertices 2188 mass = 'zero' 2189 width = 'zero' 2190 pow_part = 0 2191 else: 2192 particle = particle_dict[leg.get('id')] 2193 # Get mass 2194 if particle.get('mass').lower() == 'zero': 2195 mass = particle.get('mass') 2196 else: 2197 mass = "abs(%s)" % particle.get('mass') 2198 # Get width 2199 if particle.get('width').lower() == 'zero': 2200 width = particle.get('width') 2201 else: 2202 width = "abs(%s)" % particle.get('width') 2203 2204 pow_part = 1 + int(particle.is_boson()) 2205 2206 lines.append("prmass(%d,%d) = %s" % \ 2207 (leg.get('number'), iconf + 1, mass)) 2208 lines.append("prwidth(%d,%d) = %s" % \ 2209 (leg.get('number'), iconf + 1, width)) 2210 lines.append("pow(%d,%d) = %d" % \ 2211 (leg.get('number'), iconf + 1, pow_part)) 2212 2213 # Write the file 2214 writer.writelines(lines) 2215 2216 return True
2217 2218 #=========================================================================== 2219 # write_processes_file 2220 #===========================================================================
2221 - def write_processes_file(self, writer, subproc_group):
2222 """Write the processes.dat file with info about the subprocesses 2223 in this group.""" 2224 2225 lines = [] 2226 2227 for ime, me in \ 2228 enumerate(subproc_group.get('matrix_elements')): 2229 lines.append("%s %s" % (str(ime+1) + " " * (7-len(str(ime+1))), 2230 ",".join(p.base_string() for p in \ 2231 me.get('processes')))) 2232 if me.get('has_mirror_process'): 2233 mirror_procs = [copy.copy(p) for p in me.get('processes')] 2234 for proc in mirror_procs: 2235 legs = copy.copy(proc.get('legs')) 2236 legs.insert(0, legs.pop(1)) 2237 proc.set("legs", legs) 2238 lines.append("mirror %s" % ",".join(p.base_string() for p in \ 2239 mirror_procs)) 2240 else: 2241 lines.append("mirror none") 2242 2243 # Write the file 2244 writer.write("\n".join(lines)) 2245 2246 return True
2247 2248 #=========================================================================== 2249 # write_symswap_file 2250 #===========================================================================
2251 - def write_symswap_file(self, writer, ident_perms):
2252 """Write the file symswap.inc for MG4 by comparing diagrams using 2253 the internal matrix element value functionality.""" 2254 2255 lines = [] 2256 2257 # Write out lines for symswap.inc file (used to permute the 2258 # external leg momenta 2259 for iperm, perm in enumerate(ident_perms): 2260 lines.append("data (isym(i,%d),i=1,nexternal)/%s/" % \ 2261 (iperm+1, ",".join([str(i+1) for i in perm]))) 2262 lines.append("data nsym/%d/" % len(ident_perms)) 2263 2264 # Write the file 2265 writer.writelines(lines) 2266 2267 return True
2268 2269 #=========================================================================== 2270 # write_symfact_file 2271 #===========================================================================
2272 - def write_symfact_file(self, writer, symmetry):
2273 """Write the files symfact.dat for MG4 by comparing diagrams using 2274 the internal matrix element value functionality.""" 2275 2276 pos = max(2, int(math.ceil(math.log10(len(symmetry))))) 2277 form = "%"+str(pos)+"r %"+str(pos+1)+"r" 2278 # Write out lines for symswap.inc file (used to permute the 2279 # external leg momenta 2280 lines = [ form %(i+1, s) for i,s in enumerate(symmetry) if s != 0] 2281 2282 # Write the file 2283 writer.writelines(lines) 2284 2285 return True
2286 2287 #=========================================================================== 2288 # write_symperms_file 2289 #===========================================================================
2290 - def write_symperms_file(self, writer, perms):
2291 """Write the symperms.inc file for subprocess group, used for 2292 symmetric configurations""" 2293 2294 lines = [] 2295 for iperm, perm in enumerate(perms): 2296 lines.append("data (perms(i,%d),i=1,nexternal)/%s/" % \ 2297 (iperm+1, ",".join([str(i+1) for i in perm]))) 2298 2299 # Write the file 2300 writer.writelines(lines) 2301 2302 return True
2303 2304 #=========================================================================== 2305 # write_subproc 2306 #===========================================================================
2307 - def write_subproc(self, writer, subprocdir):
2308 """Append this subprocess to the subproc.mg file for MG4""" 2309 2310 # Write line to file 2311 writer.write(subprocdir + "\n") 2312 2313 return True
2314 2315 2316 2317 #=============================================================================== 2318 # ProcessExporterFortranMEGroup 2319 #===============================================================================
2320 -class ProcessExporterFortranMEGroup(ProcessExporterFortranME):
2321 """Class to take care of exporting a set of matrix elements to 2322 MadEvent subprocess group format.""" 2323 2324 matrix_file = "matrix_madevent_group_v4.inc" 2325 2326 #=========================================================================== 2327 # generate_subprocess_directory_v4 2328 #===========================================================================
2329 - def generate_subprocess_directory_v4(self, subproc_group, 2330 fortran_model, 2331 group_number):
2332 """Generate the Pn directory for a subprocess group in MadEvent, 2333 including the necessary matrix_N.f files, configs.inc and various 2334 other helper files""" 2335 2336 if not isinstance(subproc_group, group_subprocs.SubProcessGroup): 2337 raise base_objects.PhysicsObject.PhysicsObjectError,\ 2338 "subproc_group object not SubProcessGroup" 2339 2340 if not self.model: 2341 self.model = subproc_group.get('matrix_elements')[0].\ 2342 get('processes')[0].get('model') 2343 2344 cwd = os.getcwd() 2345 path = os.path.join(self.dir_path, 'SubProcesses') 2346 2347 os.chdir(path) 2348 2349 pathdir = os.getcwd() 2350 2351 # Create the directory PN in the specified path 2352 subprocdir = "P%d_%s" % (subproc_group.get('number'), 2353 subproc_group.get('name')) 2354 try: 2355 os.mkdir(subprocdir) 2356 except os.error as error: 2357 logger.warning(error.strerror + " " + subprocdir) 2358 2359 try: 2360 os.chdir(subprocdir) 2361 except os.error: 2362 logger.error('Could not cd to directory %s' % subprocdir) 2363 return 0 2364 2365 logger.info('Creating files in directory %s' % subprocdir) 2366 2367 # Create the matrix.f files, auto_dsig.f files and all inc files 2368 # for all subprocesses in the group 2369 2370 maxamps = 0 2371 maxflows = 0 2372 tot_calls = 0 2373 2374 matrix_elements = subproc_group.get('matrix_elements') 2375 2376 for ime, matrix_element in \ 2377 enumerate(matrix_elements): 2378 filename = 'matrix%d.f' % (ime+1) 2379 calls, ncolor = \ 2380 self.write_matrix_element_v4(writers.FortranWriter(filename), 2381 matrix_element, 2382 fortran_model, 2383 str(ime+1), 2384 subproc_group.get('diagram_maps')[\ 2385 ime]) 2386 2387 filename = 'auto_dsig%d.f' % (ime+1) 2388 self.write_auto_dsig_file(writers.FortranWriter(filename), 2389 matrix_element, 2390 str(ime+1)) 2391 2392 # Keep track of needed quantities 2393 tot_calls += int(calls) 2394 maxflows = max(maxflows, ncolor) 2395 maxamps = max(maxamps, len(matrix_element.get('diagrams'))) 2396 2397 # Draw diagrams 2398 filename = "matrix%d.ps" % (ime+1) 2399 plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\ 2400 get('diagrams'), 2401 filename, 2402 model = \ 2403 matrix_element.get('processes')[0].\ 2404 get('model'), 2405 amplitude=True) 2406 logger.info("Generating Feynman diagrams for " + \ 2407 matrix_element.get('processes')[0].nice_string()) 2408 plot.draw() 2409 2410 # Extract number of external particles 2411 (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() 2412 2413 # Generate a list of diagrams corresponding to each configuration 2414 # [[d1, d2, ...,dn],...] where 1,2,...,n is the subprocess number 2415 # If a subprocess has no diagrams for this config, the number is 0 2416 2417 subproc_diagrams_for_config = subproc_group.get('diagrams_for_configs') 2418 2419 filename = 'auto_dsig.f' 2420 self.write_super_auto_dsig_file(writers.FortranWriter(filename), 2421 subproc_group) 2422 2423 filename = 'coloramps.inc' 2424 self.write_coloramps_file(writers.FortranWriter(filename), 2425 subproc_diagrams_for_config, 2426 maxflows, 2427 matrix_elements) 2428 2429 filename = 'get_color.f' 2430 self.write_colors_file(writers.FortranWriter(filename), 2431 matrix_elements) 2432 2433 filename = 'config_subproc_map.inc' 2434 self.write_config_subproc_map_file(writers.FortranWriter(filename), 2435 subproc_diagrams_for_config) 2436 2437 filename = 'configs.inc' 2438 nconfigs, s_and_t_channels = self.write_configs_file(\ 2439 writers.FortranWriter(filename), 2440 subproc_group, 2441 subproc_diagrams_for_config) 2442 2443 filename = 'decayBW.inc' 2444 self.write_decayBW_file(writers.FortranWriter(filename), 2445 s_and_t_channels) 2446 2447 filename = 'dname.mg' 2448 self.write_dname_file(writers.FortranWriter(filename), 2449 subprocdir) 2450 2451 filename = 'iproc.dat' 2452 self.write_iproc_file(writers.FortranWriter(filename), 2453 group_number) 2454 2455 filename = 'leshouche.inc' 2456 self.write_leshouche_file(writers.FortranWriter(filename), 2457 subproc_group) 2458 2459 filename = 'maxamps.inc' 2460 self.write_maxamps_file(writers.FortranWriter(filename), 2461 maxamps, 2462 maxflows, 2463 max([len(me.get('processes')) for me in \ 2464 matrix_elements]), 2465 len(matrix_elements)) 2466 2467 # Note that mg.sym is not relevant for this case 2468 filename = 'mg.sym' 2469 self.write_default_mg_sym_file(writers.FortranWriter(filename)) 2470 2471 filename = 'mirrorprocs.inc' 2472 self.write_mirrorprocs(writers.FortranWriter(filename), 2473 subproc_group) 2474 2475 filename = 'ncombs.inc' 2476 self.write_ncombs_file(writers.FortranWriter(filename), 2477 nexternal) 2478 2479 filename = 'nexternal.inc' 2480 self.write_nexternal_file(writers.FortranWriter(filename), 2481 nexternal, ninitial) 2482 2483 filename = 'ngraphs.inc' 2484 self.write_ngraphs_file(writers.FortranWriter(filename), 2485 nconfigs) 2486 2487 filename = 'pmass.inc' 2488 self.write_pmass_file(writers.FortranWriter(filename), 2489 matrix_element) 2490 2491 filename = 'props.inc' 2492 self.write_props_file(writers.FortranWriter(filename), 2493 matrix_element, 2494 s_and_t_channels) 2495 2496 filename = 'processes.dat' 2497 files.write_to_file(filename, 2498 self.write_processes_file, 2499 subproc_group) 2500 2501 # Find config symmetries and permutations 2502 symmetry, perms, ident_perms = \ 2503 diagram_symmetry.find_symmetry(subproc_group) 2504 2505 filename = 'symswap.inc' 2506 self.write_symswap_file(writers.FortranWriter(filename), 2507 ident_perms) 2508 2509 filename = 'symfact.dat' 2510 self.write_symfact_file(writers.FortranWriter(filename), 2511 symmetry) 2512 2513 filename = 'symperms.inc' 2514 self.write_symperms_file(writers.FortranWriter(filename), 2515 perms) 2516 2517 # Generate jpgs -> pass in make_html 2518 #os.system(os.path.join('..', '..', 'bin', 'gen_jpeg-pl')) 2519 2520 linkfiles = ['addmothers.f', 2521 'cluster.f', 2522 'cluster.inc', 2523 'coupl.inc', 2524 'cuts.f', 2525 'cuts.inc', 2526 'driver.f', 2527 'genps.f', 2528 'genps.inc', 2529 'initcluster.f', 2530 'makefile', 2531 'message.inc', 2532 'myamp.f', 2533 'reweight.f', 2534 'run.inc', 2535 'maxconfigs.inc', 2536 'maxparticles.inc', 2537 'setcuts.f', 2538 'setscales.f', 2539 'sudakov.inc', 2540 'symmetry.f', 2541 'unwgt.f'] 2542 2543 for file in linkfiles: 2544 ln('../' + file , '.') 2545 2546 #import nexternal/leshouch in Source 2547 ln('nexternal.inc', '../../Source', log=False) 2548 ln('leshouche.inc', '../../Source', log=False) 2549 ln('maxamps.inc', '../../Source', log=False) 2550 2551 # Return to SubProcesses dir 2552 os.chdir(pathdir) 2553 2554 # Add subprocess to subproc.mg 2555 filename = 'subproc.mg' 2556 files.append_to_file(filename, 2557 self.write_subproc, 2558 subprocdir) 2559 2560 # Generate info page 2561 gen_infohtml.make_info_html(os.path.pardir) 2562 2563 # Return to original dir 2564 os.chdir(cwd) 2565 2566 if not tot_calls: 2567 tot_calls = 0 2568 return tot_calls
2569 2570 #=========================================================================== 2571 # write_super_auto_dsig_file 2572 #===========================================================================
2573 - def write_super_auto_dsig_file(self, writer, subproc_group):
2574 """Write the auto_dsig.f file selecting between the subprocesses 2575 in subprocess group mode""" 2576 2577 replace_dict = {} 2578 2579 # Extract version number and date from VERSION file 2580 info_lines = self.get_mg5_info_lines() 2581 replace_dict['info_lines'] = info_lines 2582 2583 matrix_elements = subproc_group.get('matrix_elements') 2584 2585 # Extract process info lines 2586 process_lines = '\n'.join([self.get_process_info_lines(me) for me in \ 2587 matrix_elements]) 2588 replace_dict['process_lines'] = process_lines 2589 2590 nexternal, ninitial = matrix_elements[0].get_nexternal_ninitial() 2591 replace_dict['nexternal'] = nexternal 2592 2593 replace_dict['nsprocs'] = 2*len(matrix_elements) 2594 2595 # Generate dsig definition line 2596 dsig_def_line = "DOUBLE PRECISION " + \ 2597 ",".join(["DSIG%d" % (iproc + 1) for iproc in \ 2598 range(len(matrix_elements))]) 2599 replace_dict["dsig_def_line"] = dsig_def_line 2600 2601 # Generate dsig process lines 2602 call_dsig_proc_lines = [] 2603 for iproc in range(len(matrix_elements)): 2604 call_dsig_proc_lines.append(\ 2605 "IF(IPROC.EQ.%(num)d) DSIGPROC=DSIG%(num)d(P1,WGT,IMODE) ! %(proc)s" % \ 2606 {"num": iproc + 1, 2607 "proc": matrix_elements[iproc].get('processes')[0].base_string()}) 2608 replace_dict['call_dsig_proc_lines'] = "\n".join(call_dsig_proc_lines) 2609 2610 file = open(os.path.join(_file_path, \ 2611 'iolibs/template_files/super_auto_dsig_group_v4.inc')).read() 2612 file = file % replace_dict 2613 2614 # Write the file 2615 writer.writelines(file)
2616 2617 #=========================================================================== 2618 # write_mirrorprocs 2619 #===========================================================================
2620 - def write_mirrorprocs(self, writer, subproc_group):
2621 """Write the mirrorprocs.inc file determining which processes have 2622 IS mirror process in subprocess group mode.""" 2623 2624 lines = [] 2625 bool_dict = {True: '.true.', False: '.false.'} 2626 matrix_elements = subproc_group.get('matrix_elements') 2627 lines.append("DATA (MIRRORPROCS(I),I=1,%d)/%s/" % \ 2628 (len(matrix_elements), 2629 ",".join([bool_dict[me.get('has_mirror_process')] for \ 2630 me in matrix_elements]))) 2631 # Write the file 2632 writer.writelines(lines)
2633 2634 #=========================================================================== 2635 # write_coloramps_file 2636 #===========================================================================
2637 - def write_coloramps_file(self, writer, diagrams_for_config, maxflows, 2638 matrix_elements):
2639 """Write the coloramps.inc file for MadEvent in Subprocess group mode""" 2640 2641 # Create a map from subprocess (matrix element) to a list of 2642 # the diagrams corresponding to each config 2643 2644 lines = [] 2645 2646 subproc_to_confdiag = {} 2647 for config in diagrams_for_config: 2648 for subproc, diag in enumerate(config): 2649 try: 2650 subproc_to_confdiag[subproc].append(diag) 2651 except KeyError: 2652 subproc_to_confdiag[subproc] = [diag] 2653 2654 for subproc in sorted(subproc_to_confdiag.keys()): 2655 lines.extend(self.get_icolamp_lines(subproc_to_confdiag[subproc], 2656 matrix_elements[subproc], 2657 subproc + 1)) 2658 2659 lines.insert(0, "logical icolamp(%d,%d,%d)" % \ 2660 (maxflows, 2661 len(diagrams_for_config), 2662 len(matrix_elements))) 2663 2664 # Write the file 2665 writer.writelines(lines) 2666 2667 return True
2668 2669 #=========================================================================== 2670 # write_config_subproc_map_file 2671 #===========================================================================
2672 - def write_config_subproc_map_file(self, writer, config_subproc_map):
2673 """Write the config_subproc_map.inc file for subprocess groups""" 2674 2675 lines = [] 2676 # Output only configs that have some corresponding diagrams 2677 iconfig = 0 2678 for config in config_subproc_map: 2679 if set(config) == set([0]): 2680 continue 2681 lines.append("DATA (CONFSUB(i,%d),i=1,%d)/%s/" % \ 2682 (iconfig + 1, len(config), 2683 ",".join([str(i) for i in config]))) 2684 iconfig += 1 2685 # Write the file 2686 writer.writelines(lines) 2687 2688 return True
2689 2690 #=========================================================================== 2691 # write_configs_file 2692 #===========================================================================
2693 - def write_configs_file(self, writer, subproc_group, diagrams_for_config):
2694 """Write the configs.inc file with topology information for a 2695 subprocess group. Use the first subprocess with a diagram for each 2696 configuration.""" 2697 2698 matrix_elements = subproc_group.get('matrix_elements') 2699 model = matrix_elements[0].get('processes')[0].get('model') 2700 2701 diagrams = [] 2702 config_numbers = [] 2703 for iconfig, config in enumerate(diagrams_for_config): 2704 # Check if any diagrams correspond to this config 2705 if set(config) == set([0]): 2706 continue 2707 subproc_diags = [] 2708 for s,d in enumerate(config): 2709 if d: 2710 subproc_diags.append(matrix_elements[s].\ 2711 get('diagrams')[d-1]) 2712 else: 2713 subproc_diags.append(None) 2714 diagrams.append(subproc_diags) 2715 config_numbers.append(iconfig + 1) 2716 2717 # Extract number of external particles 2718 (nexternal, ninitial) = subproc_group.get_nexternal_ninitial() 2719 2720 return len(diagrams), \ 2721 self.write_configs_file_from_diagrams(writer, diagrams, 2722 config_numbers, 2723 nexternal, ninitial, 2724 model)
2725 2726 #=========================================================================== 2727 # write_run_configs_file 2728 #===========================================================================
2729 - def write_run_config_file(self, writer):
2730 """Write the run_configs.inc file for MadEvent""" 2731 2732 path = os.path.join(_file_path,'iolibs','template_files','madevent_run_config.inc') 2733 text = open(path).read() % {'chanperjob':'2'} 2734 writer.write(text) 2735 return True
2736 2737 2738 #=========================================================================== 2739 # write_leshouche_file 2740 #===========================================================================
2741 - def write_leshouche_file(self, writer, subproc_group):
2742 """Write the leshouche.inc file for MG4""" 2743 2744 all_lines = [] 2745 2746 for iproc, matrix_element in \ 2747 enumerate(subproc_group.get('matrix_elements')): 2748 all_lines.extend(self.get_leshouche_lines(matrix_element, 2749 iproc)) 2750 2751 # Write the file 2752 writer.writelines(all_lines) 2753 2754 return True
2755 2756 #=============================================================================== 2757 # UFO_model_to_mg4 2758 #=============================================================================== 2759 python_to_fortran = lambda x: parsers.UFOExpressionParserFortran().parse(x) 2760
2761 -class UFO_model_to_mg4(object):
2762 """ A converter of the UFO-MG5 Model to the MG4 format """ 2763
2764 - def __init__(self, model, output_path):
2765 """ initialization of the objects """ 2766 2767 self.model = model 2768 self.model_name = model['name'] 2769 self.dir_path = output_path 2770 2771 self.coups_dep = [] # (name, expression, type) 2772 self.coups_indep = [] # (name, expression, type) 2773 self.params_dep = [] # (name, expression, type) 2774 self.params_indep = [] # (name, expression, type) 2775 self.params_ext = [] # external parameter 2776 self.p_to_f = parsers.UFOExpressionParserFortran()
2777
2779 """modify the parameter if some of them are identical up to the case""" 2780 2781 lower_dict={} 2782 duplicate = set() 2783 keys = self.model['parameters'].keys() 2784 for key in keys: 2785 for param in self.model['parameters'][key]: 2786 lower_name = param.name.lower() 2787 try: 2788 lower_dict[lower_name].append(param) 2789 except KeyError: 2790 lower_dict[lower_name] = [param] 2791 else: 2792 duplicate.add(lower_name) 2793 2794 if not duplicate: 2795 return 2796 2797 re_expr = r'''\b(%s)\b''' 2798 to_change = [] 2799 change={} 2800 for value in duplicate: 2801 for i, var in enumerate(lower_dict[value][1:]): 2802 to_change.append(var.name) 2803 change[var.name] = '%s__%s' %( var.name.lower(), i+2) 2804 var.name = '%s__%s' %( var.name.lower(), i+2) 2805 2806 replace = lambda match_pattern: change[match_pattern.groups()[0]] 2807 2808 rep_pattern = re.compile(re_expr % '|'.join(to_change)) 2809 2810 # change parameters 2811 for key in keys: 2812 if key == ('external',): 2813 continue 2814 for param in self.model['parameters'][key]: 2815 param.expr = rep_pattern.sub(replace, param.expr) 2816 2817 # change couplings 2818 for key in self.model['couplings'].keys(): 2819 for coup in self.model['couplings'][key]: 2820 coup.expr = rep_pattern.sub(replace, coup.expr)
2821 2822 2823
2824 - def refactorize(self, wanted_couplings = []):
2825 """modify the couplings to fit with MG4 convention """ 2826 2827 # Keep only separation in alphaS 2828 keys = self.model['parameters'].keys() 2829 keys.sort(key=len) 2830 for key in keys: 2831 if key == ('external',): 2832 self.params_ext += self.model['parameters'][key] 2833 elif 'aS' in key: 2834 self.params_dep += self.model['parameters'][key] 2835 else: 2836 self.params_indep += self.model['parameters'][key] 2837 # same for couplings 2838 keys = self.model['couplings'].keys() 2839 keys.sort(key=len) 2840 for key, coup_list in self.model['couplings'].items(): 2841 if 'aS' in key: 2842 self.coups_dep += [c for c in coup_list if 2843 (not wanted_couplings or c.name in \ 2844 wanted_couplings)] 2845 else: 2846 self.coups_indep += [c for c in coup_list if 2847 (not wanted_couplings or c.name in \ 2848 wanted_couplings)] 2849 2850 # MG4 use G and not aS as it basic object for alphas related computation 2851 #Pass G in the independant list 2852 if 'G' in self.params_dep: 2853 index = self.params_dep.index('G') 2854 G = self.params_dep.pop(index)
2855 # G.expr = '2*cmath.sqrt(as*pi)' 2856 # self.params_indep.insert(0, self.params_dep.pop(index)) 2857 # No need to add it if not defined 2858 2859
2860 - def build(self, wanted_couplings = [], full=True):
2861 """modify the couplings to fit with MG4 convention and creates all the 2862 different files""" 2863 2864 self.pass_parameter_to_case_insensitive() 2865 self.refactorize(wanted_couplings) 2866 2867 # write the files 2868 if full: 2869 self.write_all()
2870
2871 - def open(self, name, comment='c', format='default'):
2872 """ Open the file name in the correct directory and with a valid 2873 header.""" 2874 2875 file_path = os.path.join(self.dir_path, name) 2876 2877 if format == 'fortran': 2878 fsock = writers.FortranWriter(file_path, 'w') 2879 else: 2880 fsock = open(file_path, 'w') 2881 2882 file.writelines(fsock, comment * 77 + '\n') 2883 file.writelines(fsock,'%(comment)s written by the UFO converter\n' % \ 2884 {'comment': comment + (6 - len(comment)) * ' '}) 2885 file.writelines(fsock, comment * 77 + '\n\n') 2886 return fsock
2887 2888
2889 - def write_all(self):
2890 """ write all the files """ 2891 #write the part related to the external parameter 2892 self.create_ident_card() 2893 self.create_param_read() 2894 2895 #write the definition of the parameter 2896 self.create_input() 2897 self.create_intparam_def() 2898 2899 2900 # definition of the coupling. 2901 self.create_coupl_inc() 2902 self.create_write_couplings() 2903 self.create_couplings() 2904 2905 # the makefile 2906 self.create_makeinc() 2907 self.create_param_write() 2908 2909 # The param_card.dat 2910 self.create_param_card() 2911 2912 2913 # All the standard files 2914 self.copy_standard_file()
2915 2916 ############################################################################ 2917 ## ROUTINE CREATING THE FILES ############################################ 2918 ############################################################################ 2919
2920 - def copy_standard_file(self):
2921 """Copy the standard files for the fortran model.""" 2922 2923 2924 #copy the library files 2925 file_to_link = ['formats.inc', 'lha_read.f', 'makefile','printout.f', \ 2926 'rw_para.f', 'testprog.f', 'rw_para.f'] 2927 2928 for filename in file_to_link: 2929 cp( MG5DIR + '/models/template_files/fortran/' + filename, self.dir_path)
2930
2931 - def create_coupl_inc(self):
2932 """ write coupling.inc """ 2933 2934 fsock = self.open('coupl.inc', format='fortran') 2935 2936 # Write header 2937 header = """double precision G 2938 common/strong/ G 2939 2940 double complex gal(2) 2941 common/weak/ gal 2942 2943 """ 2944 fsock.writelines(header) 2945 2946 # Write the Mass definition/ common block 2947 masses = set() 2948 widths = set() 2949 for particle in self.model.get('particles'): 2950 #find masses 2951 one_mass = particle.get('mass') 2952 if one_mass.lower() != 'zero': 2953 masses.add(one_mass) 2954 2955 # find width 2956 one_width = particle.get('width') 2957 if one_width.lower() != 'zero': 2958 widths.add(one_width) 2959 2960 2961 fsock.writelines('double precision '+','.join(masses)+'\n') 2962 fsock.writelines('common/masses/ '+','.join(masses)+'\n\n') 2963 if widths: 2964 fsock.writelines('double precision '+','.join(widths)+'\n') 2965 fsock.writelines('common/widths/ '+','.join(widths)+'\n\n') 2966 2967 # Write the Couplings 2968 coupling_list = [coupl.name for coupl in self.coups_dep + self.coups_indep] 2969 fsock.writelines('double complex '+', '.join(coupling_list)+'\n') 2970 fsock.writelines('common/couplings/ '+', '.join(coupling_list)+'\n')
2971
2972 - def create_write_couplings(self):
2973 """ write the file coupl_write.inc """ 2974 2975 fsock = self.open('coupl_write.inc', format='fortran') 2976 2977 fsock.writelines("""write(*,*) ' Couplings of %s' 2978 write(*,*) ' ---------------------------------' 2979 write(*,*) ' '""" % self.model_name) 2980 def format(coupl): 2981 return 'write(*,2) \'%(name)s = \', %(name)s' % {'name': coupl.name}
2982 2983 # Write the Couplings 2984 lines = [format(coupl) for coupl in self.coups_dep + self.coups_indep] 2985 fsock.writelines('\n'.join(lines)) 2986 2987
2988 - def create_input(self):
2989 """create input.inc containing the definition of the parameters""" 2990 2991 fsock = self.open('input.inc', format='fortran') 2992 2993 #find mass/ width since they are already define 2994 already_def = set() 2995 for particle in self.model.get('particles'): 2996 already_def.add(particle.get('mass').lower()) 2997 already_def.add(particle.get('width').lower()) 2998 2999 is_valid = lambda name: name !='G' and name.lower() not in already_def 3000 3001 real_parameters = [param.name for param in self.params_dep + 3002 self.params_indep if param.type == 'real' 3003 and is_valid(param.name)] 3004 3005 real_parameters += [param.name for param in self.params_ext 3006 if param.type == 'real'and 3007 is_valid(param.name)] 3008 3009 fsock.writelines('double precision '+','.join(real_parameters)+'\n') 3010 fsock.writelines('common/params_R/ '+','.join(real_parameters)+'\n\n') 3011 3012 complex_parameters = [param.name for param in self.params_dep + 3013 self.params_indep if param.type == 'complex'] 3014 3015 3016 fsock.writelines('double complex '+','.join(complex_parameters)+'\n') 3017 fsock.writelines('common/params_C/ '+','.join(complex_parameters)+'\n\n')
3018 3019 3020
3021 - def create_intparam_def(self):
3022 """ create intparam_definition.inc """ 3023 3024 fsock = self.open('intparam_definition.inc', format='fortran') 3025 3026 fsock.write_comments(\ 3027 "Parameters that should not be recomputed event by event.\n") 3028 fsock.writelines("if(readlha) then\n") 3029 fsock.writelines("G = 2 * DSQRT(AS*PI) ! for the first init\n") 3030 for param in self.params_indep: 3031 if param.name == 'ZERO': 3032 continue 3033 fsock.writelines("%s = %s\n" % (param.name, 3034 self.p_to_f.parse(param.expr))) 3035 3036 fsock.writelines('endif') 3037 3038 fsock.write_comments('\nParameters that should be recomputed at an event by even basis.\n') 3039 fsock.writelines("aS = G**2/4/pi\n") 3040 for param in self.params_dep: 3041 fsock.writelines("%s = %s\n" % (param.name, 3042 self.p_to_f.parse(param.expr))) 3043 3044 fsock.write_comments("\nDefinition of the EW coupling used in the write out of aqed\n") 3045 fsock.writelines(""" gal(1) = 1d0 3046 gal(2) = 1d0 3047 """)
3048 3049
3050 - def create_couplings(self):
3051 """ create couplings.f and all couplingsX.f """ 3052 3053 nb_def_by_file = 25 3054 3055 self.create_couplings_main(nb_def_by_file) 3056 nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 3057 nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 3058 3059 for i in range(nb_coup_indep): 3060 data = self.coups_indep[nb_def_by_file * i: 3061 min(len(self.coups_indep), nb_def_by_file * (i+1))] 3062 self.create_couplings_part(i + 1, data) 3063 3064 for i in range(nb_coup_dep): 3065 data = self.coups_dep[nb_def_by_file * i: 3066 min(len(self.coups_dep), nb_def_by_file * (i+1))] 3067 self.create_couplings_part( i + 1 + nb_coup_indep , data)
3068 3069
3070 - def create_couplings_main(self, nb_def_by_file=25):
3071 """ create couplings.f """ 3072 3073 fsock = self.open('couplings.f', format='fortran') 3074 3075 fsock.writelines("""subroutine coup() 3076 3077 implicit none 3078 double precision PI 3079 logical READLHA 3080 parameter (PI=3.141592653589793d0) 3081 3082 include \'input.inc\' 3083 include \'coupl.inc\' 3084 READLHA = .true. 3085 include \'intparam_definition.inc\'\n\n 3086 """) 3087 3088 nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 3089 nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 3090 3091 fsock.writelines('\n'.join(\ 3092 ['call coup%s()' % (i + 1) for i in range(nb_coup_indep)])) 3093 3094 fsock.write_comments('\ncouplings needed to be evaluated points by points\n') 3095 3096 fsock.writelines('\n'.join(\ 3097 ['call coup%s()' % (nb_coup_indep + i + 1) \ 3098 for i in range(nb_coup_dep)])) 3099 fsock.writelines('''\n return \n end\n''') 3100 3101 fsock.writelines("""subroutine update_as_param() 3102 3103 implicit none 3104 double precision PI 3105 logical READLHA 3106 parameter (PI=3.141592653589793d0) 3107 3108 include \'input.inc\' 3109 include \'coupl.inc\' 3110 READLHA = .false. 3111 include \'intparam_definition.inc\'\n\n 3112 """) 3113 3114 nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 3115 nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 3116 3117 fsock.write_comments('\ncouplings needed to be evaluated points by points\n') 3118 3119 fsock.writelines('\n'.join(\ 3120 ['call coup%s()' % (nb_coup_indep + i + 1) \ 3121 for i in range(nb_coup_dep)])) 3122 fsock.writelines('''\n return \n end\n''')
3123
3124 - def create_couplings_part(self, nb_file, data):
3125 """ create couplings[nb_file].f containing information coming from data 3126 """ 3127 3128 fsock = self.open('couplings%s.f' % nb_file, format='fortran') 3129 fsock.writelines("""subroutine coup%s() 3130 3131 implicit none 3132 double precision PI 3133 parameter (PI=3.141592653589793d0) 3134 3135 3136 include 'input.inc' 3137 include 'coupl.inc' 3138 """ % nb_file) 3139 3140 for coupling in data: 3141 fsock.writelines('%s = %s' % (coupling.name, 3142 self.p_to_f.parse(coupling.expr))) 3143 fsock.writelines('end')
3144 3145
3146 - def create_makeinc(self):
3147 """create makeinc.inc containing the file to compile """ 3148 3149 fsock = self.open('makeinc.inc', comment='#') 3150 text = 'MODEL = couplings.o lha_read.o printout.o rw_para.o ' 3151 3152 nb_coup_indep = 1 + len(self.coups_dep) // 25 3153 nb_coup_dep = 1 + len(self.coups_indep) // 25 3154 text += ' '.join(['couplings%s.o' % (i+1) \ 3155 for i in range(nb_coup_dep + nb_coup_indep) ]) 3156 fsock.writelines(text)
3157
3158 - def create_param_write(self):
3159 """ create param_write """ 3160 3161 fsock = self.open('param_write.inc', format='fortran') 3162 3163 fsock.writelines("""write(*,*) ' External Params' 3164 write(*,*) ' ---------------------------------' 3165 write(*,*) ' '""") 3166 def format(name): 3167 return 'write(*,*) \'%(name)s = \', %(name)s' % {'name': name}
3168 3169 # Write the external parameter 3170 lines = [format(param.name) for param in self.params_ext] 3171 fsock.writelines('\n'.join(lines)) 3172 3173 fsock.writelines("""write(*,*) ' Internal Params' 3174 write(*,*) ' ---------------------------------' 3175 write(*,*) ' '""") 3176 lines = [format(data.name) for data in self.params_indep 3177 if data.name != 'ZERO'] 3178 fsock.writelines('\n'.join(lines)) 3179 fsock.writelines("""write(*,*) ' Internal Params evaluated point by point' 3180 write(*,*) ' ----------------------------------------' 3181 write(*,*) ' '""") 3182 lines = [format(data.name) for data in self.params_dep] 3183 3184 fsock.writelines('\n'.join(lines)) 3185 3186 3187
3188 - def create_ident_card(self):
3189 """ create the ident_card.dat """ 3190 3191 def format(parameter): 3192 """return the line for the ident_card corresponding to this parameter""" 3193 colum = [parameter.lhablock.lower()] + \ 3194 [str(value) for value in parameter.lhacode] + \ 3195 [parameter.name] 3196 return ' '.join(colum)+'\n'
3197 3198 fsock = self.open('ident_card.dat') 3199 3200 external_param = [format(param) for param in self.params_ext] 3201 fsock.writelines('\n'.join(external_param)) 3202 3203
3204 - def create_param_read(self):
3205 """create param_read""" 3206 3207 def format_line(parameter): 3208 """return the line for the ident_card corresponding to this 3209 parameter""" 3210 template = \ 3211 """ call LHA_get_real(npara,param,value,'%(name)s',%(name)s,%(value)s)""" \ 3212 % {'name': parameter.name, 3213 'value': self.p_to_f.parse(str(parameter.value.real))} 3214 3215 return template 3216 3217 fsock = self.open('param_read.inc', format='fortran') 3218 res_strings = [format_line(param) \ 3219 for param in self.params_ext] 3220 3221 # Correct width sign for Majorana particles (where the width 3222 # and mass need to have the same sign) 3223 for particle in self.model.get('particles'): 3224 if particle.is_fermion() and particle.get('self_antipart') and \ 3225 particle.get('width').lower() != 'zero': 3226 3227 res_strings.append('%(width)s = sign(%(width)s,%(mass)s)' % \ 3228 {'width': particle.get('width'), 'mass': particle.get('mass')}) 3229 3230 3231 fsock.writelines('\n'.join(res_strings)) 3232
3233 - def create_param_card(self):
3234 """ create the param_card.dat """ 3235 3236 out_path = os.path.join(self.dir_path, 'param_card.dat') 3237 param_writer.ParamCardWriter(self.model, out_path) 3238 out_path2 = None 3239 if hasattr(self.model, 'rule_card'): 3240 out_path2 = os.path.join(self.dir_path, 'param_card_rule.dat') 3241 self.model.rule_card.write_file(out_path2) 3242 3243 # IF MSSM convert the card to SLAH1 3244 if self.model_name == 'mssm' or self.model_name.startswith('mssm-'): 3245 import models.check_param_card as translator 3246 3247 # Check the format of the param_card for Pythia and make it correct 3248 if out_path2: 3249 translator.make_valid_param_card(out_path, out_path2) 3250 translator.convert_to_slha1(out_path)
3251