Computing the interference of loop-induced diagrams with a tree-level background with MadEvent in MG5aMC

This page describes the technical modifications necessary to compute the interference of loop-induced diagrams with tree-level diagrams using MadEvent and the Loop-Induced module of MG5aMC.

Notice that there is as of now

/!\ no official support for such computation in MG5aMC /!\

, so that the procedure below must be performed with care and should be cross-checks in some limits since

/!\ we cannot guarantee the same level of reliability as with other MG5aMC functionalities /!\

It was however already successfully used and tested for 'g g > h > t t~ g' X 'g g > t t~ g'.

You can try to adapt this recipe to the latest version of MG5aMC or use directly v2.3.3 for which you should have a one-to-one correspondance with the intstructions.

---

a) First download the latest version of MG5_aMC, either using the tarball link:

​https://launchpad.net/mg5amcnlo/2.0/2.3.x/+download/MG5_aMC_v2.3.3.tar.gz

b) Untar the 'special_loop_sm' UFO loop model file attached to this page and place it in the 'model' directory of your MG distribution. Notice that this model is a simple copy of the default 'loop_sm' model but contains modifications in the file 'CT_vertices' which now contains special 'UVtree' counterterms, which are precisely those which will emulate the QCD background you are interfering against. Have a look at them; they are placed under the comment:

#########################################################################
# Fake UVCT vertices for loop-induced and tree interference computation #
#########################################################################

Depending on the tree-level diagrams you will want to interfere against you might need to place additional copies of the regular tree-level vertices here, with the type UVTree and involving couplings defined in 'CT_couplings.py' with the coupling_order type 'BKGQCD'. You would for example need to add copies of the electroweak vertices for example if they played a role in your tree-level background. In your case however, the tree-level diags are g g > t t~, so this model should be fine.

c) Bring the following changes to the MG5_aMC python source code:

c.1) Anywhere in the class 'LoopUVCTDiagram' of the file '<MG_root_path>/madgraph/loop/loop_base_objects.py', add the following two functions:

    def get_contracted_loop_diagram(self, model, struct_rep=None):
        return copy.copy(self)    
    def get_contracted_loop_diagram_without_tag(self, struct_rep=None):
        return copy.copy(self)

c.2) In the file '<MG_root_path>/madgraph/loop/loop_diagram_generation.py' change the folllowing:

c.2.i) Around line 740, change:

        if self['process']['has_born']:
            self.set_Born_CT()

into

        if True:
            self.set_Born_CT()

c.2.ii) Around line 1190, change

        UVCTsuccessful, UVCTdiagrams = \
          super(LoopAmplitude, self).generate_diagrams(True)

into

        # Back up of the original process characteristics
        bu = copy.copy(self['process']['required_s_channels'])
        bu_orders = copy.copy(self['process']['orders'])
        self['process']['required_s_channels'] = []
        for order in self['process']['model'].get_coupling_orders():
            self['process']['orders'][order] = 0
        self['process']['orders']['UVCT_SPECIAL'] = 999
        # Specify here the coupling orders desired for the background. By default, all.
        self['process']['orders']['BKGQCD'] = 999
        UVCTsuccessful, UVCTdiagrams = \
          super(LoopAmplitude, self).generate_diagrams(True)
        # Restore backed up properties
        self['process']['required_s_channels'] = copy.copy(bu)        
        for order, value in bu_orders.items():
            self['process']['orders'][order] = value

(Note that you can add above further modifications of self['process']['orders'] and self['process']['required_s_channels'] to have a finer selection of your tree-level interfering diagrams, but for your purposes that shouldn't be needed. Notice that thanks to the line self['process']['required_s_channels'] = [], the '> h >' constraint on the original process definition will not apply to the interfering tree-level diagrams.)

c.2.iii) Around line 1212, change

            if UVCTdiag.get_order('UVCT_SPECIAL')==1:

into

            if True:

c.3) Around line 1900 of the file '<MG_root_path>/madgraph/loop/loop_helas_objects.py' change:

            ref_orders = [lao[0] for lao in loop_orders+ct_amp_orders]

into

            ref_orders = [lao[0] for lao in loop_orders+ct_amp_orders+uvct_amp_orders]

c.4) On MG5aMC versions below or equal to 2.5.0: around line 480 of the file '<MG_root_path>/madgraph/madevent/gen_ximprove.py' change:

            ratio = th_maxwgt[-1][0]/th_maxwgt[-2][0]

to

            ratio = th_maxwgt[-1][0]/max(th_maxwgt[-2][0],1.0e-99)

d) Now that all the necessary modifications have been performed, you should be able to run the code for simulation of the loop-induced vs tree interference as follows:

d.1) First start the interface by running the script './bin/mg5' from within the directory <MG_root_path>. Make sure to always generate a new process from a freshly started interface, since the above changes introduce border effects affecting subsequent process generations from the same MG5aMC interactive session.

d.2) Whenever you plan on running event generation (i.e. not necessary when doing MadLoop standalone), make sure that the grouping of subprocess is disabled. Note that disabling the subprocess grouping optimization is no longer a requirement when starting this procedure from MG5aMC v2.5.0 and beyond.

    MG5_aMC>set group_subprocesses False

d.3) then load the model with

    MG5_aMC>import model special_loop_sm

This loads the special model for this sort of computation (with the fake 'background' UVCT vertices). (Here I show an example for the process g g > h > t t~ g, so that it is clear that it is more general that what you wanted to do in the first place, namely g g > h > t t~)

d.4) Generate the process in one of the two following ways (feel free to change the process definition at will):

If you want to run MadLoop in standalone to evaluate this interference for local phase-space points for cross-checks, run:

    MG5_aMC>generate g g > h > t t~ g [virt=QCD BKGQCD]

(You will automatically be able to see the individual contribution of the signal squared, interference term and background squared) Instead, for the generation of the process which will be output with MadEvent to generate events, use:

    MG5_aMC>generate g g > h > t t~ g [noborn=QCD BKGQCD] BKGQCD^2==3

for the generation of the process which will be output with MadEvent to generate events for this interference only.

Notice that in the second syntax, you can substitute the initial state gluons with 'p' and the final-state ones with 'j' to include the quark channels as well. This works.

d.5) You can now output the process to a directory of a chosen name 'MyOutputDirName'

MG5_aMC>output MyOutputDirName

d.6) And finally launch it with

MG5_aMC>launch

or alternatively, you can specify the directory name if you are returning to the interface and want to launch an existing output:

MG5_aMC>launch MyOutputDirName

d.7) MadGraph5_aMC@NLO will then ask you questions related to the details of either the MadLoop standalone run you are doing or the event generation run.

---

Optional

---

You can decide to use a normal model instead of the special one suggested here. Then, one can use an ad-hoc new command in MG5aMC to create a copy of all base interaction and orders so as to generate their 'BKG_' equivalent. This can be done by adding the following in the class 'MadGraphCmd' in '<MG_root_path>/madgraph/interface/madgraph_interface.py':

    def do_prepare_model_for_loopInducedXTrees(self, line):
        """Commands for adding the necessary TREE UV interactions to the current model."""
        
        # Short-hand to the current model that will be modified
        model = self._curr_model

        # First add the necessary BKG coupling_orders to the available coupling_orders 
        # in the model
        if model['order_hierarchy']:
            for coupling_order in model['order_hierarchy'].keys():
                model['order_hierarchy']['BKG_%s'%coupling_order]=model['order_hierarchy'][coupling_order]
                model['perturbation_couplings'].append('BKG_%s'%coupling_order)

        new_order_hierarchy = dict(model['order_hierarchy'])
        new_perturbation_couplings = list(model['perturbation_couplings'])

        # Dubplicate base interactions to make UV ones
        new_interactions=[]
        id_offset = 1000000
        for inter in model['interactions'].get_type('base'):
            new_interactions.append(copy.deepcopy(inter))
            # Deepcopy intorrectly treats colorstring object. Need to do it by hand.
            new_interactions[-1].set('color',[cs.create_copy() for cs in inter.get('color')])
            id_offset += 1
            new_interactions[-1].set('id',id_offset)            
            new_interactions[-1].set('type','UVtree')
            new_interactions[-1].set('loop_particles',[[]])
            new_orders = {}
            for order in new_interactions[-1]['orders']:
                new_orders['BKG_%s'%order]=new_interactions[-1]['orders'][order]
            new_interactions[-1].set('orders',new_orders)


        model.set('interactions',base_objects.InteractionList(model.get('interactions')+new_interactions))
        # Reset the model dictionaries to sync changes
        model.reset_dictionaries()
        # Refresh dictionaries manually (would be done automatically anyway)
        model.actualize_dictionaries()
       
        model.set('order_hierarchy', new_order_hierarchy)
        model.set('perturbation_couplings', new_perturbation_couplings)

        logger.info('Model %s successfully patched for (finite) Loops x Trees interference computation.'%model.get('name'))
        return

And then add the following in the class 'Switcher' of '<MG_root_path>/madgraph/interface/master_interface.py':

    def do_prepare_model_for_loopInducedXTrees(self, *args, **opts):
        return self.cmd.do_prepare_model_for_loopInducedXTrees(self, *args, **opts)

You can then simply run the command 'prepare_model_for_loopInducedXTrees' to preprocess your model before running loop-induced x tree computation as instructed above.

Notice that if the "background" tree-level structure is very different than the loop-induced ones you are interfering against, then it might be a good idea to add the corresponding tree-level topologies to the multi-channel. This can be achieved with the following small modifications. Around line 104 of the file '<MG_root_path>/madgraph/various/diagram_symmetry.py' , to the lines

        base_diagrams = base_objects.DiagramList(
                   [(d.get_contracted_loop_diagram(base_model,FDStructRepo) if
                   isinstance(d,loop_base_objects.LoopDiagram) else d) for d in
               matrix_element.get('base_amplitude').get('loop_diagrams') \
                                                            if d.get('type')>0])
        diagrams = matrix_element.get_loop_diagrams()

Add the following ones right below:

        base_diagrams.extend([(d.get_contracted_loop_diagram(base_model,FDStructRepo) if
                   isinstance(d,loop_base_objects.LoopUVCTDiagram) else d) for d in
               matrix_element.get('base_amplitude').get('loop_UVCT_diagrams')])
        diagrams.extend(matrix_element.get_loop_UVCT_diagrams())

And then around the line 317 of '<MG_root_path>/madgraph/iolibs/group_subprocs.py' , to the line

       diagrams = [(d.get_contracted_loop_diagram(model,FDStructRepo) if
                   isinstance(d,loop_base_objects.LoopDiagram) else d) for d in
                 me.get('base_amplitude').get('loop_diagrams') if d.get('type')>0]

Add the following one, right below:

      diagrams.extend([(d.get_contracted_loop_diagram(model,FDStructRepo) if
                   isinstance(d,loop_base_objects.LoopUVCTDiagram) else d) for d in
                  me.get('base_amplitude').get('loop_UVCT_diagrams')])

Finally, when using these modifications and in order to activate all these extra integration channels, you must set

    MG5_aMC>set group_subprocesses False

before the generation of your process (irrespectively of your MG5aMC distribution this time).