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