Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:23 

0001 //module for producing a TTree with jet information for doing jet validation studies
0002 // for questions/bugs please contact Virginia Bailey vbailey13@gsu.edu
0003 
0004 #include "JetValidation.h"
0005 
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/PHTFileServer.h>
0008 
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011 
0012 #include <jetbase/JetMap.h>
0013 #include <jetbase/Jetv2.h>
0014 #include <jetbase/Jetv1.h>
0015 #include <jetbase/JetContainer.h>
0016 
0017 #include <centrality/CentralityInfo.h>
0018 #include <globalvertex/GlobalVertex.h>
0019 #include <globalvertex/GlobalVertexMap.h>
0020 
0021 #include <calobase/RawTower.h>
0022 #include <calobase/RawTowerContainer.h>
0023 #include <calobase/RawTowerGeom.h>
0024 #include <calobase/RawTowerGeomContainer.h>
0025 #include <calobase/TowerInfoContainer.h>
0026 #include <calobase/TowerInfo.h>
0027 
0028 #include <jetbackground/TowerBackground.h>
0029 
0030 #include <TTree.h>
0031 
0032 //____________________________________________________________________________..
0033 JetValidation::JetValidation(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
0034  SubsysReco("JetValidation_" + recojetname + "_" + truthjetname)
0035  , m_recoJetName(recojetname)
0036  , m_truthJetName(truthjetname)
0037  , m_outputFileName(outputfilename)
0038  , m_etaRange(-1, 1)
0039  , m_ptRange(5, 100)
0040  , m_doTruthJets(0)
0041  , m_doSeeds(0)
0042  , m_doUnsubJet(0)
0043  , m_T(nullptr)
0044  , m_event(-1)
0045  , m_nTruthJet(-1)
0046  , m_nJet(-1)
0047  , m_id()
0048  , m_nComponent()
0049  , m_eta()
0050  , m_phi()
0051  , m_e()
0052  , m_pt()
0053  , m_sub_et()
0054  , m_truthID()
0055  , m_truthNComponent()
0056  , m_truthEta()
0057  , m_truthPhi()
0058  , m_truthE()
0059  , m_truthPt()
0060  , m_eta_rawseed()
0061  , m_phi_rawseed()
0062  , m_pt_rawseed()
0063  , m_e_rawseed()
0064  , m_rawseed_cut()
0065  , m_eta_subseed()
0066  , m_phi_subseed()
0067  , m_pt_subseed()
0068  , m_e_subseed()
0069  , m_subseed_cut()
0070 {
0071   std::cout << "JetValidation::JetValidation(const std::string &name) Calling ctor" << std::endl;
0072 }
0073 
0074 //____________________________________________________________________________..
0075 JetValidation::~JetValidation()
0076 {
0077   std::cout << "JetValidation::~JetValidation() Calling dtor" << std::endl;
0078 }
0079 
0080 //____________________________________________________________________________..
0081 int JetValidation::Init(PHCompositeNode *topNode)
0082 {
0083   std::cout << "JetValidation::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0084   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0085   std::cout << "JetValidation::Init - Output to " << m_outputFileName << std::endl;
0086 
0087   // configure Tree
0088   m_T = new TTree("T", "MyJetAnalysis Tree");
0089   m_T->Branch("m_event", &m_event, "event/I");
0090   m_T->Branch("nJet", &m_nJet, "nJet/I");
0091   m_T->Branch("cent", &m_centrality);
0092   m_T->Branch("zvtx", &m_zvtx);
0093   m_T->Branch("b", &m_impactparam);
0094   m_T->Branch("id", &m_id);
0095   m_T->Branch("nComponent", &m_nComponent);
0096 
0097   m_T->Branch("eta", &m_eta);
0098   m_T->Branch("phi", &m_phi);
0099   m_T->Branch("e", &m_e);
0100   m_T->Branch("pt", &m_pt);
0101   if(m_doUnsubJet)
0102     {
0103       m_T->Branch("pt_unsub", &m_unsub_pt);
0104       m_T->Branch("subtracted_et", &m_sub_et);
0105     }
0106   if(m_doTruthJets){
0107     m_T->Branch("nTruthJet", &m_nTruthJet);
0108     m_T->Branch("truthID", &m_truthID);
0109     m_T->Branch("truthNComponent", &m_truthNComponent);
0110     m_T->Branch("truthEta", &m_truthEta);
0111     m_T->Branch("truthPhi", &m_truthPhi);
0112     m_T->Branch("truthE", &m_truthE);
0113     m_T->Branch("truthPt", &m_truthPt);
0114   }
0115 
0116   if(m_doSeeds){
0117     m_T->Branch("rawseedEta", &m_eta_rawseed);
0118     m_T->Branch("rawseedPhi", &m_phi_rawseed);
0119     m_T->Branch("rawseedPt", &m_pt_rawseed);
0120     m_T->Branch("rawseedE", &m_e_rawseed);
0121     m_T->Branch("rawseedCut", &m_rawseed_cut);
0122     m_T->Branch("subseedEta", &m_eta_subseed);
0123     m_T->Branch("subseedPhi", &m_phi_subseed);
0124     m_T->Branch("subseedPt", &m_pt_subseed);
0125     m_T->Branch("subseedE", &m_e_subseed);
0126     m_T->Branch("subseedCut", &m_subseed_cut);
0127   }
0128 
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131 
0132 //____________________________________________________________________________..
0133 int JetValidation::InitRun(PHCompositeNode *topNode)
0134 {
0135   std::cout << "JetValidation::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0136   return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138 
0139 //____________________________________________________________________________..
0140 int JetValidation::process_event(PHCompositeNode *topNode)
0141 {
0142   //std::cout << "JetValidation::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0143   ++m_event;
0144 
0145   // interface to reco jets
0146   JetMap* jets = findNode::getClass<JetMap>(topNode, m_recoJetName);
0147   if (!jets)
0148     {
0149       std::cout
0150     << "MyJetAnalysis::process_event - Error can not find DST Reco JetMap node "
0151     << m_recoJetName << std::endl;
0152       exit(-1);
0153     }
0154 
0155   //interface to truth jets
0156   JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0157   if (!jetsMC && m_doTruthJets)
0158     {
0159       std::cout
0160     << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0161     << m_truthJetName << std::endl;
0162       exit(-1);
0163     }
0164   
0165   // interface to jet seeds
0166   JetMap* seedjetsraw = findNode::getClass<JetMap>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0167   if (!seedjetsraw && m_doSeeds)
0168     {
0169       std::cout
0170     << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
0171     << std::endl;
0172       exit(-1);
0173     }
0174 
0175   JetMap* seedjetssub = findNode::getClass<JetMap>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0176   if (!seedjetssub && m_doSeeds)
0177     {
0178       std::cout
0179     << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
0180     << std::endl;
0181       exit(-1);
0182     }
0183 
0184   //centrality
0185   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0186   if (!cent_node)
0187     {
0188       std::cout
0189         << "MyJetAnalysis::process_event - Error can not find centrality node "
0190         << std::endl;
0191       exit(-1);
0192     }
0193 
0194   //zvertex
0195   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0196   GlobalVertex *vtx = vertexmap->begin()->second;
0197   m_zvtx = vtx->get_z();
0198 
0199 
0200   //calorimeter towers
0201   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0202   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0203   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0204   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0205   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0206   if(!towersEM3 || !towersIH3 || !towersOH3){
0207     std::cout
0208       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0209       << std::endl;
0210     exit(-1);
0211   }
0212 
0213   if(!tower_geom || !tower_geomOH){
0214     std::cout
0215       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0216       << std::endl;
0217     exit(-1);
0218   }
0219 
0220   //underlying event
0221   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0222   if(!background){
0223     std::cout<<"Can't get background. Exiting"<<std::endl;
0224     return Fun4AllReturnCodes::EVENT_OK;
0225   }
0226 
0227   //get the event centrality/impact parameter from HIJING
0228   m_centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0229   m_impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0230 
0231   //get reco jets
0232   m_nJet = 0;
0233   float background_v2 = 0;
0234   float background_Psi2 = 0;
0235   if(m_doUnsubJet)
0236     {
0237       background_v2 = background->get_v2();
0238       background_Psi2 = background->get_Psi2();
0239     }
0240   for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0241     {
0242 
0243       Jet* jet = iter->second;
0244 
0245       if(jet->get_pt() < 1) continue; // to remove noise jets
0246 
0247       m_id.push_back(jet->get_id());
0248       m_nComponent.push_back(jet->size_comp());
0249       m_eta.push_back(jet->get_eta());
0250       m_phi.push_back(jet->get_phi());
0251       m_e.push_back(jet->get_e());
0252       m_pt.push_back(jet->get_pt());
0253 
0254       if(m_doUnsubJet)
0255     {
0256       Jet* unsubjet = new Jetv1();
0257       float totalPx = 0;
0258       float totalPy = 0;
0259       float totalPz = 0;
0260       float totalE = 0;
0261       int nconst = 0;
0262 
0263       for (Jet::ConstIter comp = jet->begin_comp(); comp != jet->end_comp(); ++comp)
0264         {
0265           TowerInfo *tower;
0266           nconst++;
0267           unsigned int channel = (*comp).second;
0268           
0269           if ((*comp).first == 15 ||  (*comp).first == 30)
0270         {
0271           tower = towersIH3->get_tower_at_channel(channel);
0272           if(!tower || !tower_geom){
0273             continue;
0274           }
0275           unsigned int calokey = towersIH3->encode_key(channel);
0276           int ieta = towersIH3->getTowerEtaBin(calokey);
0277           int iphi = towersIH3->getTowerPhiBin(calokey);
0278           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0279           float UE = background->get_UE(1).at(ieta);
0280           float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
0281           float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
0282 
0283           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0284           totalE += tower->get_energy() + UE;
0285           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0286           totalPx += pt * cos(tower_phi);
0287           totalPy += pt * sin(tower_phi);
0288           totalPz += pt * sinh(tower_eta);
0289         }
0290           else if ((*comp).first == 16 || (*comp).first == 31)
0291         {
0292           tower = towersOH3->get_tower_at_channel(channel);
0293           if(!tower || !tower_geomOH)
0294             {
0295               continue;
0296             }
0297           
0298           unsigned int calokey = towersOH3->encode_key(channel);
0299           int ieta = towersOH3->getTowerEtaBin(calokey);
0300           int iphi = towersOH3->getTowerPhiBin(calokey);
0301           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0302           float UE = background->get_UE(2).at(ieta);
0303           float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0304           float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0305 
0306           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0307           totalE +=tower->get_energy() + UE;
0308           double pt = (tower->get_energy() + UE) / cosh(tower_eta);
0309           totalPx += pt * cos(tower_phi);
0310           totalPy += pt * sin(tower_phi);
0311           totalPz += pt * sinh(tower_eta);
0312         }
0313           else if ((*comp).first == 14 || (*comp).first == 29)
0314         {
0315           tower = towersEM3->get_tower_at_channel(channel);
0316           if(!tower || !tower_geom)
0317             {
0318               continue;
0319             }
0320 
0321           unsigned int calokey = towersEM3->encode_key(channel);
0322           int ieta = towersEM3->getTowerEtaBin(calokey);
0323           int iphi = towersEM3->getTowerPhiBin(calokey);
0324           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0325           float UE = background->get_UE(0).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         }
0337         }
0338       //get unsubtracted jet
0339       unsubjet->set_px(totalPx);
0340       unsubjet->set_py(totalPy);
0341       unsubjet->set_pz(totalPz);
0342       unsubjet->set_e(totalE);
0343       m_unsub_pt.push_back(unsubjet->get_pt());
0344       m_sub_et.push_back(unsubjet->get_et() - jet->get_et());
0345     }
0346 
0347       m_nJet++;
0348     }
0349 
0350   //get truth jets
0351   if(m_doTruthJets)
0352     {
0353       m_nTruthJet = 0;
0354       for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter)
0355     {
0356       Jet* truthjet = iter->second;
0357 
0358       bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0359       bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0360       if ((not eta_cut) or (not pt_cut)) continue;
0361       m_truthID.push_back(truthjet->get_id());
0362       m_truthNComponent.push_back(truthjet->size_comp());
0363       m_truthEta.push_back(truthjet->get_eta());
0364       m_truthPhi.push_back(truthjet->get_phi());
0365       m_truthE.push_back(truthjet->get_e());
0366       m_truthPt.push_back(truthjet->get_pt());
0367       m_nTruthJet++;
0368     }
0369     }
0370 
0371   //get seed jets
0372   if(m_doSeeds)
0373     {
0374       for (JetMap::Iter iter = seedjetsraw->begin(); iter != seedjetsraw->end(); ++iter)
0375     {
0376       Jet* jet = iter->second;
0377       int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
0378       m_eta_rawseed.push_back(jet->get_eta());
0379       m_phi_rawseed.push_back(jet->get_phi());
0380       m_e_rawseed.push_back(jet->get_e());
0381       m_pt_rawseed.push_back(jet->get_pt());
0382       m_rawseed_cut.push_back(passesCut);
0383     }
0384 
0385       for (JetMap::Iter iter = seedjetssub->begin(); iter != seedjetssub->end(); ++iter)
0386     {
0387       Jet* jet = iter->second;
0388       int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
0389       m_eta_subseed.push_back(jet->get_eta());
0390       m_phi_subseed.push_back(jet->get_phi());
0391       m_e_subseed.push_back(jet->get_e());
0392       m_pt_subseed.push_back(jet->get_pt());
0393       m_subseed_cut.push_back(passesCut);
0394     }
0395     }
0396   
0397   //fill the tree
0398   m_T->Fill();
0399 
0400   return Fun4AllReturnCodes::EVENT_OK;
0401 }
0402 
0403 //____________________________________________________________________________..
0404 int JetValidation::ResetEvent(PHCompositeNode *topNode)
0405 {
0406   //std::cout << "JetValidation::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0407   m_id.clear();
0408   m_nComponent.clear();
0409   m_eta.clear();
0410   m_phi.clear();
0411   m_e.clear();
0412   m_pt.clear();
0413   m_unsub_pt.clear();
0414   m_sub_et.clear();
0415 
0416   m_truthID.clear();
0417   m_truthNComponent.clear();
0418   m_truthEta.clear();
0419   m_truthPhi.clear();
0420   m_truthE.clear();
0421   m_truthPt.clear();
0422   m_truthdR.clear();
0423 
0424   m_eta_subseed.clear();
0425   m_phi_subseed.clear();
0426   m_e_subseed.clear();
0427   m_pt_subseed.clear();
0428   m_subseed_cut.clear();
0429 
0430   m_eta_rawseed.clear();
0431   m_phi_rawseed.clear();
0432   m_e_rawseed.clear();
0433   m_pt_rawseed.clear();
0434   m_rawseed_cut.clear();
0435   return Fun4AllReturnCodes::EVENT_OK;
0436 }
0437 
0438 //____________________________________________________________________________..
0439 int JetValidation::EndRun(const int runnumber)
0440 {
0441   std::cout << "JetValidation::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0442   return Fun4AllReturnCodes::EVENT_OK;
0443 }
0444 
0445 //____________________________________________________________________________..
0446 int JetValidation::End(PHCompositeNode *topNode)
0447 {
0448   std::cout << "JetValidation::End - Output to " << m_outputFileName << std::endl;
0449   PHTFileServer::get().cd(m_outputFileName);
0450 
0451   m_T->Write();
0452   std::cout << "JetValidation::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0453   return Fun4AllReturnCodes::EVENT_OK;
0454 }
0455 
0456 //____________________________________________________________________________..
0457 int JetValidation::Reset(PHCompositeNode *topNode)
0458 {
0459  std::cout << "JetValidation::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0460   return Fun4AllReturnCodes::EVENT_OK;
0461 }
0462 
0463 //____________________________________________________________________________..
0464 void JetValidation::Print(const std::string &what) const
0465 {
0466   std::cout << "JetValidation::Print(const std::string &what) const Printing info for " << what << std::endl;
0467 }
0468
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!