steering file for the Zp to invisible

steering_files/Reconstruct_ZpToInvisible_Phase3.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#######################################################
+#
+# Reconstruction of events:
+# e+ e- -> mu+ mu- Zp
+#                  Zp -> invisible
+#
+# Giacomo De Pietro (2018)
+# Laura Zani (2019)
+######################################################
+
+import random
+import string
+import variables.utils as vu
+import variables.collections as vc
+from variables import variables as var
+import modularAnalysis as ma
+import basf2 as b2
+import sys
+
+data = True
+
+# Uncomment the following line when local MC data are used
+inputFile = sys.argv[1]
+outputFile = sys.argv[2]
+
+# Uncomment the following lines when grid is used
+# inputFile = 'input.root'
+# outputFile = 'output.root'
+
+
+b2.set_log_level(b2.LogLevel.ERROR)
+
+
+def getRandomId(size=6, chars=string.ascii_uppercase + string.digits):
+    return ''.join(random.choice(chars) for _ in range(size))
+
+
+# Create the main path
+main_path = b2.create_path()
+
+# Define several aliases
+var.addAlias('p_CMS', 'useCMSFrame(p)')
+var.addAlias('E_CMS', 'useCMSFrame(E)')
+var.addAlias('pt_CMS', 'useCMSFrame(pt)')
+var.addAlias('px_CMS', 'useCMSFrame(px)')
+var.addAlias('py_CMS', 'useCMSFrame(py)')
+var.addAlias('pz_CMS', 'useCMSFrame(pz)')
+var.addAlias('theta_CMS', 'useCMSFrame(theta)')
+var.addAlias('phi_CMS', 'useCMSFrame(phi)')
+var.addAlias('clusterE_NaN', 'ifNANgiveX(clusterE, -1)')
+var.addAlias('EoverP', 'formula(ifNANgiveX(clusterE, -1) / p)')
+var.addAlias('ROE_extraE', 'roeEextra(goodGamma)')
+var.addAlias('ROE_nTracks', 'nROE_Charged(goodGamma)')
+var.addAlias('nCleanedTracks', 'nCleanedTracks([abs(dr) < 2] and [abs(dz) < 5])')
+var.addAlias('daughterAngle_CMS', 'useCMSFrame(daughterAngle(0,1))')
+var.addAlias('daughterAngle_LAB', 'daughterAngle(0,1)')
+
+trigger_bits = [
+    'a',
+    'aa',
+    'aaa',
+    'aaao',
+    'aab',
+    'aao',
+    'aay',
+    'beklm',
+    'bf',
+    'bffo',
+    'bfyo',
+    'bg',
+    'bha3d',
+    'bha3dv',
+    'bha_id0',
+    'bha_id1',
+    'bha_id2',
+    'bha_id3',
+    'bhabha',
+    'bhapur',
+    'bhie',
+    'bs',
+    'bwd_seklm',
+    'by',
+    'bz',
+    'c1hie',
+    'c2',
+    'c2gev1',
+    'c2gev2',
+    'c2gev3',
+    'c2gev4',
+    'c2hie',
+    'c3',
+    'c4',
+    'c4v',
+    'c5',
+    'cdcbb',
+    'cdcecl1',
+    'cdcecl2',
+    'cdcecl3',
+    'cdcecl4',
+    'cdcklm1',
+    'cdcklm2',
+    'cdctiming',
+    'cdctop1',
+    'cdctop2',
+    'cdctop3',
+    'cdctop4',
+    'dpee',
+    'eclbst',
+    'ecleklm1',
+    'eclmumu',
+    'ecloflo',
+    'ecltaub2b',
+    'ecltaub2b2',
+    'ecltiming',
+    'eed',
+    'eedy',
+    'eklm2',
+    'eklmhit',
+    'f',
+    'fed',
+    'ff',
+    'ff30',
+    'ffb',
+    'fff',
+    'fffo',
+    'fffov',
+    'fffv',
+    'ffo',
+    'ffoc',
+    'ffs',
+    'ffy',
+    'ffyo',
+    'ffz',
+    'fioiecl1',
+    'fp',
+    'fs',
+    'fs30',
+    'fsb',
+    'fso',
+    'fwd_seklm',
+    'fy',
+    'fy30',
+    'fyb',
+    'fyo',
+    'fyy',
+    'fz',
+    'fzb',
+    'fzo',
+    'fzz',
+    'g_c1',
+    'g_high',
+    'gg',
+    'ggbrl',
+    'ggsel',
+    'hade',
+    'hie',
+    'hie1',
+    'hie2',
+    'hie3',
+    'hiev',
+    'iecl1',
+    'ieklm1',
+    'ioiecl1',
+    'ioiecl2',
+    'itsf_b2b',
+    'klm2',
+    'klmhit',
+    'lml0',
+    'lml1',
+    'lml10',
+    'lml12',
+    'lml13',
+    'lml14',
+    'lml15',
+    'lml16',
+    'lml2',
+    'lml3',
+    'lml4',
+    'lml5',
+    'lml6',
+    'lml7',
+    'lml8',
+    'lml9',
+    'lowe',
+    'lume',
+    'lumev',
+    'mu_b2b',
+    'mu_eb2b',
+    'mu_epair',
+    'mu_pair',
+    'mu_pairv',
+    'nima01',
+    'ohem',
+    'pls',
+    'poisson',
+    'poissonv',
+    'random',
+    'revolution',
+    's',
+    'secl1',
+    'seklm1',
+    'seklm2',
+    'shem',
+    'ss',
+    'ss30',
+    'ssb',
+    'sso',
+    'sss',
+    'stt',
+    'stt4',
+    'stt5',
+    'stt6',
+    'sttecl',
+    'sttv',
+    'syb',
+    'sybecl',
+    'syo',
+    'syoecl',
+    'toptiming',
+    'vetout',
+    'y',
+    'ycdcklm1',
+    'ycdcklm2',
+    'yecl',
+    'yed',
+    'yioiecl1',
+    'yp',
+    'yy',
+    'yyv',
+    'yyy',
+    'z',
+    'zz',
+    'zzv',
+    'zzz']
+# We have some trigger bits with same names of existing variables:
+# let's treat these cases accordingly
+for bit in trigger_bits:
+    if bit not in ['random', 'y', 'z']:
+        var.addAlias(bit, f'L1PSNM({bit})')
+    else:
+        var.addAlias(f'{bit}_trigger', f'L1PSNM({bit})')
+
+# Now you can use the user-defined variables
+main_path.add_module("CustomVariables")
+var.addAlias('ptMax_CMS', 'useCMSFrame(ptmax)')
+var.addAlias('ptMin_CMS', 'useCMSFrame(ptmin)')
+var.addAlias('plongMax_CMS', 'useCMSFrame(plong_max)')
+var.addAlias('plongMin_CMS', 'useCMSFrame(plong_min)')
+var.addAlias('deltaPhiL1_LAB', 'useLabFrame(deltaPhiL1)')
+
+
+eventExtraInfoVariables = []
+
+# Function to store the KLM clusters information in the ntuples
+# Here we set how many KLM clusters we want to store
+numberOfKLMs = 10
+# Here we set the variables we want to store (only in the LAB frame)
+KLMVariables = [
+    'klmClusterInnermostLayer',
+    'klmClusterLayers',
+    'klmClusterTheta',
+    'klmClusterPhi',
+    'klmClusterTiming',
+    'klmClusterTrackDistance',
+    'nKLMClusterTrackMatches']
+
+
+def writeKLMs(
+    numberOfKLMs,
+    KLMVariables,
+    path
+):
+    """
+    Add various variables of KLM clusters to the eventExtraInfo field.
+
+    @param numberOfKLMs  Number of KLM clusters to store
+    @param KLMVariables  List of KLM clusters variables to store
+    @param path          Modules are added to this path
+    """
+
+    if (numberOfKLMs < 0):
+        b2.B2FATAL('Number of KLM clusters must be >= 0')
+    if (numberOfKLMs > 100):
+        b2.B2FATAL('Number of KLM clusters must be < 100')
+
+    KLMList = 'K_L0:klm_extraInfo'
+    ma.fillParticleList(decayString=KLMList,
+                        cut='[isFromKLM > 0]',
+                        path=path)
+    ma.rankByLowest(particleList=KLMList,
+                    variable='klmClusterLayers',
+                    numBest=numberOfKLMs,
+                    path=path)
+    for rank in range(1, numberOfKLMs + 1):
+        ma.cutAndCopyList(outputListName=f'{KLMList}_{rank}',
+                          inputListName=f'{KLMList}',
+                          cut=f'[extraInfo(klmClusterLayers_rank) == {rank}]',
+                          path=path)
+        for variable in KLMVariables:
+            ma.variablesToEventExtraInfo(particleList=f'{KLMList}_{rank}',
+                                         variables={f'{variable}': f'KLM_{rank}_{variable}'},
+                                         path=path)
+
+
+# Aliases
+for rank in range(1, numberOfKLMs + 1):
+    for variable in KLMVariables:
+        var.addAlias(f'KLM_{rank}_{variable}', f'eventExtraInfo(KLM_{rank}_{variable})')
+        eventExtraInfoVariables += [f'KLM_{rank}_{variable}', f'KLM_{rank}_{variable}']
+
+
+# Function to store the photon information in the ntuples
+# Here we set how many photons we want to store
+numberOfPhotons = 5
+# Here we set the variables we want to store
+photonVariables = ['E', 'px', 'py', 'pz']
+# And here we set other variables we want to store only in LAB frame (e.g. mdst index, etc.)
+photonOtherVariables = ['genParticleID', 'isSignal', 'mcPDG']
+
+
+def writePhotons(
+    photonSelection,
+    numberOfPhotons,
+    photonVariables,
+    photonOtherVariables,
+    path
+):
+    """
+    Add various variables of photons to the eventExtraInfo field.
+
+    @param photonSelection      Selection for the photons
+    @param numberOfPhotons      Number of photons to store
+    @param photonVariables      List of photon variables to store
+    @param photonOtherVariables List of other variables to store: the frame will not be set
+    @param path                 Modules are added to this path
+    """
+
+    if (numberOfPhotons < 0):
+        b2.B2FATAL('Number of photons must be >= 0')
+    if (numberOfPhotons > 100):
+        b2.B2FATAL('Number of photons must be < 100')
+
+    photonList = 'gamma:extraInfo'
+    ma.fillParticleList(decayString=photonList,
+                        cut=photonSelection,
+                        path=path)
+    ma.rankByHighest(particleList=photonList,
+                     variable='E',
+                     numBest=numberOfPhotons,
+                     path=path)
+    for rank in range(1, numberOfPhotons + 1):
+        ma.cutAndCopyList(outputListName=f'{photonList}_{rank}',
+                          inputListName=f'{photonList}',
+                          cut=f'[extraInfo(E_rank) == {rank}]',
+                          path=path)
+        for variable in photonVariables:
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'useLabFrame({variable})': f'photon_{rank}_{variable}'},
+                                         path=path)
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'useCMSFrame({variable})': f'photon_{rank}_{variable}_CMS'},
+                                         path=path)
+        for variable in photonOtherVariables:
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'{variable}': f'photon_{rank}_{variable}'},
+                                         path=path)
+
+
+# Aliases
+for rank in range(1, numberOfPhotons + 1):
+    for variable in photonVariables:
+        var.addAlias(f'photon_{rank}_{variable}', f'eventExtraInfo(photon_{rank}_{variable})')
+        var.addAlias(f'photon_{rank}_{variable}_CMS', f'eventExtraInfo(photon_{rank}_{variable}_CMS)')
+        eventExtraInfoVariables += [f'photon_{rank}_{variable}', f'photon_{rank}_{variable}_CMS']
+    for variable in photonOtherVariables:
+        var.addAlias(f'photon_{rank}_{variable}', f'eventExtraInfo(photon_{rank}_{variable})')
+        eventExtraInfoVariables += [f'photon_{rank}_{variable}']
+
+
+# Function to store the photon information in the ntuples
+# Here we set how many photons we want to store
+numberOfMCPhotons = 7
+# Here we set the variables we want to store
+photonMCVariables = ['E', 'px', 'py', 'pz']
+# And here we set other variables we want to store only in LAB frame (e.g. mdst index, etc.)
+photonOtherMCVariables = ['mdstIndex']
+
+
+def writeMCPhotons(
+    photonSelection,
+    numberOfPhotons,
+    photonVariables,
+    photonOtherVariables,
+    path
+):
+    """
+    Add various variables of MC photons to the eventExtraInfo field.
+
+    @param photonSelection      Selection for the photons
+    @param numberOfPhotons      Number of photons to store
+    @param photonVariables      List of photon variables to store
+    @param photonOtherVariables List of other variables to store: the frame will not be set
+    @param path                 Modules are added to this path
+    """
+
+    if (numberOfMCPhotons < 0):
+        b2.B2FATAL('Number of photons must be >= 0')
+    if (numberOfMCPhotons > 100):
+        b2.B2FATAL('Number of photons must be < 100')
+
+    photonList = 'gamma:extraInfoMC'
+    ma.fillParticleListFromMC(decayString=photonList,
+                              cut=photonSelection,
+                              path=path)
+    ma.rankByHighest(particleList=photonList,
+                     variable='E',
+                     numBest=numberOfMCPhotons,
+                     path=path)
+    for rank in range(1, numberOfMCPhotons + 1):
+        ma.cutAndCopyList(outputListName=f'{photonList}_{rank}',
+                          inputListName=f'{photonList}',
+                          cut=f'[extraInfo(E_rank) == {rank}]',
+                          path=path)
+        for variable in photonVariables:
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'useLabFrame({variable})': f'photonMC_{rank}_{variable}'},
+                                         path=path)
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'useCMSFrame({variable})': f'photonMC_{rank}_{variable}_CMS'},
+                                         path=path)
+        for variable in photonOtherVariables:
+            ma.variablesToEventExtraInfo(particleList=f'{photonList}_{rank}',
+                                         variables={f'{variable}': f'photonMC_{rank}_{variable}'},
+                                         path=path)
+
+
+# Aliases
+if not data:
+    for rank in range(1, numberOfMCPhotons + 1):
+        for variable in photonMCVariables:
+            var.addAlias(f'photonMC_{rank}_{variable}', f'eventExtraInfo(photonMC_{rank}_{variable})')
+            var.addAlias(f'photonMC_{rank}_{variable}_CMS', f'eventExtraInfo(photonMC_{rank}_{variable}_CMS)')
+            eventExtraInfoVariables += [f'photonMC_{rank}_{variable}', f'photonMC_{rank}_{variable}_CMS']
+        for variable in photonOtherMCVariables:
+            var.addAlias(f'photonMC_{rank}_{variable}', f'eventExtraInfo(photonMC_{rank}_{variable})')
+            eventExtraInfoVariables += [f'photonMC_{rank}_{variable}']
+
+
+# Function for the closest photon to the signal
+def writeClosestPhotonInfo(
+    particleList,
+    photonSelection,
+    path
+):
+    """
+    Add various variables of the closest photon to the particle extraInfo field.
+
+    @particleList          List of the input particles
+    @param photonSelection Selection for the ROE photon
+    @param path            Modules are added to this path
+    """
+
+    # Build rest of event
+    ma.buildRestOfEvent(target_list_name=particleList,
+                        path=path)
+    # Create new path for ROE
+    roe_path = b2.create_path()
+    # Get random list names (in case we run this function multiple times)
+    signalList = particleList + getRandomId()
+    photonList = 'gamma:writeClosestPhotonInfo' + getRandomId()
+    pairList = 'vpho:writeClosestPhotonInfo' + getRandomId()
+    # Pair the signal side particle with the photons
+    ma.fillSignalSideParticleList(outputListName=signalList,
+                                  decayString=f'^{particleList}',
+                                  path=roe_path)
+    ma.fillParticleList(decayString=photonList,
+                        cut=f'[isInRestOfEvent == 1] and {photonSelection}',
+                        path=roe_path)
+    ma.reconstructDecay(decayString=f'{pairList} -> {signalList} {photonList}',
+                        cut='',
+                        path=roe_path)
+    # Keep the closest photon
+    ma.rankByLowest(particleList=pairList,
+                    variable='daughterAngle(0, 1)',
+                    numBest=1,
+                    path=roe_path)
+    # Add new variables to the signal side particle
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'useLabFrame(daughterAngle(0, 1))': 'closestPhoton_angle'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'useCMSFrame(daughterAngle(0, 1))': 'closestPhoton_angle_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useLabFrame(theta))': 'closestPhoton_theta'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useCMSFrame(theta))': 'closestPhoton_theta_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useLabFrame(phi))': 'closestPhoton_phi'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useCMSFrame(phi))': 'closestPhoton_phi_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useLabFrame(E))': 'closestPhoton_E'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useCMSFrame(E))': 'closestPhoton_E_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useLabFrame(pt))': 'closestPhoton_pt'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=pairList,
+                                     varToExtraInfo={'daughter(1, useCMSFrame(pt))': 'closestPhoton_pt_CMS'},
+                                     path=roe_path)
+    # Execute the roe_path
+    path.for_each('RestOfEvent', 'RestOfEvents', roe_path)
+
+
+# Aliases for closest photon variables
+var.addAlias('closestPhoton_angle', 'extraInfo(closestPhoton_angle)')
+var.addAlias('closestPhoton_angle_CMS', 'extraInfo(closestPhoton_angle_CMS)')
+var.addAlias('closestPhoton_theta', 'extraInfo(closestPhoton_theta)')
+var.addAlias('closestPhoton_theta_CMS', 'extraInfo(closestPhoton_theta_CMS)')
+var.addAlias('closestPhoton_phi', 'extraInfo(closestPhoton_phi)')
+var.addAlias('closestPhoton_phi_CMS', 'extraInfo(closestPhoton_phi_CMS)')
+var.addAlias('closestPhoton_E', 'extraInfo(closestPhoton_E)')
+var.addAlias('closestPhoton_E_CMS', 'extraInfo(closestPhoton_E_CMS)')
+var.addAlias('closestPhoton_pt', 'extraInfo(closestPhoton_pt)')
+var.addAlias('closestPhoton_pt_CMS', 'extraInfo(closestPhoton_pt_CMS)')
+
+# Function for most energetic photon in the event
+
+
+def writeMostEnergeticPhotonInfo(
+    particleList,
+    photonSelection,
+    path
+):
+    """
+    Add various variables of the most energetic photon to the particle extraInfo field.
+
+    @particleList          List of the input particles
+    @param photonSelection Selection for the ROE photon
+    @param path            Modules are added to this path
+    """
+
+    # Build rest of event
+    ma.buildRestOfEvent(target_list_name=particleList,
+                        path=path)
+    # Create new path for ROE
+    roe_path = b2.create_path()
+    # Get random list names (in case we run this function multiple times)
+    signalList = particleList + getRandomId()
+    photonList = 'gamma:writeMostEnergeticPhotonInfo' + getRandomId()
+    # Pair the signal side particle with the photons
+    ma.fillSignalSideParticleList(outputListName=signalList,
+                                  decayString=f'^{particleList}',
+                                  path=roe_path)
+    ma.fillParticleList(decayString=photonList,
+                        cut=f'[isInRestOfEvent == 1] and {photonSelection}',
+                        path=roe_path)
+    # Add the ranking per highest photon energy
+    ma.rankByHighest(particleList=photonList,
+                     variable='E_CMS',
+                     numBest=1,
+                     path=roe_path)
+    # Most energetic photon: Add new variables to the signal side particle
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'theta': 'mostEnergeticPhoton_theta'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'theta_CMS': 'mostEnergeticPhoton_theta_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'phi': 'mostEnergeticPhoton_phi'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'phi_CMS': 'mostEnergeticPhoton_phi_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'E': 'mostEnergeticPhoton_E'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'E_CMS': 'mostEnergeticPhoton_E_CMS'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'pt': 'mostEnergeticPhoton_pt'},
+                                     path=roe_path)
+    ma.variableToSignalSideExtraInfo(particleList=photonList,
+                                     varToExtraInfo={'pt_CMS': 'mostEnergeticPhoton_pt_CMS'},
+                                     path=roe_path)
+    # Execute the roe_path
+    path.for_each('RestOfEvent', 'RestOfEvents', roe_path)
+
+
+# Aliases for most energetic photon variables
+var.addAlias('mostEnergeticPhoton_theta', 'extraInfo(mostEnergeticPhoton_theta)')
+var.addAlias('mostEnergeticPhoton_theta_CMS', 'extraInfo(mostEnergeticPhoton_theta_CMS)')
+var.addAlias('mostEnergeticPhoton_phi', 'extraInfo(mostEnergeticPhoton_phi)')
+var.addAlias('mostEnergeticPhoton_phi_CMS', 'extraInfo(mostEnergeticPhoton_phi_CMS)')
+var.addAlias('mostEnergeticPhoton_E', 'extraInfo(mostEnergeticPhoton_E)')
+var.addAlias('mostEnergeticPhoton_E_CMS', 'extraInfo(mostEnergeticPhoton_E_CMS)')
+var.addAlias('mostEnergeticPhoton_pt', 'extraInfo(mostEnergeticPhoton_pt)')
+var.addAlias('mostEnergeticPhoton_pt_CMS', 'extraInfo(mostEnergeticPhoton_pt_CMS)')
+
+# Function for the pi0 veto
+
+
+def writePi0Veto(
+    particleList,
+    photonSelection,
+    path
+):
+    """
+    Add various variables of the MostEnergetic photon to the particle extraInfo field.
+
+    @particleList          List of the input particles
+    @param photonSelection Selection for the ROE photo
+    @param path            Modules are added to this path
+    """
+
+    # Build rest of event
+    ma.buildRestOfEvent(target_list_name=particleList,
+                        path=path)
+    # Create new path for ROE
+    roe_path = b2.create_path()
+    # Get random list names (in case we run this function multiple times)
+    signalList = particleList + getRandomId()
+    photonList = 'gamma:writePi0Veto' + getRandomId()
+    pi0List = 'pi0:writePi0Veto' + getRandomId()
+    # Pair the signal side particle with the photons
+    ma.fillSignalSideParticleList(outputListName=signalList,
+                                  decayString=f'^{particleList}',
+                                  path=roe_path)
+    ma.fillParticleList(decayString=photonList,
+                        cut=f'[isInRestOfEvent == 1] and {photonSelection}',
+                        path=roe_path)
+    ma.reconstructDecay(decayString=f'{pi0List} -> {photonList} {photonList}',
+                        cut='0.080 < M < 0.200',
+                        path=roe_path)
+    # Add the ranking per closest mass to pi0 mass
+    ma.rankByLowest(particleList=pi0List,
+                    variable='abs(dM)',
+                    numBest=1,
+                    path=roe_path)
+    # Most energetic photon: Add new variables to the signal side particle
+    ma.variableToSignalSideExtraInfo(particleList=pi0List,
+                                     varToExtraInfo={'M': 'pi0veto_M'},
+                                     path=roe_path)
+    # Execute the roe_path
+    path.for_each('RestOfEvent', 'RestOfEvents', roe_path)
+
+
+var.addAlias('pi0veto_M', 'extraInfo(pi0veto_M)')
+
+# Import input files
+main_path.add_module('RootInput', inputFileNames=inputFile)
+
+main_path.add_module('Progress')
+
+# Create "mu+:selected" ParticleList and apply some cuts
+ma.fillParticleList(decayString='mu+:selected',
+                    cut='[abs(dr) < 2] and [abs(dz) < 5] and [muonID_noSVD > 0.1]',
+                    path=main_path)
+
+# Correct data applying the momentumm scale due to the B field
+if data:
+    ma.scaleTrackMomenta(inputListNames='mu+:selected',
+                         scale=1.00056,
+                         path=main_path)
+
+# Create muon pairs
+ma.reconstructDecay(decayString='Z0:mumu -> mu+:selected mu-:selected',
+                    cut='[nCleanedTracks < 4]',
+                    path=main_path)
+
+# Perform the MC matching
+ma.matchMCTruth(list_name='Z0:mumu',
+                path=main_path)
+
+# Define some sets of cuts
+trackSelection = ''
+gammaSelection = '[theta > 0.296706] and [theta < 2.61799] and [clusterErrorTiming < 1e6] and [clusterE1E9 > 0.4 or E > 0.075]'
+
+# Reconstruct the recoil
+ma.reconstructRecoil(decayString='Z0:recoil -> Z0:mumu',
+                     cut='[M < 12.]',
+                     path=main_path)
+
+# Store the information of the closest photon wrt the recoil
+writeClosestPhotonInfo(particleList='Z0:recoil',
+                       photonSelection=f'{gammaSelection} and [E >= 0.1]',
+                       path=main_path)
+
+# Store the information of the most energetic photon
+writeMostEnergeticPhotonInfo(particleList='Z0:recoil',
+                             photonSelection=f'{gammaSelection} and [E >= 0.1]',
+                             path=main_path)
+
+# Store the information of KLM clusters
+writeKLMs(numberOfKLMs=numberOfKLMs,
+          KLMVariables=KLMVariables,
+          path=main_path)
+
+# Store the information of the photons
+writePhotons(photonSelection=f'{gammaSelection} and [E >= 0.1]',
+             numberOfPhotons=numberOfPhotons,
+             photonVariables=photonVariables,
+             photonOtherVariables=photonOtherVariables,
+             path=main_path)
+
+# Store the information of the MC photons
+if not data:
+    writeMCPhotons(photonSelection='',
+                   numberOfPhotons=numberOfMCPhotons,
+                   photonVariables=photonMCVariables,
+                   photonOtherVariables=photonOtherMCVariables,
+                   path=main_path)
+
+# Store the information of the "best" pi0 in the ROE
+writePi0Veto(particleList='Z0:recoil',
+             photonSelection=gammaSelection,
+             path=main_path)
+
+# Build the ROE to "check" the extra energy in the ECL
+ma.buildRestOfEvent(target_list_name='Z0:recoil',
+                    path=main_path)
+
+# Create a ROE mask for extra energy calculation
+ma.appendROEMask(list_name='Z0:recoil',
+                 mask_name='goodGamma',
+                 trackSelection=trackSelection,
+                 eclClusterSelection=f'{gammaSelection} and [E >= 0.1]',
+                 path=main_path)
+
+# Create a brand new photon list for an alternative ROE and extraEnergy definition
+ma.fillParticleList('gamma:forROE',
+                    f'{gammaSelection} and [E >= 0.1]',
+                    True,
+                    path=main_path)
+
+# Define the list of variables to be stored in the output .root file
+kinematics_variables = vc.kinematics + ['theta',
+                                        'phi'] + ['E_CMS',
+                                                  'p_CMS',
+                                                  'pt_CMS',
+                                                  'px_CMS',
+                                                  'py_CMS',
+                                                  'pz_CMS',
+                                                  'theta_CMS',
+                                                  'phi_CMS']
+
+closestPhoton_variables = ['closestPhoton_angle',
+                           'closestPhoton_angle_CMS',
+                           'closestPhoton_theta',
+                           'closestPhoton_theta_CMS',
+                           'closestPhoton_phi',
+                           'closestPhoton_phi_CMS',
+                           'closestPhoton_E',
+                           'closestPhoton_E_CMS',
+                           'closestPhoton_pt',
+                           'closestPhoton_pt_CMS']
+
+mostEnergeticPhoton_variables = ['mostEnergeticPhoton_theta',
+                                 'mostEnergeticPhoton_theta_CMS',
+                                 'mostEnergeticPhoton_phi',
+                                 'mostEnergeticPhoton_phi_CMS',
+                                 'mostEnergeticPhoton_E',
+                                 'mostEnergeticPhoton_E_CMS',
+                                 'mostEnergeticPhoton_pt',
+                                 'mostEnergeticPhoton_pt_CMS']
+
+recoil_variables = vc.inv_mass + ['M2'] + kinematics_variables + \
+    closestPhoton_variables + mostEnergeticPhoton_variables + ['pi0veto_M']
+recoil_variables += ['ROE_extraE',
+                     'ROE_nTracks',
+                     'nTracks',
+                     'nCleanedTracks']
+recoil_variables += ['ptMax_CMS',
+                     'ptMin_CMS',
+                     'plongMax_CMS',
+                     'plongMin_CMS',
+                     'deltaPhiL1_LAB']
+
+dimuon_variables = vc.inv_mass
+dimuon_variables += ['M2',
+                     'mRecoil',
+                     'm2Recoil']
+dimuon_variables += kinematics_variables
+dimuon_variables += ['daughterAngle_LAB',
+                     'daughterAngle_CMS']
+recoil_variables += vu.create_aliases_for_selected(list_of_variables=dimuon_variables,
+                                                   decay_string='Z0:recoil -> ^Z0:mumu')
+
+muon_variables = kinematics_variables
+muon_variables += vc.track
+muon_variables += vc.track_hits
+muon_variables += vc.pid + ['muonID_noSVD']
+muon_variables += ['clusterE_NaN',
+                   'EoverP']
+muon_variables += vc.mc_truth
+muon_variables += ['isSignal',
+                   'charge']
+recoil_variables += vu.create_aliases_for_selected(list_of_variables=muon_variables,
+                                                   decay_string='Z0:recoil -> [Z0:mumu -> ^mu+:selected ^mu-:selected]')
+recoil_variables += trigger_bits
+recoil_variables += eventExtraInfoVariables
+
+# Rest of event and photons
+var.addAlias('nPhotons', 'countInList(gamma:forROE)')
+var.addAlias('sumPhotonsE', 'totalEnergyOfParticlesInList(gamma:forROE)')
+var.addAlias('sumPhotonsPx', 'totalPxOfParticlesInList(gamma:forROE)')
+var.addAlias('sumPhotonsPy', 'totalPyOfParticlesInList(gamma:forROE)')
+var.addAlias('sumPhotonsPz', 'totalPzOfParticlesInList(gamma:forROE)')
+recoil_variables += ['nPhotons', 'sumPhotonsE', 'sumPhotonsPx', 'sumPhotonsPy', 'sumPhotonsPz']
+
+# Print all the aliases
+var.printAliases()
+
+# Saving the variables into a .root file
+ma.variablesToNtuple(decayString='Z0:recoil',
+                     variables=recoil_variables,
+                     filename=outputFile,
+                     treename='Z0_recoil',
+                     path=main_path)
+
+# Process the events
+b2.process(main_path)
+
+# print out the summary
+print(b2.statistics)