From a0ff223f6d80c370c1e60705423cc3a054eb3ee5 Mon Sep 17 00:00:00 2001 From: Birgit Stapf Date: Fri, 20 Sep 2024 17:39:42 +0200 Subject: [PATCH] first versions of FCChh example final and plotting scripts --- examples/FCChh/ggHH_bbyy/analysis_final.py | 39 ++++ examples/FCChh/ggHH_bbyy/analysis_plots.py | 41 +++++ examples/FCChh/ggHH_bbyy/analysis_stage1.py | 53 +++--- .../ggHH_bbyy/analysis_stage1_fullfile.py | 173 ++++++++++++++++++ 4 files changed, 276 insertions(+), 30 deletions(-) create mode 100644 examples/FCChh/ggHH_bbyy/analysis_final.py create mode 100644 examples/FCChh/ggHH_bbyy/analysis_plots.py create mode 100644 examples/FCChh/ggHH_bbyy/analysis_stage1_fullfile.py diff --git a/examples/FCChh/ggHH_bbyy/analysis_final.py b/examples/FCChh/ggHH_bbyy/analysis_final.py new file mode 100644 index 0000000000..fd62d5345c --- /dev/null +++ b/examples/FCChh/ggHH_bbyy/analysis_final.py @@ -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 /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 +} \ No newline at end of file diff --git a/examples/FCChh/ggHH_bbyy/analysis_plots.py b/examples/FCChh/ggHH_bbyy/analysis_plots.py new file mode 100644 index 0000000000..0d87651332 --- /dev/null +++ b/examples/FCChh/ggHH_bbyy/analysis_plots.py @@ -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' \ No newline at end of file diff --git a/examples/FCChh/ggHH_bbyy/analysis_stage1.py b/examples/FCChh/ggHH_bbyy/analysis_stage1.py index 845301345e..23eabc1cfe 100644 --- a/examples/FCChh/ggHH_bbyy/analysis_stage1.py +++ b/examples/FCChh/ggHH_bbyy/analysis_stage1.py @@ -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' @@ -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]") @@ -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)") @@ -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)") @@ -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 @@ -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 diff --git a/examples/FCChh/ggHH_bbyy/analysis_stage1_fullfile.py b/examples/FCChh/ggHH_bbyy/analysis_stage1_fullfile.py new file mode 100644 index 0000000000..16d0810216 --- /dev/null +++ b/examples/FCChh/ggHH_bbyy/analysis_stage1_fullfile.py @@ -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: + ## '':{'fraction':, 'chunks':, 'output': }, + # # - 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