Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-18 09:20:45 

0001 /*!
0002  *  \file       PHG4TrackFastSimEval.cc
0003  *  \brief      Evaluation module for PHG4TrackFastSim output
0004  *  \details    input: PHG4TruthInfoContainer, SvtxTrackMap with SvtxTrack_FastSim inside
0005  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0006  */
0007 
0008 #include "PHG4TrackFastSimEval.h"
0009 
0010 #include <globalvertex/SvtxVertex.h>  // for SvtxVertex
0011 #include <globalvertex/SvtxVertexMap.h>
0012 
0013 #include <trackbase_historic/SvtxTrack.h>
0014 #include <trackbase_historic/SvtxTrackMap.h>
0015 #include <trackbase_historic/SvtxTrack_FastSim.h>
0016 
0017 #include <g4main/PHG4Hit.h>
0018 #include <g4main/PHG4HitContainer.h>
0019 #include <g4main/PHG4Particle.h>
0020 #include <g4main/PHG4TruthInfoContainer.h>
0021 #include <g4main/PHG4VtxPoint.h>
0022 
0023 #include <pdbcalbase/PdbParameterMap.h>
0024 
0025 #include <phparameter/PHParameters.h>
0026 
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028 #include <fun4all/PHTFileServer.h>
0029 #include <fun4all/SubsysReco.h>  // for SubsysReco
0030 
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>
0033 
0034 #include <TH2.h>
0035 #include <TSystem.h>
0036 #include <TTree.h>
0037 #include <TVector3.h>
0038 
0039 #include <cassert>
0040 #include <cmath>
0041 #include <iostream>
0042 #include <map>      // for _Rb_tree_const_ite...
0043 #include <utility>  // for pair
0044 
0045 #define LogDebug(exp) \
0046   (std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl)
0047 
0048 #define LogError(exp) \
0049   (std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl)
0050 
0051 #define LogWarning(exp) \
0052   (std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl)
0053 
0054 const std::string xyzt[] = {"x", "y", "z", "t"};
0055 
0056 //----------------------------------------------------------------------------//
0057 //-- Constructor:
0058 //--  simple initialization
0059 //----------------------------------------------------------------------------//
0060 PHG4TrackFastSimEval::PHG4TrackFastSimEval(const std::string &name, const std::string &filename, const std::string &trackmapname)
0061   : SubsysReco(name)
0062   , m_TruthInfoContainer(nullptr)
0063   , m_TrackMap(nullptr)
0064   , m_VertexMap(nullptr)
0065   , m_TracksEvalTree(nullptr)
0066   , m_VertexEvalTree(nullptr)
0067   , m_H2D_DeltaMomVsTruthMom(nullptr)
0068   , m_H2D_DeltaMomVsTruthEta(nullptr)
0069   , m_EventCounter(0)
0070   , m_OutFileName(filename)
0071   , m_TrackMapName(trackmapname)
0072 {
0073   reset_variables();
0074 }
0075 
0076 //----------------------------------------------------------------------------//
0077 //-- Init():
0078 //--   Intialize all histograms, trees, and ntuples
0079 //----------------------------------------------------------------------------//
0080 int PHG4TrackFastSimEval::Init(PHCompositeNode * /*topNode*/)
0081 {
0082   return Fun4AllReturnCodes::EVENT_OK;
0083 }
0084 
0085 //----------------------------------------------------------------------------//
0086 //-- InitRun():
0087 //--   add related hit object
0088 //----------------------------------------------------------------------------//
0089 int PHG4TrackFastSimEval::InitRun(PHCompositeNode *topNode)
0090 {
0091   if (Verbosity())
0092   {
0093     std::cout << PHWHERE << " Openning file " << m_OutFileName << std::endl;
0094   }
0095   PHTFileServer::open(m_OutFileName, "RECREATE");
0096 
0097   // create TTree
0098   m_TracksEvalTree = new TTree("tracks", "FastSim Eval => tracks");
0099   m_TracksEvalTree->Branch("event", &m_TTree_Event, "event/I");
0100   m_TracksEvalTree->Branch("gtrackID", &m_TTree_gTrackID, "gtrackID/I");
0101   m_TracksEvalTree->Branch("gflavor", &m_TTree_gFlavor, "gflavor/I");
0102   m_TracksEvalTree->Branch("gpx", &m_TTree_gpx, "gpx/F");
0103   m_TracksEvalTree->Branch("gpy", &m_TTree_gpy, "gpy/F");
0104   m_TracksEvalTree->Branch("gpz", &m_TTree_gpz, "gpz/F");
0105   m_TracksEvalTree->Branch("gvx", &m_TTree_gvx, "gvx/F");
0106   m_TracksEvalTree->Branch("gvy", &m_TTree_gvy, "gvy/F");
0107   m_TracksEvalTree->Branch("gvz", &m_TTree_gvz, "gvz/F");
0108   m_TracksEvalTree->Branch("gvt", &m_TTree_gvt, "gvt/F");
0109   m_TracksEvalTree->Branch("trackID", &m_TTree_TrackID, "trackID/I");
0110   m_TracksEvalTree->Branch("charge", &m_TTree_Charge, "charge/I");
0111   m_TracksEvalTree->Branch("nhits", &m_TTree_nHits, "nhits/I");
0112   m_TracksEvalTree->Branch("px", &m_TTree_px, "px/F");
0113   m_TracksEvalTree->Branch("py", &m_TTree_py, "py/F");
0114   m_TracksEvalTree->Branch("pz", &m_TTree_pz, "pz/F");
0115   m_TracksEvalTree->Branch("pcax", &m_TTree_pcax, "pcax/F");
0116   m_TracksEvalTree->Branch("pcay", &m_TTree_pcay, "pcay/F");
0117   m_TracksEvalTree->Branch("pcaz", &m_TTree_pcaz, "pcaz/F");
0118   m_TracksEvalTree->Branch("dca2d", &m_TTree_dca2d, "dca2d/F");
0119 
0120   // next a stat. on hits
0121   PHParameters PHG4TrackFastSim_Parameter("PHG4TrackFastSim");
0122 
0123   PdbParameterMap *nodeparams = findNode::getClass<PdbParameterMap>(topNode,
0124                                                                     "PHG4TrackFastSim_Parameter");
0125   if (!nodeparams)
0126   {
0127     std::cout << __PRETTY_FUNCTION__ << " : Warning, missing PHG4TrackFastSim_Parameter node and skip saving hits"
0128               << std::endl;
0129   }
0130   else
0131   {
0132     PHG4TrackFastSim_Parameter.FillFrom(nodeparams);
0133     if (Verbosity())
0134     {
0135       std::cout << __PRETTY_FUNCTION__ << " PHG4TrackFastSim_Parameter : ";
0136       PHG4TrackFastSim_Parameter.Print();
0137     }
0138 
0139     auto range = PHG4TrackFastSim_Parameter.get_all_int_params();
0140     for (auto iter = range.first; iter != range.second; ++iter)
0141     {
0142       const std::string &phg4hit_node_name = iter->first;
0143       const int &phg4hit_node_id = iter->second;
0144 
0145       std::cout << __PRETTY_FUNCTION__ << " Prepare PHG4Hit node name " << phg4hit_node_name
0146                 << " with ID = " << phg4hit_node_id << std::endl;
0147 
0148       std::string branch_name = std::string("nHit_") + phg4hit_node_name;
0149       m_TracksEvalTree->Branch(branch_name.c_str(),
0150                                &m_TTree_HitContainerID_nHits_map[phg4hit_node_id],
0151                                (branch_name + "/I").c_str());
0152     }
0153   }
0154 
0155   m_H2D_DeltaMomVsTruthEta = new TH2D("DeltaMomVsTruthEta",
0156                                       "#frac{#Delta p}{truth p} vs. truth #eta", 54, -4.5, +4.5, 1000, -1,
0157                                       1);
0158 
0159   m_H2D_DeltaMomVsTruthMom = new TH2D("DeltaMomVsTruthMom",
0160                                       "#frac{#Delta p}{truth p} vs. truth p", 41, -0.5, 40.5, 1000, -1,
0161                                       1);
0162 
0163   // create TTree - vertex
0164   m_VertexEvalTree = new TTree("vertex", "FastSim Eval => vertces");
0165   m_VertexEvalTree->Branch("event", &m_TTree_Event, "event/I");
0166   m_VertexEvalTree->Branch("gvx", &m_TTree_gvx, "gvx/F");
0167   m_VertexEvalTree->Branch("gvy", &m_TTree_gvy, "gvy/F");
0168   m_VertexEvalTree->Branch("gvz", &m_TTree_gvz, "gvz/F");
0169   m_VertexEvalTree->Branch("gvt", &m_TTree_gvt, "gvt/F");
0170   m_VertexEvalTree->Branch("vx", &m_TTree_vx, "vx/F");
0171   m_VertexEvalTree->Branch("vy", &m_TTree_vy, "vy/F");
0172   m_VertexEvalTree->Branch("vz", &m_TTree_vz, "vz/F");
0173   m_VertexEvalTree->Branch("deltavx", &m_TTree_DeltaVx, "deltavx/F");
0174   m_VertexEvalTree->Branch("deltavy", &m_TTree_DeltaVy, "deltavy/F");
0175   m_VertexEvalTree->Branch("deltavz", &m_TTree_DeltaVz, "deltavz/F");
0176   m_VertexEvalTree->Branch("gID", &m_TTree_gTrackID, "gID/I");
0177   m_VertexEvalTree->Branch("ID", &m_TTree_TrackID, "ID/I");
0178   m_VertexEvalTree->Branch("ntracks", &m_TTree_nTracks, "ntracks/I");
0179   m_VertexEvalTree->Branch("n_from_truth", &m_TTree_nFromTruth, "n_from_truth/I");
0180 
0181   for (std::map<std::string, unsigned int>::const_iterator iter = m_ProjectionNameMap.begin(); iter != m_ProjectionNameMap.end(); ++iter)
0182   {
0183     for (int i = 0; i < 4; i++)
0184     {
0185       std::string bname = iter->first + "_proj_" + xyzt[i];
0186       std::string bdef = bname + "/F";
0187 
0188       // fourth element is the path length
0189       if (i == 3)
0190       {
0191         bdef = iter->first + "_proj_path_length" + "/F";
0192       }
0193 
0194       m_TracksEvalTree->Branch(bname.c_str(), &m_TTree_proj_vec[iter->second][i], bdef.c_str());
0195     }
0196 
0197     for (int i = 0; i < 3; i++)
0198     {
0199       std::string bname = iter->first + "_proj_p" + xyzt[i];
0200       std::string bdef = bname + "/F";
0201       m_TracksEvalTree->Branch(bname.c_str(), &m_TTree_proj_p_vec[iter->second][i], bdef.c_str());
0202     }
0203     std::string nodename = "G4HIT_" + iter->first;
0204     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0205     if (hits)
0206     {
0207       for (int i = 0; i < 4; i++)
0208       {
0209         std::string bname = iter->first + "_" + xyzt[i];
0210         std::string bdef = bname + "/F";
0211         m_TracksEvalTree->Branch(bname.c_str(), &m_TTree_ref_vec[iter->second][i], bdef.c_str());
0212       }
0213       for (int i = 0; i < 3; i++)
0214       {
0215         std::string bname = iter->first + "_p" + xyzt[i];
0216         std::string bdef = bname + "/F";
0217 
0218         m_TracksEvalTree->Branch(bname.c_str(), &m_TTree_ref_p_vec[iter->second][i], bdef.c_str());
0219       }
0220     }
0221     if (!hits && Verbosity() > 0)
0222     {
0223       std::cout << "InitRun: could not find " << nodename << std::endl;
0224     }
0225   }
0226 
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229 
0230 //----------------------------------------------------------------------------//
0231 //-- process_event():
0232 //--   Call user instructions for every event.
0233 //--   This function contains the analysis structure.
0234 //----------------------------------------------------------------------------//
0235 int PHG4TrackFastSimEval::process_event(PHCompositeNode *topNode)
0236 {
0237   m_EventCounter++;
0238   if (Verbosity() >= 2 && m_EventCounter % 1000 == 0)
0239   {
0240     std::cout << PHWHERE << "Events processed: " << m_EventCounter << std::endl;
0241   }
0242 
0243   // std::cout << "Opening nodes" << std::endl;
0244   GetNodes(topNode);
0245 
0246   // std::cout << "Filling trees" << std::endl;
0247   fill_track_tree(topNode);
0248   fill_vertex_tree(topNode);
0249   // std::cout << "DONE" << std::endl;
0250 
0251   return Fun4AllReturnCodes::EVENT_OK;
0252 }
0253 
0254 //----------------------------------------------------------------------------//
0255 //-- End():
0256 //--   End method, wrap everything up
0257 //----------------------------------------------------------------------------//
0258 int PHG4TrackFastSimEval::End(PHCompositeNode * /*topNode*/)
0259 {
0260   PHTFileServer::cd(m_OutFileName);
0261 
0262   m_TracksEvalTree->Write();
0263   m_VertexEvalTree->Write();
0264 
0265   m_H2D_DeltaMomVsTruthEta->Write();
0266   m_H2D_DeltaMomVsTruthMom->Write();
0267 
0268   // PHTFileServer::get().close();
0269 
0270   return Fun4AllReturnCodes::EVENT_OK;
0271 }
0272 
0273 //----------------------------------------------------------------------------//
0274 //-- fill_tree():
0275 //--   Fill the trees with truth, track fit, and cluster information
0276 //----------------------------------------------------------------------------//
0277 void PHG4TrackFastSimEval::fill_track_tree(PHCompositeNode *topNode)
0278 {
0279   // Make sure to reset all the TTree variables before trying to set them.
0280 
0281   if (!m_TruthInfoContainer)
0282   {
0283     LogError("m_TruthInfoContainer not found!");
0284     return;
0285   }
0286 
0287   if (!m_TrackMap)
0288   {
0289     LogError("m_TrackMap not found!");
0290     return;
0291   }
0292 
0293   PHG4TruthInfoContainer::ConstRange range =
0294       m_TruthInfoContainer->GetPrimaryParticleRange();
0295   for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first;
0296        truth_itr != range.second; ++truth_itr)
0297   {
0298     reset_variables();
0299     m_TTree_Event = m_EventCounter;
0300 
0301     PHG4Particle *g4particle = truth_itr->second;
0302     if (!g4particle)
0303     {
0304       LogDebug("");
0305       continue;
0306     }
0307 
0308     SvtxTrack_FastSim *track = nullptr;
0309 
0310     if (Verbosity())
0311     {
0312       std::cout << __PRETTY_FUNCTION__ << "TRACKmap size " << m_TrackMap->size() << std::endl;
0313     }
0314     for (SvtxTrackMap::ConstIter track_itr = m_TrackMap->begin();
0315          track_itr != m_TrackMap->end();
0316          track_itr++)
0317     {
0318       // std::cout << "TRACK * " << track_itr->first << std::endl;
0319       SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(track_itr->second);
0320       if (!temp)
0321       {
0322         if (Verbosity())
0323         {
0324           std::cout << "PHG4TrackFastSimEval::fill_track_tree - ignore track that is not a SvtxTrack_FastSim:";
0325           track_itr->second->identify();
0326         }
0327         continue;
0328       }
0329       if (Verbosity())
0330       {
0331         std::cout << __PRETTY_FUNCTION__ << " PARTICLE!" << std::endl;
0332       }
0333 
0334       if ((temp->get_truth_track_id() - g4particle->get_track_id()) == 0)
0335       {
0336         track = temp;
0337       }
0338     }
0339 
0340     // std::cout << "B2" << std::endl;
0341     m_TTree_gTrackID = g4particle->get_track_id();
0342     m_TTree_gFlavor = g4particle->get_pid();
0343 
0344     m_TTree_gpx = g4particle->get_px();
0345     m_TTree_gpy = g4particle->get_py();
0346     m_TTree_gpz = g4particle->get_pz();
0347 
0348     m_TTree_gvx = std::numeric_limits<float>::quiet_NaN();
0349     m_TTree_gvy = std::numeric_limits<float>::quiet_NaN();
0350     m_TTree_gvz = std::numeric_limits<float>::quiet_NaN();
0351     m_TTree_gvt = std::numeric_limits<float>::quiet_NaN();
0352     PHG4VtxPoint *vtx = m_TruthInfoContainer->GetVtx(g4particle->get_vtx_id());
0353     if (vtx)
0354     {
0355       m_TTree_gvx = vtx->get_x();
0356       m_TTree_gvy = vtx->get_y();
0357       m_TTree_gvz = vtx->get_z();
0358       m_TTree_gvt = vtx->get_t();
0359     }
0360 
0361     if (track)
0362     {
0363       // std::cout << "C1" << std::endl;
0364       m_TTree_TrackID = track->get_id();
0365       m_TTree_Charge = track->get_charge();
0366       m_TTree_nHits = track->size_clusters();
0367 
0368       m_TTree_px = track->get_px();
0369       m_TTree_py = track->get_py();
0370       m_TTree_pz = track->get_pz();
0371       m_TTree_pcax = track->get_x();
0372       m_TTree_pcay = track->get_y();
0373       m_TTree_pcaz = track->get_z();
0374       m_TTree_dca2d = track->get_dca2d();
0375 
0376       TVector3 truth_mom(m_TTree_gpx, m_TTree_gpy, m_TTree_gpz);
0377       TVector3 reco_mom(m_TTree_px, m_TTree_py, m_TTree_pz);
0378       // std::cout << "C2" << std::endl;
0379 
0380       m_H2D_DeltaMomVsTruthMom->Fill(truth_mom.Mag(), (reco_mom.Mag() - truth_mom.Mag()) / truth_mom.Mag());
0381       m_H2D_DeltaMomVsTruthEta->Fill(truth_mom.Eta(), (reco_mom.Mag() - truth_mom.Mag()) / truth_mom.Mag());
0382       // find projections
0383       for (SvtxTrack::ConstStateIter trkstates = track->begin_states();
0384            trkstates != track->end_states();
0385            ++trkstates)
0386       {
0387         if (Verbosity())
0388         {
0389           std::cout << __PRETTY_FUNCTION__ << " checking " << trkstates->second->get_name() << std::endl;
0390         }
0391         std::map<std::string, unsigned int>::const_iterator iter = m_ProjectionNameMap.find(trkstates->second->get_name());
0392         if (iter != m_ProjectionNameMap.end())
0393         {
0394           if (Verbosity())
0395           {
0396             std::cout << __PRETTY_FUNCTION__ << " found " << trkstates->second->get_name() << std::endl;
0397           }
0398           // setting the projection (xyz and pxpypz)
0399           for (int i = 0; i < 3; i++)
0400           {
0401             m_TTree_proj_vec[iter->second][i] = trkstates->second->get_pos(i);
0402             m_TTree_proj_p_vec[iter->second][i] = trkstates->second->get_mom(i);
0403           }
0404           // fourth element is the path length
0405           m_TTree_proj_vec[iter->second][3] = trkstates->first;
0406 
0407           std::string nodename = "G4HIT_" + trkstates->second->get_name();
0408           PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0409           if (!hits)
0410           {
0411             if (Verbosity())
0412             {
0413               std::cout << __PRETTY_FUNCTION__ << " could not find " << nodename << std::endl;
0414             }
0415             continue;
0416           }
0417           if (Verbosity())
0418           {
0419             std::cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << std::endl;
0420           }
0421           PHG4HitContainer::ConstRange hit_range = hits->getHits();
0422           for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0423           {
0424             if (Verbosity())
0425             {
0426               std::cout << __PRETTY_FUNCTION__ << " checking hit id " << hit_iter->second->get_trkid() << " against " << track->get_truth_track_id() << std::endl;
0427             }
0428             if (hit_iter->second->get_trkid() - track->get_truth_track_id() == 0)
0429             {
0430               if (Verbosity())
0431               {
0432                 std::cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << std::endl;
0433               }
0434               if (iter->second > m_ProjectionNameMap.size())
0435               {
0436                 std::cout << "bad index: " << iter->second << std::endl;
0437                 gSystem->Exit(1);
0438               }
0439               m_TTree_ref_vec[iter->second][0] = hit_iter->second->get_x(0);
0440               m_TTree_ref_vec[iter->second][1] = hit_iter->second->get_y(0);
0441               m_TTree_ref_vec[iter->second][2] = hit_iter->second->get_z(0);
0442               m_TTree_ref_vec[iter->second][3] = hit_iter->second->get_t(0);
0443 
0444               m_TTree_ref_p_vec[iter->second][0] = hit_iter->second->get_px(0);
0445               m_TTree_ref_p_vec[iter->second][1] = hit_iter->second->get_py(0);
0446               m_TTree_ref_p_vec[iter->second][2] = hit_iter->second->get_pz(0);
0447             }
0448           }
0449         }
0450       }  // find projections
0451 
0452       // find number of hits
0453       for (const auto &g4hit_id_hitset : track->g4hit_ids())
0454       {
0455         const int &g4hit_id = g4hit_id_hitset.first;
0456         const std::set<PHG4HitDefs::keytype> &g4hit_set = g4hit_id_hitset.second;
0457         //
0458         auto nhit_iter = m_TTree_HitContainerID_nHits_map.find(g4hit_id);
0459         assert(nhit_iter != m_TTree_HitContainerID_nHits_map.end());
0460         //
0461         nhit_iter->second = g4hit_set.size();
0462       }
0463 
0464     }  //     if (track)
0465 
0466     m_TracksEvalTree->Fill();
0467   }  // PHG4TruthInfoContainer::ConstRange range =   m_TruthInfoContainer->GetPrimaryParticleRange();
0468 
0469   return;
0470 }
0471 
0472 //----------------------------------------------------------------------------//
0473 //-- fill_tree():
0474 //--   Fill the trees with truth, track fit, and cluster information
0475 //----------------------------------------------------------------------------//
0476 void PHG4TrackFastSimEval::fill_vertex_tree(PHCompositeNode * /*topNode*/)
0477 {
0478   if (!m_TruthInfoContainer)
0479   {
0480     LogError("m_TruthInfoContainer not found!");
0481     return;
0482   }
0483 
0484   if (!m_TrackMap)
0485   {
0486     LogError("m_TrackMap not found!");
0487     return;
0488   }
0489 
0490   if (!m_VertexMap)
0491   {
0492     return;
0493   }
0494 
0495   for (auto &iter : *m_VertexMap)
0496   {
0497     SvtxVertex *vertex = iter.second;
0498 
0499     // Make sure to reset all the TTree variables before trying to set them.
0500     reset_variables();
0501     // std::cout << "A1" << std::endl;
0502     m_TTree_Event = m_EventCounter;
0503 
0504     if (!vertex)
0505     {
0506       LogDebug("");
0507       continue;
0508     }
0509 
0510     // std::cout << "C1" << std::endl;
0511     m_TTree_TrackID = vertex->get_id();
0512     m_TTree_nTracks = vertex->size_tracks();
0513 
0514     m_TTree_vx = vertex->get_x();
0515     m_TTree_vy = vertex->get_y();
0516     m_TTree_vz = vertex->get_z();
0517     m_TTree_DeltaVx = std::sqrt(vertex->get_error(1, 1));
0518     m_TTree_DeltaVy = std::sqrt(vertex->get_error(2, 2));
0519     m_TTree_DeltaVz = std::sqrt(vertex->get_error(3, 3));
0520 
0521     // best matched vertex
0522     PHG4VtxPoint *best_vtx = nullptr;
0523     int best_n_match = -1;
0524     std::map<PHG4VtxPoint *, int> vertex_match_map;
0525     for (auto iterA = vertex->begin_tracks(); iterA != vertex->end_tracks(); ++iterA)
0526     {
0527       const auto &trackid = *iterA;
0528       const auto trackIter = m_TrackMap->find(trackid);
0529 
0530       if (trackIter == m_TrackMap->end())
0531       {
0532         continue;
0533       }
0534 
0535       SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(trackIter->second);
0536 
0537       if (!temp)
0538       {
0539         continue;
0540       }
0541 
0542       const auto g4trackID = temp->get_truth_track_id();
0543       const PHG4Particle *g4particle = m_TruthInfoContainer->GetParticle(g4trackID);
0544       assert(g4particle);
0545       PHG4VtxPoint *vtx = m_TruthInfoContainer->GetVtx(g4particle->get_vtx_id());
0546 
0547       int n_match = ++vertex_match_map[vtx];
0548 
0549       if (n_match > best_n_match)
0550       {
0551         best_n_match = n_match;
0552         best_vtx = vtx;
0553       }
0554     }
0555     if (best_vtx)
0556     {
0557       m_TTree_gvx = best_vtx->get_x();
0558       m_TTree_gvy = best_vtx->get_y();
0559       m_TTree_gvz = best_vtx->get_z();
0560       m_TTree_gvt = best_vtx->get_t();
0561 
0562       m_TTree_nFromTruth = best_n_match;
0563       m_TTree_gTrackID = best_vtx->get_id();
0564     }
0565     m_VertexEvalTree->Fill();
0566   }
0567   // std::cout << "B3" << std::endl;
0568 
0569   return;
0570 }
0571 
0572 //----------------------------------------------------------------------------//
0573 //-- reset_variables():
0574 //--   Reset all the tree variables to their default values.
0575 //--   Needs to be called at the start of every event
0576 //----------------------------------------------------------------------------//
0577 void PHG4TrackFastSimEval::reset_variables()
0578 {
0579   m_TTree_Event = -9999;
0580 
0581   //-- truth
0582   m_TTree_gTrackID = -9999;
0583   m_TTree_gFlavor = -9999;
0584   m_TTree_gpx = std::numeric_limits<float>::quiet_NaN();
0585   m_TTree_gpy = std::numeric_limits<float>::quiet_NaN();
0586   m_TTree_gpz = std::numeric_limits<float>::quiet_NaN();
0587 
0588   m_TTree_gvx = std::numeric_limits<float>::quiet_NaN();
0589   m_TTree_gvy = std::numeric_limits<float>::quiet_NaN();
0590   m_TTree_gvz = std::numeric_limits<float>::quiet_NaN();
0591   m_TTree_gvt = std::numeric_limits<float>::quiet_NaN();
0592 
0593   //-- reco
0594   m_TTree_TrackID = -9999;
0595   m_TTree_Charge = -9999;
0596   m_TTree_nHits = -9999;
0597   m_TTree_px = std::numeric_limits<float>::quiet_NaN();
0598   m_TTree_py = std::numeric_limits<float>::quiet_NaN();
0599   m_TTree_pz = std::numeric_limits<float>::quiet_NaN();
0600   m_TTree_pcax = std::numeric_limits<float>::quiet_NaN();
0601   m_TTree_pcay = std::numeric_limits<float>::quiet_NaN();
0602   m_TTree_pcaz = std::numeric_limits<float>::quiet_NaN();
0603   m_TTree_dca2d = std::numeric_limits<float>::quiet_NaN();
0604 
0605   m_TTree_vx = std::numeric_limits<float>::quiet_NaN();
0606   m_TTree_vy = std::numeric_limits<float>::quiet_NaN();
0607   m_TTree_vz = std::numeric_limits<float>::quiet_NaN();
0608   m_TTree_DeltaVx = std::numeric_limits<float>::quiet_NaN();
0609   m_TTree_DeltaVy = std::numeric_limits<float>::quiet_NaN();
0610   m_TTree_DeltaVz = std::numeric_limits<float>::quiet_NaN();
0611   m_TTree_nTracks = -9999;
0612   m_TTree_nFromTruth = -9999;
0613   for (auto &elem : m_TTree_proj_vec)
0614   {
0615     std::fill(elem.begin(), elem.end(), -9999);
0616   }
0617   for (auto &elem : m_TTree_proj_p_vec)
0618   {
0619     std::fill(elem.begin(), elem.end(), -9999);
0620   }
0621   for (auto &elem : m_TTree_ref_vec)
0622   {
0623     std::fill(elem.begin(), elem.end(), -9999);
0624   }
0625   for (auto &elem : m_TTree_ref_p_vec)
0626   {
0627     std::fill(elem.begin(), elem.end(), -9999);
0628   }
0629   for (auto &pair : m_TTree_HitContainerID_nHits_map)
0630   {
0631     pair.second = 0;
0632   }
0633 }
0634 
0635 //----------------------------------------------------------------------------//
0636 //-- GetNodes():
0637 //--   Get all the all the required nodes off the node tree
0638 //----------------------------------------------------------------------------//
0639 int PHG4TrackFastSimEval::GetNodes(PHCompositeNode *topNode)
0640 {
0641   // DST objects
0642   // Truth container
0643   m_TruthInfoContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0644                                                                     "G4TruthInfo");
0645   if (!m_TruthInfoContainer && m_EventCounter < 2)
0646   {
0647     std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
0648               << std::endl;
0649     return Fun4AllReturnCodes::ABORTEVENT;
0650   }
0651 
0652   m_TrackMap = findNode::getClass<SvtxTrackMap>(topNode,
0653                                                 m_TrackMapName);
0654   // std::cout << m_TrackMapName << std::endl;
0655   if (!m_TrackMap)
0656   {
0657     std::cout << PHWHERE << "SvtxTrackMap node with name "
0658               << m_TrackMapName
0659               << " not found on node tree"
0660               << std::endl;
0661     return Fun4AllReturnCodes::ABORTEVENT;
0662   }
0663 
0664   m_VertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0665   if (!m_VertexMap && Verbosity())
0666   {
0667     std::cout << PHWHERE << "SvtxTrackMap node with name SvtxVertexMap not found on node tree. Will not build the vertex eval tree"
0668               << std::endl;
0669   }
0670 
0671   return Fun4AllReturnCodes::EVENT_OK;
0672 }
0673 
0674 void PHG4TrackFastSimEval::AddProjection(const std::string &name)
0675 {
0676   std::vector<float> floatvec{-9999, -9999, -9999, -9999};
0677   m_TTree_proj_vec.push_back(floatvec);
0678   m_TTree_proj_p_vec.push_back(floatvec);
0679   m_TTree_ref_vec.push_back(floatvec);
0680   m_TTree_ref_p_vec.push_back(floatvec);
0681   // using m_ProjectionNameMap.size() makes sure it starts with 0 and then increments by 1
0682   // for each additional projection
0683   m_ProjectionNameMap.insert(std::make_pair(name, m_ProjectionNameMap.size()));
0684   return;
0685 }
[ 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!