Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:19 

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