Updating FCC-ee mH recoil mumu example to work with edm4hep v1 format

.github/workflows/test.yml

@@ -52,4 +52,4 @@ jobs:
         docker exec CI_container /bin/bash -c 'cd ./Package; \
           source ${{ matrix.STACK }}; \
           source ./setup.sh; \
-          fccanalysis run examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py --output myoutput.root --files-list root://eospublic.cern.ch//eos/experiment/fcc/ee/generation/DelphesEvents/spring2021/IDEA/p8_ee_Zbb_ecm91_EvtGen_Bc2TauNuTAUHADNU/events_131527278.root'
+          fccanalysis run examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py --output myoutput.root --test'

examples/FCCee/higgs/mH-recoil/histmaker_recoil.py

@@ -0,0 +1,185 @@
+# list of processes (mandatory)
+processList = {
+    #'p8_ee_ZZ_ecm240':{'fraction':1},
+    #'p8_ee_WW_ecm240':{'fraction':1}, 
+    #'wzp6_ee_mumuH_ecm240':{'fraction':1},
+    # 'p8_ee_WW_mumu_ecm240':    {'fraction':1, 'crossSection': 0.25792}, 
+    # 'p8_ee_ZZ_mumubb_ecm240':  {'fraction':1, 'crossSection': 2 * 1.35899 * 0.034 * 0.152},
+    # 'p8_ee_ZH_Zmumu_ecm240':   {'fraction':1, 'crossSection': 0.201868 * 0.034},
+    'p8_ee_ZZ_ecm240':    {'fraction':1, 'crossSection': 0.25792}, 
+    'p8_ee_WW_ecm240':  {'fraction':1, 'crossSection': 2 * 1.35899 * 0.034 * 0.152},
+    'p8_ee_ZH_ecm240':   {'fraction':1, 'crossSection': 0.201868 * 0.034},
+}
+
+# Production tag when running over EDM4Hep centrally produced events, this points to the yaml files for getting sample statistics (mandatory)
+#prodTag     = "FCCee/winter2023/IDEA/"
+
+# Link to the dictonary that contains all the cross section informations etc... (mandatory)
+procDict = "FCCee_procDict_winter2023_IDEA.json"
+
+# additional/custom C++ functions, defined in header files (optional)
+includePaths = ["functions.h"]
+
+# Define the input dir (optional)
+#inputDir    = "outputs/FCCee/higgs/mH-recoil/mumu/stage1"
+inputDir    = "/eos/experiment/fcc/hh/tutorials/edm4hep_tutorial_data/"
+
+#Optional: output directory, default is local running directory
+outputDir   = "./histmaker_outputs/"
+
+
+# optional: ncpus, default is 4, -1 uses all cores available
+nCPUS       = -1
+
+# scale the histograms with the cross-section and integrated luminosity
+doScale = True
+intLumi = 5000000 # 5 /ab
+
+
+# define some binning for various histograms
+bins_p_mu = (2000, 0, 200) # 100 MeV bins
+bins_m_ll = (2000, 0, 200) # 100 MeV bins
+bins_p_ll = (2000, 0, 200) # 100 MeV bins
+bins_recoil = (200000, 0, 200) # 1 MeV bins 
+bins_cosThetaMiss = (10000, 0, 1)
+
+bins_theta = (500, -5, 5)
+bins_eta = (600, -3, 3)
+bins_phi = (500, -5, 5)
+
+bins_count = (50, 0, 50)
+bins_charge = (10, -5, 5)
+bins_iso = (500, 0, 5)
+
+
+
+# build_graph function that contains the analysis logic, cuts and histograms (mandatory)
+def build_graph(df, dataset):
+
+    results = []
+    df = df.Define("weight", "1.0")
+    weightsum = df.Sum("weight")
+
+    # define some aliases to be used later on
+    df = df.Alias("Particle0", "_Particle_daughters.index")
+    df = df.Alias("Particle1", "_Particle_parents.index")
+
+    df = df.Alias("MCRecoAssociations0", "_MCRecoAssociations_from.index")
+    df = df.Alias("MCRecoAssociations1", "_MCRecoAssociations_to.index")
+
+    df = df.Alias("Muon0", "Muon_objIdx.index")
+
+
+
+    # get all the leptons from the collection
+    df = df.Define("muons_all", "FCCAnalyses::ReconstructedParticle::get(Muon0, ReconstructedParticles)")
+
+
+    # select leptons with momentum > 20 GeV
+    df = df.Define("muons", "FCCAnalyses::ReconstructedParticle::sel_p(20)(muons_all)")
+
+
+    df = df.Define("muons_p", "FCCAnalyses::ReconstructedParticle::get_p(muons)")
+    df = df.Define("muons_theta", "FCCAnalyses::ReconstructedParticle::get_theta(muons)")
+    df = df.Define("muons_phi", "FCCAnalyses::ReconstructedParticle::get_phi(muons)")
+    df = df.Define("muons_q", "FCCAnalyses::ReconstructedParticle::get_charge(muons)")
+    df = df.Define("muons_no", "FCCAnalyses::ReconstructedParticle::get_n(muons)")
+
+    # compute the muon isolation and store muons with an isolation cut of 0.25 in a separate column muons_sel_iso
+    df = df.Define("muons_iso", "FCCAnalyses::ZHfunctions::coneIsolation(0.01, 0.5)(muons, ReconstructedParticles)")
+    df = df.Define("muons_sel_iso", "FCCAnalyses::ZHfunctions::sel_iso(0.25)(muons, muons_iso)")
+
+
+    # baseline histograms, before any selection cuts (store with _cut0)
+    results.append(df.Histo1D(("muons_p_cut0", "", *bins_p_mu), "muons_p"))
+    results.append(df.Histo1D(("muons_theta_cut0", "", *bins_theta), "muons_theta"))
+    results.append(df.Histo1D(("muons_phi_cut0", "", *bins_phi), "muons_phi"))
+    results.append(df.Histo1D(("muons_q_cut0", "", *bins_charge), "muons_q"))
+    results.append(df.Histo1D(("muons_no_cut0", "", *bins_count), "muons_no"))
+    results.append(df.Histo1D(("muons_iso_cut0", "", *bins_iso), "muons_iso"))
+
+
+    #########
+    ### CUT 0: all events
+    #########
+    df = df.Define("cut0", "0")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut0"))
+
+
+    #########
+    ### CUT 1: at least 1 muon with at least one isolated one
+    #########
+    df = df.Filter("muons_no >= 1 && muons_sel_iso.size() > 0")
+    df = df.Define("cut1", "1")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut1"))
+
+
+    #########
+    ### CUT 2 :at least 2 opposite-sign (OS) leptons
+    #########
+    df = df.Filter("muons_no >= 2 && abs(Sum(muons_q)) < muons_q.size()")
+    df = df.Define("cut2", "2")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut2"))
+
+    # now we build the Z resonance based on the available leptons.
+    # the function resonanceBuilder_mass_recoil returns the best lepton pair compatible with the Z mass (91.2 GeV) and recoil at 125 GeV
+    # the argument 0.4 gives a weight to the Z mass and the recoil mass in the chi2 minimization
+    # technically, it returns a ReconstructedParticleData object with index 0 the di-lepton system, index and 2 the leptons of the pair
+    df = df.Define("zbuilder_result", "FCCAnalyses::ZHfunctions::resonanceBuilder_mass_recoil(91.2, 125, 0.4, 240, false)(muons, MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle, Particle0, Particle1)")
+    df = df.Define("zmumu", "Vec_rp{zbuilder_result[0]}") # the Z
+    df = df.Define("zmumu_muons", "Vec_rp{zbuilder_result[1],zbuilder_result[2]}") # the leptons 
+    df = df.Define("zmumu_m", "FCCAnalyses::ReconstructedParticle::get_mass(zmumu)[0]") # Z mass
+    df = df.Define("zmumu_p", "FCCAnalyses::ReconstructedParticle::get_p(zmumu)[0]") # momentum of the Z
+    df = df.Define("zmumu_recoil", "FCCAnalyses::ReconstructedParticle::recoilBuilder(240)(zmumu)") # compute the recoil based on the reconstructed Z
+    df = df.Define("zmumu_recoil_m", "FCCAnalyses::ReconstructedParticle::get_mass(zmumu_recoil)[0]") # recoil mass
+    df = df.Define("zmumu_muons_p", "FCCAnalyses::ReconstructedParticle::get_p(zmumu_muons)") # get the momentum of the 2 muons from the Z resonance
+
+
+
+    #########
+    ### CUT 3: Z mass window
+    #########  
+    df = df.Filter("zmumu_m > 86 && zmumu_m < 96")
+    df = df.Define("cut3", "3")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut3"))
+
+
+    #########
+    ### CUT 4: Z momentum
+    #########  
+    df = df.Filter("zmumu_p > 20 && zmumu_p < 70")
+    df = df.Define("cut4", "4")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut4"))
+
+
+    #########
+    ### CUT 5: cosThetaMiss
+    #########  
+    df = df.Define("missingEnergy", "FCCAnalyses::ZHfunctions::missingEnergy(240., ReconstructedParticles)")
+    #df = df.Define("cosTheta_miss", "FCCAnalyses::get_cosTheta_miss(missingEnergy)")
+    df = df.Define("cosTheta_miss", "FCCAnalyses::ZHfunctions::get_cosTheta_miss(MissingET)")
+    results.append(df.Histo1D(("cosThetaMiss_cut4", "", *bins_cosThetaMiss), "cosTheta_miss")) # plot it before the cut
+
+    df = df.Filter("cosTheta_miss < 0.98")
+    df = df.Define("cut5", "5")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut5"))
+
+
+    #########
+    ### CUT 6: recoil mass window
+    #########  
+    df = df.Filter("zmumu_recoil_m < 140 && zmumu_recoil_m > 120")
+    df = df.Define("cut6", "6")
+    results.append(df.Histo1D(("cutFlow", "", *bins_count), "cut6"))
+
+
+    ########################
+    # Final histograms
+    ########################
+    results.append(df.Histo1D(("zmumu_m", "", *bins_m_ll), "zmumu_m"))
+    results.append(df.Histo1D(("zmumu_recoil_m", "", *bins_recoil), "zmumu_recoil_m"))
+    results.append(df.Histo1D(("zmumu_p", "", *bins_p_ll), "zmumu_p"))
+    results.append(df.Histo1D(("zmumu_muons_p", "", *bins_p_mu), "zmumu_muons_p"))
+
+
+    return results, weightsum

examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py

@@ -23,22 +23,23 @@ def __init__(self, cmdline_args):
 
         # Mandatory: List of processes to run over
         self.process_list = {
-            # Run the full statistics in one output file named
-            # <outputDir>/p8_ee_ZZ_ecm240.root
-            'p8_ee_ZZ_ecm240': {'fraction': 0.005},
-            # Run 50% of the statistics with output into two files named
-            # <outputDir>/p8_ee_WW_ecm240/chunk<N>.root
-            'p8_ee_WW_ecm240': {'fraction': 0.5, 'chunks': 2},
-            # Run 20% of the statistics in one file named
-            # <outputDir>/p8_ee_ZH_ecm240_out.root (example on how to change
-            # the output name)
+            # # Run the full statistics in one output file named
+            # # <outputDir>/p8_ee_ZZ_ecm240.root
+            'p8_ee_ZZ_ecm240': {'fraction': 1.},
+            # # Run 50% of the statistics with output into two files named
+            # # <outputDir>/p8_ee_WW_ecm240/chunk<N>.root
+            'p8_ee_WW_ecm240': {'fraction': 0.5, 'chunks': 2}, #this doesn't work? 
+            # # Run 20% of the statistics in one file named
+            # # <outputDir>/p8_ee_ZH_ecm240_out_f02.root (example on how to change
+            # # the output name)
             'p8_ee_ZH_ecm240': {'fraction': 0.2,
-                                'output': 'p8_ee_ZH_ecm240_out'}
+                                'output': 'p8_ee_ZH_ecm240_out_f02'}
         }
 
         # Mandatory: Production tag when running over the centrally produced
         # samples, this points to the yaml files for getting sample statistics
-        self.prod_tag = 'FCCee/spring2021/IDEA/'
+        self.input_dir = '/eos/experiment/fcc/hh/tutorials/edm4hep_tutorial_data/'
+        # self.prod_tag = 'FCCee/spring2021/IDEA/'
 
         # Optional: output directory, default is local running directory
         self.output_dir = 'outputs/FCCee/higgs/mH-recoil/mumu/' \
@@ -54,9 +55,10 @@ def __init__(self, cmdline_args):
         # self.run_batch = False
 
         # Optional: test file
-        self.test_file = 'root://eospublic.cern.ch//eos/experiment/fcc/ee/' \
-                         'generation/DelphesEvents/spring2021/IDEA/' \
-                         'p8_ee_ZH_ecm240/events_101027117.root'
+        self.test_file = 'root://eospublic.cern.ch//eos/experiment/fcc/hh/' \
+                         'tutorials/edm4hep_tutorial_data/' \
+                         'p8_ee_ZH_ecm240.root'
+
 
     # Mandatory: analyzers function to define the analysis graph, please make
     # sure you return the dataframe, in this example it is dframe2
@@ -70,7 +72,7 @@ def analyzers(self, dframe):
         dframe2 = (
             dframe
             # define an alias for muon index collection
-            .Alias('Muon0', 'Muon#0.index')
+            .Alias('Muon0', 'Muon_objIdx.index')
             # define the muon collection
             .Define(
                 'muons',

examples/FCCee/higgs/mH-recoil/plots_recoil.py

@@ -0,0 +1,107 @@
+import ROOT
+
+# global parameters
+intLumi        = 1.
+intLumiLabel   = "L = 5 ab^{-1}"
+ana_tex        = 'e^{+}e^{-} #rightarrow ZH #rightarrow #mu^{+}#mu^{-} + X'
+delphesVersion = '3.4.2'
+energy         = 240.0
+collider       = 'FCC-ee'
+formats        = ['png','pdf']
+
+outdir         = './histmaker_outputs/' 
+inputDir       = './histmaker_outputs/' 
+
+plotStatUnc    = True
+
+colors = {}
+colors['ZH'] = ROOT.kRed
+colors['WW'] = ROOT.kBlue+1
+colors['ZZ'] = ROOT.kGreen+2
+
+#procs = {}
+#procs['signal'] = {'ZH':['wzp6_ee_mumuH_ecm240']}
+#procs['backgrounds'] =  {'WW':['p8_ee_WW_ecm240'], 'ZZ':['p8_ee_ZZ_ecm240']}
+procs = {}
+procs['signal'] = {'ZH':['p8_ee_ZH_ecm240']}
+# procs['signal'] = {'ZH':['p8_ee_ZH_Zmumu_ecm240']}
+#procs['backgrounds'] =  {'WW':['p8_ee_WW_mumu_ecm240'], 'ZZ':['p8_ee_ZZ_mumubb_ecm240']}
+procs['backgrounds'] =  {'ZZ':['p8_ee_ZZ_ecm240']}
+# procs['backgrounds'] =  {'ZZ':['p8_ee_ZZ_mumubb_ecm240']}
+
+
+legend = {}
+legend['ZH'] = 'ZH'
+legend['WW'] = 'WW'
+legend['ZZ'] = 'ZZ'
+
+
+
+hists = {}
+
+hists["zmumu_recoil_m"] = {
+    "output":   "zmumu_recoil_m",
+    "logy":     False,
+    "stack":    True,
+    "rebin":    100,
+    "xmin":     120,
+    "xmax":     140,
+    "ymin":     0,
+    "ymax":     2500,
+    "xtitle":   "Recoil (GeV)",
+    "ytitle":   "Events / 100 MeV",
+}
+
+hists["zmumu_p"] = {
+    "output":   "zmumu_p",
+    "logy":     False,
+    "stack":    True,
+    "rebin":    2,
+    "xmin":     0,
+    "xmax":     80,
+    "ymin":     0,
+    "ymax":     2000,
+    "xtitle":   "p(#mu^{#plus}#mu^{#minus}) (GeV)",
+    "ytitle":   "Events ",
+}
+
+hists["zmumu_m"] = {
+    "output":   "zmumu_m",
+    "logy":     False,
+    "stack":    True,
+    "rebin":    2,
+    "xmin":     86,
+    "xmax":     96,
+    "ymin":     0,
+    "ymax":     3000,
+    "xtitle":   "m(#mu^{#plus}#mu^{#minus}) (GeV)",
+    "ytitle":   "Events ",
+}
+
+hists["cosThetaMiss_cut4"] = {
+    "output":   "cosThetaMiss_cut4",
+    "logy":     True,
+    "stack":    True,
+    "rebin":    10,
+    "xmin":     0,
+    "xmax":     1,
+    "ymin":     10,
+    "ymax":     100000,
+    "xtitle":   "cos(#theta_{miss})",
+    "ytitle":   "Events ",
+    "extralab": "Before cos(#theta_{miss}) cut",
+}
+
+
+# hists["cutFlow"] = {
+#     "output":   "cutFlow",
+#     "logy":     True,
+#     "stack":    True,
+#     "xmin":     0,
+#     "xmax":     6,
+#     "ymin":     1e4,
+#     "ymax":     1e11,
+#     "xtitle":   ["All events", "#geq 1 #mu^{#pm} + ISO", "#geq 2 #mu^{#pm} + OS", "86 < m_{#mu^{+}#mu^{#minus}} < 96", "20 < p_{#mu^{+}#mu^{#minus}} < 70", "|cos#theta_{miss}| < 0.98", "120 < m_{rec} < 140"],
+#     "ytitle":   "Events ",
+#     "scaleSig": 10
+# }