Skip to content

Commit

Permalink
Updating FCC-ee mH recoil mumu example to work with edm4hep v1 format (
Browse files Browse the repository at this point in the history
…#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 <[email protected]>
  • Loading branch information
bistapf and Birgit Stapf authored Sep 12, 2024
1 parent fa672d4 commit 54e8495
Show file tree
Hide file tree
Showing 4 changed files with 310 additions and 16 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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'
185 changes: 185 additions & 0 deletions examples/FCCee/higgs/mH-recoil/histmaker_recoil.py
Original file line number Diff line number Diff line change
@@ -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
32 changes: 17 additions & 15 deletions examples/FCCee/higgs/mH-recoil/mumu/analysis_stage1.py
Original file line number Diff line number Diff line change
Expand Up @@ -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/' \
Expand All @@ -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
Expand All @@ -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',
Expand Down
107 changes: 107 additions & 0 deletions examples/FCCee/higgs/mH-recoil/plots_recoil.py
Original file line number Diff line number Diff line change
@@ -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
# }

0 comments on commit 54e8495

Please sign in to comment.