57 #define CHECK( ARG ) \
59 const bool result = ARG; \
61 printf("Failed to execute: %s \n", #ARG);\
70 printf(
"Running options:\n");
71 printf(
"YOU HAVE TO ADAPT THE OPTIONS TO FFJETSMEARINGCORRECTION\n");
72 printf(
" --help : To get the help you're reading\n");
73 printf(
" --jetColl= : Specify the jet collection (TrimmedLargeR)\n");
74 printf(
" --MassDef= : Specify the kind of jet mass used (Calo, TA, Comb)\n");
75 printf(
" --sample= : Specify input xAOD\n");
76 printf(
" Example: FFJetSmearingTool_MyExample --truth_jetColl=AntiKt10TruthTrimmedPtFrac5SmallR20Jets --reco_jetColl=AntiKt10LCTopoTrimmedPtFrac5SmallR20Jets --MassDef=Comb --sample=DAOD_JETM6.16317590._000003.pool.root.1 --output=file_name.root --ConfigFile=rel21/Spring2020/FFJetSmearingTool_TestOnly_JMR.config --DebugTool=true \n");
96 std::string truth_jetColl =
"";
97 std::string reco_jetColl =
"";
98 std::string string_kindofmass =
"";
100 std::string string_configfile_path =
"";
101 std::string string_debugtool =
"";
108 std::string
opt(
argv[
i]); std::vector< std::string >
v;
110 std::istringstream iss(
opt);
115 while (std::getline(iss,
item, delim)){
119 if (
opt.find(
"--help") != std::string::npos ) {
123 if (
opt.find(
"--sample=") != std::string::npos )
sample =
v[1];
124 if (
opt.find(
"--truth_jetColl=") != std::string::npos ) truth_jetColl =
v[1];
125 if (
opt.find(
"--reco_jetColl=") != std::string::npos ) reco_jetColl =
v[1];
126 if (
opt.find(
"--MassDef=") != std::string::npos ) string_kindofmass =
v[1];
127 if (
opt.find(
"--ConfigFile=") != std::string::npos ) string_configfile_path =
v[1];
128 if (
opt.find(
"--DebugTool=") != std::string::npos ) string_debugtool =
v[1];
129 if (
opt.find(
"--output=") != std::string::npos )
output =
v[1];
134 printf(
"No input xAOD file specified, exiting\n");
137 if(truth_jetColl==
""){
138 printf(
"No truth jet collection specified, exiting\n");
141 if(reco_jetColl==
""){
142 printf(
"No truth jet collection specified, exiting\n");
145 if(string_kindofmass==
""){
146 printf(
"No kind of jet mass specified, exiting\n");
149 if(string_configfile_path==
""){
150 printf(
"No ConfigFile specified, exiting\n");
153 if(string_debugtool==
""){
154 printf(
"No debugtool specified, exiting\n");
158 printf(
"output not specified, exiting\n");
161 printf(
"%s %s %s %s %s \n",
sample.c_str(), truth_jetColl.c_str(), reco_jetColl.c_str(),
output.c_str(), string_debugtool.c_str());
164 TString kindofmass = string_kindofmass;
166 bool want_to_debug =
false;
168 if(string_debugtool ==
"true"){want_to_debug =
true;}
169 if(string_debugtool ==
"false"){want_to_debug =
false;}
178 std::unique_ptr< TFile >
ifile( TFile::Open(
sample.c_str(),
"READ" ) );
191 const std::string name_FFJetSmearingTool =
"FFJetSmearing_Example";
212 (ffjetsmearingtool.setProperty(
"MassDef", kindofmass )).
ignore();
214 (ffjetsmearingtool.setProperty(
"ConfigFile", string_configfile_path )).
ignore();
217 if(!(ffjetsmearingtool.
initialize().isSuccess())){
218 printf(
"Initialization of FFJetSmearingTool failed, exiting\n");
224 ffjetsmearingtool.setProperty (
"OutputLevel", MSG::ERROR);
227 ffjetsmearingtool.setProperty (
"OutputLevel",
MSG::VERBOSE);
233 printf(
"Recommended systematics: \n");
234 for(
auto sysItr = recommendedSysts.
begin();
235 sysItr != recommendedSysts.
end(); ++sysItr){
236 printf(
"%s \n", sysItr->name().c_str());
240 std::vector<CP::SystematicSet> sysList;
249 CHECK(m_JetTruthLabelingTool.setProperty(
"TruthLabelName",
"R10TruthLabel_R21Consolidated"));
250 CHECK(m_JetTruthLabelingTool.setProperty(
"UseTRUTH3",
false));
251 CHECK(m_JetTruthLabelingTool.setProperty(
"TruthParticleContainerName",
"TruthParticles"));
261 const std::string name_JetCalibTools =
"JetCalibTool_LC";
266 jetCalibrationTool.setProperty(
"JetCollection",
"AntiKt10LCTopoTrimmedPtFrac5SmallR20");
268 if(kindofmass==
"Comb"){
269 jetCalibrationTool.setProperty(
"ConfigFile",
"JES_MC16recommendation_FatJet_Trimmed_JMS_comb_17Oct2018.config");
271 if(kindofmass==
"Calo"){
272 jetCalibrationTool.setProperty(
"ConfigFile",
"JES_MC16recommendation_FatJet_Trimmed_JMS_calo_12Oct2018.config");
274 if(kindofmass==
"TA"){
275 jetCalibrationTool.setProperty(
"ConfigFile",
"JES_MC16recommendation_FatJet_Trimmed_JMS_TA_12Oct2018.config");
279 jetCalibrationTool.setProperty(
"CalibSequence",
"EtaJES_JMS");
280 jetCalibrationTool.setProperty(
"IsData",
false);
281 jetCalibrationTool.setProperty(
"CalibArea",
"00-04-82");
282 jetCalibrationTool.setProperty(
"DEVmode",
false);
284 if(!(jetCalibrationTool.
initialize().isSuccess())){
285 printf(
"Initialization of JetCalibTools failed, exiting \n");
293 printf(
"\n=============SYSTEMATICS CHECK NOW \n");
294 for (
auto sys : sysList) {
297 printf(
"Error, Cannot configure calibration tool for systematics \n");
300 printf(
"\nWe are using the systematic %s \n",
sys.name().c_str());
307 static constexpr
float MeVtoGeV = 1.e-3;
313 TH1F* reco_jet_mass_hist=
new TH1F(
"reco_jet_mass_hist",
"reco_jet_mass_hist", 100, 0, 200);
314 TH1F* matched_truth_jet_mass_hist=
new TH1F(
"matched_truth_jet_mass_hist",
"matched_truth_jet_mass_hist", 100, 0, 200);
315 TH1F* smeared_reco_jet_mass_hist=
new TH1F(
"smeared_reco_jet_mass_hist",
"smeared_reco_jet_mass_hist", 100, 0, 200);
317 reco_jet_mass_hist->Sumw2();
318 matched_truth_jet_mass_hist->Sumw2();
319 smeared_reco_jet_mass_hist->Sumw2();
326 int upperlimit1 = 600;
327 int upperlimit2 = 1000;
329 int numBinsMass = 120;
336 std::unique_ptr<TH3F> hist_jet_mass_scale_change_3D;
337 hist_jet_mass_scale_change_3D = std::make_unique<TH3F>(
"hist_jet_mass_scale_change_3D",
"hist_jet_mass_scale_change_3D",numBinsPt,0,upperlimit2,numBinsMass,0,upperlimit1,numBinsMass,0,upperlimit1);
340 float lowerlimit3 = -0.5;
341 float upperlimit3 = 0.5;
342 int numBinsDiff = 100;
345 std::unique_ptr<TH3F> hist_jet_mass_resolution_change_3D;
346 hist_jet_mass_resolution_change_3D = std::make_unique<TH3F>(
"hist_jet_mass_resolution_change_3D",
"hist_jet_mass_resolution_change_3D",numBinsPt,0,upperlimit2,numBinsMass,0,upperlimit1,numBinsDiff,lowerlimit3,upperlimit3);
355 const Long64_t
nevents =
event.getEntries();
362 for(Long64_t ievent = 0; ievent <
nevents; ++ievent){
366 printf(
"Failed to load entry %lld \n", ievent);
372 if(ievent % 100==0) printf(
"Event %lld of %lld \n", ievent,
nevents);
380 printf(
"===>>> start processing event #%llu, run #%i, %lld events processed so far <<<=== \n", ei->
eventNumber(), ei->
runNumber(), ievent);
387 printf(
"Number of truth jets: %lu \n", jets_truth->
size());
390 for(
const xAOD::Jet* jet_truth : *jets_truth){
392 printf(
"Truth Jet: pt = %f, mass = %f, eta = %f \n", jet_truth->pt()*
MeVtoGeV, jet_truth->m()*
MeVtoGeV, jet_truth->eta());
401 printf(
"Number of reco jets: %lu \n", jets_reco->
size());
404 for(
const xAOD::Jet* jet_reco : *jets_reco){
407 printf(
"Reco Jet: pt = %f, mass = %f, eta = %f \n", jet_reco->pt()*
MeVtoGeV, jet_reco->m()*
MeVtoGeV, jet_reco->eta());
425 m_JetTruthLabelingTool.
modify(*(jets_shallowCopy.first));
427 if(want_to_debug) printf(
"Start the loop over the jets \n");
429 bool lead_jet =
true;
432 for(
xAOD::Jet* jet_reco : *( jets_shallowCopy.first ) ) {
443 jet_reco->getAttribute<
xAOD::JetFourMom_t>(
"JetJMSScaleMomentumCombQCD",jet_reco_Comb_FourMom);
444 jet_reco->getAttribute<
xAOD::JetFourMom_t>(
"JetJMSScaleMomentumCalo",jet_reco_CALO_FourMom);
445 jet_reco->getAttribute<
xAOD::JetFourMom_t>(
"JetJMSScaleMomentumTA",jet_reco_TA_FourMom);
447 if(want_to_debug && kindofmass ==
"Comb"){
448 printf(
"Comb jet mass = %f \n", jet_reco_Comb_FourMom.mass());
449 printf(
"Calo jet mass = %f \n", jet_reco_CALO_FourMom.mass());
450 printf(
"TA jet mass = %f \n", jet_reco_TA_FourMom.mass());
453 std::unique_ptr<xAOD::Jet> jet_reco_Comb = std::make_unique<xAOD::Jet>();
456 jet_reco_Comb->
setJetP4(jet_reco_Comb_FourMom);
457 std::unique_ptr<xAOD::Jet> jet_reco_CALO = std::make_unique<xAOD::Jet>();
459 jet_reco_CALO->
setJetP4(jet_reco_CALO_FourMom);
460 std::unique_ptr<xAOD::Jet> jet_reco_TA = std::make_unique<xAOD::Jet>();
462 jet_reco_TA->
setJetP4(jet_reco_TA_FourMom);
464 if(want_to_debug && kindofmass ==
"Comb"){
465 printf(
"NEW Comb jet mass = %f \n", jet_reco_Comb->
m());
466 printf(
"NEW CALO jet mass = %f \n", jet_reco_CALO->
m());
467 printf(
"NEW TA jet mass = %f \n", jet_reco_TA->
m());
472 if(want_to_debug && kindofmass==
"Calo"){ printf(
"CALO jet mass = %f \n", jet_reco->m()); }
473 if(want_to_debug && kindofmass==
"TA"){ printf(
"TA jet mass = %f \n", jet_reco_TA->
m()); }
475 if(kindofmass==
"TA"){m_JetTruthLabelingTool.modifyJet(*(jet_reco_TA)); };
478 if(!(ffjetsmearingtool.
getMatchedTruthJet(*jet_reco, jet_truth_matched).isSuccess())){
continue;}
480 if(lead_jet ==
true){reco_jet_mass_hist->Fill(jet_reco->m()*
MeVtoGeV); matched_truth_jet_mass_hist->Fill(jet_truth_matched.
m()*
MeVtoGeV); }
482 if(kindofmass==
"TA"){
484 Double_t aux_original_jet_mass = jet_reco->m()*
MeVtoGeV;
487 if(lead_jet ==
true && aux_original_jet_mass > 0){
488 smeared_reco_jet_mass_hist->Fill(jet_reco_TA->
m()*
MeVtoGeV); lead_jet=
false;
490 hist_jet_mass_scale_change_3D->Fill(jet_reco->pt()*
MeVtoGeV, aux_original_jet_mass, (jet_reco_TA->
m()*
MeVtoGeV));
494 Double_t aux_original_jet_mass = jet_reco->m()*
MeVtoGeV;
498 if(lead_jet ==
true && aux_original_jet_mass > 0){
499 smeared_reco_jet_mass_hist->Fill(jet_reco->m()*
MeVtoGeV); lead_jet=
false;
501 hist_jet_mass_scale_change_3D->Fill(jet_reco->pt()*
MeVtoGeV, aux_original_jet_mass, (jet_reco->m()*
MeVtoGeV));
503 hist_jet_mass_resolution_change_3D->Fill(jet_reco->pt()*
MeVtoGeV, aux_original_jet_mass, TMath::Abs((jet_reco->m()*
MeVtoGeV) - (aux_original_jet_mass))/ aux_original_jet_mass);
509 if(kindofmass==
"Comb"){ printf(
"Comb jet mass after smearing = %f \n", jet_reco->m()); }
510 if(kindofmass==
"Calo"){ printf(
"CALO jet mass after smearing = %f \n", jet_reco->m()); }
511 if(kindofmass==
"TA"){ printf(
"TA jet mass after smearing = %f \n", jet_reco_TA->
m()); }
517 if( ! ffjetsmearingtool.
correctedCopy(*jet_reco, jet_reco_copy) ){
518 if(want_to_debug){ printf(
"Failed to fill the corrected copy of the jet \n"); }
521 if(want_to_debug){ printf(
"Corrected Copy Comb jet mass after smearing = %f \n", jet_reco_copy->
m()); }
524 delete jet_reco_copy;
530 delete jets_shallowCopy.first;
531 delete jets_shallowCopy.second;
535 TString output_path =
"output/" +
sys.name() +
"_" +
output ;
537 TFile *jetmass_histograms =
new TFile(output_path,
"recreate");
540 reco_jet_mass_hist->Write();
541 matched_truth_jet_mass_hist->Write();
542 smeared_reco_jet_mass_hist->Write();
549 TH2F* hist_jet_mass_scale_change_2D =
new TH2F(
"hist_jet_mass_scale_change_2D",
"hist_jet_mass_scale_change_2D",numBinsPt,0,upperlimit2,numBinsMass,0,upperlimit1);
551 for (
int i=1;
i<=hist_jet_mass_scale_change_3D->GetNbinsX();
i++)
553 for (
int j=1; j<= hist_jet_mass_scale_change_3D->GetNbinsY(); j++)
555 TH1F* hold =
new TH1F(
"",
"",numBinsMass,0,upperlimit1);
556 for (
int k=1;
k<= hist_jet_mass_scale_change_3D->GetNbinsZ();
k++)
560 hist_jet_mass_scale_change_2D->
SetBinContent(
i,j,hold->GetMean()/hist_jet_mass_scale_change_3D->GetYaxis()->GetBinCenter(j));
567 TH2F* hist_jet_mass_resolution_change_2D =
new TH2F(
"hist_jet_mass_resolution_change_2D",
"hist_jet_mass_resolution_change_2D",numBinsPt,0,upperlimit2,numBinsMass,0,upperlimit1);
569 for (
int i=1;
i<=hist_jet_mass_resolution_change_3D->GetNbinsX();
i++)
571 for (
int j=1; j<= hist_jet_mass_resolution_change_3D->GetNbinsY(); j++)
573 TH1F* hold =
new TH1F(
"",
"",numBinsDiff,lowerlimit3,upperlimit3);
574 for (
int k=1;
k<= hist_jet_mass_resolution_change_3D->GetNbinsZ();
k++)
578 hist_jet_mass_resolution_change_2D->
SetBinContent(
i,j,hold->GetMean()/hist_jet_mass_resolution_change_3D->GetYaxis()->GetBinCenter(j));
587 TCanvas *
c1 =
new TCanvas(
"c1",
"c1",500,500);
590 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitleOffset(1.5);
591 hist_jet_mass_scale_change_2D->GetXaxis()->SetNdivisions(5);
592 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitleOffset(1.5);
593 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitleOffset(1.5);
594 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitle(
"Initial Reco Mass [GeV]");
595 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitle(
"Reco p_{T} [GeV]");
596 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitle(
"Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
598 hist_jet_mass_scale_change_2D->GetZaxis()->SetRangeUser(0.95,1.05);
601 hist_jet_mass_scale_change_2D->Draw(
"colz");
603 TString output_path_scale_debug =
"output/debug_plots/scale_variations/" +
sys.name() +
"_scaleDebug.pdf" ;
605 c1->Print(output_path_scale_debug);
607 TCanvas *
c2 =
new TCanvas(
"c2",
"c2",500,500);
609 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitleOffset(1.5);
610 hist_jet_mass_resolution_change_2D->GetXaxis()->SetNdivisions(5);
611 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitleOffset(1.5);
612 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitleOffset(1.5);
613 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitle(
"Initial Reco Mass [GeV]");
614 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitle(
"Reco p_{T} [GeV]");
615 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitle(
"Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
617 hist_jet_mass_resolution_change_2D->GetZaxis()->SetRangeUser(-0.1,0.1);
619 hist_jet_mass_resolution_change_2D->Draw(
"colz");
621 TString output_path_resolution_debug =
"output/debug_plots/resolution_variations/" +
sys.name() +
"_resolutionDebug.pdf" ;
623 c2->Print(output_path_resolution_debug);
625 delete hist_jet_mass_scale_change_2D;
627 delete hist_jet_mass_resolution_change_2D;
638 jetmass_histograms->Close();
640 delete jetmass_histograms;
643 delete reco_jet_mass_hist;
644 delete matched_truth_jet_mass_hist;
645 delete smeared_reco_jet_mass_hist;