Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationJet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationJet.cc
1 #include "QAG4SimulationJet.h"
2 
3 #include <qautils/QAHistManagerDef.h>
4 
5 #include <g4eval/JetEvalStack.h>
6 #include <g4eval/JetRecoEval.h>
7 #include <g4eval/JetTruthEval.h>
8 
9 #include <jetbase/Jet.h>
10 #include <jetbase/JetContainer.h>
11 
12 #include <g4main/PHG4HitDefs.h>
13 #include <g4main/PHG4Shower.h>
14 
15 #include <fun4all/Fun4AllBase.h>
16 #include <fun4all/Fun4AllHistoManager.h>
17 #include <fun4all/Fun4AllReturnCodes.h>
18 #include <fun4all/SubsysReco.h> // for SubsysReco
19 
20 #include <phool/getClass.h>
21 
22 #include <TAxis.h>
23 #include <TH1.h>
24 #include <TH2.h>
25 #include <TNamed.h>
26 #include <TString.h> // for operator+, TString, Form
27 
28 #include <cassert>
29 #include <cmath>
30 #include <cstdlib>
31 #include <iostream>
32 #include <set>
33 
34 QAG4SimulationJet::QAG4SimulationJet(const std::string& truth_jet,
35  enu_flags flags)
36  : SubsysReco("QAG4SimulationJet_" + truth_jet)
37  , _jetevalstacks()
38  , _truth_jet(truth_jet)
39  , _reco_jets()
40  , _flags(flags)
41  , eta_range(-1, 1)
42  , _jet_match_dEta(.1)
43  , _jet_match_dPhi(.1)
44  , _jet_match_dE_Ratio(.5)
45 {
46 }
47 
48 int QAG4SimulationJet::InitRun(PHCompositeNode* topNode)
49 {
50  if (flag(kProcessTruthMatching) || flag(kProcessRecoSpectrum))
51  {
52  for (const auto& reco_jet : _reco_jets)
53  {
54  jetevalstacks_map::iterator it_jetevalstack = _jetevalstacks.find(
55  reco_jet);
56 
57  if (it_jetevalstack == _jetevalstacks.end())
58  {
59  _jetevalstacks[reco_jet] = std::shared_ptr<JetEvalStack>(new JetEvalStack(topNode, reco_jet, _truth_jet));
60  assert(_jetevalstacks[reco_jet]);
61  _jetevalstacks[reco_jet]->set_strict(true);
62  _jetevalstacks[reco_jet]->set_verbosity(Verbosity() + 1);
63  }
64  else
65  {
66  assert(it_jetevalstack->second);
67  }
68  }
69  }
70 
71  if (flag(kProcessTruthSpectrum))
72  {
73  if (not _jettrutheval)
74  {
75  _jettrutheval = std::shared_ptr<JetTruthEval>(new JetTruthEval(topNode, _truth_jet));
76  }
77 
78  assert(_jettrutheval);
79  _jettrutheval->set_strict(true);
80  _jettrutheval->set_verbosity(Verbosity() + 1);
81  }
82 
83  return Fun4AllReturnCodes::EVENT_OK;
84 }
85 
86 int QAG4SimulationJet::Init(PHCompositeNode* topNode)
87 {
88  Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
89  assert(hm);
90 
91  if (flag(kProcessTruthSpectrum))
92  {
93  if (Verbosity() >= 1)
94  {
95  std::cout << "QAG4SimulationJet::Init - Process TruthSpectrum " << _truth_jet
96  << std::endl;
97  }
98  Init_Spectrum(topNode, _truth_jet);
99  }
100 
101  if (flag(kProcessRecoSpectrum))
102  {
103  for (const auto& reco_jet : _reco_jets)
104  {
105  if (Verbosity() >= 1)
106  {
107  std::cout << "QAG4SimulationJet::Init - Process Reco jet spectrum "
108  << reco_jet << std::endl;
109  }
110  Init_Spectrum(topNode, reco_jet);
111  }
112  }
113 
114  if (flag(kProcessTruthMatching))
115  {
116  for (const auto& reco_jet : _reco_jets)
117  {
118  if (Verbosity() >= 1)
119  {
120  std::cout << "QAG4SimulationJet::Init - Process Reco jet spectrum "
121  << reco_jet << std::endl;
122  }
123  Init_TruthMatching(topNode, reco_jet);
124  }
125  }
126 
127  return Fun4AllReturnCodes::EVENT_OK;
128 }
129 
130 int QAG4SimulationJet::process_event(PHCompositeNode* topNode)
131 {
132  if (Verbosity() > 2)
133  {
134  std::cout << "QAG4SimulationJet::process_event() entered" << std::endl;
135  }
136 
137  for (jetevalstacks_map::iterator it_jetevalstack = _jetevalstacks.begin();
138  it_jetevalstack != _jetevalstacks.end(); ++it_jetevalstack)
139  {
140  assert(it_jetevalstack->second);
141  it_jetevalstack->second->next_event(topNode);
142  }
143  if (_jettrutheval)
144  {
145  _jettrutheval->next_event(topNode);
146  }
147 
148  if (flag(kProcessTruthSpectrum))
149  {
150  if (Verbosity() >= 1)
151  {
152  std::cout << "QAG4SimulationJet::process_event - Process TruthSpectrum "
153  << _truth_jet << std::endl;
154  }
155  process_Spectrum(topNode, _truth_jet, false);
156  }
157 
158  if (flag(kProcessRecoSpectrum))
159  {
160  for (const auto& reco_jet : _reco_jets)
161  {
162  if (Verbosity() >= 1)
163  {
164  std::cout
165  << "QAG4SimulationJet::process_event - Process Reco jet spectrum "
166  << reco_jet << std::endl;
167  }
168  process_Spectrum(topNode, reco_jet, true);
169  }
170  }
171 
172  if (flag(kProcessTruthMatching))
173  {
174  for (const auto& reco_jet : _reco_jets)
175  {
176  if (Verbosity() >= 1)
177  {
178  std::cout
179  << "QAG4SimulationJet::process_event - Process Reco jet spectrum "
180  << reco_jet << std::endl;
181  }
182  process_TruthMatching(topNode, reco_jet);
183  }
184  }
185 
186  return Fun4AllReturnCodes::EVENT_OK;
187 }
188 
190 void QAG4SimulationJet::set_eta_range(double low, double high)
191 {
192  if (low > high)
193  {
194  std::swap(low, high);
195  }
196  assert(low < high); // eliminate zero range
197 
198  eta_range.first = low;
199  eta_range.second = high;
200 }
201 
204 TString
205 QAG4SimulationJet::get_eta_range_str(const char* eta_name) const
206 {
207  assert(eta_name);
208  return TString(Form("%.1f < %s < %.1f", eta_range.first, eta_name, eta_range.second));
209 }
210 
212 bool QAG4SimulationJet::jet_acceptance_cut(const Jet* jet) const
213 {
214  assert(jet);
215  bool eta_cut = (jet->get_eta() >= eta_range.first) && (jet->get_eta() <= eta_range.second);
216  return eta_cut;
217 }
218 
219 std::string
220 QAG4SimulationJet::get_histo_prefix(const std::string& src_jet_name,
221  const std::string& reco_jet_name)
222 {
223  std::string histo_prefix = "h_QAG4SimJet_";
224 
225  if (src_jet_name.length() > 0)
226  {
227  histo_prefix += src_jet_name;
228  histo_prefix += "_";
229  }
230  if (reco_jet_name.length() > 0)
231  {
232  histo_prefix += reco_jet_name;
233  histo_prefix += "_";
234  }
235 
236  return histo_prefix;
237 }
238 
239 int QAG4SimulationJet::Init_Spectrum(PHCompositeNode* /*topNode*/,
240  const std::string& jet_name)
241 {
242  Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
243  assert(hm);
244 
245  // normalization plot with counts
246  const int norm_size = 3;
247 
248  TH1D* h_norm = new TH1D(
249  TString(get_histo_prefix(jet_name)) + "Normalization",
250  " Normalization;Item;Count", norm_size, .5, norm_size + .5);
251  int i = 1;
252 
253  h_norm->GetXaxis()->SetBinLabel(i++, "Event");
254  h_norm->GetXaxis()->SetBinLabel(i++, "Inclusive Jets");
255  h_norm->GetXaxis()->SetBinLabel(i++, "Leading Jets");
256  assert(norm_size >= i - 1);
257  h_norm->GetXaxis()->LabelsOption("v");
258  hm->registerHisto(h_norm);
259 
260  hm->registerHisto(
261  new TH1F(
262  TString(get_histo_prefix(jet_name)) + "Leading_Et", //
263  TString(jet_name) + " leading jet Et, " + get_eta_range_str() + ";E_{T} (GeV)", 100, 0, 100));
264  hm->registerHisto(
265  new TH1F(
266  TString(get_histo_prefix(jet_name)) + "Leading_eta", //
267  TString(jet_name) + " leading jet #eta, " + get_eta_range_str() + ";#eta", 50, -1, 1));
268  hm->registerHisto(
269  new TH1F(
270  TString(get_histo_prefix(jet_name)) + "Leading_phi", //
271  TString(jet_name) + " leading jet #phi, " + get_eta_range_str() + ";#phi", 50, -M_PI, M_PI));
272 
273  TH1F* lcomp = new TH1F(
274  TString(get_histo_prefix(jet_name)) + "Leading_CompSize", //
275  TString(jet_name) + " leading jet # of component, " + get_eta_range_str() + ";Number of component;", 100, 1, 3000);
276  QAHistManagerDef::useLogBins(lcomp->GetXaxis());
277  hm->registerHisto(lcomp);
278 
279  hm->registerHisto(
280  new TH1F(
281  TString(get_histo_prefix(jet_name)) + "Leading_Mass", //
282  TString(jet_name) + " leading jet mass, " + get_eta_range_str() + ";Jet Mass (GeV);", 100, 0, 20));
283  hm->registerHisto(
284  new TH1F(
285  TString(get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio", //
286  TString(jet_name) + " leading jet EMCal ratio, " + get_eta_range_str() + ";Energy ratio CEMC/Total;", 100, 0, 1.01));
287  hm->registerHisto(
288  new TH1F(
289  TString(get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio", //
290  TString(jet_name) + " leading jet EMCal+HCal_{IN} ratio, " + get_eta_range_str() + ";Energy ratio (CEMC + HCALIN)/Total;",
291  100, 0, 1.01));
292 
293  // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
294  // It is only available for truth jets
295  hm->registerHisto(
296  new TH1F(
297  TString(get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio", //
298  TString(jet_name) + " leading jet leakage ratio, " + get_eta_range_str() + ";Energy ratio, Back leakage/Total;", 100,
299  0, 1.01));
300 
301  TH1F* h = new TH1F(
302  TString(get_histo_prefix(jet_name)) + "Inclusive_E", //
303  TString(jet_name) + " inclusive jet E, " + get_eta_range_str() + ";Total jet energy (GeV)", 100, 1e-3, 100);
304  QAHistManagerDef::useLogBins(h->GetXaxis());
305  hm->registerHisto(h);
306  hm->registerHisto(
307  new TH1F(
308  TString(get_histo_prefix(jet_name)) + "Inclusive_eta", //
309  TString(jet_name) + " inclusive jet #eta, " + get_eta_range_str() + ";#eta;Jet energy density", 50, -1, 1));
310  hm->registerHisto(
311  new TH1F(
312  TString(get_histo_prefix(jet_name)) + "Inclusive_phi", //
313  TString(jet_name) + " inclusive jet #phi, " + get_eta_range_str() + ";#phi;Jet energy density", 50, -M_PI, M_PI));
314 
315  return Fun4AllReturnCodes::EVENT_OK;
316 }
317 
318 int QAG4SimulationJet::process_Spectrum(PHCompositeNode* topNode,
319  const std::string& jet_name, const bool is_reco_jet)
320 {
321  JetContainer* jets = findNode::getClass<JetContainer>(topNode, jet_name.c_str());
322  if (!jets)
323  {
324  std::cout
325  << "QAG4SimulationJet::process_Spectrum - Error can not find DST JetContainer node "
326  << jet_name << std::endl;
327  exit(-1);
328  }
329 
330  Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
331  assert(hm);
332 
333  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto(get_histo_prefix(jet_name) + "Normalization"));
334  assert(h_norm);
335  h_norm->Fill("Event", 1);
336  h_norm->Fill("Inclusive Jets", jets->size());
337 
338  TH1F* ie = dynamic_cast<TH1F*>(hm->getHisto(
339  (get_histo_prefix(jet_name)) + "Inclusive_E" //
340  ));
341  assert(ie);
342  TH1F* ieta = dynamic_cast<TH1F*>(hm->getHisto(
343  (get_histo_prefix(jet_name)) + "Inclusive_eta" //
344  ));
345  assert(ieta);
346  TH1F* iphi = dynamic_cast<TH1F*>(hm->getHisto(
347  (get_histo_prefix(jet_name)) + "Inclusive_phi" //
348  ));
349  assert(iphi);
350 
351  Jet* leading_jet = nullptr;
352  double max_et = 0;
353  for (auto jet : *jets)
354  {
355  assert(jet);
356 
357  if (not jet_acceptance_cut(jet))
358  {
359  continue;
360  }
361 
362  if (jet->get_et() > max_et)
363  {
364  leading_jet = jet;
365  max_et = jet->get_et();
366  }
367 
368  ie->Fill(jet->get_e());
369  ieta->Fill(jet->get_eta());
370  iphi->Fill(jet->get_phi());
371  }
372 
373  if (leading_jet)
374  {
375  if (Verbosity() > 0)
376  {
377  std::cout
378  << "QAG4SimulationJet::process_Spectrum - processing leading jet with # comp = "
379  << leading_jet->size_comp() << std::endl;
380  leading_jet->identify();
381  }
382 
383  h_norm->Fill("Leading Jets", 1);
384 
385  TH1F* let = dynamic_cast<TH1F*>(hm->getHisto(
386  (get_histo_prefix(jet_name)) + "Leading_Et" //
387  ));
388  assert(let);
389  TH1F* leta = dynamic_cast<TH1F*>(hm->getHisto(
390  (get_histo_prefix(jet_name)) + "Leading_eta" //
391  ));
392  assert(leta);
393  TH1F* lphi = dynamic_cast<TH1F*>(hm->getHisto(
394  (get_histo_prefix(jet_name)) + "Leading_phi" //
395  ));
396  assert(lphi);
397 
398  TH1F* lcomp = dynamic_cast<TH1F*>(hm->getHisto(
399  (get_histo_prefix(jet_name)) + "Leading_CompSize" //
400  ));
401  assert(lcomp);
402  TH1F* lmass = dynamic_cast<TH1F*>(hm->getHisto(
403  (get_histo_prefix(jet_name)) + "Leading_Mass" //
404  ));
405  assert(lmass);
406  TH1F* lcemcr = dynamic_cast<TH1F*>(hm->getHisto(
407  (get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio" //
408  ));
409  assert(lcemcr);
410  TH1F* lemchcalr = dynamic_cast<TH1F*>(hm->getHisto(
411  (get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio" //
412  ));
413  assert(lemchcalr);
414  TH1F* lleak = dynamic_cast<TH1F*>(hm->getHisto(
415  (get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio" //
416  ));
417  assert(lleak);
418 
419  let->Fill(leading_jet->get_et());
420  leta->Fill(leading_jet->get_eta());
421  lphi->Fill(leading_jet->get_phi());
422  lcomp->Fill(leading_jet->size_comp());
423  lmass->Fill(leading_jet->get_mass());
424 
425  if (is_reco_jet)
426  { // this is a reco jet
427  jetevalstacks_map::iterator it_stack = _jetevalstacks.find(jet_name);
428  assert(it_stack != _jetevalstacks.end());
429  std::shared_ptr<JetEvalStack> eval_stack = it_stack->second;
430  assert(eval_stack);
431  JetRecoEval* recoeval = eval_stack->get_reco_eval();
432  assert(recoeval);
433 
434  if (Verbosity() >= VERBOSITY_A_LOT)
435  {
436  std::cout << __PRETTY_FUNCTION__ << "Leading Jet " << jet_name << ": ";
437  // leading_jet->identify();
438  std::cout << "CEMC_TOWER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) << std::endl;
439  std::cout << "CEMC_CLUSTER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::CEMC_CLUSTER) << std::endl;
440  std::cout << "HCALIN_TOWER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::HCALIN_TOWER) << std::endl;
441  std::cout << "HCALIN_CLUSTER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::HCALIN_CLUSTER) << std::endl;
442  std::cout << "leading_jet->get_e() = " << leading_jet->get_e() << std::endl;
443  }
444 
445  lcemcr->Fill( //
446  (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
447  + //
448  recoeval->get_energy_contribution(leading_jet,
449  Jet::CEMC_CLUSTER) //
450  ) /
451  leading_jet->get_e());
452 
453  lemchcalr->Fill( //
454  (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) //
455  + //
456  recoeval->get_energy_contribution(leading_jet,
457  Jet::CEMC_CLUSTER) //
458  + //
459  recoeval->get_energy_contribution(leading_jet,
460  Jet::HCALIN_TOWER) //
461  + //
462  recoeval->get_energy_contribution(leading_jet,
463  Jet::HCALIN_CLUSTER) //
464  ) /
465  leading_jet->get_e());
466  // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
467  // It is only available for truth jets
468  }
469  else
470  { // this is a truth jet
471  assert(_jettrutheval);
472 
473  double cemc_e = 0;
474  double hcalin_e = 0;
475  double bh_e = 0;
476 
477  std::set<PHG4Shower*> showers = _jettrutheval->all_truth_showers(leading_jet);
478 
479  for (auto shower : showers)
480  {
481  if (Verbosity() >= VERBOSITY_A_LOT)
482  {
483  std::cout << __PRETTY_FUNCTION__ << "Leading Truth Jet shower : ";
484  shower->identify();
485  }
486 
487  cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC"));
488  cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"));
489  cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"));
490 
491  hcalin_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_HCALIN"));
492  hcalin_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"));
493 
494  bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"));
495  bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"));
496  bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"));
497 
498  if (Verbosity() >= VERBOSITY_A_LOT)
499  {
500  // leading_jet->identify();
501  std::cout << "Shower cemc_e sum = "
502  << shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC"))
503  << " + "
504  << shower->get_edep(
505  PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"))
506  << " + "
507  << shower->get_edep(
508  PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"))
509  << std::endl;
510  std::cout << "Shower hcalin_e sum = "
511  << shower->get_edep(
512  PHG4HitDefs::get_volume_id("G4HIT_HCALIN"))
513  << " + "
514  << shower->get_edep(
515  PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"))
516  << std::endl;
517  std::cout << "Shower bh_e sum = "
518  << shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"))
519  << " + "
520  << shower->get_edep(
521  PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"))
522  << " + "
523  << shower->get_edep(
524  PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"))
525  << std::endl;
526  }
527  }
528 
529  if (Verbosity() >= VERBOSITY_A_LOT)
530  {
531  std::cout << "cemc_e sum = " << cemc_e << std::endl;
532  std::cout << "hcalin_e sum = " << hcalin_e << std::endl;
533  std::cout << "leading_jet->get_e() = " << leading_jet->get_e() << std::endl;
534  }
535 
536  lcemcr->Fill(cemc_e / leading_jet->get_e());
537  lemchcalr->Fill((cemc_e + hcalin_e) / leading_jet->get_e());
538  lleak->Fill(bh_e / leading_jet->get_e());
539  }
540  }
541 
542  return Fun4AllReturnCodes::EVENT_OK;
543 }
544 
545 int QAG4SimulationJet::Init_TruthMatching(PHCompositeNode* /*topNode*/,
546  const std::string& reco_jet_name)
547 {
548  Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
549  assert(hm);
550 
551  TH2F* h = new TH2F(
552  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Truth_Et", //
553  TString(reco_jet_name) + " Matching Count, " + get_eta_range_str() + ";E_{T, Truth} (GeV)", 20, 0, 100, 3, 0.5, 3.5);
554  h->GetYaxis()->SetBinLabel(1, "Total");
555  h->GetYaxis()->SetBinLabel(2, "Matched");
556  h->GetYaxis()->SetBinLabel(3, "Unique Matched");
557  hm->registerHisto(h);
558 
559  h = new TH2F(
560  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Reco_Et", //
561  TString(reco_jet_name) + " Matching Count, " + get_eta_range_str() + ";E_{T, Reco} (GeV)", 20, 0, 100, 3, 0.5, 3.5);
562  h->GetYaxis()->SetBinLabel(1, "Total");
563  h->GetYaxis()->SetBinLabel(2, "Matched");
564  h->GetYaxis()->SetBinLabel(3, "Unique Matched");
565  hm->registerHisto(h);
566 
567  h = new TH2F(
568  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEt", //
569  TString(reco_jet_name) + " E_{T} difference, " + get_eta_range_str() + ";E_{T, Truth} (GeV);E_{T, Reco} / E_{T, Truth}", 20, 0, 100, 100,
570  0, 2);
571  hm->registerHisto(h);
572 
573  h = new TH2F(
574  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dE", //
575  TString(reco_jet_name) + " Jet Energy Difference, " + get_eta_range_str() + ";E_{Truth} (GeV);E_{Reco} / E_{Truth}", 20, 0, 100, 100, 0, 2);
576  hm->registerHisto(h);
577 
578  h = new TH2F(
579  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEta", //
580  TString(reco_jet_name) + " #eta difference, " + get_eta_range_str() + ";E_{T, Truth} (GeV);#eta_{Reco} - #eta_{Truth}", 20, 0, 100, 200,
581  -.1, .1);
582  hm->registerHisto(h);
583 
584  h = new TH2F(
585  TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dPhi", //
586  TString(reco_jet_name) + " #phi difference, " + get_eta_range_str() + ";E_{T, Truth} (GeV);#phi_{Reco} - #phi_{Truth} (rad)", 20, 0, 100,
587  200, -.1, .1);
588  hm->registerHisto(h);
589 
590  return Fun4AllReturnCodes::EVENT_OK;
591 }
592 
593 int QAG4SimulationJet::process_TruthMatching(PHCompositeNode* topNode,
594  const std::string& reco_jet_name)
595 {
596  assert(_jet_match_dPhi > 0);
597  assert(_jet_match_dEta > 0);
598  assert(_jet_match_dE_Ratio > 0);
599 
600  Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
601  assert(hm);
602 
603  TH2F* Matching_Count_Truth_Et = dynamic_cast<TH2F*>(hm->getHisto(
604  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Truth_Et" //
605  ));
606  assert(Matching_Count_Truth_Et);
607  TH2F* Matching_Count_Reco_Et = dynamic_cast<TH2F*>(hm->getHisto(
608  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Reco_Et" //
609  ));
610  assert(Matching_Count_Reco_Et);
611  TH2F* Matching_dEt = dynamic_cast<TH2F*>(hm->getHisto(
612  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEt" //
613  ));
614  assert(Matching_dEt);
615  TH2F* Matching_dE = dynamic_cast<TH2F*>(hm->getHisto(
616  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dE" //
617  ));
618  assert(Matching_dE);
619  TH2F* Matching_dEta = dynamic_cast<TH2F*>(hm->getHisto(
620  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEta" //
621  ));
622  assert(Matching_dEta);
623  TH2F* Matching_dPhi = dynamic_cast<TH2F*>(hm->getHisto(
624  (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dPhi" //
625  ));
626  assert(Matching_dPhi);
627 
628  jetevalstacks_map::iterator it_stack = _jetevalstacks.find(reco_jet_name);
629  assert(it_stack != _jetevalstacks.end());
630  std::shared_ptr<JetEvalStack> eval_stack = it_stack->second;
631  assert(eval_stack);
632  JetRecoEval* recoeval = eval_stack->get_reco_eval();
633  assert(recoeval);
634 
635  // iterate over truth jets
636  JetContainer* truthjets = findNode::getClass<JetContainer>(topNode, _truth_jet);
637  if (!truthjets)
638  {
639  std::cout
640  << "QAG4SimulationJet::process_TruthMatching - Error can not find DST JetContainer node "
641  << _truth_jet << std::endl;
642  exit(-1);
643  }
644 
645  // search for leading truth
646  Jet* truthjet = nullptr;
647  double max_et = 0;
648  for (auto jet : *truthjets)
649  {
650  assert(jet);
651 
652  if (not jet_acceptance_cut(jet))
653  {
654  continue;
655  }
656 
657  if (jet->get_et() > max_et)
658  {
659  truthjet = jet;
660  max_et = jet->get_et();
661  }
662  }
663 
664 
665 
666  // match leading truth
667  if (truthjet)
668  {
669  if (Verbosity() > 1)
670  {
671  std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
672  << " process truth jet ";
673  truthjet->identify();
674  }
675 
676 
677  Matching_Count_Truth_Et->Fill(truthjet->get_et(), "Total", 1);
678  { // inclusive best energy match
679 
680  const Jet* recojet = recoeval->best_jet_from(truthjet);
681  if (Verbosity() > 1)
682  {
683  std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
684  << " inclusively matched with best reco jet: ";
685  recojet->identify();
686  }
687 
688 
689  if (recojet)
690  {
691  const double dPhi = recojet->get_phi() - truthjet->get_phi();
692  Matching_dPhi->Fill(truthjet->get_et(), dPhi);
693 
694  if (fabs(dPhi) < _jet_match_dPhi)
695  {
696  const double dEta = recojet->get_eta() - truthjet->get_eta();
697  Matching_dEta->Fill(truthjet->get_et(), dEta);
698 
699  if (fabs(dEta) < _jet_match_dEta)
700  {
701  const double Et_r = recojet->get_et() / (truthjet->get_et() + 1e-9);
702  const double E_r = recojet->get_e() / (truthjet->get_e() + 1e-9);
703  Matching_dEt->Fill(truthjet->get_et(), Et_r);
704  Matching_dE->Fill(truthjet->get_et(), E_r);
705 
706  if (fabs(E_r - 1) < _jet_match_dE_Ratio)
707  {
708  // matched in eta, phi and energy
709 
710  Matching_Count_Truth_Et->Fill(truthjet->get_et(),
711  "Matched", 1);
712  }
713 
714  } // if (fabs(dEta) < 0.1)
715 
716  } // if (fabs(dPhi) < 0.1)
717 
718  } // if (recojet)
719  } // inclusive best energy match
720  { // unique match
721 
722  const Jet* recojet = recoeval->unique_reco_jet_from_truth(truthjet);
723  if (recojet)
724  {
725  if (Verbosity() > 1)
726  {
727  std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
728  << " uniquely matched with reco jet: ";
729  recojet->identify();
730  }
731 
732  const double dPhi = recojet->get_phi() - truthjet->get_phi();
733 
734  if (fabs(dPhi) < _jet_match_dPhi)
735  {
736  const double dEta = recojet->get_eta() - truthjet->get_eta();
737 
738  if (fabs(dEta) < _jet_match_dEta)
739  {
740  const double E_r = recojet->get_e() / (truthjet->get_e() + 1e-9);
741  if (fabs(E_r - 1) < _jet_match_dE_Ratio)
742  {
743  // matched in eta, phi and energy
744 
745  Matching_Count_Truth_Et->Fill(truthjet->get_et(),
746  "Unique Matched", 1);
747  }
748 
749  } // if (fabs(dEta) < 0.1)
750 
751  } // if (fabs(dPhi) < 0.1)
752 
753  } // if (recojet)
754  } // unique match
755 
756  } // if (truthjet)
757 
758  // next for reco jets
759  JetContainer* recojets = findNode::getClass<JetContainer>(topNode, reco_jet_name);
760  if (!recojets)
761  {
762  std::cout
763  << "QAG4SimulationJet::process_TruthMatching - Error can not find DST JetContainer node "
764  << reco_jet_name << std::endl;
765  exit(-1);
766  }
767 
768  // search for leading reco jet
769  Jet* recojet = nullptr;
770  max_et = 0;
771  for (auto jet : *recojets)
772  {
773  assert(jet);
774 
775  if (not jet_acceptance_cut(jet))
776  {
777  continue;
778  }
779 
780  if (jet->get_et() > max_et)
781  {
782  recojet = jet;
783  max_et = jet->get_et();
784  }
785  }
786 
787  // match leading reco jet
788  if (recojet)
789  {
790  if (Verbosity() > 1)
791  {
792  std::cout << "QAG4SimulationJet::process_TruthMatching - " << reco_jet_name
793  << " process reco jet ";
794  recojet->identify();
795  }
796 
797  Matching_Count_Reco_Et->Fill(recojet->get_et(), "Total", 1);
798 
799  { // inclusive best energy match
800  Jet* truthjet1 = recoeval->max_truth_jet_by_energy(recojet);
801  if (truthjet1)
802  {
803  const double dPhi = recojet->get_phi() - truthjet1->get_phi();
804  if (fabs(dPhi) < _jet_match_dPhi)
805  {
806  const double dEta = recojet->get_eta() - truthjet1->get_eta();
807  if (fabs(dEta) < _jet_match_dEta)
808  {
809  const double E_r = recojet->get_e() / (truthjet1->get_e() + 1e-9);
810 
811  if (fabs(E_r - 1) < _jet_match_dE_Ratio)
812  {
813  // matched in eta, phi and energy
814 
815  Matching_Count_Reco_Et->Fill(recojet->get_et(),
816  "Unique Matched", 1);
817  }
818 
819  } // if (fabs(dEta) < 0.1)
820 
821  } // if (fabs(dPhi) < 0.1)
822 
823  } // if (truthjet1)
824  } // inclusive best energy match
825 
826  { // unique match
827  Jet* truthjet2 = recoeval->unique_truth_jet_from_reco(recojet);
828  if (truthjet2)
829  {
830  const double dPhi = recojet->get_phi() - truthjet2->get_phi();
831  if (fabs(dPhi) < _jet_match_dPhi)
832  {
833  const double dEta = recojet->get_eta() - truthjet2->get_eta();
834  if (fabs(dEta) < _jet_match_dEta)
835  {
836  const double E_r = recojet->get_e() / (truthjet2->get_e() + 1e-9);
837 
838  if (fabs(E_r - 1) < _jet_match_dE_Ratio)
839  {
840  // matched in eta, phi and energy
841 
842  Matching_Count_Reco_Et->Fill(recojet->get_et(),
843  "Matched", 1);
844  }
845 
846  } // if (fabs(dEta) < 0.1)
847 
848  } // if (fabs(dPhi) < 0.1)
849 
850  } // if (truthjet2)
851  } // unique match
852 
853  } // if (recojet)
854 
855  return Fun4AllReturnCodes::EVENT_OK;
856 }