File indexing completed on 2025-12-18 09:15:12
0001
0002
0003 #include <fun4all/Fun4AllBase.h>
0004 #include <JetValidation.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <fun4all/PHTFileServer.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009 #include <jetbase/JetMap.h>
0010 #include <jetbase/JetContainer.h>
0011 #include <jetbase/Jetv2.h>
0012 #include <jetbase/Jetv1.h>
0013 #include <centrality/CentralityInfo.h>
0014 #include <globalvertex/GlobalVertex.h>
0015 #include <globalvertex/GlobalVertexMap.h>
0016 #include <calobase/RawTower.h>
0017 #include <calobase/RawTowerContainer.h>
0018 #include <calobase/RawTowerGeom.h>
0019 #include <calobase/RawTowerGeomContainer.h>
0020 #include <calobase/TowerInfoContainer.h>
0021 #include <calobase/TowerInfo.h>
0022 #include <ffarawobjects/Gl1Packet.h>
0023 #include <jetbackground/TowerBackground.h>
0024 #include <calobase/RawCluster.h>
0025 #include <calobase/RawClusterContainer.h>
0026 #include <calobase/RawClusterUtility.h>
0027 #include <calobase/RawTowerDefs.h>
0028
0029 #include <TTree.h>
0030 #include <iostream>
0031 #include <sstream>
0032 #include <iomanip>
0033 #include <cmath>
0034 #include <vector>
0035
0036 #include "TH1.h"
0037 #include "TH2.h"
0038 #include "TFile.h"
0039 #include <TTree.h>
0040
0041
0042 JetValidation::JetValidation(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
0043 SubsysReco("JetValidation_" + recojetname + "_" + truthjetname)
0044 , m_recoJetName(recojetname)
0045 , m_truthJetName(truthjetname)
0046 , m_outputFileName(outputfilename)
0047 , m_etaRange(-1, 1)
0048 , m_ptRange(5, 100)
0049 , m_doTruthJets(0)
0050 , m_doSeeds(0)
0051 , m_doUnsubJet(0)
0052 , m_T(nullptr)
0053 , m_event(-1)
0054 , m_nTruthJet(-1)
0055 , m_nJet(-1)
0056 , m_id()
0057 , m_nComponent()
0058 , m_eta()
0059 , m_phi()
0060 , m_e()
0061 , m_pt()
0062 , m_cleta()
0063 , m_clphi()
0064 , m_cle()
0065 , m_clecore()
0066 , m_clpt()
0067 , m_clprob()
0068 , m_sub_et()
0069 , m_truthID()
0070 , m_truthNComponent()
0071 , m_truthEta()
0072 , m_truthPhi()
0073 , m_truthE()
0074 , m_truthPt()
0075 , m_eta_rawseed()
0076 , m_phi_rawseed()
0077 , m_pt_rawseed()
0078 , m_e_rawseed()
0079 , m_rawseed_cut()
0080 , m_eta_subseed()
0081 , m_phi_subseed()
0082 , m_pt_subseed()
0083 , m_e_subseed()
0084 , m_subseed_cut()
0085 {
0086 std::cout << "JetValidation::JetValidation(const std::string &name) Calling ctor" << std::endl;
0087 }
0088
0089
0090 JetValidation::~JetValidation()
0091 {
0092 std::cout << "JetValidation::~JetValidation() Calling dtor" << std::endl;
0093 }
0094
0095
0096
0097 int JetValidation::Init(PHCompositeNode *topNode)
0098 {
0099 std::cout << "JetValidation::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0100 PHTFileServer::get().open(m_outputFileName, "RECREATE");
0101 std::cout << "JetValidation::Init - Output to " << m_outputFileName << std::endl;
0102
0103
0104 m_T = new TTree("T", "MyJetAnalysis Tree");
0105 m_T->Branch("m_event", &m_event, "event/I");
0106 m_T->Branch("nJet", &m_nJet, "nJet/I");
0107 m_T->Branch("cent", &m_centrality);
0108 m_T->Branch("zvtx", &m_zvtx);
0109 m_T->Branch("b", &m_impactparam);
0110 m_T->Branch("id", &m_id);
0111 m_T->Branch("nComponent", &m_nComponent);
0112 m_T->Branch("triggerVector", &m_triggerVector);
0113
0114 m_T->Branch("eta", &m_eta);
0115 m_T->Branch("phi", &m_phi);
0116 m_T->Branch("e", &m_e);
0117 m_T->Branch("pt", &m_pt);
0118
0119 m_T->Branch("cleta", &m_cleta);
0120 m_T->Branch("clphi", &m_clphi);
0121 m_T->Branch("cle", &m_cle);
0122 m_T->Branch("clecore", &m_clecore);
0123 m_T->Branch("clpt", &m_clpt);
0124 m_T->Branch("clprob", &m_clprob);
0125
0126 if(m_doUnsubJet)
0127 {
0128 m_T->Branch("pt_unsub", &m_unsub_pt);
0129 m_T->Branch("subtracted_et", &m_sub_et);
0130 }
0131 if(m_doTruthJets){
0132 m_T->Branch("nTruthJet", &m_nTruthJet);
0133 m_T->Branch("truthID", &m_truthID);
0134 m_T->Branch("truthNComponent", &m_truthNComponent);
0135 m_T->Branch("truthEta", &m_truthEta);
0136 m_T->Branch("truthPhi", &m_truthPhi);
0137 m_T->Branch("truthE", &m_truthE);
0138 m_T->Branch("truthPt", &m_truthPt);
0139 }
0140
0141 if(m_doSeeds){
0142 m_T->Branch("rawseedEta", &m_eta_rawseed);
0143 m_T->Branch("rawseedPhi", &m_phi_rawseed);
0144 m_T->Branch("rawseedPt", &m_pt_rawseed);
0145 m_T->Branch("rawseedE", &m_e_rawseed);
0146 m_T->Branch("rawseedCut", &m_rawseed_cut);
0147 m_T->Branch("subseedEta", &m_eta_subseed);
0148 m_T->Branch("subseedPhi", &m_phi_subseed);
0149 m_T->Branch("subseedPt", &m_pt_subseed);
0150 m_T->Branch("subseedE", &m_e_subseed);
0151 m_T->Branch("subseedCut", &m_subseed_cut);
0152 }
0153
0154 return Fun4AllReturnCodes::EVENT_OK;
0155 }
0156
0157
0158 int JetValidation::InitRun(PHCompositeNode *topNode)
0159 {
0160 std::cout << "JetValidation::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0161
0162 return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165
0166 int JetValidation::process_event(PHCompositeNode *topNode)
0167 {
0168
0169 ++m_event;
0170
0171
0172 JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0173 if (!jets)
0174 {
0175 std::cout
0176 << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0177 << m_recoJetName << std::endl;
0178 exit(-1);
0179 }
0180
0181
0182
0183 JetContainer* jetsMC = findNode::getClass<JetContainer>(topNode, m_truthJetName);
0184 if (!jetsMC && m_doTruthJets)
0185 {
0186 std::cout
0187 << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0188 << m_truthJetName << std::endl;
0189 exit(-1);
0190 }
0191
0192
0193 JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0194 if (!seedjetsraw && m_doSeeds)
0195 {
0196 std::cout
0197 << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
0198 << std::endl;
0199 exit(-1);
0200 }
0201
0202 JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0203 if (!seedjetssub && m_doSeeds)
0204 {
0205 std::cout
0206 << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
0207 << std::endl;
0208 exit(-1);
0209 }
0210
0211
0212 CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0213 if (!cent_node)
0214 {
0215 std::cout
0216 << "MyJetAnalysis::process_event - Error can not find centrality node "
0217 << std::endl;
0218 exit(-1);
0219 }
0220
0221
0222 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0223 if (!vertexmap)
0224 {
0225 std::cout
0226 << "MyJetAnalysis::process_event - Error can not find global vertex node "
0227 << std::endl;
0228 exit(-1);
0229 }
0230 if (vertexmap->empty())
0231 {
0232 std::cout
0233 << "MyJetAnalysis::process_event - global vertex node is empty "
0234 << std::endl;
0235 return Fun4AllReturnCodes::ABORTEVENT;
0236 }
0237 else
0238 {
0239 GlobalVertex *vtx = vertexmap->begin()->second;
0240 m_zvtx = vtx->get_z();
0241 }
0242
0243
0244 TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0245 TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0246 TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0247 RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0248 RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0249 if(!towersEM3 || !towersIH3 || !towersOH3){
0250 std::cout
0251 <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0252 << std::endl;
0253 exit(-1);
0254 }
0255
0256 if(!tower_geom || !tower_geomOH){
0257 std::cout
0258 <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0259 << std::endl;
0260 exit(-1);
0261 }
0262
0263
0264 TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0265 if(m_doUnsubJet) {
0266 if(!background){
0267 std::cout<<"Can't get background. Exiting"<<std::endl;
0268 return Fun4AllReturnCodes::EVENT_OK;
0269 }
0270 }
0271
0272
0273
0274 m_centrality = (int)(100*cent_node->get_centile(CentralityInfo::PROP::mbd_NS));
0275 m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
0276
0277
0278 m_nJet = 0;
0279 float background_v2 = 0;
0280 float background_Psi2 = 0;
0281 if(m_doUnsubJet)
0282 {
0283 background_v2 = background->get_v2();
0284 background_Psi2 = background->get_Psi2();
0285 }
0286
0287 for (auto jet : *jets)
0288 {
0289
0290 if(jet->get_pt() < 1) continue;
0291
0292 m_id.push_back(jet->get_id());
0293 m_nComponent.push_back(jet->size_comp());
0294 m_eta.push_back(jet->get_eta());
0295 m_phi.push_back(jet->get_phi());
0296 m_e.push_back(jet->get_e());
0297 m_pt.push_back(jet->get_pt());
0298
0299 if(m_doUnsubJet)
0300 {
0301 Jet* unsubjet = new Jetv1();
0302
0303 float totalPx = 0;
0304 float totalPy = 0;
0305 float totalPz = 0;
0306 float totalE = 0;
0307 int nconst = 0;
0308
0309 for (auto comp: jet->get_comp_vec())
0310 {
0311 TowerInfo *tower;
0312 nconst++;
0313 unsigned int channel = comp.second;
0314
0315 if (comp.first == 15 || comp.first == 30)
0316 {
0317 tower = towersIH3->get_tower_at_channel(channel);
0318 if(!tower || !tower_geom){
0319 continue;
0320 }
0321 unsigned int calokey = towersIH3->encode_key(channel);
0322 int ieta = towersIH3->getTowerEtaBin(calokey);
0323 int iphi = towersIH3->getTowerPhiBin(calokey);
0324 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0325 float UE = background->get_UE(1).at(ieta);
0326 float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0327 float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0328
0329 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0330 totalE += tower->get_energy() + UE;
0331 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0332 totalPx += pt * cos(tower_phi);
0333 totalPy += pt * sin(tower_phi);
0334 totalPz += pt * sinh(tower_eta);
0335 }
0336 else if (comp.first == 16 || comp.first == 31)
0337 {
0338 tower = towersOH3->get_tower_at_channel(channel);
0339 if(!tower || !tower_geomOH)
0340 {
0341 continue;
0342 }
0343
0344 unsigned int calokey = towersOH3->encode_key(channel);
0345 int ieta = towersOH3->getTowerEtaBin(calokey);
0346 int iphi = towersOH3->getTowerPhiBin(calokey);
0347 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0348 float UE = background->get_UE(2).at(ieta);
0349 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0350 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0351
0352 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0353 totalE +=tower->get_energy() + UE;
0354 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0355 totalPx += pt * cos(tower_phi);
0356 totalPy += pt * sin(tower_phi);
0357 totalPz += pt * sinh(tower_eta);
0358 }
0359 else if (comp.first == 14 || comp.first == 29)
0360 {
0361 tower = towersEM3->get_tower_at_channel(channel);
0362 if(!tower || !tower_geom)
0363 {
0364 continue;
0365 }
0366
0367 unsigned int calokey = towersEM3->encode_key(channel);
0368 int ieta = towersEM3->getTowerEtaBin(calokey);
0369 int iphi = towersEM3->getTowerPhiBin(calokey);
0370 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0371 float UE = background->get_UE(0).at(ieta);
0372 float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0373 float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0374
0375 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0376 totalE +=tower->get_energy() + UE;
0377 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0378 totalPx += pt * cos(tower_phi);
0379 totalPy += pt * sin(tower_phi);
0380 totalPz += pt * sinh(tower_eta);
0381
0382 }
0383 }
0384
0385 unsubjet->set_px(totalPx);
0386 unsubjet->set_py(totalPy);
0387 unsubjet->set_pz(totalPz);
0388 unsubjet->set_e(totalE);
0389 m_unsub_pt.push_back(unsubjet->get_pt());
0390 m_sub_et.push_back(unsubjet->get_et() - jet->get_et());
0391 }
0392
0393
0394 m_nJet++;
0395 }
0396
0397
0398 if(m_doTruthJets)
0399 {
0400 m_nTruthJet = 0;
0401
0402 for (auto truthjet : *jetsMC)
0403 {
0404
0405
0406 bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0407 bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0408 if ((not eta_cut) or (not pt_cut)) continue;
0409 m_truthID.push_back(truthjet->get_id());
0410 m_truthNComponent.push_back(truthjet->size_comp());
0411 m_truthEta.push_back(truthjet->get_eta());
0412 m_truthPhi.push_back(truthjet->get_phi());
0413 m_truthE.push_back(truthjet->get_e());
0414 m_truthPt.push_back(truthjet->get_pt());
0415 m_nTruthJet++;
0416 }
0417 }
0418
0419
0420 if(m_doSeeds)
0421 {
0422 for (auto jet : *seedjetsraw)
0423 {
0424 int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0425 m_eta_rawseed.push_back(jet->get_eta());
0426 m_phi_rawseed.push_back(jet->get_phi());
0427 m_e_rawseed.push_back(jet->get_e());
0428 m_pt_rawseed.push_back(jet->get_pt());
0429 m_rawseed_cut.push_back(passesCut);
0430 }
0431
0432 for (auto jet : *seedjetssub)
0433 {
0434 int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0435 m_eta_subseed.push_back(jet->get_eta());
0436 m_phi_subseed.push_back(jet->get_phi());
0437 m_e_subseed.push_back(jet->get_e());
0438 m_pt_subseed.push_back(jet->get_pt());
0439 m_subseed_cut.push_back(passesCut);
0440 }
0441 }
0442
0443
0444 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "14001");
0445 if (!gl1PacketInfo)
0446 {
0447 std::cout << PHWHERE << "caloTreeGen::process_event: GL1Packet node is missing. Output related to this node will be empty" << std::endl;
0448 }
0449
0450 if (gl1PacketInfo)
0451 {
0452 uint64_t triggervec = gl1PacketInfo->getTriggerVector();
0453 for (int i = 0; i < 64; i++)
0454 {
0455 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0456 m_triggerVector.push_back(trig_decision);
0457 triggervec = (triggervec >> 1U) & 0xffffffffU;
0458 }
0459 }
0460
0461
0462 RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0463 if (!clusterContainer)
0464 {
0465 std::cout << PHWHERE << "caloTreeGen::process_event: CLUSTERINFO_CEMC node is missing. Output related to this node will be empty" << std::endl;
0466 return 0;
0467 }
0468 RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0469 RawClusterContainer::ConstIterator clusterIter;
0470 for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0471 {
0472 RawCluster *recoCluster = clusterIter->second;
0473
0474 CLHEP::Hep3Vector vertex(0, 0, 0);
0475 if (m_zvtx != -9999)
0476 {
0477 vertex.setZ(m_zvtx);
0478 }
0479 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0480 CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, vertex);
0481
0482 float clusE = E_vec_cluster_Full.mag();
0483 float clusEcore = E_vec_cluster.mag();
0484 float clus_eta = E_vec_cluster.pseudoRapidity();
0485 float clus_phi = E_vec_cluster.phi();
0486 float clus_pt = E_vec_cluster.perp();
0487 float clus_prob = recoCluster->get_prob();
0488
0489 m_cle.push_back(clusE);
0490 m_clecore.push_back(clusEcore);
0491 m_cleta.push_back(clus_eta);
0492 m_clphi.push_back(clus_phi);
0493 m_clpt.push_back(clus_pt);
0494 m_clprob.push_back(clus_prob);
0495 }
0496
0497
0498 m_T->Fill();
0499
0500 return Fun4AllReturnCodes::EVENT_OK;
0501 }
0502
0503
0504 int JetValidation::ResetEvent(PHCompositeNode *topNode)
0505 {
0506
0507 m_id.clear();
0508 m_nComponent.clear();
0509 m_eta.clear();
0510 m_phi.clear();
0511 m_e.clear();
0512 m_pt.clear();
0513 m_unsub_pt.clear();
0514 m_sub_et.clear();
0515
0516 m_truthID.clear();
0517 m_truthNComponent.clear();
0518 m_truthEta.clear();
0519 m_truthPhi.clear();
0520 m_truthE.clear();
0521 m_truthPt.clear();
0522 m_truthdR.clear();
0523
0524 m_eta_subseed.clear();
0525 m_phi_subseed.clear();
0526 m_e_subseed.clear();
0527 m_pt_subseed.clear();
0528 m_subseed_cut.clear();
0529
0530 m_eta_rawseed.clear();
0531 m_phi_rawseed.clear();
0532 m_e_rawseed.clear();
0533 m_pt_rawseed.clear();
0534 m_rawseed_cut.clear();
0535
0536 m_triggerVector.clear();
0537
0538 m_cle.clear();
0539 m_clecore.clear();
0540 m_cleta.clear();
0541 m_clphi.clear();
0542 m_clpt.clear();
0543 m_clprob.clear();
0544 return Fun4AllReturnCodes::EVENT_OK;
0545 }
0546
0547
0548 int JetValidation::EndRun(const int runnumber)
0549 {
0550 std::cout << "JetValidation::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0551 return Fun4AllReturnCodes::EVENT_OK;
0552 }
0553
0554
0555 int JetValidation::End(PHCompositeNode *topNode)
0556 {
0557 std::cout << "JetValidation::End - Output to " << m_outputFileName << std::endl;
0558 PHTFileServer::get().cd(m_outputFileName);
0559
0560 m_T->Write();
0561 std::cout << "JetValidation::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0562 return Fun4AllReturnCodes::EVENT_OK;
0563 }
0564
0565
0566 int JetValidation::Reset(PHCompositeNode *topNode)
0567 {
0568 std::cout << "JetValidation::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0569 return Fun4AllReturnCodes::EVENT_OK;
0570 }
0571
0572
0573 void JetValidation::Print(const std::string &what) const
0574 {
0575 std::cout << "JetValidation::Print(const std::string &what) const Printing info for " << what << std::endl;
0576 }
