JetRTrack/src/jetrtrack.cc

0001 #include "jetrtrack.h"
0002 #include <fun4all/SubsysReco.h>  // for SubsysReco
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <TFile.h>
0006 #include <TH1F.h>
0007 #include <fun4all/Fun4AllHistoManager.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <g4jets/JetMap.h>
0011 #include <fun4all/Fun4AllServer.h>
0012 #include <fun4all/PHTFileServer.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014 #include <TTree.h>
0015 #include <phool/PHCompositeNode.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>    // for PHIODataNode
0019 #include <phool/PHNode.h>          // for PHNode
0020 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0021 #include <phool/PHObject.h>        // for PHObject
0022 #include <phool/getClass.h>
0023
0024 #include <iostream>
0025 #include <sstream>
0026 #include <string>
0027 #include <phool/phool.h>                    // for PHWHERE
0028
0029 #include <g4main/PHG4Particle.h>
0030 #include <g4main/PHG4TruthInfoContainer.h>
0031 #include <centrality/CentralityInfo.h>
0032
0033 #include <trackbase/TrkrCluster.h>
0034 #include <trackbase/TrkrClusterv3.h>
0035 #include <trackbase/TrkrClusterv4.h>
0036 #include <trackbase/TrkrHit.h>
0037 #include <trackbase/TrkrHitSetContainer.h>
0038 #include <trackbase/TrkrDefs.h>
0039 #include <trackbase/ActsGeometry.h>
0040 #include <trackbase/ClusterErrorPara.h>
0041
0042 #include <trackbase/TpcDefs.h>
0043
0044 #include <trackbase/TrkrClusterContainer.h>
0045 #include <trackbase/TrkrHitSet.h>
0046 #include <trackbase/TrkrClusterHitAssoc.h>
0047 #include <trackbase/TrkrClusterIterationMapv1.h>
0048
0049 #include <trackbase_historic/SvtxTrack.h>
0050 #include <trackbase_historic/SvtxTrackMap.h>
0051 #include <trackbase_historic/SvtxVertex.h>
0052 #include <trackbase_historic/SvtxVertexMap.h>
0053 #include <trackbase_historic/ActsTransformations.h>
0054 #include <trackbase_historic/TrackSeed.h>
0055 #include <g4main/PHG4Particle.h>
0056
0057 #include <g4main/PHG4TruthInfoContainer.h>
0058
0059
0060
0061 //____________________________________________________________________________..
0062 jetrtrack::jetrtrack(const std::string &name):
0063  SubsysReco(name)
0064 {
0065   std::cout << "jetrtrack::jetrtrack(const std::string &name) Calling ctor" << std::endl;
0066 }
0067
0068 //____________________________________________________________________________..
0069 jetrtrack::~jetrtrack()
0070 {
0071   std::cout << "jetrtrack::~jetrtrack() Calling dtor" << std::endl;
0072 }
0073 //____________________________________________________________________________..
0074
0075 int jetrtrack::Init(PHCompositeNode *topNode)
0076 {
0077   //-------------------------------------------------------------
0078   //print the list of nodes available in memory
0079   //-------------------------------------------------------------
0080   topNode->print();
0081   return Fun4AllReturnCodes::EVENT_OK;
0082 }
0083
0084 Fun4AllHistoManager *jetrtrack::get_HistoManager()  //This enables the handling of the output object and eventual writing of it to file
0085 {
0086   Fun4AllServer *se = Fun4AllServer::instance();
0087   Fun4AllHistoManager *hm = se->getHistoManager("jetrtrack_HISTOS");
0088
0089   if (not hm)
0090   {
0091     std:: cout
0092         << "jetrtrack::get_HistoManager - Making Fun4AllHistoManager jetrtrack_HISTOS"
0093         << std::endl;
0094     hm = new Fun4AllHistoManager("jetrtrack_HISTOS");
0095     se->registerHistoManager(hm);
0096   }
0097   assert(hm);
0098   return hm;
0099 }
0100
0101 //____________________________________________________________________________..
0102 int jetrtrack::InitRun(PHCompositeNode *topNode)
0103 {
0104   std::cout << "jetrtrack::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0105
0106   //-------------------------------
0107   //Create the output file
0108   //-------------------------------
0109   PHTFileServer::get().open(m_outfilename.c_str(), "RECREATE");
0110   // PHTFileServer::get().open("jettree.root", "RECREATE");
0111
0112
0113   Fun4AllHistoManager *hm = get_HistoManager();
0114   assert(hm);
0115
0116   //------------------------------------------------------------------------------------------------
0117   // create the output tree where jets, truth particles, and tracks will
0118   // be stored and register it to the histogram manager
0119   //------------------------------------------------------------------------------------------------
0120
0121   TTree* tree=  new TTree("tree","tree");
0122   tree->Print();
0123
0124   // Truth jet kinematics
0125   tree->Branch("tjet_pt",&m_tjet_pt);
0126   tree->Branch("tjet_phi",&m_tjet_phi);
0127   tree->Branch("tjet_eta",&m_tjet_eta);
0128   tree->Branch("tjet_jcpt",&m_tjet_jcpt);
0129
0130   //Reconstructed jet kinematics
0131   tree->Branch("rjet_pt",&m_rjet_pt);
0132   tree->Branch("rjet_phi",&m_rjet_phi);
0133   tree->Branch("rjet_eta",&m_rjet_eta);
0134
0135   //Reconstructed track kinematics
0136   tree->Branch("trk_pt",&m_trk_pt);
0137   tree->Branch("trk_eta",&m_trk_eta);
0138   tree->Branch("trk_phi",&m_trk_phi);
0139   tree->Branch("trk_quality",&m_trk_qual);
0140   tree->Branch("nmvtx_hits", &m_nmvtxhits);
0141
0142   //global variables
0143   tree->Branch("cent", &m_centrality);
0144   tree->Branch("zvtx", &m_zvtx);
0145
0146   //truth charged particle information
0147   tree->Branch("tp_pt",&m_tp_pt);
0148   tree->Branch("tp_px",&m_tp_px);
0149   tree->Branch("tp_py",&m_tp_py);
0150   tree->Branch("tp_pz",&m_tp_pz);
0151   tree->Branch("tp_phi",&m_tp_phi);
0152   tree->Branch("tp_eta",&m_tp_eta);
0153   tree->Branch("tp_pid",&m_tp_pid);
0154
0155   //Jet constituent information for truth jets
0156   tree->Branch("jc_pt",&m_jc_pt);
0157   tree->Branch("jc_phi",&m_jc_phi);
0158   tree->Branch("jc_eta",&m_jc_eta);
0159   tree->Branch("jc_tjet_index",&m_jc_index);
0160
0161   hm->registerHisto(tree);
0162
0163   return Fun4AllReturnCodes::EVENT_OK;
0164 }
0165
0166 //____________________________________________________________________________..
0167 int jetrtrack::process_event(PHCompositeNode *topNode)
0168 {
0169   double pi = TMath::Pi();
0170
0171   //-----------------------------------------------------------------------------------------------------------------
0172   // load in the containers which hold the truth and reconstructed information
0173   //-----------------------------------------------------------------------------------------------------------------
0174
0175   JetMap* tjets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r04");
0176   if (!tjets)
0177   {
0178     std::cout
0179         << "Jetrtrack::process_event - Error can not find DST JetMap node "
0180         << std::endl;
0181     exit(-1);
0182   }
0183
0184
0185   JetMap* rjets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r04_Sub1");
0186   if (!rjets)
0187   {
0188     std::cout
0189         << "Jetrtrack::process_event - Error can not find DST JetMap node "
0190         << std::endl;
0191     exit(-1);
0192   }
0193
0194   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0195   if (!trackmap)
0196   {
0197     std::cout
0198         << "Jetrtrack::process_event - Error can not find DST SvtxTrackMap node "
0199         << std::endl;
0200     exit(-1);
0201   }
0202
0203
0204   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0205   if (!cent_node)
0206     {
0207       std::cout
0208         << "Jetrtrack::process_event - Error can not find centrality node "
0209         << std::endl;
0210       exit(-1);
0211     }
0212
0213   m_centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);  //Centrality as defined by the impact parameter
0214
0215
0216
0217   SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0218  if (!vertexmap)
0219     {
0220       std::cout
0221         << "Jetrtrack::process_event - Error can not find vertex map node "
0222         << std::endl;
0223       exit(-1);
0224     }
0225
0226   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0227   if (!truthinfo)
0228     {
0229       std::cout
0230         << "Jetrtrack::process_event - Error can not find G4TruthInfo map node "
0231         << std::endl;
0232       exit(-1);
0233     }
0234
0235
0236
0237   //-------------------------------------------------------------
0238   //Call the histogram manager in order to
0239   //pull the tree into the perview of the event
0240   //-------------------------------------------------------------
0241
0242   Fun4AllHistoManager *hm = get_HistoManager();
0243   assert(hm);
0244
0245
0246
0247   //----------------------------------------------------
0248   // loop over vertex information and
0249   // extract z-vertex
0250   //----------------------------------------------------
0251
0252   if (vertexmap)
0253     {
0254       if (!vertexmap->empty())
0255     {
0256         SvtxVertex* vertex = (vertexmap->begin()->second);
0257         m_zvtx = vertex->get_z();
0258       }
0259     }
0260
0261
0262
0263   //---------------------------------------------------
0264   // loop over all truth jets in the event
0265   //---------------------------------------------------
0266   int chargedparticlespids[] = {11,211,321,2212}; // corresponds to electrons, charged pions, charged kaons, and protons
0267
0268   int jetnumber = 0;
0269   //---------------------------------
0270   // loop over all truth jets
0271   //---------------------------------
0272   for (JetMap::Iter iter = tjets->begin(); iter != tjets->end(); ++iter)
0273   {
0274     Jet* tjet = iter->second;
0275     assert(tjet);
0276     float tjet_phi = tjet->get_phi();
0277     float tjet_eta = tjet->get_eta();
0278     float tjet_pt = tjet->get_pt();
0279     if (tjet_pt < 5){continue;} //Filter out  low pT jets to reduce file size
0280
0281
0282     //---------------------------------------------------------
0283     //loop over truth jet constituent particles
0284     //extract the jet constitutent kinematics
0285     // from charged particles
0286     //---------------------------------------------------------
0287
0288     float sumjcpt = 0;
0289     for (Jet::ConstIter comp = tjet->begin_comp(); comp != tjet->end_comp(); ++comp)
0290       {
0291     PHG4Particle* truth = truthinfo->GetParticle((*comp).second);
0292     bool reject_particle = true;
0293     for (int k = 0 ; k < 4;k++)
0294        {
0295          if (abs(truth->get_pid()) == chargedparticlespids[k])
0296            {
0297          reject_particle = false;
0298            }
0299        }
0300     if (reject_particle) {continue;}
0301     float m_truthpx = truth->get_px();
0302     float m_truthpy = truth->get_py();
0303     float m_truthpz = truth->get_pz();
0304     float m_truthenergy = truth->get_e();
0305     float m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
0306     float m_truthphi = atan2(m_truthpy , m_truthpx);
0307     float m_trutheta = atanh(m_truthpz / m_truthenergy);
0308     m_jc_pt.push_back(m_truthpt);
0309     m_jc_phi.push_back(m_truthphi);
0310     m_jc_eta.push_back(m_trutheta);
0311     m_jc_index.push_back(jetnumber);
0312     sumjcpt += m_truthpt;
0313       }
0314     jetnumber++;
0315
0316
0317
0318     float drmin = 0.3;
0319     float matched_pt = -999;
0320     float matched_eta = -999;
0321     float matched_phi = -999;
0322     //---------------------------------------------------------------------------------------------------
0323     //loop over reconstructed jets and match each truth jet to the closest
0324     //reconstructed jet
0325     //----------------------------------------------------------------------------------------------------
0326     for (JetMap::Iter riter = rjets->begin(); riter != rjets->end(); ++riter)
0327       {
0328     Jet* rjet = riter->second;
0329     float rjet_phi = rjet->get_phi();
0330     float rjet_eta = rjet->get_eta();
0331     float rjet_pt = rjet->get_pt();
0332
0333         float deta = fabs(rjet_eta - tjet_eta);
0334     float dphi = fabs(rjet_phi - tjet_phi);
0335     if (dphi > pi)
0336       {
0337         dphi = 2*pi-dphi;
0338       }
0339
0340     float dr = TMath::Sqrt(dphi*dphi + deta*deta);
0341     if (dr < drmin)//truth-reco jet matching
0342       {
0343         matched_pt = rjet_pt;
0344         matched_eta = rjet_eta;
0345         matched_phi = rjet_phi;
0346         drmin = dr;
0347       }
0348       }
0349     //-------------------------------------------------------------------------
0350     //append the truth and matched reconstructed jet
0351     // to the end of the corresponding vectors
0352     //-------------------------------------------------------------------------
0353     if (tjet_pt > 0)
0354       {
0355
0356     m_tjet_pt.push_back(tjet_pt);
0357     m_tjet_jcpt.push_back(sumjcpt);
0358     m_tjet_phi.push_back(tjet_phi);
0359     m_tjet_eta.push_back(tjet_eta);
0360     m_rjet_pt.push_back(matched_pt);
0361     m_rjet_phi.push_back(matched_phi);
0362     m_rjet_eta.push_back(matched_eta);
0363     m_dr.push_back(drmin);
0364       }
0365     }
0366
0367   //-------------------------------------------------------------
0368   //loop over reconstructed tracks to extract
0369   //reconstructed track information
0370   //-------------------------------------------------------------
0371
0372  if (trackmap)
0373     {
0374       for (SvtxTrackMap::Iter iter = trackmap->begin();
0375            iter != trackmap->end();
0376            ++iter)
0377     {
0378       SvtxTrack* track = iter->second;
0379       float quality = track->get_quality();
0380       auto silicon_seed = track->get_silicon_seed();
0381       int nmvtxhits = 0;
0382       if (silicon_seed)
0383         {
0384           nmvtxhits = silicon_seed->size_cluster_keys();
0385         }
0386       if (quality < 10) //select on some minimal track quality as reccomended by tracking group
0387       {
0388         if (track->get_pt() < 1) {continue;} //exclude  low pT tracks to reduce file size
0389         m_trk_pt.push_back(track->get_pt());
0390         m_trk_eta.push_back(track->get_eta());
0391         m_trk_phi.push_back(track->get_phi());
0392         m_trk_qual.push_back(quality);
0393         m_nmvtxhits.push_back(nmvtxhits);
0394       }
0395       }
0396     }
0397
0398
0399  //------------------------------------------------------------------------------------------------------
0400  // loop over primary truth particles to extract truth particle information
0401  //------------------------------------------------------------------------------------------------------
0402  if (truthinfo)
0403    {
0404      PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0405      /// Loop over the G4 truth (stable) particles
0406      for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0407       iter != range.second;
0408       ++iter)
0409        {
0410      /// Get this truth particle
0411      const PHG4Particle *truth = iter->second;
0412
0413
0414      /// Get this particles momentum, etc.
0415      float m_truthpx = truth->get_px();
0416      float m_truthpy = truth->get_py();
0417      float m_truthpz = truth->get_pz();
0418      float m_truthenergy = truth->get_e();
0419
0420      float m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
0421      float m_truthphi = atan2(m_truthpy , m_truthpx);
0422      float m_trutheta = atanh(m_truthpz / m_truthenergy);
0423      int  m_truthpid = truth->get_pid();
0424
0425      if (m_truthpt < 1){continue;}
0426      //--------------------------------------------------------------------------------------------------
0427      //Filter truth particles to select those which are charged and stablish
0428      //--------------------------------------------------------------------------------------------------
0429      bool reject_particle = true;
0430      for (int k = 0 ; k < 4;k++)
0431        {
0432          if (abs(truth->get_pid()) == chargedparticlespids[k])
0433            {
0434          reject_particle = false;
0435            }
0436        }
0437      if (reject_particle) {continue;}
0438
0439      m_tp_pt.push_back(m_truthpt);
0440      m_tp_px.push_back(m_truthpx);
0441      m_tp_py.push_back(m_truthpy);
0442      m_tp_pz.push_back(m_truthpz);
0443      m_tp_phi.push_back(m_truthphi);
0444      m_tp_eta.push_back(m_trutheta);
0445      m_tp_pid.push_back(m_truthpid);
0446        }
0447    }
0448   //----------------------------------------------------------
0449   //Record the vectors of jet kinematic info
0450   //to the ttree.
0451   //----------------------------------------------------------
0452
0453   TTree *tree = dynamic_cast<TTree *>(hm->getHisto("tree"));
0454   assert(tree);
0455   tree->Fill();
0456
0457
0458   return Fun4AllReturnCodes::EVENT_OK;
0459 }
0460
0461 //____________________________________________________________________________..
0462 int jetrtrack::ResetEvent(PHCompositeNode *topNode)
0463 {
0464   //-----------------------------------------------------
0465   //Clear every vector inbetween events
0466   // this is CRITICAL!! Otherwise you will
0467   // have runaway behaviors
0468   //------------------------------------------------------
0469   m_tjet_pt.clear();
0470   m_tjet_jcpt.clear();
0471   m_tjet_phi.clear();
0472   m_tjet_eta.clear();
0473   m_nmvtxhits.clear();
0474
0475   m_jc_index.clear();
0476   m_tp_pt.clear();
0477   m_tp_px.clear();
0478   m_tp_py.clear();
0479   m_tp_pz.clear();
0480   m_tp_phi.clear();
0481   m_tp_eta.clear();
0482   m_tp_pid.clear();
0483
0484
0485   m_jc_pt.clear();
0486   m_jc_phi.clear();
0487   m_jc_eta.clear();
0488
0489   m_rjet_pt.clear();
0490   m_rjet_phi.clear();
0491   m_rjet_eta.clear();
0492   m_dr.clear();
0493
0494   m_trk_pt.clear();
0495   m_trk_phi.clear();
0496   m_trk_eta.clear();
0497   m_trk_qual.clear();
0498
0499   return Fun4AllReturnCodes::EVENT_OK;
0500 }
0501
0502 //____________________________________________________________________________..
0503 int jetrtrack::EndRun(const int runnumber)
0504 {
0505   std::cout << "jetrtrack::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0506   return Fun4AllReturnCodes::EVENT_OK;
0507 }
0508
0509 //____________________________________________________________________________..
0510 int jetrtrack::End(PHCompositeNode *topNode)
0511 {
0512   std::cout << "jetrtrack::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0513   //----------------------------------------------------------------------
0514   //Enter the created output file and write the ttree
0515   //----------------------------------------------------------------------
0516
0517   PHTFileServer::get().cd(m_outfilename.c_str());
0518   // PHTFileServer::get().cd("jettree.root");
0519
0520   Fun4AllHistoManager *hm = get_HistoManager();
0521   assert(hm);
0522
0523   for (unsigned int i = 0; i < hm->nHistos(); i++)
0524     hm->getHisto(i)->Write();
0525
0526
0527   return Fun4AllReturnCodes::EVENT_OK;
0528 }
0529
0530 //____________________________________________________________________________..
0531 int jetrtrack::Reset(PHCompositeNode *topNode)
0532 {
0533  std::cout << "jetrtrack::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0534   return Fun4AllReturnCodes::EVENT_OK;
0535 }
0536
0537 //____________________________________________________________________________..
0538 void jetrtrack::Print(const std::string &what) const
0539 {
0540   std::cout << "jetrtrack::Print(const std::string &what) const Printing info for " << what << std::endl;
0541 }