################################################################################ # # Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors # # This file is a part of the MadGraph5_aMC@NLO project, an application which # automatically generates Feynman diagrams and matrix elements for arbitrary # high-energy processes in the Standard Model and beyond. # # It is subject to the MadGraph5_aMC@NLO license which should accompany this # distribution. # # For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch # ################################################################################ """Methods and classes to export matrix elements to fks format.""" from distutils import dir_util import glob import logging import os import re import shutil import subprocess import string import copy import platform import madgraph.core.color_algebra as color import madgraph.core.helas_objects as helas_objects import madgraph.core.base_objects as base_objects import madgraph.fks.fks_helas_objects as fks_helas_objects import madgraph.fks.fks_base as fks import madgraph.fks.fks_common as fks_common import madgraph.iolibs.drawing_eps as draw import madgraph.iolibs.gen_infohtml as gen_infohtml import madgraph.iolibs.files as files import madgraph.various.misc as misc import madgraph.iolibs.file_writers as writers import madgraph.iolibs.template_files as template_files import madgraph.iolibs.ufo_expression_parsers as parsers import madgraph.iolibs.export_v4 as export_v4 import madgraph.loop.loop_exporters as loop_exporters import madgraph.various.q_polynomial as q_polynomial import madgraph.various.banner as banner_mod import aloha.create_aloha as create_aloha import models.write_param_card as write_param_card import models.check_param_card as check_param_card from madgraph import MadGraph5Error, MG5DIR, InvalidCmd from madgraph.iolibs.files import cp, ln, mv pjoin = os.path.join _file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/' logger = logging.getLogger('madgraph.export_fks') def make_jpeg_async(args): Pdir = args[0] old_pos = args[1] dir_path = args[2] devnull = os.open(os.devnull, os.O_RDWR) os.chdir(Pdir) subprocess.call([os.path.join(old_pos, dir_path, 'bin', 'internal', 'gen_jpeg-pl')], stdout = devnull) os.chdir(os.path.pardir) #================================================================================= # Class for used of the (non-optimized) Loop process #================================================================================= class ProcessExporterFortranFKS(loop_exporters.LoopProcessExporterFortranSA): """Class to take care of exporting a set of matrix elements to Fortran (v4) format.""" #=============================================================================== # copy the Template in a new directory. #=============================================================================== def copy_fkstemplate(self): """create the directory run_name as a copy of the MadEvent Template, and clean the directory For now it is just the same as copy_v4template, but it will be modified """ mgme_dir = self.mgme_dir dir_path = self.dir_path clean =self.opt['clean'] #First copy the full template tree if dir_path doesn't exit if not os.path.isdir(dir_path): if not mgme_dir: raise MadGraph5Error, \ "No valid MG_ME path given for MG4 run directory creation." logger.info('initialize a new directory: %s' % \ os.path.basename(dir_path)) shutil.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True) # distutils.dir_util.copy_tree since dir_path already exists dir_util.copy_tree(pjoin(self.mgme_dir, 'Template', 'Common'),dir_path) # Copy plot_card for card in ['plot_card']: if os.path.isfile(pjoin(self.dir_path, 'Cards',card + '.dat')): try: shutil.copy(pjoin(self.dir_path, 'Cards', card + '.dat'), pjoin(self.dir_path, 'Cards', card + '_default.dat')) except IOError: logger.warning("Failed to move " + card + ".dat to default") elif not os.path.isfile(os.path.join(dir_path, 'TemplateVersion.txt')): if not mgme_dir: raise MadGraph5Error, \ "No valid MG_ME path given for MG4 run directory creation." try: shutil.copy(os.path.join(mgme_dir, 'MGMEVersion.txt'), dir_path) except IOError: MG5_version = misc.get_pkg_info() open(os.path.join(dir_path, 'MGMEVersion.txt'), 'w').write( \ "5." + MG5_version['version']) #Ensure that the Template is clean if clean: logger.info('remove old information in %s' % os.path.basename(dir_path)) if os.environ.has_key('MADGRAPH_BASE'): subprocess.call([os.path.join('bin', 'internal', 'clean_template'), '--web'],cwd=dir_path) else: try: subprocess.call([os.path.join('bin', 'internal', 'clean_template')], \ cwd=dir_path) except Exception, why: raise MadGraph5Error('Failed to clean correctly %s: \n %s' \ % (os.path.basename(dir_path),why)) #Write version info MG_version = misc.get_pkg_info() open(os.path.join(dir_path, 'SubProcesses', 'MGVersion.txt'), 'w').write( MG_version['version']) # We must link the CutTools to the Library folder of the active Template self.link_CutTools(dir_path) link_tir_libs=[] tir_libs=[] os.remove(os.path.join(self.dir_path,'SubProcesses','makefile_loop.inc')) dirpath = os.path.join(self.dir_path, 'SubProcesses') filename = pjoin(self.dir_path, 'SubProcesses','makefile_loop') calls = self.write_makefile_TIR(writers.MakefileWriter(filename), link_tir_libs,tir_libs) os.remove(os.path.join(self.dir_path,'Source','make_opts.inc')) filename = pjoin(self.dir_path, 'Source','make_opts') calls = self.write_make_opts(writers.MakefileWriter(filename), link_tir_libs,tir_libs) # Duplicate run_card and FO_analyse_card for card in ['FO_analyse_card', 'shower_card']: try: shutil.copy(pjoin(self.dir_path, 'Cards', card + '.dat'), pjoin(self.dir_path, 'Cards', card + '_default.dat')) except IOError: logger.warning("Failed to copy " + card + ".dat to default") cwd = os.getcwd() dirpath = os.path.join(self.dir_path, 'SubProcesses') try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 # We add here the user-friendly MadLoop option setter. cpfiles= ["SubProcesses/MadLoopParamReader.f", "Cards/MadLoopParams.dat", "SubProcesses/MadLoopParams.inc"] for file in cpfiles: shutil.copy(os.path.join(self.loop_dir,'StandAlone/', file), os.path.join(self.dir_path, file)) shutil.copy(pjoin(self.dir_path, 'Cards','MadLoopParams.dat'), pjoin(self.dir_path, 'Cards','MadLoopParams_default.dat')) if os.path.exists(pjoin(self.dir_path, 'Cards', 'MadLoopParams.dat')): self.MadLoopparam = banner_mod.MadLoopParam(pjoin(self.dir_path, 'Cards', 'MadLoopParams.dat')) # write the output file self.MadLoopparam.write(pjoin(self.dir_path,"SubProcesses", "MadLoopParams.dat")) # We need minimal editing of MadLoopCommons.f MadLoopCommon = open(os.path.join(self.loop_dir,'StandAlone', "SubProcesses","MadLoopCommons.inc")).read() writer = writers.FortranWriter(os.path.join(self.dir_path, "SubProcesses","MadLoopCommons.f")) writer.writelines(MadLoopCommon%{ 'print_banner_commands':self.MadLoop_banner}, context={'collier_available':False}) writer.close() # Write the cts_mpc.h and cts_mprec.h files imported from CutTools self.write_mp_files(writers.FortranWriter('cts_mprec.h'),\ writers.FortranWriter('cts_mpc.h')) # Finally make sure to turn off MC over Hel for the default mode. FKS_card_path = pjoin(self.dir_path,'Cards','FKS_params.dat') FKS_card_file = open(FKS_card_path,'r') FKS_card = FKS_card_file.read() FKS_card_file.close() FKS_card = re.sub(r"#NHelForMCoverHels\n-?\d+", "#NHelForMCoverHels\n-1", FKS_card) FKS_card_file = open(FKS_card_path,'w') FKS_card_file.write(FKS_card) FKS_card_file.close() # Return to original PWD os.chdir(cwd) # Copy the different python files in the Template self.copy_python_files() # We need to create the correct open_data for the pdf self.write_pdf_opendata() # I put it here not in optimized one, because I want to use the same makefile_loop.inc # Also, we overload this function (i.e. it is already defined in # LoopProcessExporterFortranSA) because the path of the template makefile # is different. def write_makefile_TIR(self, writer, link_tir_libs,tir_libs,tir_include=[]): """ Create the file makefile_loop which links to the TIR libraries.""" file = open(os.path.join(self.mgme_dir,'Template','NLO', 'SubProcesses','makefile_loop.inc')).read() replace_dict={} replace_dict['link_tir_libs']=' '.join(link_tir_libs) replace_dict['tir_libs']=' '.join(tir_libs) replace_dict['dotf']='%.f' replace_dict['doto']='%.o' replace_dict['tir_include']=' '.join(tir_include) file=file%replace_dict if writer: writer.writelines(file) else: return file # I put it here not in optimized one, because I want to use the same make_opts.inc def write_make_opts(self, writer, link_tir_libs,tir_libs): """ Create the file make_opts which links to the TIR libraries.""" file = open(os.path.join(self.mgme_dir,'Template','NLO', 'Source','make_opts.inc')).read() replace_dict={} replace_dict['link_tir_libs']=' '.join(link_tir_libs) replace_dict['tir_libs']=' '.join(tir_libs) replace_dict['dotf']='%.f' replace_dict['doto']='%.o' file=file%replace_dict if writer: writer.writelines(file) else: return file #=========================================================================== # copy_python_files #=========================================================================== def copy_python_files(self): """copy python files required for the Template""" files_to_copy = [ \ pjoin('interface','amcatnlo_run_interface.py'), pjoin('interface','extended_cmd.py'), pjoin('interface','common_run_interface.py'), pjoin('interface','coloring_logging.py'), pjoin('various','misc.py'), pjoin('various','shower_card.py'), pjoin('various','FO_analyse_card.py'), pjoin('various','histograms.py'), pjoin('various','banner.py'), pjoin('various','cluster.py'), pjoin('various','systematics.py'), pjoin('various','lhe_parser.py'), pjoin('madevent','sum_html.py'), pjoin('madevent','gen_crossxhtml.py'), pjoin('iolibs','files.py'), pjoin('iolibs','save_load_object.py'), pjoin('iolibs','file_writers.py'), pjoin('..','models','check_param_card.py'), pjoin('__init__.py') ] cp(_file_path+'/interface/.mg5_logging.conf', self.dir_path+'/bin/internal/me5_logging.conf') for cp_file in files_to_copy: cp(pjoin(_file_path,cp_file), pjoin(self.dir_path,'bin','internal',os.path.basename(cp_file))) def convert_model(self, model, wanted_lorentz = [], wanted_couplings = []): super(ProcessExporterFortranFKS,self).convert_model(model, wanted_lorentz, wanted_couplings) IGNORE_PATTERNS = ('*.pyc','*.dat','*.py~') try: shutil.rmtree(pjoin(self.dir_path,'bin','internal','ufomodel')) except OSError as error: pass model_path = model.get('modelpath') shutil.copytree(model_path, pjoin(self.dir_path,'bin','internal','ufomodel'), ignore=shutil.ignore_patterns(*IGNORE_PATTERNS)) if hasattr(model, 'restrict_card'): out_path = pjoin(self.dir_path, 'bin', 'internal','ufomodel', 'restrict_default.dat') if isinstance(model.restrict_card, check_param_card.ParamCard): model.restrict_card.write(out_path) else: files.cp(model.restrict_card, out_path) #=========================================================================== # write_maxparticles_file #=========================================================================== def write_maxparticles_file(self, writer, maxparticles): """Write the maxparticles.inc file for MadEvent""" lines = "integer max_particles, max_branch\n" lines += "parameter (max_particles=%d) \n" % maxparticles lines += "parameter (max_branch=max_particles-1)" # Write the file writer.writelines(lines) return True #=========================================================================== # write_maxconfigs_file #=========================================================================== def write_maxconfigs_file(self, writer, maxconfigs): """Write the maxconfigs.inc file for MadEvent""" lines = "integer lmaxconfigs\n" lines += "parameter (lmaxconfigs=%d)" % maxconfigs # Write the file writer.writelines(lines) return True #=============================================================================== # write a procdef_mg5 (an equivalent of the MG4 proc_card.dat) #=============================================================================== def write_procdef_mg5(self, file_pos, modelname, process_str): """ write an equivalent of the MG4 proc_card in order that all the Madevent Perl script of MadEvent4 are still working properly for pure MG5 run.""" proc_card_template = template_files.mg4_proc_card.mg4_template process_template = template_files.mg4_proc_card.process_template process_text = '' coupling = '' new_process_content = [] # First find the coupling and suppress the coupling from process_str #But first ensure that coupling are define whithout spaces: process_str = process_str.replace(' =', '=') process_str = process_str.replace('= ', '=') process_str = process_str.replace(',',' , ') #now loop on the element and treat all the coupling for info in process_str.split(): if '=' in info: coupling += info + '\n' else: new_process_content.append(info) # Recombine the process_str (which is the input process_str without coupling #info) process_str = ' '.join(new_process_content) #format the SubProcess process_text += process_template.substitute({'process': process_str, \ 'coupling': coupling}) text = proc_card_template.substitute({'process': process_text, 'model': modelname, 'multiparticle':''}) ff = open(file_pos, 'w') ff.write(text) ff.close() #=============================================================================== # write a initial states map, useful for the fast PDF NLO interface #=============================================================================== def write_init_map(self, file_pos, initial_states): """ Write an initial state process map. Each possible PDF combination gets an unique identifier.""" text='' for i,e in enumerate(initial_states): text=text+str(i+1)+' '+str(len(e)) for t in e: text=text+' ' try: for p in t: text=text+' '+str(p) except TypeError: text=text+' '+str(t) text=text+'\n' ff = open(file_pos, 'w') ff.write(text) ff.close() def get_ME_identifier(self, matrix_element, *args, **opts): """ A function returning a string uniquely identifying the matrix element given in argument so that it can be used as a prefix to all MadLoop5 subroutines and common blocks related to it. This allows to compile several processes into one library as requested by the BLHA (Binoth LesHouches Accord) guidelines. The MadFKS design necessitates that there is no process prefix.""" return '' #=============================================================================== # write_coef_specs #=============================================================================== def write_coef_specs_file(self, virt_me_list): """writes the coef_specs.inc in the DHELAS folder. Should not be called in the non-optimized mode""" raise fks_common.FKSProcessError(), \ "write_coef_specs should be called only in the loop-optimized mode" #=============================================================================== # generate_directories_fks #=============================================================================== def generate_directories_fks(self, matrix_element, fortran_model, me_number, me_ntot, path=os.getcwd(),OLP='MadLoop'): """Generate the Pxxxxx_i directories for a subprocess in MadFKS, including the necessary matrix.f and various helper files""" proc = matrix_element.born_matrix_element['processes'][0] if not self.model: self.model = matrix_element.get('processes')[0].get('model') cwd = os.getcwd() try: os.chdir(path) except OSError, error: error_msg = "The directory %s should exist in order to be able " % path + \ "to \"export\" in it. If you see this error message by " + \ "typing the command \"export\" please consider to use " + \ "instead the command \"output\". " raise MadGraph5Error, error_msg calls = 0 self.fksdirs = [] #first make and cd the direcrory corresponding to the born process: borndir = "P%s" % \ (matrix_element.get('processes')[0].shell_string()) os.mkdir(borndir) os.chdir(borndir) logger.info('Writing files in %s (%d / %d)' % (borndir, me_number + 1, me_ntot)) ## write the files corresponding to the born process in the P* directory self.generate_born_fks_files(matrix_element, fortran_model, me_number, path) # With NJET you want to generate the order file per subprocess and most # likely also generate it for each subproc. if OLP=='NJET': filename = 'OLE_order.lh' self.write_lh_order(filename, [matrix_element.born_matrix_element.get('processes')[0]], OLP) if matrix_element.virt_matrix_element: calls += self.generate_virt_directory( \ matrix_element.virt_matrix_element, \ fortran_model, \ os.path.join(path, borndir)) #write the infortions for the different real emission processes self.write_real_matrix_elements(matrix_element, fortran_model) self.write_pdf_calls(matrix_element, fortran_model) filename = 'nFKSconfigs.inc' self.write_nfksconfigs_file(writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'iproc.dat' self.write_iproc_file(writers.FortranWriter(filename), me_number) filename = 'fks_info.inc' self.write_fks_info_file(writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'leshouche_info.dat' nfksconfs,maxproc,maxflow,nexternal=\ self.write_leshouche_info_file(filename,matrix_element) # if no corrections are generated ([LOonly] mode), get # these variables from the born if nfksconfs == maxproc == maxflow == 0: nfksconfs = 1 (dummylines, maxproc, maxflow) = self.get_leshouche_lines( matrix_element.born_matrix_element, 1) filename = 'leshouche_decl.inc' self.write_leshouche_info_declarations( writers.FortranWriter(filename), nfksconfs,maxproc,maxflow,nexternal, fortran_model) filename = 'genps.inc' ngraphs = matrix_element.born_matrix_element.get_number_of_amplitudes() ncolor = max(1,len(matrix_element.born_matrix_element.get('color_basis'))) self.write_genps(writers.FortranWriter(filename),maxproc,ngraphs,\ ncolor,maxflow,fortran_model) filename = 'configs_and_props_info.dat' nconfigs,max_leg_number,nfksconfs=self.write_configs_and_props_info_file( filename, matrix_element) filename = 'configs_and_props_decl.inc' self.write_configs_and_props_info_declarations( writers.FortranWriter(filename), nconfigs,max_leg_number,nfksconfs, fortran_model) filename = 'real_from_born_configs.inc' self.write_real_from_born_configs( writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'ngraphs.inc' self.write_ngraphs_file(writers.FortranWriter(filename), nconfigs) #write the wrappers filename = 'real_me_chooser.f' self.write_real_me_wrapper(writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'parton_lum_chooser.f' self.write_pdf_wrapper(writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'get_color.f' self.write_colors_file(writers.FortranWriter(filename), matrix_element) filename = 'nexternal.inc' (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() self.write_nexternal_file(writers.FortranWriter(filename), nexternal, ninitial) self.proc_characteristic['ninitial'] = ninitial self.proc_characteristic['nexternal'] = max(self.proc_characteristic['nexternal'], nexternal) filename = 'pmass.inc' try: self.write_pmass_file(writers.FortranWriter(filename), matrix_element.real_processes[0].matrix_element) except IndexError: self.write_pmass_file(writers.FortranWriter(filename), matrix_element.born_matrix_element) #draw the diagrams self.draw_feynman_diagrams(matrix_element) linkfiles = ['BinothLHADummy.f', 'check_poles.f', 'MCmasses_HERWIG6.inc', 'MCmasses_HERWIGPP.inc', 'MCmasses_PYTHIA6Q.inc', 'MCmasses_PYTHIA6PT.inc', 'MCmasses_PYTHIA8.inc', 'add_write_info.f', 'coupl.inc', 'cuts.f', 'FKS_params.dat', 'initial_states_map.dat', 'OLE_order.olc', 'FKSParams.inc', 'FKSParamReader.f', 'cuts.inc', 'unlops.inc', 'pythia_unlops.f', 'driver_mintMC.f', 'driver_mintFO.f', 'appl_interface.cc', 'appl_interface_dummy.f', 'appl_common.inc', 'reweight_appl.inc', 'fastjetfortran_madfks_core.cc', 'fastjetfortran_madfks_full.cc', 'fjcore.cc', 'fastjet_wrapper.f', 'fjcore.hh', 'fks_Sij.f', 'fks_powers.inc', 'fks_singular.f', 'veto_xsec.f', 'veto_xsec.inc', 'c_weight.inc', 'fks_inc_chooser.f', 'leshouche_inc_chooser.f', 'configs_and_props_inc_chooser.f', 'genps_fks.f', 'boostwdir2.f', 'madfks_mcatnlo.inc', 'open_output_files.f', 'open_output_files_dummy.f', 'HwU_dummy.f', 'madfks_plot.f', 'analysis_dummy.f', 'analysis_lhe.f', 'mint-integrator2.f', 'MC_integer.f', 'mint.inc', 'montecarlocounter.f', 'q_es.inc', 'recluster.cc', 'Boosts.h', 'reweight.inc', 'reweight0.inc', 'reweight1.inc', 'reweightNLO.inc', 'reweight_all.inc', 'reweight_xsec.f', 'reweight_xsec_events.f', 'reweight_xsec_events_pdf_dummy.f', 'iproc_map.f', 'run.inc', 'run_card.inc', 'setcuts.f', 'setscales.f', 'symmetry_fks_test_MC.f', 'symmetry_fks_test_ME.f', 'symmetry_fks_v3.f', 'trapfpe.c', 'vegas2.for', 'write_ajob.f', 'handling_lhe_events.f', 'write_event.f', 'fill_MC_mshell.f', 'maxparticles.inc', 'message.inc', 'initcluster.f', 'cluster.inc', 'cluster.f', 'reweight.f', 'randinit', 'sudakov.inc', 'maxconfigs.inc', 'timing_variables.inc'] for file in linkfiles: ln('../' + file , '.') os.system("ln -s ../../Cards/param_card.dat .") #copy the makefile os.system("ln -s ../makefile_fks_dir ./makefile") if matrix_element.virt_matrix_element: os.system("ln -s ../BinothLHA.f ./BinothLHA.f") elif OLP!='MadLoop': os.system("ln -s ../BinothLHA_OLP.f ./BinothLHA.f") else: os.system("ln -s ../BinothLHA_user.f ./BinothLHA.f") #import nexternal/leshouches in Source # ln('nexternal.inc', '../../Source', log=False) # ln('born_leshouche.inc', '../../Source', log=False) # Return to SubProcesses dir os.chdir(os.path.pardir) # Add subprocess to subproc.mg filename = 'subproc.mg' files.append_to_file(filename, self.write_subproc, borndir) os.chdir(cwd) # Generate info page gen_infohtml.make_info_html_nlo(self.dir_path) return calls #=========================================================================== # create the run_card #=========================================================================== def create_run_card(self, processes, history): """ """ run_card = banner_mod.RunCardNLO() run_card.create_default_for_process(self.proc_characteristic, history, processes) run_card.write(pjoin(self.dir_path, 'Cards', 'run_card_default.dat')) run_card.write(pjoin(self.dir_path, 'Cards', 'run_card.dat')) def pass_information_from_cmd(self, cmd): """pass information from the command interface to the exporter. Please do not modify any object of the interface from the exporter. """ self.proc_defs = cmd._curr_proc_defs if hasattr(cmd,'born_processes'): self.born_processes = cmd.born_processes else: self.born_processes = [] return def finalize(self, matrix_elements, history, mg5options, flaglist): """Finalize FKS directory by creating jpeg diagrams, html pages,proc_card_mg5.dat and madevent.tar.gz and create the MA5 card if necessary.""" devnull = os.open(os.devnull, os.O_RDWR) try: res = misc.call([mg5options['lhapdf'], '--version'], \ stdout=subprocess.PIPE, stderr=subprocess.PIPE) except Exception: res = 1 if res != 0: logger.info('The value for lhapdf in the current configuration does not ' + \ 'correspond to a valid executable.\nPlease set it correctly either in ' + \ 'input/mg5_configuration or with "set lhapdf /path/to/lhapdf-config" ' + \ 'and regenrate the process. \nTo avoid regeneration, edit the ' + \ ('%s/Cards/amcatnlo_configuration.txt file.\n' % self.dir_path ) + \ 'Note that you can still compile and run aMC@NLO with the built-in PDFs\n') compiler_dict = {'fortran': mg5options['fortran_compiler'], 'cpp': mg5options['cpp_compiler'], 'f2py': mg5options['f2py_compiler']} if 'nojpeg' in flaglist: makejpg = False else: makejpg = True output_dependencies = mg5options['output_dependencies'] self.proc_characteristic['grouped_matrix'] = False self.proc_characteristic['complex_mass_scheme'] = mg5options['complex_mass_scheme'] self.create_proc_charac() self.create_run_card(matrix_elements.get_processes(), history) # modelname = self.model.get('name') # if modelname == 'mssm' or modelname.startswith('mssm-'): # param_card = os.path.join(self.dir_path, 'Cards','param_card.dat') # mg5_param = os.path.join(self.dir_path, 'Source', 'MODEL', 'MG5_param.dat') # check_param_card.convert_to_mg5card(param_card, mg5_param) # check_param_card.check_valid_param_card(mg5_param) # # write the model functions get_mass/width_from_id filename = os.path.join(self.dir_path,'Source','MODEL','get_mass_width_fcts.f') makeinc = os.path.join(self.dir_path,'Source','MODEL','makeinc.inc') self.write_get_mass_width_file(writers.FortranWriter(filename), makeinc, self.model) # # Write maxconfigs.inc based on max of ME's/subprocess groups filename = os.path.join(self.dir_path,'Source','maxconfigs.inc') self.write_maxconfigs_file(writers.FortranWriter(filename), matrix_elements.get_max_configs()) # # Write maxparticles.inc based on max of ME's/subprocess groups filename = os.path.join(self.dir_path,'Source','maxparticles.inc') self.write_maxparticles_file(writers.FortranWriter(filename), matrix_elements.get_max_particles()) # Touch "done" file os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses')) # Check for compiler fcompiler_chosen = self.set_fortran_compiler(compiler_dict) ccompiler_chosen = self.set_cpp_compiler(compiler_dict['cpp']) old_pos = os.getcwd() os.chdir(os.path.join(self.dir_path, 'SubProcesses')) P_dir_list = [proc for proc in os.listdir('.') if os.path.isdir(proc) and \ proc[0] == 'P'] devnull = os.open(os.devnull, os.O_RDWR) # Convert the poscript in jpg files (if authorize) if makejpg: logger.info("Generate jpeg diagrams") for Pdir in P_dir_list: os.chdir(Pdir) subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')], stdout = devnull) os.chdir(os.path.pardir) # logger.info("Generate web pages") # Create the WebPage using perl script subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], \ stdout = devnull) os.chdir(os.path.pardir) # # obj = gen_infohtml.make_info_html(self.dir_path) # [mv(name, './HTML/') for name in os.listdir('.') if \ # (name.endswith('.html') or name.endswith('.jpg')) and \ # name != 'index.html'] # if online: # nb_channel = obj.rep_rule['nb_gen_diag'] # open(os.path.join('./Online'),'w').write(str(nb_channel)) # Write command history as proc_card_mg5 if os.path.isdir('Cards'): output_file = os.path.join('Cards', 'proc_card_mg5.dat') history.write(output_file) # Duplicate run_card and FO_analyse_card for card in ['run_card', 'FO_analyse_card', 'shower_card']: try: shutil.copy(pjoin(self.dir_path, 'Cards', card + '.dat'), pjoin(self.dir_path, 'Cards', card + '_default.dat')) except IOError: logger.warning("Failed to copy " + card + ".dat to default") subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], stdout = devnull) # Run "make" to generate madevent.tar.gz file if os.path.exists(pjoin('SubProcesses', 'subproc.mg')): if os.path.exists('amcatnlo.tar.gz'): os.remove('amcatnlo.tar.gz') subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'make_amcatnlo_tar')], stdout = devnull) # subprocess.call([os.path.join(old_pos, self.dir_path, 'bin', 'internal', 'gen_cardhtml-pl')], stdout = devnull) #return to the initial dir os.chdir(old_pos) # Setup stdHep # Find the correct fortran compiler base_compiler= ['FC=g77','FC=gfortran'] StdHep_path = pjoin(MG5DIR, 'vendor', 'StdHEP') if output_dependencies == 'external': # check if stdhep has to be compiled (only the first time) if not os.path.exists(pjoin(MG5DIR, 'vendor', 'StdHEP', 'lib', 'libstdhep.a')) or \ not os.path.exists(pjoin(MG5DIR, 'vendor', 'StdHEP', 'lib', 'libFmcfio.a')): if 'FC' not in os.environ or not os.environ['FC']: path = os.path.join(StdHep_path, 'src', 'make_opts') text = open(path).read() for base in base_compiler: text = text.replace(base,'FC=%s' % fcompiler_chosen) open(path, 'w').writelines(text) logger.info('Compiling StdHEP. This has to be done only once.') misc.compile(cwd = pjoin(MG5DIR, 'vendor', 'StdHEP')) logger.info('Done.') #then link the libraries in the exported dir files.ln(pjoin(StdHep_path, 'lib', 'libstdhep.a'), \ pjoin(self.dir_path, 'MCatNLO', 'lib')) files.ln(pjoin(StdHep_path, 'lib', 'libFmcfio.a'), \ pjoin(self.dir_path, 'MCatNLO', 'lib')) elif output_dependencies == 'internal': StdHEP_internal_path = pjoin(self.dir_path,'Source','StdHEP') shutil.copytree(StdHep_path,StdHEP_internal_path, symlinks=True) # Create the links to the lib folder linkfiles = ['libstdhep.a', 'libFmcfio.a'] for file in linkfiles: ln(pjoin(os.path.pardir,os.path.pardir,'Source','StdHEP','lib',file), os.path.join(self.dir_path, 'MCatNLO', 'lib')) if 'FC' not in os.environ or not os.environ['FC']: path = pjoin(StdHEP_internal_path, 'src', 'make_opts') text = open(path).read() for base in base_compiler: text = text.replace(base,'FC=%s' % fcompiler_chosen) open(path, 'w').writelines(text) # To avoid compiler version conflicts, we force a clean here misc.compile(['clean'],cwd = StdHEP_internal_path) elif output_dependencies == 'environment_paths': # Here the user chose to define the dependencies path in one of # his environmental paths libStdHep = misc.which_lib('libstdhep.a') libFmcfio = misc.which_lib('libFmcfio.a') if not libStdHep is None and not libFmcfio is None: logger.info('MG5_aMC is using StdHep installation found at %s.'%\ os.path.dirname(libStdHep)) ln(pjoin(libStdHep),pjoin(self.dir_path, 'MCatNLO', 'lib'),abspath=True) ln(pjoin(libFmcfio),pjoin(self.dir_path, 'MCatNLO', 'lib'),abspath=True) else: raise InvalidCmd("Could not find the location of the files"+\ " libstdhep.a and libFmcfio.a in you environment paths.") else: raise MadGraph5Error, 'output_dependencies option %s not recognized'\ %output_dependencies # Create the default MadAnalysis5 cards if 'madanalysis5_path' in self.opt and not \ self.opt['madanalysis5_path'] is None and not self.proc_defs is None: # When using processes = sum([me.get('processes') if not isinstance(me, str) else [] \ for me in matrix_elements.get('matrix_elements')],[]) # Try getting the processes from the generation info directly if no ME are # available (as it is the case for parallel generation if len(processes)==0: processes = self.born_processes if len(processes)==0: logger.warning( """MG5aMC could not provide to Madanalysis5 the list of processes generated. As a result, the default card will not be tailored to the process generated. This typically happens when using the 'low_mem_multicore_nlo_generation' NLO generation mode.""") # For now, simply assign all processes to each proc_defs. # That shouldn't really affect the default analysis card created by MA5 self.create_default_madanalysis5_cards( history, self.proc_defs, [processes,]*len(self.proc_defs), self.opt['madanalysis5_path'], pjoin(self.dir_path,'Cards'), levels =['hadron']) def write_real_from_born_configs(self, writer, matrix_element, fortran_model): """Writes the real_from_born_configs.inc file that contains the mapping to go for a given born configuration (that is used e.g. in the multi-channel phase-space integration to the corresponding real-emission diagram, i.e. the real emission diagram in which the combined ij is split in i_fks and j_fks.""" lines=[] lines2=[] max_links=0 born_me=matrix_element.born_matrix_element for iFKS, conf in enumerate(matrix_element.get_fks_info_list()): iFKS=iFKS+1 links=conf['fks_info']['rb_links'] max_links=max(max_links,len(links)) for i,diags in enumerate(links): if not i == diags['born_conf']: print links raise MadGraph5Error, "born_conf should be canonically ordered" real_configs=', '.join(['%d' % int(diags['real_conf']+1) for diags in links]) lines.append("data (real_from_born_conf(irfbc,%d),irfbc=1,%d) /%s/" \ % (iFKS,len(links),real_configs)) lines2.append("integer irfbc") lines2.append("integer real_from_born_conf(%d,%d)" \ % (max_links,len(matrix_element.get_fks_info_list()))) # Write the file writer.writelines(lines2+lines) #=============================================================================== # write_get_mass_width_file #=============================================================================== #test written def write_get_mass_width_file(self, writer, makeinc, model): """Write the get_mass_width_file.f file for MG4. Also update the makeinc.inc file """ mass_particles = [p for p in model['particles'] if p['mass'].lower() != 'zero'] width_particles = [p for p in model['particles'] if p['width'].lower() != 'zero'] iflines_mass = '' iflines_width = '' for i, part in enumerate(mass_particles): if i == 0: ifstring = 'if' else: ifstring = 'else if' if part['self_antipart']: iflines_mass += '%s (id.eq.%d) then\n' % \ (ifstring, part.get_pdg_code()) else: iflines_mass += '%s (id.eq.%d.or.id.eq.%d) then\n' % \ (ifstring, part.get_pdg_code(), part.get_anti_pdg_code()) iflines_mass += 'get_mass_from_id=abs(%s)\n' % part.get('mass') for i, part in enumerate(width_particles): if i == 0: ifstring = 'if' else: ifstring = 'else if' if part['self_antipart']: iflines_width += '%s (id.eq.%d) then\n' % \ (ifstring, part.get_pdg_code()) else: iflines_width += '%s (id.eq.%d.or.id.eq.%d) then\n' % \ (ifstring, part.get_pdg_code(), part.get_anti_pdg_code()) iflines_width += 'get_width_from_id=abs(%s)\n' % part.get('width') replace_dict = {'iflines_mass' : iflines_mass, 'iflines_width' : iflines_width} file = open(os.path.join(_file_path, \ 'iolibs/template_files/get_mass_width_fcts.inc')).read() file = file % replace_dict # Write the file writer.writelines(file) # update the makeinc makeinc_content = open(makeinc).read() makeinc_content = makeinc_content.replace('MODEL = ', 'MODEL = get_mass_width_fcts.o ') open(makeinc, 'w').write(makeinc_content) return def write_configs_and_props_info_declarations(self, writer, max_iconfig, max_leg_number, nfksconfs, fortran_model): """writes the declarations for the variables relevant for configs_and_props """ lines = [] lines.append("integer ifr,lmaxconfigs_used,max_branch_used") lines.append("parameter (lmaxconfigs_used=%4d)" % max_iconfig) lines.append("parameter (max_branch_used =%4d)" % -max_leg_number) lines.append("integer mapconfig_d(%3d,0:lmaxconfigs_used)" % nfksconfs) lines.append("integer iforest_d(%3d,2,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) lines.append("integer sprop_d(%3d,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) lines.append("integer tprid_d(%3d,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) lines.append("double precision pmass_d(%3d,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) lines.append("double precision pwidth_d(%3d,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) lines.append("integer pow_d(%3d,-max_branch_used:-1,lmaxconfigs_used)" % nfksconfs) writer.writelines(lines) def write_configs_and_props_info_file(self, filename, matrix_element): """writes the configs_and_props_info.inc file that cointains all the (real-emission) configurations (IFOREST) as well as the masses and widths of intermediate particles""" lines = [] lines.append("# C -> MAPCONFIG_D") lines.append("# F/D -> IFOREST_D") lines.append("# S -> SPROP_D") lines.append("# T -> TPRID_D") lines.append("# M -> PMASS_D/PWIDTH_D") lines.append("# P -> POW_D") lines2 = [] nconfs = len(matrix_element.get_fks_info_list()) (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() max_iconfig=0 max_leg_number=0 for iFKS, conf in enumerate(matrix_element.get_fks_info_list()): iFKS=iFKS+1 iconfig = 0 s_and_t_channels = [] mapconfigs = [] fks_matrix_element=matrix_element.real_processes[conf['n_me'] - 1].matrix_element base_diagrams = fks_matrix_element.get('base_amplitude').get('diagrams') model = fks_matrix_element.get('base_amplitude').get('process').get('model') minvert = min([max([len(vert.get('legs')) for vert in \ diag.get('vertices')]) for diag in base_diagrams]) lines.append("# ") lines.append("# nFKSprocess %d" % iFKS) for idiag, diag in enumerate(base_diagrams): if any([len(vert.get('legs')) > minvert for vert in diag.get('vertices')]): # Only 3-vertices allowed in configs.inc continue iconfig = iconfig + 1 helas_diag = fks_matrix_element.get('diagrams')[idiag] mapconfigs.append(helas_diag.get('number')) lines.append("# Diagram %d for nFKSprocess %d" % \ (helas_diag.get('number'),iFKS)) # Correspondance between the config and the amplitudes lines.append("C %4d %4d %4d " % (iFKS,iconfig, helas_diag.get('number'))) # Need to reorganize the topology so that we start with all # final state external particles and work our way inwards schannels, tchannels = helas_diag.get('amplitudes')[0].\ get_s_and_t_channels(ninitial, model, 990) s_and_t_channels.append([schannels, tchannels]) # Write out propagators for s-channel and t-channel vertices allchannels = schannels if len(tchannels) > 1: # Write out tchannels only if there are any non-trivial ones allchannels = schannels + tchannels for vert in allchannels: daughters = [leg.get('number') for leg in vert.get('legs')[:-1]] last_leg = vert.get('legs')[-1] lines.append("F %4d %4d %4d %4d" % \ (iFKS,last_leg.get('number'), iconfig, len(daughters))) for d in daughters: lines.append("D %4d" % d) if vert in schannels: lines.append("S %4d %4d %4d %10d" % \ (iFKS,last_leg.get('number'), iconfig, last_leg.get('id'))) elif vert in tchannels[:-1]: lines.append("T %4d %4d %4d %10d" % \ (iFKS,last_leg.get('number'), iconfig, abs(last_leg.get('id')))) # update what the array sizes (mapconfig,iforest,etc) will be max_leg_number = min(max_leg_number,last_leg.get('number')) max_iconfig = max(max_iconfig,iconfig) # Write out number of configs lines.append("# Number of configs for nFKSprocess %d" % iFKS) lines.append("C %4d %4d %4d" % (iFKS,0,iconfig)) # write the props.inc information lines2.append("# ") particle_dict = fks_matrix_element.get('processes')[0].get('model').\ get('particle_dict') for iconf, configs in enumerate(s_and_t_channels): for vertex in configs[0] + configs[1][:-1]: leg = vertex.get('legs')[-1] if leg.get('id') not in particle_dict: # Fake propagator used in multiparticle vertices pow_part = 0 else: particle = particle_dict[leg.get('id')] pow_part = 1 + int(particle.is_boson()) lines2.append("M %4d %4d %4d %10d " % \ (iFKS,leg.get('number'), iconf + 1, leg.get('id'))) lines2.append("P %4d %4d %4d %4d " % \ (iFKS,leg.get('number'), iconf + 1, pow_part)) # Write the file open(filename,'w').write('\n'.join(lines+lines2)) return max_iconfig, max_leg_number, nconfs def write_leshouche_info_declarations(self, writer, nfksconfs, maxproc, maxflow, nexternal, fortran_model): """writes the declarations for the variables relevant for leshouche_info """ lines = [] lines.append('integer maxproc_used, maxflow_used') lines.append('parameter (maxproc_used = %d)' % maxproc) lines.append('parameter (maxflow_used = %d)' % maxflow) lines.append('integer idup_d(%d,%d,maxproc_used)' % (nfksconfs, nexternal)) lines.append('integer mothup_d(%d,%d,%d,maxproc_used)' % (nfksconfs, 2, nexternal)) lines.append('integer icolup_d(%d,%d,%d,maxflow_used)' % (nfksconfs, 2, nexternal)) lines.append('integer niprocs_d(%d)' % (nfksconfs)) writer.writelines(lines) def write_genps(self, writer, maxproc,ngraphs,ncolor,maxflow, fortran_model): """writes the genps.inc file """ lines = [] lines.append("include 'maxparticles.inc'") lines.append("include 'maxconfigs.inc'") lines.append("integer maxproc,ngraphs,ncolor,maxflow") lines.append("parameter (maxproc=%d,ngraphs=%d,ncolor=%d,maxflow=%d)" % \ (maxproc,ngraphs,ncolor,maxflow)) writer.writelines(lines) def write_leshouche_info_file(self, filename, matrix_element): """writes the leshouche_info.inc file which contains the LHA informations for all the real emission processes """ lines = [] lines.append("# I -> IDUP_D") lines.append("# M -> MOTHUP_D") lines.append("# C -> ICOLUP_D") nfksconfs = len(matrix_element.get_fks_info_list()) (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() maxproc = 0 maxflow = 0 for i, conf in enumerate(matrix_element.get_fks_info_list()): # for i, real in enumerate(matrix_element.real_processes): (newlines, nprocs, nflows) = self.get_leshouche_lines( matrix_element.real_processes[conf['n_me'] - 1].matrix_element, i + 1) lines.extend(newlines) maxproc = max(maxproc, nprocs) maxflow = max(maxflow, nflows) # Write the file open(filename,'w').write('\n'.join(lines)) return nfksconfs, maxproc, maxflow, nexternal def write_pdf_wrapper(self, writer, matrix_element, fortran_model): """writes the wrapper which allows to chose among the different real matrix elements""" file = \ """double precision function dlum() implicit none include 'timing_variables.inc' integer nfksprocess common/c_nfksprocess/nfksprocess call cpu_time(tbefore) """ if matrix_element.real_processes: for n, info in enumerate(matrix_element.get_fks_info_list()): file += \ """if (nfksprocess.eq.%(n)d) then call dlum_%(n_me)d(dlum) else""" % {'n': n + 1, 'n_me' : info['n_me']} file += \ """ write(*,*) 'ERROR: invalid n in dlum :', nfksprocess stop endif call cpu_time(tAfter) tPDF = tPDF + (tAfter-tBefore) return end """ else: file+= \ """call dlum_0(dlum) call cpu_time(tAfter) tPDF = tPDF + (tAfter-tBefore) return end """ # Write the file writer.writelines(file) return 0 def write_real_me_wrapper(self, writer, matrix_element, fortran_model): """writes the wrapper which allows to chose among the different real matrix elements""" file = \ """subroutine smatrix_real(p, wgt) implicit none include 'nexternal.inc' double precision p(0:3, nexternal) double precision wgt integer nfksprocess common/c_nfksprocess/nfksprocess """ for n, info in enumerate(matrix_element.get_fks_info_list()): file += \ """if (nfksprocess.eq.%(n)d) then call smatrix_%(n_me)d(p, wgt) else""" % {'n': n + 1, 'n_me' : info['n_me']} if matrix_element.real_processes: file += \ """ write(*,*) 'ERROR: invalid n in real_matrix :', nfksprocess stop endif return end """ else: file += \ """ wgt=0d0 return end """ # Write the file writer.writelines(file) return 0 def draw_feynman_diagrams(self, matrix_element): """Create the ps files containing the feynman diagrams for the born process, as well as for all the real emission processes""" filename = 'born.ps' plot = draw.MultiEpsDiagramDrawer(matrix_element.born_matrix_element.\ get('base_amplitude').get('diagrams'), filename, model=matrix_element.born_matrix_element.\ get('processes')[0].get('model'), amplitude=True, diagram_type='born') plot.draw() for n, fksreal in enumerate(matrix_element.real_processes): filename = 'matrix_%d.ps' % (n + 1) plot = draw.MultiEpsDiagramDrawer(fksreal.matrix_element.\ get('base_amplitude').get('diagrams'), filename, model=fksreal.matrix_element.\ get('processes')[0].get('model'), amplitude=True, diagram_type='real') plot.draw() def write_real_matrix_elements(self, matrix_element, fortran_model): """writes the matrix_i.f files which contain the real matrix elements""" for n, fksreal in enumerate(matrix_element.real_processes): filename = 'matrix_%d.f' % (n + 1) self.write_matrix_element_fks(writers.FortranWriter(filename), fksreal.matrix_element, n + 1, fortran_model) def write_pdf_calls(self, matrix_element, fortran_model): """writes the parton_lum_i.f files which contain the real matrix elements. If no real emission existst, write the one for the born""" if matrix_element.real_processes: for n, fksreal in enumerate(matrix_element.real_processes): filename = 'parton_lum_%d.f' % (n + 1) self.write_pdf_file(writers.FortranWriter(filename), fksreal.matrix_element, n + 1, fortran_model) else: filename = 'parton_lum_0.f' self.write_pdf_file(writers.FortranWriter(filename), matrix_element.born_matrix_element, 0, fortran_model) def generate_born_fks_files(self, matrix_element, fortran_model, me_number, path): """generates the files needed for the born amplitude in the P* directory, which will be needed by the P* directories""" pathdir = os.getcwd() filename = 'born.f' calls_born, ncolor_born = \ self.write_born_fks(writers.FortranWriter(filename),\ matrix_element, fortran_model) filename = 'born_hel.f' self.write_born_hel(writers.FortranWriter(filename),\ matrix_element, fortran_model) filename = 'born_conf.inc' nconfigs, mapconfigs, s_and_t_channels = \ self.write_configs_file( writers.FortranWriter(filename), matrix_element.born_matrix_element, fortran_model) filename = 'born_props.inc' self.write_props_file(writers.FortranWriter(filename), matrix_element.born_matrix_element, fortran_model, s_and_t_channels) filename = 'born_decayBW.inc' self.write_decayBW_file(writers.FortranWriter(filename), s_and_t_channels) filename = 'born_leshouche.inc' nflows = self.write_leshouche_file(writers.FortranWriter(filename), matrix_element.born_matrix_element, fortran_model) filename = 'born_nhel.inc' self.write_born_nhel_file(writers.FortranWriter(filename), matrix_element.born_matrix_element, nflows, fortran_model, ncolor_born) filename = 'born_ngraphs.inc' self.write_ngraphs_file(writers.FortranWriter(filename), matrix_element.born_matrix_element.get_number_of_amplitudes()) filename = 'ncombs.inc' self.write_ncombs_file(writers.FortranWriter(filename), matrix_element.born_matrix_element, fortran_model) filename = 'born_maxamps.inc' maxamps = len(matrix_element.get('diagrams')) maxflows = ncolor_born self.write_maxamps_file(writers.FortranWriter(filename), maxamps, maxflows, max([len(matrix_element.get('processes')) for me in \ matrix_element.born_matrix_element]),1) filename = 'config_subproc_map.inc' self.write_config_subproc_map_file(writers.FortranWriter(filename), s_and_t_channels) filename = 'coloramps.inc' self.write_coloramps_file(writers.FortranWriter(filename), mapconfigs, matrix_element.born_matrix_element, fortran_model) #write the sborn_sf.f and the b_sf_files filename = ['sborn_sf.f', 'sborn_sf_dum.f'] for i, links in enumerate([matrix_element.color_links, []]): self.write_sborn_sf(writers.FortranWriter(filename[i]), links, fortran_model) self.color_link_files = [] for i in range(len(matrix_element.color_links)): filename = 'b_sf_%3.3d.f' % (i + 1) self.color_link_files.append(filename) self.write_b_sf_fks(writers.FortranWriter(filename), matrix_element, i, fortran_model) def generate_virtuals_from_OLP(self,process_list,export_path, OLP): """Generates the library for computing the loop matrix elements necessary for this process using the OLP specified.""" # Start by writing the BLHA order file virtual_path = pjoin(export_path,'OLP_virtuals') if not os.path.exists(virtual_path): os.makedirs(virtual_path) filename = os.path.join(virtual_path,'OLE_order.lh') self.write_lh_order(filename, process_list, OLP) fail_msg='Generation of the virtuals with %s failed.\n'%OLP+\ 'Please check the virt_generation.log file in %s.'\ %str(pjoin(virtual_path,'virt_generation.log')) # Perform some tasks specific to certain OLP's if OLP=='GoSam': cp(pjoin(self.mgme_dir,'Template','loop_material','OLP_specifics', 'GoSam','makevirt'),pjoin(virtual_path,'makevirt')) cp(pjoin(self.mgme_dir,'Template','loop_material','OLP_specifics', 'GoSam','gosam.rc'),pjoin(virtual_path,'gosam.rc')) ln(pjoin(export_path,'Cards','param_card.dat'),virtual_path) # Now generate the process logger.info('Generating the loop matrix elements with %s...'%OLP) virt_generation_log = \ open(pjoin(virtual_path,'virt_generation.log'), 'w') retcode = subprocess.call(['./makevirt'],cwd=virtual_path, stdout=virt_generation_log, stderr=virt_generation_log) virt_generation_log.close() # Check what extension is used for the share libraries on this system possible_other_extensions = ['so','dylib'] shared_lib_ext='so' for ext in possible_other_extensions: if os.path.isfile(pjoin(virtual_path,'Virtuals','lib', 'libgolem_olp.'+ext)): shared_lib_ext = ext # Now check that everything got correctly generated files_to_check = ['olp_module.mod',str(pjoin('lib', 'libgolem_olp.'+shared_lib_ext))] if retcode != 0 or any([not os.path.exists(pjoin(virtual_path, 'Virtuals',f)) for f in files_to_check]): raise fks_common.FKSProcessError(fail_msg) # link the library to the lib folder ln(pjoin(virtual_path,'Virtuals','lib','libgolem_olp.'+shared_lib_ext), pjoin(export_path,'lib')) # Specify in make_opts the right library necessitated by the OLP make_opts_content=open(pjoin(export_path,'Source','make_opts')).read() make_opts=open(pjoin(export_path,'Source','make_opts'),'w') if OLP=='GoSam': if platform.system().lower()=='darwin': # On mac the -rpath is not supported and the path of the dynamic # library is automatically wired in the executable make_opts_content=make_opts_content.replace('libOLP=', 'libOLP=-Wl,-lgolem_olp') else: # On other platforms the option , -rpath= path to libgolem.so is necessary # Using a relative path is not ideal because the file libgolem.so is not # copied on the worker nodes. # make_opts_content=make_opts_content.replace('libOLP=', # 'libOLP=-Wl,-rpath=../$(LIBDIR) -lgolem_olp') # Using the absolute path is working in the case where the disk of the # front end machine is mounted on all worker nodes as well. make_opts_content=make_opts_content.replace('libOLP=', 'libOLP=-Wl,-rpath='+str(pjoin(export_path,'lib'))+' -lgolem_olp') make_opts.write(make_opts_content) make_opts.close() # A priori this is generic to all OLP's # Parse the contract file returned and propagate the process label to # the include of the BinothLHA.f file proc_to_label = self.parse_contract_file( pjoin(virtual_path,'OLE_order.olc')) self.write_BinothLHA_inc(process_list,proc_to_label,\ pjoin(export_path,'SubProcesses')) # Link the contract file to within the SubProcess directory ln(pjoin(virtual_path,'OLE_order.olc'),pjoin(export_path,'SubProcesses')) def write_BinothLHA_inc(self, processes, proc_to_label, SubProcPath): """ Write the file Binoth_proc.inc in each SubProcess directory so as to provide the right process_label to use in the OLP call to get the loop matrix element evaluation. The proc_to_label is the dictionary of the format of the one returned by the function parse_contract_file.""" for proc in processes: name = "P%s"%proc.shell_string() proc_pdgs=(tuple([leg.get('id') for leg in proc.get('legs') if \ not leg.get('state')]), tuple([leg.get('id') for leg in proc.get('legs') if \ leg.get('state')])) incFile = open(pjoin(SubProcPath, name,'Binoth_proc.inc'),'w') try: incFile.write( """ INTEGER PROC_LABEL PARAMETER (PROC_LABEL=%d)"""%(proc_to_label[proc_pdgs])) except KeyError: raise fks_common.FKSProcessError('Could not found the target'+\ ' process %s > %s in '%(str(proc_pdgs[0]),str(proc_pdgs[1]))+\ ' the proc_to_label argument in write_BinothLHA_inc.') incFile.close() def parse_contract_file(self, contract_file_path): """ Parses the BLHA contract file, make sure all parameters could be understood by the OLP and return a mapping of the processes (characterized by the pdg's of the initial and final state particles) to their process label. The format of the mapping is {((in_pdgs),(out_pdgs)):proc_label}. """ proc_def_to_label = {} if not os.path.exists(contract_file_path): raise fks_common.FKSProcessError('Could not find the contract file'+\ ' OLE_order.olc in %s.'%str(contract_file_path)) comment_re=re.compile(r"^\s*#") proc_def_re=re.compile( r"^(?P(\s*-?\d+\s*)+)->(?P(\s*-?\d+\s*)+)\|"+ r"\s*(?P\d+)\s*(?P\d+)\s*$") line_OK_re=re.compile(r"^.*\|\s*OK") for line in file(contract_file_path): # Ignore comments if not comment_re.match(line) is None: continue # Check if it is a proc definition line proc_def = proc_def_re.match(line) if not proc_def is None: if int(proc_def.group('proc_class'))!=1: raise fks_common.FKSProcessError( 'aMCatNLO can only handle loop processes generated by the OLP which have only '+\ ' process class attribute. Found %s instead in: \n%s'\ %(proc_def.group('proc_class'),line)) in_pdgs=tuple([int(in_pdg) for in_pdg in \ proc_def.group('in_pdgs').split()]) out_pdgs=tuple([int(out_pdg) for out_pdg in \ proc_def.group('out_pdgs').split()]) proc_def_to_label[(in_pdgs,out_pdgs)]=\ int(proc_def.group('proc_label')) continue # For the other types of line, just make sure they end with | OK if line_OK_re.match(line) is None: raise fks_common.FKSProcessError( 'The OLP could not process the following line: \n%s'%line) return proc_def_to_label def generate_virt_directory(self, loop_matrix_element, fortran_model, dir_name): """writes the V**** directory inside the P**** directories specified in dir_name""" cwd = os.getcwd() matrix_element = loop_matrix_element # Create the MadLoop5_resources directory if not already existing dirpath = os.path.join(dir_name, 'MadLoop5_resources') try: os.mkdir(dirpath) except os.error as error: logger.warning(error.strerror + " " + dirpath) # Create the directory PN_xx_xxxxx in the specified path name = "V%s" % matrix_element.get('processes')[0].shell_string() dirpath = os.path.join(dir_name, name) try: os.mkdir(dirpath) except os.error as error: logger.warning(error.strerror + " " + dirpath) try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 logger.info('Creating files in directory %s' % name) # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() calls=self.write_loop_matrix_element_v4(None,matrix_element,fortran_model) # The born matrix element, if needed filename = 'born_matrix.f' calls = self.write_bornmatrix( writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'nexternal.inc' self.write_nexternal_file(writers.FortranWriter(filename), nexternal, ninitial) filename = 'pmass.inc' self.write_pmass_file(writers.FortranWriter(filename), matrix_element) filename = 'ngraphs.inc' self.write_ngraphs_file(writers.FortranWriter(filename), len(matrix_element.get_all_amplitudes())) filename = "loop_matrix.ps" plot = draw.MultiEpsDiagramDrawer(base_objects.DiagramList( matrix_element.get('base_amplitude').get('loop_diagrams')[:1000]), filename, model=matrix_element.get('processes')[0].get('model'), amplitude='') logger.info("Drawing loop Feynman diagrams for " + \ matrix_element.get('processes')[0].nice_string(print_weighted=False)) plot.draw() filename = "born_matrix.ps" plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\ get('born_diagrams'),filename,model=matrix_element.get('processes')[0].\ get('model'),amplitude='') logger.info("Generating born Feynman diagrams for " + \ matrix_element.get('processes')[0].nice_string(print_weighted=False)) plot.draw() # We also need to write the overall maximum quantities for this group # of processes in 'global_specs.inc'. In aMCatNLO, there is always # only one process, so this is trivial self.write_global_specs(matrix_element, output_path=pjoin(dirpath,'global_specs.inc')) open('unique_id.inc','w').write( """ integer UNIQUE_ID parameter(UNIQUE_ID=1)""") linkfiles = ['coupl.inc', 'mp_coupl.inc', 'mp_coupl_same_name.inc', 'cts_mprec.h', 'cts_mpc.h', 'MadLoopParamReader.f', 'MadLoopCommons.f','MadLoopParams.inc'] # We should move to MadLoop5_resources directory from the SubProcesses ln(pjoin(os.path.pardir,os.path.pardir,'MadLoopParams.dat'), pjoin('..','MadLoop5_resources')) for file in linkfiles: ln('../../%s' % file) os.system("ln -s ../../makefile_loop makefile") linkfiles = ['mpmodule.mod'] for file in linkfiles: ln('../../../lib/%s' % file) linkfiles = ['coef_specs.inc'] for file in linkfiles: ln('../../../Source/DHELAS/%s' % file) # Return to original PWD os.chdir(cwd) if not calls: calls = 0 return calls def get_qed_qcd_orders_from_weighted(self, nexternal, weighted): """computes the QED/QCD orders from the knowledge of the n of ext particles and of the weighted orders""" # n vertices = nexternal - 2 =QED + QCD # weighted = 2*QED + QCD QED = weighted - nexternal + 2 QCD = weighted - 2 * QED return QED, QCD #=============================================================================== # write_lh_order #=============================================================================== #test written def write_lh_order(self, filename, process_list, OLP='MadLoop'): """Creates the OLE_order.lh file. This function should be edited according to the OLP which is used. For now it is generic.""" if len(process_list)==0: raise fks_common.FKSProcessError('No matrix elements provided to '+\ 'the function write_lh_order.') return # We assume the orders to be common to all Subprocesses orders = process_list[0].get('orders') if 'QED' in orders.keys() and 'QCD' in orders.keys(): QED=orders['QED'] QCD=orders['QCD'] elif 'QED' in orders.keys(): QED=orders['QED'] QCD=0 elif 'QCD' in orders.keys(): QED=0 QCD=orders['QCD'] else: QED, QCD = self.get_qed_qcd_orders_from_weighted(\ len(process_list[0].get('legs')), orders['WEIGHTED']) replace_dict = {} replace_dict['mesq'] = 'CHaveraged' replace_dict['corr'] = ' '.join(process_list[0].\ get('perturbation_couplings')) replace_dict['irreg'] = 'CDR' replace_dict['aspow'] = QCD replace_dict['aepow'] = QED replace_dict['modelfile'] = './param_card.dat' replace_dict['params'] = 'alpha_s' proc_lines=[] for proc in process_list: proc_lines.append('%s -> %s' % \ (' '.join(str(l['id']) for l in proc['legs'] if not l['state']), ' '.join(str(l['id']) for l in proc['legs'] if l['state']))) replace_dict['pdgs'] = '\n'.join(proc_lines) replace_dict['symfin'] = 'Yes' content = \ "#OLE_order written by MadGraph5_aMC@NLO\n\ \n\ MatrixElementSquareType %(mesq)s\n\ CorrectionType %(corr)s\n\ IRregularisation %(irreg)s\n\ AlphasPower %(aspow)d\n\ AlphaPower %(aepow)d\n\ NJetSymmetrizeFinal %(symfin)s\n\ ModelFile %(modelfile)s\n\ Parameters %(params)s\n\ \n\ # process\n\ %(pdgs)s\n\ " % replace_dict file = open(filename, 'w') file.write(content) file.close return #=============================================================================== # write_born_fks #=============================================================================== # test written def write_born_fks(self, writer, fksborn, fortran_model): """Export a matrix element to a born.f file in MadFKS format""" matrix_element = fksborn.born_matrix_element if not matrix_element.get('processes') or \ not matrix_element.get('diagrams'): return 0 if not isinstance(writer, writers.FortranWriter): raise writers.FortranWriter.FortranWriterError(\ "writer not FortranWriter") # Set lowercase/uppercase Fortran code writers.FortranWriter.downcase = False replace_dict = {} # Extract version number and date from VERSION file info_lines = self.get_mg5_info_lines() replace_dict['info_lines'] = info_lines # Extract process info lines process_lines = self.get_process_info_lines(matrix_element) replace_dict['process_lines'] = process_lines # Extract ncomb ncomb = matrix_element.get_helicity_combinations() replace_dict['ncomb'] = ncomb # Extract helicity lines helicity_lines = self.get_helicity_lines(matrix_element) replace_dict['helicity_lines'] = helicity_lines # Extract IC line ic_line = self.get_ic_line(matrix_element) replace_dict['ic_line'] = ic_line # Extract overall denominator # Averaging initial state color, spin, and identical FS particles #den_factor_line = get_den_factor_line(matrix_element) # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() replace_dict['ngraphs'] = ngraphs # Extract nwavefuncs nwavefuncs = matrix_element.get_number_of_wavefunctions() replace_dict['nwavefuncs'] = nwavefuncs # Extract ncolor ncolor = max(1, len(matrix_element.get('color_basis'))) replace_dict['ncolor'] = ncolor # Extract color data lines color_data_lines = self.get_color_data_lines(matrix_element) replace_dict['color_data_lines'] = "\n".join(color_data_lines) # Extract helas calls helas_calls = fortran_model.get_matrix_element_calls(\ matrix_element) replace_dict['helas_calls'] = "\n".join(helas_calls) # Extract amp2 lines amp2_lines = self.get_amp2_lines(matrix_element) replace_dict['amp2_lines'] = '\n'.join(amp2_lines) # Extract JAMP lines jamp_lines = self.get_JAMP_lines(matrix_element) replace_dict['jamp_lines'] = '\n'.join(jamp_lines) # Set the size of Wavefunction if not self.model or any([p.get('spin') in [4,5] for p in self.model.get('particles') if p]): replace_dict['wavefunctionsize'] = 20 else: replace_dict['wavefunctionsize'] = 8 # Extract glu_ij_lines ij_lines = self.get_ij_lines(fksborn) replace_dict['ij_lines'] = '\n'.join(ij_lines) # Extract den_factor_lines den_factor_lines = self.get_den_factor_lines(fksborn) replace_dict['den_factor_lines'] = '\n'.join(den_factor_lines) # Extract the number of FKS process replace_dict['nconfs'] = max(len(fksborn.get_fks_info_list()),1) file = open(os.path.join(_file_path, \ 'iolibs/template_files/born_fks.inc')).read() file = file % replace_dict # Write the file writer.writelines(file) return len(filter(lambda call: call.find('#') != 0, helas_calls)), ncolor def write_born_hel(self, writer, fksborn, fortran_model): """Export a matrix element to a born_hel.f file in MadFKS format""" matrix_element = fksborn.born_matrix_element if not matrix_element.get('processes') or \ not matrix_element.get('diagrams'): return 0 if not isinstance(writer, writers.FortranWriter): raise writers.FortranWriter.FortranWriterError(\ "writer not FortranWriter") # Set lowercase/uppercase Fortran code writers.FortranWriter.downcase = False replace_dict = {} # Extract version number and date from VERSION file info_lines = self.get_mg5_info_lines() replace_dict['info_lines'] = info_lines # Extract process info lines process_lines = self.get_process_info_lines(matrix_element) replace_dict['process_lines'] = process_lines # Extract ncomb ncomb = matrix_element.get_helicity_combinations() replace_dict['ncomb'] = ncomb # Extract helicity lines helicity_lines = self.get_helicity_lines(matrix_element) replace_dict['helicity_lines'] = helicity_lines # Extract IC line ic_line = self.get_ic_line(matrix_element) replace_dict['ic_line'] = ic_line # Extract overall denominator # Averaging initial state color, spin, and identical FS particles #den_factor_line = get_den_factor_line(matrix_element) # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() replace_dict['ngraphs'] = ngraphs # Extract nwavefuncs nwavefuncs = matrix_element.get_number_of_wavefunctions() replace_dict['nwavefuncs'] = nwavefuncs # Extract ncolor ncolor = max(1, len(matrix_element.get('color_basis'))) replace_dict['ncolor'] = ncolor # Extract color data lines color_data_lines = self.get_color_data_lines(matrix_element) replace_dict['color_data_lines'] = "\n".join(color_data_lines) # Extract amp2 lines amp2_lines = self.get_amp2_lines(matrix_element) replace_dict['amp2_lines'] = '\n'.join(amp2_lines) # Extract JAMP lines jamp_lines = self.get_JAMP_lines(matrix_element) replace_dict['jamp_lines'] = '\n'.join(jamp_lines) # Extract den_factor_lines den_factor_lines = self.get_den_factor_lines(fksborn) replace_dict['den_factor_lines'] = '\n'.join(den_factor_lines) # Extract the number of FKS process replace_dict['nconfs'] = len(fksborn.get_fks_info_list()) file = open(os.path.join(_file_path, \ 'iolibs/template_files/born_fks_hel.inc')).read() file = file % replace_dict # Write the file writer.writelines(file) return #=============================================================================== # write_born_sf_fks #=============================================================================== #test written def write_sborn_sf(self, writer, color_links, fortran_model): """Creates the sborn_sf.f file, containing the calls to the different color linked borns""" replace_dict = {} nborns = len(color_links) ifkss = [] iborns = [] mms = [] nns = [] iflines = "\n" #header for the sborn_sf.f file file = """subroutine sborn_sf(p_born,m,n,wgt) implicit none include "nexternal.inc" double precision p_born(0:3,nexternal-1),wgt double complex wgt1(2) integer m,n \n""" if nborns > 0: for i, c_link in enumerate(color_links): iborn = i+1 iff = {True : 'if', False : 'elseif'}[i==0] m, n = c_link['link'] if m != n: iflines += \ "c b_sf_%(iborn)3.3d links partons %(m)d and %(n)d \n\ %(iff)s ((m.eq.%(m)d .and. n.eq.%(n)d).or.(m.eq.%(n)d .and. n.eq.%(m)d)) then \n\ call sb_sf_%(iborn)3.3d(p_born,wgt)\n\n" \ %{'m':m, 'n': n, 'iff': iff, 'iborn': iborn} else: iflines += \ "c b_sf_%(iborn)3.3d links partons %(m)d and %(n)d \n\ %(iff)s (m.eq.%(m)d .and. n.eq.%(n)d) then \n\ call sb_sf_%(iborn)3.3d(p_born,wgt)\n\n" \ %{'m':m, 'n': n, 'iff': iff, 'iborn': iborn} file += iflines + \ """else wgt = 0d0 endif return end""" elif nborns == 0: #write a dummy file file+=""" c This is a dummy function because c this subdir has no soft singularities wgt = 0d0 return end""" # Write the end of the file writer.writelines(file) #=============================================================================== # write_b_sf_fks #=============================================================================== #test written def write_b_sf_fks(self, writer, fksborn, i, fortran_model): """Create the b_sf_xxx.f file for the soft linked born in MadFKS format""" matrix_element = copy.copy(fksborn.born_matrix_element) if not matrix_element.get('processes') or \ not matrix_element.get('diagrams'): return 0 if not isinstance(writer, writers.FortranWriter): raise writers.FortranWriter.FortranWriterError(\ "writer not FortranWriter") # Set lowercase/uppercase Fortran code writers.FortranWriter.downcase = False iborn = i + 1 link = fksborn.color_links[i] replace_dict = {} replace_dict['iborn'] = iborn # Extract version number and date from VERSION file info_lines = self.get_mg5_info_lines() replace_dict['info_lines'] = info_lines # Extract process info lines process_lines = self.get_process_info_lines(matrix_element) replace_dict['process_lines'] = process_lines + \ "\nc spectators: %d %d \n" % tuple(link['link']) # Extract ncomb ncomb = matrix_element.get_helicity_combinations() replace_dict['ncomb'] = ncomb # Extract helicity lines helicity_lines = self.get_helicity_lines(matrix_element) replace_dict['helicity_lines'] = helicity_lines # Extract IC line ic_line = self.get_ic_line(matrix_element) replace_dict['ic_line'] = ic_line # Extract den_factor_lines den_factor_lines = self.get_den_factor_lines(fksborn) replace_dict['den_factor_lines'] = '\n'.join(den_factor_lines) # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() replace_dict['ngraphs'] = ngraphs # Extract nwavefuncs nwavefuncs = matrix_element.get_number_of_wavefunctions() replace_dict['nwavefuncs'] = nwavefuncs # Extract ncolor ncolor1 = max(1, len(link['orig_basis'])) replace_dict['ncolor1'] = ncolor1 ncolor2 = max(1, len(link['link_basis'])) replace_dict['ncolor2'] = ncolor2 # Extract color data lines color_data_lines = self.get_color_data_lines_from_color_matrix(\ link['link_matrix']) replace_dict['color_data_lines'] = "\n".join(color_data_lines) # Extract amp2 lines amp2_lines = self.get_amp2_lines(matrix_element) replace_dict['amp2_lines'] = '\n'.join(amp2_lines) # Extract JAMP lines jamp_lines = self.get_JAMP_lines(matrix_element) new_jamp_lines = [] for line in jamp_lines: line = string.replace(line, 'JAMP', 'JAMP1') new_jamp_lines.append(line) replace_dict['jamp1_lines'] = '\n'.join(new_jamp_lines) matrix_element.set('color_basis', link['link_basis'] ) jamp_lines = self.get_JAMP_lines(matrix_element) new_jamp_lines = [] for line in jamp_lines: line = string.replace(line, 'JAMP', 'JAMP2') new_jamp_lines.append(line) replace_dict['jamp2_lines'] = '\n'.join(new_jamp_lines) # Extract the number of FKS process replace_dict['nconfs'] = len(fksborn.get_fks_info_list()) file = open(os.path.join(_file_path, \ 'iolibs/template_files/b_sf_xxx_fks.inc')).read() file = file % replace_dict # Write the file writer.writelines(file) return 0 , ncolor1 #=============================================================================== # write_born_nhel_file #=============================================================================== #test written def write_born_nhel_file(self, writer, matrix_element, nflows, fortran_model, ncolor): """Write the born_nhel.inc file for MG4.""" ncomb = matrix_element.get_helicity_combinations() file = " integer max_bhel, max_bcol \n" file = file + "parameter (max_bhel=%d)\nparameter(max_bcol=%d)" % \ (ncomb, nflows) # Write the file writer.writelines(file) return True #=============================================================================== # write_fks_info_file #=============================================================================== def write_nfksconfigs_file(self, writer, fksborn, fortran_model): """Writes the content of nFKSconfigs.inc, which just gives the total FKS dirs as a parameter. nFKSconfigs is always >=1 (use a fake configuration for LOonly)""" replace_dict = {} replace_dict['nconfs'] = max(len(fksborn.get_fks_info_list()), 1) content = \ """ INTEGER FKS_CONFIGS PARAMETER (FKS_CONFIGS=%(nconfs)d) """ % replace_dict writer.writelines(content) #=============================================================================== # write_fks_info_file #=============================================================================== def write_fks_info_file(self, writer, fksborn, fortran_model): #test_written """Writes the content of fks_info.inc, which lists the informations on the possible splittings of the born ME. nconfs is always >=1 (use a fake configuration for LOonly). The fake configuration use an 'antigluon' (id -21, color=8) as i_fks and the last colored particle as j_fks.""" replace_dict = {} fks_info_list = fksborn.get_fks_info_list() replace_dict['nconfs'] = max(len(fks_info_list), 1) # this is for processes with 'real' or 'all' as NLO mode if len(fks_info_list) > 0: fks_i_values = ', '.join(['%d' % info['fks_info']['i'] \ for info in fks_info_list]) fks_j_values = ', '.join(['%d' % info['fks_info']['j'] \ for info in fks_info_list]) col_lines = [] pdg_lines = [] charge_lines = [] fks_j_from_i_lines = [] for i, info in enumerate(fks_info_list): col_lines.append( \ 'DATA (PARTICLE_TYPE_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /' \ % (i + 1, ', '.join('%d' % col for col in fksborn.real_processes[info['n_me']-1].colors) )) pdg_lines.append( \ 'DATA (PDG_TYPE_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /' \ % (i + 1, ', '.join('%d' % pdg for pdg in info['pdgs']))) charge_lines.append(\ 'DATA (PARTICLE_CHARGE_D(%d, IPOS), IPOS=1, NEXTERNAL) / %s /'\ % (i + 1, ', '.join('%19.15fd0' % charg\ for charg in fksborn.real_processes[info['n_me']-1].charges) )) fks_j_from_i_lines.extend(self.get_fks_j_from_i_lines(fksborn.real_processes[info['n_me']-1],\ i + 1)) else: # this is for 'LOonly', generate a fake FKS configuration with # - i_fks = nexternal, pdg type = -21 and color =8 # - j_fks = the last colored particle bornproc = fksborn.born_matrix_element.get('processes')[0] pdgs = [l.get('id') for l in bornproc.get('legs')] + [-21] colors = [l.get('color') for l in bornproc.get('legs')] + [8] charges = [0.] * len(colors) fks_i = len(colors) # use the first colored particle if it exists, or # just the first fks_j=1 for cpos, col in enumerate(colors[:-1]): if col != 1: fks_j = cpos+1 fks_i_values = str(fks_i) fks_j_values = str(fks_j) col_lines = ['DATA (PARTICLE_TYPE_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \ % ', '.join([str(col) for col in colors])] pdg_lines = ['DATA (PDG_TYPE_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \ % ', '.join([str(pdg) for pdg in pdgs])] charge_lines = ['DATA (PARTICLE_CHARGE_D(1, IPOS), IPOS=1, NEXTERNAL) / %s /' \ % ', '.join('%19.15fd0' % charg for charg in charges)] fks_j_from_i_lines = ['DATA (FKS_J_FROM_I_D(1, %d, JPOS), JPOS = 0, 1) / 1, %d /' \ % (fks_i, fks_j)] replace_dict['fks_i_line'] = "data fks_i_D / %s /" % fks_i_values replace_dict['fks_j_line'] = "data fks_j_D / %s /" % fks_j_values replace_dict['col_lines'] = '\n'.join(col_lines) replace_dict['pdg_lines'] = '\n'.join(pdg_lines) replace_dict['charge_lines'] = '\n'.join(charge_lines) replace_dict['fks_j_from_i_lines'] = '\n'.join(fks_j_from_i_lines) content = \ """ INTEGER IPOS, JPOS INTEGER FKS_I_D(%(nconfs)d), FKS_J_D(%(nconfs)d) INTEGER FKS_J_FROM_I_D(%(nconfs)d, NEXTERNAL, 0:NEXTERNAL) INTEGER PARTICLE_TYPE_D(%(nconfs)d, NEXTERNAL), PDG_TYPE_D(%(nconfs)d, NEXTERNAL) REAL*8 PARTICLE_CHARGE_D(%(nconfs)d, NEXTERNAL) %(fks_i_line)s %(fks_j_line)s %(fks_j_from_i_lines)s C C Particle type: C octet = 8, triplet = 3, singlet = 1 %(col_lines)s C C Particle type according to PDG: C %(pdg_lines)s C C Particle charge: C charge is set 0. with QCD corrections, which is irrelevant %(charge_lines)s """ % replace_dict if not isinstance(writer, writers.FortranWriter): raise writers.FortranWriter.FortranWriterError(\ "writer not FortranWriter") # Set lowercase/uppercase Fortran code writers.FortranWriter.downcase = False writer.writelines(content) return True #=============================================================================== # write_matrix_element_fks #=============================================================================== #test written def write_matrix_element_fks(self, writer, matrix_element, n, fortran_model): """Export a matrix element to a matrix.f file in MG4 madevent format""" if not matrix_element.get('processes') or \ not matrix_element.get('diagrams'): return 0,0 if not isinstance(writer, writers.FortranWriter): raise writers.FortranWriter.FortranWriterError(\ "writer not FortranWriter") # Set lowercase/uppercase Fortran code writers.FortranWriter.downcase = False replace_dict = {} replace_dict['N_me'] = n # Extract version number and date from VERSION file info_lines = self.get_mg5_info_lines() replace_dict['info_lines'] = info_lines # Extract process info lines process_lines = self.get_process_info_lines(matrix_element) replace_dict['process_lines'] = process_lines # Extract ncomb ncomb = matrix_element.get_helicity_combinations() replace_dict['ncomb'] = ncomb # Extract helicity lines helicity_lines = self.get_helicity_lines(matrix_element) replace_dict['helicity_lines'] = helicity_lines # Extract IC line ic_line = self.get_ic_line(matrix_element) replace_dict['ic_line'] = ic_line # Extract overall denominator # Averaging initial state color, spin, and identical FS particles den_factor_line = self.get_den_factor_line(matrix_element) replace_dict['den_factor_line'] = den_factor_line # Extract ngraphs ngraphs = matrix_element.get_number_of_amplitudes() replace_dict['ngraphs'] = ngraphs # Extract ncolor ncolor = max(1, len(matrix_element.get('color_basis'))) replace_dict['ncolor'] = ncolor # Extract color data lines color_data_lines = self.get_color_data_lines(matrix_element) replace_dict['color_data_lines'] = "\n".join(color_data_lines) # Extract helas calls helas_calls = fortran_model.get_matrix_element_calls(\ matrix_element) replace_dict['helas_calls'] = "\n".join(helas_calls) # Extract nwavefuncs (important to place after get_matrix_element_calls # so that 'me_id' is set) nwavefuncs = matrix_element.get_number_of_wavefunctions() replace_dict['nwavefuncs'] = nwavefuncs # Extract amp2 lines amp2_lines = self.get_amp2_lines(matrix_element) replace_dict['amp2_lines'] = '\n'.join(amp2_lines) # Set the size of Wavefunction if not self.model or any([p.get('spin') in [4,5] for p in self.model.get('particles') if p]): replace_dict['wavefunctionsize'] = 20 else: replace_dict['wavefunctionsize'] = 8 # Extract JAMP lines jamp_lines = self.get_JAMP_lines(matrix_element) replace_dict['jamp_lines'] = '\n'.join(jamp_lines) realfile = open(os.path.join(_file_path, \ 'iolibs/template_files/realmatrix_fks.inc')).read() realfile = realfile % replace_dict # Write the file writer.writelines(realfile) return len(filter(lambda call: call.find('#') != 0, helas_calls)), ncolor #=============================================================================== # write_pdf_file #=============================================================================== def write_pdf_file(self, writer, matrix_element, n, fortran_model): #test written """Write the auto_dsig.f file for MadFKS, which contains pdf call information""" if not matrix_element.get('processes') or \ not matrix_element.get('diagrams'): return 0 nexternal, ninitial = matrix_element.get_nexternal_ninitial() if ninitial < 1 or ninitial > 2: raise writers.FortranWriter.FortranWriterError, \ """Need ninitial = 1 or 2 to write auto_dsig file""" replace_dict = {} replace_dict['N_me'] = n # Extract version number and date from VERSION file info_lines = self.get_mg5_info_lines() replace_dict['info_lines'] = info_lines # Extract process info lines process_lines = self.get_process_info_lines(matrix_element) replace_dict['process_lines'] = process_lines pdf_vars, pdf_data, pdf_lines = \ self.get_pdf_lines_mir(matrix_element, ninitial, False, False) replace_dict['pdf_vars'] = pdf_vars replace_dict['pdf_data'] = pdf_data replace_dict['pdf_lines'] = pdf_lines pdf_vars_mirr, pdf_data_mirr, pdf_lines_mirr = \ self.get_pdf_lines_mir(matrix_element, ninitial, False, True) replace_dict['pdf_lines_mirr'] = pdf_lines_mirr file = open(os.path.join(_file_path, \ 'iolibs/template_files/parton_lum_n_fks.inc')).read() file = file % replace_dict # Write the file writer.writelines(file) #=============================================================================== # write_coloramps_file #=============================================================================== #test written def write_coloramps_file(self, writer, mapconfigs, matrix_element, fortran_model): """Write the coloramps.inc file for MadEvent""" lines = [] lines.append( "logical icolamp(%d,%d,1)" % \ (max(len(matrix_element.get('color_basis').keys()), 1), len(mapconfigs))) lines += self.get_icolamp_lines(mapconfigs, matrix_element, 1) # Write the file writer.writelines(lines) return True #=============================================================================== # write_leshouche_file #=============================================================================== #test written def write_leshouche_file(self, writer, matrix_element, fortran_model): """Write the leshouche.inc file for MG4""" # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() lines = [] for iproc, proc in enumerate(matrix_element.get('processes')): legs = proc.get_legs_with_decays() lines.append("DATA (IDUP(i,%d),i=1,%d)/%s/" % \ (iproc + 1, nexternal, ",".join([str(l.get('id')) for l in legs]))) for i in [1, 2]: lines.append("DATA (MOTHUP(%d,i,%3r),i=1,%2r)/%s/" % \ (i, iproc + 1, nexternal, ",".join([ "%3r" % 0 ] * ninitial + \ [ "%3r" % i ] * (nexternal - ninitial)))) # Here goes the color connections corresponding to the JAMPs # Only one output, for the first subproc! if iproc == 0: # If no color basis, just output trivial color flow if not matrix_element.get('color_basis'): for i in [1, 2]: lines.append("DATA (ICOLUP(%d,i, 1),i=1,%2r)/%s/" % \ (i, nexternal, ",".join([ "%3r" % 0 ] * nexternal))) color_flow_list = [] else: # First build a color representation dictionnary repr_dict = {} for l in legs: repr_dict[l.get('number')] = \ proc.get('model').get_particle(l.get('id')).get_color()\ * (-1)**(1+l.get('state')) # Get the list of color flows color_flow_list = \ matrix_element.get('color_basis').color_flow_decomposition(repr_dict, ninitial) # And output them properly for cf_i, color_flow_dict in enumerate(color_flow_list): for i in [0, 1]: lines.append("DATA (ICOLUP(%d,i,%3r),i=1,%2r)/%s/" % \ (i + 1, cf_i + 1, nexternal, ",".join(["%3r" % color_flow_dict[l.get('number')][i] \ for l in legs]))) # Write the file writer.writelines(lines) return len(color_flow_list) #=============================================================================== # write_configs_file #=============================================================================== #test_written def write_configs_file(self, writer, matrix_element, fortran_model): """Write the configs.inc file for MadEvent""" # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() lines = [] iconfig = 0 s_and_t_channels = [] mapconfigs = [] model = matrix_element.get('processes')[0].get('model') # new_pdg = model.get_first_non_pdg() base_diagrams = matrix_element.get('base_amplitude').get('diagrams') model = matrix_element.get('base_amplitude').get('process').get('model') minvert = min([max([len(vert.get('legs')) for vert in \ diag.get('vertices')]) for diag in base_diagrams]) for idiag, diag in enumerate(base_diagrams): if any([len(vert.get('legs')) > minvert for vert in diag.get('vertices')]): # Only 3-vertices allowed in configs.inc continue iconfig = iconfig + 1 helas_diag = matrix_element.get('diagrams')[idiag] mapconfigs.append(helas_diag.get('number')) lines.append("# Diagram %d, Amplitude %d" % \ (helas_diag.get('number'),helas_diag.get('amplitudes')[0]['number'])) # Correspondance between the config and the amplitudes lines.append("data mapconfig(%4d)/%4d/" % (iconfig, helas_diag.get('amplitudes')[0]['number'])) # Need to reorganize the topology so that we start with all # final state external particles and work our way inwards schannels, tchannels = helas_diag.get('amplitudes')[0].\ get_s_and_t_channels(ninitial, model, 990) s_and_t_channels.append([schannels, tchannels]) # Write out propagators for s-channel and t-channel vertices allchannels = schannels if len(tchannels) > 1: # Write out tchannels only if there are any non-trivial ones allchannels = schannels + tchannels for vert in allchannels: daughters = [leg.get('number') for leg in vert.get('legs')[:-1]] last_leg = vert.get('legs')[-1] lines.append("data (iforest(i,%3d,%4d),i=1,%d)/%s/" % \ (last_leg.get('number'), iconfig, len(daughters), ",".join(["%3d" % d for d in daughters]))) if vert in schannels: lines.append("data sprop(%4d,%4d)/%8d/" % \ (last_leg.get('number'), iconfig, last_leg.get('id'))) elif vert in tchannels[:-1]: lines.append("data tprid(%4d,%4d)/%8d/" % \ (last_leg.get('number'), iconfig, abs(last_leg.get('id')))) # Write out number of configs lines.append("# Number of configs") lines.append("data mapconfig(0)/%4d/" % iconfig) # Write the file writer.writelines(lines) return iconfig, mapconfigs, s_and_t_channels #=============================================================================== # write_decayBW_file #=============================================================================== #test written def write_decayBW_file(self, writer, s_and_t_channels): """Write the decayBW.inc file for MadEvent""" lines = [] booldict = {False: ".false.", True: ".false."} ####Changed by MZ 2011-11-23!!!! for iconf, config in enumerate(s_and_t_channels): schannels = config[0] for vertex in schannels: # For the resulting leg, pick out whether it comes from # decay or not, as given by the from_group flag leg = vertex.get('legs')[-1] lines.append("data gForceBW(%d,%d)/%s/" % \ (leg.get('number'), iconf + 1, booldict[leg.get('from_group')])) # Write the file writer.writelines(lines) return True #=============================================================================== # write_dname_file #=============================================================================== def write_dname_file(self, writer, matrix_element, fortran_model): """Write the dname.mg file for MG4""" line = "DIRNAME=P%s" % \ matrix_element.get('processes')[0].shell_string() # Write the file writer.write(line + "\n") return True #=============================================================================== # write_iproc_file #=============================================================================== def write_iproc_file(self, writer, me_number): """Write the iproc.dat file for MG4""" line = "%d" % (me_number + 1) # Write the file for line_to_write in writer.write_line(line): writer.write(line_to_write) return True #=============================================================================== # Helper functions #=============================================================================== #=============================================================================== # get_fks_j_from_i_lines #=============================================================================== def get_fks_j_from_i_lines(self, me, i = 0): #test written """generate the lines for fks.inc describing initializating the fks_j_from_i array""" lines = [] if not me.isfinite: for ii, js in me.fks_j_from_i.items(): if js: lines.append('DATA (FKS_J_FROM_I_D(%d, %d, JPOS), JPOS = 0, %d) / %d, %s /' \ % (i, ii, len(js), len(js), ', '.join(["%d" % j for j in js]))) else: lines.append('DATA (FKS_J_FROM_I_D(%d, JPOS), JPOS = 0, %d) / %d, %s /' \ % (2, 1, 1, '1')) lines.append('') return lines #=============================================================================== # get_leshouche_lines #=============================================================================== def get_leshouche_lines(self, matrix_element, ime): #test written """Write the leshouche.inc file for MG4""" # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() lines = [] for iproc, proc in enumerate(matrix_element.get('processes')): legs = proc.get_legs_with_decays() lines.append("I %4d %4d %s" % \ (ime, iproc + 1, " ".join([str(l.get('id')) for l in legs]))) for i in [1, 2]: lines.append("M %4d %4d %4d %s" % \ (ime, i, iproc + 1, " ".join([ "%3d" % 0 ] * ninitial + \ [ "%3d" % i ] * (nexternal - ninitial)))) # Here goes the color connections corresponding to the JAMPs # Only one output, for the first subproc! if iproc == 0: # If no color basis, just output trivial color flow if not matrix_element.get('color_basis'): for i in [1, 2]: lines.append("C %4d %4d 1 %s" % \ (ime, i, " ".join([ "%3d" % 0 ] * nexternal))) color_flow_list = [] nflow = 1 else: # First build a color representation dictionnary repr_dict = {} for l in legs: repr_dict[l.get('number')] = \ proc.get('model').get_particle(l.get('id')).get_color()\ * (-1)**(1+l.get('state')) # Get the list of color flows color_flow_list = \ matrix_element.get('color_basis').color_flow_decomposition(repr_dict, ninitial) # And output them properly for cf_i, color_flow_dict in enumerate(color_flow_list): for i in [0, 1]: lines.append("C %4d %4d %4d %s" % \ (ime, i + 1, cf_i + 1, " ".join(["%3d" % color_flow_dict[l.get('number')][i] \ for l in legs]))) nflow = len(color_flow_list) nproc = len(matrix_element.get('processes')) return lines, nproc, nflow #=============================================================================== # get_den_factor_lines #=============================================================================== def get_den_factor_lines(self, fks_born): """returns the lines with the information on the denominator keeping care of the identical particle factors in the various real emissions""" lines = [] info_list = fks_born.get_fks_info_list() if info_list: # if the reals have been generated, fill with the corresponding average factor lines.append('INTEGER IDEN_VALUES(%d)' % len(info_list)) lines.append('DATA IDEN_VALUES /' + \ ', '.join(['%d' % ( fks_born.born_matrix_element.get_denominator_factor() ) \ for info in info_list]) + '/') else: # otherwise use the born lines.append('INTEGER IDEN_VALUES(1)') lines.append('DATA IDEN_VALUES / %d /' \ % fks_born.born_matrix_element.get_denominator_factor()) return lines #=============================================================================== # get_ij_lines #=============================================================================== def get_ij_lines(self, fks_born): """returns the lines with the information on the particle number of the born that splits""" info_list = fks_born.get_fks_info_list() lines = [] if info_list: # if the reals have been generated, fill with the corresponding value of ij lines.append('INTEGER IJ_VALUES(%d)' % len(info_list)) lines.append('DATA IJ_VALUES /' + \ ', '.join(['%d' % info['fks_info']['ij'] for info in info_list]) + '/') else: #otherwise just put the first leg lines.append('INTEGER IJ_VALUES(1)') lines.append('DATA IJ_VALUES / 1 /') return lines def get_pdf_lines_mir(self, matrix_element, ninitial, subproc_group = False,\ mirror = False): #test written """Generate the PDF lines for the auto_dsig.f file""" processes = matrix_element.get('processes') model = processes[0].get('model') pdf_definition_lines = "" pdf_data_lines = "" pdf_lines = "" if ninitial == 1: pdf_lines = "PD(0) = 0d0\nIPROC = 0\n" for i, proc in enumerate(processes): process_line = proc.base_string() pdf_lines = pdf_lines + "IPROC=IPROC+1 ! " + process_line pdf_lines = pdf_lines + "\nPD(IPROC) = 1d0\n" pdf_lines = pdf_lines + "\nPD(0)=PD(0)+PD(IPROC)\n" else: # Pick out all initial state particles for the two beams initial_states = [sorted(list(set([p.get_initial_pdg(1) for \ p in processes]))), sorted(list(set([p.get_initial_pdg(2) for \ p in processes])))] # Prepare all variable names pdf_codes = dict([(p, model.get_particle(p).get_name()) for p in \ sum(initial_states,[])]) for key,val in pdf_codes.items(): pdf_codes[key] = val.replace('~','x').replace('+','p').replace('-','m') # Set conversion from PDG code to number used in PDF calls pdgtopdf = {21: 0, 22: 7} # Fill in missing entries of pdgtopdf for pdg in sum(initial_states,[]): if not pdg in pdgtopdf and not pdg in pdgtopdf.values(): pdgtopdf[pdg] = pdg elif pdg not in pdgtopdf and pdg in pdgtopdf.values(): # If any particle has pdg code 7, we need to use something else pdgtopdf[pdg] = 6000000 + pdg # Get PDF variable declarations for all initial states for i in [0,1]: pdf_definition_lines += "DOUBLE PRECISION " + \ ",".join(["%s%d" % (pdf_codes[pdg],i+1) \ for pdg in \ initial_states[i]]) + \ "\n" # Get PDF data lines for all initial states for i in [0,1]: pdf_data_lines += "DATA " + \ ",".join(["%s%d" % (pdf_codes[pdg],i+1) \ for pdg in initial_states[i]]) + \ "/%d*1D0/" % len(initial_states[i]) + \ "\n" # Get PDF values for the different initial states for i, init_states in enumerate(initial_states): if not mirror: ibeam = i + 1 else: ibeam = 2 - i if subproc_group: pdf_lines = pdf_lines + \ "IF (ABS(LPP(IB(%d))).GE.1) THEN\nLP=SIGN(1,LPP(IB(%d)))\n" \ % (ibeam, ibeam) else: pdf_lines = pdf_lines + \ "IF (ABS(LPP(%d)) .GE. 1) THEN\nLP=SIGN(1,LPP(%d))\n" \ % (ibeam, ibeam) for initial_state in init_states: if initial_state in pdf_codes.keys(): if subproc_group: if abs(pdgtopdf[initial_state]) <= 7: pdf_lines = pdf_lines + \ ("%s%d=PDG2PDF(ABS(LPP(IB(%d))),%d*LP," + \ "XBK(IB(%d)),DSQRT(Q2FACT(%d)))\n") % \ (pdf_codes[initial_state], i + 1, ibeam, pdgtopdf[initial_state], ibeam, ibeam) else: # setting other partons flavours outside quark, gluon, photon to be 0d0 pdf_lines = pdf_lines + \ ("c settings other partons flavours outside quark, gluon, photon to 0d0\n" + \ "%s%d=0d0\n") % \ (pdf_codes[initial_state],i + 1) else: if abs(pdgtopdf[initial_state]) <= 7: pdf_lines = pdf_lines + \ ("%s%d=PDG2PDF(ABS(LPP(%d)),%d*LP," + \ "XBK(%d),DSQRT(Q2FACT(%d)))\n") % \ (pdf_codes[initial_state], i + 1, ibeam, pdgtopdf[initial_state], ibeam, ibeam) else: # setting other partons flavours outside quark, gluon, photon to be 0d0 pdf_lines = pdf_lines + \ ("c settings other partons flavours outside quark, gluon, photon to 0d0\n" + \ "%s%d=0d0\n") % \ (pdf_codes[initial_state],i + 1) pdf_lines = pdf_lines + "ENDIF\n" # Add up PDFs for the different initial state particles pdf_lines = pdf_lines + "PD(0) = 0d0\nIPROC = 0\n" for proc in processes: process_line = proc.base_string() pdf_lines = pdf_lines + "IPROC=IPROC+1 ! " + process_line pdf_lines = pdf_lines + "\nPD(IPROC) = " for ibeam in [1, 2]: initial_state = proc.get_initial_pdg(ibeam) if initial_state in pdf_codes.keys(): pdf_lines = pdf_lines + "%s%d*" % \ (pdf_codes[initial_state], ibeam) else: pdf_lines = pdf_lines + "1d0*" # Remove last "*" from pdf_lines pdf_lines = pdf_lines[:-1] + "\n" # Remove last line break from pdf_lines return pdf_definition_lines[:-1], pdf_data_lines[:-1], pdf_lines[:-1] #test written def get_color_data_lines_from_color_matrix(self, color_matrix, n=6): """Return the color matrix definition lines for the given color_matrix. Split rows in chunks of size n.""" if not color_matrix: return ["DATA Denom(1)/1/", "DATA (CF(i,1),i=1,1) /1/"] else: ret_list = [] my_cs = color.ColorString() for index, denominator in \ enumerate(color_matrix.get_line_denominators()): # First write the common denominator for this color matrix line ret_list.append("DATA Denom(%i)/%i/" % (index + 1, denominator)) # Then write the numerators for the matrix elements num_list = color_matrix.get_line_numerators(index, denominator) for k in xrange(0, len(num_list), n): ret_list.append("DATA (CF(i,%3r),i=%3r,%3r) /%s/" % \ (index + 1, k + 1, min(k + n, len(num_list)), ','.join(["%5r" % i for i in num_list[k:k + n]]))) return ret_list #=========================================================================== # write_maxamps_file #=========================================================================== def write_maxamps_file(self, writer, maxamps, maxflows, maxproc,maxsproc): """Write the maxamps.inc file for MG4.""" file = " integer maxamps, maxflow, maxproc, maxsproc\n" file = file + "parameter (maxamps=%d, maxflow=%d)\n" % \ (maxamps, maxflows) file = file + "parameter (maxproc=%d, maxsproc=%d)" % \ (maxproc, maxsproc) # Write the file writer.writelines(file) return True #=============================================================================== # write_ncombs_file #=============================================================================== def write_ncombs_file(self, writer, matrix_element, fortran_model): # #test written """Write the ncombs.inc file for MadEvent.""" # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() # ncomb (used for clustering) is 2^(nexternal) file = " integer n_max_cl\n" file = file + "parameter (n_max_cl=%d)" % (2 ** (nexternal+1)) # Write the file writer.writelines(file) return True #=========================================================================== # write_config_subproc_map_file #=========================================================================== def write_config_subproc_map_file(self, writer, s_and_t_channels): """Write a dummy config_subproc.inc file for MadEvent""" lines = [] for iconfig in range(len(s_and_t_channels)): lines.append("DATA CONFSUB(1,%d)/1/" % \ (iconfig + 1)) # Write the file writer.writelines(lines) return True #=========================================================================== # write_colors_file #=========================================================================== def write_colors_file(self, writer, matrix_element): """Write the get_color.f file for MadEvent, which returns color for all particles used in the matrix element.""" try: matrix_elements=matrix_element.real_processes[0].matrix_element except IndexError: matrix_elements=[matrix_element.born_matrix_element] if isinstance(matrix_elements, helas_objects.HelasMatrixElement): matrix_elements = [matrix_elements] model = matrix_elements[0].get('processes')[0].get('model') # We need the both particle and antiparticle wf_ids, since the identity # depends on the direction of the wf. # loop on the real emissions wf_ids = set(sum([sum([sum([sum([[wf.get_pdg_code(),wf.get_anti_pdg_code()] \ for wf in d.get('wavefunctions')],[]) \ for d in me.get('diagrams')],[]) \ for me in [real_proc.matrix_element]],[])\ for real_proc in matrix_element.real_processes],[])) # and also on the born wf_ids = wf_ids.union(set(sum([sum([[wf.get_pdg_code(),wf.get_anti_pdg_code()] \ for wf in d.get('wavefunctions')],[]) \ for d in matrix_element.born_matrix_element.get('diagrams')],[]))) # loop on the real emissions leg_ids = set(sum([sum([sum([[l.get('id') for l in \ p.get_legs_with_decays()] for p in \ me.get('processes')], []) for me in \ [real_proc.matrix_element]], []) for real_proc in \ matrix_element.real_processes],[])) # and also on the born leg_ids = leg_ids.union(set(sum([[l.get('id') for l in \ p.get_legs_with_decays()] for p in \ matrix_element.born_matrix_element.get('processes')], []))) particle_ids = sorted(list(wf_ids.union(leg_ids))) lines = """function get_color(ipdg) implicit none integer get_color, ipdg if(ipdg.eq.%d)then get_color=%d return """ % (particle_ids[0], model.get_particle(particle_ids[0]).get_color()) for part_id in particle_ids[1:]: lines += """else if(ipdg.eq.%d)then get_color=%d return """ % (part_id, model.get_particle(part_id).get_color()) # Dummy particle for multiparticle vertices with pdg given by # first code not in the model lines += """else if(ipdg.eq.%d)then c This is dummy particle used in multiparticle vertices get_color=2 return """ % model.get_first_non_pdg() lines += """else write(*,*)'Error: No color given for pdg ',ipdg get_color=0 return endif end """ # Write the file writer.writelines(lines) return True #=============================================================================== # write_props_file #=============================================================================== #test_written def write_props_file(self, writer, matrix_element, fortran_model, s_and_t_channels): """Write the props.inc file for MadEvent. Needs input from write_configs_file. With respect to the parent routine, it has some more specific formats that allow the props.inc file to be read by the link program""" lines = [] particle_dict = matrix_element.get('processes')[0].get('model').\ get('particle_dict') for iconf, configs in enumerate(s_and_t_channels): for vertex in configs[0] + configs[1][:-1]: leg = vertex.get('legs')[-1] if leg.get('id') not in particle_dict: # Fake propagator used in multiparticle vertices mass = 'zero' width = 'zero' pow_part = 0 else: particle = particle_dict[leg.get('id')] # Get mass if particle.get('mass').lower() == 'zero': mass = particle.get('mass') else: mass = "abs(%s)" % particle.get('mass') # Get width if particle.get('width').lower() == 'zero': width = particle.get('width') else: width = "abs(%s)" % particle.get('width') pow_part = 1 + int(particle.is_boson()) lines.append("pmass(%3d,%4d) = %s" % \ (leg.get('number'), iconf + 1, mass)) lines.append("pwidth(%3d,%4d) = %s" % \ (leg.get('number'), iconf + 1, width)) lines.append("pow(%3d,%4d) = %d" % \ (leg.get('number'), iconf + 1, pow_part)) # Write the file writer.writelines(lines) return True #=========================================================================== # write_subproc #=========================================================================== def write_subproc(self, writer, subprocdir): """Append this subprocess to the subproc.mg file for MG4""" # Write line to file writer.write(subprocdir + "\n") return True #================================================================================= # Class for using the optimized Loop process #================================================================================= class ProcessOptimizedExporterFortranFKS(loop_exporters.LoopProcessOptimizedExporterFortranSA,\ ProcessExporterFortranFKS): """Class to take care of exporting a set of matrix elements to Fortran (v4) format.""" def finalize(self, *args, **opts): ProcessExporterFortranFKS.finalize(self, *args, **opts) #export_v4.ProcessExporterFortranSA.finalize(self, *args, **opts) #=============================================================================== # copy the Template in a new directory. #=============================================================================== def copy_fkstemplate(self): """create the directory run_name as a copy of the MadEvent Template, and clean the directory For now it is just the same as copy_v4template, but it will be modified """ mgme_dir = self.mgme_dir dir_path = self.dir_path clean =self.opt['clean'] #First copy the full template tree if dir_path doesn't exit if not os.path.isdir(dir_path): if not mgme_dir: raise MadGraph5Error, \ "No valid MG_ME path given for MG4 run directory creation." logger.info('initialize a new directory: %s' % \ os.path.basename(dir_path)) shutil.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True) # distutils.dir_util.copy_tree since dir_path already exists dir_util.copy_tree(pjoin(self.mgme_dir, 'Template', 'Common'), dir_path) # Copy plot_card for card in ['plot_card']: if os.path.isfile(pjoin(self.dir_path, 'Cards',card + '.dat')): try: shutil.copy(pjoin(self.dir_path, 'Cards', card + '.dat'), pjoin(self.dir_path, 'Cards', card + '_default.dat')) except IOError: logger.warning("Failed to copy " + card + ".dat to default") elif not os.path.isfile(os.path.join(dir_path, 'TemplateVersion.txt')): if not mgme_dir: raise MadGraph5Error, \ "No valid MG_ME path given for MG4 run directory creation." try: shutil.copy(os.path.join(mgme_dir, 'MGMEVersion.txt'), dir_path) except IOError: MG5_version = misc.get_pkg_info() open(os.path.join(dir_path, 'MGMEVersion.txt'), 'w').write( \ "5." + MG5_version['version']) #Ensure that the Template is clean if clean: logger.info('remove old information in %s' % os.path.basename(dir_path)) if os.environ.has_key('MADGRAPH_BASE'): subprocess.call([os.path.join('bin', 'internal', 'clean_template'), '--web'], cwd=dir_path) else: try: subprocess.call([os.path.join('bin', 'internal', 'clean_template')], \ cwd=dir_path) except Exception, why: raise MadGraph5Error('Failed to clean correctly %s: \n %s' \ % (os.path.basename(dir_path),why)) #Write version info MG_version = misc.get_pkg_info() open(os.path.join(dir_path, 'SubProcesses', 'MGVersion.txt'), 'w').write( MG_version['version']) # We must link the CutTools to the Library folder of the active Template self.link_CutTools(dir_path) # We must link the TIR to the Library folder of the active Template link_tir_libs=[] tir_libs=[] tir_include=[] for tir in self.all_tir: tir_dir="%s_dir"%tir libpath=getattr(self,tir_dir) libpath = self.link_TIR(os.path.join(self.dir_path, 'lib'), libpath,"lib%s.a"%tir,tir_name=tir) setattr(self,tir_dir,libpath) if libpath != "": if tir in ['pjfry','ninja','golem', 'samurai','collier']: # We should link dynamically when possible, so we use the original # location of these libraries. link_tir_libs.append('-L%s/ -l%s'%(libpath,tir)) tir_libs.append('%s/lib%s.$(libext)'%(libpath,tir)) # For Ninja, we must also link against OneLoop. if tir in ['ninja']: if not any(os.path.isfile(pjoin(libpath,'libavh_olo.%s'%ext)) for ext in ['a','dylib','so']): raise MadGraph5Error( "The OneLOop library 'libavh_olo.(a|dylib|so)' could no be found in path '%s'. Please place a symlink to it there."%libpath) link_tir_libs.append('-L%s/ -l%s'%(libpath,'avh_olo')) tir_libs.append('%s/lib%s.$(libext)'%(libpath,'avh_olo')) # We must add the corresponding includes for these TIR if tir in ['golem','samurai','ninja','collier']: trg_path = pjoin(os.path.dirname(libpath),'include') if os.path.isdir(trg_path): to_include = misc.find_includes_path(trg_path, self.include_names[tir]) else: to_include = None # Special possible location for collier if to_include is None and tir=='collier': to_include = misc.find_includes_path( pjoin(libpath,'modules'),self.include_names[tir]) if to_include is None: logger.error( 'Could not find the include directory for %s, looking in %s.\n' % (tir ,str(trg_path))+ 'Generation carries on but you will need to edit the include path by hand in the makefiles.') to_include = '' tir_include.append('-I %s'%to_include) else: link_tir_libs.append('-l%s'%tir) tir_libs.append('$(LIBDIR)lib%s.$(libext)'%tir) os.remove(os.path.join(self.dir_path,'SubProcesses','makefile_loop.inc')) cwd = os.getcwd() dirpath = os.path.join(self.dir_path, 'SubProcesses') try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 filename = 'makefile_loop' calls = self.write_makefile_TIR(writers.MakefileWriter(filename), link_tir_libs,tir_libs,tir_include=tir_include) os.remove(os.path.join(self.dir_path,'Source','make_opts.inc')) dirpath = os.path.join(self.dir_path, 'Source') try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 filename = 'make_opts' calls = self.write_make_opts(writers.MakefileWriter(filename), link_tir_libs,tir_libs) # Return to original PWD os.chdir(cwd) cwd = os.getcwd() dirpath = os.path.join(self.dir_path, 'SubProcesses') try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 # We add here the user-friendly MadLoop option setter. cpfiles= ["SubProcesses/MadLoopParamReader.f", "Cards/MadLoopParams.dat", "SubProcesses/MadLoopParams.inc"] for file in cpfiles: shutil.copy(os.path.join(self.loop_dir,'StandAlone/', file), os.path.join(self.dir_path, file)) shutil.copy(pjoin(self.dir_path, 'Cards','MadLoopParams.dat'), pjoin(self.dir_path, 'Cards','MadLoopParams_default.dat')) if os.path.exists(pjoin(self.dir_path, 'Cards', 'MadLoopParams.dat')): self.MadLoopparam = banner_mod.MadLoopParam(pjoin(self.dir_path, 'Cards', 'MadLoopParams.dat')) # write the output file self.MadLoopparam.write(pjoin(self.dir_path,"SubProcesses", "MadLoopParams.dat")) # We need minimal editing of MadLoopCommons.f MadLoopCommon = open(os.path.join(self.loop_dir,'StandAlone', "SubProcesses","MadLoopCommons.inc")).read() writer = writers.FortranWriter(os.path.join(self.dir_path, "SubProcesses","MadLoopCommons.f")) writer.writelines(MadLoopCommon%{ 'print_banner_commands':self.MadLoop_banner}, context={'collier_available':self.tir_available_dict['collier']}) writer.close() # link the files from the MODEL model_path = self.dir_path + '/Source/MODEL/' # Note that for the [real=] mode, these files are not present if os.path.isfile(os.path.join(model_path,'mp_coupl.inc')): ln(model_path + '/mp_coupl.inc', self.dir_path + '/SubProcesses') if os.path.isfile(os.path.join(model_path,'mp_coupl_same_name.inc')): ln(model_path + '/mp_coupl_same_name.inc', \ self.dir_path + '/SubProcesses') # Write the cts_mpc.h and cts_mprec.h files imported from CutTools self.write_mp_files(writers.FortranWriter('cts_mprec.h'),\ writers.FortranWriter('cts_mpc.h'),) self.copy_python_files() # We need to create the correct open_data for the pdf self.write_pdf_opendata() # Return to original PWD os.chdir(cwd) def generate_virt_directory(self, loop_matrix_element, fortran_model, dir_name): """writes the V**** directory inside the P**** directories specified in dir_name""" cwd = os.getcwd() matrix_element = loop_matrix_element # Create the MadLoop5_resources directory if not already existing dirpath = os.path.join(dir_name, 'MadLoop5_resources') try: os.mkdir(dirpath) except os.error as error: logger.warning(error.strerror + " " + dirpath) # Create the directory PN_xx_xxxxx in the specified path name = "V%s" % matrix_element.get('processes')[0].shell_string() dirpath = os.path.join(dir_name, name) try: os.mkdir(dirpath) except os.error as error: logger.warning(error.strerror + " " + dirpath) try: os.chdir(dirpath) except os.error: logger.error('Could not cd to directory %s' % dirpath) return 0 logger.info('Creating files in directory %s' % name) # Extract number of external particles (nexternal, ninitial) = matrix_element.get_nexternal_ninitial() calls=self.write_loop_matrix_element_v4(None,matrix_element,fortran_model) # We need a link to coefs.inc from DHELAS ln(pjoin(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc'), abspath=False, cwd=None) # The born matrix element, if needed filename = 'born_matrix.f' calls = self.write_bornmatrix( writers.FortranWriter(filename), matrix_element, fortran_model) filename = 'nexternal.inc' self.write_nexternal_file(writers.FortranWriter(filename), nexternal, ninitial) filename = 'pmass.inc' self.write_pmass_file(writers.FortranWriter(filename), matrix_element) filename = 'ngraphs.inc' self.write_ngraphs_file(writers.FortranWriter(filename), len(matrix_element.get_all_amplitudes())) filename = "loop_matrix.ps" writers.FortranWriter(filename).writelines("""C Post-helas generation loop-drawing is not ready yet.""") plot = draw.MultiEpsDiagramDrawer(base_objects.DiagramList( matrix_element.get('base_amplitude').get('loop_diagrams')[:1000]), filename, model=matrix_element.get('processes')[0].get('model'), amplitude='') logger.info("Drawing loop Feynman diagrams for " + \ matrix_element.get('processes')[0].nice_string(\ print_weighted=False)) plot.draw() filename = "born_matrix.ps" plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\ get('born_diagrams'), filename, model=matrix_element.get('processes')[0].\ get('model'), amplitude='') logger.info("Generating born Feynman diagrams for " + \ matrix_element.get('processes')[0].nice_string(\ print_weighted=False)) plot.draw() # We also need to write the overall maximum quantities for this group # of processes in 'global_specs.inc'. In aMCatNLO, there is always # only one process, so this is trivial self.write_global_specs(matrix_element, output_path=pjoin(dirpath,'global_specs.inc')) open('unique_id.inc','w').write( """ integer UNIQUE_ID parameter(UNIQUE_ID=1)""") linkfiles = ['coupl.inc', 'mp_coupl.inc', 'mp_coupl_same_name.inc', 'cts_mprec.h', 'cts_mpc.h', 'MadLoopParamReader.f', 'MadLoopParams.inc','MadLoopCommons.f'] for file in linkfiles: ln('../../%s' % file) os.system("ln -s ../../makefile_loop makefile") # We should move to MadLoop5_resources directory from the SubProcesses ln(pjoin(os.path.pardir,os.path.pardir,'MadLoopParams.dat'), pjoin('..','MadLoop5_resources')) linkfiles = ['mpmodule.mod'] for file in linkfiles: ln('../../../lib/%s' % file) linkfiles = ['coef_specs.inc'] for file in linkfiles: ln('../../../Source/DHELAS/%s' % file) # Return to original PWD os.chdir(cwd) if not calls: calls = 0 return calls #=============================================================================== # write_coef_specs #=============================================================================== def write_coef_specs_file(self, max_loop_vertex_ranks): """ writes the coef_specs.inc in the DHELAS folder. Should not be called in the non-optimized mode""" filename = os.path.join(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc') replace_dict = {} replace_dict['max_lwf_size'] = 4 replace_dict['vertex_max_coefs'] = max(\ [q_polynomial.get_number_of_coefs_for_rank(n) for n in max_loop_vertex_ranks]) IncWriter=writers.FortranWriter(filename,'w') IncWriter.writelines("""INTEGER MAXLWFSIZE PARAMETER (MAXLWFSIZE=%(max_lwf_size)d) INTEGER VERTEXMAXCOEFS PARAMETER (VERTEXMAXCOEFS=%(vertex_max_coefs)d)"""\ % replace_dict) IncWriter.close()