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
0489 /*
0490 The following function matches tracks to calo clusters. As of 7/1/2025, this only extends to the EMCal. HCal matching is in development.
0491
0492 To run EMCal matching, DST_CALO files must be read into the Fun4All server.
0493 */
0494 void KFParticle_truthAndDetTools::fillCaloBranch(PHCompositeNode *topNode,
0495                                                  TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id, bool &isTrackEMCalmatch)
0496 {
0497   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
0498   if (!dst_trackmap)
0499   {
0500     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
0501     gSystem->Exit(1);
0502     exit(1);
0503   }
0504
0505   track = getTrack(daughter.Id(), dst_trackmap);
0506
0507   if (!clustersEM)
0508   {
0509     clustersEM = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0510     if (!clustersEM)
0511     {
0512       std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTERINFO_CEMC" << std::endl;
0513     }
0514   }
0515   // if (!EMCalGeo)
0516   // {
0517   //   EMCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0518   //   if (!EMCalGeo)
0519   //   {
0520   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_CEMC" << std::endl;
0521   //     // return Fun4AllReturnCodes::ABORTEVENT;
0522   //   }
0523   // }
0524   // if(!_towersEM)
0525   // {
0526   //   _towersEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0527   //   if(!_towersEM)
0528   //   {
0529   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_CEMC" << std::endl;
0530   //     //return Fun4AllReturnCodes::ABORTEVENT;
0531   //   }
0532   // }
0533   if (!clustersIH)
0534   {
0535     clustersIH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALIN");
0536     // if (!clustersIH)
0537     // {
0538     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTER_HCALIN" << std::endl;
0539     //   //return Fun4AllReturnCodes::ABORTEVENT;
0540     // }
0541   }
0542   if (!IHCalGeo)
0543   {
0544     IHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0545     // if(!IHCalGeo)
0546     // {
0547     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_HCALIN" << std::endl;
0548     //   //return Fun4AllReturnCodes::ABORTEVENT;
0549     // }
0550   }
0551   // if(!_towersIH)
0552   // {
0553   //   _towersIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0554   //   if(!_towersIH)
0555   //   {
0556   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_HCALIN" << std::endl;
0557   //     //return Fun4AllReturnCodes::ABORTEVENT;
0558   //   }
0559   // }
0560   if (!clustersOH)
0561   {
0562     clustersOH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALOUT");
0563     // if (!clustersOH)
0564     // {
0565     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "CLUSTER_HCALOUT" << std::endl;
0566     //   //return Fun4AllReturnCodes::ABORTEVENT;
0567     // }
0568   }
0569   if (!OHCalGeo)
0570   {
0571     OHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0572     // if(!OHCalGeo)
0573     // {
0574     //   std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWERGEOM_HCALOUT" << std::endl;
0575     //   //return Fun4AllReturnCodes::ABORTEVENT;
0576     // }
0577   }
0578   // if(!_towersOH)
0579   // {
0580   //   _towersOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0581   //   if(!_towersOH)
0582   //   {
0583   //     std::cout << __FILE__ << "::" << __func__ << " : FATAL ERROR, cannot find cluster container " << "TOWER_CALIB_HCALOUT" << std::endl;
0584   //     //return Fun4AllReturnCodes::ABORTEVENT;
0585   //   }
0586   // }
0587
0588   // Radii for track projections
0589   double caloRadiusEMCal;
0590   caloRadiusEMCal = m_emcal_radius_user;  //Use function set_emcal_radius_user(float set_variable) to set this in your Fun4All macro
0591   // double caloRadiusIHCal;
0592   // double caloRadiusOHCal;
0593   // caloRadiusEMCal = 100.70;
0594   // caloRadiusEMCal = EMCalGeo->get_radius(); //This requires DST_CALOFITTING
0595   // caloRadiusOHCal = OHCalGeo->get_radius();
0596   // caloRadiusIHCal = IHCalGeo->get_radius();
0597
0598   // Create variables and containers etc.
0599   bool is_match;
0600   int index;
0601   float radius_scale;
0602   RawCluster *cluster;
0603
0604   // EMCAL******************************************************
0605   //  std::cout << "Starting EMCAL-track matching!" << std::endl;
0606
0607   SvtxTrackState *thisState = nullptr;
0608   thisState = track->get_state(caloRadiusEMCal);
0609   // thisState->identify();
0610   // std::cout << "size states " << (size_t)track->size_states() << std::endl;
0611   // for (auto state_iter = track->begin_states();
0612   //       state_iter != track->end_states();
0613   //       ++state_iter)
0614   //   {
0615   //     SvtxTrackState* tstate = state_iter->second;
0616   //     tstate->identify();
0617   //     std::cout<< "radius " << sqrt((tstate->get_x())*(tstate->get_x()) + (tstate->get_y())*(tstate->get_y())) << std::endl;
0618   //   }
0619
0620   // assert(thisState);
0621   float _track_phi_emc = std::numeric_limits<float>::quiet_NaN();
0622   float _track_eta_emc = std::numeric_limits<float>::quiet_NaN();
0623   float _track_x_emc = std::numeric_limits<float>::quiet_NaN();
0624   float _track_y_emc = std::numeric_limits<float>::quiet_NaN();
0625   float _track_z_emc = std::numeric_limits<float>::quiet_NaN();
0626
0627   // EMCal variables and vectors
0628   float _emcal_phi = std::numeric_limits<float>::quiet_NaN();
0629   float _emcal_eta = std::numeric_limits<float>::quiet_NaN();
0630   float _emcal_x = std::numeric_limits<float>::quiet_NaN();
0631   float _emcal_y = std::numeric_limits<float>::quiet_NaN();
0632   float _emcal_z = std::numeric_limits<float>::quiet_NaN();
0633   radius_scale = std::numeric_limits<float>::quiet_NaN();
0634   float _emcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0635   float _emcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0636   float _emcal_clusE = std::numeric_limits<float>::quiet_NaN();
0637   std::vector<float> v_emcal_phi;
0638   std::vector<float> v_emcal_eta;
0639   std::vector<float> v_emcal_x;
0640   std::vector<float> v_emcal_y;
0641   std::vector<float> v_emcal_z;
0642   std::vector<float> v_emcal_dphi;
0643   std::vector<float> v_emcal_deta;
0644   std::vector<float> v_emcal_dr;
0645   std::vector<float> v_emcal_dz;
0646   std::vector<float> v_emcal_3x3;
0647   std::vector<float> v_emcal_5x5;
0648   std::vector<float> v_emcal_clusE;
0649
0650   // Set variables for matching
0651   is_match = false;
0652   index = -1;
0653   // int ijk = 0; // nothing is being done with this variable in the end
0654
0655   //clustersEM->identify();
0656
0657   if (thisState != nullptr)
0658   {
0659     _track_x_emc = thisState->get_x();
0660     _track_y_emc = thisState->get_y();
0661     _track_z_emc = thisState->get_z();
0662     _track_phi_emc = std::atan2(_track_y_emc, _track_x_emc);
0663     _track_eta_emc = std::asinh(_track_z_emc / std::sqrt((_track_x_emc * _track_x_emc) + (_track_y_emc * _track_y_emc)));
0664
0665     // Create objects, containers, iterators for clusters
0666     cluster = nullptr;
0667     RawClusterContainer::Range begin_end_EMC = clustersEM->getClusters();
0668     RawClusterContainer::Iterator clusIter_EMC;
0669
0670     // Loop over the EMCal clusters
0671     for (clusIter_EMC = begin_end_EMC.first; clusIter_EMC != begin_end_EMC.second; ++clusIter_EMC)
0672     {
0673       // Minimum energy cut
0674       cluster = clusIter_EMC->second;
0675       float cluster_energy = cluster->get_energy();
0676
0677       if (cluster_energy < m_emcal_e_low_cut)
0678       {
0679         // ijk++;
0680         continue;
0681       }
0682
0683       // Get cluster information
0684       _emcal_x = cluster->get_x();
0685       _emcal_y = cluster->get_y();
0686       _emcal_z = cluster->get_z();
0687       radius_scale = m_emcal_radius_user / std::sqrt((_emcal_x * _emcal_x) + (_emcal_y * _emcal_y));
0688       _emcal_x *= radius_scale;
0689       _emcal_y *= radius_scale;
0690       _emcal_z *= radius_scale;
0691       _emcal_phi = std::atan2(_emcal_y, _emcal_y);
0692       _emcal_eta = std::asinh(_emcal_z / std::sqrt((_emcal_x * _emcal_x) + (_emcal_y * _emcal_y)));
0693       // _emcal_3x3 = get_e3x3(cluster, _towersEM, 0); //0 for emcal
0694       // _emcal_5x5 = get_e5x5(cluster, _towersEM, 0); //0 for emcal
0695       _emcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0696       _emcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0697       _emcal_clusE = cluster_energy;
0698
0699       // Variables to determine potential matches
0700       float dphi = PiRange(_track_phi_emc - _emcal_phi);
0701       float dz = _track_z_emc - _emcal_z;
0702       float deta = _track_eta_emc - _emcal_eta;
0703       float tmparg = caloRadiusEMCal * dphi;
0704       float dr = std::sqrt((tmparg * tmparg) + (dz * dz));  // sqrt((R*dphi)^2 + (dz)^2
0705       // float dr = sqrt((dphi*dphi + deta*deta)); //previous version
0706
0707       // Requirements for a possible match
0708       if (dz > m_dz_cut_high || dz < m_dz_cut_low)
0709       {
0710         continue;
0711       }
0712       if (dphi > m_dphi_cut_high || dphi < m_dphi_cut_low)
0713       {
0714         continue;
0715       }
0716
0717       // std::cout << "**********DELTA INFORMATION************" << std::endl;
0718       // std::cout << "dphi = " << dphi << std::endl;
0719       // std::cout << "deta = " << deta << std::endl;
0720       // std::cout << "dz =   " << dz <<   std::endl;
0721       // std::cout << "dr =   " << dr <<   std::endl;
0722       // std::cout << "****************************************" << std::endl;
0723       // Add potential match's information to vectors
0724       v_emcal_phi.push_back(_emcal_phi);
0725       v_emcal_eta.push_back(_emcal_eta);
0726       v_emcal_x.push_back(_emcal_x);
0727       v_emcal_y.push_back(_emcal_y);
0728       v_emcal_z.push_back(_emcal_z);
0729       v_emcal_dphi.push_back(dphi);
0730       v_emcal_dz.push_back(dz);
0731       v_emcal_deta.push_back(deta);
0732       v_emcal_dr.push_back(dr);
0733       v_emcal_3x3.push_back(_emcal_3x3);
0734       v_emcal_5x5.push_back(_emcal_5x5);
0735       v_emcal_clusE.push_back(_emcal_clusE);
0736       is_match = true;
0737       // ijk++;
0738     }
0739
0740     // Find the closest match from all potential matches
0741     if (is_match == true)
0742     {
0743       float tmp = 99999;
0744       for (long unsigned int i = 0; i < v_emcal_dr.size(); i++)
0745       {
0746         if (v_emcal_dr[i] < tmp)
0747         {
0748           index = i;
0749           tmp = v_emcal_dr[i];
0750         }
0751       }
0752     }
0753   }
0754
0755   /*
0756   // Print out statements
0757   if (index != -1)
0758   {
0759     std::cout << "matched tracks!!!" << std::endl;
0760     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;
0761     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;
0762     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;
0763   }
0764   */
0765
0766   // Save values to the branches!
0767   if (index == -1)
0768   {
0769     detector_emcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0770     detector_emcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0771     detector_emcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0772     detector_emcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0773     detector_emcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0774     detector_emcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0775     isTrackEMCalmatch = false;
0776   }
0777   else
0778   {
0779     detector_emcal_deltaphi[daughter_id] = v_emcal_dphi[index];
0780     detector_emcal_deltaeta[daughter_id] = v_emcal_deta[index];
0781     detector_emcal_deltaz[daughter_id] = v_emcal_dz[index];
0782     detector_emcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0783     detector_emcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0784     detector_emcal_cluster_energy[daughter_id] = v_emcal_clusE[index];
0785     isTrackEMCalmatch = true;
0786   }
0787
0788   // HCAL*******************************************************
0789
0790   // INNER
0791   /*
0792   std::cout << "Starting IHCAL-track matching!" << std::endl;
0793
0794   // Track projection
0795   thisState = nullptr;
0796   thisState = track->get_state(caloRadiusIHCal);
0797   float _track_phi_ihc = std::numeric_limits<float>::quiet_NaN();
0798   float _track_eta_ihc = std::numeric_limits<float>::quiet_NaN();
0799   float _track_x_ihc = std::numeric_limits<float>::quiet_NaN();
0800   float _track_y_ihc = std::numeric_limits<float>::quiet_NaN();
0801   float _track_z_ihc = std::numeric_limits<float>::quiet_NaN();
0802
0803   // IHCal variables and vectors
0804   float _ihcal_phi = std::numeric_limits<float>::quiet_NaN();
0805   float _ihcal_eta = std::numeric_limits<float>::quiet_NaN();
0806   float _ihcal_x = std::numeric_limits<float>::quiet_NaN();
0807   float _ihcal_y = std::numeric_limits<float>::quiet_NaN();
0808   float _ihcal_z = std::numeric_limits<float>::quiet_NaN();
0809   float _ihcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0810   //float _ihcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0811   float _ihcal_clusE = std::numeric_limits<float>::quiet_NaN();
0812   radius_scale = std::numeric_limits<float>::quiet_NaN();
0813   std::vector<float> v_ihcal_phi;
0814   std::vector<float> v_ihcal_eta;
0815   std::vector<float> v_ihcal_x;
0816   std::vector<float> v_ihcal_y;
0817   std::vector<float> v_ihcal_z;
0818   std::vector<float> v_ihcal_dphi;
0819   std::vector<float> v_ihcal_deta;
0820   std::vector<float> v_ihcal_dr;
0821   std::vector<float> v_ihcal_3x3;
0822   //std::vector<float> v_ihcal_5x5;
0823   std::vector<float> v_ihcal_clusE;
0824
0825   if(thisState != nullptr){
0826
0827   // Reset variables for matching
0828   is_match = false;
0829   index = -1;
0830
0831   _track_phi_ihc = atan2(thisState->get_y(), thisState->get_x());
0832   _track_eta_ihc = asinh(thisState->get_z()/sqrt(thisState->get_x()*thisState->get_x() + thisState->get_y()*thisState->get_y()));
0833   _track_x_ihc = thisState->get_x();
0834   _track_y_ihc = thisState->get_y();
0835   _track_z_ihc = thisState->get_z();
0836
0837
0838   // Create objects, containers, iterators for clusters
0839   cluster = nullptr;
0840   RawClusterContainer::Range begin_end_IHC = clustersIH->getClusters();
0841   RawClusterContainer::Iterator clusIter_IHC;
0842
0843   // Loop over the IHCal clusters
0844   for (clusIter_IHC = begin_end_IHC.first; clusIter_IHC != begin_end_IHC.second; ++clusIter_IHC)
0845   {
0846
0847   // Minimum energy cut
0848   cluster = clusIter_IHC->second;
0849   if(cluster->get_energy() < m_ihcal_e_low_cut)
0850   {
0851   continue;
0852   }
0853
0854   std::cout << "Found a cluster about threshhold energy!" << std::endl;
0855
0856   // Get cluster information
0857   _ihcal_phi = atan2(cluster->get_y(), cluster->get_x());
0858   _ihcal_eta = asinh(cluster->get_z()/sqrt(cluster->get_x()*cluster->get_x() + cluster->get_y()*cluster->get_y()));
0859   _ihcal_x = cluster->get_x();
0860   _ihcal_y = cluster->get_y();
0861   radius_scale = m_ihcal_radius_user / sqrt(_ihcal_x*_ihcal_x+_ihcal_y*_ihcal_y);
0862   _ihcal_z = radius_scale*cluster->get_z();
0863   // _ihcal_3x3 = get_e3x3(cluster, _towersIH, 0); //1 for ihcal
0864   // _ihcal_5x5 = get_e5x5(cluster, _towersIH, 0); //1 for ihcal
0865   _ihcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0866   _ihcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0867   _ihcal_clusE = cluster->get_energy();
0868
0869
0870   // Variables to determine potential matches
0871   float dphi = abs(PiRange(_track_phi_ihc - _ihcal_phi));
0872   float dz = abs(_track_z_ihc - _ihcal_z);
0873   float deta = abs(_ihcal_eta - _track_eta_ihc);
0874   float dr = sqrt((dphi*dphi + deta*deta));
0875
0876   // Requirements for a possible match
0877   if(dphi<m_dphi_cut && dz<m_dz_cut)
0878   {
0879   //Add potential match's information to vectors
0880   v_ihcal_phi.push_back(_ihcal_phi);
0881   v_ihcal_eta.push_back(_ihcal_eta);
0882   v_ihcal_x.push_back(_ihcal_x);
0883   v_ihcal_y.push_back(_ihcal_y);
0884   v_ihcal_z.push_back(_ihcal_z);
0885   v_ihcal_dphi.push_back(dphi);
0886   v_ihcal_deta.push_back(deta);
0887   v_ihcal_dr.push_back(dr);
0888   v_ihcal_3x3.push_back(_ihcal_3x3);
0889   v_ihcal_clusE.push_back(_ihcal_clusE);
0890
0891   is_match = true;
0892   }
0893   }
0894
0895   // Find the closest match from all potential matches
0896   if (is_match == true)
0897   {
0898   float tmp = 99999;
0899   for(long unsigned int i = 0; i < v_ihcal_dr.size(); i++){
0900   if(v_ihcal_dr[i] < tmp){
0901   index = i;
0902   tmp = v_ihcal_dr[i];
0903   }
0904   }
0905   }
0906   }
0907
0908   // Print out statihents
0909   if(index != -1){
0910   std::cout<<"matched tracks!!!"<<std::endl;
0911   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;
0912   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;
0913   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;
0914   }
0915   */
0916
0917   // Save values to the branches!
0918   // if(index == -1){
0919   detector_ihcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0920   detector_ihcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0921   detector_ihcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0922   detector_ihcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0923   detector_ihcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0924   // }
0925   // else{
0926   // detector_ihcal_deltaphi[daughter_id] = v_ihcal_dphi[index];
0927   // detector_ihcal_deltaeta[daughter_id] = v_ihcal_deta[index];
0928   // detector_ihcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0929   // detector_ihcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
0930   // detector_ihcal_cluster_energy[daughter_id] = v_ihcal_clusE[index];
0931   // }
0932
0933   // std::cout << "IHCAL CLLUSTERS MATCHED TO TRACK" << std::endl;
0934
0935   // OUTER
0936
0937   // std::cout << "Starting OHCAL-track matching!" << std::endl;
0938
0939   /*
0940
0941   // Track projection
0942   thisState = nullptr;
0943   thisState = track->get_state(caloRadiusOHCal);
0944   float _track_phi_ohc = std::numeric_limits<float>::quiet_NaN();
0945   float _track_eta_ohc = std::numeric_limits<float>::quiet_NaN();
0946   float _track_x_ohc = std::numeric_limits<float>::quiet_NaN();
0947   float _track_y_ohc = std::numeric_limits<float>::quiet_NaN();
0948   float _track_z_ohc = std::numeric_limits<float>::quiet_NaN();
0949
0950   // OHCal variables and vectors
0951   float _ohcal_phi = std::numeric_limits<float>::quiet_NaN();
0952   float _ohcal_eta = std::numeric_limits<float>::quiet_NaN();
0953   float _ohcal_x = std::numeric_limits<float>::quiet_NaN();
0954   float _ohcal_y = std::numeric_limits<float>::quiet_NaN();
0955   float _ohcal_z = std::numeric_limits<float>::quiet_NaN();
0956   float _ohcal_3x3 = std::numeric_limits<float>::quiet_NaN();
0957   //float _ohcal_5x5 = std::numeric_limits<float>::quiet_NaN();
0958   float _ohcal_clusE = std::numeric_limits<float>::quiet_NaN();
0959   radius_scale = std::numeric_limits<float>::quiet_NaN();
0960   std::vector<float> v_ohcal_phi;
0961   std::vector<float> v_ohcal_eta;
0962   std::vector<float> v_ohcal_x;
0963   std::vector<float> v_ohcal_y;
0964   std::vector<float> v_ohcal_z;
0965   std::vector<float> v_ohcal_dphi;
0966   std::vector<float> v_ohcal_deta;
0967   std::vector<float> v_ohcal_dr;
0968   std::vector<float> v_ohcal_3x3;
0969   //std::vector<float> v_ohcal_5x5;
0970   std::vector<float> v_ohcal_clusE;
0971
0972   // Reset variables for matching
0973   is_match = false;
0974   index = -1;
0975
0976   if(thisState != nullptr)
0977   {
0978   _track_phi_ohc = atan2(thisState->get_y(), thisState->get_x());
0979   _track_eta_ohc = asinh(thisState->get_z()/sqrt(thisState->get_x()*thisState->get_x() + thisState->get_y()*thisState->get_y()));
0980   _track_x_ohc = thisState->get_x();
0981   _track_y_ohc = thisState->get_y();
0982   _track_z_ohc = thisState->get_z();
0983
0984
0985   // Create objects, containers, iterators for clusters
0986   cluster = nullptr;
0987   RawClusterContainer::Range begin_end_OHC = clustersOH->getClusters();
0988   RawClusterContainer::Iterator clusIter_OHC;
0989
0990   // Loop over the OHCal clusters
0991   for (clusIter_OHC = begin_end_OHC.first; clusIter_OHC != begin_end_OHC.second; ++clusIter_OHC)
0992   {
0993
0994   // Minimum energy cut
0995   cluster = clusIter_OHC->second;
0996   if(cluster->get_energy() < m_ohcal_e_low_cut)
0997   {
0998   continue;
0999   }
1000
1001   // Get cluster information
1002   _ohcal_phi = atan2(cluster->get_y(), cluster->get_x());
1003   _ohcal_eta = asinh(cluster->get_z()/sqrt(cluster->get_x()*cluster->get_x() + cluster->get_y()*cluster->get_y()));
1004   _ohcal_x = cluster->get_x();
1005   _ohcal_y = cluster->get_y();
1006   radius_scale = m_ohcal_radius_user / sqrt(_ohcal_x*_ohcal_x+_ohcal_y*_ohcal_y);
1007   _ohcal_z = radius_scale*cluster->get_z();
1008   // _ohcal_3x3 = get_e3x3(cluster, _towersOH, 0); //2 for ohcal
1009   // _ohcal_5x5 = get_e5x5(cluster, _towersOH, 0); //2 for ohcal
1010   _ohcal_3x3 = std::numeric_limits<float>::quiet_NaN();
1011   _ohcal_5x5 = std::numeric_limits<float>::quiet_NaN();
1012   _ohcal_clusE = cluster->get_energy();
1013
1014   // Variables to determine potential matches
1015   float dphi = abs(PiRange(_track_phi_ohc - _ohcal_phi));
1016   float dz = abs(_track_z_ohc - _ohcal_z);
1017   float deta = abs(_ohcal_eta - _track_eta_ohc);
1018   float dr = sqrt((dphi*dphi + deta*deta));
1019
1020   // Requirohents for a possible match
1021   if(dphi<m_dphi_cut && dz<m_dz_cut)
1022   {
1023   //Add potential match's information to vectors
1024   v_ohcal_phi.push_back(_ohcal_phi);
1025   v_ohcal_eta.push_back(_ohcal_eta);
1026   v_ohcal_x.push_back(_ohcal_x);
1027   v_ohcal_y.push_back(_ohcal_y);
1028   v_ohcal_z.push_back(_ohcal_z);
1029   v_ohcal_dphi.push_back(dphi);
1030   v_ohcal_deta.push_back(deta);
1031   v_ohcal_dr.push_back(dr);
1032   v_ohcal_3x3.push_back(_ohcal_3x3);
1033   v_ohcal_clusE.push_back(_ohcal_clusE);
1034
1035   is_match = true;
1036   }
1037   }
1038
1039
1040   // Find the closest match from all potential matches
1041   if (is_match == true)
1042   {
1043   float tmp = 99999;
1044   for(long unsigned int i = 0; i < v_ohcal_dr.size(); i++){
1045   if(v_ohcal_dr[i] < tmp){
1046   index = i;
1047   tmp = v_ohcal_dr[i];
1048   }
1049   }
1050   }
1051   }
1052
1053   std::cout << "match identified" << std::endl;
1054
1055   // Print out statements
1056   if(index != -1){
1057   std::cout<<"matched tracks!!!"<<std::endl;
1058   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;
1059   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;
1060   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;
1061   }
1062   */
1063
1064   // Save values to the branches!
1065   // if(index == -1){
1066   detector_ohcal_deltaphi[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1067   detector_ohcal_deltaeta[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1068   detector_ohcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1069   detector_ohcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1070   detector_ohcal_cluster_energy[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1071   // }
1072   // else{
1073   // detector_ohcal_deltaphi[daughter_id] = v_ohcal_dphi[index];
1074   // detector_ohcal_deltaeta[daughter_id] = v_ohcal_deta[index];
1075   // detector_ohcal_energy_3x3[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1076   // detector_ohcal_energy_5x5[daughter_id] = std::numeric_limits<float>::quiet_NaN();
1077   // detector_ohcal_cluster_energy[daughter_id] = v_ohcal_clusE[index];
1078   // }
1079   // std::cout << "OHCAL CLLUSTERS MATCHED TO TRACK" << std::endl;
1080 }
1081
1082 void KFParticle_truthAndDetTools::initializeDetectorBranches(TTree *m_tree, int daughter_id, const std::string &daughter_number)
1083 {
1084   m_tree->Branch((daughter_number + "_residual_x").c_str(), &residual_x[daughter_id]);
1085   m_tree->Branch((daughter_number + "_residual_y").c_str(), &residual_y[daughter_id]);
1086   m_tree->Branch((daughter_number + "_residual_z").c_str(), &residual_z[daughter_id]);
1087   m_tree->Branch((daughter_number + "_layer").c_str(), &detector_layer[daughter_id]);
1088
1089   for (auto const &subdetector : Use)
1090   {
1091     if (subdetector.second)
1092     {
1093       initializeSubDetectorBranches(m_tree, subdetector.first, daughter_id, daughter_number);
1094     }
1095   }
1096 }
1097
1098 void KFParticle_truthAndDetTools::initializeSubDetectorBranches(TTree *m_tree, const std::string &detectorName, int daughter_id, const std::string &daughter_number)
1099 {
1100   if (detectorName == "MVTX")
1101   {
1102     m_tree->Branch((daughter_number + "_" + detectorName + "_staveID").c_str(), &mvtx_staveID[daughter_id]);
1103     m_tree->Branch((daughter_number + "_" + detectorName + "_chipID").c_str(), &mvtx_chipID[daughter_id]);
1104     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_MVTX[daughter_id]);
1105     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_MVTX[daughter_id]);
1106   }
1107   if (detectorName == "INTT")
1108   {
1109     m_tree->Branch((daughter_number + "_" + detectorName + "_ladderZID").c_str(), &intt_ladderZID[daughter_id]);
1110     m_tree->Branch((daughter_number + "_" + detectorName + "_ladderPhiID").c_str(), &intt_ladderPhiID[daughter_id]);
1111     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_INTT[daughter_id]);
1112     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_INTT[daughter_id]);
1113   }
1114   if (detectorName == "TPC")
1115   {
1116     m_tree->Branch((daughter_number + "_" + detectorName + "_sectorID").c_str(), &tpc_sectorID[daughter_id]);
1117     m_tree->Branch((daughter_number + "_" + detectorName + "_side").c_str(), &tpc_side[daughter_id]);
1118     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_TPC[daughter_id]);
1119     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_TPC[daughter_id]);
1120   }
1121   if (detectorName == "TPOT")
1122   {
1123     m_tree->Branch((daughter_number + "_" + detectorName + "_nHits").c_str(), &detector_nHits_TPOT[daughter_id]);
1124     m_tree->Branch((daughter_number + "_" + detectorName + "_nStates").c_str(), &detector_nStates_TPOT[daughter_id]);
1125   }
1126 }
1127
1128 void KFParticle_truthAndDetTools::fillDetectorBranch(PHCompositeNode *topNode,
1129                                                      TTree * /*m_tree*/, const KFParticle &daughter, int daughter_id)
1130 {
1131   dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name_nTuple);
1132   if (!dst_trackmap)
1133   {
1134     std::cout << "KFParticle truth matching: " << m_trk_map_node_name_nTuple << " does not exist" << std::endl;
1135   }
1136
1137   std::string geoName = "ActsGeometry";
1138   geometry = findNode::getClass<ActsGeometry>(topNode, geoName);
1139   if (!geometry)
1140   {
1141     std::cout << "KFParticle detector info: " << geoName << " does not exist" << std::endl;
1142   }
1143
1144   dst_clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1145   if (!dst_clustermap)
1146   {
1147     std::cout << "KFParticle detector info: TRKR_CLUSTER does not exist" << std::endl;
1148   }
1149
1150   track = getTrack(daughter.Id(), dst_trackmap);
1151   detector_nHits_MVTX[daughter_id] = 0;
1152   detector_nHits_INTT[daughter_id] = 0;
1153   detector_nHits_TPC[daughter_id] = 0;
1154   detector_nHits_TPOT[daughter_id] = 0;
1155
1156   TrackSeed *silseed = track->get_silicon_seed();
1157   TrackSeed *tpcseed = track->get_tpc_seed();
1158
1159   if (silseed)
1160   {
1161     for (auto cluster_iter = silseed->begin_cluster_keys(); cluster_iter != silseed->end_cluster_keys(); ++cluster_iter)
1162     {
1163       const auto &cluster_key = *cluster_iter;
1164       const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
1165
1166       detector_layer[daughter_id].push_back(TrkrDefs::getLayer(cluster_key));
1167
1168       unsigned int staveId;
1169       unsigned int chipId;
1170       unsigned int ladderZId;
1171       unsigned int ladderPhiId;
1172       unsigned int sectorId;
1173       unsigned int side;
1174       staveId = chipId = ladderZId = ladderPhiId = sectorId = side = std::numeric_limits<unsigned int>::quiet_NaN();
1175
1176       if (Use["MVTX"] && trackerID == TrkrDefs::mvtxId)
1177       {
1178         staveId = MvtxDefs::getStaveId(cluster_key);
1179         chipId = MvtxDefs::getChipId(cluster_key);
1180         ++detector_nHits_MVTX[daughter_id];
1181       }
1182       else if (Use["INTT"] && trackerID == TrkrDefs::inttId)
1183       {
1184         ladderZId = InttDefs::getLadderZId(cluster_key);
1185         ladderPhiId = InttDefs::getLadderPhiId(cluster_key);
1186         ++detector_nHits_INTT[daughter_id];
1187       }
1188
1189       mvtx_staveID[daughter_id].push_back(staveId);
1190       mvtx_chipID[daughter_id].push_back(chipId);
1191       intt_ladderZID[daughter_id].push_back(ladderZId);
1192       intt_ladderPhiID[daughter_id].push_back(ladderPhiId);
1193       tpc_sectorID[daughter_id].push_back(sectorId);
1194       tpc_side[daughter_id].push_back(side);
1195     }
1196   }
1197
1198   if (tpcseed)
1199   {
1200     for (auto cluster_iter = tpcseed->begin_cluster_keys(); cluster_iter != tpcseed->end_cluster_keys(); ++cluster_iter)
1201     {
1202       const auto &cluster_key = *cluster_iter;
1203       const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
1204
1205       detector_layer[daughter_id].push_back(TrkrDefs::getLayer(cluster_key));
1206
1207       unsigned int staveId;
1208       unsigned int chipId;
1209       unsigned int ladderZId;
1210       unsigned int ladderPhiId;
1211       unsigned int sectorId;
1212       unsigned int side;
1213       staveId = chipId = ladderZId = ladderPhiId = sectorId = side = std::numeric_limits<unsigned int>::quiet_NaN();
1214
1215       if (Use["TPC"] && trackerID == TrkrDefs::tpcId)
1216       {
1217         sectorId = TpcDefs::getSectorId(cluster_key);
1218         side = TpcDefs::getSide(cluster_key);
1219         ++detector_nHits_TPC[daughter_id];
1220       }
1221       else if (Use["TPOT"] && trackerID == TrkrDefs::micromegasId)
1222       {
1223         ++detector_nHits_TPOT[daughter_id];
1224       }
1225
1226       mvtx_staveID[daughter_id].push_back(staveId);
1227       mvtx_chipID[daughter_id].push_back(chipId);
1228       intt_ladderZID[daughter_id].push_back(ladderZId);
1229       intt_ladderPhiID[daughter_id].push_back(ladderPhiId);
1230       tpc_sectorID[daughter_id].push_back(sectorId);
1231       tpc_side[daughter_id].push_back(side);
1232     }
1233   }
1234
1235   for (auto state_iter = track->begin_states();
1236        state_iter != track->end_states();
1237        ++state_iter)
1238   {
1239     SvtxTrackState *tstate = state_iter->second;
1240     if (tstate->get_pathlength() != 0)  // The first track state is an extrapolation so has no cluster
1241     {
1242       auto stateckey = tstate->get_cluskey();
1243       TrkrCluster *cluster = dst_clustermap->findCluster(stateckey);
1244       if (!cluster)
1245       {
1246     // do not have associated cluster, could be track states projected to calo system
1247         continue;
1248       }
1249       auto global = geometry->getGlobalPosition(stateckey, cluster);
1250
1251       residual_x[daughter_id].push_back(global.x() - tstate->get_x());
1252       residual_y[daughter_id].push_back(global.y() - tstate->get_y());
1253       residual_z[daughter_id].push_back(global.z() - tstate->get_z());
1254
1255       uint8_t id = TrkrDefs::getTrkrId(stateckey);
1256
1257       switch (id)
1258       {
1259       case TrkrDefs::mvtxId:
1260         ++detector_nStates_MVTX[daughter_id];
1261         break;
1262       case TrkrDefs::inttId:
1263         ++detector_nStates_INTT[daughter_id];
1264         break;
1265       case TrkrDefs::tpcId:
1266         ++detector_nStates_TPC[daughter_id];
1267         break;
1268       case TrkrDefs::micromegasId:
1269         ++detector_nStates_TPOT[daughter_id];
1270         break;
1271       default:
1272         //std::cout << "Cluster key doesnt match a tracking system, could be related with projected track state to calorimeter system" << std::endl;
1273         break;
1274       }
1275     }
1276   }
1277 }
1278
1279 int KFParticle_truthAndDetTools::getPVID(PHCompositeNode *topNode, const KFParticle &kfpvertex)
1280 {
1281   if (m_dont_use_global_vertex_truth)
1282   {
1283     if (m_use_mbd_vertex_truth)
1284     {
1285       dst_mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
1286       if (dst_mbdvertexmap)
1287       {
1288         MbdVertex *m_dst_vertex = dst_mbdvertexmap->get(kfpvertex.Id());
1289         return m_dst_vertex->get_beam_crossing();
1290       }
1291
1292       std::cout << "KFParticle vertex matching: " << m_vtx_map_node_name_nTuple << " does not exist" << std::endl;
1293     }
1294     else
1295     {
1296       dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vtx_map_node_name_nTuple);
1297       if (dst_vertexmap)
1298       {
1299         SvtxVertex *m_dst_vertex = dst_vertexmap->get(kfpvertex.Id());
1300         return m_dst_vertex->get_beam_crossing();
1301       }
1302
1303       std::cout << "KFParticle vertex matching: " << m_vtx_map_node_name_nTuple << " does not exist" << std::endl;
1304     }
1305   }
1306   else
1307   {
1308     auto globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1309     if (!globalvertexmap)
1310     {
1311       return -100;
1312     }
1313
1314     GlobalVertex *gvertex = globalvertexmap->get(kfpvertex.Id());
1315     return gvertex->get_beam_crossing();
1316   }
1317
1318   return -100;
1319 }
1320
1321 void KFParticle_truthAndDetTools::allPVInfo(PHCompositeNode *topNode,
1322                                             TTree * /*m_tree*/,
1323                                             const KFParticle &motherParticle,
1324                                             std::vector<KFParticle> daughters,
1325                                             std::vector<KFParticle> intermediates)
1326 {
1327   KFParticle_Tools kfpTupleTools;
1328   kfpTupleTools.set_dont_use_global_vertex(m_dont_use_global_vertex_truth);
1329   std::vector<KFParticle> primaryVertices = kfpTupleTools.makeAllPrimaryVertices(topNode, m_vtx_map_node_name_nTuple);
1330
1331   for (auto &primaryVertice : primaryVertices)
1332   {
1333     allPV_x.push_back(primaryVertice.GetX());
1334     allPV_y.push_back(primaryVertice.GetY());
1335     allPV_z.push_back(primaryVertice.GetZ());
1336
1337     allPV_mother_IP.push_back(motherParticle.GetDistanceFromVertex(primaryVertice));
1338     allPV_mother_IPchi2.push_back(motherParticle.GetDeviationFromVertex(primaryVertice));
1339
1340     for (unsigned int j = 0; j < daughters.size(); ++j)
1341     {
1342       allPV_daughter_IP[j].push_back(daughters[j].GetDistanceFromVertex(primaryVertice));
1343       allPV_daughter_IPchi2[j].push_back(daughters[j].GetDeviationFromVertex(primaryVertice));
1344     }
1345
1346     for (unsigned int j = 0; j < intermediates.size(); ++j)
1347     {
1348       allPV_intermediates_IP[j].push_back(intermediates[j].GetDistanceFromVertex(primaryVertice));
1349       allPV_intermediates_IPchi2[j].push_back(intermediates[j].GetDeviationFromVertex(primaryVertice));
1350     }
1351   }
1352 }
1353
1354 void KFParticle_truthAndDetTools::clearVectors()
1355 {
1356   for (int i = 0; i < m_num_tracks_nTuple; ++i)
1357   {
1358     // Truth vectors
1359     m_true_daughter_track_history_PDG_ID[i].clear();
1360     m_true_daughter_track_history_PDG_mass[i].clear();
1361     m_true_daughter_track_history_px[i].clear();
1362     m_true_daughter_track_history_py[i].clear();
1363     m_true_daughter_track_history_pz[i].clear();
1364     m_true_daughter_track_history_pE[i].clear();
1365     m_true_daughter_track_history_pT[i].clear();
1366
1367     // Detector vectors
1368     residual_x[i].clear();
1369     residual_y[i].clear();
1370     residual_z[i].clear();
1371     detector_layer[i].clear();
1372     mvtx_staveID[i].clear();
1373     mvtx_chipID[i].clear();
1374     intt_ladderZID[i].clear();
1375     intt_ladderPhiID[i].clear();
1376     tpc_sectorID[i].clear();
1377     tpc_side[i].clear();
1378
1379     detector_nStates_MVTX[i] = 0;
1380     detector_nStates_INTT[i] = 0;
1381     detector_nStates_TPC[i] = 0;
1382     detector_nStates_TPOT[i] = 0;
1383
1384     // PV vectors
1385     allPV_daughter_IP[i].clear();
1386     allPV_daughter_IPchi2[i].clear();
1387   }
1388
1389   allPV_x.clear();
1390   allPV_y.clear();
1391   allPV_z.clear();
1392   allPV_z.clear();
1393
1394   allPV_mother_IP.clear();
1395   allPV_mother_IPchi2.clear();
1396
1397   for (int i = 0; i < m_num_intermediate_states_nTuple; ++i)
1398   {
1399     allPV_intermediates_IP[i].clear();
1400     allPV_intermediates_IPchi2[i].clear();
1401   }
1402 }