File indexing completed on 2025-08-05 08:13:20
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(!background){
0266 std::cout<<"Can't get background. Exiting"<<std::endl;
0267 return Fun4AllReturnCodes::EVENT_OK;
0268 }
0269
0270
0271
0272 m_centrality = (int)(100*cent_node->get_centile(CentralityInfo::PROP::mbd_NS));
0273 m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
0274
0275
0276 m_nJet = 0;
0277 float background_v2 = 0;
0278 float background_Psi2 = 0;
0279 if(m_doUnsubJet)
0280 {
0281 background_v2 = background->get_v2();
0282 background_Psi2 = background->get_Psi2();
0283 }
0284
0285 for (auto jet : *jets)
0286 {
0287
0288 if(jet->get_pt() < 1) continue;
0289
0290 m_id.push_back(jet->get_id());
0291 m_nComponent.push_back(jet->size_comp());
0292 m_eta.push_back(jet->get_eta());
0293 m_phi.push_back(jet->get_phi());
0294 m_e.push_back(jet->get_e());
0295 m_pt.push_back(jet->get_pt());
0296
0297 if(m_doUnsubJet)
0298 {
0299 Jet* unsubjet = new Jetv1();
0300
0301 float totalPx = 0;
0302 float totalPy = 0;
0303 float totalPz = 0;
0304 float totalE = 0;
0305 int nconst = 0;
0306
0307 for (auto comp: jet->get_comp_vec())
0308 {
0309 TowerInfo *tower;
0310 nconst++;
0311 unsigned int channel = comp.second;
0312
0313 if (comp.first == 15 || comp.first == 30)
0314 {
0315 tower = towersIH3->get_tower_at_channel(channel);
0316 if(!tower || !tower_geom){
0317 continue;
0318 }
0319 unsigned int calokey = towersIH3->encode_key(channel);
0320 int ieta = towersIH3->getTowerEtaBin(calokey);
0321 int iphi = towersIH3->getTowerPhiBin(calokey);
0322 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0323 float UE = background->get_UE(1).at(ieta);
0324 float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0325 float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0326
0327 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0328 totalE += tower->get_energy() + UE;
0329 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0330 totalPx += pt * cos(tower_phi);
0331 totalPy += pt * sin(tower_phi);
0332 totalPz += pt * sinh(tower_eta);
0333 }
0334 else if (comp.first == 16 || comp.first == 31)
0335 {
0336 tower = towersOH3->get_tower_at_channel(channel);
0337 if(!tower || !tower_geomOH)
0338 {
0339 continue;
0340 }
0341
0342 unsigned int calokey = towersOH3->encode_key(channel);
0343 int ieta = towersOH3->getTowerEtaBin(calokey);
0344 int iphi = towersOH3->getTowerPhiBin(calokey);
0345 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0346 float UE = background->get_UE(2).at(ieta);
0347 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0348 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0349
0350 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0351 totalE +=tower->get_energy() + UE;
0352 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0353 totalPx += pt * cos(tower_phi);
0354 totalPy += pt * sin(tower_phi);
0355 totalPz += pt * sinh(tower_eta);
0356 }
0357 else if (comp.first == 14 || comp.first == 29)
0358 {
0359 tower = towersEM3->get_tower_at_channel(channel);
0360 if(!tower || !tower_geom)
0361 {
0362 continue;
0363 }
0364
0365 unsigned int calokey = towersEM3->encode_key(channel);
0366 int ieta = towersEM3->getTowerEtaBin(calokey);
0367 int iphi = towersEM3->getTowerPhiBin(calokey);
0368 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0369 float UE = background->get_UE(0).at(ieta);
0370 float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0371 float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0372
0373 UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0374 totalE +=tower->get_energy() + UE;
0375 double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0376 totalPx += pt * cos(tower_phi);
0377 totalPy += pt * sin(tower_phi);
0378 totalPz += pt * sinh(tower_eta);
0379
0380 }
0381 }
0382
0383 unsubjet->set_px(totalPx);
0384 unsubjet->set_py(totalPy);
0385 unsubjet->set_pz(totalPz);
0386 unsubjet->set_e(totalE);
0387 m_unsub_pt.push_back(unsubjet->get_pt());
0388 m_sub_et.push_back(unsubjet->get_et() - jet->get_et());
0389 }
0390
0391
0392 m_nJet++;
0393 }
0394
0395
0396 if(m_doTruthJets)
0397 {
0398 m_nTruthJet = 0;
0399
0400 for (auto truthjet : *jetsMC)
0401 {
0402
0403
0404 bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0405 bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0406 if ((not eta_cut) or (not pt_cut)) continue;
0407 m_truthID.push_back(truthjet->get_id());
0408 m_truthNComponent.push_back(truthjet->size_comp());
0409 m_truthEta.push_back(truthjet->get_eta());
0410 m_truthPhi.push_back(truthjet->get_phi());
0411 m_truthE.push_back(truthjet->get_e());
0412 m_truthPt.push_back(truthjet->get_pt());
0413 m_nTruthJet++;
0414 }
0415 }
0416
0417
0418 if(m_doSeeds)
0419 {
0420 for (auto jet : *seedjetsraw)
0421 {
0422 int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0423 m_eta_rawseed.push_back(jet->get_eta());
0424 m_phi_rawseed.push_back(jet->get_phi());
0425 m_e_rawseed.push_back(jet->get_e());
0426 m_pt_rawseed.push_back(jet->get_pt());
0427 m_rawseed_cut.push_back(passesCut);
0428 }
0429
0430 for (auto jet : *seedjetssub)
0431 {
0432 int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0433 m_eta_subseed.push_back(jet->get_eta());
0434 m_phi_subseed.push_back(jet->get_phi());
0435 m_e_subseed.push_back(jet->get_e());
0436 m_pt_subseed.push_back(jet->get_pt());
0437 m_subseed_cut.push_back(passesCut);
0438 }
0439 }
0440
0441
0442 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0443 if (!gl1PacketInfo)
0444 {
0445 std::cout << PHWHERE << "caloTreeGen::process_event: GL1Packet node is missing. Output related to this node will be empty" << std::endl;
0446 }
0447
0448 if (gl1PacketInfo)
0449 {
0450 uint64_t triggervec = gl1PacketInfo->getTriggerVector();
0451 for (int i = 0; i < 64; i++)
0452 {
0453 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0454 m_triggerVector.push_back(trig_decision);
0455 triggervec = (triggervec >> 1U) & 0xffffffffU;
0456 }
0457 }
0458
0459
0460 RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0461 if (!clusterContainer)
0462 {
0463 std::cout << PHWHERE << "caloTreeGen::process_event: CLUSTERINFO_CEMC node is missing. Output related to this node will be empty" << std::endl;
0464 return 0;
0465 }
0466 RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0467 RawClusterContainer::ConstIterator clusterIter;
0468 for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0469 {
0470 RawCluster *recoCluster = clusterIter->second;
0471
0472 CLHEP::Hep3Vector vertex(0, 0, 0);
0473 if (m_zvtx != -9999)
0474 {
0475 vertex.setZ(m_zvtx);
0476 }
0477 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0478 CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, vertex);
0479
0480 float clusE = E_vec_cluster_Full.mag();
0481 float clusEcore = E_vec_cluster.mag();
0482 float clus_eta = E_vec_cluster.pseudoRapidity();
0483 float clus_phi = E_vec_cluster.phi();
0484 float clus_pt = E_vec_cluster.perp();
0485 float clus_prob = recoCluster->get_prob();
0486
0487 m_cle.push_back(clusE);
0488 m_clecore.push_back(clusEcore);
0489 m_cleta.push_back(clus_eta);
0490 m_clphi.push_back(clus_phi);
0491 m_clpt.push_back(clus_pt);
0492 m_clprob.push_back(clus_prob);
0493 }
0494
0495
0496 m_T->Fill();
0497
0498 return Fun4AllReturnCodes::EVENT_OK;
0499 }
0500
0501
0502 int JetValidation::ResetEvent(PHCompositeNode *topNode)
0503 {
0504
0505 m_id.clear();
0506 m_nComponent.clear();
0507 m_eta.clear();
0508 m_phi.clear();
0509 m_e.clear();
0510 m_pt.clear();
0511 m_unsub_pt.clear();
0512 m_sub_et.clear();
0513
0514 m_truthID.clear();
0515 m_truthNComponent.clear();
0516 m_truthEta.clear();
0517 m_truthPhi.clear();
0518 m_truthE.clear();
0519 m_truthPt.clear();
0520 m_truthdR.clear();
0521
0522 m_eta_subseed.clear();
0523 m_phi_subseed.clear();
0524 m_e_subseed.clear();
0525 m_pt_subseed.clear();
0526 m_subseed_cut.clear();
0527
0528 m_eta_rawseed.clear();
0529 m_phi_rawseed.clear();
0530 m_e_rawseed.clear();
0531 m_pt_rawseed.clear();
0532 m_rawseed_cut.clear();
0533
0534 m_triggerVector.clear();
0535
0536 m_cle.clear();
0537 m_clecore.clear();
0538 m_cleta.clear();
0539 m_clphi.clear();
0540 m_clpt.clear();
0541 m_clprob.clear();
0542 return Fun4AllReturnCodes::EVENT_OK;
0543 }
0544
0545
0546 int JetValidation::EndRun(const int runnumber)
0547 {
0548 std::cout << "JetValidation::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0549 return Fun4AllReturnCodes::EVENT_OK;
0550 }
0551
0552
0553 int JetValidation::End(PHCompositeNode *topNode)
0554 {
0555 std::cout << "JetValidation::End - Output to " << m_outputFileName << std::endl;
0556 PHTFileServer::get().cd(m_outputFileName);
0557
0558 m_T->Write();
0559 std::cout << "JetValidation::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0560 return Fun4AllReturnCodes::EVENT_OK;
0561 }
0562
0563
0564 int JetValidation::Reset(PHCompositeNode *topNode)
0565 {
0566 std::cout << "JetValidation::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0567 return Fun4AllReturnCodes::EVENT_OK;
0568 }
0569
0570
0571 void JetValidation::Print(const std::string &what) const
0572 {
0573 std::cout << "JetValidation::Print(const std::string &what) const Printing info for " << what << std::endl;
0574 }
