Skip to content

Commit

Permalink
first versions of FCChh example final and plotting scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
bistapf committed Sep 20, 2024
1 parent a178c42 commit a0ff223
Show file tree
Hide file tree
Showing 4 changed files with 276 additions and 30 deletions.
39 changes: 39 additions & 0 deletions examples/FCChh/ggHH_bbyy/analysis_final.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
#Input directory where the files produced at the pre-selection level are
inputDir = "outputs/FCChh/ggHH_bbyy/presel/"

#Input directory where the files produced at the pre-selection level are
outputDir = "outputs/FCChh/ggHH_bbyy/final/"

processList = {
'pwp8_pp_hh_5f_hhbbyy':{},#Run over the full statistics from stage2 input file <inputDir>/p8_ee_ZZ_ecm240.root. Keep the same output name as input
}

#Link to the dictonary that contains all the cross section informations etc...
procDict = "FCCee_procDict_spring2021_IDEA.json" # will need an FCC-hh one!

#Add MySample_p8_ee_ZH_ecm240 as it is not an offical process TO UPDATE
procDictAdd={"pwp8_pp_hh_5f_hhbbyy":{"numberOfEvents": 10000000, "sumOfWeights": 10000000, "crossSection": 1.0, "kfactor": 1.0, "matchingEfficiency": 1.0}}

#Number of CPUs to use
nCPUS = 2

#produces ROOT TTrees, default is False
doTree = True

saveTabular = True

###Dictionnay of the list of cuts. The key is the name of the selection that will be added to the output file
cutList = {
"sel0_myy":"m_yy[0] > 100. && m_yy[0] < 180.",
"sel1_mbb":"(m_yy[0] > 100. && m_yy[0] < 180.) && (m_bb[0] > 80. && m_bb[0] < 200.)",
}


#Dictionary for the ouput variable/hitograms. The key is the name of the variable in the output files. "name" is the name of the variable in the input file, "title" is the x-axis label of the histogram, "bin" the number of bins of the histogram, "xmin" the minimum x-axis value and "xmax" the maximum x-axis value.
histoList = {
"myy":{"name":"m_yy","title":"m_{#gamma#gamma} [GeV]","bin":50,"xmin":0,"xmax":200},
"myy_zoom":{"name":"m_yy","title":"m_{#gamma#gamma} [GeV]","bin":50,"xmin":100,"xmax":180},
"y1_pT":{"name":"g1_pt","title":"pT_{#gamma1} [GeV]","bin":50,"xmin":0.,"xmax":200.},
"y2_pT":{"name":"g2_pt","title":"pT_{#gamma2} [GeV]","bin":50,"xmin":0.,"xmax":200.},
"pT_y1_vs_y2_2D":{"cols":["g1_pt", "g2_pt"],"title":"m_{Z} - leptonic recoil [GeV]", "bins": [(40,80,100), (100,120,140)]}, # 2D histogram
}
41 changes: 41 additions & 0 deletions examples/FCChh/ggHH_bbyy/analysis_plots.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
import ROOT

# global parameters
intLumi = 30e+06 #in pb-1
ana_tex = 'pp #rightarrow HH #rightarrow b #bar{b} #gamma #gamma '
delphesVersion = '3.4.2'
energy = 100
collider = 'FCC-hh'
inputDir = 'outputs/FCChh/ggHH_bbyy/final/'
formats = ['png','pdf']
yaxis = ['lin','log']
stacksig = ['nostack']
# stacksig = ['stack','nostack']
outdir = 'outputs/FCChh/ggHH_bbyy/plots/'
plotStatUnc = True

variables = ['myy','myy_zoom','y1_pT','y2_pT']
# variables = ['myy','myy_zoom','y1_pT','y2_pT','pT_y1_vs_y2_2D']
# rebin = [1, 1, 1, 1, 2] # uniform rebin per variable (optional)

### Dictionary with the analysis name as a key, and the list of selections to be plotted for this analysis. The name of the selections should be the same than in the final selection
selections = {}
selections['bbyy_analysis'] = ["sel0_myy","sel1_mbb"]

extralabel = {}
extralabel['sel0_myy'] = "Sel.: 100 < m_{#gamma#gamma} < 180 GeV"
extralabel['sel1_mbb'] = "Sel.: 100 < m_{#gamma#gamma} < 180 GeV and 80 < m_{bb} < 200 GeV"

colors = {}
colors['bbyy_signal'] = ROOT.kRed
# colors['WW'] = ROOT.kBlue+1
# colors['ZZ'] = ROOT.kGreen+2
# colors['VV'] = ROOT.kGreen+3

plots = {}
plots['bbyy_analysis'] = {'signal':{'bbyy_signal':['pwp8_pp_hh_5f_hhbbyy']},
}


legend = {}
legend['bbyy_signal'] = 'HH signal'
53 changes: 23 additions & 30 deletions examples/FCChh/ggHH_bbyy/analysis_stage1.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,7 @@ def __init__(self, cmdline_args):
# self.prod_tag = 'FCCee/spring2021/IDEA/'

# Optional: output directory, default is local running directory
self.output_dir = 'outputs/FCChh/ggHH_bbyy/' \
# f'stage1_{self.ana_args.muon_pt}'
self.output_dir = 'outputs/FCChh/ggHH_bbyy/presel/'

# Optional: analysisName, default is ''
# self.analysis_name = 'My Analysis'
Expand Down Expand Up @@ -72,7 +71,7 @@ def analyzers(self, dframe):
.Define("gamma", "FCCAnalyses::ReconstructedParticle::get(Photon_objIdx.index, ReconstructedParticles)")
.Define("selpt_gamma", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(gamma)")
.Define("sel_gamma_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_gamma)")
.Define("sel_gamma", "AnalysisFCChh::SortParticleCollection(sel_gamma_unsort)")
.Define("sel_gamma", "AnalysisFCChh::SortParticleCollection(sel_gamma_unsort)") #sort by pT

.Define("ngamma", "FCCAnalyses::ReconstructedParticle::get_n(sel_gamma)")
.Define("g1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_gamma)[0]")
Expand All @@ -85,16 +84,14 @@ def analyzers(self, dframe):
.Define("g2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_gamma)[1]")

#H(yy) if it exists, if there are no 2 selected photons, doesnt get filled
.Define("yy_pairs_unmerged", "AnalysisFCChh::getPairs(sel_gamma)")
.Define("yy_pairs", "AnalysisFCChh::merge_pairs(yy_pairs_unmerged)")
.Define("yy_pairs_unmerged", "AnalysisFCChh::getPairs(sel_gamma)") # retrieves the leading pT pair of all possible
.Define("yy_pairs", "AnalysisFCChh::merge_pairs(yy_pairs_unmerged)") # merge pair into one object to access inv masses etc
.Define("m_yy", "FCCAnalyses::ReconstructedParticle::get_mass(yy_pairs)")


########################################### JETS ###########################################

# jets after overlap removal is performed between jets and isolated electrons, muons and photons

# selected jets above a pT threshold of 30 GeV, eta < 4, tight ID
# selected jets above a pT threshold of 30 GeV, eta < 4
.Define("selpt_jets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(Jet)")
.Define("sel_jets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_jets)")
.Define("sel_jets", "AnalysisFCChh::SortParticleCollection(sel_jets_unsort)")
Expand All @@ -110,12 +107,8 @@ def analyzers(self, dframe):
.Define("j2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[1]")

#b-tagged jets:
# .Alias("Jet3","Jet#3.index")
.Alias("Jet3","_Jet_particles.index")
#LOOSE WP : bit 1 = medium WP, bit 2 = tight WP
#b tagged jets
.Define("bjets", "AnalysisFCChh::get_tagged_jets(ReconstructedParticles, ParticleIDs, _ParticleIDs_particle, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl
# .Define("bjets", "AnalysisFCChh::get_tagged_jets(Jet, Jet3, ParticleIDs, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl
.Define("bjets", "AnalysisFCChh::get_tagged_jets(ReconstructedParticles, ParticleIDs, _ParticleIDs_particle, _ParticleIDs_parameters, 1)") #bit 1 = medium WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl
.Define("selpt_bjets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(bjets)")
.Define("sel_bjets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_bjets)")
.Define("sel_bjets", "AnalysisFCChh::SortParticleCollection(sel_bjets_unsort)")
Expand All @@ -129,12 +122,20 @@ def analyzers(self, dframe):
.Define("b2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_bjets)[1]")
.Define("b2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_bjets)[1]")

#H(bb) system - using the medium WP jets - if there are no 2 b-tagged jets these variable don't get filled
.Define("bb_pairs_unmerged", "AnalysisFCChh::getPairs(sel_bjets)") # retrieves the leading pT pair of all possible
.Define("bb_pairs", "AnalysisFCChh::merge_pairs(bb_pairs_unmerged)") # merge pair into one object to access inv masses etc
.Define("m_bb", "FCCAnalyses::ReconstructedParticle::get_mass(bb_pairs)")

########################################### APPLY PRE-SELECTION ###########################################

# #H(bb) system - using the loose WP b-jets -> if the events has < 2 bjets these variables do not get filled!
# .Define("bb_pairs_unmerged", "AnalysisFCChh::getPairs(sel_bjets)") #currently gets only leading pT pair, as a RecoParticlePair
# #then merge the bb pair into one object and get its kinematic properties
# .Define("bb_pairs", "AnalysisFCChh::merge_pairs(bb_pairs_unmerged)") #merge into one object to access inv masses etc
# .Define("m_bb", "FCCAnalyses::ReconstructedParticle::get_mass(bb_pairs)")
#require at least two b-jets and two photons, both with invariant masses compatible with the Higgs mass
.Filter("sel_bjets.size()>1")
.Filter("sel_gamma.size()>1")
# .Filter("m_bb[0] < 200.")
# .Filter("m_bb[0] > 80.")
# .Filter("m_yy[0] < 180.")
# .Filter("m_yy[0] > 100.")

)
return dframe2
Expand All @@ -146,17 +147,9 @@ def output(self):
Output variables which will be saved to output root file.
'''
branch_list = [
'nbjets',
# 'm_bb',
'ngamma',
'm_yy',
# 'selected_muons_pt',
# 'selected_muons_y',
# 'selected_muons_p',
# 'selected_muons_e',
# 'zed_leptonic_pt',
# 'zed_leptonic_m',
# 'zed_leptonic_charge',
# 'zed_leptonic_recoil_m'
# Photons and H(yy) system:
'ngamma', 'g1_pt', 'g2_pt', 'g1_eta', 'g2_eta', 'm_yy',
# b-jets and H(bb) system:
'nbjets', 'b1_pt', 'b2_pt', 'b1_eta', 'b2_eta', 'm_bb',
]
return branch_list
173 changes: 173 additions & 0 deletions examples/FCChh/ggHH_bbyy/analysis_stage1_fullfile.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
'''
Analysis example for FCC-hh, using gg->HH->bbyy di-Higgs production events
'''
from argparse import ArgumentParser


# Mandatory: Analysis class where the user defines the operations on the
# dataframe.
class Analysis():
'''
Di-Higgs analysis in bbyy.
'''
def __init__(self, cmdline_args):
parser = ArgumentParser(
description='Additional analysis arguments',
usage='Provide additional arguments after analysis script path')
# parser.add_argument('--bjet-pt', default='10.', type=float,
# help='Minimal pT of the selected b-jets.')
# Parse additional arguments not known to the FCCAnalyses parsers
# All command line arguments know to fccanalysis are provided in the
# `cmdline_arg` dictionary.
self.ana_args, _ = parser.parse_known_args(cmdline_args['unknown'])

# Mandatory: List of processes to run over
self.process_list = {
# # Add your processes like this:
## '<name of process>':{'fraction':<fraction of events to run over>, 'chunks':<number of chunks to split the output into>, 'output':<name of the output file> },
# # - <name of process> needs to correspond either the name of the input .root file, or the name of a directory containing root files
# # If you want to process only part of the events, split the output into chunks or give a different name to the output use the optional arguments
# # or leave blank to use defaults = run the full statistics in one output file named the same as the process:
'pwp8_pp_hh_lambda100_5f_hhbbaa': {},
}

# Mandatory: Input directory where to find the samples, or a production tag when running over the centrally produced
# samples (this points to the yaml files for getting sample statistics)
self.input_dir = '/eos/experiment/fcc/hh/generation/DelphesEvents/fcc_v06_scenarioI/'
# self.prod_tag = 'FCCee/spring2021/IDEA/'

# Optional: output directory, default is local running directory
self.output_dir = 'fcc_v06_bbyy/' \
# f'stage1_{self.ana_args.muon_pt}'

# Optional: analysisName, default is ''
# self.analysis_name = 'My Analysis'

# Optional: number of threads to run on, default is 'all available'
# self.n_threads = 4

# Optional: running on HTCondor, default is False
# self.run_batch = False

# Optional: test file that is used if you run with the --test argument (fccanalysis run ./examples/FCChh/ggHH_bbyy/analysis_stage1.py --test)
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
def analyzers(self, dframe):
'''
Analysis graph.
'''

dframe2 = (
dframe

.Define("weight", "EventHeader.weight")

########################################### PHOTONS ###########################################

.Define("gamma", "FCCAnalyses::ReconstructedParticle::get(Photon_objIdx.index, ReconstructedParticles)")
.Define("selpt_gamma", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(gamma)")
.Define("sel_gamma_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_gamma)")
.Define("sel_gamma", "AnalysisFCChh::SortParticleCollection(sel_gamma_unsort)")

.Define("ngamma", "FCCAnalyses::ReconstructedParticle::get_n(sel_gamma)")
.Define("g1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_gamma)[0]")
.Define("g1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_gamma)[0]")
.Define("g1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_gamma)[0]")
.Define("g1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_gamma)[0]")
.Define("g2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_gamma)[1]")
.Define("g2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_gamma)[1]")
.Define("g2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_gamma)[1]")
.Define("g2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_gamma)[1]")

#H(yy) if it exists, if there are no 2 selected photons, doesnt get filled
.Define("yy_pairs_unmerged", "AnalysisFCChh::getPairs(sel_gamma)")
.Define("yy_pairs", "AnalysisFCChh::merge_pairs(yy_pairs_unmerged)")
.Define("m_yy", "FCCAnalyses::ReconstructedParticle::get_mass(yy_pairs)")


########################################### JETS ###########################################

# jets after overlap removal is performed between jets and isolated electrons, muons and photons

# selected jets above a pT threshold of 30 GeV, eta < 4, tight ID
.Define("selpt_jets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(Jet)")
.Define("sel_jets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_jets)")
.Define("sel_jets", "AnalysisFCChh::SortParticleCollection(sel_jets_unsort)")
.Define("njets", "FCCAnalyses::ReconstructedParticle::get_n(sel_jets)")
.Define("j1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[0]")
.Define("j1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[0]")
.Define("j1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[0]")
.Define("j1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[0]")

.Define("j2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_jets)[1]")
.Define("j2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_jets)[1]")
.Define("j2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_jets)[1]")
.Define("j2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_jets)[1]")

#b-tagged jets:
# .Alias("Jet3","Jet#3.index")
.Alias("Jet3","_Jet_particles.index")
#LOOSE WP : bit 1 = medium WP, bit 2 = tight WP
#b tagged jets
.Define("bjets", "AnalysisFCChh::get_tagged_jets(ReconstructedParticles, ParticleIDs, _ParticleIDs_particle, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl
# .Define("bjets", "AnalysisFCChh::get_tagged_jets(Jet, Jet3, ParticleIDs, _ParticleIDs_parameters, 1)") #bit 0 = loose WP, see: https://github.com/delphes/delphes/blob/master/cards/FCC/scenarios/FCChh_I.tcl
# .Define("selpt_bjets", "FCCAnalyses::ReconstructedParticle::sel_pt(30.)(bjets)")
# .Define("sel_bjets_unsort", "FCCAnalyses::ReconstructedParticle::sel_eta(4)(selpt_bjets)")
# .Define("sel_bjets", "AnalysisFCChh::SortParticleCollection(sel_bjets_unsort)")
.Alias("sel_bjets", "bjets")
.Define("nbjets", "FCCAnalyses::ReconstructedParticle::get_n(sel_bjets)")
.Define("b1_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_bjets)[0]")
.Define("b1_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_bjets)[0]")
.Define("b1_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_bjets)[0]")
.Define("b1_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_bjets)[0]")
.Define("b2_e", "FCCAnalyses::ReconstructedParticle::get_e(sel_bjets)[1]")
.Define("b2_pt", "FCCAnalyses::ReconstructedParticle::get_pt(sel_bjets)[1]")
.Define("b2_eta", "FCCAnalyses::ReconstructedParticle::get_eta(sel_bjets)[1]")
.Define("b2_phi", "FCCAnalyses::ReconstructedParticle::get_phi(sel_bjets)[1]")


# #H(bb) system - using the loose WP b-jets -> if the events has < 2 bjets these variables do not get filled!
.Define("bb_pairs_unmerged", "AnalysisFCChh::getPairs(sel_bjets)") #currently gets only leading pT pair, as a RecoParticlePair
#then merge the bb pair into one object and get its kinematic properties
.Define("bb_pairs", "AnalysisFCChh::merge_pairs(bb_pairs_unmerged)") #merge into one object to access inv masses etc
.Define("m_bb", "FCCAnalyses::ReconstructedParticle::get_mass(bb_pairs)")

# Filter at least one candidate
.Filter("sel_bjets.size()>1")
# .Filter("sel_gamma.size()>1")
# .Filter("m_bb[0] < 200.")
# .Filter("m_bb[0] > 80.")
# .Filter("m_yy[0] < 180.")
# .Filter("m_yy[0] > 100.")

)
return dframe2

# Mandatory: output function, please make sure you return the branch list
# as a python list
def output(self):
'''
Output variables which will be saved to output root file.
'''
branch_list = [
'nbjets',
# 'm_bb',
'ngamma',
'm_yy',
'm_bb',
'b1_pt', 'b2_pt', 'b1_eta', 'b2_eta',
# 'selected_muons_pt',
# 'selected_muons_y',
# 'selected_muons_p',
# 'selected_muons_e',
# 'zed_leptonic_pt',
# 'zed_leptonic_m',
# 'zed_leptonic_charge',
# 'zed_leptonic_recoil_m'
]
return branch_list

0 comments on commit a0ff223

Please sign in to comment.