packages/KFParticle_sPHENIX/KFParticle_truthAndDetTools.cc

0001 #include "KFParticle_truthAndDetTools.h"
0002 #include "KFParticle_Tools.h"  // for KFParticle_Tools
0003
0004 #include <g4eval/SvtxEvalStack.h>   // for SvtxEvalStack
0005 #include <g4eval/SvtxTrackEval.h>   // for SvtxTrackEval
0006 #include <g4eval/SvtxTruthEval.h>   // for SvtxTruthEval
0007 #include <g4eval/SvtxVertexEval.h>  // for SvtxVertexEval
0008
0009 #include <globalvertex/SvtxVertex.h>         // for SvtxVertex
0010 #include <globalvertex/SvtxVertexMap.h>      // for SvtxVertexMap, SvtxVer...
0011 #include <trackbase/InttDefs.h>              // for getLadderPhiId, getLad...
0012 #include <trackbase/MvtxDefs.h>              // for getChipId, getStaveId
0013 #include <trackbase/TpcDefs.h>               // for getSectorId, getSide
0014 #include <trackbase/TrkrCluster.h>           // for TrkrCluster
0015 #include <trackbase/TrkrClusterContainer.h>  // for TrkrClusterContainer
0016 #include <trackbase/TrkrDefs.h>              // for getLayer, getTrkrId
0017 #include <trackbase_historic/SvtxPHG4ParticleMap.h>
0018 #include <trackbase_historic/SvtxTrack.h>     // for SvtxTrack, SvtxTrack::...
0019 #include <trackbase_historic/SvtxTrackMap.h>  // for SvtxTrackMap, SvtxTrac...
0020
0021 #include <globalvertex/GlobalVertex.h>
0022 #include <globalvertex/GlobalVertexMap.h>
0023
0024 #include <g4main/PHG4Particle.h>            // for PHG4Particle
0025 #include <g4main/PHG4TruthInfoContainer.h>  // for PHG4TruthInfoContainer
0026 #include <g4main/PHG4VtxPoint.h>            // for PHG4VtxPoint
0027
0028 #include <phhepmc/PHHepMCGenEvent.h>     // for PHHepMCGenEvent
0029 #include <phhepmc/PHHepMCGenEventMap.h>  // for PHHepMCGenEventMap
0030
0031 #include <phool/getClass.h>  // for getClass
0032
0033 #include <KFParticle.h>  // for KFParticle
0034
0035 #include <TSystem.h>     // for gSystem->Exit()
0036 #include <TTree.h>       // for TTree
0037
0038 #include <HepMC/GenEvent.h>       // for GenEvent::particle_con...
0039 #include <HepMC/GenParticle.h>    // for GenParticle
0040 #include <HepMC/GenVertex.h>      // for GenVertex::particle_it...
0041 #include <HepMC/IteratorRange.h>  // for parents
0042 #include <HepMC/SimpleVector.h>   // for FourVector
0043
0044 #include <algorithm>  // for max, find
0045 #include <cmath>      // for pow, sqrt
0046 #include <cstdlib>    // for NULL, abs
0047 #include <iostream>   // for operator<<, endl, basi...
0048 #include <iterator>   // for end, begin
0049 #include <map>        // for _Rb_tree_iterator, map
0050 #include <memory>     // for allocator_traits<>::va...
0051 #include <utility>    // for pair
0052
0053 std::map<std::string, int> Use =
0054     {
0055         {"MVTX", 1},
0056         {"INTT", 1},
0057         {"TPC", 1},
0058         {"TPOT", 1},
0059         {"EMCAL", 0},
0060         {"OHCAL", 0},
0061         {"IHCAL", 0}};
0062
0063 SvtxTrack *KFParticle_truthAndDetTools::getTrack(unsigned int track_id, SvtxTrackMap *trackmap)
0064 {
0065   SvtxTrack *matched_track = nullptr;
0066
0067   for (auto &iter : *trackmap)
0068   {
0069     if (iter.first == track_id)
0070     {
0071       matched_track = iter.second;
0072     }
0073   }
0074
0075   return matched_track;
0076 }
0077
0078 GlobalVertex *KFParticle_truthAndDetTools::getVertex(unsigned int vertex_id, GlobalVertexMap *vertexmap)
0079 {
0080   GlobalVertex *matched_vertex = vertexmap->get(vertex_id);
0081
0082   return matched_vertex;
0083 }
0084
0085 PHG4Particle *KFParticle_truthAndDetTools::getTruthTrack(SvtxTrack *thisTrack, PHCompositeNode *topNode)
0086 {
0087   /*
0088    * There are two methods for getting the truth rack from the reco track
0089    * 1. (recommended) Use the reco -> truth tables (requires SvtxPHG4ParticleMap). Introduced Summer of 2022
0090    * 2. Get truth track via nClusters. Older method and will work with older DSTs
0091    */
0092
0093   PHG4Particle *particle = nullptr;
0094
0095   SvtxPHG4ParticleMap *dst_reco_truth_map = findNode::getClass<SvtxPHG4ParticleMap>(topNode, "SvtxPHG4ParticleMap");
0096   if (dst_reco_truth_map)
0097   {
0098     m_truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0099     if (!m_truthinfo)
0100     {
0101       std::cout << "KFParticle truth matching: G4TruthInfo does not exist" << std::endl;
0102     }
0103
0104     std::map<float, std::set<int>> truth_set = dst_reco_truth_map->get(thisTrack->get_id());
0105     if (!truth_set.empty())
0106     {
0107       std::pair<float, std::set<int>> best_weight = *truth_set.rbegin();
0108       int best_truth_id = *best_weight.second.rbegin();
0109       particle = m_truthinfo->GetParticle(best_truth_id);
0110     }
0111   }
0112   else
0113   {
0114     std::cout << __FILE__ << ": SvtxPHG4ParticleMap not found, reverting to max_truth_particle_by_nclusters()" << std::endl;
0115
0116     if (!m_svtx_evalstack)
0117     {
0118       m_svtx_evalstack = new SvtxEvalStack(topNode);
0119       // clustereval = m_svtx_evalstack->get_cluster_eval();
0120       // hiteval = m_svtx_evalstack->get_hit_eval();
0121       trackeval = m_svtx_evalstack->get_track_eval();
0122       trutheval = m_svtx_evalstack->get_truth_eval();
0123       vertexeval = m_svtx_evalstack->get_vertex_eval();
0124     }
0125
0126     m_svtx_evalstack->next_event(topNode);
0127
0128     particle = trackeval->max_truth_particle_by_nclusters(thisTrack);
0129   }
0130   return particle;
0131 }
0132
0133 void KFParticle_truthAndDetTools::initializeTruthBranches(TTree *m_tree, int daughter_id, const std::string &daughter_number, bool m_constrain_to_vertex_truthMatch)
0134 {
0135   m_tree->Branch((daughter_number + "_true_ID").c_str(), &m_true_daughter_id[daughter_id], (daughter_number + "_true_ID/I").c_str());
0136   if (m_constrain_to_vertex_truthMatch)
0137   {
0138     m_tree->Branch((daughter_number + "_true_IP").c_str(), &m_true_daughter_ip[daughter_id], (daughter_number + "_true_IP/F").c_str());
0139     m_tree->Branch((daughter_number + "_true_IP_xy").c_str(), &m_true_daughter_ip_xy[daughter_id], (daughter_number + "_true_IP_xy/F").c_str());
0140   }
0141   m_tree->Branch((daughter_number + "_true_px").c_str(), &m_true_daughter_px[daughter_id], (daughter_number + "_true_px/F").c_str());
0142   m_tree->Branch((daughter_number + "_true_py").c_str(), &m_true_daughter_py[daughter_id], (daughter_number + "_true_py/F").c_str());
0143   m_tree->Branch((daughter_number + "_true_pz").c_str(), &m_true_daughter_pz[daughter_id], (daughter_number + "_true_pz/F").c_str());
0144   m_tree->Branch((daughter_number + "_true_p").c_str(), &m_true_daughter_p[daughter_id], (daughter_number + "_true_p/F").c_str());
0145   m_tree->Branch((daughter_number + "_true_pT").c_str(), &m_true_daughter_pt[daughter_id], (daughter_number + "_true_pT/F").c_str());
0146   m_tree->Branch((daughter_number + "_true_EV_x").c_str(), &m_true_daughter_vertex_x[daughter_id], (daughter_number + "_true_EV_x/F").c_str());
0147   m_tree->Branch((daughter_number + "_true_EV_y").c_str(), &m_true_daughter_vertex_y[daughter_id], (daughter_number + "_true_EV_y/F").c_str());
0148   m_tree->Branch((daughter_number + "_true_EV_z").c_str(), &m_true_daughter_vertex_z[daughter_id], (daughter_number + "_true_EV_z/F").c_str());
0149   if (m_constrain_to_vertex_truthMatch)
0150   {
0151     m_tree->Branch((daughter_number + "_true_PV_x").c_str(), &m_true_daughter_pv_x[daughter_id], (daughter_number + "_true_PV_x/F").c_str());
0152     m_tree->Branch((daughter_number + "_true_PV_y").c_str(), &m_true_daughter_pv_y[daughter_id], (daughter_number + "_true_PV_y/F").c_str());
0153     m_tree->Branch((daughter_number + "_true_PV_z").c_str(), &m_true_daughter_pv_z[daughter_id], (daughter_number + "_true_PV_z/F").c_str());
0154   }
0155   m_tree->Branch((daughter_number + "_true_track_history_PDG_ID").c_str(), &m_true_daughter_track_history_PDG_ID[daughter_id]);
0156   m_tree->Branch((daughter_number + "_true_track_history_PDG_mass").c_str(), &m_true_daughter_track_history_PDG_mass[daughter_id]);
0157   m_tree->Branch((daughter_number + "_true_track_history_px").c_str(), &m_true_daughter_track_history_px[daughter_id]);
0158   m_tree->Branch((daughter_number + "_true_track_history_py").c_str(), &m_true_daughter_track_history_py[daughter_id]);
0159   m_tree->Branch((daughter_number + "_true_track_history_pz").c_str(), &m_true_daughter_track_history_pz[daughter_id]);
0160   m_tree->Branch((daughter_number + "_true_track_history_pE").c_str(), &m_true_daughter_track_history_pE[daughter_id]);
0161   m_tree->Branch((daughter_number + "_true_track_history_pT").c_str(), &m_true_daughter_track_history_pT[daughter_id]);
0162 }
0163
0164 void KFParticle_truthAndDetTools::fillTruthBranch(PHCompositeNode *topNode, TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id, const KFParticle &kfvertex, bool m_constrain_to_vertex_truthMatch)
0165 {
0166   float true_px;
0167   float true_py;
0168   float true_pz;
0169   float true_p;
0170   float true_pt;
0171
0172   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
0173   if (!dst_trackmap)
0174   {
0175     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
0176   }
0177
0178   if (m_use_mbd_vertex_truth)
0179   {
0180     dst_mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0181     if (!dst_mbdvertexmap)
0182     {
0183       std::cout << "KFParticle truth matching: MbdVertexMap does not exist" << std::endl;
0184     }
0185   }
0186   else
0187   {
0188     dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vtx_map_node_name_nTuple);
0189     if (!dst_vertexmap)
0190     {
0191       std::cout << "KFParticle truth matching: " << m_vtx_map_node_name_nTuple << " does not exist" << std::endl;
0192     }
0193   }
0194
0195   auto *globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0196   if (!globalvertexmap)
0197   {
0198     std::cout << "KFParticle truth matching: GlobalVertexMap does not exist" << std::endl;
0199   }
0200
0201   track = getTrack(daughter.Id(), dst_trackmap);
0202   g4particle = getTruthTrack(track, topNode);
0203
0204   bool isParticleValid = g4particle == nullptr ? false : true;
0205
0206   true_px = isParticleValid ? (Float_t) g4particle->get_px() : 0.;
0207   true_py = isParticleValid ? (Float_t) g4particle->get_py() : 0.;
0208   true_pz = isParticleValid ? (Float_t) g4particle->get_pz() : 0.;
0209   true_p = sqrt(pow(true_px, 2) + pow(true_py, 2) + pow(true_pz, 2));
0210   true_pt = sqrt(pow(true_px, 2) + pow(true_py, 2));
0211
0212   m_true_daughter_px[daughter_id] = true_px;
0213   m_true_daughter_py[daughter_id] = true_py;
0214   m_true_daughter_pz[daughter_id] = true_pz;
0215   m_true_daughter_p[daughter_id] = true_p;
0216   m_true_daughter_pt[daughter_id] = true_pt;
0217   m_true_daughter_id[daughter_id] = isParticleValid ? g4particle->get_pid() : 0;
0218
0219   if (!m_svtx_evalstack)
0220   {
0221     m_svtx_evalstack = new SvtxEvalStack(topNode);
0222     // clustereval = m_svtx_evalstack->get_cluster_eval();
0223     // hiteval = m_svtx_evalstack->get_hit_eval();
0224     trackeval = m_svtx_evalstack->get_track_eval();
0225     trutheval = m_svtx_evalstack->get_truth_eval();
0226     vertexeval = m_svtx_evalstack->get_vertex_eval();
0227   }
0228
0229   if (isParticleValid)
0230   {
0231     g4vertex_point = trutheval->get_vertex(g4particle);
0232   }
0233
0234   m_true_daughter_vertex_x[daughter_id] = isParticleValid ? g4vertex_point->get_x() : 0.;
0235   m_true_daughter_vertex_y[daughter_id] = isParticleValid ? g4vertex_point->get_y() : 0.;
0236   m_true_daughter_vertex_z[daughter_id] = isParticleValid ? g4vertex_point->get_z() : 0.;
0237
0238   if (m_constrain_to_vertex_truthMatch)
0239   {
0240     // Calculate true DCA
0241     GlobalVertex *recoVertex = getVertex(kfvertex.Id(), globalvertexmap);
0242     GlobalVertex::VTXTYPE whichVtx = m_use_mbd_vertex_truth ? GlobalVertex::MBD : GlobalVertex::SVTX;
0243     auto svtxviter = recoVertex->find_vertexes(whichVtx);
0244
0245     // check that it contains a track vertex
0246     if (svtxviter == recoVertex->end_vertexes())
0247     {
0248       std::string vtxType = m_use_mbd_vertex_truth ? "MBD" : "silicon";
0249       std::cout << "Have a global vertex with no " << vtxType << " vertex... shouldn't happen in KFParticle_truthAndDetTools::fillTruthBranch..." << std::endl;
0250     }
0251
0252     auto svtxvertexvector = svtxviter->second;
0253     MbdVertex *mbdvertex = nullptr;
0254     SvtxVertex *svtxvertex = nullptr;
0255
0256     for (auto &vertex_iter : svtxvertexvector)
0257     {
0258       if (m_use_mbd_vertex_truth)
0259       {
0260         mbdvertex = dst_mbdvertexmap->find(vertex_iter->get_id())->second;
0261       }
0262       else
0263       {
0264         svtxvertex = dst_vertexmap->find(vertex_iter->get_id())->second;
0265       }
0266     }
0267
0268     PHG4VtxPoint *truePoint = nullptr;
0269     if (m_use_mbd_vertex_truth)
0270     {
0271       std::set<PHG4VtxPoint *> truePointSet = vertexeval->all_truth_points(mbdvertex);
0272       truePoint = *truePointSet.begin();
0273     }
0274     else
0275     {
0276       truePoint = vertexeval->max_truth_point_by_ntracks(svtxvertex);
0277     }
0278
0279     if (truePoint == nullptr && isParticleValid)
0280     {
0281       // PHG4Particle *g4mother = m_truthinfo->GetParticle(g4particle->get_parent_id());
0282       PHG4Particle *g4mother = m_truthinfo->GetPrimaryParticle(g4particle->get_parent_id());
0283       if (!g4mother)
0284       {
0285         std::cout << "KFParticle truth matching: True mother not found!\n";
0286         std::cout << "Your truth track DCA will be measured wrt a reconstructed vertex!" << std::endl;
0287         truePoint = nullptr;
0288       }
0289       else
0290       {
0291         truePoint = m_truthinfo->GetVtx(g4mother->get_vtx_id());  // Note, this may not be the PV for a decay with tertiaries
0292       }
0293     }
0294
0295     KFParticle trueKFParticleVertex;
0296
0297     float f_vertexParameters[6] = {0};
0298
0299     if (truePoint == nullptr)
0300     {
0301       std::cout << "KFParticle truth matching: This event has no PHG4VtxPoint information!\n";
0302       std::cout << "Your truth track DCA will be measured wrt a reconstructed vertex!" << std::endl;
0303
0304       f_vertexParameters[0] = recoVertex->get_x();
0305       f_vertexParameters[1] = recoVertex->get_y();
0306       f_vertexParameters[2] = recoVertex->get_z();
0307     }
0308     else
0309     {
0310       f_vertexParameters[0] = truePoint->get_x();
0311       f_vertexParameters[1] = truePoint->get_y();
0312       f_vertexParameters[2] = truePoint->get_z();
0313     }
0314
0315     float f_vertexCovariance[21] = {0};
0316
0317     trueKFParticleVertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
0318
0319     KFParticle trueKFParticle;
0320
0321     float f_trackParameters[6] = {m_true_daughter_vertex_x[daughter_id],
0322                                   m_true_daughter_vertex_y[daughter_id],
0323                                   m_true_daughter_vertex_z[daughter_id],
0324                                   true_px,
0325                                   true_py,
0326                                   true_pz};
0327
0328     float f_trackCovariance[21] = {0};
0329
0330     trueKFParticle.Create(f_trackParameters, f_trackCovariance, 1, -1);
0331
0332     m_true_daughter_ip[daughter_id] = trueKFParticle.GetDistanceFromVertex(trueKFParticleVertex);
0333     m_true_daughter_ip_xy[daughter_id] = trueKFParticle.GetDistanceFromVertexXY(trueKFParticleVertex);
0334
0335     m_true_daughter_pv_x[daughter_id] = truePoint == nullptr ? -99. : truePoint->get_x();
0336     m_true_daughter_pv_y[daughter_id] = truePoint == nullptr ? -99. : truePoint->get_y();
0337     m_true_daughter_pv_z[daughter_id] = truePoint == nullptr ? -99. : truePoint->get_z();
0338   }
0339 }
0340
0341 void KFParticle_truthAndDetTools::fillGeant4Branch(PHG4Particle *particle, int daughter_id)
0342 {
0343   Float_t pT = sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2));
0344
0345   m_true_daughter_track_history_PDG_ID[daughter_id].push_back(particle->get_pid());
0346   m_true_daughter_track_history_PDG_mass[daughter_id].push_back(0);
0347   m_true_daughter_track_history_px[daughter_id].push_back((Float_t) particle->get_px());
0348   m_true_daughter_track_history_py[daughter_id].push_back((Float_t) particle->get_py());
0349   m_true_daughter_track_history_pz[daughter_id].push_back((Float_t) particle->get_pz());
0350   m_true_daughter_track_history_pE[daughter_id].push_back((Float_t) particle->get_e());
0351   m_true_daughter_track_history_pT[daughter_id].push_back(pT);
0352 }
0353
0354 void KFParticle_truthAndDetTools::fillHepMCBranch(HepMC::GenParticle *particle, int daughter_id)
0355 {
0356   const HepMC::FourVector &myFourVector = particle->momentum();
0357
0358   m_true_daughter_track_history_PDG_ID[daughter_id].push_back(particle->pdg_id());
0359   m_true_daughter_track_history_PDG_mass[daughter_id].push_back((Float_t) particle->generatedMass());
0360   m_true_daughter_track_history_px[daughter_id].push_back((Float_t) myFourVector.px());
0361   m_true_daughter_track_history_py[daughter_id].push_back((Float_t) myFourVector.py());
0362   m_true_daughter_track_history_pz[daughter_id].push_back((Float_t) myFourVector.pz());
0363   m_true_daughter_track_history_pE[daughter_id].push_back((Float_t) myFourVector.e());
0364   m_true_daughter_track_history_pT[daughter_id].push_back((Float_t) myFourVector.perp());
0365 }
0366
0367 int KFParticle_truthAndDetTools::getHepMCInfo(PHCompositeNode *topNode, TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id)
0368 {
0369   // Make dummy particle for null pointers and missing nodes
0370   HepMC::GenParticle *dummyParticle = new HepMC::GenParticle();
0371   HepMC::FourVector dummyFourVector(0, 0, 0, 0);
0372   dummyParticle->set_momentum(dummyFourVector);
0373   dummyParticle->set_pdg_id(0);
0374   dummyParticle->set_generated_mass(0.);
0375
0376   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
0377   if (!dst_trackmap)
0378   {
0379     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
0380     gSystem->Exit(1);
0381     exit(1);
0382   }
0383
0384   track = getTrack(daughter.Id(), dst_trackmap);
0385   g4particle = getTruthTrack(track, topNode);
0386
0387   bool isParticleValid = g4particle == nullptr ? false : true;
0388
0389   if (!isParticleValid)
0390   {
0391     std::cout << "KFParticle truth matching: this track is a ghost" << std::endl;
0392     fillHepMCBranch(dummyParticle, daughter_id);
0393     return 0;
0394   }
0395
0396   m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0397   if (!m_geneventmap)
0398   {
0399     std::cout << "KFParticle truth matching: Missing node PHHepMCGenEventMap" << std::endl;
0400     std::cout << "You will have no mother information" << std::endl;
0401     fillHepMCBranch(dummyParticle, daughter_id);
0402     return 0;
0403   }
0404
0405   m_genevt = m_geneventmap->get(1);
0406   if (!m_genevt)
0407   {
0408     std::cout << "KFParticle truth matching: Missing node PHHepMCGenEvent" << std::endl;
0409     std::cout << "You will have no mother information" << std::endl;
0410     fillHepMCBranch(dummyParticle, daughter_id);
0411     return 0;
0412   }
0413
0414   // Start by looking for our particle in the Geant record
0415   // Any decay that Geant4 handles will not be in the HepMC record
0416   // This can happen if you limit the decay volume in the generator
0417   if (g4particle->get_parent_id() != 0)
0418   {
0419     m_truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0420     if (!m_truthinfo)
0421     {
0422       std::cout << "KFParticle truth matching: G4TruthInfo does not exist" << std::endl;
0423     }
0424     while (g4particle->get_parent_id() != 0)
0425     {
0426       g4particle = m_truthinfo->GetParticle(g4particle->get_parent_id());
0427       fillGeant4Branch(g4particle, daughter_id);
0428     }
0429   }
0430
0431   HepMC::GenEvent *theEvent = m_genevt->getEvent();
0432   HepMC::GenParticle *prevParticle = nullptr;
0433
0434   // int forbiddenPDGIDs[] = {21, 22};  //Stop tracing history when we reach quarks, gluons and photons
0435   int forbiddenPDGIDs[] = {0};  // 20230921 - Request made to have gluon information to see about gluon-splitting
0436
0437   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin(); p != theEvent->particles_end(); ++p)
0438   {
0439     if (((*p)->barcode() == g4particle->get_barcode()))
0440     {
0441       prevParticle = *p;
0442       while (!prevParticle->is_beam())
0443       {
0444         bool breakOut = false;
0445         for (HepMC::GenVertex::particle_iterator mother = prevParticle->production_vertex()->particles_begin(HepMC::parents);
0446              mother != prevParticle->production_vertex()->particles_end(HepMC::parents); ++mother)
0447         {
0448           if (std::find(std::begin(forbiddenPDGIDs), std::end(forbiddenPDGIDs),
0449                         abs((*mother)->pdg_id())) != std::end(forbiddenPDGIDs))
0450           {
0451             breakOut = true;
0452             break;
0453           }
0454
0455           fillHepMCBranch((*mother), daughter_id);
0456           prevParticle = *mother;
0457         }
0458         if (breakOut)
0459         {
0460           break;
0461         }
0462       }
0463     }
0464   }
0465
0466   return 0;
0467 }  // End of function
0468
0469 void KFParticle_truthAndDetTools::initializeCaloBranches(TTree *m_tree, int daughter_id, const std::string &daughter_number)
0470 {
0471   m_tree->Branch((daughter_number + "_EMCAL_DeltaPhi").c_str(), &detector_emcal_deltaphi[daughter_id], (daughter_number + "_EMCAL_DeltaPhi/F").c_str());
0472   m_tree->Branch((daughter_number + "_EMCAL_DeltaEta").c_str(), &detector_emcal_deltaeta[daughter_id], (daughter_number + "_EMCAL_DeltaEta/F").c_str());
0473   m_tree->Branch((daughter_number + "_EMCAL_DeltaZ").c_str(), &detector_emcal_deltaz[daughter_id], (daughter_number + "_EMCAL_DeltaZ/F").c_str());
0474   m_tree->Branch((daughter_number + "_EMCAL_energy_3x3").c_str(), &detector_emcal_energy_3x3[daughter_id], (daughter_number + "_EMCAL_energy_3x3/F").c_str());
0475   m_tree->Branch((daughter_number + "_EMCAL_energy_5x5").c_str(), &detector_emcal_energy_5x5[daughter_id], (daughter_number + "_EMCAL_energy_5x5/F").c_str());
0476   m_tree->Branch((daughter_number + "_EMCAL_energy_cluster").c_str(), &detector_emcal_cluster_energy[daughter_id], (daughter_number + "_EMCAL_energy_cluster/F").c_str());
0477   m_tree->Branch((daughter_number + "_IHCAL_DeltaPhi").c_str(), &detector_ihcal_deltaphi[daughter_id], (daughter_number + "_IHCAL_DeltaPhi/F").c_str());
0478   m_tree->Branch((daughter_number + "_IHCAL_DeltaEta").c_str(), &detector_ihcal_deltaeta[daughter_id], (daughter_number + "_IHCAL_DeltaEta/F").c_str());
0479   m_tree->Branch((daughter_number + "_IHCAL_energy_3x3").c_str(), &detector_ihcal_energy_3x3[daughter_id], (daughter_number + "_IHCAL_energy_3x3/F").c_str());
0480   m_tree->Branch((daughter_number + "_IHCAL_energy_5x5").c_str(), &detector_ihcal_energy_5x5[daughter_id], (daughter_number + "_IHCAL_energy_5x5/F").c_str());
0481   m_tree->Branch((daughter_number + "_IHCAL_energy_cluster").c_str(), &detector_ihcal_cluster_energy[daughter_id], (daughter_number + "_IHCAL_energy_cluster/F").c_str());
0482   m_tree->Branch((daughter_number + "_OHCAL_DeltaPhi").c_str(), &detector_ohcal_deltaphi[daughter_id], (daughter_number + "_OHCAL_DeltaEta/F").c_str());
0483   m_tree->Branch((daughter_number + "_OHCAL_DeltaEta").c_str(), &detector_ohcal_deltaeta[daughter_id], (daughter_number + "_OHCAL_DeltaEta/F").c_str());
0484   m_tree->Branch((daughter_number + "_OHCAL_energy_3x3").c_str(), &detector_ohcal_energy_3x3[daughter_id], (daughter_number + "_OHCAL_energy_3x3/F").c_str());
0485   m_tree->Branch((daughter_number + "_OHCAL_energy_5x5").c_str(), &detector_ohcal_energy_5x5[daughter_id], (daughter_number + "_OHCAL_energy_5x5/F").c_str());
0486   m_tree->Branch((daughter_number + "_OHCAL_energy_cluster").c_str(), &detector_ohcal_cluster_energy[daughter_id], (daughter_number + "_OHCAL_energy_cluster/F").c_str());
0487
0488   //Cluster Shape Studies
0489   if(m_get_detailed_calorimetry){
0490     m_tree->Branch((daughter_number + "_EMCAL_5x5Cell_binPhi").c_str(), &detector_emcal_5x5Cell_Phi[daughter_id]);
0491     m_tree->Branch((daughter_number + "_EMCAL_5x5Cell_binEta").c_str(), &detector_emcal_5x5Cell_Eta[daughter_id]);
0492     m_tree->Branch((daughter_number + "_EMCAL_5x5Cell_E").c_str(), &detector_emcal_5x5Cell_E[daughter_id]);
0493     m_tree->Branch((daughter_number + "_EMCAL_nTowers").c_str(), &detector_emcal_ntowers[daughter_id], (daughter_number + "_EMCAL_nTowers/i").c_str());
0494     m_tree->Branch((daughter_number + "_EMCAL_Chi2").c_str(), &detector_emcal_chi2[daughter_id], (daughter_number + "_EMCAL_Chi2/F").c_str());
0495   }
0496
0497 }
0498
0499 /*
0500 The following function matches tracks to calo clusters. As of 7/1/2025, this only extends to the EMCal. HCal matching is in development.
0501
0502 To run EMCal matching, DST_CALO files must be read into the Fun4All server.
0503 */
0504 void KFParticle_truthAndDetTools::fillCaloBranch(PHCompositeNode *topNode,
0505                                                  TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id, bool &isTrackEMCalmatch)
0506 {
0507   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
0508   if (!dst_trackmap)
0509   {
0510     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
0511     gSystem->Exit(1);
0512     exit(1);
0513   }
0514
0515   track = getTrack(daughter.Id(), dst_trackmap);
0516
0517   if (!clustersEM)
0518   {
0519     clustersEM = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0520     if (!clustersEM)
0521     {
0522       std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTERINFO_CEMC" << std::endl;
0523     }
0524   }
0525   if (!EMCalGeo)
0526   {
0527     EMCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0528     if (!EMCalGeo)
0529     {
0530       std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_CEMC" << std::endl;
0531       // return Fun4AllReturnCodes::ABORTEVENT;
0532     }
0533   }
0534     if(!_towersEM)
0535   {
0536     _towersEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0537     if(!_towersEM)
0538     {
0539       std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_CEMC" << std::endl;
0540       //return Fun4AllReturnCodes::ABORTEVENT;
0541     }
0542   }
0543   if (!clustersIH)
0544   {
0545     clustersIH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALIN");
0546     // if (!clustersIH)
0547     // {
0548     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTER_HCALIN" << std::endl;
0549     //   //return Fun4AllReturnCodes::ABORTEVENT;
0550     // }
0551   }
0552   if (!IHCalGeo)
0553   {
0554     IHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0555     // if(!IHCalGeo)
0556     // {
0557     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_HCALIN" << std::endl;
0558     //   //return Fun4AllReturnCodes::ABORTEVENT;
0559     // }
0560   }
0561   // if(!_towersIH)
0562   // {
0563   //   _towersIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0564   //   if(!_towersIH)
0565   //   {
0566   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_HCALIN" << std::endl;
0567   //     //return Fun4AllReturnCodes::ABORTEVENT;
0568   //   }
0569   // }
0570   if (!clustersOH)
0571   {
0572     clustersOH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALOUT");
0573     // if (!clustersOH)
0574     // {
0575     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTER_HCALOUT" << std::endl;
0576     //   //return Fun4AllReturnCodes::ABORTEVENT;
0577     // }
0578   }
0579   if (!OHCalGeo)
0580   {
0581     OHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0582     // if(!OHCalGeo)
0583     // {
0584     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_HCALOUT" << std::endl;
0585     //   //return Fun4AllReturnCodes::ABORTEVENT;
0586     // }
0587   }
0588   // if(!_towersOH)
0589   // {
0590   //   _towersOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0591   //   if(!_towersOH)
0592   //   {
0593   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_HCALOUT" << std::endl;
0594   //     //return Fun4AllReturnCodes::ABORTEVENT;
0595   //   }
0596   // }
0597
0598   // Radii for track projections
0599   double caloRadiusEMCal;
0600   caloRadiusEMCal = m_emcal_radius_user;  //Use function set_emcal_radius_user(float set_variable) to set this in your Fun4All macro
0601   // double caloRadiusIHCal;
0602   // double caloRadiusOHCal;
0603   // caloRadiusEMCal = 100.70;
0604   // caloRadiusEMCal = EMCalGeo->get_radius(); //This requires DST_CALOFITTING
0605   // caloRadiusOHCal = OHCalGeo->get_radius();
0606   // caloRadiusIHCal = IHCalGeo->get_radius();
0607
0608   // Create variables and containers etc.
0609   bool is_match;
0610   int index;
0611   float radius_scale;
0612   RawCluster *cluster;
0613
0614   // EMCAL******************************************************
0615   //  std::cout << "Starting EMCAL-track matching!" << std::endl;
0616
0617   SvtxTrackState *thisState = nullptr;
0618   thisState = track->get_state(caloRadiusEMCal);
0619   // thisState->identify();
0620   // std::cout << "size states " << (size_t)track->size_states() << std::endl;
0621   // for (auto state_iter = track->begin_states();
0622   //       state_iter != track->end_states();
0623   //       ++state_iter)
0624   //   {
0625   //     SvtxTrackState* tstate = state_iter->second;
0626   //     tstate->identify();
0627   //     std::cout<< "radius " << sqrt((tstate->get_x())*(tstate->get_x()) + (tstate->get_y())*(tstate->get_y())) << std::endl;
0628   //   }
0629
0630   // assert(thisState);
0631   float _track_phi_emc = std::numeric_limits<float>::quiet_NaN();
0632   float _track_eta_emc = std::numeric_limits<float>::quiet_NaN();
0633   float _track_x_emc = std::numeric_limits<float>::quiet_NaN();
0634   float _track_y_emc = std::numeric_limits<float>::quiet_NaN();
0635   float _track_z_emc = std::numeric_limits<float>::quiet_NaN();
0636
0637   // EMCal variables and vectors
0638   float _emcal_phi = std::numeric_limits<float>::quiet_NaN();
0639   float _emcal_eta = std::numeric_limits<float>::quiet_NaN();
0640   float _emcal_x = std::numeric_limits<float>::quiet_NaN();
0641   float _emcal_y = std::numeric_limits<float>::quiet_NaN();
0642   float _emcal_z = std::numeric_limits<float>::quiet_NaN();
0643   radius_scale = std::numeric_limits<float>::quiet_NaN();
0644   float _emcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0645   float _emcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0646   float _emcal_clusE = std::numeric_limits<float>::quiet_NaN();
0647   std::vector<float> v_emcal_phi;
0648   std::vector<float> v_emcal_eta;
0649   std::vector<float> v_emcal_x;
0650   std::vector<float> v_emcal_y;
0651   std::vector<float> v_emcal_z;
0652   std::vector<float> v_emcal_dphi;
0653   std::vector<float> v_emcal_deta;
0654   std::vector<float> v_emcal_dr;
0655   std::vector<float> v_emcal_dz;
0656   std::vector<float> v_emcal_3x3;
0657   std::vector<float> v_emcal_5x5;
0658   std::vector<float> v_emcal_clusE;
0659   //Detailed Calo Info
0660   float _emcal_nTowers = std::numeric_limits<float>::quiet_NaN();
0661   float _emcal_chi2 = std::numeric_limits<float>::quiet_NaN();
0662   std::vector<float> v_emcal_nTowers;
0663   std::vector<float> v_emcal_chi2;
0664   RawClusterDefs::keytype _clus_key;
0665   std::vector<RawClusterDefs::keytype> v_clus_key;
0666
0667   // Set variables for matching
0668   is_match = false;
0669   index = -1;
0670   // int ijk = 0; // nothing is being done with this variable in the end
0671
0672   //clustersEM->identify();
0673
0674   if (thisState != nullptr)
0675   {
0676     _track_x_emc = thisState->get_x();
0677     _track_y_emc = thisState->get_y();
0678     _track_z_emc = thisState->get_z();
0679     _track_phi_emc = std::atan2(_track_y_emc, _track_x_emc);
0680     _track_eta_emc = std::asinh(_track_z_emc / std::sqrt((_track_x_emc * _track_x_emc) + (_track_y_emc * _track_y_emc)));
0681
0682     // Create objects, containers, iterators for clusters
0683     cluster = nullptr;
0684     RawClusterContainer::Range begin_end_EMC = clustersEM->getClusters();
0685     RawClusterContainer::Iterator clusIter_EMC;
0686
0687     // Loop over the EMCal clusters
0688     for (clusIter_EMC = begin_end_EMC.first; clusIter_EMC != begin_end_EMC.second; ++clusIter_EMC)
0689     {
0690       // Minimum energy cut
0691       cluster = clusIter_EMC->second;
0692       float cluster_energy = cluster->get_energy();
0693
0694       if (cluster_energy < m_emcal_e_low_cut)
0695       {
0696         // ijk++;
0697         continue;
0698       }
0699
0700       // Get cluster information
0701       _emcal_x = cluster->get_x();
0702       _emcal_y = cluster->get_y();
0703       _emcal_z = cluster->get_z();
0704       radius_scale = m_emcal_radius_user / std::sqrt((_emcal_x * _emcal_x) + (_emcal_y * _emcal_y));
0705       _emcal_x *= radius_scale;
0706       _emcal_y *= radius_scale;
0707       _emcal_z *= radius_scale;
0708       _emcal_phi = std::atan2(_emcal_y, _emcal_x);
0709       _emcal_eta = std::asinh(_emcal_z / std::sqrt((_emcal_x * _emcal_x) + (_emcal_y * _emcal_y)));
0710       // _emcal_3x3 = get_e3x3(cluster, _towersEM, 0); //0 for emcal
0711       // _emcal_5x5 = get_e5x5(cluster, _towersEM, 0); //0 for emcal
0712       _emcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0713       _emcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0714       _emcal_clusE = cluster_energy;
0715       //Detailed Calo Info
0716       _emcal_nTowers = cluster->getNTowers();
0717       _emcal_chi2 = cluster->get_chi2();
0718       _clus_key = cluster->get_id();
0719
0720
0721
0722   //     virtual RawClusterDefs::keytype get_id() const
0723   // {
0724   //   PHOOL_VIRTUAL_WARN("get_id()");
0725   //   return 0;
0726   // }
0727
0728       // Variables to determine potential matches
0729       float dphi = PiRange(_track_phi_emc - _emcal_phi);
0730       float dz = _track_z_emc - _emcal_z;
0731       float deta = _track_eta_emc - _emcal_eta;
0732       float tmparg = caloRadiusEMCal * dphi;
0733       float dr = std::sqrt((tmparg * tmparg) + (dz * dz));  // sqrt((R*dphi)^2 + (dz)^2
0734       // float dr = sqrt((dphi*dphi + deta*deta)); //previous version
0735
0736       // Requirements for a possible match
0737       if (dz > m_dz_cut_high || dz < m_dz_cut_low)
0738       {
0739         continue;
0740       }
0741       if (dphi > m_dphi_cut_high || dphi < m_dphi_cut_low)
0742       {
0743         continue;
0744       }
0745
0746       // std::cout << "**********DELTA INFORMATION************" << std::endl;
0747       // std::cout << "dphi = " << dphi << std::endl;
0748       // std::cout << "deta = " << deta << std::endl;
0749       // std::cout << "dz =   " << dz <<   std::endl;
0750       // std::cout << "dr =   " << dr <<   std::endl;
0751       // std::cout << "****************************************" << std::endl;
0752       // Add potential match's information to vectors
0753       v_emcal_phi.push_back(_emcal_phi);
0754       v_emcal_eta.push_back(_emcal_eta);
0755       v_emcal_x.push_back(_emcal_x);
0756       v_emcal_y.push_back(_emcal_y);
0757       v_emcal_z.push_back(_emcal_z);
0758       v_emcal_dphi.push_back(dphi);
0759       v_emcal_dz.push_back(dz);
0760       v_emcal_deta.push_back(deta);
0761       v_emcal_dr.push_back(dr);
0762       v_emcal_3x3.push_back(_emcal_3x3);
0763       v_emcal_5x5.push_back(_emcal_5x5);
0764       v_emcal_clusE.push_back(_emcal_clusE);
0765       //Detailed Calo Info
0766       v_emcal_nTowers.push_back(_emcal_nTowers);
0767       v_emcal_chi2.push_back(_emcal_chi2);
0768       v_clus_key.push_back(_clus_key);
0769
0770       is_match = true;
0771     }
0772
0773     // Find the closest match from all potential matches
0774     if (is_match == true)
0775     {
0776       float tmp = 99999;
0777       for (long unsigned int i = 0; i < v_emcal_dr.size(); i++)
0778       {
0779         if (v_emcal_dr[i] < tmp)
0780         {
0781           index = i;
0782           tmp = v_emcal_dr[i];
0783         }
0784       }
0785     }
0786
0787   }
0788
0789   /*
0790   // Print out statements
0791   if (index != -1)
0792   {
0793     std::cout << "matched tracks!!!" << std::endl;
0794     std::cout << "emcal x = " << v_emcal_x[index] << " , y = " << v_emcal_y[index] << " , z = " << v_emcal_z[index] << " , phi = " << v_emcal_phi[index] << " , eta = " << v_emcal_eta[index] << std::endl;
0795     std::cout << "track projected x = " << _track_x_emc << " , y = " << _track_y_emc << " , z = " << _track_z_emc << " , phi = " << _track_phi_emc << " , eta = " << _track_eta_emc << std::endl;
0796     std::cout << "track px = " << track->get_px() << " , py = " << track->get_py() << " , pz = " << track->get_pz() << " , pt = " << track->get_pt() << " , p = " << track->get_p() << " , charge = " << track->get_charge() << std::endl;
0797   }
0798   */
0799
0800   // Save values to the branches!
0801   if (index == -1)
0802   {
0803     detector_emcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0804     detector_emcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0805     detector_emcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0806     detector_emcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0807     detector_emcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0808     detector_emcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0809     if(m_get_detailed_calorimetry)
0810     {
0811     detector_emcal_ntowers[daughter_id] = std::numeric_limits<unsigned int>::quiet_NaN();
0812     detector_emcal_chi2[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0813     detector_emcal_5x5Cell_Phi[daughter_id].push_back(std::numeric_limits<unsigned int>::quiet_NaN());
0814     detector_emcal_5x5Cell_Eta[daughter_id].push_back(std::numeric_limits<unsigned int>::quiet_NaN());
0815     detector_emcal_5x5Cell_E[daughter_id].push_back(std::numeric_limits<float>::quiet_NaN());
0816     }
0817     isTrackEMCalmatch = false;
0818   }
0819   else
0820   {
0821     detector_emcal_deltaphi[daughter_id] = v_emcal_dphi[index];
0822     detector_emcal_deltaeta[daughter_id] = v_emcal_deta[index];
0823     detector_emcal_deltaz[daughter_id] = v_emcal_dz[index];
0824     detector_emcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0825     detector_emcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0826     detector_emcal_cluster_energy[daughter_id] = v_emcal_clusE[index];
0827     if(m_get_detailed_calorimetry)
0828     {
0829     detector_emcal_ntowers[daughter_id] = v_emcal_nTowers[index];
0830     detector_emcal_chi2[daughter_id] = v_emcal_chi2[index];
0831     KFParticle_truthAndDetTools::Get5x5CellInfo(v_clus_key[index], daughter_id); //Adds the tower info to the vectors
0832     }
0833     isTrackEMCalmatch = true;
0834   }
0835
0836
0837   // HCAL*******************************************************
0838
0839   // INNER
0840   /*
0841   std::cout << "Starting IHCAL-track matching!" << std::endl;
0842
0843   // Track projection
0844   thisState = nullptr;
0845   thisState = track->get_state(caloRadiusIHCal);
0846   float _track_phi_ihc = std::numeric_limits<float>::quiet_NaN();
0847   float _track_eta_ihc = std::numeric_limits<float>::quiet_NaN();
0848   float _track_x_ihc = std::numeric_limits<float>::quiet_NaN();
0849   float _track_y_ihc = std::numeric_limits<float>::quiet_NaN();
0850   float _track_z_ihc = std::numeric_limits<float>::quiet_NaN();
0851
0852   // IHCal variables and vectors
0853   float _ihcal_phi = std::numeric_limits<float>::quiet_NaN();
0854   float _ihcal_eta = std::numeric_limits<float>::quiet_NaN();
0855   float _ihcal_x = std::numeric_limits<float>::quiet_NaN();
0856   float _ihcal_y = std::numeric_limits<float>::quiet_NaN();
0857   float _ihcal_z = std::numeric_limits<float>::quiet_NaN();
0858   float _ihcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0859   //float _ihcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0860   float _ihcal_clusE = std::numeric_limits<float>::quiet_NaN();
0861   radius_scale = std::numeric_limits<float>::quiet_NaN();
0862   std::vector<float> v_ihcal_phi;
0863   std::vector<float> v_ihcal_eta;
0864   std::vector<float> v_ihcal_x;
0865   std::vector<float> v_ihcal_y;
0866   std::vector<float> v_ihcal_z;
0867   std::vector<float> v_ihcal_dphi;
0868   std::vector<float> v_ihcal_deta;
0869   std::vector<float> v_ihcal_dr;
0870   std::vector<float> v_ihcal_3x3;
0871   //std::vector<float> v_ihcal_5x5;
0872   std::vector<float> v_ihcal_clusE;
0873
0874   if(thisState != nullptr){
0875
0876   // Reset variables for matching
0877   is_match = false;
0878   index = -1;
0879
0880   _track_phi_ihc = atan2(thisState->get_y(), thisState->get_x());
0881   _track_eta_ihc = asinh(thisState->get_z()/sqrt(thisState->get_x()*thisState->get_x() + thisState->get_y()*thisState->get_y()));
0882   _track_x_ihc = thisState->get_x();
0883   _track_y_ihc = thisState->get_y();
0884   _track_z_ihc = thisState->get_z();
0885
0886
0887   // Create objects, containers, iterators for clusters
0888   cluster = nullptr;
0889   RawClusterContainer::Range begin_end_IHC = clustersIH->getClusters();
0890   RawClusterContainer::Iterator clusIter_IHC;
0891
0892   // Loop over the IHCal clusters
0893   for (clusIter_IHC = begin_end_IHC.first; clusIter_IHC != begin_end_IHC.second; ++clusIter_IHC)
0894   {
0895
0896   // Minimum energy cut
0897   cluster = clusIter_IHC->second;
0898   if(cluster->get_energy() < m_ihcal_e_low_cut)
0899   {
0900   continue;
0901   }
0902
0903   std::cout << "Found a cluster about threshhold energy!" << std::endl;
0904
0905   // Get cluster information
0906   _ihcal_phi = atan2(cluster->get_y(), cluster->get_x());
0907   _ihcal_eta = asinh(cluster->get_z()/sqrt(cluster->get_x()*cluster->get_x() + cluster->get_y()*cluster->get_y()));
0908   _ihcal_x = cluster->get_x();
0909   _ihcal_y = cluster->get_y();
0910   radius_scale = m_ihcal_radius_user / sqrt(_ihcal_x*_ihcal_x+_ihcal_y*_ihcal_y);
0911   _ihcal_z = radius_scale*cluster->get_z();
0912   // _ihcal_3x3 = get_e3x3(cluster, _towersIH, 0); //1 for ihcal
0913   // _ihcal_5x5 = get_e5x5(cluster, _towersIH, 0); //1 for ihcal
0914   _ihcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0915   _ihcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0916   _ihcal_clusE = cluster->get_energy();
0917
0918
0919   // Variables to determine potential matches
0920   float dphi = abs(PiRange(_track_phi_ihc - _ihcal_phi));
0921   float dz = abs(_track_z_ihc - _ihcal_z);
0922   float deta = abs(_ihcal_eta - _track_eta_ihc);
0923   float dr = sqrt((dphi*dphi + deta*deta));
0924
0925   // Requirements for a possible match
0926   if(dphi<m_dphi_cut && dz<m_dz_cut)
0927   {
0928   //Add potential match's information to vectors
0929   v_ihcal_phi.push_back(_ihcal_phi);
0930   v_ihcal_eta.push_back(_ihcal_eta);
0931   v_ihcal_x.push_back(_ihcal_x);
0932   v_ihcal_y.push_back(_ihcal_y);
0933   v_ihcal_z.push_back(_ihcal_z);
0934   v_ihcal_dphi.push_back(dphi);
0935   v_ihcal_deta.push_back(deta);
0936   v_ihcal_dr.push_back(dr);
0937   v_ihcal_3x3.push_back(_ihcal_3x3);
0938   v_ihcal_clusE.push_back(_ihcal_clusE);
0939
0940   is_match = true;
0941   }
0942   }
0943
0944   // Find the closest match from all potential matches
0945   if (is_match == true)
0946   {
0947   float tmp = 99999;
0948   for(long unsigned int i = 0; i < v_ihcal_dr.size(); i++){
0949   if(v_ihcal_dr[i] < tmp){
0950   index = i;
0951   tmp = v_ihcal_dr[i];
0952   }
0953   }
0954   }
0955   }
0956
0957   // Print out statihents
0958   if(index != -1){
0959   std::cout<<"matched tracks!!!"<<std::endl;
0960   std::cout<<"ihcal x = "<<v_ihcal_x[index]<<" , y = "<<v_ihcal_y[index]<<" , z = "<<v_ihcal_z[index]<<" , phi = "<<v_ihcal_phi[index]<<" , eta = "<<v_ihcal_eta[index]<<std::endl;
0961   std::cout<<"track projected x = "<<_track_x_ihc<<" , y = "<<_track_y_ihc<<" , z = "<<_track_z_ihc<<" , phi = "<<_track_phi_ihc<<" , eta = "<<_track_eta_ihc<<std::endl;
0962   std::cout<<"track px = "<<track->get_px()<<" , py = "<<track->get_py()<<" , pz = "<<track->get_pz()<<" , pt = "<<track->get_pt()<<" , p = "<<track->get_p()<<" , charge = "<<track->get_charge()<<std::endl;
0963   }
0964   */
0965
0966   // Save values to the branches!
0967   // if(index == -1){
0968   detector_ihcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0969   detector_ihcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0970   detector_ihcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0971   detector_ihcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0972   detector_ihcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0973   // }
0974   // else{
0975   // detector_ihcal_deltaphi[daughter_id] = v_ihcal_dphi[index];
0976   // detector_ihcal_deltaeta[daughter_id] = v_ihcal_deta[index];
0977   // detector_ihcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0978   // detector_ihcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0979   // detector_ihcal_cluster_energy[daughter_id] = v_ihcal_clusE[index];
0980   // }
0981
0982   // std::cout << "IHCAL CLLUSTERS MATCHED TO TRACK" << std::endl;
0983
0984   // OUTER
0985
0986   // std::cout << "Starting OHCAL-track matching!" << std::endl;
0987
0988   /*
0989
0990   // Track projection
0991   thisState = nullptr;
0992   thisState = track->get_state(caloRadiusOHCal);
0993   float _track_phi_ohc = std::numeric_limits<float>::quiet_NaN();
0994   float _track_eta_ohc = std::numeric_limits<float>::quiet_NaN();
0995   float _track_x_ohc = std::numeric_limits<float>::quiet_NaN();
0996   float _track_y_ohc = std::numeric_limits<float>::quiet_NaN();
0997   float _track_z_ohc = std::numeric_limits<float>::quiet_NaN();
0998
0999   // OHCal variables and vectors
1000   float _ohcal_phi = std::numeric_limits<float>::quiet_NaN();
1001   float _ohcal_eta = std::numeric_limits<float>::quiet_NaN();
1002   float _ohcal_x = std::numeric_limits<float>::quiet_NaN();
1003   float _ohcal_y = std::numeric_limits<float>::quiet_NaN();
1004   float _ohcal_z = std::numeric_limits<float>::quiet_NaN();
1005   float _ohcal_3x3 = std::numeric_limits<float>::quiet_NaN();
1006   //float _ohcal_5x5 = std::numeric_limits<float>::quiet_NaN();
1007   float _ohcal_clusE = std::numeric_limits<float>::quiet_NaN();
1008   radius_scale = std::numeric_limits<float>::quiet_NaN();
1009   std::vector<float> v_ohcal_phi;
1010   std::vector<float> v_ohcal_eta;
1011   std::vector<float> v_ohcal_x;
1012   std::vector<float> v_ohcal_y;
1013   std::vector<float> v_ohcal_z;
1014   std::vector<float> v_ohcal_dphi;
1015   std::vector<float> v_ohcal_deta;
1016   std::vector<float> v_ohcal_dr;
1017   std::vector<float> v_ohcal_3x3;
1018   //std::vector<float> v_ohcal_5x5;
1019   std::vector<float> v_ohcal_clusE;
1020
1021   // Reset variables for matching
1022   is_match = false;
1023   index = -1;
1024
1025   if(thisState != nullptr)
1026   {
1027   _track_phi_ohc = atan2(thisState->get_y(), thisState->get_x());
1028   _track_eta_ohc = asinh(thisState->get_z()/sqrt(thisState->get_x()*thisState->get_x() + thisState->get_y()*thisState->get_y()));
1029   _track_x_ohc = thisState->get_x();
1030   _track_y_ohc = thisState->get_y();
1031   _track_z_ohc = thisState->get_z();
1032
1033
1034   // Create objects, containers, iterators for clusters
1035   cluster = nullptr;
1036   RawClusterContainer::Range begin_end_OHC = clustersOH->getClusters();
1037   RawClusterContainer::Iterator clusIter_OHC;
1038
1039   // Loop over the OHCal clusters
1040   for (clusIter_OHC = begin_end_OHC.first; clusIter_OHC != begin_end_OHC.second; ++clusIter_OHC)
1041   {
1042
1043   // Minimum energy cut
1044   cluster = clusIter_OHC->second;
1045   if(cluster->get_energy() < m_ohcal_e_low_cut)
1046   {
1047   continue;
1048   }
1049
1050   // Get cluster information
1051   _ohcal_phi = atan2(cluster->get_y(), cluster->get_x());
1052   _ohcal_eta = asinh(cluster->get_z()/sqrt(cluster->get_x()*cluster->get_x() + cluster->get_y()*cluster->get_y()));
1053   _ohcal_x = cluster->get_x();
1054   _ohcal_y = cluster->get_y();
1055   radius_scale = m_ohcal_radius_user / sqrt(_ohcal_x*_ohcal_x+_ohcal_y*_ohcal_y);
1056   _ohcal_z = radius_scale*cluster->get_z();
1057   // _ohcal_3x3 = get_e3x3(cluster, _towersOH, 0); //2 for ohcal
1058   // _ohcal_5x5 = get_e5x5(cluster, _towersOH, 0); //2 for ohcal
1059   _ohcal_3x3 = std::numeric_limits<float>::quiet_NaN();
1060   _ohcal_5x5 = std::numeric_limits<float>::quiet_NaN();
1061   _ohcal_clusE = cluster->get_energy();
1062
1063   // Variables to determine potential matches
1064   float dphi = abs(PiRange(_track_phi_ohc - _ohcal_phi));
1065   float dz = abs(_track_z_ohc - _ohcal_z);
1066   float deta = abs(_ohcal_eta - _track_eta_ohc);
1067   float dr = sqrt((dphi*dphi + deta*deta));
1068
1069   // Requirohents for a possible match
1070   if(dphi<m_dphi_cut && dz<m_dz_cut)
1071   {
1072   //Add potential match's information to vectors
1073   v_ohcal_phi.push_back(_ohcal_phi);
1074   v_ohcal_eta.push_back(_ohcal_eta);
1075   v_ohcal_x.push_back(_ohcal_x);
1076   v_ohcal_y.push_back(_ohcal_y);
1077   v_ohcal_z.push_back(_ohcal_z);
1078   v_ohcal_dphi.push_back(dphi);
1079   v_ohcal_deta.push_back(deta);
1080   v_ohcal_dr.push_back(dr);
1081   v_ohcal_3x3.push_back(_ohcal_3x3);
1082   v_ohcal_clusE.push_back(_ohcal_clusE);
1083
1084   is_match = true;
1085   }
1086   }
1087
1088
1089   // Find the closest match from all potential matches
1090   if (is_match == true)
1091   {
1092   float tmp = 99999;
1093   for(long unsigned int i = 0; i < v_ohcal_dr.size(); i++){
1094   if(v_ohcal_dr[i] < tmp){
1095   index = i;
1096   tmp = v_ohcal_dr[i];
1097   }
1098   }
1099   }
1100   }
1101
1102   std::cout << "match identified" << std::endl;
1103
1104   // Print out statements
1105   if(index != -1){
1106   std::cout<<"matched tracks!!!"<<std::endl;
1107   std::cout<<"ohcal x = "<<v_ohcal_x[index]<<" , y = "<<v_ohcal_y[index]<<" , z = "<<v_ohcal_z[index]<<" , phi = "<<v_ohcal_phi[index]<<" , eta = "<<v_ohcal_eta[index]<<std::endl;
1108   std::cout<<"track projected x = "<<_track_x_ohc<<" , y = "<<_track_y_ohc<<" , z = "<<_track_z_ohc<<" , phi = "<<_track_phi_ohc<<" , eta = "<<_track_eta_ohc<<std::endl;
1109   std::cout<<"track px = "<<track->get_px()<<" , py = "<<track->get_py()<<" , pz = "<<track->get_pz()<<" , pt = "<<track->get_pt()<<" , p = "<<track->get_p()<<" , charge = "<<track->get_charge()<<std::endl;
1110   }
1111   */
1112
1113   // Save values to the branches!
1114   // if(index == -1){
1115   detector_ohcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1116   detector_ohcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1117   detector_ohcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1118   detector_ohcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1119   detector_ohcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1120   // }
1121   // else{
1122   // detector_ohcal_deltaphi[daughter_id] = v_ohcal_dphi[index];
1123   // detector_ohcal_deltaeta[daughter_id] = v_ohcal_deta[index];
1124   // detector_ohcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1125   // detector_ohcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1126   // detector_ohcal_cluster_energy[daughter_id] = v_ohcal_clusE[index];
1127   // }
1128   // std::cout << "OHCAL CLLUSTERS MATCHED TO TRACK" << std::endl;
1129 }
1130
1131
1132 void KFParticle_truthAndDetTools::Get5x5CellInfo(RawClusterDefs::keytype key_in, int daughter_id){
1133
1134   //Get cluster from cluster key
1135   RawCluster *clus = clustersEM->getCluster(key_in);
1136
1137   //Set up variables for looping over towers
1138   int center_phi = std::numeric_limits<int>::quiet_NaN();
1139   int center_eta = std::numeric_limits<int>::quiet_NaN();
1140   float center_energy = -9999;
1141   const RawCluster::TowerConstRange begin_end_Towers = clus->get_towers();
1142   RawCluster::TowerConstIterator towersIter;
1143
1144   //Loop over the towers in the cluster to find the cluster center
1145   for (towersIter = begin_end_Towers.first; towersIter != begin_end_Towers.second; ++towersIter){
1146     if(towersIter->second > center_energy){
1147       //Setup for phi and eta
1148       RawTowerGeom *tower_geom = EMCalGeo->get_tower_geometry(towersIter->first);
1149       //Get energy, phi, and eta of center tower
1150       center_energy = towersIter->second;
1151       center_phi = EMCalGeo->get_phibin(tower_geom->get_phi());
1152       center_eta = EMCalGeo->get_etabin(tower_geom->get_eta());
1153     }
1154   } // Now we have the center tower and can grab the 5x5 grid info
1155
1156   //Vectors used in next for loops
1157   std::vector<unsigned int> v_phi_tmp;
1158   std::vector<unsigned int> v_eta_tmp;
1159
1160   //Get phi bins, and account for wrap-around case
1161   for(int i = -2; i<=2; i++){
1162     int phibin = center_phi+i;
1163     //Wraparound
1164     if(phibin < 0 || phibin > 255){
1165       if(phibin == -1)  { phibin = 255;}
1166       if(phibin == -2)  { phibin = 254;}
1167       if(phibin == 256) { phibin = 0;}
1168       if(phibin == 257) { phibin = 1;}
1169       if(phibin == 258) { phibin = 2;}
1170     }
1171     v_phi_tmp.push_back(phibin);
1172   }
1173
1174   //Get eta bins, and account for end of acceptance
1175   for(int i = -2; i<=2; i++){
1176     int etabin = center_eta+i;
1177     //End of eta acceptance case
1178     if(etabin < 0 || etabin > 95){ continue; }
1179     v_eta_tmp.push_back(etabin);
1180     // else{ v_eta_tmp.push_back(etabin); }
1181   }
1182
1183   //Vector to hold keys and energies
1184   std::vector<unsigned int> v_tower_phi = {};
1185   std::vector<unsigned int> v_tower_eta = {};
1186   std::vector<float> v_tower_energy = {};
1187
1188   //Ranges for for-loops
1189   int phi_range = v_phi_tmp.size();
1190   int eta_range = v_eta_tmp.size();
1191
1192   //for loop to fill vectors which will be stored in the KFP output nTuple
1193   for(int i = 0; i<phi_range; i++){
1194     //Get ith phi bin
1195     const unsigned int iphi = v_phi_tmp[i];
1196
1197     for(int j=0; j<eta_range; j++){
1198
1199       //Get jth eta bin
1200       const unsigned int jeta = v_eta_tmp[j];
1201       //Get key at phi and eta bin
1202       unsigned int towerinfo_key = TowerInfoDefs::encode_emcal(jeta, iphi);
1203       //Get energy
1204       TowerInfo *tower = _towersEM->get_tower_at_key(towerinfo_key);
1205       float energytmp = tower->get_energy();
1206
1207       //Store in vectors
1208       detector_emcal_5x5Cell_Phi[daughter_id].push_back(iphi);
1209       detector_emcal_5x5Cell_Eta[daughter_id].push_back(jeta);
1210       detector_emcal_5x5Cell_E[daughter_id].push_back(energytmp);
1211     }
1212   }
1213
1214 }
1215
1216 void KFParticle_truthAndDetTools::initializeDetectorBranches(TTree *m_tree, int daughter_id, const std::string &daughter_number)
1217 {
1218
1219   if(m_get_detailed_tracking){
1220     m_tree->Branch((daughter_number + "_residual_x").c_str(), &residual_x[daughter_id]);
1221     m_tree->Branch((daughter_number + "_residual_y").c_str(), &residual_y[daughter_id]);
1222     m_tree->Branch((daughter_number + "_residual_z").c_str(), &residual_z[daughter_id]);
1223     m_tree->Branch((daughter_number + "_layer").c_str(), &detector_layer[daughter_id]);
1224   }
1225
1226   for (auto const &subdetector : Use)
1227   {
1228     if (subdetector.second)
1229     {
1230       initializeSubDetectorBranches(m_tree, subdetector.first, daughter_id, daughter_number);
1231     }
1232   }
1233 }
1234 void KFParticle_truthAndDetTools::initializeSubDetectorBranches(TTree *m_tree, const std::string &detectorName, int daughter_id, const std::string &daughter_number)
1235 {
1236   if (detectorName == "MVTX")
1237   {
1238     if(m_get_detailed_tracking){
1239     m_tree->Branch((daughter_number + "_" + detectorName + "_staveID").c_str(), &mvtx_staveID[daughter_id]);
1240     m_tree->Branch((daughter_number + "_" + detectorName + "_chipID").c_str(), &mvtx_chipID[daughter_id]);
1241     }
1242     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_MVTX[daughter_id]);
1243     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_MVTX[daughter_id]);
1244   }
1245   if (detectorName == "INTT")
1246   {
1247     if(m_get_detailed_tracking){
1248     m_tree->Branch((daughter_number + "_" + detectorName + "_ladderZID").c_str(), &intt_ladderZID[daughter_id]);
1249     m_tree->Branch((daughter_number + "_" + detectorName + "_ladderPhiID").c_str(), &intt_ladderPhiID[daughter_id]);
1250     }
1251     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_INTT[daughter_id]);
1252     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_INTT[daughter_id]);
1253   }
1254   if (detectorName == "TPC")
1255   {
1256     if(m_get_detailed_tracking){
1257     m_tree->Branch((daughter_number + "_" + detectorName + "_sectorID").c_str(), &tpc_sectorID[daughter_id]);
1258     m_tree->Branch((daughter_number + "_" + detectorName + "_side").c_str(), &tpc_side[daughter_id]);
1259     }
1260     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_TPC[daughter_id]);
1261     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_TPC[daughter_id]);
1262   }
1263   if (detectorName == "TPOT")
1264   {
1265     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_TPOT[daughter_id]);
1266     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_TPOT[daughter_id]);
1267   }
1268 }
1269
1270 void KFParticle_truthAndDetTools::fillDetectorBranch(PHCompositeNode *topNode,
1271                                                      TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id)
1272 {
1273   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
1274   if (!dst_trackmap)
1275   {
1276     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
1277   }
1278
1279   std::string geoName = "ActsGeometry";
1280   geometry = findNode::getClass<ActsGeometry>(topNode, geoName);
1281   if (!geometry)
1282   {
1283     std::cout << "KFParticle detector info: " << geoName << " does not exist" << std::endl;
1284   }
1285
1286   dst_clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1287   if (!dst_clustermap)
1288   {
1289     std::cout << "KFParticle detector info: TRKR_CLUSTER does not exist" << std::endl;
1290   }
1291
1292   track = getTrack(daughter.Id(), dst_trackmap);
1293   detector_nHits_MVTX[daughter_id] = 0;
1294   detector_nHits_INTT[daughter_id] = 0;
1295   detector_nHits_TPC[daughter_id] = 0;
1296   detector_nHits_TPOT[daughter_id] = 0;
1297
1298   TrackSeed *silseed = track->get_silicon_seed();
1299   TrackSeed *tpcseed = track->get_tpc_seed();
1300
1301   if (silseed)
1302   {
1303     for (auto cluster_iter = silseed->begin_cluster_keys(); cluster_iter != silseed->end_cluster_keys(); ++cluster_iter)
1304     {
1305       const auto &cluster_key = *cluster_iter;
1306       const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
1307
1308       detector_layer[daughter_id].push_back(TrkrDefs::getLayer(cluster_key));
1309
1310       unsigned int staveId;
1311       unsigned int chipId;
1312       unsigned int ladderZId;
1313       unsigned int ladderPhiId;
1314       unsigned int sectorId;
1315       unsigned int side;
1316       staveId = chipId = ladderZId = ladderPhiId = sectorId = side = std::numeric_limits<unsigned int>::quiet_NaN();
1317
1318       if (Use["MVTX"] && trackerID == TrkrDefs::mvtxId)
1319       {
1320         staveId = MvtxDefs::getStaveId(cluster_key);
1321         chipId = MvtxDefs::getChipId(cluster_key);
1322         ++detector_nHits_MVTX[daughter_id];
1323       }
1324       else if (Use["INTT"] && trackerID == TrkrDefs::inttId)
1325       {
1326         ladderZId = InttDefs::getLadderZId(cluster_key);
1327         ladderPhiId = InttDefs::getLadderPhiId(cluster_key);
1328         ++detector_nHits_INTT[daughter_id];
1329       }
1330
1331       if(m_get_detailed_tracking){
1332       mvtx_staveID[daughter_id].push_back(staveId);
1333       mvtx_chipID[daughter_id].push_back(chipId);
1334       intt_ladderZID[daughter_id].push_back(ladderZId);
1335       intt_ladderPhiID[daughter_id].push_back(ladderPhiId);
1336       tpc_sectorID[daughter_id].push_back(sectorId);
1337       tpc_side[daughter_id].push_back(side);
1338       }
1339     }
1340   }
1341
1342   if (tpcseed)
1343   {
1344     for (auto cluster_iter = tpcseed->begin_cluster_keys(); cluster_iter != tpcseed->end_cluster_keys(); ++cluster_iter)
1345     {
1346       const auto &cluster_key = *cluster_iter;
1347       const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
1348
1349       detector_layer[daughter_id].push_back(TrkrDefs::getLayer(cluster_key));
1350
1351       unsigned int staveId;
1352       unsigned int chipId;
1353       unsigned int ladderZId;
1354       unsigned int ladderPhiId;
1355       unsigned int sectorId;
1356       unsigned int side;
1357       staveId = chipId = ladderZId = ladderPhiId = sectorId = side = std::numeric_limits<unsigned int>::quiet_NaN();
1358
1359       if (Use["TPC"] && trackerID == TrkrDefs::tpcId)
1360       {
1361         sectorId = TpcDefs::getSectorId(cluster_key);
1362         side = TpcDefs::getSide(cluster_key);
1363         ++detector_nHits_TPC[daughter_id];
1364       }
1365       else if (Use["TPOT"] && trackerID == TrkrDefs::micromegasId)
1366       {
1367         ++detector_nHits_TPOT[daughter_id];
1368       }
1369
1370       if(m_get_detailed_tracking){
1371       mvtx_staveID[daughter_id].push_back(staveId);
1372       mvtx_chipID[daughter_id].push_back(chipId);
1373       intt_ladderZID[daughter_id].push_back(ladderZId);
1374       intt_ladderPhiID[daughter_id].push_back(ladderPhiId);
1375       tpc_sectorID[daughter_id].push_back(sectorId);
1376       tpc_side[daughter_id].push_back(side);
1377       }
1378     }
1379   }
1380
1381   for (auto state_iter = track->begin_states();
1382        state_iter != track->end_states();
1383        ++state_iter)
1384   {
1385     SvtxTrackState *tstate = state_iter->second;
1386     if (tstate->get_pathlength() != 0)  // The first track state is an extrapolation so has no cluster
1387     {
1388       auto stateckey = tstate->get_cluskey();
1389       TrkrCluster *cluster = dst_clustermap->findCluster(stateckey);
1390       if (!cluster)
1391       {
1392     // do not have associated cluster, could be track states projected to calo system
1393         continue;
1394       }
1395       auto global = geometry->getGlobalPosition(stateckey, cluster);
1396
1397       if(m_get_detailed_tracking){
1398       residual_x[daughter_id].push_back(global.x() - tstate->get_x());
1399       residual_y[daughter_id].push_back(global.y() - tstate->get_y());
1400       residual_z[daughter_id].push_back(global.z() - tstate->get_z());
1401       }
1402
1403       uint8_t id = TrkrDefs::getTrkrId(stateckey);
1404
1405       switch (id)
1406       {
1407       case TrkrDefs::mvtxId:
1408         ++detector_nStates_MVTX[daughter_id];
1409         break;
1410       case TrkrDefs::inttId:
1411         ++detector_nStates_INTT[daughter_id];
1412         break;
1413       case TrkrDefs::tpcId:
1414         ++detector_nStates_TPC[daughter_id];
1415         break;
1416       case TrkrDefs::micromegasId:
1417         ++detector_nStates_TPOT[daughter_id];
1418         break;
1419       default:
1420         //std::cout << "Cluster key doesnt match a tracking system, could be related with projected track state to calorimeter system" << std::endl;
1421         break;
1422       }
1423     }
1424   }
1425 }
1426
1427 int KFParticle_truthAndDetTools::getPVID(PHCompositeNode *topNode, const KFParticle &kfpvertex)
1428 {
1429   if (m_dont_use_global_vertex_truth)
1430   {
1431     if (m_use_mbd_vertex_truth)
1432     {
1433       dst_mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
1434       if (dst_mbdvertexmap)
1435       {
1436         MbdVertex *m_dst_vertex = dst_mbdvertexmap->get(kfpvertex.Id());
1437         return m_dst_vertex->get_beam_crossing();
1438       }
1439
1440       std::cout << "KFParticle vertex matching: " << m_vtx_map_node_name_nTuple << " does not exist" << std::endl;
1441     }
1442     else
1443     {
1444       dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vtx_map_node_name_nTuple);
1445       if (dst_vertexmap)
1446       {
1447         SvtxVertex *m_dst_vertex = dst_vertexmap->get(kfpvertex.Id());
1448         return m_dst_vertex->get_beam_crossing();
1449       }
1450
1451       std::cout << "KFParticle vertex matching: " << m_vtx_map_node_name_nTuple << " does not exist" << std::endl;
1452     }
1453   }
1454   else
1455   {
1456     auto globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1457     if (!globalvertexmap)
1458     {
1459       return -100;
1460     }
1461
1462     GlobalVertex *gvertex = globalvertexmap->get(kfpvertex.Id());
1463     return gvertex->get_beam_crossing();
1464   }
1465
1466   return -100;
1467 }
1468
1469 void KFParticle_truthAndDetTools::allPVInfo(PHCompositeNode *topNode,
1470                                             TTree * /*m_tree*/,
1471                                             const KFParticle &motherParticle,
1472                                             std::vector<KFParticle> daughters,
1473                                             std::vector<KFParticle> intermediates)
1474 {
1475   KFParticle_Tools kfpTupleTools;
1476   kfpTupleTools.set_dont_use_global_vertex(m_dont_use_global_vertex_truth);
1477   std::vector<KFParticle> primaryVertices = kfpTupleTools.makeAllPrimaryVertices(topNode, m_vtx_map_node_name_nTuple);
1478
1479   for (auto &primaryVertice : primaryVertices)
1480   {
1481     allPV_x.push_back(primaryVertice.GetX());
1482     allPV_y.push_back(primaryVertice.GetY());
1483     allPV_z.push_back(primaryVertice.GetZ());
1484
1485     allPV_mother_IP.push_back(motherParticle.GetDistanceFromVertex(primaryVertice));
1486     allPV_mother_IPchi2.push_back(motherParticle.GetDeviationFromVertex(primaryVertice));
1487
1488     for (unsigned int j = 0; j < daughters.size(); ++j)
1489     {
1490       allPV_daughter_IP[j].push_back(daughters[j].GetDistanceFromVertex(primaryVertice));
1491       allPV_daughter_IPchi2[j].push_back(daughters[j].GetDeviationFromVertex(primaryVertice));
1492     }
1493
1494     for (unsigned int j = 0; j < intermediates.size(); ++j)
1495     {
1496       allPV_intermediates_IP[j].push_back(intermediates[j].GetDistanceFromVertex(primaryVertice));
1497       allPV_intermediates_IPchi2[j].push_back(intermediates[j].GetDeviationFromVertex(primaryVertice));
1498     }
1499   }
1500 }
1501
1502 void KFParticle_truthAndDetTools::clearVectors()
1503 {
1504   for (int i = 0; i < m_num_tracks_nTuple; ++i)
1505   {
1506     // Truth vectors
1507     m_true_daughter_track_history_PDG_ID[i].clear();
1508     m_true_daughter_track_history_PDG_mass[i].clear();
1509     m_true_daughter_track_history_px[i].clear();
1510     m_true_daughter_track_history_py[i].clear();
1511     m_true_daughter_track_history_pz[i].clear();
1512     m_true_daughter_track_history_pE[i].clear();
1513     m_true_daughter_track_history_pT[i].clear();
1514
1515     // Detector vectors
1516     residual_x[i].clear();
1517     residual_y[i].clear();
1518     residual_z[i].clear();
1519     detector_layer[i].clear();
1520     mvtx_staveID[i].clear();
1521     mvtx_chipID[i].clear();
1522     intt_ladderZID[i].clear();
1523     intt_ladderPhiID[i].clear();
1524     tpc_sectorID[i].clear();
1525     tpc_side[i].clear();
1526
1527     detector_nStates_MVTX[i] = 0;
1528     detector_nStates_INTT[i] = 0;
1529     detector_nStates_TPC[i] = 0;
1530     detector_nStates_TPOT[i] = 0;
1531
1532     // PV vectors
1533     allPV_daughter_IP[i].clear();
1534     allPV_daughter_IPchi2[i].clear();
1535
1536     //Detailed Calo
1537     if(m_get_detailed_calorimetry){
1538     detector_emcal_5x5Cell_Phi[i].clear();
1539     detector_emcal_5x5Cell_Eta[i].clear();
1540     detector_emcal_5x5Cell_E[i].clear();
1541     }
1542
1543   }
1544
1545   allPV_x.clear();
1546   allPV_y.clear();
1547   allPV_z.clear();
1548   allPV_z.clear();
1549
1550   allPV_mother_IP.clear();
1551   allPV_mother_IPchi2.clear();
1552
1553   for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
1554   {
1555     allPV_intermediates_IP[i].clear();
1556     allPV_intermediates_IPchi2[i].clear();
1557   }
1558 }