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);
91 std::string item;
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 << truth_jetColl << reco_jetColl <<
output << string_debugtool << 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 if (!event) {
166 std::cout << "Failed to open file " << std::endl;
167 return 1;
168 }
169
171
172 const std::string name_FFJetSmearingTool = "FFJetSmearing_Example";
174 CHECK(ffjetsmearingtool.setProperty(
"MassDef", kindofmass));
175 CHECK(ffjetsmearingtool.setProperty(
"MCType", kindofmc));
176 CHECK(ffjetsmearingtool.setProperty(
"ConfigFile", string_configfile_path));
177 CHECK(ffjetsmearingtool.setProperty(
"OutputLevel", MSG::ERROR));
178 if (want_to_debug){
179 CHECK(ffjetsmearingtool.setProperty(
"OutputLevel", MSG::VERBOSE));
180 }
181 CHECK(ffjetsmearingtool.initialize());
182
184
185 const CP::SystematicSet& recommendedSysts = ffjetsmearingtool.recommendedSystematics();
186 std::cout << "\nRecommended systematics:\n" << std::endl;
187 for (
auto sysItr = recommendedSysts.
begin(); sysItr != recommendedSysts.
end(); ++sysItr){
188 std::cout << sysItr->name().c_str() << std::endl;
189 }
190
191 std::vector<CP::SystematicSet> sysList;
192
194
195 const std::string name_JetCalibTools = "JetCalibrationTool";
197
198 CHECK(m_JetCalibrationTool_handle.setProperty(
"JetCollection",
"AntiKt10UFOCSSKSoftDropBeta100Zcut10"));
199 CHECK(m_JetCalibrationTool_handle.setProperty(
"ConfigFile",
"JES_MC20PreRecommendation_R10_UFO_CSSK_SoftDrop_JMS_R21Insitu_26Nov2024.config"));
200
201 CHECK(m_JetCalibrationTool_handle.setProperty(
"CalibSequence",
"EtaJES_JMS"));
202 CHECK(m_JetCalibrationTool_handle.setProperty(
"IsData",
false));
203 CHECK(m_JetCalibrationTool_handle.setProperty(
"CalibArea",
"00-04-82"));
204 CHECK(m_JetCalibrationTool_handle.setProperty(
"DEVmode",
false));
205
206 CHECK(m_JetCalibrationTool_handle.initialize());
207
209
211
212
213
214
215
216
217 for (auto sys : recommendedSysts){
219 }
220
221 std::cout << "\n=============SYSTEMATICS CHECK NOW" << std::endl;
222 for (auto sys : sysList){
223
224 std::cout <<
"\nRunning over the systematic " <<
sys.name().c_str() << std::endl;
225 static constexpr float MeVtoGeV = 1.e-3;
226
227
228 if (ffjetsmearingtool.applySystematicVariation(sys) != StatusCode::SUCCESS){
229 std::cout << "Error, Cannot configure calibration tool for systematics" << std::endl;
230 }
231
232
233 TH1F* reco_jet_mass_hist =
new TH1F(
"reco_jet_mass_hist",
"reco_jet_mass_hist", 300, 0, 600);
234 TH1F* matched_truth_jet_mass_hist =
new TH1F(
"matched_truth_jet_mass_hist",
"matched_truth_jet_mass_hist", 300, 0, 600);
235 TH1F* smeared_reco_jet_mass_hist =
new TH1F(
"smeared_reco_jet_mass_hist",
"smeared_reco_jet_mass_hist", 300, 0, 600);
236
237 TH1F* reco_jet_pt_hist =
new TH1F(
"reco_jet_pt_hist",
"reco_jet_pt_hist", 800, 0, 4000);
238 TH1F* matched_truth_jet_pt_hist =
new TH1F(
"matched_truth_jet_pt_hist",
"matched_truth_jet_pt_hist", 800, 0, 4000);
239 TH1F* smeared_reco_jet_pt_hist =
new TH1F(
"smeared_reco_jet_pt_hist",
"smeared_reco_jet_pt_hist", 800, 0, 4000);
240
241 TH1F* reco_jet_rapidity_hist =
new TH1F(
"reco_jet_rapidity_hist",
"reco_jet_rapidity_hist", 50, -2.5, 2.5);
242 TH1F* matched_truth_jet_rapidity_hist =
new TH1F(
"matched_truth_jet_rapidity_hist",
"matched_truth_jet_rapidity_hist", 50, -2.5, 2.5);
243 TH1F* smeared_reco_jet_rapidity_hist =
new TH1F(
"smeared_reco_jet_rapidity_hist",
"smeared_reco_jet_rapidity_hist", 50, -2.5, 2.5);
244
245 reco_jet_mass_hist->Sumw2();
246 matched_truth_jet_mass_hist->Sumw2();
247 smeared_reco_jet_mass_hist->Sumw2();
248
249 reco_jet_pt_hist->Sumw2();
250 matched_truth_jet_pt_hist->Sumw2();
251 smeared_reco_jet_pt_hist->Sumw2();
252
253 reco_jet_rapidity_hist->Sumw2();
254 matched_truth_jet_rapidity_hist->Sumw2();
255 smeared_reco_jet_rapidity_hist->Sumw2();
256
257 int upperlimit1 = 600;
258 int upperlimit2 = 1000;
259
260 int numBinsMass = 120;
261 int numBinsPt = 100;
262
263 std::unique_ptr<TH3F> hist_jet_mass_scale_change_3D;
264 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);
265
266 float lowerlimit3 = -0.5;
267 float upperlimit3 = 0.5;
268 int numBinsDiff = 100;
269
270 std::unique_ptr<TH3F> hist_jet_mass_resolution_change_3D;
271 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);
272
274
275 Long64_t
nevents =
event->getEntries();
276
277 if (eventsMax < nevents){
279 }
280
281 for (Long64_t ievent = 0; ievent <
nevents; ++ievent){
282
283
284 if (
event->getEntry( ievent ) < 0 ){
285 std::cout << "Failed to load entry " << ievent << std::endl;
286 return 1;
287 }
288
289
290 if (ievent % 1000==0){
291 std::cout <<
"Event " << ievent <<
" of " <<
nevents << std::endl;
292 }
293
294
295 if (want_to_debug){
298
299 std::cout <<
"===>>> start processing event " << ei->
eventNumber() <<
", run " << ei->
runNumber() <<
" - Events processed so far: " << ievent << std::endl;
300
301
303 CHECK(
event->retrieve(jets_truth, truth_jetColl) );
304 std::cout <<
"Number of truth jets: " << jets_truth->
size() << std::endl;
305
306
307 for (
const xAOD::Jet* jet_truth : *jets_truth){
308
309 std::cout <<
"Truth Jet: pt = " << jet_truth->pt()*
MeVtoGeV <<
", mass = " << jet_truth->m()*
MeVtoGeV <<
", eta = " << jet_truth->eta() << std::endl;
310 }
311
312
314 CHECK(
event->retrieve(jets_reco, reco_jetColl) );
315 std::cout <<
"Number of reco jets: " << jets_reco->
size() << std::endl;
316
317
318 for (
const xAOD::Jet* jet_reco : *jets_reco){
319
320 std::cout <<
"Reco Jet: pt = " << jet_reco->pt()*
MeVtoGeV <<
", mass = " << jet_reco->m()*
MeVtoGeV <<
", eta = " << jet_reco->eta() << std::endl;
321 }
322 }
323
325 jet_truth_matched.makePrivateStore();
326
327
329 CHECK(
event->retrieve( jets, reco_jetColl ) );
330
331
333
334 if (want_to_debug){
335 std::cout << "Start the loop over the jets" << std::endl;
336 }
337
338 bool lead_jet = true;
339
340
341 if ( m_JetCalibrationTool_handle.applyCalibration( *(jets_shallowCopy.first) ) != StatusCode::SUCCESS ){
342 std::cout << "JetCalibration tool reported a EL::StatusCode::FAILURE" << std::endl;
343 return 1;
344 }
345
347
348 for (
xAOD::Jet* jet_reco : *(jets_shallowCopy.first) ){
349
350 double jetpt = jet_reco->pt() *
MeVtoGeV;
351 double jetm = jet_reco->m() *
MeVtoGeV;
352 double jetrapidity = jet_reco->rapidity();
353
354 if (jetpt < 200 or jetpt > 3000) continue;
355
356 if (jetm > 600) continue;
357
358 if (abs(jetrapidity) > 2) continue;
359
360 if ( ffjetsmearingtool.getMatchedTruthJet(*jet_reco, jet_truth_matched) != StatusCode::SUCCESS ){
361 continue;
362 }
363
364 double aux_original_jet_mass = jet_reco->m() *
MeVtoGeV;
365
366 if (lead_jet == true && aux_original_jet_mass > 0){
367
368 reco_jet_mass_hist->Fill(jet_reco->m() * MeVtoGeV);
369 reco_jet_pt_hist->Fill(jet_reco->pt() * MeVtoGeV);
370 reco_jet_rapidity_hist->Fill(jet_reco->rapidity());
371
372 matched_truth_jet_mass_hist->Fill(jet_truth_matched.
m() * MeVtoGeV);
373 matched_truth_jet_pt_hist->Fill(jet_truth_matched.
pt() * MeVtoGeV);
374 matched_truth_jet_rapidity_hist->Fill(jet_truth_matched.
rapidity());
375
376
378 std::cout << "FFJetSmearingTool reported a EL::StatusCode::FAILURE" << std::endl;
379 return 1;
380 }
381
382 smeared_reco_jet_mass_hist->Fill(jet_reco->m() * MeVtoGeV);
383 smeared_reco_jet_pt_hist->Fill(jet_reco->pt() * MeVtoGeV);
384 smeared_reco_jet_rapidity_hist->Fill(jet_reco->rapidity());
385
386 hist_jet_mass_scale_change_3D->Fill(jet_reco->pt() * MeVtoGeV, aux_original_jet_mass, jet_reco->m() * MeVtoGeV);
387 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);
388
389 lead_jet = false;
390
391 }
392
393 }
394
395 }
396
397 TString output_path =
"output/" +
sys.name() +
output ;
398 TFile *jetmass_histograms = new TFile(output_path,"recreate");
399
400 reco_jet_mass_hist->Write();
401 matched_truth_jet_mass_hist->Write();
402 smeared_reco_jet_mass_hist->Write();
403
404 reco_jet_pt_hist->Write();
405 matched_truth_jet_pt_hist->Write();
406 smeared_reco_jet_pt_hist->Write();
407
408 reco_jet_rapidity_hist->Write();
409 matched_truth_jet_rapidity_hist->Write();
410 smeared_reco_jet_rapidity_hist->Write();
411
413
414 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);
415 for (
int i=1;
i<=hist_jet_mass_scale_change_3D->GetNbinsX();
i++){
416 for (
int j=1;
j<= hist_jet_mass_scale_change_3D->GetNbinsY();
j++){
417 TH1F* hold =
new TH1F(
"",
"",numBinsMass,0,upperlimit1);
418 for (
int k=1;
k<= hist_jet_mass_scale_change_3D->GetNbinsZ();
k++){
419 hold->SetBinContent(k,hist_jet_mass_scale_change_3D->GetBinContent(i,j,k));
420 }
421 hist_jet_mass_scale_change_2D->SetBinContent(i,j,hold->GetMean()/hist_jet_mass_scale_change_3D->GetYaxis()->GetBinCenter(j));
422 delete hold;
423 }
424 }
425
426 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);
427 for (
int i=1;
i<=hist_jet_mass_resolution_change_3D->GetNbinsX();
i++){
428 for (
int j=1;
j<= hist_jet_mass_resolution_change_3D->GetNbinsY();
j++){
429 TH1F* hold =
new TH1F(
"",
"",numBinsDiff,lowerlimit3,upperlimit3);
430 for (
int k=1;
k<= hist_jet_mass_resolution_change_3D->GetNbinsZ();
k++){
431 hold->SetBinContent(k,hist_jet_mass_resolution_change_3D->GetBinContent(i,j,k));
432 }
433 hist_jet_mass_resolution_change_2D->SetBinContent(i,j,hold->GetMean()/hist_jet_mass_resolution_change_3D->GetYaxis()->GetBinCenter(j));
434 delete hold;
435 }
436 }
437
439
442
443 TCanvas *
c1 =
new TCanvas(
"c1",
"c1",w,
h);
444 c1->SetWindowSize(w+(w-
c1->GetWw()),
h+(
h-
c1->GetWh()));
446
447 TPad *canvas_1 = new TPad("canvas_1", "canvas_1",0.0,0.0,1.0,1.0);
448 canvas_1->SetRightMargin(0.10);
449 canvas_1->SetFillStyle(4000);
450 canvas_1->Draw();
451 canvas_1->cd();
452
453 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitleOffset(1.5);
454 hist_jet_mass_scale_change_2D->GetXaxis()->SetNdivisions(5);
455 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitleOffset(1.5);
456 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitleOffset(1.5);
457 hist_jet_mass_scale_change_2D->GetYaxis()->SetTitle("Initial Reco Mass [GeV]");
458 hist_jet_mass_scale_change_2D->GetXaxis()->SetTitle("Reco p_{T} [GeV]");
459 hist_jet_mass_scale_change_2D->GetZaxis()->SetTitle("Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
460 hist_jet_mass_scale_change_2D->GetZaxis()->SetRangeUser(0.95,1.05);
461 hist_jet_mass_scale_change_2D->GetZaxis()->SetLabelSize(0.035);
462 hist_jet_mass_scale_change_2D->Draw("colz");
463 gPad->RedrawAxis();
464
465 TString output_path_scale_debug =
"output/debug_plots/scale_variations/" +
sys.name() +
"_scaleDebug.pdf" ;
466 c1->Print(output_path_scale_debug);
467
468 delete hist_jet_mass_scale_change_2D;
470
471 TCanvas *
c2 =
new TCanvas(
"c2",
"c2",w,
h);
472 c2->SetWindowSize(w+(w-
c2->GetWw()),
h+(
h-
c2->GetWh()));
474
475 TPad *canvas_2 = new TPad("canvas_2", "canvas_2",0.0,0.0,1.0,1.0);
476 canvas_2->SetRightMargin(0.10);
477 canvas_2->SetFillStyle(4000);
478 canvas_2->Draw();
479 canvas_2->cd();
480
481 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitleOffset(1.5);
482 hist_jet_mass_resolution_change_2D->GetXaxis()->SetNdivisions(5);
483 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitleOffset(1.5);
484 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitleOffset(1.5);
485 hist_jet_mass_resolution_change_2D->GetYaxis()->SetTitle("Initial Reco Mass [GeV]");
486 hist_jet_mass_resolution_change_2D->GetXaxis()->SetTitle("Reco p_{T} [GeV]");
487 hist_jet_mass_resolution_change_2D->GetZaxis()->SetTitle("Average Mass smearing (Initial_reco_mass/smeared_reco_mass)");
488 hist_jet_mass_resolution_change_2D->GetZaxis()->SetRangeUser(-0.1,0.1);
489 hist_jet_mass_resolution_change_2D->GetZaxis()->SetLabelSize(0.035);
490 hist_jet_mass_resolution_change_2D->Draw("colz");
491 gPad->RedrawAxis();
492
493 TString output_path_resolution_debug =
"output/debug_plots/resolution_variations/" +
sys.name() +
"_resolutionDebug.pdf" ;
494 c2->Print(output_path_resolution_debug);
495
496 delete hist_jet_mass_resolution_change_2D;
498
499 jetmass_histograms->Close();
500 delete jetmass_histograms;
501
502 delete reco_jet_mass_hist;
503 delete matched_truth_jet_mass_hist;
504 delete smeared_reco_jet_mass_hist;
505
506 delete reco_jet_pt_hist;
507 delete matched_truth_jet_pt_hist;
508 delete smeared_reco_jet_pt_hist;
509
510 delete reco_jet_rapidity_hist;
511 delete matched_truth_jet_rapidity_hist;
512 delete smeared_reco_jet_rapidity_hist;
513
514 }
515
516 return 0;
517
518}
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".