From 54e849554155d86eeb0eb70d2c5ea7b753fb1dba Mon Sep 17 00:00:00 2001 From: bistapf <79460851+bistapf@users.noreply.github.com> Date: Thu, 12 Sep 2024 14:55:44 +0200 Subject: [PATCH] Updating FCC-ee mH recoil mumu example to work with edm4hep v1 format (#400) * updating mH recoil mumu example to work with edm4hep v1, adding also histmaker and plots scripts as used in FCC SW tutorial * fix test file needed by FCCAnalysis_analysis_example_run_analysis * adapt build&test to work with edm4hep v1 format with temporary changes/ommissions to be able to develop further on this branch * fix comments in test * cleanup/changes from review of PR#400 * add missing end quote --------- Co-authored-by: Birgit Stapf --- .github/workflows/test.yml | 2 +- .../FCCee/higgs/mH-recoil/histmaker_recoil.py | 185 ++++++++++++++++++ .../higgs/mH-recoil/mumu/analysis_stage1.py | 32 +-- .../FCCee/higgs/mH-recoil/plots_recoil.py | 107 ++++++++++ 4 files changed, 310 insertions(+), 16 deletions(-) create mode 100644 examples/FCCee/higgs/mH-recoil/histmaker_recoil.py create mode 100644 examples/FCCee/higgs/mH-recoil/plots_recoil.py diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 4b0fad4ffc..0fc6706602 100644 --- a/.github/workflows/test.yml +++ b/.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' \ No newline at end of file diff --git a/examples/FCCee/higgs/mH-recoil/histmaker_recoil.py b/examples/FCCee/higgs/mH-recoil/histmaker_recoil.py new file mode 100644 index 0000000000..92969c5c51 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py b/examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py index 1b0baedb7f..fd1d9f6949 100644 --- a/examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py +++ b/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 - # /p8_ee_ZZ_ecm240.root - 'p8_ee_ZZ_ecm240': {'fraction': 0.005}, - # Run 50% of the statistics with output into two files named - # /p8_ee_WW_ecm240/chunk.root - 'p8_ee_WW_ecm240': {'fraction': 0.5, 'chunks': 2}, - # Run 20% of the statistics in one file named - # /p8_ee_ZH_ecm240_out.root (example on how to change - # the output name) + # # Run the full statistics in one output file named + # # /p8_ee_ZZ_ecm240.root + 'p8_ee_ZZ_ecm240': {'fraction': 1.}, + # # Run 50% of the statistics with output into two files named + # # /p8_ee_WW_ecm240/chunk.root + 'p8_ee_WW_ecm240': {'fraction': 0.5, 'chunks': 2}, #this doesn't work? + # # Run 20% of the statistics in one file named + # # /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', diff --git a/examples/FCCee/higgs/mH-recoil/plots_recoil.py b/examples/FCCee/higgs/mH-recoil/plots_recoil.py new file mode 100644 index 0000000000..4586368467 --- /dev/null +++ b/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 +# } \ No newline at end of file