Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:15:27 

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 #include <fun4all/Fun4AllBase.h>
0004 #include <InclusiveJet.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 <globalvertex/MbdVertex.h>
0017 #include <globalvertex/MbdVertexMap.h>
0018 #include <calobase/RawTower.h>
0019 #include <calobase/RawTowerContainer.h>
0020 #include <calobase/RawTowerGeom.h>
0021 #include <calobase/RawTowerGeomContainer.h>
0022 #include <calobase/TowerInfoContainer.h>
0023 #include <calobase/TowerInfo.h>
0024 #include <ffarawobjects/Gl1Packet.h>
0025 #include <jetbackground/TowerBackground.h>
0026 #include <calobase/RawTowerDefs.h>
0027 #include <g4main/PHG4TruthInfoContainer.h>
0028 #include <g4main/PHG4Particle.h>
0029 #include <ffaobjects/EventHeaderv1.h>
0030 
0031 #include <calobase/RawCluster.h>
0032 #include <calobase/RawClusterContainer.h>
0033 #include <calobase/RawClusterUtility.h>
0034 #include <calobase/RawTowerGeomContainer.h>
0035 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0036 
0037 #include <trackbase_historic/SvtxTrack.h>
0038 #include <trackbase_historic/SvtxTrackMap.h>
0039 #include <trackbase_historic/SvtxTrackState.h>
0040 #include <trackbase/TrkrCluster.h>
0041 #include <trackbase/TrkrClusterv3.h>
0042 #include <trackbase/TrkrClusterv4.h>
0043 #include <trackbase/TrkrHit.h>
0044 #include <trackbase/TrkrHitSetContainer.h>
0045 #include <trackbase/TrkrDefs.h>
0046 #include <trackbase/ActsGeometry.h>
0047 #include <trackbase/ClusterErrorPara.h>
0048 #include <trackbase/TpcDefs.h>
0049 #include <trackbase/TrkrClusterContainer.h>
0050 #include <trackbase/TrkrHitSet.h>
0051 #include <trackbase/TrkrClusterHitAssoc.h>
0052 #include <trackbase/TrkrClusterIterationMapv1.h>
0053 
0054 #include <TTree.h>
0055 #include <iostream>
0056 #include <sstream>
0057 #include <iomanip>
0058 #include <cmath>
0059 #include <vector>
0060 
0061 #include "TH1.h"
0062 #include "TH2.h"
0063 #include "TFile.h"
0064 #include <TTree.h>
0065 #include <TVector3.h>
0066 
0067 //____________________________________________________________________________..
0068 InclusiveJet::InclusiveJet(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
0069   SubsysReco("InclusiveJet_" + recojetname + "_" + truthjetname)
0070   , m_recoJetName(recojetname)
0071   , m_truthJetName(truthjetname)
0072   , m_outputFileName(outputfilename)
0073   , m_etaRange(-0.7, 0.7)
0074   , m_ptRange(5, 100)
0075   , m_doTruthJets(0)
0076   , m_doTruth(0)
0077   , m_doTowers(0)
0078   , m_doEmcalClusters(0)
0079   , m_doTopoclusters(0)
0080   , m_doSeeds(0)
0081   , m_doTracks(0)
0082   , m_T(nullptr)
0083   , m_event(-1)
0084   , m_nTruthJet(-1)
0085   , m_nJet(-1)
0086   , m_triggerVector()
0087   , m_nComponent()
0088   , m_eta()
0089   , m_phi()
0090   , m_e()
0091   , m_pt()
0092   , m_jetEmcalE()
0093   , m_jetIhcalE()
0094   , m_jetOhcalE()
0095   , m_truthNComponent()
0096   , m_truthEta()
0097   , m_truthPhi()
0098   , m_truthE()
0099   , m_truthPt()
0100   , m_eta_rawseed()
0101   , m_phi_rawseed()
0102   , m_pt_rawseed()
0103   , m_e_rawseed()
0104   , m_rawseed_cut()
0105   , m_eta_subseed()
0106   , m_phi_subseed()
0107   , m_pt_subseed()
0108   , m_e_subseed()
0109   , m_subseed_cut()
0110 {
0111   std::cout << "InclusiveJet::InclusiveJet(const std::string &name) Calling ctor" << std::endl;
0112 }
0113 
0114 //____________________________________________________________________________..
0115 InclusiveJet::~InclusiveJet()
0116 {
0117   std::cout << "InclusiveJet::~InclusiveJet() Calling dtor" << std::endl;
0118 }
0119 
0120 //____________________________________________________________________________..
0121 //____________________________________________________________________________..
0122 int InclusiveJet::Init(PHCompositeNode *topNode)
0123 {
0124   std::cout << "InclusiveJet::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0125   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0126   std::cout << "InclusiveJet::Init - Output to " << m_outputFileName << std::endl;
0127 
0128   // configure Tree
0129   m_T = new TTree("T", "MyJetAnalysis Tree");
0130   m_T->Branch("m_event", &m_event, "event/I");
0131   m_T->Branch("nJet", &m_nJet, "nJet/I");
0132   m_T->Branch("zvtx", &m_zvtx);
0133   m_T->Branch("nComponent", &m_nComponent);
0134   m_T->Branch("triggerVector", &m_triggerVector);
0135 
0136   m_T->Branch("eta", &m_eta);
0137   m_T->Branch("phi", &m_phi);
0138   m_T->Branch("e", &m_e);
0139   m_T->Branch("pt", &m_pt);
0140 
0141   m_T->Branch("jetEmcalE", &m_jetEmcalE);
0142   m_T->Branch("jetIhcalE", &m_jetIhcalE);
0143   m_T->Branch("jetOhcalE", &m_jetOhcalE);
0144 
0145   if(m_doTruthJets){
0146     m_T->Branch("nTruthJet", &m_nTruthJet);
0147     m_T->Branch("truthNComponent", &m_truthNComponent);
0148     m_T->Branch("truthEta", &m_truthEta);
0149     m_T->Branch("truthPhi", &m_truthPhi);
0150     m_T->Branch("truthE", &m_truthE);
0151     m_T->Branch("truthPt", &m_truthPt);
0152   }
0153 
0154   if(m_doSeeds){
0155     m_T->Branch("rawseedEta", &m_eta_rawseed);
0156     m_T->Branch("rawseedPhi", &m_phi_rawseed);
0157     m_T->Branch("rawseedPt", &m_pt_rawseed);
0158     m_T->Branch("rawseedE", &m_e_rawseed);
0159     m_T->Branch("rawseedCut", &m_rawseed_cut);
0160     m_T->Branch("subseedEta", &m_eta_subseed);
0161     m_T->Branch("subseedPhi", &m_phi_subseed);
0162     m_T->Branch("subseedPt", &m_pt_subseed);
0163     m_T->Branch("subseedE", &m_e_subseed);
0164     m_T->Branch("subseedCut", &m_subseed_cut);
0165   }
0166 
0167   if(m_doTowers) {
0168     m_T->Branch("emcaln",&m_emcaln,"emcaln/I");
0169     m_T->Branch("emcale",m_emcale,"emcale[emcaln]/F");
0170     m_T->Branch("emcalchi2", m_emcalchi2,"emcalchi2[emcaln]/F");
0171     m_T->Branch("emcaleta",m_emcaleta,"emcaleta[emcaln]/F");
0172     m_T->Branch("emcalphi",m_emcalphi,"emcalphi[emcaln]/F");
0173     m_T->Branch("emcalieta",m_emcalieta,"emcalieta[emcaln]/I");
0174     m_T->Branch("emcaliphi",m_emcaliphi,"emcaliphi[emcaln]/I");
0175     m_T->Branch("emcalstatus",m_emcalstatus,"emcalstatus[emcaln]/I");
0176     //m_T->Branch("emcaltime",m_emcaltime,"emcaltime[emcaln]/F");
0177 
0178     m_T->Branch("ihcaln",&m_ihcaln,"ihcaln/I");
0179     m_T->Branch("ihcale",m_ihcale,"ihcale[ihcaln]/F");
0180     m_T->Branch("ihcalchi2", m_ihcalchi2,"ihcalchi2[ihcaln]/F");
0181     m_T->Branch("ihcaleta",m_ihcaleta,"ihcaleta[ihcaln]/F");
0182     m_T->Branch("ihcalphi",m_ihcalphi,"ihcalphi[ihcaln]/F");
0183     m_T->Branch("ihcalieta",m_ihcalieta,"ihcalieta[ihcaln]/I");
0184     m_T->Branch("ihcaliphi",m_ihcaliphi,"ihcaliphi[ihcaln]/I");
0185     m_T->Branch("ihcalstatus",m_ihcalstatus,"ihcalstatus[ihcaln]/I");
0186     //m_T->Branch("ihcaltime",m_ihcaltime,"ihcaltime[ihcaln]/F");
0187 
0188     m_T->Branch("ohcaln",&m_ohcaln,"ohcaln/I");
0189     m_T->Branch("ohcale",m_ohcale,"ohcale[ohcaln]/F");
0190     m_T->Branch("ohcalchi2", m_ohcalchi2,"ohcalchi2[ohcaln]/F");
0191     m_T->Branch("ohcaleta",m_ohcaleta,"ohcaleta[ohcaln]/F");
0192     m_T->Branch("ohcalphi",m_ohcalphi,"ohcalphi[ohcaln]/F");
0193     m_T->Branch("ohcalieta",m_ohcalieta,"ohcalieta[ohcaln]/I");
0194     m_T->Branch("ohcaliphi",m_ohcaliphi,"ohcaliphi[ohcaln]/I");
0195     m_T->Branch("ohcalstatus",m_ohcalstatus,"ohcalstatus[ohcaln]/I");
0196     //m_T->Branch("ohcaltime",m_ohcaltime,"ohcaltime[ohcaln]/F");
0197   }
0198 
0199   if(m_doEmcalClusters) {
0200     m_T->Branch("emcal_clsmult",&m_emcal_clsmult,"emcal_clsmult/I");
0201     m_T->Branch("emcal_cluster_e",m_emcal_cluster_e,"emcal_cluster_e[emcal_clsmult]/F");
0202     m_T->Branch("emcal_cluster_eta",m_emcal_cluster_eta,"emcal_cluster_eta[emcal_clsmult]/F");
0203     m_T->Branch("emcal_cluster_phi",m_emcal_cluster_phi,"emcal_cluster_phi[emcal_clsmult]/F");
0204   }
0205 
0206   if(m_doTopoclusters) {
0207     m_T->Branch("clsmult",&m_clsmult,"clsmult/I");
0208     m_T->Branch("cluster_e",m_cluster_e,"cluster_e[clsmult]/F");
0209     m_T->Branch("cluster_eta",m_cluster_eta,"cluster_eta[clsmult]/F");
0210     m_T->Branch("cluster_phi",m_cluster_phi,"cluster_phi[clsmult]/F");
0211     m_T->Branch("cluster_ntowers",m_cluster_ntowers,"cluster_ntowers[clsmult]/I");
0212     m_T->Branch("cluster_tower_e",m_cluster_tower_e,"cluster_tower_e[clsmult][500]/F");
0213     m_T->Branch("cluster_tower_calo",m_cluster_tower_calo,"cluster_tower_calo[clsmult][500]/I");
0214     m_T->Branch("cluster_tower_ieta",m_cluster_tower_ieta,"cluster_tower_ieta[clsmult][500]/I");
0215     m_T->Branch("cluster_tower_iphi",m_cluster_tower_iphi,"cluster_tower_iphi[clsmult][500]/I");
0216   }
0217 
0218   if(m_doTruth) {
0219     m_T->Branch("truthpar_n",&truthpar_n,"truthpar_n/I");
0220     m_T->Branch("truthpar_pz",truthpar_pz,"truthpar_pz[truthpar_n]/F");
0221     m_T->Branch("truthpar_pt",truthpar_pt,"truthpar_pt[truthpar_n]/F");
0222     m_T->Branch("truthpar_e",truthpar_e,"truthpar_e[truthpar_n]/F");
0223     m_T->Branch("truthpar_eta",truthpar_eta,"truthpar_eta[truthpar_n]/F");
0224     m_T->Branch("truthpar_phi",truthpar_phi,"truthpar_phi[truthpar_n]/F");
0225     m_T->Branch("truthpar_pid",truthpar_pid,"truthpar_pid[truthpar_n]/I");
0226   }
0227 
0228   if(m_doTracks) {
0229     m_T->Branch("trkmult",&m_trkmult,"trkmult/I");
0230     m_T->Branch("tr_p",m_tr_p,"tr_p[trkmult]/F");
0231     m_T->Branch("tr_pt",m_tr_pt,"tr_pt[trkmult]/F");
0232     m_T->Branch("tr_eta",m_tr_eta,"tr_eta[trkmult]/F");
0233     m_T->Branch("tr_phi",m_tr_phi,"tr_phi[trkmult]/F");
0234     m_T->Branch("tr_charge",m_tr_charge,"tr_charge[trkmult]/I");
0235     m_T->Branch("tr_chisq",m_tr_chisq,"tr_chisq[trkmult]/F");
0236     m_T->Branch("tr_ndf",m_tr_ndf,"tr_ndf[trkmult]/I");
0237     m_T->Branch("tr_nintt",m_tr_nintt,"tr_nintt[trkmult]/I");
0238     m_T->Branch("tr_nmaps",m_tr_nmaps,"tr_nmaps[trkmult]/I");
0239     m_T->Branch("tr_ntpc",m_tr_ntpc,"tr_ntpc[trkmult]/I");
0240     m_T->Branch("tr_quality",m_tr_quality,"tr_quality[trkmult]/F");
0241     m_T->Branch("tr_vertex_id",m_tr_vertex_id,"tr_vertex_id[trkmult]/I");
0242     m_T->Branch("tr_cemc_eta",m_tr_cemc_eta,"tr_cemc_eta[trkmult]/F");
0243     m_T->Branch("tr_cemc_phi",m_tr_cemc_phi,"tr_cemc_phi[trkmult]/F");
0244     m_T->Branch("tr_ihcal_eta",m_tr_ihcal_eta,"tr_ihcal_eta[trkmult]/F");
0245     m_T->Branch("tr_ihcal_phi",m_tr_ihcal_phi,"tr_ihcal_phi[trkmult]/F");
0246     m_T->Branch("tr_ohcal_eta",m_tr_ohcal_eta,"tr_ohcal_eta[trkmult]/F");
0247     m_T->Branch("tr_ohcal_phi",m_tr_ohcal_phi,"tr_ohcal_phi[trkmult]/F");
0248     m_T->Branch("tr_outer_cemc_eta",m_tr_outer_cemc_eta,"tr_outer_cemc_eta[trkmult]/F");
0249     m_T->Branch("tr_outer_cemc_phi",m_tr_outer_cemc_phi,"tr_outer_cemc_phi[trkmult]/F");
0250     m_T->Branch("tr_outer_ihcal_eta",m_tr_outer_ihcal_eta,"tr_outer_ihcal_eta[trkmult]/F");
0251     m_T->Branch("tr_outer_ihcal_phi",m_tr_outer_ihcal_phi,"tr_outer_ihcal_phi[trkmult]/F");
0252     m_T->Branch("tr_outer_ohcal_eta",m_tr_outer_ohcal_eta,"tr_outer_ohcal_eta[trkmult]/F");
0253     m_T->Branch("tr_outer_ohcal_phi",m_tr_outer_ohcal_phi,"tr_outer_ohcal_phi[trkmult]/F");
0254   }
0255   
0256   return Fun4AllReturnCodes::EVENT_OK;
0257 }
0258 
0259 //____________________________________________________________________________..
0260 int InclusiveJet::InitRun(PHCompositeNode *topNode)
0261 {
0262   std::cout << "InclusiveJet::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0263   
0264   return Fun4AllReturnCodes::EVENT_OK;
0265 }
0266 
0267 //____________________________________________________________________________..
0268 int InclusiveJet::process_event(PHCompositeNode *topNode)
0269 {
0270   //  std::cout << "InclusiveJet::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0271   ++m_event;
0272 
0273   // interface to reco jets
0274   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0275   if (!jets)
0276     {
0277       std::cout
0278     << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0279     << m_recoJetName << std::endl;
0280       exit(-1);
0281     }
0282 
0283   //interface to truth jets
0284   //JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0285   JetContainer* jetsMC = findNode::getClass<JetContainer>(topNode, m_truthJetName);
0286   if (!jetsMC && m_doTruthJets)
0287     {
0288       std::cout
0289     << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
0290     << m_truthJetName << std::endl;
0291       exit(-1);
0292     }
0293   
0294   // interface to jet seeds
0295   JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0296   if (!seedjetsraw && m_doSeeds)
0297     {
0298       std::cout
0299     << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
0300 << std::endl;
0301       exit(-1);
0302     }
0303 
0304   JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0305   if (!seedjetssub && m_doSeeds)
0306     {
0307       std::cout
0308 << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
0309 << std::endl;
0310       exit(-1);
0311     }
0312   
0313   
0314   //zvertex
0315   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0316   if (!vertexmap)
0317     {
0318       std::cout
0319         << "MyJetAnalysis::process_event - Error can not find global vertex  node "
0320         << std::endl;
0321       exit(-1);
0322     }
0323   if (vertexmap->empty())
0324     {
0325       std::cout
0326         << "MyJetAnalysis::process_event - global vertex node is empty "
0327         << std::endl;
0328       return Fun4AllReturnCodes::ABORTEVENT;
0329     }
0330   else
0331     {
0332       GlobalVertex *globalVertex = vertexmap->begin()->second;
0333       //std::cout << "Event: " << m_event << std::endl;
0334       //globalVertex->identify();
0335 
0336       auto mbdStartIter = globalVertex->find_vertexes(GlobalVertex::MBD);
0337       auto mbdEndIter = globalVertex->end_vertexes();
0338 
0339       for (auto iter = mbdStartIter; iter != mbdEndIter; ++iter)
0340       {
0341         const auto &[type, vertexVec] = *iter;
0342         if (type != GlobalVertex::MBD) continue;
0343         for (const auto *vertex : vertexVec)
0344         {
0345           if (!vertex) continue; 
0346           m_zvtx = vertex->get_z();
0347         }
0348       }
0349 
0350     }
0351 
0352     //std::cout << "using vertex " << m_zvtx << std::endl;
0353 
0354  /*
0355   MbdVertexMap *mvertexmap = findNode::getClass<MbdVertexMap>(topNode,"MbdVertexMap");
0356   if (!mvertexmap)
0357   {
0358     std::cout << "MbdVertexMap node is missing" << std::endl;
0359   }
0360   if (mvertexmap && !mvertexmap->empty())
0361   {
0362     MbdVertex *mvtx = mvertexmap->begin()->second;
0363     std::cout << "Event: " << m_event << std::endl;
0364     mvtx->identify();
0365     if (mvtx)
0366     {
0367         std::cout << "mbd vertex: " << mvtx->get_z() << std::endl;
0368     }
0369     else
0370     {
0371         std::cout << "no mbd vertex" << std::endl;
0372     }
0373   }
0374   */
0375 
0376   if (fabs(m_zvtx) > 30) {
0377     return Fun4AllReturnCodes::EVENT_OK;
0378   }
0379 
0380   //calorimeter towers
0381   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0382   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0383   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0384   RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0385   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0386   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0387   if(!towersEM3 || !towersIH3 || !towersOH3){
0388     std::cout
0389       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0390       << std::endl;
0391     exit(-1);
0392   }
0393 
0394   TowerInfoContainer *EMtowers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0395   TowerInfoContainer *IHtowers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0396   TowerInfoContainer *OHtowers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0397 
0398   if(!tower_geom || !tower_geomOH){
0399     std::cout
0400       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0401       << std::endl;
0402     exit(-1);
0403   }
0404 
0405   RawClusterContainer *topoclusters = findNode::getClass<RawClusterContainer>(topNode,"TOPOCLUSTER_ALLCALO"); 
0406   if (m_doTopoclusters && !topoclusters) {
0407         std::cout
0408     <<"MyJetAnalysis::process_event - Error can not find topoclusters "
0409     << std::endl;
0410   exit(-1);
0411   }
0412 
0413   RawClusterContainer *emcalclusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_CEMC");
0414   if (m_doEmcalClusters && !emcalclusters) {
0415       std::cout
0416     <<"MyJetAnalysis::process_event - Error can not find EMCal clusters "
0417     << std::endl;
0418   exit(-1);
0419   }
0420 
0421   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0422   if (m_doTruth && !truthinfo) {
0423     std::cout << PHWHERE << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"<< std::endl;
0424     return Fun4AllReturnCodes::EVENT_OK;
0425   }
0426 
0427   if(m_doTracks){
0428     RawTowerGeomContainer_Cylinderv1* cemcGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0429     RawTowerGeomContainer_Cylinderv1* ihcalGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0430     RawTowerGeomContainer_Cylinderv1* ohcalGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0431 
0432     if(!cemcGeomContainer){
0433       std::cout << "ERROR: TOWERGEOM_CEMC not found" << std::endl;
0434     }
0435     if(!ihcalGeomContainer){
0436       std::cout << "ERROR: TOWERGEOM_HCALIN not found" << std::endl;
0437     }
0438     if(!ohcalGeomContainer){
0439       std::cout << "ERROR: TOWERGEOM_HCALOUT not found" << std::endl;
0440     }
0441 
0442     m_cemcRadius = cemcGeomContainer->get_radius();
0443     m_ihcalRadius = ihcalGeomContainer->get_radius();
0444     m_ohcalRadius = ohcalGeomContainer->get_radius();
0445     m_cemcOuterRadius = cemcGeomContainer->get_radius() + cemcGeomContainer->get_thickness();
0446     m_ihcalOuterRadius = ihcalGeomContainer->get_radius() + ihcalGeomContainer->get_thickness();
0447     m_ohcalOuterRadius = ohcalGeomContainer->get_radius() + ohcalGeomContainer->get_thickness();
0448   }
0449 
0450   //get reco jets
0451   m_nJet = 0;
0452   float leadpt = 0;
0453   for (auto jet : *jets)
0454     {
0455       bool eta_cut = (jet->get_eta() >= m_etaRange.first) and (jet->get_eta() <= m_etaRange.second);
0456       bool pt_cut = (jet->get_pt() >= m_ptRange.first) and (jet->get_pt() <= m_ptRange.second);
0457       if ((not eta_cut) or (not pt_cut)) continue;
0458       //if(jet->get_e() < 0) {
0459       //  return Fun4AllReturnCodes::EVENT_OK; // currently applied to deal with cold EMCal IB
0460       //}
0461       m_nComponent.push_back(jet->size_comp());
0462       m_eta.push_back(jet->get_eta());
0463       m_phi.push_back(jet->get_phi());
0464       m_e.push_back(jet->get_e());
0465       m_pt.push_back(jet->get_pt());
0466 
0467       if (m_pt.back() > leadpt) { leadpt = m_pt.back(); }
0468 
0469       float emcalE = 0;
0470       float ihcalE = 0;
0471       float ohcalE = 0;
0472       for (auto comp: jet->get_comp_vec())
0473       {
0474         TowerInfo *tower;
0475         unsigned int channel = comp.second;
0476         if (comp.first == 14 || comp.first == 29 || comp.first == 28) {
0477           tower = towersEM3->get_tower_at_channel(channel);
0478           if(!tower) { continue; }
0479           emcalE += tower->get_energy();
0480         }    
0481         if (comp.first == 15 ||  comp.first == 30 || comp.first == 26)
0482         {
0483           tower = towersIH3->get_tower_at_channel(channel);
0484           if(!tower) { continue; }
0485           ihcalE += tower->get_energy();
0486         }
0487 
0488         if (comp.first == 16 ||  comp.first == 31 || comp.first == 27)
0489         {
0490           tower = towersOH3->get_tower_at_channel(channel);
0491           if(!tower) { continue; }
0492           ohcalE += tower->get_energy();
0493         }
0494       }
0495 
0496       m_jetEmcalE.push_back(emcalE);
0497       m_jetIhcalE.push_back(ihcalE);
0498       m_jetOhcalE.push_back(ohcalE);
0499 
0500       if (ohcalE/m_e.back() > 0.9) {
0501         std::cout << "event " << m_event << " emcalE " << emcalE << " ohcalE " << ohcalE << " totalE " << emcalE + ihcalE + ohcalE << " Jet E " << m_e.back() << std::endl;
0502         //return Fun4AllReturnCodes::EVENT_OK;
0503       }
0504       m_nJet++;
0505     }
0506 
0507     if (m_doLeadPtCut && leadpt < m_leadPtCut) { 
0508       return Fun4AllReturnCodes::EVENT_OK; 
0509     }
0510 
0511   //get truth jets
0512   if(m_doTruthJets)
0513     {
0514       m_nTruthJet = 0;
0515       //for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter)
0516       for (auto truthjet : *jetsMC)
0517     {
0518       //Jet* truthjet = iter->second;
0519         
0520       bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0521       bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
0522       if ((not eta_cut) or (not pt_cut)) continue;
0523       m_truthNComponent.push_back(truthjet->size_comp());
0524       m_truthEta.push_back(truthjet->get_eta());
0525       m_truthPhi.push_back(truthjet->get_phi());
0526       m_truthE.push_back(truthjet->get_e());
0527       m_truthPt.push_back(truthjet->get_pt());
0528       m_nTruthJet++;
0529     }
0530     }
0531   
0532   //get seed jets
0533   if(m_doSeeds)
0534     {
0535       for (auto jet : *seedjetsraw)
0536     {
0537       int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0538       m_eta_rawseed.push_back(jet->get_eta());
0539       m_phi_rawseed.push_back(jet->get_phi());
0540       m_e_rawseed.push_back(jet->get_e());
0541       m_pt_rawseed.push_back(jet->get_pt());
0542       m_rawseed_cut.push_back(passesCut);
0543     }
0544       
0545       for (auto jet : *seedjetssub)
0546     {
0547       int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0548       m_eta_subseed.push_back(jet->get_eta());
0549       m_phi_subseed.push_back(jet->get_phi());
0550       m_e_subseed.push_back(jet->get_e());
0551       m_pt_subseed.push_back(jet->get_pt());
0552       m_subseed_cut.push_back(passesCut);
0553     }
0554     }
0555 
0556   //grab the gl1 data
0557   Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0558   if (m_doTriggerCut && !gl1PacketInfo)
0559     {
0560       std::cout << PHWHERE << "caloTreeGen::process_event: GL1Packet node is missing. Output related to this node will be empty" << std::endl;
0561     }
0562 
0563   if (gl1PacketInfo) {
0564     bool jettrig = false;
0565     uint64_t triggervec = gl1PacketInfo->getScaledVector();
0566     for (int i = 0; i < 64; i++) {
0567         bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0568       if (trig_decision) {
0569         m_triggerVector.push_back(i);
0570         if ((i >= 16 && i <= 23) || (i >= 32 && i <= 35)) { jettrig = true; }
0571       }
0572         triggervec = triggervec >> 1U;
0573       }
0574     if (m_doTriggerCut && !jettrig) {
0575       return Fun4AllReturnCodes::EVENT_OK;
0576     }
0577   }
0578 
0579   if (m_doTowers) {
0580     if(EMtowers) { 
0581       m_emcaln = 0;
0582       int nchannels = 24576; 
0583       for(int i=0; i<nchannels; ++i) {
0584         TowerInfo *tower = EMtowers->get_tower_at_channel(i); //get EMCal tower
0585         //if(!tower->get_isGood()) { continue; }
0586         int key = EMtowers->encode_key(i);
0587         int etabin = EMtowers->getTowerEtaBin(key);
0588         int phibin = EMtowers->getTowerPhiBin(key);
0589         m_emcale[m_emcaln] = tower->get_energy(); 
0590         m_emcalchi2[m_emcaln] = tower->get_chi2();
0591         if (!tower->get_isBadChi2() && tower->get_chi2() > 10000) {
0592           std::cout << "event " << m_event << " ieta " << etabin << " iphi " << phibin << " chi2 " << tower->get_chi2() << " e " << tower->get_energy() << std::endl;
0593         }
0594         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, etabin, phibin);
0595         RawTowerGeom *geom = tower_geomEM->get_tower_geometry(geomkey); 
0596         TVector3 tower_pos;
0597         tower_pos.SetXYZ(geom->get_center_x(),geom->get_center_y(),geom->get_center_z() - m_zvtx);
0598         m_emcaleta[m_emcaln] = tower_pos.Eta();
0599         m_emcalphi[m_emcaln] = tower_pos.Phi();
0600         m_emcalieta[m_emcaln] = etabin;
0601         m_emcaliphi[m_emcaln] = phibin;
0602         m_emcalstatus[m_emcaln] = tower->get_status();
0603         //m_emcaltime[m_emcaln] = tower->get_time_float();
0604         m_emcaln++;
0605       }
0606     }
0607 
0608     if(IHtowers) { 
0609       m_ihcaln = 0;
0610       int nchannels = 1536; 
0611       for(int i=0; i<nchannels; ++i) {
0612         TowerInfo *tower = IHtowers->get_tower_at_channel(i); //get IHCal tower
0613         //if(!tower->get_isGood()) { continue; }
0614         int key = IHtowers->encode_key(i);
0615         int etabin = IHtowers->getTowerEtaBin(key);
0616         int phibin = IHtowers->getTowerPhiBin(key);
0617         m_ihcale[m_ihcaln] = tower->get_energy(); 
0618         m_ihcalchi2[m_ihcaln] = tower->get_chi2();
0619         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, etabin, phibin);
0620         RawTowerGeom *geom = tower_geom->get_tower_geometry(geomkey); 
0621         TVector3 tower_pos;
0622         tower_pos.SetXYZ(geom->get_center_x(),geom->get_center_y(),geom->get_center_z() - m_zvtx);
0623         m_ihcaleta[m_ihcaln] = tower_pos.Eta();
0624         m_ihcalphi[m_ihcaln] = tower_pos.Phi();
0625         m_ihcalieta[m_ihcaln] = etabin;
0626         m_ihcaliphi[m_ihcaln] = phibin;
0627         m_ihcalstatus[m_ihcaln] = tower->get_status();
0628         //m_ihcaltime[m_ihcaln] = tower->get_time_float();
0629         m_ihcaln++;
0630       }
0631     }
0632 
0633     if(OHtowers) { 
0634       m_ohcaln = 0;
0635       int nchannels = 1536; 
0636       for(int i=0; i<nchannels; ++i) {
0637         TowerInfo *tower = OHtowers->get_tower_at_channel(i); //get OHCal tower
0638         //if(!tower->get_isGood()) { continue; }
0639         int key = OHtowers->encode_key(i);
0640         int etabin = OHtowers->getTowerEtaBin(key);
0641         int phibin = OHtowers->getTowerPhiBin(key);
0642         m_ohcale[m_ohcaln] = tower->get_energy(); 
0643         m_ohcalchi2[m_ohcaln] = tower->get_chi2();
0644         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, etabin, phibin);
0645         RawTowerGeom *geom = tower_geomOH->get_tower_geometry(geomkey); 
0646         TVector3 tower_pos;
0647         tower_pos.SetXYZ(geom->get_center_x(),geom->get_center_y(),geom->get_center_z() - m_zvtx);
0648         m_ohcaleta[m_ohcaln] = tower_pos.Eta();
0649         m_ohcalphi[m_ohcaln] = tower_pos.Phi();
0650         m_ohcalieta[m_ohcaln] = etabin;
0651         m_ohcaliphi[m_ohcaln] = phibin;
0652         m_ohcalstatus[m_ohcaln] = tower->get_status();
0653         //m_ohcaltime[m_ohcaln] = tower->get_time_float();
0654         m_ohcaln++;
0655       }
0656     }
0657   }
0658   if (m_doEmcalClusters) {
0659     if (emcalclusters) {
0660       RawClusterContainer::Map clusterMap = emcalclusters->getClustersMap();
0661       m_emcal_clsmult = 0;
0662       for(auto entry : clusterMap){
0663         RawCluster* cluster = entry.second;
0664         CLHEP::Hep3Vector origin(0, 0, m_zvtx);
0665         m_emcal_cluster_e[m_emcal_clsmult] = cluster->get_energy();
0666         m_emcal_cluster_eta[m_emcal_clsmult] = RawClusterUtility::GetPseudorapidity(*cluster, origin);
0667         m_emcal_cluster_phi[m_emcal_clsmult] = RawClusterUtility::GetAzimuthAngle(*cluster, origin);
0668         m_emcal_clsmult++;
0669         if (m_emcal_clsmult == 10000) { break; }
0670       } 
0671     }
0672   }
0673   if (m_doTopoclusters) {
0674     if (topoclusters) {
0675       RawClusterContainer::Map clusterMap = topoclusters->getClustersMap();
0676       m_clsmult = 0;
0677       for(auto entry : clusterMap){
0678         RawCluster* cluster = entry.second;
0679         CLHEP::Hep3Vector origin(0, 0, m_zvtx);
0680         m_cluster_e[m_clsmult] = cluster->get_energy();
0681         m_cluster_eta[m_clsmult] = RawClusterUtility::GetPseudorapidity(*cluster, origin);
0682         m_cluster_phi[m_clsmult] = RawClusterUtility::GetAzimuthAngle(*cluster, origin);
0683         m_cluster_ntowers[m_clsmult] = (int)cluster->getNTowers();
0684         int m_clstower = 0;
0685         for (const auto& [tower_id, tower_e] : cluster->get_towermap())
0686         {
0687           /*
0688             if (m_cluster_e[m_clsmult] < -1.0) {
0689               std::cout << "Tower ID: " << tower_id
0690                         << " | Tower E: " << tower_e
0691                         << " | Calorimeter ID: " << static_cast<int>(RawTowerDefs::decode_caloid(tower_id))
0692                         << " | Index1: " << RawTowerDefs::decode_index1(tower_id)
0693                         << " | Index2: " << RawTowerDefs::decode_index2(tower_id)
0694                         << std::endl;
0695             }
0696             */
0697             m_cluster_tower_calo[m_clsmult][m_clstower] = static_cast<int>(RawTowerDefs::decode_caloid(tower_id));
0698             m_cluster_tower_ieta[m_clsmult][m_clstower]  = RawTowerDefs::decode_index1(tower_id);
0699             m_cluster_tower_iphi[m_clsmult][m_clstower]  = RawTowerDefs::decode_index2(tower_id);
0700             m_cluster_tower_e[m_clsmult][m_clstower] = tower_e;
0701             /*
0702             if (m_cluster_e[m_clsmult] < -1.0) {
0703               std::cout << "Tower ID: " << tower_id
0704                         << " | Tower E: " << tower_e << " " << m_cluster_tower_e[m_clsmult][m_clstower]
0705                         << " | Calorimeter ID: " << m_cluster_tower_calo[m_clsmult][m_clstower]
0706                         << " | Index1: " << m_cluster_tower_ieta[m_clsmult][m_clstower]
0707                         << " | Index2: " << m_cluster_tower_iphi[m_clsmult][m_clstower]
0708                         << std::endl;
0709             }
0710             */
0711             m_clstower++;
0712         }
0713         //std::cout << std::endl;
0714         m_clsmult++;
0715         if (m_clsmult == 10000) { break; }
0716       } 
0717     }
0718   }
0719 
0720   if (m_doTruth) {
0721     PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0722     truthpar_n = 0;
0723     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) {
0724       // Get truth particle
0725       const PHG4Particle *truth = iter->second;
0726       if (!truth) { std::cout << "missing particle" << std::endl; continue; }
0727       if (!truthinfo->is_primary(truth)) continue;
0728       /// Get this particles momentum, etc.
0729       truthpar_pt[truthpar_n] = sqrt(truth->get_px() * truth->get_px() + truth->get_py() * truth->get_py());
0730       truthpar_pz[truthpar_n] = truth->get_pz();
0731       truthpar_e[truthpar_n] = truth->get_e();
0732       truthpar_phi[truthpar_n] = atan2(truth->get_py(), truth->get_px());
0733       truthpar_eta[truthpar_n] = atanh(truth->get_pz() / sqrt(truth->get_px()*truth->get_px()+truth->get_py()*truth->get_py()+truth->get_pz()*truth->get_pz()));
0734       if (truthpar_eta[truthpar_n] != truthpar_eta[truthpar_n]) truthpar_eta[truthpar_n] = -999; // check for nans
0735       truthpar_pid[truthpar_n] = truth->get_pid();
0736       truthpar_n++;
0737       if(truthpar_n > 99999)
0738       {
0739         return Fun4AllReturnCodes::EVENT_OK;
0740       }
0741     }
0742   }
0743 
0744   if(m_doTracks) {
0745     SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0746     // Check whether SvtxTrackMap is present in the node tree or not
0747     if (!trackmap) {
0748       std::cout << PHWHERE << "SvtxTrackMap node is missing, can't collect tracks" << std::endl;
0749       return Fun4AllReturnCodes::EVENT_OK;
0750     }
0751 
0752     /*
0753     SvtxVertexMap *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0754     if (!vertexmap) {
0755       std::cout << PHWHERE << "SvtxVertexMap node is missing, can't collect tracks" << std::endl;
0756       return;
0757     }
0758     */
0759 
0760     /*
0761     // EvalStack for truth track matching
0762     if(!m_svtxEvalStack){ m_svtxEvalStack = new SvtxEvalStack(topNode); }
0763     m_svtxEvalStack->next_event(topNode);
0764     SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
0765     */
0766 
0767     m_trkmult = 0;
0768     // Looping over tracks in an event
0769     for(auto entry : *trackmap) {
0770       SvtxTrack *track = entry.second;
0771 
0772       m_tr_p[m_trkmult] = track->get_p();
0773       m_tr_pt[m_trkmult] = track->get_pt();
0774       m_tr_eta[m_trkmult] = track->get_eta();
0775       m_tr_phi[m_trkmult] = track->get_phi();
0776     
0777       m_tr_charge[m_trkmult] = track->get_charge();
0778       m_tr_chisq[m_trkmult] = track->get_chisq();
0779       m_tr_ndf[m_trkmult] = track->get_ndf();
0780       m_tr_quality[m_trkmult] = track->get_quality();
0781       m_tr_vertex_id[m_trkmult] = track->get_vertex_id();
0782 
0783       TrackSeed *silseed = track->get_silicon_seed();
0784       TrackSeed *tpcseed = track->get_tpc_seed();
0785 
0786       m_tr_nmaps[m_trkmult] = 0;
0787       m_tr_nintt[m_trkmult] = 0;
0788       m_tr_ntpc[m_trkmult] = 0;
0789 
0790       if(tpcseed) {
0791         for (TrackSeed::ConstClusterKeyIter iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter) {
0792           TrkrDefs::cluskey cluster_key = *iter;
0793           unsigned int layer = TrkrDefs::getLayer(cluster_key);
0794           if (_nlayers_maps > 0 && layer < _nlayers_maps) {
0795             m_tr_nmaps[m_trkmult]++;
0796           }
0797           if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) {
0798             m_tr_nintt[m_trkmult]++;
0799           }
0800           if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc)) {
0801             m_tr_ntpc[m_trkmult]++;
0802           }
0803         }
0804       } 
0805 
0806       if(silseed) {
0807         for (TrackSeed::ConstClusterKeyIter iter = silseed->begin_cluster_keys(); iter != silseed->end_cluster_keys(); ++iter) {
0808           TrkrDefs::cluskey cluster_key = *iter;
0809           unsigned int layer = TrkrDefs::getLayer(cluster_key);
0810           if (_nlayers_maps > 0 && layer < _nlayers_maps) {
0811             m_tr_nmaps[m_trkmult]++;
0812           }
0813           if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) {
0814             m_tr_nintt[m_trkmult]++;
0815           }
0816           if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc)) {
0817             m_tr_ntpc[m_trkmult]++;
0818             }
0819           }
0820         }
0821 
0822       // Project tracks to CEMC
0823       if(track->find_state(m_cemcRadius) != track->end_states()){
0824         m_tr_cemc_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_cemcRadius)->second);
0825         m_tr_cemc_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_cemcRadius)->second);
0826       }
0827       else{
0828         m_tr_cemc_eta[m_trkmult] = -999;
0829         m_tr_cemc_phi[m_trkmult] = -999;
0830       }
0831       
0832       // Project tracks to inner HCAL
0833       if(track->find_state(m_ihcalRadius) != track->end_states()){
0834         m_tr_ihcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ihcalRadius)->second);
0835         m_tr_ihcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ihcalRadius)->second);
0836       }
0837       else{
0838         m_tr_ihcal_eta[m_trkmult] = -999;
0839         m_tr_ihcal_phi[m_trkmult] = -999;
0840       }
0841 
0842       // Project tracks to outer HCAL
0843       if(track->find_state(m_ohcalRadius) != track->end_states()){
0844         m_tr_ohcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ohcalRadius)->second);
0845         m_tr_ohcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ohcalRadius)->second);
0846       }
0847       else{
0848         m_tr_ohcal_eta[m_trkmult] = -999;
0849         m_tr_ohcal_phi[m_trkmult] = -999;
0850       }
0851 
0852       if(track->find_state(m_cemcOuterRadius) != track->end_states()){
0853         m_tr_outer_cemc_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_cemcOuterRadius)->second);
0854         m_tr_outer_cemc_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_cemcOuterRadius)->second);
0855       }
0856       else{
0857         m_tr_outer_cemc_eta[m_trkmult] = -999;
0858         m_tr_outer_cemc_phi[m_trkmult] = -999;
0859       }
0860       
0861       // Project tracks to inner HCAL
0862       if(track->find_state(m_ihcalOuterRadius) != track->end_states()){
0863         m_tr_outer_ihcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ihcalOuterRadius)->second);
0864         m_tr_outer_ihcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ihcalOuterRadius)->second);
0865       }
0866       else{
0867         m_tr_outer_ihcal_eta[m_trkmult] = -999;
0868         m_tr_outer_ihcal_phi[m_trkmult] = -999;
0869       }
0870 
0871       // Project tracks to outer HCAL
0872       if(track->find_state(m_ohcalOuterRadius) != track->end_states()){
0873         m_tr_outer_ohcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ohcalOuterRadius)->second);
0874         m_tr_outer_ohcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ohcalOuterRadius)->second);
0875       }
0876       else{
0877         m_tr_outer_ohcal_eta[m_trkmult] = -999;
0878         m_tr_outer_ohcal_phi[m_trkmult] = -999;
0879       }
0880 
0881       /*
0882       // Matching SvtxTracks to G4 truth
0883       PHG4Particle *truthtrack = trackeval->max_truth_particle_by_nclusters(track);
0884 
0885       if(truthtrack){
0886         m_tr_truth_pid[m_trkmult] =truthtrack->get_pid();
0887         m_tr_truth_is_primary[m_trkmult] =truthinfo->is_primary(truthtrack);
0888     
0889         m_tr_truth_e[m_trkmult] =truthtrack->get_e();
0890 
0891         float px = truthtrack->get_px();
0892         float py = truthtrack->get_py();
0893         float pz = truthtrack->get_pz();
0894 
0895         m_tr_truth_pt[m_trkmult] = sqrt(px*px+py*py);
0896         m_tr_truth_phi[m_trkmult] = atan2(py, px);
0897         m_tr_truth_eta[m_trkmult] = atanh(pz/sqrt(px*px+py*py+pz*pz));
0898         m_tr_truth_track_id[m_trkmult] = truthtrack->get_track_id();
0899       } else{
0900         m_tr_truth_pid[m_trkmult] = -999;
0901         m_tr_truth_is_primary[m_trkmult] = -999;
0902         m_tr_truth_e[m_trkmult] = -999;
0903         m_tr_truth_pt[m_trkmult] = -999;
0904         m_tr_truth_eta[m_trkmult] = -999;
0905         m_tr_truth_phi[m_trkmult] = -999;
0906         m_tr_truth_track_id[m_trkmult] = -999;
0907       }
0908       */
0909       m_trkmult++;
0910     }
0911 
0912     /*
0913     m_vtxmult = 0;
0914     for(auto entry : *vertexmap)
0915     {
0916       SvtxVertex *vertex = entry.second;
0917       
0918       m_vertex_id[m_vtxmult] = vertex->get_id();
0919       m_vx[m_vtxmult] = vertex->get_x();
0920       m_vy[m_vtxmult] = vertex->get_y();
0921       m_vz[m_vtxmult] = vertex->get_z();
0922       m_vtxmult++;
0923     }
0924     */
0925   }
0926   
0927   //fill the tree
0928   m_T->Fill();
0929 
0930   return Fun4AllReturnCodes::EVENT_OK;
0931 }
0932 
0933 //____________________________________________________________________________..
0934 int InclusiveJet::ResetEvent(PHCompositeNode *topNode)
0935 {
0936   //std::cout << "InclusiveJet::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0937   m_nComponent.clear();
0938   m_eta.clear();
0939   m_phi.clear();
0940   m_e.clear();
0941   m_pt.clear();
0942 
0943   m_jetEmcalE.clear();
0944   m_jetIhcalE.clear();
0945   m_jetOhcalE.clear();
0946 
0947   m_truthNComponent.clear();
0948   m_truthEta.clear();
0949   m_truthPhi.clear();
0950   m_truthE.clear();
0951   m_truthPt.clear();
0952   m_truthdR.clear();
0953 
0954   m_eta_subseed.clear();
0955   m_phi_subseed.clear();
0956   m_e_subseed.clear();
0957   m_pt_subseed.clear();
0958   m_subseed_cut.clear();
0959 
0960   m_eta_rawseed.clear();
0961   m_phi_rawseed.clear();
0962   m_e_rawseed.clear();
0963   m_pt_rawseed.clear();
0964   m_rawseed_cut.clear();
0965   
0966   m_triggerVector.clear();
0967 
0968   m_zvtx = -9999;
0969   m_emcaln = 0;
0970   m_ihcaln = 0;
0971   m_ohcaln = 0;
0972   m_clsmult = 0;
0973   m_emcal_clsmult = 0;
0974   m_trkmult = 0;
0975 
0976   for (int i = 0; i < 24576; i++) {
0977     m_emcale[i] = 0;
0978     m_emcalieta[i] = 0;
0979     m_emcaliphi[i] = 0;
0980     m_emcaleta[i] = 0;
0981     m_emcalphi[i] = 0;
0982     m_emcalchi2[i] = 0;
0983     m_emcalstatus[i] = 0;
0984     //m_emcaltime[i] = 0;
0985   }
0986 
0987   for (int i = 0; i < 1536; i++) {
0988     m_ihcale[i] = 0;
0989     m_ihcalieta[i] = 0;
0990     m_ihcaliphi[i] = 0;
0991     m_ihcaleta[i] = 0;
0992     m_ihcalphi[i] = 0;
0993     m_ihcalchi2[i] = 0;
0994     m_ihcalstatus[i] = 0;
0995     //m_ihcaltime[i] = 0;
0996     m_ohcale[i] = 0;
0997     m_ohcalieta[i] = 0;
0998     m_ohcaliphi[i] = 0;
0999     m_ohcaleta[i] = 0;
1000     m_ohcalphi[i] = 0;
1001     m_ohcalchi2[i] = 0;
1002     m_ohcalstatus[i] = 0;
1003     //m_ohcaltime[i] = 0;
1004   }
1005 
1006   for (int i = 0; i < 100000; i++) {
1007     truthpar_pt[i] = 0;
1008     truthpar_pz[i] = 0;
1009     truthpar_e[i] = 0;
1010     truthpar_phi[i] = 0;
1011     truthpar_eta[i] = 0;
1012     truthpar_pid[i] = 0;
1013   }
1014 
1015   for (int i = 0; i < 2000; i++) {
1016     m_cluster_e[i] = 0;
1017     m_cluster_eta[i] = 0;
1018     m_cluster_phi[i] = 0;
1019     m_cluster_ntowers[i] = 0;
1020     m_emcal_cluster_e[i] = 0;
1021     m_emcal_cluster_eta[i] = 0;
1022     m_emcal_cluster_phi[i] = 0;
1023     for (int j = 0; j < 500; j++) {
1024       m_cluster_tower_calo[i][j] = 0;
1025       m_cluster_tower_ieta[i][j] = 0;
1026       m_cluster_tower_iphi[i][j] = 0;
1027       m_cluster_tower_e[i][j] = 0;
1028     }
1029   }
1030 
1031   for (int i = 0; i < 2000; i++) {
1032     m_tr_p[i] = 0;
1033     m_tr_pt[i] = 0;
1034     m_tr_eta[i] = 0;
1035     m_tr_phi[i] = 0;
1036     m_tr_charge[i] = 0;
1037     m_tr_chisq[i] = 0;
1038     m_tr_ndf[i] = 0;
1039     m_tr_nintt[i] = 0;
1040     m_tr_nmaps[i] = 0;
1041     m_tr_ntpc[i] = 0;
1042     m_tr_quality[i] = 0;
1043     m_tr_vertex_id[i] = 0;
1044     m_tr_cemc_eta[i] = 0; // Projection of track to calorimeters
1045     m_tr_cemc_phi[i] = 0;
1046     m_tr_ihcal_eta[i] = 0;
1047     m_tr_ihcal_phi[i] = 0;
1048     m_tr_ohcal_eta[i] = 0;
1049     m_tr_ohcal_phi[i] = 0;
1050     m_tr_outer_cemc_eta[i] = 0;
1051     m_tr_outer_cemc_phi[i] = 0;
1052     m_tr_outer_ihcal_eta[i] = 0;
1053     m_tr_outer_ihcal_phi[i] = 0;
1054     m_tr_outer_ohcal_eta[i] = 0;
1055     m_tr_outer_ohcal_phi[i] = 0;
1056   }
1057 
1058   return Fun4AllReturnCodes::EVENT_OK;
1059 }
1060 
1061 float InclusiveJet::calculateProjectionEta(SvtxTrackState* projectedState){
1062   float x = projectedState->get_x();// - initialState->get_x();
1063   float y = projectedState->get_y();// - initialState->get_y();
1064   float z = projectedState->get_z();// - initialState->get_z();
1065 
1066   float theta = acos(z / sqrt(x*x + y*y + z*z) );
1067 
1068   return -log( tan(theta/2.0) );
1069 }
1070 
1071 float InclusiveJet::calculateProjectionPhi(SvtxTrackState* projectedState){
1072   float x = projectedState->get_x();// - initialState->get_x();
1073   float y = projectedState->get_y();// - initialState->get_y();
1074  
1075   return atan2(y,x);
1076 }
1077 
1078 //____________________________________________________________________________..
1079 int InclusiveJet::EndRun(const int runnumber)
1080 {
1081   std::cout << "InclusiveJet::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
1082   return Fun4AllReturnCodes::EVENT_OK;
1083 }
1084 
1085 //____________________________________________________________________________..
1086 int InclusiveJet::End(PHCompositeNode *topNode)
1087 {
1088   std::cout << "InclusiveJet::End - Output to " << m_outputFileName << std::endl;
1089   PHTFileServer::get().cd(m_outputFileName);
1090 
1091   m_T->Write();
1092   std::cout << "InclusiveJet::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1093   return Fun4AllReturnCodes::EVENT_OK;
1094 }
1095 
1096 //____________________________________________________________________________..
1097 int InclusiveJet::Reset(PHCompositeNode *topNode)
1098 {
1099   std::cout << "InclusiveJet::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1100   return Fun4AllReturnCodes::EVENT_OK;
1101 }
1102 
1103 //____________________________________________________________________________..
1104 void InclusiveJet::Print(const std::string &what) const
1105 {
1106   std::cout << "InclusiveJet::Print(const std::string &what) const Printing info for " << what << std::endl;
1107 }
[ 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!