~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to madgraph/iolibs/export_v4.py

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
1
################################################################################
2
2
#
3
 
# Copyright (c) 2009 The MadGraph Development team and Contributors
 
3
# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors
4
4
#
5
 
# This file is a part of the MadGraph 5 project, an application which 
 
5
# This file is a part of the MadGraph5_aMC@NLO project, an application which 
6
6
# automatically generates Feynman diagrams and matrix elements for arbitrary
7
7
# high-energy processes in the Standard Model and beyond.
8
8
#
9
 
# It is subject to the MadGraph license which should accompany this 
 
9
# It is subject to the MadGraph5_aMC@NLO license which should accompany this 
10
10
# distribution.
11
11
#
12
 
# For more information, please visit: http://madgraph.phys.ucl.ac.be
 
12
# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
13
13
#
14
14
################################################################################
15
15
"""Methods and classes to export matrix elements to v4 format."""
16
16
 
17
17
import copy
 
18
from distutils import dir_util
18
19
import fractions
19
20
import glob
20
21
import logging
42
43
import madgraph.various.diagram_symmetry as diagram_symmetry
43
44
import madgraph.various.misc as misc
44
45
import madgraph.various.process_checks as process_checks
45
 
 
46
 
 
47
 
 
48
46
import aloha.create_aloha as create_aloha
49
47
import models.import_ufo as import_ufo
50
48
import models.write_param_card as param_writer
70
68
        self.mgme_dir = mgme_dir
71
69
        self.dir_path = dir_path
72
70
        self.model = None
 
71
 
73
72
        if opt:
74
73
            self.opt = opt
75
74
        else:
76
75
            self.opt = {'clean': False, 'complex_mass':False,
77
 
                        'export_format':'madevent'}
78
 
            
 
76
                        'export_format':'madevent', 'mp': False}
 
77
 
79
78
    #===========================================================================
80
79
    # copy the Template in a new directory.
81
80
    #===========================================================================
83
82
        """create the directory run_name as a copy of the MadEvent
84
83
        Template, and clean the directory
85
84
        """
86
 
 
 
85
        
87
86
        #First copy the full template tree if dir_path doesn't exit
88
87
        if not os.path.isdir(self.dir_path):
89
88
            assert self.mgme_dir, \
90
89
                     "No valid MG_ME path given for MG4 run directory creation."
91
90
            logger.info('initialize a new directory: %s' % \
92
91
                        os.path.basename(self.dir_path))
93
 
            shutil.copytree(pjoin(self.mgme_dir, 'Template'),
 
92
            shutil.copytree(pjoin(self.mgme_dir, 'Template/LO'),
94
93
                            self.dir_path, True)
 
94
            # distutils.dir_util.copy_tree since dir_path already exists
 
95
            dir_util.copy_tree(pjoin(self.mgme_dir, 'Template/Common'), 
 
96
                               self.dir_path)
95
97
            # Duplicate run_card and plot_card
96
98
            for card in ['run_card', 'plot_card']:
97
99
                try:
219
221
        mv(model_path + '/param_card.dat', self.dir_path + '/Cards/param_card_default.dat')
220
222
        ln(model_path + '/coupl.inc', self.dir_path + '/Source')
221
223
        ln(model_path + '/coupl.inc', self.dir_path + '/SubProcesses')
 
224
        self.make_source_links()
 
225
        
 
226
    def make_source_links(self):
 
227
        """ Create the links from the files in sources """
222
228
        ln(self.dir_path + '/Source/run.inc', self.dir_path + '/SubProcesses', log=False)
223
229
        ln(self.dir_path + '/Source/genps.inc', self.dir_path + '/SubProcesses', log=False)
224
230
        ln(self.dir_path + '/Source/maxconfigs.inc', self.dir_path + '/SubProcesses', log=False)
225
231
        ln(self.dir_path + '/Source/maxparticles.inc', self.dir_path + '/SubProcesses', log=False)
 
232
        ln(self.dir_path + '/Source/run_config.inc', self.dir_path + '/SubProcesses', log=False)
226
233
 
227
234
    #===========================================================================
228
235
    # export the helas routine
286
293
        return True
287
294
 
288
295
    #===========================================================================
 
296
    # write_nexternal_madspin
 
297
    #===========================================================================
 
298
    def write_nexternal_madspin(self, writer, nexternal, ninitial):
 
299
        """Write the nexternal_prod.inc file for madspin"""
 
300
 
 
301
        replace_dict = {}
 
302
 
 
303
        replace_dict['nexternal'] = nexternal
 
304
        replace_dict['ninitial'] = ninitial
 
305
 
 
306
        file = """ \
 
307
          integer    nexternal_prod
 
308
          parameter (nexternal_prod=%(nexternal)d)
 
309
          integer    nincoming_prod
 
310
          parameter (nincoming_prod=%(ninitial)d)""" % replace_dict
 
311
 
 
312
        # Write the file
 
313
        writer.writelines(file)
 
314
 
 
315
        return True
 
316
 
 
317
    #===========================================================================
289
318
    # write_nexternal_file
290
319
    #===========================================================================
291
320
    def write_nexternal_file(self, writer, nexternal, ninitial):
352
381
                             wanted_couplings = []):
353
382
        """ Create a full valid MG4 model from a MG5 model (coming from UFO)"""
354
383
 
 
384
        # Make sure aloha is in quadruple precision if needed
 
385
        old_aloha_mp=aloha.mp_precision
 
386
        aloha.mp_precision=self.opt['mp']
 
387
 
355
388
        # create the MODEL
356
389
        write_dir=pjoin(self.dir_path, 'Source', 'MODEL')
357
390
        model_builder = UFO_model_to_mg4(model, write_dir, self.opt)
358
391
        model_builder.build(wanted_couplings)
359
392
 
360
 
        # Create and write ALOHA Routine
361
 
        aloha_model = create_aloha.AbstractALOHAModel(model.get('name'))
 
393
        # Backup the loop mode, because it can be changed in what follows.
 
394
        old_loop_mode = aloha.loop_mode
 
395
 
 
396
        # Create the aloha model or use the existing one (for loop exporters
 
397
        # this is useful as the aloha model will be used again in the 
 
398
        # LoopHelasMatrixElements generated). We do not save the model generated
 
399
        # here if it didn't exist already because it would be a waste of
 
400
        # memory for tree level applications since aloha is only needed at the
 
401
        # time of creating the aloha fortran subroutines.
 
402
        if hasattr(self, 'aloha_model'):
 
403
            aloha_model = self.aloha_model
 
404
        else:
 
405
            aloha_model = create_aloha.AbstractALOHAModel(model.get('name'))            
362
406
        aloha_model.add_Lorentz_object(model.get('lorentz'))
 
407
 
 
408
        # Compute the subroutines
363
409
        if wanted_lorentz:
364
410
            aloha_model.compute_subset(wanted_lorentz)
365
411
        else:
366
412
            aloha_model.compute_all(save=False)
 
413
 
 
414
        # Write them out
367
415
        write_dir=pjoin(self.dir_path, 'Source', 'DHELAS')
368
416
        aloha_model.write(write_dir, 'Fortran')
369
417
 
 
418
        # Revert the original aloha loop mode
 
419
        aloha.loop_mode = old_loop_mode
 
420
 
370
421
        #copy Helas Template
371
422
        cp(MG5DIR + '/aloha/template_files/Makefile_F', write_dir+'/makefile')
372
 
        for filename in os.listdir(pjoin(MG5DIR,'aloha','template_files')):
373
 
            if not filename.lower().endswith('.f'):
374
 
                continue
375
 
            cp((MG5DIR + '/aloha/template_files/' + filename), write_dir)
 
423
        if any([any(['L' in tag for tag in d[1]]) for d in wanted_lorentz]):
 
424
            cp(MG5DIR + '/aloha/template_files/aloha_functions_loop.f', write_dir+'/aloha_functions.f')
 
425
            aloha_model.loop_mode = False
 
426
        else:
 
427
            cp(MG5DIR + '/aloha/template_files/aloha_functions.f', write_dir+'/aloha_functions.f')
376
428
        create_aloha.write_aloha_file_inc(write_dir, '.f', '.o')
377
429
 
378
430
        # Make final link in the Process
379
431
        self.make_model_symbolic_link()
 
432
        
 
433
        # Re-establish original aloha mode
 
434
        aloha.mp_precision=old_aloha_mp
380
435
 
381
436
    #===========================================================================
382
437
    # Helper functions
388
443
        info = misc.get_pkg_info()
389
444
        info_lines = ""
390
445
        if info and info.has_key('version') and  info.has_key('date'):
391
 
            info_lines = "#  Generated by MadGraph 5 v. %s, %s\n" % \
 
446
            info_lines = "#  Generated by MadGraph5_aMC@NLO v. %s, %s\n" % \
392
447
                         (info['version'], info['date'])
393
448
            info_lines = info_lines + \
394
 
                         "#  By the MadGraph Development Team\n" + \
395
 
                         "#  Please visit us at https://launchpad.net/madgraph5"
 
449
                         "#  By the MadGraph5_aMC@NLO Development Team\n" + \
 
450
                         "#  Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch"
396
451
        else:
397
 
            info_lines = "#  Generated by MadGraph 5\n" + \
398
 
                         "#  By the MadGraph Development Team\n" + \
399
 
                         "#  Please visit us at https://launchpad.net/madgraph5"        
 
452
            info_lines = "#  Generated by MadGraph5_aMC@NLO\n" + \
 
453
                         "#  By the MadGraph5_aMC@NLO Development Team\n" + \
 
454
                         "#  Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch"        
400
455
 
401
456
        return info_lines
402
457
 
407
462
                         for process in matrix_element.get('processes')])
408
463
 
409
464
 
410
 
    def get_helicity_lines(self, matrix_element):
 
465
    def get_helicity_lines(self, matrix_element,array_name='NHEL'):
411
466
        """Return the Helicity matrix definition lines for this matrix element"""
412
467
 
413
468
        helicity_line_list = []
417
472
            int_list = [i, len(helicities)]
418
473
            int_list.extend(helicities)
419
474
            helicity_line_list.append(\
420
 
                ("DATA (NHEL(I,%4r),I=1,%d) /" + \
 
475
                ("DATA ("+array_name+"(I,%4r),I=1,%d) /" + \
421
476
                 ",".join(['%2r'] * len(helicities)) + "/") % tuple(int_list))
422
477
 
423
478
        return "\n".join(helicity_line_list)
580
635
 
581
636
        return ret_lines
582
637
 
583
 
    def get_JAMP_lines(self, matrix_element):
584
 
        """Return the JAMP = sum(fermionfactor * AMP(i)) lines"""
585
 
 
586
 
        res_list = []
587
 
 
 
638
    #===========================================================================
 
639
    # Returns the data statements initializing the coeffictients for the JAMP
 
640
    # decomposition. It is used when the JAMP initialization is decided to be 
 
641
    # done through big arrays containing the projection coefficients.
 
642
    #===========================================================================    
 
643
    def get_JAMP_coefs(self, color_amplitudes, color_basis=None, tag_letter="",\
 
644
                       n=50, Nc_value=3):
 
645
        """This functions return the lines defining the DATA statement setting
 
646
        the coefficients building the JAMPS out of the AMPS. Split rows in
 
647
        bunches of size n.
 
648
        One can specify the color_basis from which the color amplitudes originates
 
649
        so that there are commentaries telling what color structure each JAMP
 
650
        corresponds to."""
 
651
        
 
652
        if(not isinstance(color_amplitudes,list) or 
 
653
           not (color_amplitudes and isinstance(color_amplitudes[0],list))):
 
654
                raise MadGraph5Error, "Incorrect col_amps argument passed to get_JAMP_coefs"
 
655
 
 
656
        res_list = []
 
657
        my_cs = color.ColorString()
 
658
        for index, coeff_list in enumerate(color_amplitudes):
 
659
            # Create the list of the complete numerical coefficient.
 
660
            coefs_list=[coefficient[0][0]*coefficient[0][1]*\
 
661
                        (fractions.Fraction(Nc_value)**coefficient[0][3]) for \
 
662
                        coefficient in coeff_list]
 
663
            # Create the list of the numbers of the contributing amplitudes.
 
664
            # Mutliply by -1 for those which have an imaginary coefficient.
 
665
            ampnumbers_list=[coefficient[1]*(-1 if coefficient[0][2] else 1) \
 
666
                              for coefficient in coeff_list]
 
667
            # Find the common denominator.      
 
668
            commondenom=abs(reduce(fractions.gcd, coefs_list).denominator)
 
669
            num_list=[(coefficient*commondenom).numerator \
 
670
                      for coefficient in coefs_list]
 
671
            res_list.append("DATA NCONTRIBAMPS%s(%i)/%i/"%(tag_letter,\
 
672
                                                         index+1,len(num_list)))
 
673
            res_list.append("DATA DENOMCCOEF%s(%i)/%i/"%(tag_letter,\
 
674
                                                         index+1,commondenom))
 
675
            if color_basis:
 
676
                my_cs.from_immutable(sorted(color_basis.keys())[index])
 
677
                res_list.append("C %s" % repr(my_cs))
 
678
            for k in xrange(0, len(num_list), n):
 
679
                res_list.append("DATA (NUMCCOEF%s(%3r,i),i=%6r,%6r) /%s/" % \
 
680
                    (tag_letter,index + 1, k + 1, min(k + n, len(num_list)),
 
681
                                 ','.join(["%6r" % i for i in num_list[k:k + n]])))
 
682
                res_list.append("DATA (AMPNUMBERS%s(%3r,i),i=%6r,%6r) /%s/" % \
 
683
                    (tag_letter,index + 1, k + 1, min(k + n, len(num_list)),
 
684
                                 ','.join(["%6r" % i for i in ampnumbers_list[k:k + n]])))
 
685
                pass
 
686
        return res_list
 
687
 
 
688
    def get_JAMP_lines(self, col_amps, basename="JAMP(", basename2="AMP(", split=-1):
 
689
        """Return the JAMP = sum(fermionfactor * AMP(i)) lines from col_amps 
 
690
        defined as a matrix element or directly as a color_amplitudes dictionary"""
 
691
 
 
692
        # Let the user call get_JAMP_lines directly from a MatrixElement or from
 
693
        # the color amplitudes lists.
 
694
        if(isinstance(col_amps,helas_objects.HelasMatrixElement)):
 
695
            color_amplitudes=col_amps.get_color_amplitudes()
 
696
        elif(isinstance(col_amps,list)):
 
697
            if(col_amps and isinstance(col_amps[0],list)):
 
698
                color_amplitudes=col_amps
 
699
            else:
 
700
                raise MadGraph5Error, "Incorrect col_amps argument passed to get_JAMP_lines"
 
701
        else:
 
702
            raise MadGraph5Error, "Incorrect col_amps argument passed to get_JAMP_lines"
 
703
 
 
704
 
 
705
        res_list = []
588
706
        for i, coeff_list in \
589
 
                enumerate(matrix_element.get_color_amplitudes()):
590
 
 
591
 
            res = "JAMP(%i)=" % (i + 1)
592
 
 
593
 
            # Optimization: if all contributions to that color basis element have
594
 
            # the same coefficient (up to a sign), put it in front
595
 
            list_fracs = [abs(coefficient[0][1]) for coefficient in coeff_list]
596
 
            common_factor = False
597
 
            diff_fracs = list(set(list_fracs))
598
 
            if len(diff_fracs) == 1 and abs(diff_fracs[0]) != 1:
599
 
                common_factor = True
600
 
                global_factor = diff_fracs[0]
601
 
                res = res + '%s(' % self.coeff(1, global_factor, False, 0)
602
 
 
603
 
            for (coefficient, amp_number) in coeff_list:
 
707
                enumerate(color_amplitudes):
 
708
            # Break the JAMP definition into 'n=split' pieces to avoid having
 
709
            # arbitrarly long lines.
 
710
            first=True
 
711
            n = (len(coeff_list)+1 if split<=0 else split) 
 
712
            while coeff_list!=[]:
 
713
                coefs=coeff_list[:n]
 
714
                coeff_list=coeff_list[n:]
 
715
                res = ((basename+"%i)=") % (i + 1)) + \
 
716
                      (((basename+"%i)") % (i + 1)) if not first and split>0 else '')
 
717
                first=False
 
718
                # Optimization: if all contributions to that color basis element have
 
719
                # the same coefficient (up to a sign), put it in front
 
720
                list_fracs = [abs(coefficient[0][1]) for coefficient in coefs]
 
721
                common_factor = False
 
722
                diff_fracs = list(set(list_fracs))
 
723
                if len(diff_fracs) == 1 and abs(diff_fracs[0]) != 1:
 
724
                    common_factor = True
 
725
                    global_factor = diff_fracs[0]
 
726
                    res = res + '%s(' % self.coeff(1, global_factor, False, 0)
 
727
    
 
728
                for (coefficient, amp_number) in coefs:
 
729
                    if common_factor:
 
730
                        res = (res + "%s" + basename2 + "%d)") % \
 
731
                                                   (self.coeff(coefficient[0],
 
732
                                                   coefficient[1] / abs(coefficient[1]),
 
733
                                                   coefficient[2],
 
734
                                                   coefficient[3]),
 
735
                                                   amp_number)
 
736
                    else:
 
737
                        res = (res + "%s" + basename2 + "%d)") % (self.coeff(coefficient[0],
 
738
                                                   coefficient[1],
 
739
                                                   coefficient[2],
 
740
                                                   coefficient[3]),
 
741
                                                   amp_number)
 
742
    
604
743
                if common_factor:
605
 
                    res = res + "%sAMP(%d)" % (self.coeff(coefficient[0],
606
 
                                               coefficient[1] / abs(coefficient[1]),
607
 
                                               coefficient[2],
608
 
                                               coefficient[3]),
609
 
                                               amp_number)
610
 
                else:
611
 
                    res = res + "%sAMP(%d)" % (self.coeff(coefficient[0],
612
 
                                               coefficient[1],
613
 
                                               coefficient[2],
614
 
                                               coefficient[3]),
615
 
                                               amp_number)
616
 
 
617
 
            if common_factor:
618
 
                res = res + ')'
619
 
 
620
 
            res_list.append(res)
 
744
                    res = res + ')'
 
745
    
 
746
                res_list.append(res)
621
747
 
622
748
        return res_list
623
749
 
727
853
        # Remove last line break from the return variables
728
854
        return pdf_definition_lines[:-1], pdf_data_lines[:-1], pdf_lines[:-1]
729
855
 
 
856
    #===========================================================================
 
857
    # write_props_file
 
858
    #===========================================================================
 
859
    def write_props_file(self, writer, matrix_element, s_and_t_channels):
 
860
        """Write the props.inc file for MadEvent. Needs input from
 
861
        write_configs_file."""
 
862
 
 
863
        lines = []
 
864
 
 
865
        particle_dict = matrix_element.get('processes')[0].get('model').\
 
866
                        get('particle_dict')
 
867
 
 
868
        for iconf, configs in enumerate(s_and_t_channels):
 
869
            for vertex in configs[0] + configs[1][:-1]:
 
870
                leg = vertex.get('legs')[-1]
 
871
                if leg.get('id') not in particle_dict:
 
872
                    # Fake propagator used in multiparticle vertices
 
873
                    mass = 'zero'
 
874
                    width = 'zero'
 
875
                    pow_part = 0
 
876
                else:
 
877
                    particle = particle_dict[leg.get('id')]
 
878
                    # Get mass
 
879
                    if particle.get('mass').lower() == 'zero':
 
880
                        mass = particle.get('mass')
 
881
                    else:
 
882
                        mass = "abs(%s)" % particle.get('mass')
 
883
                    # Get width
 
884
                    if particle.get('width').lower() == 'zero':
 
885
                        width = particle.get('width')
 
886
                    else:
 
887
                        width = "abs(%s)" % particle.get('width')
 
888
 
 
889
                    pow_part = 1 + int(particle.is_boson())
 
890
 
 
891
                lines.append("prmass(%d,%d)  = %s" % \
 
892
                             (leg.get('number'), iconf + 1, mass))
 
893
                lines.append("prwidth(%d,%d) = %s" % \
 
894
                             (leg.get('number'), iconf + 1, width))
 
895
                lines.append("pow(%d,%d) = %d" % \
 
896
                             (leg.get('number'), iconf + 1, pow_part))
 
897
 
 
898
        # Write the file
 
899
        writer.writelines(lines)
 
900
 
 
901
        return True
 
902
 
 
903
    #===========================================================================
 
904
    # write_configs_file
 
905
    #===========================================================================
 
906
    def write_configs_file(self, writer, matrix_element):
 
907
        """Write the configs.inc file for MadEvent"""
 
908
 
 
909
        # Extract number of external particles
 
910
        (nexternal, ninitial) = matrix_element.get_nexternal_ninitial()
 
911
 
 
912
        configs = [(i+1, d) for i,d in enumerate(matrix_element.get('diagrams'))]
 
913
        mapconfigs = [c[0] for c in configs]
 
914
        model = matrix_element.get('processes')[0].get('model')
 
915
        return mapconfigs, self.write_configs_file_from_diagrams(writer,
 
916
                                                            [[c[1]] for c in configs],
 
917
                                                            mapconfigs,
 
918
                                                            nexternal, ninitial,
 
919
                                                            model)
 
920
 
 
921
    #===========================================================================
 
922
    # write_configs_file_from_diagrams
 
923
    #===========================================================================
 
924
    def write_configs_file_from_diagrams(self, writer, configs, mapconfigs,
 
925
                                         nexternal, ninitial, model):
 
926
        """Write the actual configs.inc file.
 
927
        
 
928
        configs is the diagrams corresponding to configs (each
 
929
        diagrams is a list of corresponding diagrams for all
 
930
        subprocesses, with None if there is no corresponding diagrams
 
931
        for a given process).
 
932
        mapconfigs gives the diagram number for each config.
 
933
 
 
934
        For s-channels, we need to output one PDG for each subprocess in
 
935
        the subprocess group, in order to be able to pick the right
 
936
        one for multiprocesses."""
 
937
 
 
938
        lines = []
 
939
 
 
940
        s_and_t_channels = []
 
941
 
 
942
        minvert = min([max([d for d in config if d][0].get_vertex_leg_numbers()) \
 
943
                       for config in configs])
 
944
 
 
945
        # Number of subprocesses
 
946
        nsubprocs = len(configs[0])
 
947
 
 
948
        nconfigs = 0
 
949
 
 
950
        new_pdg = model.get_first_non_pdg()
 
951
 
 
952
        for iconfig, helas_diags in enumerate(configs):
 
953
            if any([vert > minvert for vert in
 
954
                    [d for d in helas_diags if d][0].get_vertex_leg_numbers()]):
 
955
                # Only 3-vertices allowed in configs.inc
 
956
                continue
 
957
            nconfigs += 1
 
958
 
 
959
            # Need s- and t-channels for all subprocesses, including
 
960
            # those that don't contribute to this config
 
961
            empty_verts = []
 
962
            stchannels = []
 
963
            for h in helas_diags:
 
964
                if h:
 
965
                    # get_s_and_t_channels gives vertices starting from
 
966
                    # final state external particles and working inwards
 
967
                    stchannels.append(h.get('amplitudes')[0].\
 
968
                                      get_s_and_t_channels(ninitial, model, new_pdg))
 
969
                else:
 
970
                    stchannels.append((empty_verts, None))
 
971
 
 
972
            # For t-channels, just need the first non-empty one
 
973
            tchannels = [t for s,t in stchannels if t != None][0]
 
974
 
 
975
            # For s_and_t_channels (to be used later) use only first config
 
976
            s_and_t_channels.append([[s for s,t in stchannels if t != None][0],
 
977
                                     tchannels])
 
978
 
 
979
            # Make sure empty_verts is same length as real vertices
 
980
            if any([s for s,t in stchannels]):
 
981
                empty_verts[:] = [None]*max([len(s) for s,t in stchannels])
 
982
 
 
983
                # Reorganize s-channel vertices to get a list of all
 
984
                # subprocesses for each vertex
 
985
                schannels = zip(*[s for s,t in stchannels])
 
986
            else:
 
987
                schannels = []
 
988
 
 
989
            allchannels = schannels
 
990
            if len(tchannels) > 1:
 
991
                # Write out tchannels only if there are any non-trivial ones
 
992
                allchannels = schannels + tchannels
 
993
 
 
994
            # Write out propagators for s-channel and t-channel vertices
 
995
 
 
996
            lines.append("# Diagram %d" % (mapconfigs[iconfig]))
 
997
            # Correspondance between the config and the diagram = amp2
 
998
            lines.append("data mapconfig(%d)/%d/" % (nconfigs,
 
999
                                                     mapconfigs[iconfig]))
 
1000
 
 
1001
            for verts in allchannels:
 
1002
                if verts in schannels:
 
1003
                    vert = [v for v in verts if v][0]
 
1004
                else:
 
1005
                    vert = verts
 
1006
                daughters = [leg.get('number') for leg in vert.get('legs')[:-1]]
 
1007
                last_leg = vert.get('legs')[-1]
 
1008
                lines.append("data (iforest(i,%d,%d),i=1,%d)/%s/" % \
 
1009
                             (last_leg.get('number'), nconfigs, len(daughters),
 
1010
                              ",".join([str(d) for d in daughters])))
 
1011
                if verts in schannels:
 
1012
                    pdgs = []
 
1013
                    for v in verts:
 
1014
                        if v:
 
1015
                            pdgs.append(v.get('legs')[-1].get('id'))
 
1016
                        else:
 
1017
                            pdgs.append(0)
 
1018
                    lines.append("data (sprop(i,%d,%d),i=1,%d)/%s/" % \
 
1019
                                 (last_leg.get('number'), nconfigs, nsubprocs,
 
1020
                                  ",".join([str(d) for d in pdgs])))
 
1021
                    lines.append("data tprid(%d,%d)/0/" % \
 
1022
                                 (last_leg.get('number'), nconfigs))
 
1023
                elif verts in tchannels[:-1]:
 
1024
                    lines.append("data tprid(%d,%d)/%d/" % \
 
1025
                                 (last_leg.get('number'), nconfigs,
 
1026
                                  abs(last_leg.get('id'))))
 
1027
                    lines.append("data (sprop(i,%d,%d),i=1,%d)/%s/" % \
 
1028
                                 (last_leg.get('number'), nconfigs, nsubprocs,
 
1029
                                  ",".join(['0'] * nsubprocs)))
 
1030
 
 
1031
        # Write out number of configs
 
1032
        lines.append("# Number of configs")
 
1033
        lines.append("data mapconfig(0)/%d/" % nconfigs)
 
1034
 
 
1035
        # Write the file
 
1036
        writer.writelines(lines)
 
1037
 
 
1038
        return s_and_t_channels
730
1039
 
731
1040
    #===========================================================================
732
1041
    # Global helper methods
748
1057
            else:
749
1058
                return '-'
750
1059
 
751
 
        res_str = '%+i' % total_coeff.numerator
 
1060
        res_str = '%+iD0' % total_coeff.numerator
752
1061
 
753
1062
        if total_coeff.denominator != 1:
754
1063
            # Check if total_coeff is an integer
755
 
            res_str = res_str + './%i.' % total_coeff.denominator
 
1064
            res_str = res_str + '/%iD0' % total_coeff.denominator
756
1065
 
757
1066
        if is_imaginary:
758
1067
            res_str = res_str + '*imag1'
759
1068
 
760
1069
        return res_str + '*'
761
1070
 
762
 
    def set_compiler(self, default_compiler):
 
1071
    def set_compiler(self, default_compiler, force=False):
763
1072
        """Set compiler based on what's available on the system"""
764
 
        
765
 
        
 
1073
                
766
1074
        # Check for compiler
767
1075
        if default_compiler and misc.which(default_compiler):
768
1076
            compiler = default_compiler
781
1089
        self.replace_make_opt_compiler(compiler)
782
1090
        # Replace also for Template but not for cluster
783
1091
        if not os.environ.has_key('MADGRAPH_DATA') and ReadWrite:
784
 
            self.replace_make_opt_compiler(compiler, pjoin(MG5DIR, 'Template'))
 
1092
            self.replace_make_opt_compiler(compiler, pjoin(MG5DIR, 'Template', 'LO'))
785
1093
 
786
1094
    def replace_make_opt_compiler(self, compiler, root_dir = ""):
787
1095
        """Set FC=compiler in Source/make_opts"""
815
1123
    """Class to take care of exporting a set of matrix elements to
816
1124
    MadGraph v4 StandAlone format."""
817
1125
 
 
1126
    def __init__(self, *args, **opts):
 
1127
        """add the format information compare to standard init"""
 
1128
        
 
1129
        if 'format' in opts:
 
1130
            self.format = opts['format']
 
1131
            del opts['format']
 
1132
        else:
 
1133
            self.format = 'standalone'
 
1134
        ProcessExporterFortran.__init__(self, *args, **opts)
 
1135
 
818
1136
    def copy_v4template(self, modelname):
819
1137
        """Additional actions needed for setup of Template
820
1138
        """
821
1139
 
822
 
        
823
1140
        #First copy the full template tree if dir_path doesn't exit
824
1141
        if os.path.isdir(self.dir_path):
825
1142
            return
826
1143
        
827
1144
        logger.info('initialize a new standalone directory: %s' % \
828
1145
                        os.path.basename(self.dir_path))
829
 
        temp_dir = pjoin(self.mgme_dir, 'Template')
830
 
        
 
1146
        temp_dir = pjoin(self.mgme_dir, 'Template/LO')
831
1147
        
832
1148
        # Create the directory structure
833
1149
        os.mkdir(self.dir_path)
857
1173
        # Add file in SubProcesses
858
1174
        shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'makefile_sa_f_sp'), 
859
1175
                    pjoin(self.dir_path, 'SubProcesses', 'makefile'))
860
 
        shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'check_sa.f'), 
861
 
                    pjoin(self.dir_path, 'SubProcesses', 'check_sa.f'))
862
1176
        
 
1177
        if self.format == 'standalone':
 
1178
            shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'check_sa.f'), 
 
1179
                    pjoin(self.dir_path, 'SubProcesses', 'check_sa.f'))
 
1180
        elif self.format == 'standalone_rw':
 
1181
            shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'driver_reweight.f'), 
 
1182
                    pjoin(self.dir_path, 'SubProcesses', 'check_sa.f'))
 
1183
                        
863
1184
        # Add file in Source
864
1185
        shutil.copy(pjoin(temp_dir, 'Source', 'make_opts'), 
865
1186
                    pjoin(self.dir_path, 'Source'))        
866
1187
        # add the makefile 
867
1188
        filename = pjoin(self.dir_path,'Source','makefile')
868
 
        self.write_source_makefile(writers.FileWriter(filename))            
 
1189
        self.write_source_makefile(writers.FileWriter(filename))          
869
1190
        
870
1191
    #===========================================================================
871
1192
    # export model files
907
1228
                              online = False, compiler='g77'):
908
1229
        """Finalize Standalone MG4 directory by generation proc_card_mg5.dat"""
909
1230
 
910
 
        self.set_compiler(compiler)
 
1231
        self.compiler_choice(compiler)
911
1232
        self.make()
912
1233
 
913
1234
        # Write command history as proc_card_mg5
918
1239
            output_file.write(text)
919
1240
            output_file.close()
920
1241
 
 
1242
    def compiler_choice(self, compiler):
 
1243
        """ Different daughter classes might want different compilers.
 
1244
        So this function is meant to be overloaded if desired."""
 
1245
        
 
1246
        self.set_compiler(compiler)
 
1247
 
921
1248
    #===========================================================================
922
1249
    # generate_subprocess_directory_v4
923
1250
    #===========================================================================
932
1259
        dirpath = pjoin(self.dir_path, 'SubProcesses', \
933
1260
                       "P%s" % matrix_element.get('processes')[0].shell_string())
934
1261
 
 
1262
        if self.opt['sa_symmetry']:
 
1263
            # avoid symmetric output
 
1264
            for proc in matrix_element.get('processes'):
 
1265
                leg0 = proc.get('legs')[0]
 
1266
                leg1 = proc.get('legs')[1]
 
1267
                if not leg1.get('state'):
 
1268
                    proc.get('legs')[0] = leg1
 
1269
                    proc.get('legs')[1] = leg0
 
1270
                    dirpath2 =  pjoin(self.dir_path, 'SubProcesses', \
 
1271
                           "P%s" % proc.shell_string())
 
1272
                    #restore original order
 
1273
                    proc.get('legs')[1] = leg1
 
1274
                    proc.get('legs')[0] = leg0                
 
1275
                    if os.path.exists(dirpath2):
 
1276
                        logger.info('Symmetric directory exists')
 
1277
                        return 0
 
1278
 
935
1279
        try:
936
1280
            os.mkdir(dirpath)
937
1281
        except os.error as error:
949
1293
        (nexternal, ninitial) = matrix_element.get_nexternal_ninitial()
950
1294
 
951
1295
        # Create the matrix.f file and the nexternal.inc file
952
 
        filename = 'matrix.f'
 
1296
        if self.opt['export_format']=='standalone_msP':
 
1297
            filename = 'matrix_prod.f'
 
1298
        else:
 
1299
            filename = 'matrix.f'
953
1300
        calls = self.write_matrix_element_v4(
954
1301
            writers.FortranWriter(filename),
955
1302
            matrix_element,
956
1303
            fortran_model)
957
1304
 
 
1305
        if self.opt['export_format']=='standalone_msP':
 
1306
            filename = 'configs_production.inc'
 
1307
            mapconfigs, s_and_t_channels = self.write_configs_file(\
 
1308
                writers.FortranWriter(filename),
 
1309
                matrix_element)
 
1310
 
 
1311
            filename = 'props_production.inc'
 
1312
            self.write_props_file(writers.FortranWriter(filename),
 
1313
                             matrix_element,
 
1314
                             s_and_t_channels)
 
1315
 
 
1316
            filename = 'nexternal_prod.inc'
 
1317
            self.write_nexternal_madspin(writers.FortranWriter(filename),
 
1318
                             nexternal, ninitial)
 
1319
 
958
1320
        filename = 'nexternal.inc'
959
1321
        self.write_nexternal_file(writers.FortranWriter(filename),
960
1322
                             nexternal, ninitial)
981
1343
 
982
1344
        linkfiles = ['check_sa.f', 'coupl.inc', 'makefile']
983
1345
 
984
 
 
985
1346
        for file in linkfiles:
986
1347
            ln('../%s' % file)
987
1348
 
1022
1383
        if not isinstance(writer, writers.FortranWriter):
1023
1384
            raise writers.FortranWriter.FortranWriterError(\
1024
1385
                "writer not FortranWriter")
 
1386
            
 
1387
        if not self.opt.has_key('sa_symmetry'):
 
1388
            self.opt['sa_symmetry']=False
1025
1389
 
1026
1390
        # Set lowercase/uppercase Fortran code
1027
1391
        writers.FortranWriter.downcase = False
1028
1392
 
1029
 
        replace_dict = {}
 
1393
        replace_dict = {'global_variable':'', 'amp2_lines':''}
1030
1394
 
1031
1395
        # Extract helas calls
1032
1396
        helas_calls = fortran_model.get_matrix_element_calls(\
1055
1419
 
1056
1420
        # Extract overall denominator
1057
1421
        # Averaging initial state color, spin, and identical FS particles
1058
 
        den_factor_line = self.get_den_factor_line(matrix_element)
1059
 
        replace_dict['den_factor_line'] = den_factor_line
 
1422
        replace_dict['den_factor_line'] = self.get_den_factor_line(matrix_element)
1060
1423
 
1061
1424
        # Extract ngraphs
1062
1425
        ngraphs = matrix_element.get_number_of_amplitudes()
1074
1437
        color_data_lines = self.get_color_data_lines(matrix_element)
1075
1438
        replace_dict['color_data_lines'] = "\n".join(color_data_lines)
1076
1439
 
 
1440
        if self.opt['export_format']=='standalone_msP':
 
1441
        # For MadSpin need to return the AMP2
 
1442
            amp2_lines = self.get_amp2_lines(matrix_element, [] )
 
1443
            replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
 
1444
            replace_dict['global_variable'] = "       Double Precision amp2(NGRAPHS)\n       common/to_amps/  amp2\n"
1077
1445
 
1078
1446
        # Extract JAMP lines
1079
1447
        jamp_lines = self.get_JAMP_lines(matrix_element)
1080
1448
        replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
1081
1449
 
1082
 
        file = open(pjoin(_file_path, \
 
1450
        if self.opt['export_format']=='standalone_msP' :
 
1451
            file = open(pjoin(_file_path, \
 
1452
                          'iolibs/template_files/matrix_standalone_msP_v4.inc')).read()
 
1453
 
 
1454
        elif self.opt['export_format']=='standalone_msF':
 
1455
            file = open(pjoin(_file_path, \
 
1456
                          'iolibs/template_files/matrix_standalone_msF_v4.inc')).read()
 
1457
 
 
1458
        else:
 
1459
            file = open(pjoin(_file_path, \
1083
1460
                          'iolibs/template_files/matrix_standalone_v4.inc')).read()
 
1461
 
1084
1462
        file = file % replace_dict
1085
1463
 
1086
1464
        # Write the file
1109
1487
        
1110
1488
        # The next file are model dependant (due to SLAH convention)
1111
1489
        self.model_name = modelname
1112
 
        # Add the combine_events.f 
1113
 
        filename = pjoin(self.dir_path,'Source','combine_events.f')
1114
 
        self.write_combine_events(writers.FortranWriter(filename)) # already formatted
1115
1490
        # Add the symmetry.f 
1116
1491
        filename = pjoin(self.dir_path,'SubProcesses','symmetry.f')
1117
1492
        self.write_symmetry(writers.FortranWriter(filename))
1139
1514
                            self.dir_path+'/bin/internal/madevent_interface.py')
1140
1515
        cp(_file_path+'/interface/extended_cmd.py',
1141
1516
                                  self.dir_path+'/bin/internal/extended_cmd.py')
 
1517
        cp(_file_path+'/interface/common_run_interface.py',
 
1518
                            self.dir_path+'/bin/internal/common_run_interface.py')
1142
1519
        cp(_file_path+'/various/misc.py', self.dir_path+'/bin/internal/misc.py')        
1143
1520
        cp(_file_path+'/iolibs/files.py', self.dir_path+'/bin/internal/files.py')
1144
1521
        cp(_file_path+'/iolibs/save_load_object.py', 
1168
1545
                                 self.dir_path+'/bin/internal/me5_logging.conf') 
1169
1546
        cp(_file_path+'/interface/coloring_logging.py', 
1170
1547
                                 self.dir_path+'/bin/internal/coloring_logging.py')
 
1548
        # shower card and FO_analyse_card. 
 
1549
        #  Although not needed, it is imported by banner.py
 
1550
        cp(_file_path+'/various/shower_card.py', 
 
1551
                                 self.dir_path+'/bin/internal/shower_card.py') 
 
1552
        cp(_file_path+'/various/FO_analyse_card.py', 
 
1553
                                 self.dir_path+'/bin/internal/FO_analyse_card.py') 
1171
1554
 
1172
1555
 
1173
1556
    def convert_model_to_mg4(self, model, wanted_lorentz = [], 
1174
1557
                                                         wanted_couplings = []):
1175
1558
         
1176
 
         super(ProcessExporterFortranME,self).convert_model_to_mg4(model, 
 
1559
        super(ProcessExporterFortranME,self).convert_model_to_mg4(model, 
1177
1560
                                               wanted_lorentz, wanted_couplings)
1178
1561
         
1179
 
         IGNORE_PATTERNS = ('*.pyc','*.dat','*.py~')
1180
 
         try:
1181
 
             shutil.rmtree(pjoin(self.dir_path,'bin','internal','ufomodel'))
1182
 
         except OSError as error:
1183
 
             pass
1184
 
         shutil.copytree(model.get('version_tag').split('##')[0], 
 
1562
        IGNORE_PATTERNS = ('*.pyc','*.dat','*.py~')
 
1563
        try:
 
1564
            shutil.rmtree(pjoin(self.dir_path,'bin','internal','ufomodel'))
 
1565
        except OSError as error:
 
1566
            pass
 
1567
        model_path = model.get('modelpath')
 
1568
        # This is not safe if there is a '##' or '-' in the path.
 
1569
        shutil.copytree(model_path, 
1185
1570
                               pjoin(self.dir_path,'bin','internal','ufomodel'),
1186
1571
                               ignore=shutil.ignore_patterns(*IGNORE_PATTERNS))
1187
 
         if hasattr(model, 'restrict_card'):
1188
 
             out_path = pjoin(self.dir_path, 'bin', 'internal','ufomodel',
 
1572
        if hasattr(model, 'restrict_card'):
 
1573
            out_path = pjoin(self.dir_path, 'bin', 'internal','ufomodel',
1189
1574
                                                         'restrict_default.dat')
1190
 
             if isinstance(model.restrict_card, check_param_card.ParamCard):
1191
 
                 model.restrict_card.write(out_path)
1192
 
             else:
1193
 
                 files.cp(model.restrict_card, out_path)
 
1575
            if isinstance(model.restrict_card, check_param_card.ParamCard):
 
1576
                model.restrict_card.write(out_path)
 
1577
            else:
 
1578
                files.cp(model.restrict_card, out_path)
1194
1579
 
1195
1580
                
1196
1581
    #===========================================================================
1280
1665
                             matrix_element)
1281
1666
 
1282
1667
        filename = 'configs.inc'
1283
 
        mapconfigs, s_and_t_channels = self.write_configs_file(\
 
1668
        mapconfigs, (s_and_t_channels, nqcd_list) = self.write_configs_file(\
1284
1669
            writers.FortranWriter(filename),
1285
1670
            matrix_element)
1286
1671
 
 
1672
        filename = 'config_nqcd.inc'
 
1673
        self.write_config_nqcd_file(writers.FortranWriter(filename),
 
1674
                               nqcd_list)
 
1675
 
1287
1676
        filename = 'config_subproc_map.inc'
1288
1677
        self.write_config_subproc_map_file(writers.FortranWriter(filename),
1289
1678
                                           s_and_t_channels)
1388
1777
                     'run.inc',
1389
1778
                     'maxconfigs.inc',
1390
1779
                     'maxparticles.inc',
 
1780
                     'run_config.inc',
1391
1781
                     'setcuts.f',
1392
1782
                     'setscales.f',
1393
1783
                     'sudakov.inc',
1434
1824
            check_param_card.convert_to_mg5card(param_card, mg5_param)
1435
1825
            check_param_card.check_valid_param_card(mg5_param)
1436
1826
 
1437
 
 
 
1827
        # Add the combine_events.f modify param_card path/number of @X
 
1828
        filename = pjoin(self.dir_path,'Source','combine_events.f')
 
1829
        try:
 
1830
            nb_proc =[p.get('id') for me in matrix_elements for m in me.get('matrix_elements') for p in m.get('processes')]
 
1831
        except AttributeError:
 
1832
            nb_proc =[p.get('id') for m in matrix_elements.get('matrix_elements') for p in m.get('processes')]
 
1833
        nb_proc = len(set(nb_proc))
 
1834
        self.write_combine_events(writers.FortranWriter(filename), nb_proc) # already formatted
1438
1835
        # Write maxconfigs.inc based on max of ME's/subprocess groups
1439
1836
        filename = pjoin(self.dir_path,'Source','maxconfigs.inc')
1440
1837
        self.write_maxconfigs_file(writers.FortranWriter(filename),
1461
1858
        if makejpg:
1462
1859
            logger.info("Generate jpeg diagrams")
1463
1860
            for Pdir in P_dir_list:
1464
 
                os.chdir(Pdir)
1465
 
                try:
1466
 
                    misc.call([pjoin(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')],
1467
 
                                stdout = devnull)
1468
 
                except:
1469
 
                    os.system('chmod +x %s ' % pjoin(old_pos, self.dir_path, 'bin', 'internal', '*'))
1470
 
                    misc.call([pjoin(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')],
1471
 
                                stdout = devnull)
1472
 
                os.chdir(os.path.pardir)
 
1861
                misc.call([pjoin(old_pos, self.dir_path, 'bin', 'internal', 'gen_jpeg-pl')],
 
1862
                                stdout = devnull, cwd=Pdir)
1473
1863
 
1474
1864
        logger.info("Generate web pages")
1475
1865
        # Create the WebPage using perl script
1720
2110
        return True
1721
2111
 
1722
2112
    #===========================================================================
1723
 
    # write_coloramps_file
 
2113
    # write_colors_file
1724
2114
    #===========================================================================
1725
2115
    def write_colors_file(self, writer, matrix_elements):
1726
2116
        """Write the get_color.f file for MadEvent, which returns color
1738
2128
                               for d in me.get('diagrams')], []) \
1739
2129
                          for me in matrix_elements], []))
1740
2130
 
1741
 
        leg_ids = set(sum([sum([[l.get('id') for l in \
1742
 
                                 p.get_legs_with_decays()] for p in \
1743
 
                                me.get('processes')], []) for me in
1744
 
                           matrix_elements], []))
 
2131
        leg_ids = set(sum([sum([sum([[l.get('id'), 
 
2132
                          model.get_particle(l.get('id')).get_anti_pdg_code()] \
 
2133
                                  for l in p.get_legs_with_decays()], []) \
 
2134
                                for p in me.get('processes')], []) \
 
2135
                           for me in matrix_elements], []))
1745
2136
        particle_ids = sorted(list(wf_ids.union(leg_ids)))
1746
2137
 
1747
2138
        lines = """function get_color(ipdg)
1779
2170
        return True
1780
2171
 
1781
2172
    #===========================================================================
 
2173
    # write_config_nqcd_file
 
2174
    #===========================================================================
 
2175
    def write_config_nqcd_file(self, writer, nqcd_list):
 
2176
        """Write the config_nqcd.inc with the number of QCD couplings
 
2177
        for each config"""
 
2178
 
 
2179
        lines = []
 
2180
        for iconf, n in enumerate(nqcd_list):
 
2181
            lines.append("data nqcd(%d)/%d/" % (iconf+1, n))
 
2182
 
 
2183
        # Write the file
 
2184
        writer.writelines(lines)
 
2185
 
 
2186
        return True
 
2187
 
 
2188
    #===========================================================================
1782
2189
    # write_maxconfigs_file
1783
2190
    #===========================================================================
1784
2191
    def write_maxconfigs_file(self, writer, matrix_elements):
1839
2246
                                                            nexternal, ninitial,
1840
2247
                                                            model)
1841
2248
 
1842
 
 
1843
 
 
1844
2249
    #===========================================================================
1845
2250
    # write_run_configs_file
1846
2251
    #===========================================================================
1874
2279
 
1875
2280
        s_and_t_channels = []
1876
2281
 
 
2282
        nqcd_list = []
 
2283
 
1877
2284
        minvert = min([max([d for d in config if d][0].get_vertex_leg_numbers()) \
1878
2285
                       for config in configs])
1879
2286
 
1933
2340
            # Correspondance between the config and the diagram = amp2
1934
2341
            lines.append("data mapconfig(%d)/%d/" % (nconfigs,
1935
2342
                                                     mapconfigs[iconfig]))
 
2343
            # Number of QCD couplings in this diagram
 
2344
            nqcd = 0
 
2345
            for h in helas_diags:
 
2346
                if h:
 
2347
                    try:
 
2348
                        nqcd = h.calculate_orders()['QCD']
 
2349
                    except KeyError:
 
2350
                        pass
 
2351
                    break
 
2352
                else:
 
2353
                    continue
 
2354
 
 
2355
            nqcd_list.append(nqcd)
1936
2356
 
1937
2357
            for verts in allchannels:
1938
2358
                if verts in schannels:
1971
2391
        # Write the file
1972
2392
        writer.writelines(lines)
1973
2393
 
1974
 
        return s_and_t_channels
 
2394
        return s_and_t_channels, nqcd_list
1975
2395
 
1976
2396
    #===========================================================================
1977
2397
    # write_config_subproc_map_file
2066
2486
    #===========================================================================
2067
2487
    # write_combine_events
2068
2488
    #===========================================================================
2069
 
    def write_combine_events(self, writer):
 
2489
    def write_combine_events(self, writer, nb_proc=100):
2070
2490
        """Write the SubProcess/driver.f file for MG4"""
2071
2491
 
2072
2492
        path = pjoin(_file_path,'iolibs','template_files','madevent_combine_events.f')
2075
2495
            card = 'Source/MODEL/MG5_param.dat'
2076
2496
        else:
2077
2497
            card = 'param_card.dat' 
2078
 
        text = open(path).read() % {'param_card_name':card} 
 
2498
        
 
2499
        #set maxpup (number of @X in the process card)
 
2500
            
 
2501
        text = open(path).read() % {'param_card_name':card, 'maxpup':nb_proc+1}
 
2502
        #the +1 is just a security. This is not needed but I feel(OM) safer with it.
2079
2503
        writer.write(text)
2080
2504
        
2081
2505
        return True
2271
2695
        return True
2272
2696
 
2273
2697
    #===========================================================================
2274
 
    # write_props_file
2275
 
    #===========================================================================
2276
 
    def write_props_file(self, writer, matrix_element, s_and_t_channels):
2277
 
        """Write the props.inc file for MadEvent. Needs input from
2278
 
        write_configs_file."""
2279
 
 
2280
 
        lines = []
2281
 
 
2282
 
        particle_dict = matrix_element.get('processes')[0].get('model').\
2283
 
                        get('particle_dict')
2284
 
 
2285
 
        for iconf, configs in enumerate(s_and_t_channels):
2286
 
            for vertex in configs[0] + configs[1][:-1]:
2287
 
                leg = vertex.get('legs')[-1]
2288
 
                if leg.get('id') not in particle_dict:
2289
 
                    # Fake propagator used in multiparticle vertices
2290
 
                    mass = 'zero'
2291
 
                    width = 'zero'
2292
 
                    pow_part = 0
2293
 
                else:
2294
 
                    particle = particle_dict[leg.get('id')]
2295
 
                    # Get mass
2296
 
                    if particle.get('mass').lower() == 'zero':
2297
 
                        mass = particle.get('mass')
2298
 
                    else:
2299
 
                        mass = "abs(%s)" % particle.get('mass')
2300
 
                    # Get width
2301
 
                    if particle.get('width').lower() == 'zero':
2302
 
                        width = particle.get('width')
2303
 
                    else:
2304
 
                        width = "abs(%s)" % particle.get('width')
2305
 
 
2306
 
                    pow_part = 1 + int(particle.is_boson())
2307
 
 
2308
 
                lines.append("prmass(%d,%d)  = %s" % \
2309
 
                             (leg.get('number'), iconf + 1, mass))
2310
 
                lines.append("prwidth(%d,%d) = %s" % \
2311
 
                             (leg.get('number'), iconf + 1, width))
2312
 
                lines.append("pow(%d,%d) = %d" % \
2313
 
                             (leg.get('number'), iconf + 1, pow_part))
2314
 
 
2315
 
        # Write the file
2316
 
        writer.writelines(lines)
2317
 
 
2318
 
        return True
2319
 
 
2320
 
    #===========================================================================
2321
2698
    # write_processes_file
2322
2699
    #===========================================================================
2323
2700
    def write_processes_file(self, writer, subproc_group):
2414
2791
 
2415
2792
        return True
2416
2793
 
2417
 
 
2418
 
 
2419
2794
#===============================================================================
2420
2795
# ProcessExporterFortranMEGroup
2421
2796
#===============================================================================
2535
2910
                                           subproc_diagrams_for_config)
2536
2911
 
2537
2912
        filename = 'configs.inc'
2538
 
        nconfigs, s_and_t_channels = self.write_configs_file(\
 
2913
        nconfigs, (s_and_t_channels, nqcd_list) = self.write_configs_file(\
2539
2914
            writers.FortranWriter(filename),
2540
2915
            subproc_group,
2541
2916
            subproc_diagrams_for_config)
2542
2917
 
 
2918
        filename = 'config_nqcd.inc'
 
2919
        self.write_config_nqcd_file(writers.FortranWriter(filename),
 
2920
                                    nqcd_list)
 
2921
 
2543
2922
        filename = 'decayBW.inc'
2544
2923
        self.write_decayBW_file(writers.FortranWriter(filename),
2545
2924
                           s_and_t_channels)
2634
3013
                     'run.inc',
2635
3014
                     'maxconfigs.inc',
2636
3015
                     'maxparticles.inc',
 
3016
                     'run_config.inc',
2637
3017
                     'setcuts.f',
2638
3018
                     'setscales.f',
2639
3019
                     'sudakov.inc',
2874
3254
 
2875
3255
class UFO_model_to_mg4(object):
2876
3256
    """ A converter of the UFO-MG5 Model to the MG4 format """
 
3257
 
 
3258
    # The list below shows the only variables the user is allowed to change by
 
3259
    # himself for each PS point. If he changes any other, then calling 
 
3260
    # UPDATE_AS_PARAM() (or equivalently MP_UPDATE_AS_PARAM()) will not
 
3261
    # correctly account for the change.
 
3262
    PS_dependent_key = ['aS','MU_R']
 
3263
    mp_complex_format = 'complex*32'
 
3264
    mp_real_format = 'real*16'
 
3265
    # Warning, it is crucial none of the couplings/parameters of the model
 
3266
    # starts with this prefix. I should add a check for this.
 
3267
    # You can change it as the global variable to check_param_card.ParamCard
 
3268
    mp_prefix = check_param_card.ParamCard.mp_prefix
2877
3269
    
2878
3270
    def __init__(self, model, output_path, opt=None):
2879
3271
        """ initialization of the objects """
2884
3276
        if opt:
2885
3277
            self.opt = opt
2886
3278
        else:
2887
 
            self.opt = {'complex_mass': False, 'export_format': 'madevent'}
 
3279
            self.opt = {'complex_mass': False, 'export_format': 'madevent', 'mp':True}
2888
3280
            
2889
3281
        self.coups_dep = []    # (name, expression, type)
2890
3282
        self.coups_indep = []  # (name, expression, type)
2892
3284
        self.params_indep = [] # (name, expression, type)
2893
3285
        self.params_ext = []   # external parameter
2894
3286
        self.p_to_f = parsers.UFOExpressionParserFortran()
 
3287
        self.mp_p_to_f = parsers.UFOExpressionParserMPFortran()            
2895
3288
    
2896
3289
    def pass_parameter_to_case_insensitive(self):
2897
3290
        """modify the parameter if some of them are identical up to the case"""
2959
3352
 
2960
3353
            if key == ('external',):
2961
3354
                self.params_ext += to_add
2962
 
            elif 'aS' in key:
 
3355
            elif any([(k in key) for k in self.PS_dependent_key]):
2963
3356
                self.params_dep += to_add
2964
3357
            else:
2965
3358
                self.params_indep += to_add
2967
3360
        keys = self.model['couplings'].keys()
2968
3361
        keys.sort(key=len)
2969
3362
        for key, coup_list in self.model['couplings'].items():
2970
 
            if 'aS' in key:
2971
 
                self.coups_dep += [c for c in coup_list if 
 
3363
            if any([(k in key) for k in self.PS_dependent_key]):
 
3364
                self.coups_dep += [c for c in coup_list if
2972
3365
                                   (not wanted_couplings or c.name in \
2973
3366
                                    wanted_couplings)]
2974
3367
            else:
3024
3417
        
3025
3418
        #write the definition of the parameter
3026
3419
        self.create_input()
3027
 
        self.create_intparam_def()
 
3420
        self.create_intparam_def(dp=True,mp=False)
 
3421
        if self.opt['mp']:
 
3422
            self.create_intparam_def(dp=False,mp=True)
3028
3423
        
3029
3424
        
3030
3425
        # definition of the coupling.
 
3426
        self.create_actualize_mp_ext_param_inc()
3031
3427
        self.create_coupl_inc()
3032
3428
        self.create_write_couplings()
3033
3429
        self.create_couplings()
3035
3431
        # the makefile
3036
3432
        self.create_makeinc()
3037
3433
        self.create_param_write()
 
3434
 
 
3435
        # The model functions
 
3436
        self.create_model_functions_inc()
 
3437
        self.create_model_functions_def()
3038
3438
        
3039
3439
        # The param_card.dat        
3040
3440
        self.create_param_card()
3052
3452
    
3053
3453
        
3054
3454
        #copy the library files
3055
 
        file_to_link = ['formats.inc', 'lha_read.f','printout.f', \
 
3455
        file_to_link = ['formats.inc','printout.f', \
3056
3456
                        'rw_para.f', 'testprog.f']
3057
3457
    
3058
 
    
3059
3458
        for filename in file_to_link:
3060
 
            cp( MG5DIR + '/models/template_files/fortran/' + filename, self.dir_path)
3061
 
 
3062
 
        if self.opt['export_format'] == 'madevent':
 
3459
            cp( MG5DIR + '/models/template_files/fortran/' + filename, \
 
3460
                                                                self.dir_path)
 
3461
            
 
3462
        file = open(os.path.join(MG5DIR,\
 
3463
                              'models/template_files/fortran/rw_para.f')).read()
 
3464
 
 
3465
        includes=["include \'coupl.inc\'","include \'input.inc\'"]
 
3466
        if self.opt['mp']:
 
3467
            includes.extend(["include \'mp_coupl.inc\'","include \'mp_input.inc\'"])
 
3468
        # In standalone and madloop we do no use the compiled param card but
 
3469
        # still parse the .dat one so we must load it.
 
3470
        if self.opt['export_format'] in ['madloop','madloop_optimized']:
 
3471
            load_card = 'call LHA_loadcard(param_name,npara,param,value)'
 
3472
            lha_read_filename='lha_read_mp.f'
 
3473
        elif self.opt['export_format'].startswith('standalone'):
 
3474
            load_card = 'call LHA_loadcard(param_name,npara,param,value)'
 
3475
            lha_read_filename='lha_read.f'
 
3476
        else:
 
3477
            load_card = ''
 
3478
            lha_read_filename='lha_read.f'
 
3479
        cp( MG5DIR + '/models/template_files/fortran/' + lha_read_filename, \
 
3480
                                       os.path.join(self.dir_path,'lha_read.f'))
 
3481
        
 
3482
        file=file%{'includes':'\n      '.join(includes),
 
3483
                   'load_card':load_card}
 
3484
        writer=open(os.path.join(self.dir_path,'rw_para.f'),'w')
 
3485
        writer.writelines(file)
 
3486
        writer.close()
 
3487
 
 
3488
        if self.opt['export_format'] in ['madevent', 'FKS5_default', 'FKS5_optimized']:
3063
3489
            cp( MG5DIR + '/models/template_files/fortran/makefile_madevent', 
3064
3490
                self.dir_path + '/makefile')
3065
 
        else:
 
3491
            if self.opt['export_format'] in ['FKS5_default', 'FKS5_optimized']:
 
3492
                path = pjoin(self.dir_path, 'makefile')
 
3493
                text = open(path).read()
 
3494
                text = text.replace('madevent','aMCatNLO')
 
3495
                open(path, 'w').writelines(text)
 
3496
 
 
3497
        elif self.opt['export_format'] in ['standalone', 'standalone_msP','standalone_msF',
 
3498
                                  'madloop','madloop_optimized', 'standalone_rw']:
3066
3499
            cp( MG5DIR + '/models/template_files/fortran/makefile_standalone', 
3067
3500
                self.dir_path + '/makefile')
3068
 
            text = open(pjoin(self.dir_path, 'rw_para.f')).read()
3069
 
            text = re.sub(r'c\s*call LHA_loadcard','       call LHA_loadcard',text, re.I)
3070
 
            fsock = open(pjoin(self.dir_path, 'rw_para.f'), 'w')
3071
 
            fsock.write(text)
3072
 
            fsock.close()
3073
 
            
3074
 
     
3075
 
 
 
3501
        else:
 
3502
            raise MadGraph5Error('Unknown format')
3076
3503
 
3077
3504
    def create_coupl_inc(self):
3078
3505
        """ write coupling.inc """
3079
3506
        
3080
3507
        fsock = self.open('coupl.inc', format='fortran')
3081
 
        
 
3508
        if self.opt['mp']:
 
3509
            mp_fsock = self.open('mp_coupl.inc', format='fortran')
 
3510
            mp_fsock_same_name = self.open('mp_coupl_same_name.inc',\
 
3511
                                            format='fortran')
 
3512
 
3082
3513
        # Write header
3083
3514
        header = """double precision G
3084
3515
                common/strong/ G
3087
3518
                common/weak/ gal
3088
3519
 
3089
3520
                """        
 
3521
        if self.model.get('expansion_order'):
 
3522
            header=header+"""double precision MU_R
 
3523
                common/rscale/ MU_R
 
3524
 
 
3525
                """
 
3526
        header = header+"""double precision Nf
 
3527
                parameter(Nf=%d)
 
3528
                """ % self.model.get_nflav()
 
3529
                
3090
3530
        fsock.writelines(header)
3091
3531
        
 
3532
        if self.opt['mp']:
 
3533
            header = """%(real_mp_format)s %(mp_prefix)sG
 
3534
                    common/MP_strong/ %(mp_prefix)sG
 
3535
                     
 
3536
                    %(complex_mp_format)s %(mp_prefix)sgal(2)
 
3537
                    common/MP_weak/ %(mp_prefix)sgal
 
3538
    
 
3539
                    """        
 
3540
            if self.model.get('expansion_order'):
 
3541
                header=header+"""%(complex_mp_format)s %(mp_prefix)sMU_R
 
3542
                    common/MP_rscale/ %(mp_prefix)sMU_R
 
3543
    
 
3544
                    """            
 
3545
            mp_fsock.writelines(header%{'real_mp_format':self.mp_real_format,
 
3546
                                  'complex_mp_format':self.mp_complex_format,
 
3547
                                  'mp_prefix':self.mp_prefix})
 
3548
            mp_fsock_same_name.writelines(header%{'real_mp_format':self.mp_real_format,
 
3549
                                  'complex_mp_format':self.mp_complex_format,
 
3550
                                  'mp_prefix':''})
 
3551
 
3092
3552
        # Write the Mass definition/ common block
3093
3553
        masses = set()
3094
3554
        widths = set()
3107
3567
                widths.add(one_width)
3108
3568
                if self.opt['complex_mass'] and one_mass.lower() != 'zero':
3109
3569
                    complex_mass.add('CMASS_%s' % one_mass)
3110
 
        
3111
 
        fsock.writelines('double precision '+','.join(masses)+'\n')
3112
 
        fsock.writelines('common/masses/ '+','.join(masses)+'\n\n')
 
3570
            
 
3571
        if masses:
 
3572
            fsock.writelines('double precision '+','.join(masses)+'\n')
 
3573
            fsock.writelines('common/masses/ '+','.join(masses)+'\n\n')
 
3574
            if self.opt['mp']:
 
3575
                mp_fsock_same_name.writelines(self.mp_real_format+' '+\
 
3576
                                                          ','.join(masses)+'\n')
 
3577
                mp_fsock_same_name.writelines('common/MP_masses/ '+\
 
3578
                                                        ','.join(masses)+'\n\n')                
 
3579
                mp_fsock.writelines(self.mp_real_format+' '+','.join([\
 
3580
                                        self.mp_prefix+m for m in masses])+'\n')
 
3581
                mp_fsock.writelines('common/MP_masses/ '+\
 
3582
                            ','.join([self.mp_prefix+m for m in masses])+'\n\n')                
 
3583
 
3113
3584
        if widths:
3114
3585
            fsock.writelines('double precision '+','.join(widths)+'\n')
3115
3586
            fsock.writelines('common/widths/ '+','.join(widths)+'\n\n')
 
3587
            if self.opt['mp']:
 
3588
                mp_fsock_same_name.writelines(self.mp_real_format+' '+\
 
3589
                                                          ','.join(widths)+'\n')
 
3590
                mp_fsock_same_name.writelines('common/MP_widths/ '+\
 
3591
                                                        ','.join(widths)+'\n\n')                
 
3592
                mp_fsock.writelines(self.mp_real_format+' '+','.join([\
 
3593
                                        self.mp_prefix+w for w in widths])+'\n')
 
3594
                mp_fsock.writelines('common/MP_widths/ '+\
 
3595
                            ','.join([self.mp_prefix+w for w in widths])+'\n\n')
3116
3596
        
3117
3597
        # Write the Couplings
3118
3598
        coupling_list = [coupl.name for coupl in self.coups_dep + self.coups_indep]       
3119
3599
        fsock.writelines('double complex '+', '.join(coupling_list)+'\n')
3120
3600
        fsock.writelines('common/couplings/ '+', '.join(coupling_list)+'\n')
 
3601
        if self.opt['mp']:
 
3602
            mp_fsock_same_name.writelines(self.mp_complex_format+' '+\
 
3603
                                                   ','.join(coupling_list)+'\n')
 
3604
            mp_fsock_same_name.writelines('common/MP_couplings/ '+\
 
3605
                                                 ','.join(coupling_list)+'\n\n')                
 
3606
            mp_fsock.writelines(self.mp_complex_format+' '+','.join([\
 
3607
                                 self.mp_prefix+c for c in coupling_list])+'\n')
 
3608
            mp_fsock.writelines('common/MP_couplings/ '+\
 
3609
                     ','.join([self.mp_prefix+c for c in coupling_list])+'\n\n')            
3121
3610
        
3122
3611
        # Write complex mass for complex mass scheme (if activated)
3123
 
        if self.opt['complex_mass']:
 
3612
        if self.opt['complex_mass'] and complex_mass:
3124
3613
            fsock.writelines('double complex '+', '.join(complex_mass)+'\n')
3125
 
            fsock.writelines('common/complex_mass/ '+', '.join(complex_mass)+'\n')            
 
3614
            fsock.writelines('common/complex_mass/ '+', '.join(complex_mass)+'\n')
 
3615
            if self.opt['mp']:
 
3616
                mp_fsock_same_name.writelines(self.mp_complex_format+' '+\
 
3617
                                                    ','.join(complex_mass)+'\n')
 
3618
                mp_fsock_same_name.writelines('common/MP_complex_mass/ '+\
 
3619
                                                  ','.join(complex_mass)+'\n\n')                
 
3620
                mp_fsock.writelines(self.mp_complex_format+' '+','.join([\
 
3621
                                self.mp_prefix+cm for cm in complex_mass])+'\n')
 
3622
                mp_fsock.writelines('common/MP_complex_mass/ '+\
 
3623
                    ','.join([self.mp_prefix+cm for cm in complex_mass])+'\n\n')                       
3126
3624
        
3127
3625
    def create_write_couplings(self):
3128
3626
        """ write the file coupl_write.inc """
3144
3642
        """create input.inc containing the definition of the parameters"""
3145
3643
        
3146
3644
        fsock = self.open('input.inc', format='fortran')
3147
 
 
 
3645
        if self.opt['mp']:
 
3646
            mp_fsock = self.open('mp_input.inc', format='fortran')
 
3647
                    
3148
3648
        #find mass/ width since they are already define
3149
3649
        already_def = set()
3150
3650
        for particle in self.model.get('particles'):
3153
3653
            if self.opt['complex_mass']:
3154
3654
                already_def.add('cmass_%s' % particle.get('mass').lower())
3155
3655
 
3156
 
        is_valid = lambda name: name !='G' and name.lower() not in already_def
 
3656
        is_valid = lambda name: name!='G' and name!='MU_R' and name.lower() not in already_def
3157
3657
        
3158
3658
        real_parameters = [param.name for param in self.params_dep + 
3159
3659
                            self.params_indep if param.type == 'real'
3165
3665
        
3166
3666
        fsock.writelines('double precision '+','.join(real_parameters)+'\n')
3167
3667
        fsock.writelines('common/params_R/ '+','.join(real_parameters)+'\n\n')
 
3668
        if self.opt['mp']:
 
3669
            mp_fsock.writelines(self.mp_real_format+' '+','.join([\
 
3670
                              self.mp_prefix+p for p in real_parameters])+'\n')
 
3671
            mp_fsock.writelines('common/MP_params_R/ '+','.join([\
 
3672
                            self.mp_prefix+p for p in real_parameters])+'\n\n')        
3168
3673
        
3169
3674
        complex_parameters = [param.name for param in self.params_dep + 
3170
3675
                            self.params_indep if param.type == 'complex' and
3171
3676
                            is_valid(param.name)]
3172
3677
 
3173
 
 
3174
3678
        fsock.writelines('double complex '+','.join(complex_parameters)+'\n')
3175
 
        fsock.writelines('common/params_C/ '+','.join(complex_parameters)+'\n\n')                
3176
 
    
3177
 
    
3178
 
 
3179
 
    def create_intparam_def(self):
3180
 
        """ create intparam_definition.inc """
3181
 
 
3182
 
        fsock = self.open('intparam_definition.inc', format='fortran')
 
3679
        fsock.writelines('common/params_C/ '+','.join(complex_parameters)+'\n\n')
 
3680
        if self.opt['mp']:
 
3681
            mp_fsock.writelines(self.mp_complex_format+' '+','.join([\
 
3682
                            self.mp_prefix+p for p in complex_parameters])+'\n')
 
3683
            mp_fsock.writelines('common/MP_params_C/ '+','.join([\
 
3684
                          self.mp_prefix+p for p in complex_parameters])+'\n\n')                
 
3685
    
 
3686
    
 
3687
 
 
3688
    def create_intparam_def(self, dp=True, mp=False):
 
3689
        """ create intparam_definition.inc setting the internal parameters.
 
3690
        Output the double precision and/or the multiple precision parameters
 
3691
        depending on the parameters dp and mp. If mp only, then the file names
 
3692
        get the 'mp_' prefix.
 
3693
         """
 
3694
 
 
3695
        fsock = self.open('%sintparam_definition.inc'%
 
3696
                             ('mp_' if mp and not dp else ''), format='fortran')
3183
3697
        
3184
3698
        fsock.write_comments(\
3185
3699
                "Parameters that should not be recomputed event by event.\n")
3186
3700
        fsock.writelines("if(readlha) then\n")
3187
 
        fsock.writelines("G = 2 * DSQRT(AS*PI) ! for the first init\n")
 
3701
        if dp:        
 
3702
            fsock.writelines("G = 2 * DSQRT(AS*PI) ! for the first init\n")
 
3703
        if mp:
 
3704
            fsock.writelines("MP__G = 2 * SQRT(MP__AS*MP__PI) ! for the first init\n")
3188
3705
        for param in self.params_indep:
3189
3706
            if param.name == 'ZERO':
3190
3707
                continue
3191
 
            fsock.writelines("%s = %s\n" % (param.name,
 
3708
            if dp:
 
3709
                fsock.writelines("%s = %s\n" % (param.name,
3192
3710
                                            self.p_to_f.parse(param.expr)))
3193
 
        
 
3711
            if mp:
 
3712
                fsock.writelines("%s%s = %s\n" % (self.mp_prefix,param.name,
 
3713
                                            self.mp_p_to_f.parse(param.expr)))    
 
3714
 
3194
3715
        fsock.writelines('endif')
3195
3716
        
3196
3717
        fsock.write_comments('\nParameters that should be recomputed at an event by even basis.\n')
3197
 
        fsock.writelines("aS = G**2/4/pi\n")
 
3718
        if dp:        
 
3719
            fsock.writelines("aS = G**2/4/pi\n")
 
3720
        if mp:
 
3721
            fsock.writelines("MP__aS = MP__G**2/4/MP__PI\n")
3198
3722
        for param in self.params_dep:
3199
 
            fsock.writelines("%s = %s\n" % (param.name,
 
3723
            if dp:
 
3724
                fsock.writelines("%s = %s\n" % (param.name,
3200
3725
                                            self.p_to_f.parse(param.expr)))
 
3726
            elif mp:
 
3727
                fsock.writelines("%s%s = %s\n" % (self.mp_prefix,param.name,
 
3728
                                            self.mp_p_to_f.parse(param.expr)))
3201
3729
 
3202
3730
        fsock.write_comments("\nDefinition of the EW coupling used in the write out of aqed\n")
3203
3731
        if ('aEWM1',) in self.model['parameters']:
3204
 
            fsock.writelines(""" gal(1) = 3.5449077018110318 / DSQRT(aEWM1)
 
3732
            if dp:
 
3733
                fsock.writelines(""" gal(1) = 3.5449077018110318 / DSQRT(aEWM1)
3205
3734
                                 gal(2) = 1d0
3206
 
                         """)            
 
3735
                         """)
 
3736
            elif mp:
 
3737
                fsock.writelines(""" %(mp_prefix)sgal(1) = 2 * SQRT(MP__PI/MP__aEWM1)
 
3738
                                 %(mp_prefix)sgal(2) = 1d0 
 
3739
                                 """ %{'mp_prefix':self.mp_prefix})
 
3740
                pass
3207
3741
        else:
3208
 
            logger.warning('$RED aEWM1 not define in MODEL. AQED will not be written correcty in LHE FILE')
3209
 
            fsock.writelines(""" gal(1) = 1d0
3210
 
                             gal(2) = 1d0
3211
 
                         """)
 
3742
            if dp:
 
3743
                logger.warning('$RED aEWM1 not define in MODEL. AQED will not be written correcty in LHE FILE')
 
3744
                fsock.writelines(""" gal(1) = 1d0
 
3745
                                 gal(2) = 1d0
 
3746
                             """)
 
3747
            elif mp:
 
3748
                fsock.writelines(""" %(mp_prefix)sgal(1) = 1e0_16
 
3749
                                 %(mp_prefix)sgal(2) = 1e0_16
 
3750
                             """%{'mp_prefix':self.mp_prefix})
3212
3751
 
3213
3752
    
3214
3753
    def create_couplings(self):
3221
3760
        nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 
3222
3761
        
3223
3762
        for i in range(nb_coup_indep):
 
3763
            # For the independent couplings, we compute the double and multiple
 
3764
            # precision ones together
3224
3765
            data = self.coups_indep[nb_def_by_file * i: 
3225
3766
                             min(len(self.coups_indep), nb_def_by_file * (i+1))]
3226
 
            self.create_couplings_part(i + 1, data)
 
3767
            self.create_couplings_part(i + 1, data, dp=True, mp=self.opt['mp'])
3227
3768
            
3228
3769
        for i in range(nb_coup_dep):
 
3770
            # For the dependent couplings, we compute the double and multiple
 
3771
            # precision ones in separate subroutines.
3229
3772
            data = self.coups_dep[nb_def_by_file * i: 
3230
3773
                               min(len(self.coups_dep), nb_def_by_file * (i+1))]
3231
 
            self.create_couplings_part( i + 1 + nb_coup_indep , data)        
 
3774
            self.create_couplings_part( i + 1 + nb_coup_indep , data, 
 
3775
                                                               dp=True,mp=False)
 
3776
            if self.opt['mp']:
 
3777
                self.create_couplings_part( i + 1 + nb_coup_indep , data, 
 
3778
                                                              dp=False,mp=True)
3232
3779
        
3233
3780
        
3234
3781
    def create_couplings_main(self, nb_def_by_file=25):
3239
3786
        fsock.writelines("""subroutine coup()
3240
3787
 
3241
3788
                            implicit none
3242
 
                            double precision PI
 
3789
                            double precision PI, ZERO
3243
3790
                            logical READLHA
3244
3791
                            parameter  (PI=3.141592653589793d0)
3245
 
                            
3246
 
                            include \'input.inc\'
 
3792
                            parameter  (ZERO=0d0)""")
 
3793
        if self.opt['mp']:
 
3794
            fsock.writelines("""%s MP__PI, MP__ZERO
 
3795
                                parameter (MP__PI=3.1415926535897932384626433832795e0_16)
 
3796
                                parameter (MP__ZERO=0e0_16)
 
3797
                                include \'mp_input.inc\'
 
3798
                                include \'mp_coupl.inc\'
 
3799
                        """%self.mp_real_format) 
 
3800
        fsock.writelines("""include \'input.inc\'
3247
3801
                            include \'coupl.inc\'
3248
3802
                            READLHA = .true.
3249
 
                            include \'intparam_definition.inc\'\n\n
3250
 
                         """)
 
3803
                            include \'intparam_definition.inc\'""")
 
3804
        if self.opt['mp']:
 
3805
            fsock.writelines("""include \'mp_intparam_definition.inc\'\n""")
3251
3806
        
3252
3807
        nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 
3253
3808
        nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 
3260
3815
        fsock.writelines('\n'.join(\
3261
3816
                    ['call coup%s()' %  (nb_coup_indep + i + 1) \
3262
3817
                      for i in range(nb_coup_dep)]))
 
3818
        if self.opt['mp']:
 
3819
            fsock.writelines('\n'.join(\
 
3820
                    ['call mp_coup%s()' %  (nb_coup_indep + i + 1) \
 
3821
                      for i in range(nb_coup_dep)]))
3263
3822
        fsock.writelines('''\n return \n end\n''')
3264
3823
 
3265
3824
        fsock.writelines("""subroutine update_as_param()
3266
3825
 
3267
3826
                            implicit none
3268
 
                            double precision PI
 
3827
                            double precision PI, ZERO
3269
3828
                            logical READLHA
3270
 
                            parameter  (PI=3.141592653589793d0)
3271
 
                            
3272
 
                            include \'input.inc\'
 
3829
                            parameter  (PI=3.141592653589793d0)            
 
3830
                            parameter  (ZERO=0d0)""")
 
3831
        fsock.writelines("""include \'input.inc\'
3273
3832
                            include \'coupl.inc\'
3274
 
                            READLHA = .false.
3275
 
                            include \'intparam_definition.inc\'\n\n
 
3833
                            READLHA = .false.""")
 
3834
        fsock.writelines("""    
 
3835
                            include \'intparam_definition.inc\'\n
3276
3836
                         """)
3277
 
        
 
3837
            
3278
3838
        nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 
3279
3839
        nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 
3280
3840
                
3285
3845
                      for i in range(nb_coup_dep)]))
3286
3846
        fsock.writelines('''\n return \n end\n''')
3287
3847
 
3288
 
    def create_couplings_part(self, nb_file, data):
3289
 
        """ create couplings[nb_file].f containing information coming from data
 
3848
        if self.opt['mp']:
 
3849
            fsock.writelines("""subroutine mp_update_as_param()
 
3850
    
 
3851
                                implicit none
 
3852
                                logical READLHA""")
 
3853
            fsock.writelines("""%s MP__PI, MP__ZERO
 
3854
                                    parameter (MP__PI=3.1415926535897932384626433832795e0_16)
 
3855
                                    parameter (MP__ZERO=0e0_16)
 
3856
                                    include \'mp_input.inc\'
 
3857
                                    include \'mp_coupl.inc\'
 
3858
                            """%self.mp_real_format)
 
3859
            fsock.writelines("""include \'input.inc\'
 
3860
                                include \'coupl.inc\'
 
3861
                                include \'actualize_mp_ext_params.inc\'
 
3862
                                READLHA = .false.
 
3863
                                include \'mp_intparam_definition.inc\'\n
 
3864
                             """)
 
3865
            
 
3866
            nb_coup_indep = 1 + len(self.coups_indep) // nb_def_by_file 
 
3867
            nb_coup_dep = 1 + len(self.coups_dep) // nb_def_by_file 
 
3868
                    
 
3869
            fsock.write_comments('\ncouplings needed to be evaluated points by points\n')
 
3870
    
 
3871
            fsock.writelines('\n'.join(\
 
3872
                        ['call mp_coup%s()' %  (nb_coup_indep + i + 1) \
 
3873
                          for i in range(nb_coup_dep)]))
 
3874
            fsock.writelines('''\n return \n end\n''')
 
3875
 
 
3876
    def create_couplings_part(self, nb_file, data, dp=True, mp=False):
 
3877
        """ create couplings[nb_file].f containing information coming from data.
 
3878
        Outputs the computation of the double precision and/or the multiple
 
3879
        precision couplings depending on the parameters dp and mp.
 
3880
        If mp is True and dp is False, then the prefix 'MP_' is appended to the
 
3881
        filename and subroutine name.
3290
3882
        """
3291
3883
        
3292
 
        fsock = self.open('couplings%s.f' % nb_file, format='fortran')
3293
 
        fsock.writelines("""subroutine coup%s()
3294
 
        
3295
 
          implicit none
3296
 
          double precision PI
3297
 
          parameter  (PI=3.141592653589793d0)
3298
 
      
3299
 
      
3300
 
          include 'input.inc'
3301
 
          include 'coupl.inc'
3302
 
                        """ % nb_file)
3303
 
        
3304
 
        for coupling in data:            
3305
 
            fsock.writelines('%s = %s' % (coupling.name,
 
3884
        fsock = self.open('%scouplings%s.f' %('mp_' if mp and not dp else '',
 
3885
                                                     nb_file), format='fortran')
 
3886
        fsock.writelines("""subroutine %scoup%s()
 
3887
          
 
3888
          implicit none"""%('mp_' if mp and not dp else '',nb_file))
 
3889
        if dp:
 
3890
            fsock.writelines("""
 
3891
              double precision PI
 
3892
              parameter  (PI=3.141592653589793d0)
 
3893
              include 'input.inc'
 
3894
              include 'coupl.inc'""")
 
3895
        if mp:
 
3896
            fsock.writelines("""%s MP__PI
 
3897
                                parameter (MP__PI=3.1415926535897932384626433832795e0_16)
 
3898
                                include \'mp_input.inc\'
 
3899
                                include \'mp_coupl.inc\'
 
3900
                        """%self.mp_real_format) 
 
3901
        fsock.writelines("""
 
3902
          include 'model_functions.inc'""")
 
3903
        for coupling in data:
 
3904
            if dp:            
 
3905
                fsock.writelines('%s = %s' % (coupling.name,
3306
3906
                                          self.p_to_f.parse(coupling.expr)))
 
3907
            if mp:
 
3908
                fsock.writelines('%s%s = %s' % (self.mp_prefix,coupling.name,
 
3909
                                          self.mp_p_to_f.parse(coupling.expr)))
3307
3910
        fsock.writelines('end')
3308
3911
 
 
3912
    def create_model_functions_inc(self):
 
3913
        """ Create model_functions.inc which contains the various declarations
 
3914
        of auxiliary functions which might be used in the couplings expressions
 
3915
        """
 
3916
        
 
3917
        fsock = self.open('model_functions.inc', format='fortran')
 
3918
        fsock.writelines("""double complex cond
 
3919
          double complex reglog""")
 
3920
        if self.opt['mp']:
 
3921
            fsock.writelines("""%(complex_mp_format)s mp_cond
 
3922
          %(complex_mp_format)s mp_reglog"""\
 
3923
          %{'complex_mp_format':self.mp_complex_format})
 
3924
 
 
3925
    def create_model_functions_def(self):
 
3926
        """ Create model_functions.f which contains the various definitions
 
3927
        of auxiliary functions which might be used in the couplings expressions
 
3928
        Add the functions.f functions for formfactors support
 
3929
        """
 
3930
 
 
3931
        fsock = self.open('model_functions.f', format='fortran')
 
3932
        fsock.writelines("""double complex function cond(condition,truecase,falsecase)
 
3933
          implicit none
 
3934
          double complex condition,truecase,falsecase
 
3935
          if(condition.eq.(0.0d0,0.0d0)) then
 
3936
             cond=truecase
 
3937
          else
 
3938
             cond=falsecase
 
3939
          endif
 
3940
          end
 
3941
          
 
3942
          double complex function reglog(arg)
 
3943
          implicit none
 
3944
          double complex arg
 
3945
          if(arg.eq.(0.0d0,0.0d0)) then
 
3946
             reglog=(0.0d0,0.0d0)
 
3947
          else
 
3948
             reglog=log(arg)
 
3949
          endif
 
3950
          end""")
 
3951
        if self.opt['mp']:
 
3952
            fsock.writelines("""
 
3953
              
 
3954
              %(complex_mp_format)s function mp_cond(condition,truecase,falsecase)
 
3955
              implicit none
 
3956
              %(complex_mp_format)s condition,truecase,falsecase
 
3957
              if(condition.eq.(0.0e0_16,0.0e0_16)) then
 
3958
                 mp_cond=truecase
 
3959
              else
 
3960
                 mp_cond=falsecase
 
3961
              endif
 
3962
              end
 
3963
              
 
3964
              %(complex_mp_format)s function mp_reglog(arg)
 
3965
              implicit none
 
3966
              %(complex_mp_format)s arg
 
3967
              if(arg.eq.(0.0e0_16,0.0e0_16)) then
 
3968
                 mp_reglog=(0.0e0_16,0.0e0_16)
 
3969
              else
 
3970
                 mp_reglog=log(arg)
 
3971
              endif
 
3972
              end"""%{'complex_mp_format':self.mp_complex_format})            
 
3973
 
 
3974
        #check for the file functions.f
 
3975
        model_path = self.model.get('modelpath')
 
3976
        if os.path.exists(pjoin(model_path,'Fortran','functions.f')):
 
3977
            fsock.write_comment_line(' USER DEFINE FUNCTIONS ')
 
3978
            input = pjoin(self.model_path,'Fortran','functions.f')
 
3979
            fsock.writelines(open(input).read())
 
3980
            fsock.write_comment_line(' END USER DEFINE FUNCTIONS ')
3309
3981
 
3310
3982
    def create_makeinc(self):
3311
3983
        """create makeinc.inc containing the file to compile """
3312
3984
        
3313
3985
        fsock = self.open('makeinc.inc', comment='#')
3314
 
        text = 'MODEL = couplings.o lha_read.o printout.o rw_para.o '
 
3986
        text = 'MODEL = couplings.o lha_read.o printout.o rw_para.o'
 
3987
        text += ' model_functions.o '
3315
3988
        
3316
3989
        nb_coup_indep = 1 + len(self.coups_dep) // 25 
3317
3990
        nb_coup_dep = 1 + len(self.coups_indep) // 25
3318
 
        text += ' '.join(['couplings%s.o' % (i+1) \
3319
 
                                  for i in range(nb_coup_dep + nb_coup_indep) ])
 
3991
        couplings_files=['couplings%s.o' % (i+1) \
 
3992
                                for i in range(nb_coup_dep + nb_coup_indep) ]
 
3993
        if self.opt['mp']:
 
3994
            couplings_files+=['mp_couplings%s.o' % (i+1) for i in \
 
3995
                               range(nb_coup_dep,nb_coup_dep + nb_coup_indep) ]
 
3996
        text += ' '.join(couplings_files)
3320
3997
        fsock.writelines(text)
3321
3998
        
3322
3999
    def create_param_write(self):
3366
4043
        external_param = [format(param) for param in self.params_ext]
3367
4044
        fsock.writelines('\n'.join(external_param))
3368
4045
 
3369
 
        
 
4046
    def create_actualize_mp_ext_param_inc(self):
 
4047
        """ create the actualize_mp_ext_params.inc code """
 
4048
        
 
4049
        # In principle one should actualize all external, but for now, it is
 
4050
        # hardcoded that only AS and MU_R can by dynamically changed by the user
 
4051
        # so that we only update those ones.
 
4052
        # Of course, to be on the safe side, one could decide to update all
 
4053
        # external parameters.
 
4054
        update_params_list=[p for p in self.params_ext if p.name in 
 
4055
                                                          self.PS_dependent_key]
 
4056
        
 
4057
        res_strings = ["%(mp_prefix)s%(name)s=%(name)s"\
 
4058
                        %{'mp_prefix':self.mp_prefix,'name':param.name}\
 
4059
                                                for param in update_params_list]
 
4060
        # When read_lha is false, it is G which is taken in input and not AS, so
 
4061
        # this is what should be reset here too.
 
4062
        if 'aS' in [param.name for param in update_params_list]:
 
4063
            res_strings.append("%(mp_prefix)sG=G"%{'mp_prefix':self.mp_prefix})
 
4064
            
 
4065
        fsock = self.open('actualize_mp_ext_params.inc', format='fortran')
 
4066
        fsock.writelines('\n'.join(res_strings))
 
4067
 
3370
4068
    def create_param_read(self):    
3371
4069
        """create param_read"""
3372
4070
        
3373
 
        if self.opt['export_format'] == 'madevent':
 
4071
        if self.opt['export_format'] in ['madevent', 'FKS5_default', 'FKS5_optimized']:
3374
4072
            fsock = self.open('param_read.inc', format='fortran')
3375
4073
            fsock.writelines(' include \'../param_card.inc\'')
3376
4074
            return
3382
4080
            """ call LHA_get_real(npara,param,value,'%(name)s',%(name)s,%(value)s)""" \
3383
4081
                % {'name': parameter.name,
3384
4082
                   'value': self.p_to_f.parse(str(parameter.value.real))}
3385
 
        
 
4083
            if self.opt['mp']:
 
4084
                template = template+ \
 
4085
                ("\n call MP_LHA_get_real(npara,param,value,'%(name)s',"+
 
4086
                 "%(mp_prefix)s%(name)s,%(value)s)") \
 
4087
                % {'name': parameter.name,'mp_prefix': self.mp_prefix,
 
4088
                   'value': self.mp_p_to_f.parse(str(parameter.value.real))}    
3386
4089
            return template        
3387
4090
    
3388
4091
        fsock = self.open('param_read.inc', format='fortran')
3397
4100
                
3398
4101
                res_strings.append('%(width)s = sign(%(width)s,%(mass)s)' % \
3399
4102
                 {'width': particle.get('width'), 'mass': particle.get('mass')})
3400
 
                
3401
 
        
 
4103
                if self.opt['mp']:
 
4104
                    res_strings.append(\
 
4105
                      ('%(mp_pref)s%(width)s = sign(%(mp_pref)s%(width)s,'+\
 
4106
                       '%(mp_pref)s%(mass)s)')%{'width': particle.get('width'),\
 
4107
                       'mass': particle.get('mass'),'mp_pref':self.mp_prefix})
 
4108
 
3402
4109
        fsock.writelines('\n'.join(res_strings))
3403
4110
        
3404
4111
    def create_param_card(self):
3420
4127
                translator.make_valid_param_card(out_path, out_path2)
3421
4128
            translator.convert_to_slha1(out_path)
3422
4129
        
3423
 
def ExportV4Factory(cmd, noclean):
 
4130
def ExportV4Factory(cmd, noclean, output_type='default'):
3424
4131
    """ Determine which Export_v4 class is required. cmd is the command 
3425
 
        interface containing all potential usefull information."""
 
4132
        interface containing all potential usefull information.
 
4133
        The output_type argument specifies from which context the output
 
4134
        is called. It is 'madloop' for MadLoop5, 'amcatnlo' for FKS5 output
 
4135
        and 'default' for tree-level outputs."""
3426
4136
 
3427
4137
    group_subprocesses = cmd.options['group_subprocesses']
3428
 
    #check if we need to group processes
3429
 
    if cmd.options['group_subprocesses'] == 'Auto':
3430
 
        if cmd._curr_amps[0].get_ninitial()  == 2:
3431
 
            group_subprocesses = True
3432
 
        else:
3433
 
            group_subprocesses = False
3434
 
 
3435
 
    assert group_subprocesses in [True, False]
3436
 
    
3437
 
    
3438
 
    opt = {'clean': not noclean, 
3439
 
           'complex_mass': cmd.options['complex_mass_scheme'],
3440
 
           'export_format':cmd._export_format,
3441
 
           'model': cmd._curr_model.get('name')}
3442
 
    
3443
 
    if cmd._export_format in ['standalone', 'matrix']:
3444
 
        return ProcessExporterFortranSA(cmd._mgme_dir, cmd._export_dir, opt)
3445
 
    elif cmd._export_format in ['madevent'] and group_subprocesses:
3446
 
        return  ProcessExporterFortranMEGroup(cmd._mgme_dir, cmd._export_dir,
3447
 
                                                                        opt)
3448
 
    elif cmd._export_format in ['madevent']:
3449
 
        return ProcessExporterFortranME(cmd._mgme_dir, cmd._export_dir,opt)
3450
 
    
 
4138
    
 
4139
    # First treat the MadLoop5 standalone case
 
4140
    if output_type=='madloop':
 
4141
        import madgraph.loop.loop_exporters as loop_exporters
 
4142
        if os.path.isdir(os.path.join(cmd._mgme_dir, 'Template/loop_material')):
 
4143
            ExporterClass=None
 
4144
            options = {'clean': not noclean, 
 
4145
              'complex_mass':cmd.options['complex_mass_scheme'],
 
4146
              'export_format':'madloop', 
 
4147
              'mp':True,
 
4148
              'loop_dir': os.path.join(cmd._mgme_dir, 'Template/loop_material'),
 
4149
              'cuttools_dir': cmd._cuttools_dir,
 
4150
              'fortran_compiler':cmd.options['fortran_compiler']}
 
4151
 
 
4152
            if not cmd.options['loop_optimized_output']:
 
4153
                ExporterClass=loop_exporters.LoopProcessExporterFortranSA
 
4154
            else:
 
4155
                if all([amp['process']['has_born'] for amp in cmd._curr_amps]):
 
4156
                    ExporterClass=loop_exporters.LoopProcessOptimizedExporterFortranSA
 
4157
                    options['export_format'] = 'madloop_optimized'
 
4158
                else:
 
4159
                    logger.warning('ML5 can only exploit the optimized output for '+\
 
4160
                                   ' processes with born diagrams. The optimization '+\
 
4161
                                   ' option is therefore turned off for this process.')
 
4162
                    ExporterClass=loop_exporters.LoopProcessExporterFortranSA
 
4163
                    options['export_format'] = 'madloop_default'
 
4164
            return ExporterClass(cmd._mgme_dir, cmd._export_dir, options)
 
4165
        else:
 
4166
            raise MadGraph5Error('MG5_aMC cannot find the \'loop_material\' directory'+\
 
4167
                                 ' in %s'%str(cmd._mgme_dir))
 
4168
 
 
4169
    # Then treat the aMC@NLO output     
 
4170
    elif output_type=='amcatnlo':
 
4171
        import madgraph.iolibs.export_fks as export_fks
 
4172
        ExporterClass=None
 
4173
        options = {'clean': not noclean, 
 
4174
              'complex_mass':cmd.options['complex_mass_scheme'],
 
4175
              'export_format':'madloop', 
 
4176
              #use MP for HELAS only if there are virtual amps 
 
4177
              'mp':len(cmd._fks_multi_proc.get_virt_amplitudes()) > 0,
 
4178
              'loop_dir': os.path.join(cmd._mgme_dir,'Template','loop_material'),
 
4179
              'cuttools_dir': cmd._cuttools_dir,
 
4180
              'fortran_compiler':cmd.options['fortran_compiler']}
 
4181
        if not cmd.options['loop_optimized_output']:
 
4182
            logger.info("Writing out the aMC@NLO code")
 
4183
            ExporterClass = export_fks.ProcessExporterFortranFKS
 
4184
            options['export_format']='FKS5_default'
 
4185
        else:
 
4186
            logger.info("Writing out the aMC@NLO code, using optimized Loops")
 
4187
            ExporterClass = export_fks.ProcessOptimizedExporterFortranFKS
 
4188
            options['export_format']='FKS5_optimized'
 
4189
        return ExporterClass(cmd._mgme_dir, cmd._export_dir, options)
 
4190
 
 
4191
    # Then the default tree-level output
 
4192
    elif output_type=='default':
 
4193
        #check if we need to group processes
 
4194
        if cmd.options['group_subprocesses'] == 'Auto':
 
4195
            if cmd._curr_amps[0].get_ninitial()  == 2:
 
4196
                group_subprocesses = True
 
4197
            else:
 
4198
                group_subprocesses = False
 
4199
 
 
4200
        assert group_subprocesses in [True, False]
 
4201
        
 
4202
        opt = {'clean': not noclean,
 
4203
               'complex_mass': cmd.options['complex_mass_scheme'],
 
4204
               'export_format':cmd._export_format,
 
4205
               'mp': False,  
 
4206
               'sa_symmetry':False, 
 
4207
               'model': cmd._curr_model.get('name') }
 
4208
 
 
4209
        format = cmd._export_format #shortcut
 
4210
 
 
4211
        if format in ['standalone_msP', 'standalone_msF', 'standalone_rw']:
 
4212
            opt['sa_symmetry'] = True        
 
4213
    
 
4214
        if format == 'matrix' or format.startswith('standalone'):
 
4215
            return ProcessExporterFortranSA(cmd._mgme_dir, cmd._export_dir, opt,
 
4216
                                            format=format)
 
4217
        
 
4218
        elif format in ['madevent'] and group_subprocesses:
 
4219
            return  ProcessExporterFortranMEGroup(cmd._mgme_dir, cmd._export_dir,
 
4220
                                                                            opt)
 
4221
        elif format in ['madevent']:
 
4222
            return ProcessExporterFortranME(cmd._mgme_dir, cmd._export_dir,opt)
 
4223
        
 
4224
        else:
 
4225
            raise Exception, 'Wrong export_v4 format'
3451
4226
    else:
3452
 
        raise Exception, 'Wrong export_v4 format'
3453
 
        
3454
 
    
3455
 
    
 
4227
        raise MadGraph5Error, 'Output type %s not reckognized in ExportV4Factory.'
3456
4228
    
3457
4229
            
 
4230