71 {
72
74
77 std::string truth_jetColl = "";
78 std::string reco_jetColl = "";
79 std::string string_kindofmass = "";
80 std::string string_kindofmc = "";
82 std::string string_configfile_path = "";
83 std::string string_debugtool = "";
84 int eventsMax = 30000;
85
87 for (
int i=1;
i<
argc;
i++){
88
89 std::string
opt(argv[i]); std::vector< std::string >
v;
90 std::istringstream iss(opt);
92 char delim = '=';
93
94 while (std::getline(iss, item, delim)){
96 }
97
98 if (
opt.find(
"--help") != std::string::npos ){
100 return 0;
101 }
102
103 if (
opt.find(
"--sample=") != std::string::npos )
sample =
v[1];
104 if (
opt.find(
"--truth_jetColl=") != std::string::npos ) truth_jetColl =
v[1];
105 if (
opt.find(
"--reco_jetColl=") != std::string::npos ) reco_jetColl =
v[1];
106 if (
opt.find(
"--MassDef=") != std::string::npos ) string_kindofmass =
v[1];
107 if (
opt.find(
"--MCType=") != std::string::npos ) string_kindofmc =
v[1];
108 if (
opt.find(
"--ConfigFile=") != std::string::npos ) string_configfile_path =
v[1];
109 if (
opt.find(
"--DebugTool=") != std::string::npos ) string_debugtool =
v[1];
110 if (
opt.find(
"--output=") != std::string::npos )
output =
v[1];
111 if (
opt.find(
"--eventsMax=") != std::string::npos ) eventsMax = std::atoi(v[1].
data());
112
113 }
114
115 if (sample==""){
116 std::cout << "No input xAOD file specified, exiting" << std::endl;
117 return 1;
118 }
119 if (truth_jetColl==""){
120 std::cout << "No truth jet collection specified, exiting" << std::endl;
121 return 1;
122 }
123 if (reco_jetColl==""){
124 std::cout << "No truth jet collection specified, exiting" << std::endl;
125 return 1;
126 }
127 if (string_kindofmass==""){
128 std::cout << "No kind of jet mass specified, exiting" << std::endl;
129 return 1;
130 }
131 if (string_kindofmc==""){
132 std::cout << "No kind of MC specified, exiting" << std::endl;
133 return 1;
134 }
135 if (string_configfile_path==""){
136 std::cout << "No ConfigFile specified, exiting" << std::endl;
137 return 1;
138 }
139 if (string_debugtool==""){
140 std::cout << "No debugtool specified, exiting" << std::endl;
141 return 1;
142 }
143 if (output==""){
144 std::cout << "Output not specified, exiting" << std::endl;
145 return 1;
146 }
147 std::cout <<
sample.c_str() << truth_jetColl.c_str() << reco_jetColl.c_str() <<
output.c_str() << string_debugtool.c_str() << std::endl;
148
149 TString kindofmass = string_kindofmass;
150 TString kindofmc = string_kindofmc;
151
152 bool want_to_debug = false;
153 if (string_debugtool == "true"){
154 want_to_debug = true;
155 }else if (string_debugtool == "false"){
156 want_to_debug = false;
157 }
158
160
161 std::unique_ptr< TFile >
ifile( TFile::Open(
sample.c_str(),
"READ" ) );
162
163
165
167
169
170 const std::string name_FFJetSmearingTool = "FFJetSmearing_Example";
172 CHECK(ffjetsmearingtool.setProperty(
"MassDef", kindofmass));
173 CHECK(ffjetsmearingtool.setProperty(
"MCType", kindofmc));
174 CHECK(ffjetsmearingtool.setProperty(
"ConfigFile", string_configfile_path));
175 CHECK(ffjetsmearingtool.setProperty(
"OutputLevel", MSG::ERROR));
176 if (want_to_debug){
177 CHECK(ffjetsmearingtool.setProperty(
"OutputLevel", MSG::VERBOSE));
178 }
179 CHECK(ffjetsmearingtool.initialize());
180
182
183 const CP::SystematicSet& recommendedSysts = ffjetsmearingtool.recommendedSystematics();
184 std::cout << "\nRecommended systematics:\n" << std::endl;
185 for (
auto sysItr = recommendedSysts.
begin(); sysItr != recommendedSysts.
end(); ++sysItr){
186 std::cout << sysItr->name().c_str() << std::endl;
187 }
188
189 std::vector<CP::SystematicSet> sysList;
190
192
193
194
195
196
197
198
199
200
202
203 const std::string name_JetCalibTools = "JetCalibrationTool";
205
206 CHECK(m_JetCalibrationTool_handle.setProperty(
"JetCollection",
"AntiKt10UFOCSSKSoftDropBeta100Zcut10"));
207 CHECK(m_JetCalibrationTool_handle.setProperty(
"ConfigFile",
"JES_MC20PreRecommendation_R10_UFO_CSSK_SoftDrop_JMS_R21Insitu_26Nov2024.config"));
208
209 CHECK(m_JetCalibrationTool_handle.setProperty(
"CalibSequence",
"EtaJES_JMS"));
210 CHECK(m_JetCalibrationTool_handle.setProperty(
"IsData",
false));
211 CHECK(m_JetCalibrationTool_handle.setProperty(
"CalibArea",
"00-04-82"));
212 CHECK(m_JetCalibrationTool_handle.setProperty(
"DEVmode",
false));
213
214 CHECK(m_JetCalibrationTool_handle.initialize());
215
217
219
220
221
222
223
224
225 for (auto sys : recommendedSysts){
227 }
228
229 std::cout << "\n=============SYSTEMATICS CHECK NOW" << std::endl;
230 for (auto sys : sysList){
231
232 std::cout <<
"\nRunning over the systematic " <<
sys.name().c_str() << std::endl;
233 static constexpr float MeVtoGeV = 1.e-3;
234
235
236 if (ffjetsmearingtool.applySystematicVariation(sys) != StatusCode::SUCCESS){
237 std::cout << "Error, Cannot configure calibration tool for systematics" << std::endl;
238 }
239
240
241 TH1F* reco_jet_mass_hist =
new TH1F(
"reco_jet_mass_hist",
"reco_jet_mass_hist", 300, 0, 600);
242 TH1F* matched_truth_jet_mass_hist =
new TH1F(
"matched_truth_jet_mass_hist",
"matched_truth_jet_mass_hist", 300, 0, 600);
243 TH1F* smeared_reco_jet_mass_hist =
new TH1F(
"smeared_reco_jet_mass_hist",
"smeared_reco_jet_mass_hist", 300, 0, 600);
244
245 TH1F* reco_jet_pt_hist =
new TH1F(
"reco_jet_pt_hist",
"reco_jet_pt_hist", 800, 0, 4000);
246 TH1F* matched_truth_jet_pt_hist =
new TH1F(
"matched_truth_jet_pt_hist",
"matched_truth_jet_pt_hist", 800, 0, 4000);
247 TH1F* smeared_reco_jet_pt_hist =
new TH1F(
"smeared_reco_jet_pt_hist",
"smeared_reco_jet_pt_hist", 800, 0, 4000);
248
249 TH1F* reco_jet_rapidity_hist =
new TH1F(
"reco_jet_rapidity_hist",
"reco_jet_rapidity_hist", 50, -2.5, 2.5);
250 TH1F* matched_truth_jet_rapidity_hist =
new TH1F(
"matched_truth_jet_rapidity_hist",
"matched_truth_jet_rapidity_hist", 50, -2.5, 2.5);
251 TH1F* smeared_reco_jet_rapidity_hist =
new TH1F(
"smeared_reco_jet_rapidity_hist",
"smeared_reco_jet_rapidity_hist", 50, -2.5, 2.5);
252
253 reco_jet_mass_hist->Sumw2();
254 matched_truth_jet_mass_hist->Sumw2();
255 smeared_reco_jet_mass_hist->Sumw2();
256
257 reco_jet_pt_hist->Sumw2();
258 matched_truth_jet_pt_hist->Sumw2();
259 smeared_reco_jet_pt_hist->Sumw2();
260
261 reco_jet_rapidity_hist->Sumw2();
262 matched_truth_jet_rapidity_hist->Sumw2();
263 smeared_reco_jet_rapidity_hist->Sumw2();
264
265 int upperlimit1 = 600;
266 int upperlimit2 = 1000;
267
268 int numBinsMass = 120;
269 int numBinsPt = 100;
270
271 std::unique_ptr<TH3F> hist_jet_mass_scale_change_3D;
272 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);
273
274 float lowerlimit3 = -0.5;
275 float upperlimit3 = 0.5;
276 int numBinsDiff = 100;
277
278 std::unique_ptr<TH3F> hist_jet_mass_resolution_change_3D;
279 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);
280
282
283 Long64_t
nevents =
event.getEntries();
284
285 if (eventsMax < nevents){
287 }
288
289 for (Long64_t ievent = 0; ievent <
nevents; ++ievent){
290
291
292 if (
event.getEntry( ievent ) < 0 ){
293 std::cout << "Failed to load entry " << ievent << std::endl;
294 return 1;
295 }
296
297
298 if (ievent % 1000==0){
299 std::cout <<
"Event " << ievent <<
" of " <<
nevents << std::endl;
300 }
301
302
303 if (want_to_debug){
306
307 std::cout <<
"===>>> start processing event " << ei->
eventNumber() <<
", run " << ei->
runNumber() <<
" - Events processed so far: " << ievent << std::endl;
308
309
311 CHECK(
event.retrieve(jets_truth, truth_jetColl) );
312 std::cout <<
"Number of truth jets: " << jets_truth->
size() << std::endl;
313
314
315 for (
const xAOD::Jet* jet_truth : *jets_truth){
316
317 std::cout <<
"Truth Jet: pt = " << jet_truth->pt()*
MeVtoGeV <<
", mass = " << jet_truth->m()*
MeVtoGeV <<
", eta = " << jet_truth->eta() << std::endl;
318 }
319
320
322 CHECK(
event.retrieve(jets_reco, reco_jetColl) );
323 std::cout <<
"Number of reco jets: " << jets_reco->
size() << std::endl;
324
325
326 for (
const xAOD::Jet* jet_reco : *jets_reco){
327
328 std::cout <<
"Reco Jet: pt = " << jet_reco->pt()*
MeVtoGeV <<
", mass = " << jet_reco->m()*
MeVtoGeV <<
", eta = " << jet_reco->eta() << std::endl;
329 }
330 }
331
334
335
337 CHECK(
event.retrieve( jets, reco_jetColl ) );
338
339
341
343
344
345
346 if (want_to_debug){
347 std::cout << "Start the loop over the jets" << std::endl;
348 }
349
350 bool lead_jet = true;
351
352
353 if ( m_JetCalibrationTool_handle.applyCalibration( *(jets_shallowCopy.first) ) != StatusCode::SUCCESS ){
354 std::cout << "JetCalibration tool reported a EL::StatusCode::FAILURE" << std::endl;
355 return 1;
356 }
357
359
360 for (
xAOD::Jet* jet_reco : *(jets_shallowCopy.first) ){
361
362 double jetpt = jet_reco->pt() *
MeVtoGeV;
363 double jetm = jet_reco->m() *
MeVtoGeV;
364 double jetrapidity = jet_reco->rapidity();
365
366 if (jetpt < 200 or jetpt > 3000) continue;
367
368 if (jetm > 600) continue;
369
370 if (abs(jetrapidity) > 2) continue;
371
372 if ( ffjetsmearingtool.getMatchedTruthJet(*jet_reco, jet_truth_matched) != StatusCode::SUCCESS ){
373 continue;
374 }
375
376 double aux_original_jet_mass = jet_reco->m() *
MeVtoGeV;
377
378 if (lead_jet == true && aux_original_jet_mass > 0){
379
380 reco_jet_mass_hist->Fill(jet_reco->m() * MeVtoGeV);
381 reco_jet_pt_hist->Fill(jet_reco->pt() * MeVtoGeV);
382 reco_jet_rapidity_hist->Fill(jet_reco->rapidity());
383
384 matched_truth_jet_mass_hist->Fill(jet_truth_matched.
m() * MeVtoGeV);
385 matched_truth_jet_pt_hist->Fill(jet_truth_matched.
pt() * MeVtoGeV);
386 matched_truth_jet_rapidity_hist->Fill(jet_truth_matched.
rapidity());
387
388
390 std::cout << "FFJetSmearingTool reported a EL::StatusCode::FAILURE" << std::endl;
391 return 1;
392 }
393
394 smeared_reco_jet_mass_hist->Fill(jet_reco->m() * MeVtoGeV);
395 smeared_reco_jet_pt_hist->Fill(jet_reco->pt() * MeVtoGeV);
396 smeared_reco_jet_rapidity_hist->Fill(jet_reco->rapidity());
397
398 hist_jet_mass_scale_change_3D->Fill(jet_reco->pt() * MeVtoGeV, aux_original_jet_mass, jet_reco->m() * MeVtoGeV);
399 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);
400
401 lead_jet = false;
402
403 }
404
405 }
406
407 delete jets_shallowCopy.first;
408 delete jets_shallowCopy.second;
409
410 }
411
412 TString output_path =
"output/" +
sys.name() +
output ;
413 TFile *jetmass_histograms = new TFile(output_path,"recreate");
414
415 reco_jet_mass_hist->Write();
416 matched_truth_jet_mass_hist->Write();
417 smeared_reco_jet_mass_hist->Write();
418
419 reco_jet_pt_hist->Write();
420 matched_truth_jet_pt_hist->Write();
421 smeared_reco_jet_pt_hist->Write();
422
423 reco_jet_rapidity_hist->Write();
424 matched_truth_jet_rapidity_hist->Write();
425 smeared_reco_jet_rapidity_hist->Write();
426
428
429 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);
430 for (
int i=1;
i<=hist_jet_mass_scale_change_3D->GetNbinsX();
i++){
431 for (int j=1; j<= hist_jet_mass_scale_change_3D->GetNbinsY(); j++){
432 TH1F* hold =
new TH1F(
"",
"",numBinsMass,0,upperlimit1);
433 for (
int k=1;
k<= hist_jet_mass_scale_change_3D->GetNbinsZ();
k++){
434 hold->SetBinContent(k,hist_jet_mass_scale_change_3D->GetBinContent(i,j,k));
435 }
436 hist_jet_mass_scale_change_2D->SetBinContent(i,j,hold->GetMean()/hist_jet_mass_scale_change_3D->GetYaxis()->GetBinCenter(j));
437 delete hold;
438 }
439 }
440
441 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);
442 for (
int i=1;
i<=hist_jet_mass_resolution_change_3D->GetNbinsX();
i++){
443 for (int j=1; j<= hist_jet_mass_resolution_change_3D->GetNbinsY(); j++){
444 TH1F* hold =
new TH1F(
"",
"",numBinsDiff,lowerlimit3,upperlimit3);
445 for (
int k=1;
k<= hist_jet_mass_resolution_change_3D->GetNbinsZ();
k++){
446 hold->SetBinContent(k,hist_jet_mass_resolution_change_3D->GetBinContent(i,j,k));
447 }
448 hist_jet_mass_resolution_change_2D->SetBinContent(i,j,hold->GetMean()/hist_jet_mass_resolution_change_3D->GetYaxis()->GetBinCenter(j));
449 delete hold;
450 }
451 }
452
454
457
458 TCanvas *
c1 =
new TCanvas(
"c1",
"c1",w,
h);
459 c1->SetWindowSize(w+(w-
c1->GetWw()),
h+(
h-
c1->GetWh()));
461
462 TPad *canvas_1 = new TPad("canvas_1", "canvas_1",0.0,0.0,1.0,1.0);
463 canvas_1->SetRightMargin(0.10);
464 canvas_1->SetFillStyle(4000);
465 canvas_1->Draw();
466 canvas_1->cd();
467
468 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitleOffset(1.5);
469 hist_jet_mass_scale_change_2D->GetXaxis()->SetNdivisions(5);
470 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitleOffset(1.5);
471 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitleOffset(1.5);
472 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitle("Initial Reco Mass [GeV]");
473 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitle("Reco p_{T} [GeV]");
474 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitle("Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
475 hist_jet_mass_scale_change_2D->GetZaxis()->SetRangeUser(0.95,1.05);
476 hist_jet_mass_scale_change_2D->GetZaxis()->SetLabelSize(0.035);
477 hist_jet_mass_scale_change_2D->Draw("colz");
478 gPad->RedrawAxis();
479
480 TString output_path_scale_debug =
"output/debug_plots/scale_variations/" +
sys.name() +
"_scaleDebug.pdf" ;
481 c1->Print(output_path_scale_debug);
482
483 delete hist_jet_mass_scale_change_2D;
485
486 TCanvas *
c2 =
new TCanvas(
"c2",
"c2",w,
h);
487 c2->SetWindowSize(w+(w-
c2->GetWw()),
h+(
h-
c2->GetWh()));
489
490 TPad *canvas_2 = new TPad("canvas_2", "canvas_2",0.0,0.0,1.0,1.0);
491 canvas_2->SetRightMargin(0.10);
492 canvas_2->SetFillStyle(4000);
493 canvas_2->Draw();
494 canvas_2->cd();
495
496 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitleOffset(1.5);
497 hist_jet_mass_resolution_change_2D->GetXaxis()->SetNdivisions(5);
498 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitleOffset(1.5);
499 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitleOffset(1.5);
500 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitle("Initial Reco Mass [GeV]");
501 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitle("Reco p_{T} [GeV]");
502 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitle("Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
503 hist_jet_mass_resolution_change_2D->GetZaxis()->SetRangeUser(-0.1,0.1);
504 hist_jet_mass_resolution_change_2D->GetZaxis()->SetLabelSize(0.035);
505 hist_jet_mass_resolution_change_2D->Draw("colz");
506 gPad->RedrawAxis();
507
508 TString output_path_resolution_debug =
"output/debug_plots/resolution_variations/" +
sys.name() +
"_resolutionDebug.pdf" ;
509 c2->Print(output_path_resolution_debug);
510
511 delete hist_jet_mass_resolution_change_2D;
513
514 jetmass_histograms->Close();
515 delete jetmass_histograms;
516
517 delete reco_jet_mass_hist;
518 delete matched_truth_jet_mass_hist;
519 delete smeared_reco_jet_mass_hist;
520
521 delete reco_jet_pt_hist;
522 delete matched_truth_jet_pt_hist;
523 delete smeared_reco_jet_pt_hist;
524
525 delete reco_jet_rapidity_hist;
526 delete matched_truth_jet_rapidity_hist;
527 delete smeared_reco_jet_rapidity_hist;
528
529 }
530
531 return 0;
532
533}
char data[hepevt_bytes_allocation_ATLAS]
Header file for AthHistogramAlgorithm.
@ Error
Some error happened during the object correction.
Class to wrap a set of SystematicVariations.
const_iterator end() const
description: const iterator to the end of the set
const_iterator begin() const
description: const iterator to the beginning of the set
size_type size() const noexcept
Returns the number of elements in the collection.
void makePrivateStore()
Create a new (empty) private store for this object.
uint32_t runNumber() const
The current event's run number.
uint64_t eventNumber() const
The current event's event number.
virtual double pt() const
The transverse momentum ( ) of the particle.
virtual double rapidity() const
The true rapidity (y) of the particle.
virtual double m() const
The invariant mass of the particle.
Tool for accessing xAOD files outside of Athena.
@ kClassAccess
Access auxiliary data using the aux containers.
TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Jet_v1 Jet
Definition of the current "jet version".
EventInfo_v1 EventInfo
Definition of the latest event info version.
std::pair< std::unique_ptr< T >, std::unique_ptr< ShallowAuxContainer > > shallowCopyContainer(const T &cont, const EventContext &ctx)
Function making a shallow copy of a constant container.
JetContainer_v1 JetContainer
Definition of the current "jet container version".