File indexing completed on 2025-08-05 08:11:14
0001 #include "JetCalibration.h"
0002 #include "FindDijets.h"
0003
0004 #include <TH1.h>
0005 #include <TH2.h>
0006 #include <TH2D.h>
0007 #include <TH1D.h>
0008 #include <TEfficiency.h>
0009 #include <g4main/PHG4HitContainer.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011 #include <g4main/PHG4Particle.h>
0012 #include <g4main/PHG4Hit.h>
0013 #include <g4main/PHG4Shower.h>
0014 #include <g4main/PHG4HitDefs.h>
0015 #include <TMath.h>
0016 #include <jetbase/Jet.h>
0017 #include <jetbase/JetContainer.h>
0018 #include <jetbase/JetContainerv1.h>
0019 #include <jetbase/Jetv2.h>
0020 #include <fun4all/Fun4AllHistoManager.h>
0021 #include <HepMC/SimpleVector.h>
0022
0023 #include <vector>
0024
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/PHIODataNode.h>
0028 #include <phool/PHNode.h>
0029 #include <phool/PHNodeIterator.h>
0030 #include <phool/PHObject.h>
0031 #include <phool/getClass.h>
0032
0033
0034 #include <iostream>
0035
0036 #include <map>
0037
0038
0039 JetCalibration::JetCalibration(const std::string &name, const std::string &outfilename):
0040 SubsysReco(name)
0041
0042 {
0043 _foutname = outfilename;
0044
0045 }
0046
0047
0048 JetCalibration::~JetCalibration()
0049 {
0050
0051 }
0052
0053
0054 int JetCalibration::Init(PHCompositeNode *topNode)
0055 {
0056 _i_event = 0;
0057 hm = new Fun4AllHistoManager("TRIGGER_HISTOGRAMS");
0058
0059 TH1D *h;
0060 TH2D *h2;
0061
0062 h = new TH1D("h_truth_reco", "", 200, 0, 2);
0063 hm->registerHisto(h);
0064
0065 h = new TH1D("h_match_dR", "", 100, 0, 0.2);
0066 hm->registerHisto(h);
0067
0068 h = new TH1D("h_truth_iso", "", 100, 0, TMath::Pi());
0069 hm->registerHisto(h);
0070
0071 TEfficiency *hp = new TEfficiency("h_match_eff", "", 80, 0, 40);
0072 hm->registerHisto(hp);
0073
0074 h = new TH1D("h_truth_em", "", 100, 0, 1);
0075 hm->registerHisto(h);
0076
0077 h2 = new TH2D("h_truth_em_v_truth", "", 60, 0, 60, 100, 0, 1);
0078 hm->registerHisto(h2);
0079
0080 h2 = new TH2D("h_truth_reco_v_truth", "", 60, 0, 60, 200, 0, 2);
0081 hm->registerHisto(h2);
0082
0083 h2 = new TH2D("h_truth_reco_v_truth_em", "", 100, 0, 1, 200, 0, 2);
0084 hm->registerHisto(h2);
0085
0086 h2 = new TH2D("h_truth_reco_v_truth_em_shower", "", 100, 0, 1, 200, 0, 2);
0087 hm->registerHisto(h2);
0088
0089 return 0;
0090 }
0091
0092
0093 int JetCalibration::InitRun(PHCompositeNode *topNode)
0094 {
0095 return Fun4AllReturnCodes::EVENT_OK;
0096 }
0097
0098 int JetCalibration::process_event(PHCompositeNode *topNode)
0099 {
0100
0101 _i_event++;
0102
0103 std::string recoJetName = "AntiKt_Tower_r04_Sub1";
0104
0105 JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0106
0107 recoJetName = "AntiKt_Truth_r04";
0108
0109 JetContainer *jets_truth_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0110 PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0111 PHG4TruthInfoContainer::ShowerRange range = truthinfo->GetPrimaryShowerRange();
0112
0113
0114 PHG4HitContainer *hitinfo;
0115 hitinfo = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_CEMC");
0116 if (!hitinfo)
0117 {
0118 return Fun4AllReturnCodes::ABORTRUN;
0119 }
0120
0121 m_HitContainerMap[hitinfo->GetID()] = hitinfo;
0122
0123 hitinfo = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_HCALOUT");
0124 m_HitContainerMap[hitinfo->GetID()] = hitinfo;
0125
0126 hitinfo = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_HCALIN");
0127 m_HitContainerMap[hitinfo->GetID()] = hitinfo;
0128
0129 std::vector<struct jetty> truth_jets;
0130 std::vector<struct jetty> reco_jets;
0131 for (auto jet : *jets_truth_4)
0132 {
0133 if (jet->get_pt() < 4) continue;
0134 if (fabs(jet->get_eta()) > .7) continue;
0135 struct jetty myjet;
0136 myjet.pt = jet->get_pt();
0137 myjet.et = 0;
0138 myjet.eta = jet->get_eta();
0139 myjet.phi = jet->get_phi();
0140 myjet.emcal = 0;
0141 myjet.ohcal = 0;
0142 float jet_em = 0;
0143 float jet_h = 0;
0144 float jet_edep = 0;
0145 for (const auto& comp : jet->get_comp_vec())
0146 {
0147 unsigned int index = comp.second;
0148 PHG4Particle* truth_particle = truthinfo->GetParticle(index);
0149
0150 assert(truth_particle);
0151
0152 PHG4Shower* shower = nullptr;
0153 for (PHG4TruthInfoContainer::ShowerIterator iter = range.first;
0154 iter != range.second;
0155 ++iter)
0156 {
0157 shower = iter->second;
0158 if (shower->get_parent_particle_id() == truth_particle->get_track_id())
0159 {
0160 break;
0161 }
0162 }
0163
0164 if (!shower)
0165 {
0166 continue;
0167 }
0168
0169
0170 float edep = 0.0;
0171 float edep_e = 0.0;
0172 float edep_h = 0.0;
0173
0174
0175 for (std::map<int, std::set<PHG4HitDefs::keytype> >::iterator iter = shower->begin_g4hit_id();
0176 iter != shower->end_g4hit_id();
0177 ++iter)
0178 {
0179 int g4hitmap_id = iter->first;
0180 std::map<int, PHG4HitContainer*>::iterator mapiter = m_HitContainerMap.find(g4hitmap_id);
0181 if (mapiter == m_HitContainerMap.end())
0182 {
0183 continue;
0184 }
0185
0186 PHG4HitContainer* hits = mapiter->second;
0187
0188
0189 for (unsigned long long g4hit_id : iter->second)
0190 {
0191 PHG4Hit* g4hit = hits->findHit(g4hit_id);
0192 if (!g4hit)
0193 {
0194 continue;
0195 }
0196
0197 PHG4Particle* particle = truthinfo->GetParticle(g4hit->get_trkid());
0198 if (!particle)
0199 {
0200 continue;
0201 }
0202
0203 if (!isnan(g4hit->get_edep()))
0204 {
0205 edep += g4hit->get_edep();
0206 if (abs(particle->get_pid()) == 11)
0207 {
0208 edep_e += g4hit->get_edep();
0209 }
0210 else
0211 {
0212 edep_h += g4hit->get_edep();
0213 }
0214 }
0215
0216 }
0217
0218 }
0219 int pid = truth_particle->get_pid();
0220 double energy = sqrt(TMath::Power(truth_particle->get_px(), 2) + TMath::Power(truth_particle->get_py(), 2));
0221
0222 if (isEM(pid))
0223 {
0224 myjet.emcal += energy;
0225 }
0226 else
0227 {
0228 myjet.ohcal += energy;
0229 }
0230 myjet.et += energy;
0231
0232 jet_em += edep_e;
0233 jet_h += edep_h;
0234 jet_edep += edep;
0235 }
0236 if (jet_edep == 0.0) continue;
0237 myjet.edep = jet_edep;
0238 myjet.edep_em = jet_em;
0239 myjet.edep_hd = jet_h;
0240 truth_jets.push_back(myjet);
0241
0242 }
0243 int size = truth_jets.size();
0244 for (int ijet = 0 ; ijet < size ; ijet++)
0245 {
0246 struct jetty jet1 = truth_jets.at(ijet);
0247 float min_dR = 99;
0248 for (int jjet = 0; jjet < size; jjet++)
0249 {
0250 if (jjet == ijet) continue;
0251 struct jetty jet2 = truth_jets.at(jjet);
0252 double dR = getDr(jet1, jet2);
0253 if (dR < min_dR)
0254 {
0255 min_dR = dR;
0256 }
0257 }
0258 truth_jets.at(ijet).iso = min_dR;
0259 }
0260 auto h_truth_reco = dynamic_cast<TH1*>(hm->getHisto("h_truth_reco"));
0261 auto h_truth_em = dynamic_cast<TH1*>(hm->getHisto("h_truth_em"));
0262 auto h_match_dR = dynamic_cast<TH1*>(hm->getHisto("h_match_dR"));
0263 auto h_match_eff = dynamic_cast<TEfficiency*>(hm->getHisto("h_match_eff"));
0264 auto h_truth_em_v_truth = dynamic_cast<TH2*>(hm->getHisto("h_truth_em_v_truth"));
0265 auto h_truth_reco_v_truth_em = dynamic_cast<TH2*>(hm->getHisto("h_truth_reco_v_truth_em"));
0266 auto h_truth_reco_v_truth_em_shower = dynamic_cast<TH2*>(hm->getHisto("h_truth_reco_v_truth_em_shower"));
0267 auto h_truth_reco_v_truth = dynamic_cast<TH2*>(hm->getHisto("h_truth_reco_v_truth"));
0268 auto h_truth_iso = dynamic_cast<TH1*>(hm->getHisto("h_truth_iso"));
0269 for (auto jet : *jets_4)
0270 {
0271
0272 struct jetty myjet;
0273 if (jet->get_pt() < 1) continue;
0274 myjet.pt = jet->get_pt();
0275 myjet.eta = jet->get_eta();
0276 myjet.phi = jet->get_phi();
0277
0278 reco_jets.push_back(myjet);
0279 }
0280
0281 std::sort(reco_jets.begin(), reco_jets.end(), [] (auto a, auto b) { return a.pt > b.pt; });
0282 std::sort(truth_jets.begin(), truth_jets.end(), [] (auto a, auto b) { return a.pt > b.pt; });
0283
0284 for (auto &jet : truth_jets)
0285 {
0286 h_truth_iso->Fill(jet.iso);
0287 if (jet.iso < 1.0) continue;
0288 bool matched = false;
0289 for (auto &jetreco : reco_jets)
0290 {
0291 double dR = getDr(jet, jetreco);
0292 if (dR < 0.2)
0293 {
0294 matched = true;
0295 h_truth_reco->Fill(jetreco.pt/jet.pt);
0296 h_truth_em->Fill(jet.emcal/jet.pt);
0297 h_match_dR->Fill(dR);
0298 h_truth_em_v_truth->Fill(jet.pt, jet.emcal/jet.pt);
0299 if (jet.pt > 20 && jet.pt < 30) h_truth_reco_v_truth_em->Fill(jet.emcal/jet.pt, jetreco.pt/jet.pt);
0300 if (jet.pt > 20 && jet.pt < 30) h_truth_reco_v_truth_em_shower->Fill(jet.edep_em/jet.edep, jetreco.pt/jet.pt);
0301 h_truth_reco_v_truth->Fill(jet.pt, jetreco.pt/jet.pt);
0302 break;
0303 }
0304 }
0305 h_match_eff->Fill(matched, jet.pt);
0306 }
0307
0308
0309 return Fun4AllReturnCodes::EVENT_OK;
0310 }
0311
0312 double JetCalibration::getDr(struct jetty jet1, struct jetty jet2)
0313 {
0314 float dphi = jet1.phi - jet2.phi;
0315
0316 if (dphi > TMath::Pi())
0317 {
0318 dphi = dphi - 2.*TMath::Pi();
0319 }
0320 if (dphi <= -1*TMath::Pi())
0321 {
0322 dphi = dphi + 2. * TMath::Pi();
0323 }
0324
0325 return sqrt(TMath::Power(dphi, 2) + TMath::Power(jet1.eta - jet2.eta, 2));
0326 }
0327 bool JetCalibration::isEM(int pid)
0328 {
0329 if (pid == 11 || pid == -11 || pid == 22 || pid == 111 || pid == 221)
0330 {
0331 return true;
0332 }
0333 return false;
0334 }
0335
0336 void JetCalibration::GetNodes(PHCompositeNode* topNode)
0337 {
0338
0339
0340 }
0341
0342 int JetCalibration::ResetEvent(PHCompositeNode *topNode)
0343 {
0344
0345 return Fun4AllReturnCodes::EVENT_OK;
0346 }
0347
0348
0349 int JetCalibration::EndRun(const int runnumber)
0350 {
0351 return Fun4AllReturnCodes::EVENT_OK;
0352 }
0353
0354
0355 int JetCalibration::End(PHCompositeNode *topNode)
0356 {
0357 if (Verbosity() > 0)
0358 {
0359 std::cout << "JetCalibration::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0360 }
0361 std::cout<<"Total events: "<<_i_event<<std::endl;
0362 hm->dumpHistos(_foutname.c_str());
0363 return Fun4AllReturnCodes::EVENT_OK;
0364 }
0365
0366
0367 int JetCalibration::Reset(PHCompositeNode *topNode)
0368 {
0369 if (Verbosity() > 0)
0370 {
0371 std::cout << "JetCalibration::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0372 }
0373 return Fun4AllReturnCodes::EVENT_OK;
0374 }
