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