QA/SimulationModules/QAG4SimulationJet.cc

0001 #include "QAG4SimulationJet.h"
0002
0003 #include <TString.h>
0004 #include <qautils/QAHistManagerDef.h>
0005
0006 #include <g4eval/JetEvalStack.h>
0007 #include <g4eval/JetRecoEval.h>
0008 #include <g4eval/JetTruthEval.h>
0009
0010 #include <jetbase/Jet.h>
0011 #include <jetbase/JetContainer.h>
0012
0013 #include <g4main/PHG4HitDefs.h>
0014 #include <g4main/PHG4Shower.h>
0015
0016 #include <fun4all/Fun4AllBase.h>
0017 #include <fun4all/Fun4AllHistoManager.h>
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019 #include <fun4all/SubsysReco.h>  // for SubsysReco
0020
0021 #include <phool/getClass.h>
0022
0023 #include <boost/format.hpp>
0024
0025 #include <TAxis.h>
0026 #include <TH1.h>
0027 #include <TH2.h>
0028
0029 #include <cassert>
0030 #include <cmath>
0031 #include <cstdlib>
0032 #include <iostream>
0033 #include <memory>
0034 #include <set>
0035 #include <string>
0036 #include <utility>
0037
0038 QAG4SimulationJet::QAG4SimulationJet(const std::string& truth_jet,
0039                                      enu_flags flags)
0040   : SubsysReco("QAG4SimulationJet_" + truth_jet)
0041   , _jetevalstacks()
0042   , _truth_jet(truth_jet)
0043   , _reco_jets()
0044   , _flags(flags)
0045   , eta_range(-1, 1)
0046   , _jet_match_dEta(.1)
0047   , _jet_match_dPhi(.1)
0048   , _jet_match_dE_Ratio(.5)
0049 {
0050 }
0051
0052 int QAG4SimulationJet::InitRun(PHCompositeNode* topNode)
0053 {
0054   if (flag(kProcessTruthMatching) || flag(kProcessRecoSpectrum))
0055   {
0056     for (const auto& reco_jet : _reco_jets)
0057     {
0058       jetevalstacks_map::iterator const it_jetevalstack = _jetevalstacks.find(
0059           reco_jet);
0060
0061       if (it_jetevalstack == _jetevalstacks.end())
0062       {
0063         _jetevalstacks[reco_jet] = std::shared_ptr<JetEvalStack>(new JetEvalStack(topNode, reco_jet, _truth_jet));
0064         assert(_jetevalstacks[reco_jet]);
0065         _jetevalstacks[reco_jet]->set_strict(true);
0066         _jetevalstacks[reco_jet]->set_verbosity(Verbosity() + 1);
0067       }
0068       else
0069       {
0070         assert(it_jetevalstack->second);
0071       }
0072     }
0073   }
0074
0075   if (flag(kProcessTruthSpectrum))
0076   {
0077     if (not _jettrutheval)
0078     {
0079       _jettrutheval = std::shared_ptr<JetTruthEval>(new JetTruthEval(topNode, _truth_jet));
0080     }
0081
0082     assert(_jettrutheval);
0083     _jettrutheval->set_strict(true);
0084     _jettrutheval->set_verbosity(Verbosity() + 1);
0085   }
0086
0087   return Fun4AllReturnCodes::EVENT_OK;
0088 }
0089
0090 int QAG4SimulationJet::Init(PHCompositeNode* topNode)
0091 {
0092   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0093   assert(hm);
0094
0095   if (flag(kProcessTruthSpectrum))
0096   {
0097     if (Verbosity() >= 1)
0098     {
0099       std::cout << "QAG4SimulationJet::Init - Process TruthSpectrum " << _truth_jet
0100                 << std::endl;
0101     }
0102     Init_Spectrum(topNode, _truth_jet);
0103   }
0104
0105   if (flag(kProcessRecoSpectrum))
0106   {
0107     for (const auto& reco_jet : _reco_jets)
0108     {
0109       if (Verbosity() >= 1)
0110       {
0111         std::cout << "QAG4SimulationJet::Init - Process Reco jet spectrum "
0112                   << reco_jet << std::endl;
0113       }
0114       Init_Spectrum(topNode, reco_jet);
0115     }
0116   }
0117
0118   if (flag(kProcessTruthMatching))
0119   {
0120     for (const auto& reco_jet : _reco_jets)
0121     {
0122       if (Verbosity() >= 1)
0123       {
0124         std::cout << "QAG4SimulationJet::Init - Process Reco jet spectrum "
0125                   << reco_jet << std::endl;
0126       }
0127       Init_TruthMatching(topNode, reco_jet);
0128     }
0129   }
0130
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133
0134 int QAG4SimulationJet::process_event(PHCompositeNode* topNode)
0135 {
0136   if (Verbosity() > 2)
0137   {
0138     std::cout << "QAG4SimulationJet::process_event() entered" << std::endl;
0139   }
0140
0141   for (jetevalstacks_map::iterator it_jetevalstack = _jetevalstacks.begin();
0142        it_jetevalstack != _jetevalstacks.end(); ++it_jetevalstack)
0143   {
0144     assert(it_jetevalstack->second);
0145     it_jetevalstack->second->next_event(topNode);
0146   }
0147   if (_jettrutheval)
0148   {
0149     _jettrutheval->next_event(topNode);
0150   }
0151
0152   if (flag(kProcessTruthSpectrum))
0153   {
0154     if (Verbosity() >= 1)
0155     {
0156       std::cout << "QAG4SimulationJet::process_event - Process TruthSpectrum "
0157                 << _truth_jet << std::endl;
0158     }
0159     process_Spectrum(topNode, _truth_jet, false);
0160   }
0161
0162   if (flag(kProcessRecoSpectrum))
0163   {
0164     for (const auto& reco_jet : _reco_jets)
0165     {
0166       if (Verbosity() >= 1)
0167       {
0168         std::cout
0169             << "QAG4SimulationJet::process_event - Process Reco jet spectrum "
0170             << reco_jet << std::endl;
0171       }
0172       process_Spectrum(topNode, reco_jet, true);
0173     }
0174   }
0175
0176   if (flag(kProcessTruthMatching))
0177   {
0178     for (const auto& reco_jet : _reco_jets)
0179     {
0180       if (Verbosity() >= 1)
0181       {
0182         std::cout
0183             << "QAG4SimulationJet::process_event - Process Reco jet spectrum "
0184             << reco_jet << std::endl;
0185       }
0186       process_TruthMatching(topNode, reco_jet);
0187     }
0188   }
0189
0190   return Fun4AllReturnCodes::EVENT_OK;
0191 }
0192
0193 //! set eta range
0194 void QAG4SimulationJet::set_eta_range(double low, double high)
0195 {
0196   if (low > high)
0197   {
0198     std::swap(low, high);
0199   }
0200   assert(low < high);  // eliminate zero range
0201
0202   eta_range.first = low;
0203   eta_range.second = high;
0204 }
0205
0206 //! string description of eta range
0207 //! @return TString as ROOT likes
0208 TString
0209 QAG4SimulationJet::get_eta_range_str(const char* eta_name) const
0210 {
0211   assert(eta_name);
0212   return TString((boost::format("%.1f < %s < %.1f") % eta_range.first % eta_name % eta_range.second).str().c_str());
0213 }
0214
0215 //! acceptance cut on jet object
0216 bool QAG4SimulationJet::jet_acceptance_cut(const Jet* jet) const
0217 {
0218   assert(jet);
0219   bool const eta_cut = (jet->get_eta() >= eta_range.first) && (jet->get_eta() <= eta_range.second);
0220   return eta_cut;
0221 }
0222
0223 std::string
0224 QAG4SimulationJet::get_histo_prefix(const std::string& src_jet_name,
0225                                     const std::string& reco_jet_name)
0226 {
0227   std::string histo_prefix = "h_QAG4SimJet_";
0228
0229   if (src_jet_name.length() > 0)
0230   {
0231     histo_prefix += src_jet_name;
0232     histo_prefix += "_";
0233   }
0234   if (reco_jet_name.length() > 0)
0235   {
0236     histo_prefix += reco_jet_name;
0237     histo_prefix += "_";
0238   }
0239
0240   return histo_prefix;
0241 }
0242
0243 int QAG4SimulationJet::Init_Spectrum(PHCompositeNode* /*topNode*/,
0244                                      const std::string& jet_name)
0245 {
0246   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0247   assert(hm);
0248
0249   // normalization plot with counts
0250   const int norm_size = 3;
0251
0252   TH1D* h_norm = new TH1D(
0253       TString(get_histo_prefix(jet_name)) + "Normalization",
0254       " Normalization;Item;Count", norm_size, .5, norm_size + .5);
0255   int i = 1;
0256
0257   h_norm->GetXaxis()->SetBinLabel(i++, "Event");
0258   h_norm->GetXaxis()->SetBinLabel(i++, "Inclusive Jets");
0259   h_norm->GetXaxis()->SetBinLabel(i++, "Leading Jets");
0260   assert(norm_size >= i - 1);
0261   h_norm->GetXaxis()->LabelsOption("v");
0262   hm->registerHisto(h_norm);
0263
0264   hm->registerHisto(
0265       new TH1F(
0266           TString(get_histo_prefix(jet_name)) + "Leading_Et",  //
0267           TString(jet_name) + " leading jet Et, " + get_eta_range_str() + ";E_{T} (GeV)", 100, 0, 100));
0268   hm->registerHisto(
0269       new TH1F(
0270           TString(get_histo_prefix(jet_name)) + "Leading_eta",  //
0271           TString(jet_name) + " leading jet #eta, " + get_eta_range_str() + ";#eta", 50, -1, 1));
0272   hm->registerHisto(
0273       new TH1F(
0274           TString(get_histo_prefix(jet_name)) + "Leading_phi",  //
0275           TString(jet_name) + " leading jet #phi, " + get_eta_range_str() + ";#phi", 50, -M_PI, M_PI));
0276
0277   TH1F* lcomp = new TH1F(
0278       TString(get_histo_prefix(jet_name)) + "Leading_CompSize",  //
0279       TString(jet_name) + " leading jet # of component, " + get_eta_range_str() + ";Number of component;", 100, 1, 3000);
0280   QAHistManagerDef::useLogBins(lcomp->GetXaxis());
0281   hm->registerHisto(lcomp);
0282
0283   hm->registerHisto(
0284       new TH1F(
0285           TString(get_histo_prefix(jet_name)) + "Leading_Mass",  //
0286           TString(jet_name) + " leading jet mass, " + get_eta_range_str() + ";Jet Mass (GeV);", 100, 0, 20));
0287   hm->registerHisto(
0288       new TH1F(
0289           TString(get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio",  //
0290           TString(jet_name) + " leading jet EMCal ratio, " + get_eta_range_str() + ";Energy ratio CEMC/Total;", 100, 0, 1.01));
0291   hm->registerHisto(
0292       new TH1F(
0293           TString(get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio",  //
0294           TString(jet_name) + " leading jet EMCal+HCal_{IN} ratio, " + get_eta_range_str() + ";Energy ratio (CEMC + HCALIN)/Total;",
0295           100, 0, 1.01));
0296
0297   // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
0298   // It is only available for truth jets
0299   hm->registerHisto(
0300       new TH1F(
0301           TString(get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio",  //
0302           TString(jet_name) + " leading jet leakage ratio, " + get_eta_range_str() + ";Energy ratio, Back leakage/Total;", 100,
0303           0, 1.01));
0304
0305   TH1F* h = new TH1F(
0306       TString(get_histo_prefix(jet_name)) + "Inclusive_E",  //
0307       TString(jet_name) + " inclusive jet E, " + get_eta_range_str() + ";Total jet energy (GeV)", 100, 1e-3, 100);
0308   QAHistManagerDef::useLogBins(h->GetXaxis());
0309   hm->registerHisto(h);
0310   hm->registerHisto(
0311       new TH1F(
0312           TString(get_histo_prefix(jet_name)) + "Inclusive_eta",  //
0313           TString(jet_name) + " inclusive jet #eta, " + get_eta_range_str() + ";#eta;Jet energy density", 50, -1, 1));
0314   hm->registerHisto(
0315       new TH1F(
0316           TString(get_histo_prefix(jet_name)) + "Inclusive_phi",  //
0317           TString(jet_name) + " inclusive jet #phi, " + get_eta_range_str() + ";#phi;Jet energy density", 50, -M_PI, M_PI));
0318
0319   return Fun4AllReturnCodes::EVENT_OK;
0320 }
0321
0322 int QAG4SimulationJet::process_Spectrum(PHCompositeNode* topNode,
0323                                         const std::string& jet_name, const bool is_reco_jet)
0324 {
0325   JetContainer* jets = findNode::getClass<JetContainer>(topNode, jet_name.c_str());
0326   if (!jets)
0327   {
0328     std::cout
0329         << "QAG4SimulationJet::process_Spectrum - Error can not find DST JetContainer node "
0330         << jet_name << std::endl;
0331     exit(-1);
0332   }
0333
0334   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0335   assert(hm);
0336
0337   TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto(get_histo_prefix(jet_name) + "Normalization"));
0338   assert(h_norm);
0339   h_norm->Fill("Event", 1);
0340   h_norm->Fill("Inclusive Jets", jets->size());
0341
0342   TH1F* ie = dynamic_cast<TH1F*>(hm->getHisto(
0343       (get_histo_prefix(jet_name)) + "Inclusive_E"  //
0344       ));
0345   assert(ie);
0346   TH1F* ieta = dynamic_cast<TH1F*>(hm->getHisto(
0347       (get_histo_prefix(jet_name)) + "Inclusive_eta"  //
0348       ));
0349   assert(ieta);
0350   TH1F* iphi = dynamic_cast<TH1F*>(hm->getHisto(
0351       (get_histo_prefix(jet_name)) + "Inclusive_phi"  //
0352       ));
0353   assert(iphi);
0354
0355   Jet* leading_jet = nullptr;
0356   double max_et = 0;
0357   for (auto jet : *jets)
0358   {
0359     assert(jet);
0360
0361     if (not jet_acceptance_cut(jet))
0362     {
0363       continue;
0364     }
0365
0366     if (jet->get_et() > max_et)
0367     {
0368       leading_jet = jet;
0369       max_et = jet->get_et();
0370     }
0371
0372     ie->Fill(jet->get_e());
0373     ieta->Fill(jet->get_eta());
0374     iphi->Fill(jet->get_phi());
0375   }
0376
0377   if (leading_jet)
0378   {
0379     if (Verbosity() > 0)
0380     {
0381       std::cout
0382           << "QAG4SimulationJet::process_Spectrum - processing leading jet with # comp = "
0383           << leading_jet->size_comp() << std::endl;
0384       leading_jet->identify();
0385     }
0386
0387     h_norm->Fill("Leading Jets", 1);
0388
0389     TH1F* let = dynamic_cast<TH1F*>(hm->getHisto(
0390         (get_histo_prefix(jet_name)) + "Leading_Et"  //
0391         ));
0392     assert(let);
0393     TH1F* leta = dynamic_cast<TH1F*>(hm->getHisto(
0394         (get_histo_prefix(jet_name)) + "Leading_eta"  //
0395         ));
0396     assert(leta);
0397     TH1F* lphi = dynamic_cast<TH1F*>(hm->getHisto(
0398         (get_histo_prefix(jet_name)) + "Leading_phi"  //
0399         ));
0400     assert(lphi);
0401
0402     TH1F* lcomp = dynamic_cast<TH1F*>(hm->getHisto(
0403         (get_histo_prefix(jet_name)) + "Leading_CompSize"  //
0404         ));
0405     assert(lcomp);
0406     TH1F* lmass = dynamic_cast<TH1F*>(hm->getHisto(
0407         (get_histo_prefix(jet_name)) + "Leading_Mass"  //
0408         ));
0409     assert(lmass);
0410     TH1F* lcemcr = dynamic_cast<TH1F*>(hm->getHisto(
0411         (get_histo_prefix(jet_name)) + "Leading_CEMC_Ratio"  //
0412         ));
0413     assert(lcemcr);
0414     TH1F* lemchcalr = dynamic_cast<TH1F*>(hm->getHisto(
0415         (get_histo_prefix(jet_name)) + "Leading_CEMC_HCalIN_Ratio"  //
0416         ));
0417     assert(lemchcalr);
0418     TH1F* lleak = dynamic_cast<TH1F*>(hm->getHisto(
0419         (get_histo_prefix(jet_name)) + "Leading_Leakage_Ratio"  //
0420         ));
0421     assert(lleak);
0422
0423     let->Fill(leading_jet->get_et());
0424     leta->Fill(leading_jet->get_eta());
0425     lphi->Fill(leading_jet->get_phi());
0426     lcomp->Fill(leading_jet->size_comp());
0427     lmass->Fill(leading_jet->get_mass());
0428
0429     if (is_reco_jet)
0430     {  // this is a reco jet
0431       jetevalstacks_map::iterator const it_stack = _jetevalstacks.find(jet_name);
0432       assert(it_stack != _jetevalstacks.end());
0433       std::shared_ptr<JetEvalStack> const eval_stack = it_stack->second;
0434       assert(eval_stack);
0435       JetRecoEval* recoeval = eval_stack->get_reco_eval();
0436       assert(recoeval);
0437
0438       if (Verbosity() >= VERBOSITY_A_LOT)
0439       {
0440         std::cout << __PRETTY_FUNCTION__ << "Leading Jet " << jet_name << ": ";
0441         //            leading_jet->identify();
0442         std::cout << "CEMC_TOWER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER) << std::endl;
0443         std::cout << "CEMC_CLUSTER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::CEMC_CLUSTER) << std::endl;
0444         std::cout << "HCALIN_TOWER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::HCALIN_TOWER) << std::endl;
0445         std::cout << "HCALIN_CLUSTER sum = " << recoeval->get_energy_contribution(leading_jet, Jet::HCALIN_CLUSTER) << std::endl;
0446         std::cout << "leading_jet->get_e() = " << leading_jet->get_e() << std::endl;
0447       }
0448
0449       lcemcr->Fill(                                                         //
0450           (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER)  //
0451            +                                                                //
0452            recoeval->get_energy_contribution(leading_jet,
0453                                              Jet::CEMC_CLUSTER)  //
0454            ) /
0455           leading_jet->get_e());
0456
0457       lemchcalr->Fill(                                                      //
0458           (recoeval->get_energy_contribution(leading_jet, Jet::CEMC_TOWER)  //
0459            +                                                                //
0460            recoeval->get_energy_contribution(leading_jet,
0461                                              Jet::CEMC_CLUSTER)  //
0462            +                                                     //
0463            recoeval->get_energy_contribution(leading_jet,
0464                                              Jet::HCALIN_TOWER)  //
0465            +                                                     //
0466            recoeval->get_energy_contribution(leading_jet,
0467                                              Jet::HCALIN_CLUSTER)  //
0468            ) /
0469           leading_jet->get_e());
0470       // reco jet has no definition for leakages, since leakage is not reconstructed as part of jet energy.
0471       // It is only available for truth jets
0472     }
0473     else
0474     {  // this is a truth jet
0475       assert(_jettrutheval);
0476
0477       double cemc_e = 0;
0478       double hcalin_e = 0;
0479       double bh_e = 0;
0480
0481       std::set<PHG4Shower*> const showers = _jettrutheval->all_truth_showers(leading_jet);
0482
0483       for (auto shower : showers)
0484       {
0485         if (Verbosity() >= VERBOSITY_A_LOT)
0486         {
0487           std::cout << __PRETTY_FUNCTION__ << "Leading Truth Jet shower : ";
0488           shower->identify();
0489         }
0490
0491         cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC"));
0492         cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"));
0493         cemc_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"));
0494
0495         hcalin_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_HCALIN"));
0496         hcalin_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"));
0497
0498         bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"));
0499         bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"));
0500         bh_e += shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"));
0501
0502         if (Verbosity() >= VERBOSITY_A_LOT)
0503         {
0504           //            leading_jet->identify();
0505           std::cout << "Shower cemc_e sum = "
0506                     << shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_CEMC"))
0507                     << " + "
0508                     << shower->get_edep(
0509                            PHG4HitDefs::get_volume_id("G4HIT_CEMC_ELECTRONICS"))
0510                     << " + "
0511                     << shower->get_edep(
0512                            PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_CEMC"))
0513                     << std::endl;
0514           std::cout << "Shower hcalin_e sum = "
0515                     << shower->get_edep(
0516                            PHG4HitDefs::get_volume_id("G4HIT_HCALIN"))
0517                     << " + "
0518                     << shower->get_edep(
0519                            PHG4HitDefs::get_volume_id("G4HIT_ABSORBER_HCALIN"))
0520                     << std::endl;
0521           std::cout << "Shower bh_e sum = "
0522                     << shower->get_edep(PHG4HitDefs::get_volume_id("G4HIT_BH_1"))
0523                     << " + "
0524                     << shower->get_edep(
0525                            PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_PLUS"))
0526                     << " + "
0527                     << shower->get_edep(
0528                            PHG4HitDefs::get_volume_id("G4HIT_BH_FORWARD_NEG"))
0529                     << std::endl;
0530         }
0531       }
0532
0533       if (Verbosity() >= VERBOSITY_A_LOT)
0534       {
0535         std::cout << "cemc_e sum = " << cemc_e << std::endl;
0536         std::cout << "hcalin_e sum = " << hcalin_e << std::endl;
0537         std::cout << "leading_jet->get_e() = " << leading_jet->get_e() << std::endl;
0538       }
0539
0540       lcemcr->Fill(cemc_e / leading_jet->get_e());
0541       lemchcalr->Fill((cemc_e + hcalin_e) / leading_jet->get_e());
0542       lleak->Fill(bh_e / leading_jet->get_e());
0543     }
0544   }
0545
0546   return Fun4AllReturnCodes::EVENT_OK;
0547 }
0548
0549 int QAG4SimulationJet::Init_TruthMatching(PHCompositeNode* /*topNode*/,
0550                                           const std::string& reco_jet_name)
0551 {
0552   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0553   assert(hm);
0554
0555   TH2F* h = new TH2F(
0556       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Truth_Et",  //
0557       TString(reco_jet_name) + " Matching Count, " + get_eta_range_str() + ";E_{T, Truth} (GeV)", 20, 0, 100, 3, 0.5, 3.5);
0558   h->GetYaxis()->SetBinLabel(1, "Total");
0559   h->GetYaxis()->SetBinLabel(2, "Matched");
0560   h->GetYaxis()->SetBinLabel(3, "Unique Matched");
0561   hm->registerHisto(h);
0562
0563   h = new TH2F(
0564       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Reco_Et",  //
0565       TString(reco_jet_name) + " Matching Count, " + get_eta_range_str() + ";E_{T, Reco} (GeV)", 20, 0, 100, 3, 0.5, 3.5);
0566   h->GetYaxis()->SetBinLabel(1, "Total");
0567   h->GetYaxis()->SetBinLabel(2, "Matched");
0568   h->GetYaxis()->SetBinLabel(3, "Unique Matched");
0569   hm->registerHisto(h);
0570
0571   h = new TH2F(
0572       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEt",  //
0573       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,
0574       0, 2);
0575   hm->registerHisto(h);
0576
0577   h = new TH2F(
0578       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dE",  //
0579       TString(reco_jet_name) + " Jet Energy Difference, " + get_eta_range_str() + ";E_{Truth} (GeV);E_{Reco} / E_{Truth}", 20, 0, 100, 100, 0, 2);
0580   hm->registerHisto(h);
0581
0582   h = new TH2F(
0583       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEta",  //
0584       TString(reco_jet_name) + " #eta difference, " + get_eta_range_str() + ";E_{T, Truth} (GeV);#eta_{Reco} - #eta_{Truth}", 20, 0, 100, 200,
0585       -.1, .1);
0586   hm->registerHisto(h);
0587
0588   h = new TH2F(
0589       TString(get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dPhi",  //
0590       TString(reco_jet_name) + " #phi difference, " + get_eta_range_str() + ";E_{T, Truth} (GeV);#phi_{Reco} - #phi_{Truth} (rad)", 20, 0, 100,
0591       200, -.1, .1);
0592   hm->registerHisto(h);
0593
0594   return Fun4AllReturnCodes::EVENT_OK;
0595 }
0596
0597 int QAG4SimulationJet::process_TruthMatching(PHCompositeNode* topNode,
0598                                              const std::string& reco_jet_name)
0599 {
0600   assert(_jet_match_dPhi > 0);
0601   assert(_jet_match_dEta > 0);
0602   assert(_jet_match_dE_Ratio > 0);
0603
0604   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0605   assert(hm);
0606
0607   TH2F* Matching_Count_Truth_Et = dynamic_cast<TH2F*>(hm->getHisto(
0608       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Truth_Et"  //
0609       ));
0610   assert(Matching_Count_Truth_Et);
0611   TH2F* Matching_Count_Reco_Et = dynamic_cast<TH2F*>(hm->getHisto(
0612       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_Count_Reco_Et"  //
0613       ));
0614   assert(Matching_Count_Reco_Et);
0615   TH2F* Matching_dEt = dynamic_cast<TH2F*>(hm->getHisto(
0616       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEt"  //
0617       ));
0618   assert(Matching_dEt);
0619   TH2F* Matching_dE = dynamic_cast<TH2F*>(hm->getHisto(
0620       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dE"  //
0621       ));
0622   assert(Matching_dE);
0623   TH2F* Matching_dEta = dynamic_cast<TH2F*>(hm->getHisto(
0624       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dEta"  //
0625       ));
0626   assert(Matching_dEta);
0627   TH2F* Matching_dPhi = dynamic_cast<TH2F*>(hm->getHisto(
0628       (get_histo_prefix(_truth_jet, reco_jet_name)) + "Matching_dPhi"  //
0629       ));
0630   assert(Matching_dPhi);
0631
0632   jetevalstacks_map::iterator const it_stack = _jetevalstacks.find(reco_jet_name);
0633   assert(it_stack != _jetevalstacks.end());
0634   std::shared_ptr<JetEvalStack> const eval_stack = it_stack->second;
0635   assert(eval_stack);
0636   JetRecoEval* recoeval = eval_stack->get_reco_eval();
0637   assert(recoeval);
0638
0639   // iterate over truth jets
0640   JetContainer* truthjets = findNode::getClass<JetContainer>(topNode, _truth_jet);
0641   if (!truthjets)
0642   {
0643     std::cout
0644         << "QAG4SimulationJet::process_TruthMatching - Error can not find DST JetContainer node "
0645         << _truth_jet << std::endl;
0646     exit(-1);
0647   }
0648
0649   // search for leading truth
0650   Jet* truthjet = nullptr;
0651   double max_et = 0;
0652   for (auto jet : *truthjets)
0653   {
0654     assert(jet);
0655
0656     if (not jet_acceptance_cut(jet))
0657     {
0658       continue;
0659     }
0660
0661     if (jet->get_et() > max_et)
0662     {
0663       truthjet = jet;
0664       max_et = jet->get_et();
0665     }
0666   }
0667
0668   // match leading truth
0669   if (truthjet)
0670   {
0671     if (Verbosity() > 1)
0672     {
0673       std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
0674                 << " process truth jet ";
0675       truthjet->identify();
0676     }
0677
0678     Matching_Count_Truth_Et->Fill(truthjet->get_et(), "Total", 1);
0679     {  // inclusive best energy match
0680
0681       const Jet* recojet = recoeval->best_jet_from(truthjet);
0682       if (Verbosity() > 1)
0683       {
0684         std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
0685                   << " inclusively matched with best reco jet: ";
0686         recojet->identify();
0687       }
0688
0689       if (recojet)
0690       {
0691         const double dPhi = recojet->get_phi() - truthjet->get_phi();
0692         Matching_dPhi->Fill(truthjet->get_et(), dPhi);
0693
0694         if (fabs(dPhi) < _jet_match_dPhi)
0695         {
0696           const double dEta = recojet->get_eta() - truthjet->get_eta();
0697           Matching_dEta->Fill(truthjet->get_et(), dEta);
0698
0699           if (fabs(dEta) < _jet_match_dEta)
0700           {
0701             const double Et_r = recojet->get_et() / (truthjet->get_et() + 1e-9);
0702             const double E_r = recojet->get_e() / (truthjet->get_e() + 1e-9);
0703             Matching_dEt->Fill(truthjet->get_et(), Et_r);
0704             Matching_dE->Fill(truthjet->get_et(), E_r);
0705
0706             if (fabs(E_r - 1) < _jet_match_dE_Ratio)
0707             {
0708               // matched in eta, phi and energy
0709
0710               Matching_Count_Truth_Et->Fill(truthjet->get_et(),
0711                                             "Matched", 1);
0712             }
0713
0714           }  //  if (fabs(dEta) < 0.1)
0715
0716         }  // if (fabs(dPhi) < 0.1)
0717
0718       }  //       if (recojet)
0719     }  // inclusive best energy match
0720     {  // unique match
0721
0722       const Jet* recojet = recoeval->unique_reco_jet_from_truth(truthjet);
0723       if (recojet)
0724       {
0725         if (Verbosity() > 1)
0726         {
0727           std::cout << "QAG4SimulationJet::process_TruthMatching - " << _truth_jet
0728                     << " uniquely matched with reco jet: ";
0729           recojet->identify();
0730         }
0731
0732         const double dPhi = recojet->get_phi() - truthjet->get_phi();
0733
0734         if (fabs(dPhi) < _jet_match_dPhi)
0735         {
0736           const double dEta = recojet->get_eta() - truthjet->get_eta();
0737
0738           if (fabs(dEta) < _jet_match_dEta)
0739           {
0740             const double E_r = recojet->get_e() / (truthjet->get_e() + 1e-9);
0741             if (fabs(E_r - 1) < _jet_match_dE_Ratio)
0742             {
0743               // matched in eta, phi and energy
0744
0745               Matching_Count_Truth_Et->Fill(truthjet->get_et(),
0746                                             "Unique Matched", 1);
0747             }
0748
0749           }  //  if (fabs(dEta) < 0.1)
0750
0751         }  // if (fabs(dPhi) < 0.1)
0752
0753       }  //       if (recojet)
0754     }  // unique match
0755
0756   }  //  if (truthjet)
0757
0758   // next for reco jets
0759   JetContainer* recojets = findNode::getClass<JetContainer>(topNode, reco_jet_name);
0760   if (!recojets)
0761   {
0762     std::cout
0763         << "QAG4SimulationJet::process_TruthMatching - Error can not find DST JetContainer node "
0764         << reco_jet_name << std::endl;
0765     exit(-1);
0766   }
0767
0768   // search for leading reco jet
0769   Jet* recojet = nullptr;
0770   max_et = 0;
0771   for (auto jet : *recojets)
0772   {
0773     assert(jet);
0774
0775     if (not jet_acceptance_cut(jet))
0776     {
0777       continue;
0778     }
0779
0780     if (jet->get_et() > max_et)
0781     {
0782       recojet = jet;
0783       max_et = jet->get_et();
0784     }
0785   }
0786
0787   // match leading reco jet
0788   if (recojet)
0789   {
0790     if (Verbosity() > 1)
0791     {
0792       std::cout << "QAG4SimulationJet::process_TruthMatching - " << reco_jet_name
0793                 << " process reco jet ";
0794       recojet->identify();
0795     }
0796
0797     Matching_Count_Reco_Et->Fill(recojet->get_et(), "Total", 1);
0798
0799     {  // inclusive best energy match
0800       Jet* truthjet1 = recoeval->max_truth_jet_by_energy(recojet);
0801       if (truthjet1)
0802       {
0803         const double dPhi = recojet->get_phi() - truthjet1->get_phi();
0804         if (fabs(dPhi) < _jet_match_dPhi)
0805         {
0806           const double dEta = recojet->get_eta() - truthjet1->get_eta();
0807           if (fabs(dEta) < _jet_match_dEta)
0808           {
0809             const double E_r = recojet->get_e() / (truthjet1->get_e() + 1e-9);
0810
0811             if (fabs(E_r - 1) < _jet_match_dE_Ratio)
0812             {
0813               // matched in eta, phi and energy
0814
0815               Matching_Count_Reco_Et->Fill(recojet->get_et(),
0816                                            "Unique Matched", 1);
0817             }
0818
0819           }  //  if (fabs(dEta) < 0.1)
0820
0821         }  // if (fabs(dPhi) < 0.1)
0822
0823       }  //      if (truthjet1)
0824     }  // inclusive best energy match
0825
0826     {  // unique match
0827       Jet* truthjet2 = recoeval->unique_truth_jet_from_reco(recojet);
0828       if (truthjet2)
0829       {
0830         const double dPhi = recojet->get_phi() - truthjet2->get_phi();
0831         if (fabs(dPhi) < _jet_match_dPhi)
0832         {
0833           const double dEta = recojet->get_eta() - truthjet2->get_eta();
0834           if (fabs(dEta) < _jet_match_dEta)
0835           {
0836             const double E_r = recojet->get_e() / (truthjet2->get_e() + 1e-9);
0837
0838             if (fabs(E_r - 1) < _jet_match_dE_Ratio)
0839             {
0840               // matched in eta, phi and energy
0841
0842               Matching_Count_Reco_Et->Fill(recojet->get_et(),
0843                                            "Matched", 1);
0844             }
0845
0846           }  //  if (fabs(dEta) < 0.1)
0847
0848         }  // if (fabs(dPhi) < 0.1)
0849
0850       }  //      if (truthjet2)
0851     }  // unique match
0852
0853   }  // if (recojet)
0854
0855   return Fun4AllReturnCodes::EVENT_OK;
0856 }