QA/SimulationModules/QAG4SimulationJet.cc

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