Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:20:17 

0001 /*
0002  * This file is part of KFParticle package
0003  * Copyright ( C ) 2007-2019 FIAS Frankfurt Institute for Advanced Studies
0004  *               2007-2019 Goethe University of Frankfurt
0005  *               2007-2019 Ivan Kisel <I.Kisel@compeng.uni-frankfurt.de>
0006  *               2007-2019 Maksym Zyzak
0007  *
0008  * KFParticle is free software: you can redistribute it and/or modify
0009  * it under the terms of the GNU General Public License as published by
0010  * the Free Software Foundation, either version 3 of the License, or
0011  * ( at your option ) any later version.
0012  *
0013  * KFParticle is distributed in the hope that it will be useful,
0014  * but WITHOUT ANY WARRANTY; without even the implied warranty of
0015  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0016  * GNU General Public License for more details.
0017  *
0018  * You should have received a copy of the GNU General Public License
0019  * along with this program.  If not, see <https://www.gnu.org/licenses/>.
0020  */
0021 
0022 /*****************/
0023 /* Cameron Dean  */
0024 /*   LANL 2020   */
0025 /* cdean@bnl.gov */
0026 /*****************/
0027 
0028 #include "KFParticle_Tools.h"
0029 
0030 #include <trackbase_historic/SvtxTrack.h>
0031 #include <trackbase_historic/SvtxTrackMap.h>
0032 #include <trackbase_historic/TrackAnalysisUtils.h>
0033 
0034 #include <trackbase/TrkrClusterContainer.h>
0035 #include <trackbase/TrkrCluster.h>
0036 #include <trackbase/TrkrDefs.h> 
0037 #include <g4detectors/PHG4TpcGeom.h>
0038 #include <g4detectors/PHG4TpcGeomContainer.h>
0039 
0040 #include <globalvertex/GlobalVertex.h>
0041 #include <globalvertex/GlobalVertexMap.h>
0042 #include <globalvertex/SvtxVertex.h>
0043 #include <globalvertex/SvtxVertexMap.h>
0044 
0045 #include <phool/getClass.h>
0046 
0047 #include <ffamodules/CDBInterface.h>
0048 
0049 // KFParticle stuff
0050 #include <KFParticle.h>
0051 #include <KFVertex.h>
0052 
0053 #include <Rtypes.h>
0054 #include <TDatabasePDG.h>
0055 #include <TMatrixD.h>
0056 #include "KFParticle_truthAndDetTools.h"
0057 
0058 #include <TFile.h>
0059 #include <TMatrixDfwd.h>  // for TMatrixD
0060 #include <TMatrixT.h>     // for TMatrixT, operator*
0061 
0062 #include <Eigen/Dense>
0063 
0064 #include <algorithm>  // for max, remove, minmax_el...
0065 #include <cmath>      // for sqrt, pow, M_PI
0066 #include <cstdlib>    // for abs, NULL
0067 #include <iostream>   // for operator<<, basic_ostream
0068 #include <iterator>   // for end
0069 #include <map>        // for _Rb_tree_iterator, map
0070 #include <memory>     // for allocator_traits<>::va...
0071 
0072 KFParticle_truthAndDetTools toolSet;
0073 
0074 /// KFParticle constructor
0075 KFParticle_Tools::KFParticle_Tools()
0076   : m_has_intermediates(false)
0077   , m_min_mass(0)
0078   , m_max_mass(0)
0079   , m_min_decayTime(-1 * FLT_MAX)
0080   , m_max_decayTime(FLT_MAX)
0081   , m_min_decayLength(-1 * FLT_MAX)
0082   , m_max_decayLength(FLT_MAX)
0083   , m_track_min_pt(0.)
0084   , m_track_max_pt(5e3)
0085   , m_track_ptchi2(FLT_MAX)
0086   , m_track_ip_xy(-100.)
0087   , m_track_ipchi2_xy(-1)
0088   , m_track_ip(-1.)
0089   , m_track_ipchi2(-1)
0090   , m_track_chi2ndof(100.)
0091   , m_nMVTXStates(3)
0092   , m_nINTTStates(1)
0093   , m_nTPCStates(20)
0094   , m_nTPOTStates(0)
0095   , m_comb_DCA_xy(1)
0096   , m_comb_DCA(0.5)
0097   , m_vertex_chi2ndof(20.)
0098   , m_fdchi2(0.)
0099   , m_dira_min(-1.01)
0100   , m_dira_max(1.01)
0101   , m_mother_pt(0.)
0102   , m_mother_ipchi2(FLT_MAX)
0103   , m_get_charge_conjugate(false)
0104   , m_extrapolateTracksToSV(true)
0105   , m_vtx_map_node_name("SvtxVertexMap")
0106   , m_trk_map_node_name("SvtxTrackMap")
0107   , m_dst_mbdvertexmap()
0108   , m_dst_mbdvertex()
0109   , m_dst_trackmap()
0110   , m_dst_track()
0111   , m_dst_vertexmap()
0112   , m_dst_vertex()
0113   , m_cluster_map()
0114   , m_geom_container()
0115 {
0116 }
0117 
0118 KFParticle KFParticle_Tools::makeVertex(PHCompositeNode * /*topNode*/)
0119 {
0120   float vtxX = m_use_mbd_vertex ? 0 : m_dst_vertex->get_x();
0121   float vtxY = m_use_mbd_vertex ? 0 : m_dst_vertex->get_y();
0122   float vtxZ = m_use_mbd_vertex ? m_dst_mbdvertex->get_z() : m_dst_vertex->get_z();
0123 
0124   float f_vertexParameters[6] = {vtxX, vtxY, vtxZ, 0, 0, 0};
0125 
0126   float f_vertexCovariance[21] = {0};
0127   unsigned int iterate = 0;
0128   if (m_use_mbd_vertex)
0129   {
0130     f_vertexCovariance[5] = m_dst_mbdvertex->get_z_err();
0131   }
0132   else
0133   {
0134     for (unsigned int i = 0; i < 3; ++i)
0135     {
0136       for (unsigned int j = 0; j <= i; ++j)
0137       {
0138         f_vertexCovariance[iterate] = m_dst_vertex->get_error(i, j);
0139         ++iterate;
0140       }
0141     }
0142   }
0143 
0144   KFParticle kfp_vertex;
0145   kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
0146   kfp_vertex.NDF() = m_use_mbd_vertex ? 0 : m_dst_vertex->get_ndof();
0147   kfp_vertex.Chi2() = m_use_mbd_vertex ? 0 : m_dst_vertex->get_chisq();
0148 
0149   return kfp_vertex;
0150 }
0151 
0152 std::vector<KFParticle> KFParticle_Tools::makeAllPrimaryVertices(PHCompositeNode *topNode, const std::string &vertexMapName)
0153 {
0154   std::string vtxMN;
0155 
0156   unsigned int vertexID = 0;
0157 
0158   if (vertexMapName.empty())
0159   {
0160     vtxMN = m_vtx_map_node_name;
0161   }
0162   else
0163   {
0164     vtxMN = vertexMapName;
0165   }
0166 
0167   std::vector<KFParticle> primaryVertices;
0168 
0169   if (m_use_mbd_vertex)
0170   {
0171     m_dst_mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0172   }
0173   else
0174   {
0175     m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, vtxMN);
0176   }
0177 
0178   if (m_dont_use_global_vertex)
0179   {
0180     if (m_use_mbd_vertex)
0181     {
0182       for (MbdVertexMap::ConstIter iter = m_dst_mbdvertexmap->begin(); iter != m_dst_mbdvertexmap->end(); ++iter)
0183       {
0184         m_dst_mbdvertex = iter->second;
0185         primaryVertices.push_back(makeVertex(topNode));
0186         primaryVertices[vertexID].SetId(iter->first);
0187         ++vertexID;
0188       }
0189     }
0190     else
0191     {
0192        for (SvtxVertexMap::ConstIter iter = m_dst_vertexmap->begin(); iter != m_dst_vertexmap->end(); ++iter)
0193       {
0194         m_dst_vertex = iter->second;
0195         primaryVertices.push_back(makeVertex(topNode));
0196         primaryVertices[vertexID].SetId(iter->first);
0197         ++vertexID;
0198       }
0199     }
0200 
0201     return primaryVertices;
0202   }
0203 
0204   m_dst_globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0205   if (!m_dst_globalvertexmap)
0206   {
0207     std::cout << "Can't continue in KFParticle_Tools::makeAllPrimaryVertices" << std::endl;
0208     return primaryVertices;
0209   }
0210 
0211   for (GlobalVertexMap::ConstIter iter = m_dst_globalvertexmap->begin(); iter != m_dst_globalvertexmap->end(); ++iter)
0212   {
0213     m_dst_globalvertex = iter->second;
0214     
0215     GlobalVertex::VTXTYPE whichVtx = m_use_mbd_vertex ? GlobalVertex::MBD : GlobalVertex::SVTX; 
0216 
0217     auto svtxiter = m_dst_globalvertex->find_vertexes(whichVtx);
0218     // check that it contains a track vertex
0219     if (svtxiter == m_dst_globalvertex->end_vertexes())
0220     {
0221       continue;
0222     }
0223 
0224     auto svtxvertexvector = svtxiter->second;
0225 
0226     for (auto &vertex : svtxvertexvector)
0227     {
0228       if (m_use_mbd_vertex)
0229       {
0230         m_dst_mbdvertex = m_dst_mbdvertexmap->find(vertex->get_id())->second;
0231       }
0232       else
0233       {
0234         m_dst_vertex = m_dst_vertexmap->find(vertex->get_id())->second;
0235       }
0236 
0237       primaryVertices.push_back(makeVertex(topNode));
0238       primaryVertices[vertexID].SetId(m_dst_globalvertex->get_id());
0239       ++vertexID;
0240     }
0241   }
0242 
0243   return primaryVertices;
0244 }
0245 
0246 KFParticle KFParticle_Tools::makeParticle(PHCompositeNode * /*topNode*/)  /// Return a KFPTrack from track vector and covariance matrix. No mass or vertex constraints
0247 {
0248   KFParticle kfp_particle;
0249 
0250   float f_trackParameters[6] = {m_dst_track->get_x(),
0251                                 m_dst_track->get_y(),
0252                                 m_dst_track->get_z(),
0253                                 m_dst_track->get_px(),
0254                                 m_dst_track->get_py(),
0255                                 m_dst_track->get_pz()};
0256 
0257   float f_trackCovariance[21];
0258   unsigned int iterate = 0;
0259   for (unsigned int i = 0; i < 6; ++i)
0260   {
0261     for (unsigned int j = 0; j <= i; ++j)
0262     {
0263       f_trackCovariance[iterate] = m_dst_track->get_error(i, j);
0264       ++iterate;
0265     }
0266   }
0267 
0268   kfp_particle.Create(f_trackParameters, f_trackCovariance, (Int_t) m_dst_track->get_charge(), -1);
0269   kfp_particle.NDF() = m_dst_track->get_ndf();
0270   kfp_particle.Chi2() = m_dst_track->get_chisq();
0271   kfp_particle.SetId(m_dst_track->get_id());
0272 
0273   return kfp_particle;
0274 }
0275 
0276 std::vector<KFParticle> KFParticle_Tools::makeAllDaughterParticles(PHCompositeNode *topNode)
0277 {
0278   std::vector<KFParticle> daughterParticles;
0279   m_dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name.c_str());
0280   unsigned int trackID = 0;
0281 
0282   for (auto &iter : *m_dst_trackmap)
0283   {
0284     m_dst_track = iter.second;
0285 
0286     if (m_bunch_crossing_zero_only && (m_dst_track->get_crossing() != 0))
0287     {
0288       continue;
0289     }
0290 
0291     //Now lets check if we have the right number of tracker states
0292     int MVTX_states = 0;
0293     int INTT_states = 0;
0294     int TPC_states = 0;
0295     int TPOT_states = 0;
0296 
0297     for (auto state_iter = m_dst_track->begin_states();
0298          state_iter != m_dst_track->end_states();
0299          ++state_iter)
0300     {
0301       SvtxTrackState* tstate = state_iter->second;
0302       if (tstate->get_pathlength() != 0) //The first track state is an extrapolation so has no cluster
0303       {
0304         auto stateckey = tstate->get_cluskey(); 
0305         uint8_t id = TrkrDefs::getTrkrId(stateckey);
0306       
0307         switch (id)
0308         {
0309           case TrkrDefs::mvtxId:
0310             ++MVTX_states;
0311             break;
0312           case TrkrDefs::inttId:
0313             ++INTT_states;
0314             break;
0315           case TrkrDefs::tpcId:
0316             ++TPC_states;
0317             break;
0318           case TrkrDefs::micromegasId:
0319             ++TPOT_states;
0320             break;
0321           default:
0322             //std::cout << "Cluster key doesnt match a tracking system, could be related with projected track state to calorimeter system" << std::endl;
0323             break; 
0324         }
0325       }
0326     }
0327 
0328     if (MVTX_states < m_nMVTXStates)
0329     {
0330       continue;
0331     }
0332     if (INTT_states < m_nINTTStates)
0333     {
0334       continue;
0335     }
0336     if (TPC_states < m_nTPCStates)
0337     {
0338       continue;
0339     }
0340     if (TPOT_states < m_nTPOTStates)
0341     {
0342       continue;
0343     }
0344  
0345     daughterParticles.push_back(makeParticle(topNode));  /// Turn all dst tracks in KFP tracks
0346     daughterParticles[trackID].SetId(iter.first);
0347     ++trackID;
0348   }
0349 
0350   return daughterParticles;
0351 }
0352 
0353 void KFParticle_Tools::getTracksFromBC(PHCompositeNode *topNode, const int &bunch_crossing, const std::string &vertexMapName, int &nTracks, int &nPVs)
0354 {
0355   if (m_use_mbd_vertex) //If you're using the MBD vertex then there is no way to know which tracks are associated to it
0356   {
0357     return;
0358   }
0359 
0360   std::string vtxMN;
0361   if (vertexMapName.empty())
0362   {
0363     vtxMN = m_vtx_map_node_name;
0364   }
0365   else
0366   {
0367     vtxMN = vertexMapName;
0368   }
0369 
0370   nTracks = 0;
0371   nPVs = 0;
0372   m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, vtxMN);
0373   for (SvtxVertexMap::ConstIter iter = m_dst_vertexmap->begin(); iter != m_dst_vertexmap->end(); ++iter)
0374   {
0375     m_dst_vertex = iter->second;
0376     if ((int) m_dst_vertex->get_beam_crossing() == bunch_crossing)
0377     {
0378       nTracks += m_dst_vertex->size_tracks();
0379       ++nPVs;
0380     }
0381   }
0382 }
0383 
0384 int KFParticle_Tools::getTracksFromVertex(PHCompositeNode *topNode, const KFParticle &vertex, const std::string &vertexMapName)
0385 {
0386   if (m_use_mbd_vertex) //If you're using the MBD vertex then there is no way to know which tracks are associated to it
0387   {
0388     return 0;
0389   }
0390 
0391   std::string vtxMN;
0392   if (vertexMapName.empty())
0393   {
0394     vtxMN = m_vtx_map_node_name;
0395   }
0396   else
0397   {
0398     vtxMN = vertexMapName;
0399   }
0400 
0401   SvtxVertex* associated_vertex = nullptr;
0402   if (m_dont_use_global_vertex)
0403   {
0404     m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, vtxMN);
0405     associated_vertex = m_dst_vertexmap->get(vertex.Id());
0406   }
0407   else
0408   {
0409     m_dst_globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0410     if (!m_dst_globalvertexmap)
0411     {
0412       std::cout << "Can't continue in KFParticle_Tools::makeAllPrimaryVertices" << std::endl;
0413       return 0;
0414     }
0415     auto associated_gvertex = m_dst_globalvertexmap->get(vertex.Id());
0416 
0417     auto svtxiter = associated_gvertex->find_vertexes(GlobalVertex::SVTX);
0418     auto svtxvertexvector = svtxiter->second;
0419 
0420     for (auto &gvertex : svtxvertexvector)
0421     {
0422       associated_vertex = m_dst_vertexmap->find(gvertex->get_id())->second;
0423     }
0424   }
0425   if (associated_vertex)
0426   {
0427     return associated_vertex->size_tracks();
0428   }
0429   else
0430   {
0431     return 0;
0432   }
0433 }
0434 
0435 /*const*/ bool KFParticle_Tools::isGoodTrack(const KFParticle &particle, const std::vector<KFParticle> &primaryVertices)
0436 {
0437   bool goodTrack = false;
0438 
0439   float min_ip = 0;
0440   float min_ipchi2 = 0;
0441   float min_ip_xy = 0;
0442   float min_ipchi2_xy = 0;
0443 
0444   float pt = particle.GetPt();
0445   float pterr = particle.GetErrPt();
0446   float ptchi2 = pow(pterr / pt, 2);
0447   float trackchi2ndof = particle.GetChi2() / particle.GetNDF();
0448   calcMinIP(particle, primaryVertices, min_ip, min_ipchi2);
0449   calcMinIP(particle, primaryVertices, min_ip_xy, min_ipchi2_xy, false);
0450 
0451   if (isInRange(m_track_min_pt, pt, m_track_max_pt)
0452    && ptchi2 <= m_track_ptchi2
0453    && min_ip >= m_track_ip
0454    && min_ipchi2 >= m_track_ipchi2 
0455    && min_ip_xy >= m_track_ip_xy
0456    && min_ipchi2_xy >= m_track_ipchi2_xy 
0457    && trackchi2ndof <= m_track_chi2ndof)
0458   {
0459     goodTrack = true;
0460   }
0461   return goodTrack;
0462 }
0463 
0464 int KFParticle_Tools::calcMinIP(const KFParticle &track, const std::vector<KFParticle> &PVs,
0465                                 float &minimumIP, float &minimumIPchi2, bool do3D)
0466 {
0467   std::vector<float> ip, ipchi2;
0468 
0469   for (const auto &PV : PVs)
0470   {
0471     float thisIPchi2 = 0;
0472 
0473     if (do3D) 
0474     {
0475       ip.push_back(track.GetDistanceFromVertex(PV));
0476       track.GetDeviationFromVertex(PV);
0477     }
0478     else
0479     {
0480       ip.push_back(abs(track.GetDistanceFromVertexXY(PV)));
0481       track.GetDeviationFromVertexXY(PV);
0482     }
0483 
0484     if (thisIPchi2 < 0)
0485     {
0486       thisIPchi2 = 0;
0487     }
0488     ipchi2.push_back(thisIPchi2);  //Τhere are times where the IPchi2 calc fails
0489   }
0490 
0491   auto minmax_ip = minmax_element(ip.begin(), ip.end());  // Order the IP from small to large
0492   minimumIP = *minmax_ip.first;
0493   auto minmax_ipchi2 = minmax_element(ipchi2.begin(), ipchi2.end());  // Order the IP chi2 from small to large
0494   minimumIPchi2 = *minmax_ipchi2.first;
0495 
0496   return 0;
0497 }
0498 
0499 std::vector<int> KFParticle_Tools::findAllGoodTracks(const std::vector<KFParticle> &daughterParticles, const std::vector<KFParticle> &primaryVertices)
0500 {
0501   std::vector<int> goodTrackIndex;
0502 
0503   for (unsigned int i_parts = 0; i_parts < daughterParticles.size(); ++i_parts)
0504   {
0505     if (isGoodTrack(daughterParticles[i_parts], primaryVertices))
0506     {
0507       goodTrackIndex.push_back(i_parts);
0508     }
0509   }
0510 
0511   removeDuplicates(goodTrackIndex);
0512 
0513   return goodTrackIndex;
0514 }
0515 
0516 std::vector<std::vector<int>> KFParticle_Tools::findTwoProngs(std::vector<KFParticle> daughterParticles, std::vector<int> goodTrackIndex, int nTracks)
0517 {
0518   std::vector<std::vector<int>> goodTracksThatMeet;
0519 
0520   for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
0521   {
0522     for (std::vector<int>::iterator j_it = goodTrackIndex.begin(); j_it != goodTrackIndex.end(); ++j_it)
0523     {
0524       if (i_it < j_it)
0525       {
0526         float dca = daughterParticles[*i_it].GetDistanceFromParticle(daughterParticles[*j_it]);
0527         float dca_xy = abs(daughterParticles[*i_it].GetDistanceFromParticleXY(daughterParticles[*j_it]));
0528 
0529         if (dca <= m_comb_DCA && dca_xy <= m_comb_DCA_xy)
0530         {
0531           KFVertex twoParticleVertex;
0532           twoParticleVertex += daughterParticles[*i_it];
0533           twoParticleVertex += daughterParticles[*j_it];
0534           float vertexchi2ndof = twoParticleVertex.GetChi2() / twoParticleVertex.GetNDF();
0535           float sv_radial_position = sqrt(pow(twoParticleVertex.GetX(), 2) + pow(twoParticleVertex.GetY(), 2));
0536           std::vector<int> combination = {*i_it, *j_it};
0537 
0538           if (nTracks == 2 && vertexchi2ndof > m_vertex_chi2ndof)
0539           {
0540             continue;
0541           }
0542           else
0543           {
0544             if (nTracks == 2 && sv_radial_position < m_min_radial_SV)
0545             {
0546               continue;
0547             }
0548             else
0549             {
0550               goodTracksThatMeet.push_back(combination);
0551             }
0552           }
0553         }
0554       }
0555     }
0556   }
0557 
0558   return goodTracksThatMeet;
0559 }
0560 
0561 std::vector<std::vector<int>> KFParticle_Tools::findNProngs(std::vector<KFParticle> daughterParticles,
0562                                                             const std::vector<int> &goodTrackIndex,
0563                                                             std::vector<std::vector<int>> goodTracksThatMeet,
0564                                                             int nRequiredTracks, unsigned int nProngs)
0565 {
0566   unsigned int nGoodProngs = goodTracksThatMeet.size();
0567 
0568   for (auto &i_it : goodTrackIndex)
0569   {
0570     for (unsigned int i_prongs = 0; i_prongs < nGoodProngs; ++i_prongs)
0571     {
0572       bool trackNotUsedAlready = true;
0573       for (unsigned int i_trackCheck = 0; i_trackCheck < nProngs - 1; ++i_trackCheck)
0574       {
0575         if (i_it == goodTracksThatMeet[i_prongs][i_trackCheck])
0576         {
0577           trackNotUsedAlready = false;
0578         }
0579       }
0580       if (trackNotUsedAlready)
0581       {
0582         bool dcaMet = true;
0583         for (unsigned int i = 0; i < nProngs - 1; ++i)
0584         {
0585           float dca = daughterParticles[i_it].GetDistanceFromParticle(daughterParticles[goodTracksThatMeet[i_prongs][i]]);
0586           float dca_xy = abs(daughterParticles[i_it].GetDistanceFromParticleXY(daughterParticles[goodTracksThatMeet[i_prongs][i]]));
0587 
0588           if (dca > m_comb_DCA || dca_xy > m_comb_DCA_xy)
0589           {
0590             dcaMet = false;
0591           }
0592         }
0593 
0594         if (dcaMet)
0595         {
0596           KFVertex particleVertex;
0597           particleVertex += daughterParticles[i_it];
0598           std::vector<int> combination;
0599           combination.push_back(i_it);
0600           for (unsigned int i = 0; i < nProngs - 1; ++i)
0601           {
0602             particleVertex += daughterParticles[goodTracksThatMeet[i_prongs][i]];
0603             combination.push_back(goodTracksThatMeet[i_prongs][i]);
0604           }
0605           float vertexchi2ndof = particleVertex.GetChi2() / particleVertex.GetNDF();
0606           float sv_radial_position = sqrt(pow(particleVertex.GetX(), 2) + pow(particleVertex.GetY(), 2));
0607 
0608           if ((unsigned int) nRequiredTracks == nProngs && vertexchi2ndof > m_vertex_chi2ndof)
0609           {
0610             continue;
0611           }
0612           else
0613           {
0614             if ((unsigned int) nRequiredTracks == nProngs && sv_radial_position < m_min_radial_SV)
0615             {
0616               continue;
0617             }
0618             else
0619             {
0620               goodTracksThatMeet.push_back(combination);
0621             }
0622           }
0623         }
0624       }
0625     }
0626   }
0627 
0628   goodTracksThatMeet.erase(goodTracksThatMeet.begin(), goodTracksThatMeet.begin() + nGoodProngs);
0629   for (auto &i : goodTracksThatMeet)
0630   {
0631     sort(i.begin(), i.end());
0632   }
0633   removeDuplicates(goodTracksThatMeet);
0634 
0635   return goodTracksThatMeet;
0636 }
0637 
0638 std::vector<std::vector<int>> KFParticle_Tools::appendTracksToIntermediates(KFParticle intermediateResonances[], const std::vector<KFParticle> &daughterParticles, const std::vector<int> &goodTrackIndex, int num_remaining_tracks)
0639 {
0640   std::vector<std::vector<int>> goodTracksThatMeet, goodTracksThatMeetIntermediates;  //, vectorOfGoodTracks;
0641   if (num_remaining_tracks == 1)
0642   {
0643     for (auto &i_it : goodTrackIndex)
0644     {
0645       std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances + m_num_intermediate_states);
0646       std::vector<std::vector<int>> dummyTrackList;
0647       std::vector<int> dummyTrackID;  // I already have the track ids stored in goodTracksThatMeet[i]
0648       v_intermediateResonances.insert(end(v_intermediateResonances), daughterParticles[i_it]);
0649       for (unsigned int k = 0; k < v_intermediateResonances.size(); ++k)
0650       {
0651         dummyTrackID.push_back(k);
0652       }
0653       dummyTrackList = findTwoProngs(v_intermediateResonances, dummyTrackID, (int) v_intermediateResonances.size());
0654       if (v_intermediateResonances.size() > 2)
0655       {
0656         for (unsigned int p = 3; p <= v_intermediateResonances.size(); ++p)
0657         {
0658           dummyTrackList = findNProngs(v_intermediateResonances,
0659                                        dummyTrackID, dummyTrackList,
0660                                        (int) v_intermediateResonances.size(), (int) p);
0661         }
0662       }
0663 
0664       if (dummyTrackList.size() != 0)
0665       {
0666         std::vector<int> goodTrack{i_it};
0667         goodTracksThatMeetIntermediates.push_back(goodTrack);
0668       }
0669     }
0670   }
0671   else
0672   {
0673     goodTracksThatMeet = findTwoProngs(daughterParticles, goodTrackIndex, num_remaining_tracks);
0674 
0675     for (int p = 3; p <= num_remaining_tracks; ++p)
0676     {
0677       goodTracksThatMeet = findNProngs(daughterParticles, goodTrackIndex, goodTracksThatMeet, num_remaining_tracks, p);
0678     }
0679 
0680     for (auto &i : goodTracksThatMeet)
0681     {
0682       std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances + m_num_intermediate_states);
0683       std::vector<std::vector<int>> dummyTrackList;
0684       std::vector<int> dummyTrackID;  // I already have the track ids stored in goodTracksThatMeet[i]
0685       for (int j : i)
0686       {
0687         v_intermediateResonances.push_back(daughterParticles[i[j]]);
0688       }
0689       for (unsigned int k = 0; k < v_intermediateResonances.size(); ++k)
0690       {
0691         dummyTrackID.push_back(k);
0692       }
0693       dummyTrackList = findTwoProngs(v_intermediateResonances, dummyTrackID, (int) v_intermediateResonances.size());
0694       for (unsigned int p = 3; p <= v_intermediateResonances.size(); ++p)
0695       {
0696         dummyTrackList = findNProngs(v_intermediateResonances, dummyTrackID, dummyTrackList, (int) v_intermediateResonances.size(), (int) p);
0697       }
0698 
0699       if (dummyTrackList.size() != 0)
0700       {
0701         goodTracksThatMeetIntermediates.push_back(i);
0702       }
0703     }
0704   }
0705 
0706   return goodTracksThatMeetIntermediates;
0707 }
0708 
0709 float KFParticle_Tools::eventDIRA(const KFParticle &particle, const KFParticle &vertex, bool do3D)
0710 {
0711   const int nDimensions = do3D ? 3 : 2;
0712   TMatrixD flightVector(nDimensions, 1);
0713   TMatrixD momVector(nDimensions, 1);
0714   flightVector(0, 0) = particle.GetX() - vertex.GetX();
0715   flightVector(1, 0) = particle.GetY() - vertex.GetY();
0716 
0717   momVector(0, 0) = particle.GetPx();
0718   momVector(1, 0) = particle.GetPy();
0719   
0720   if (do3D)
0721   {
0722     flightVector(2, 0) = particle.GetZ() - vertex.GetZ();
0723     momVector(2, 0) = particle.GetPz();
0724   }
0725 
0726   TMatrixD momDotFD(1, 1);  // Calculate momentum dot flight distance
0727   momDotFD = TMatrixD(momVector, TMatrixD::kTransposeMult, flightVector);
0728   float f_momDotFD = momDotFD(0, 0);
0729 
0730   TMatrixD sizeOfMom(1, 1);  // Calculates the size of the momentum vector
0731   sizeOfMom = TMatrixD(momVector, TMatrixD::kTransposeMult, momVector);
0732   float f_sizeOfMom = sqrt(sizeOfMom(0, 0));
0733 
0734   TMatrixD sizeOfFD(1, 1);  // Calculates the size of the flight distance vector
0735   sizeOfFD = TMatrixD(flightVector, TMatrixD::kTransposeMult, flightVector);
0736   float f_sizeOfFD = sqrt(sizeOfFD(0, 0));
0737 
0738   return f_momDotFD / (f_sizeOfMom * f_sizeOfFD);
0739 }
0740 
0741 float KFParticle_Tools::flightDistanceChi2(const KFParticle &particle, const KFParticle &vertex)
0742 {
0743   TMatrixD flightVector(3, 1);
0744   TMatrixD flightDistanceCovariance(3, 3);
0745 
0746   const KFParticle &kfp_vertex(vertex);
0747 
0748   flightVector(0, 0) = particle.GetX() - kfp_vertex.GetX();
0749   flightVector(1, 0) = particle.GetY() - kfp_vertex.GetY();
0750   flightVector(2, 0) = particle.GetZ() - kfp_vertex.GetZ();
0751 
0752   for (int i = 0; i < 3; i++)
0753   {
0754     for (int j = 0; j < 3; j++)
0755     {
0756       flightDistanceCovariance(i, j) = particle.GetCovariance(i, j) + kfp_vertex.GetCovariance(i, j);
0757     }
0758   }
0759 
0760   TMatrixD anInverseMatrix(3, 3);
0761   anInverseMatrix = flightDistanceCovariance.Invert();
0762   TMatrixD m_chi2Value(1, 1);
0763   m_chi2Value = TMatrixD(flightVector, TMatrixD::kTransposeMult, anInverseMatrix * flightVector);
0764   return m_chi2Value(0, 0);
0765 }
0766 
0767 std::tuple<KFParticle, bool> KFParticle_Tools::buildMother(KFParticle vDaughters[], int daughterOrder[],
0768                                                            bool isIntermediate, int intermediateNumber, int nTracks,
0769                                                            bool constrainMass, float required_vertexID, PHCompositeNode* topNode)
0770 {
0771   KFParticle mother;
0772   KFParticle *inputTracks = new KFParticle[nTracks];
0773 
0774   mother.SetConstructMethod(2);
0775 
0776   bool daughterMassCheck = true;
0777   int particlesWithPID[] = {11, 211, 321, 2212};
0778   float unique_vertexID = 0;
0779 
0780   // Figure out if the decay has reco. tracks mixed with resonances
0781   int num_tracks_used_by_intermediates = 0;
0782   for (int i = 0; i < m_num_intermediate_states; ++i)
0783   {
0784     num_tracks_used_by_intermediates += m_num_tracks_from_intermediate[i];
0785   }
0786   int num_remaining_tracks = m_num_tracks - num_tracks_used_by_intermediates;
0787 
0788   for (int i = 0; i < nTracks; ++i)
0789   {
0790     float daughterMass = 0;
0791 
0792     if ((Int_t) vDaughters[i].GetQ() != 0)
0793     {
0794       // For charged particle, like p+/-, pi+/-, Sigma+/-, etc...
0795       // different charged particle has different PDGID
0796       // just to protect if they have different mass for different charge
0797       // but in EvtGen, there is no C-violation...so this is just a protection
0798       daughterMass = constrainMass ? getParticleMass((Int_t) vDaughters[i].GetQ() * daughterOrder[i]) : vDaughters[i].GetMass();
0799     }
0800     else if ((Int_t) vDaughters[i].GetQ() == 0)
0801     {
0802       // For neutral particle, like pi0, eta, J/psi, etc... who do not have an anti-particle with anti-PDGID
0803       // and other neutral particle, like Lambda0/anti-Lambda0 ... who have an anti-particle with anti-PDGID
0804       // avoid charge*PDGID=0 case and getting wrong mass
0805       daughterMass = constrainMass ? getParticleMass(daughterOrder[i]) : vDaughters[i].GetMass();
0806     }
0807 
0808     if ((num_remaining_tracks > 0 && i >= m_num_intermediate_states) || isIntermediate)
0809     {
0810       if ((Int_t) vDaughters[i].GetQ() != 0)
0811       {
0812         daughterMass = getParticleMass((Int_t) vDaughters[i].GetQ() * daughterOrder[i]);
0813       }
0814       else if ((Int_t) vDaughters[i].GetQ() == 0)
0815       {
0816         daughterMass = getParticleMass(daughterOrder[i]);
0817       }
0818 
0819     }
0820     inputTracks[i].Create(vDaughters[i].Parameters(),
0821                           vDaughters[i].CovarianceMatrix(),
0822                           (Int_t) vDaughters[i].GetQ(),
0823                           daughterMass);
0824 
0825     //Run PID check
0826     if (m_use_PID)
0827     {
0828       int track_PDG_ID = (Int_t) vDaughters[i].GetQ()*daughterOrder[i];
0829       if (std::find(std::begin(particlesWithPID), std::end(particlesWithPID), std::abs(track_PDG_ID)) != std::end(particlesWithPID))
0830       {
0831         float calculated_dEdx_value = get_dEdx(topNode, vDaughters[i]);
0832         double expected_dEdx_value = get_dEdx_fitValue((Int_t) vDaughters[i].GetQ() * vDaughters[i].GetP(), track_PDG_ID);
0833         bool accept_dEdx = isInRange((1-m_dEdx_band_width)*expected_dEdx_value, calculated_dEdx_value, (1+m_dEdx_band_width)*expected_dEdx_value);
0834         if(!accept_dEdx)
0835         {
0836          delete [] inputTracks;
0837          return std::make_tuple(mother, false);
0838         }
0839       }
0840     }
0841 
0842     mother.AddDaughter(inputTracks[i]);
0843     mother.AddDaughterId(vDaughters[i].Id());
0844     unique_vertexID += (Int_t) vDaughters[i].GetQ() * getParticleMass(daughterOrder[i]);
0845   }
0846 
0847   if (isIntermediate)
0848   {
0849     mother.SetPDG(getParticleID(m_intermediate_name[intermediateNumber].c_str()));
0850   }
0851   if (!isIntermediate && !m_mother_name_Tools.empty())
0852   {
0853     mother.SetPDG(getParticleID(m_mother_name_Tools));
0854   }
0855 
0856   bool chargeCheck;
0857   if (m_get_charge_conjugate)
0858   {
0859     chargeCheck = std::abs(unique_vertexID) == std::abs(required_vertexID) ? true : false;
0860   }
0861   else
0862   {
0863     chargeCheck = unique_vertexID == required_vertexID ? true : false;
0864   }
0865 
0866   for (int j = 0; j < nTracks; ++j)
0867   {
0868     if (m_extrapolateTracksToSV)
0869     {
0870       inputTracks[j].SetProductionVertex(mother);
0871     }
0872     if (!m_allowZeroMassTracks)
0873     {
0874       if (inputTracks[j].GetMass() == 0)
0875       {
0876         daughterMassCheck = false;
0877       }
0878     }
0879   }
0880 
0881 
0882   float calculated_mass, calculated_mass_err;
0883   mother.GetMass(calculated_mass, calculated_mass_err);
0884   float calculated_pt = mother.GetPt();
0885 
0886   float min_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].first : m_min_mass;
0887   float max_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].second : m_max_mass;
0888   float min_pt = isIntermediate ? m_intermediate_min_pt[intermediateNumber] : m_mother_pt;
0889 
0890   float max_vertex_volume = isIntermediate ? m_intermediate_vertex_volume[intermediateNumber] : m_mother_vertex_volume;
0891 
0892   bool goodCandidate = false;
0893 
0894   if (calculated_mass >= min_mass && calculated_mass <= max_mass &&
0895       calculated_pt >= min_pt && daughterMassCheck && chargeCheck && calculateEllipsoidVolume(mother) <= max_vertex_volume)
0896   {
0897     goodCandidate = true;
0898   }
0899 
0900   if (goodCandidate && m_require_bunch_crossing_match)
0901   {
0902     std::vector<int> crossings;
0903     for (int i = 0; i < nTracks; ++i)
0904     {
0905       SvtxTrack *thisTrack = toolSet.getTrack(vDaughters[i].Id(), m_dst_trackmap);
0906       if (thisTrack)//This protects against intermediates which have no track but I need a way to assign the bunch crossing to an interemdiate as this was already checked when it was actually built
0907       {
0908         crossings.push_back(thisTrack->get_crossing());
0909       }
0910     }
0911 
0912     removeDuplicates(crossings);
0913 
0914     if (crossings.size() != 1)
0915     {
0916       goodCandidate = false;
0917     }
0918   }
0919 
0920   // Check the requirements of an intermediate states against this mother and re-do goodCandidate
0921   if (goodCandidate && m_has_intermediates && !isIntermediate)  // The decay has intermediate states and we are now looking at the mother
0922   {
0923     for (int k = 0; k < m_num_intermediate_states; ++k)
0924     {
0925       float intermediate_DIRA = eventDIRA(vDaughters[k], mother);
0926       float intermediate_FDchi2 = flightDistanceChi2(vDaughters[k], mother);
0927       if (intermediate_DIRA < m_intermediate_min_dira[k] ||
0928           intermediate_FDchi2 < m_intermediate_min_fdchi2[k])
0929       {
0930         goodCandidate = false;
0931       }
0932     }
0933   }
0934   delete [] inputTracks;
0935   return std::make_tuple(mother, goodCandidate);
0936 }
0937 
0938 void KFParticle_Tools::constrainToVertex(KFParticle &particle, bool &goodCandidate, KFParticle &vertex)
0939 {
0940   KFParticle particleCopy = particle;
0941   particleCopy.SetProductionVertex(vertex);
0942   particleCopy.TransportToDecayVertex();
0943 
0944   float calculated_decayTime;
0945   float calculated_decayTimeErr;
0946   float calculated_decayLength;
0947   float calculated_decayLengthErr;
0948   float calculated_decayLength_xy;
0949   float calculated_decayLengthErr_xy;
0950 
0951   particleCopy.GetLifeTime(calculated_decayTime, calculated_decayTimeErr);
0952   particleCopy.GetDecayLength(calculated_decayLength, calculated_decayLengthErr);
0953   particleCopy.GetDecayLengthXY(calculated_decayLength_xy, calculated_decayLengthErr_xy);
0954 
0955   float calculated_decayTime_xy = particleCopy.GetPseudoProperDecayTime(vertex, particleCopy.GetMass());
0956 
0957   float calculated_fdchi2 = flightDistanceChi2(particle, vertex);
0958 
0959   float calculated_ip_xy = abs(particle.GetDistanceFromVertexXY(vertex));
0960   float calculated_ipchi2_xy = particle.GetDeviationFromVertexXY(vertex);
0961   float calculated_dira_xy = eventDIRA(particle, vertex, false);
0962   float calculated_ip = particle.GetDistanceFromVertex(vertex);
0963   float calculated_ipchi2 = particle.GetDeviationFromVertex(vertex);
0964   float calculated_dira = eventDIRA(particle, vertex);
0965 
0966   float calculated_decay_time_significance = calculated_decayTime/calculated_decayTimeErr;
0967   float calculated_decay_length_significance = calculated_decayLength/calculated_decayLengthErr;
0968   float calculated_decay_length_xy_significance = calculated_decayLength_xy/calculated_decayLengthErr_xy;
0969 
0970   goodCandidate = false;
0971 
0972   const float speed = 2.99792458e-2;
0973   calculated_decayTime /= speed;
0974 
0975   if (calculated_fdchi2 >= m_fdchi2 
0976    && calculated_ip <= m_mother_ip 
0977    && calculated_ipchi2 <= m_mother_ipchi2 
0978    && calculated_ip_xy <= m_mother_ip_xy
0979    && calculated_ipchi2_xy <= m_mother_ipchi2_xy
0980    && calculated_decay_time_significance >= m_mother_min_decay_time_significance
0981    && calculated_decay_length_significance >= m_mother_min_decay_length_significance
0982    && calculated_decay_length_xy_significance >= m_mother_min_decay_length_xy_significance
0983    && isInRange(m_dira_min, calculated_dira, m_dira_max) 
0984    && isInRange(m_dira_xy_min, calculated_dira_xy, m_dira_xy_max) 
0985    && isInRange(m_min_decayTime, calculated_decayTime, m_max_decayTime) 
0986    && isInRange(m_min_decayTime_xy, calculated_decayTime_xy, m_max_decayTime_xy) 
0987    && isInRange(m_min_decayLength, calculated_decayLength, m_max_decayLength)
0988    && isInRange(m_min_decayLength_xy, calculated_decayLength_xy, m_max_decayLength_xy))
0989   {
0990     goodCandidate = true;
0991   }
0992 }
0993 
0994 std::tuple<KFParticle, bool> KFParticle_Tools::getCombination(KFParticle vDaughters[], int daughterOrder[], KFParticle vertex, bool constrain_to_vertex, bool isIntermediate, int intermediateNumber, int nTracks, bool constrainMass, float required_vertexID, PHCompositeNode* topNode)
0995 {
0996   KFParticle candidate;
0997   bool isGoodCandidate;
0998 
0999   std::tie(candidate, isGoodCandidate) = buildMother(vDaughters, daughterOrder, isIntermediate, intermediateNumber, nTracks, constrainMass, required_vertexID, topNode);
1000 
1001   if (constrain_to_vertex && isGoodCandidate && !isIntermediate)
1002   {
1003     if (m_require_track_and_vertex_match)
1004     {
1005       isGoodCandidate = checkTrackAndVertexMatch(vDaughters, nTracks, vertex);
1006     }
1007 
1008     if (isGoodCandidate)
1009     {
1010       constrainToVertex(candidate, isGoodCandidate, vertex);
1011     }
1012   }
1013   return std::make_tuple(candidate, isGoodCandidate);
1014 }
1015 
1016 std::vector<std::vector<int>> KFParticle_Tools::findUniqueDaughterCombinations(int start, int end)
1017 {
1018   std::vector<int> vect_permutations;
1019   std::vector<std::vector<int>> uniqueCombinations;
1020   std::map<int, int> daughterMap;
1021   for (int i = start; i < end; i++)
1022   {
1023     daughterMap.insert(std::pair<int, int>(i, abs(getParticleID(m_daughter_name[i].c_str()))));
1024     vect_permutations.push_back(i);
1025   }
1026   int *permutations = &vect_permutations[0];
1027 
1028   do
1029   {
1030     std::vector<int> combination;
1031     combination.reserve((end - start));
1032     for (int i = 0; i < (end - start); i++)
1033     {
1034       combination.push_back(daughterMap.find(permutations[i])->second);
1035     }
1036     uniqueCombinations.push_back(combination);
1037   } while (std::next_permutation(permutations, permutations + vect_permutations.size()));
1038 
1039   removeDuplicates(uniqueCombinations);
1040 
1041   return uniqueCombinations;
1042 }
1043 
1044 double KFParticle_Tools::calculateEllipsoidRadius(int posOrNeg, double sigma_ii, double sigma_jj, double sigma_ij)
1045 {  // Note - Only works for a 2D ellipsoid OR rotated nD ellipsoid to avoid projections
1046   if (std::abs(posOrNeg) != 1)
1047   {
1048     std::cout << "You have set posOrNeg to " << posOrNeg << ". This value must be  +/- 1! Skipping" << std::endl;
1049     return 0;
1050   }
1051 
1052   double r_ij = sqrt((sigma_ii + sigma_jj) / 2 + posOrNeg * (sqrt(pow((sigma_ii - sigma_jj) / 2, 2) + pow(sigma_ij, 2))));
1053 
1054   return r_ij;
1055 }
1056 
1057 float KFParticle_Tools::calculateEllipsoidVolume(const KFParticle &particle)
1058 {
1059   TMatrixD cov_matrix(3, 3);
1060 
1061   for (int i = 0; i < 3; ++i)
1062   {
1063     for (int j = 0; j < 3; ++j)
1064     {
1065       cov_matrix(i, j) = particle.GetCovariance(i, j);
1066     }
1067   }
1068 
1069   float volume;
1070   if (cov_matrix(0, 0) * cov_matrix(1, 1) * cov_matrix(2, 2) == 0)
1071   {
1072     volume = 0;
1073   }
1074   else
1075   {
1076     volume = (4. / 3.) * M_PI * sqrt((std::abs(cov_matrix.Determinant())));  // The covariance matrix is error-squared
1077   }
1078 
1079   return volume;
1080 }
1081 
1082 float KFParticle_Tools::calculateJT(const KFParticle &mother, const KFParticle &daughter)
1083 {
1084   Eigen::Vector3f motherP = Eigen::Vector3f(mother.GetPx(), mother.GetPy(), mother.GetPz());
1085   Eigen::Vector3f daughterP = Eigen::Vector3f(daughter.GetPx(), daughter.GetPy(), daughter.GetPz());
1086 
1087   Eigen::Vector3f motherP_X_daughterP = motherP.cross(daughterP);
1088 
1089   float jT = (motherP_X_daughterP.norm()) / motherP.norm();
1090 
1091   return jT;
1092 }
1093 
1094 bool KFParticle_Tools::isInRange(float min, float value, float max)
1095 {
1096   return min <= value && value <= max;
1097 }
1098 
1099 void KFParticle_Tools::removeDuplicates(std::vector<double> &v)
1100 {
1101   auto end = v.end();
1102   for (auto it = v.begin(); it != end; ++it)
1103   {
1104     end = remove(it + 1, end, *it);
1105   }
1106   v.erase(end, v.end());
1107 }
1108 
1109 void KFParticle_Tools::removeDuplicates(std::vector<int> &v)
1110 {
1111   auto end = v.end();
1112   for (auto it = v.begin(); it != end; ++it)
1113   {
1114     end = remove(it + 1, end, *it);
1115   }
1116   v.erase(end, v.end());
1117 }
1118 
1119 void KFParticle_Tools::removeDuplicates(std::vector<std::vector<int>> &v)
1120 {
1121   auto end = v.end();
1122   for (auto it = v.begin(); it != end; ++it)
1123   {
1124     end = remove(it + 1, end, *it);
1125   }
1126   v.erase(end, v.end());
1127 }
1128 
1129 void KFParticle_Tools::removeDuplicates(std::vector<std::vector<std::string>> &v)
1130 {
1131   auto end = v.end();
1132   for (auto it = v.begin(); it != end; ++it)
1133   {
1134     end = remove(it + 1, end, *it);
1135   }
1136   v.erase(end, v.end());
1137 }
1138 
1139 bool KFParticle_Tools::findParticle(const std::string &particle)
1140 {
1141   bool particleFound = true;
1142   if (!TDatabasePDG::Instance()->GetParticle(particle.c_str()))
1143   {
1144     std::cout << "The particle, " << particle << " is not recognized" << std::endl;
1145     std::cout << "Check TDatabasePDG for a list of available particles" << std::endl;
1146     particleFound = false;
1147   }
1148 
1149   return particleFound;
1150 }
1151 
1152 int KFParticle_Tools::getParticleID(const std::string &particle)
1153 {
1154   return TDatabasePDG::Instance()->GetParticle(particle.c_str())->PdgCode();
1155 }
1156 
1157 float KFParticle_Tools::getParticleMass(const std::string &particle)
1158 {
1159   return TDatabasePDG::Instance()->GetParticle(particle.c_str())->Mass();
1160 }
1161 
1162 float KFParticle_Tools::getParticleMass(const int PDGID)
1163 {
1164   return TDatabasePDG::Instance()->GetParticle(PDGID)->Mass();
1165 }
1166 
1167 void KFParticle_Tools::identify(const KFParticle &particle)
1168 {
1169   std::cout << "Track ID: " << particle.Id() << std::endl;
1170   std::cout << "PDG ID: " << particle.GetPDG() << ", charge: " << (int) particle.GetQ() << ", mass: " << particle.GetMass() << " GeV" << std::endl;
1171   std::cout << "(px,py,pz) = (" << particle.GetPx() << " +/- " << std::sqrt(particle.GetCovariance(3, 3)) << ", ";
1172   std::cout << particle.GetPy() << " +/- " << std::sqrt(particle.GetCovariance(4, 4)) << ", ";
1173   std::cout << particle.GetPz() << " +/- " << std::sqrt(particle.GetCovariance(5, 5)) << ") GeV" << std::endl;
1174   std::cout << "(x,y,z) = (" << particle.GetX() << " +/- " << std::sqrt(particle.GetCovariance(0, 0)) << ", ";
1175   std::cout << particle.GetY() << " +/- " << std::sqrt(particle.GetCovariance(1, 1)) << ", ";
1176   std::cout << particle.GetZ() << " +/- " << std::sqrt(particle.GetCovariance(2, 2)) << ") cm\n"
1177             << std::endl;
1178 }
1179 
1180 float KFParticle_Tools::get_dEdx(PHCompositeNode *topNode, const KFParticle &daughter)
1181 {
1182   m_dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name.c_str());
1183   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1184   m_geom_container = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
1185   auto geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1186   if(!m_cluster_map || !m_geom_container || !geometry)
1187   {
1188     return -1.0;
1189   }
1190 
1191   SvtxTrack *daughter_track = toolSet.getTrack(daughter.Id(), m_dst_trackmap);
1192   TrackSeed *tpcseed = daughter_track->get_tpc_seed();
1193   float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
1194   // These are randomly chosen layer thicknesses for the TPC, to get the
1195   // correct region thicknesses in an easy to pass way to the helper fxn
1196   layerThicknesses[0] = m_geom_container->GetLayerCellGeom(7)->get_thickness();
1197   layerThicknesses[1] = m_geom_container->GetLayerCellGeom(8)->get_thickness();
1198   layerThicknesses[2] = m_geom_container->GetLayerCellGeom(27)->get_thickness();
1199   layerThicknesses[3] = m_geom_container->GetLayerCellGeom(50)->get_thickness();
1200 
1201   return TrackAnalysisUtils::calc_dedx(tpcseed, m_cluster_map, geometry, layerThicknesses);
1202     
1203 }
1204 
1205 void KFParticle_Tools::init_dEdx_fits()
1206 {
1207   std::string dedx_fitparams = CDBInterface::instance()->getUrl("TPC_DEDX_FITPARAM");
1208   TFile *filefit = TFile::Open(dedx_fitparams.c_str());
1209 
1210   if (!filefit->IsOpen())
1211   {
1212       std::cerr << "Error opening filefit!" << std::endl;
1213       return;
1214   }
1215 
1216   filefit->GetObject("f_piband", f_pion_plus);
1217   filefit->GetObject("f_Kband", f_kaon_plus);
1218   filefit->GetObject("f_pband", f_proton_plus);
1219   filefit->GetObject("f_piminus_band", f_pion_minus);
1220   filefit->GetObject("f_Kminus_band", f_kaon_minus);
1221   filefit->GetObject("f_pbar_band", f_proton_minus);
1222 
1223   pidMap.insert(std::pair<int, TF1*>( -11,  f_pion_plus));
1224   pidMap.insert(std::pair<int, TF1*>( 211,  f_pion_plus));
1225   pidMap.insert(std::pair<int, TF1*>( 321,  f_kaon_plus));
1226   pidMap.insert(std::pair<int, TF1*>( 2212, f_proton_plus));
1227   pidMap.insert(std::pair<int, TF1*>(  11,  f_pion_minus));
1228   pidMap.insert(std::pair<int, TF1*>(-211,  f_pion_minus));
1229   pidMap.insert(std::pair<int, TF1*>(-321,  f_kaon_minus));
1230   pidMap.insert(std::pair<int, TF1*>(-2212, f_proton_minus));
1231 }
1232 
1233 double KFParticle_Tools::get_dEdx_fitValue(float momentum, int PID)
1234 {
1235   return pidMap[PID]->Eval(momentum);
1236 }
1237 
1238 bool KFParticle_Tools::checkTrackAndVertexMatch(KFParticle vDaughters[], int nTracks, KFParticle vertex)
1239 {
1240   bool vertexAndTrackMatch = true;
1241 
1242   int vertexCrossing = 1e5;
1243 
1244   if (m_dont_use_global_vertex)
1245   {
1246     if (m_use_mbd_vertex)
1247     {
1248       m_dst_mbdvertex = m_dst_mbdvertexmap->get(vertex.Id());
1249       vertexCrossing = m_dst_mbdvertex->get_beam_crossing();
1250     }
1251     else
1252     {
1253       m_dst_vertex = m_dst_vertexmap->get(vertex.Id());
1254       vertexCrossing = m_dst_vertex->get_beam_crossing();
1255     }
1256   }
1257   else
1258   {
1259     m_dst_globalvertex = m_dst_globalvertexmap->get(vertex.Id());
1260 
1261     GlobalVertex::VTXTYPE whichVtx = m_use_mbd_vertex ? GlobalVertex::MBD : GlobalVertex::SVTX; 
1262 
1263     auto svtxiter = m_dst_globalvertex->find_vertexes(whichVtx);
1264     auto svtxvertexvector = svtxiter->second;
1265 
1266     for (auto &gvertex : svtxvertexvector)
1267     {
1268       if (m_use_mbd_vertex)
1269       {
1270         m_dst_mbdvertex = m_dst_mbdvertexmap->find(gvertex->get_id())->second;
1271         vertexCrossing = m_dst_mbdvertex->get_beam_crossing();
1272       }
1273       else
1274       {
1275         m_dst_vertex = m_dst_vertexmap->find(gvertex->get_id())->second;
1276         vertexCrossing = m_dst_vertex->get_beam_crossing();
1277       }
1278     }
1279   }
1280 
1281   for (int i = 0; i < nTracks; ++i)
1282   {
1283     SvtxTrack *thisTrack = toolSet.getTrack(vDaughters[i].Id(), m_dst_trackmap);
1284     if (thisTrack)//This protects against intermediates which have no track
1285     {
1286       int trackCrossing = thisTrack->get_crossing();
1287       vertexAndTrackMatch = trackCrossing != vertexCrossing ? false : vertexAndTrackMatch;
1288     }
1289   }
1290    
1291   return vertexAndTrackMatch;
1292 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!