Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:47 

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_eventReconstruction.h"
0029 #include "KFParticle_Tools.h"
0030 #include "KFParticle_truthAndDetTools.h"
0031 
0032 //sPHENIX stuff
0033 #include <trackbase_historic/SvtxTrack.h>
0034 
0035 // KFParticle stuff
0036 #include <KFParticle.h>
0037 
0038 #include <algorithm>
0039 #include <cassert>
0040 #include <iterator>  // for begin, distance, end
0041 #include <memory>   // for allocator_traits<>::value_type
0042 #include <string>   // for string
0043 #include <tuple>    // for tie, tuple
0044 
0045 #include <iostream>
0046 
0047 /// Create necessary objects
0048 KFParticle_Tools kfp_Tools_evtReco;
0049 
0050 /// KFParticle constructor
0051 KFParticle_eventReconstruction::KFParticle_eventReconstruction()
0052   : m_constrain_to_vertex(false)
0053   , m_constrain_int_mass(false)
0054   , m_use_fake_pv(false)
0055 {
0056 }
0057 
0058 void KFParticle_eventReconstruction::createDecay(PHCompositeNode* topNode, std::vector<KFParticle>& selectedMother, std::vector<KFParticle>& selectedVertex,
0059                                                  std::vector<std::vector<KFParticle>>& selectedDaughters,
0060                                                  std::vector<std::vector<KFParticle>>& selectedIntermediates,
0061                                                  int& nPVs)
0062 {
0063   std::vector<KFParticle> primaryVertices;
0064   if (m_use_fake_pv)
0065   {
0066     primaryVertices.push_back(createFakePV());
0067   }
0068   else
0069   {
0070     primaryVertices = makeAllPrimaryVertices(topNode, m_vtx_map_node_name);
0071   }
0072 
0073   std::vector<KFParticle> daughterParticles = makeAllDaughterParticles(topNode);
0074 
0075   nPVs = primaryVertices.size();
0076 
0077   std::vector<int> goodTrackIndex = findAllGoodTracks(daughterParticles, primaryVertices);
0078 
0079   if (!m_has_intermediates)
0080   {
0081     buildBasicChain(selectedMother, selectedVertex, selectedDaughters, daughterParticles, goodTrackIndex, primaryVertices, topNode);
0082   }
0083   else
0084   {
0085     buildChain(selectedMother, selectedVertex, selectedDaughters, selectedIntermediates, daughterParticles, goodTrackIndex, primaryVertices, topNode);
0086   }
0087 }
0088 
0089 /*
0090  *  This function is used to build a basic n-body decay without any intermediate particles such as D's or J/psi's
0091  */
0092 void KFParticle_eventReconstruction::buildBasicChain(std::vector<KFParticle>& selectedMotherBasic,
0093                                                      std::vector<KFParticle>& selectedVertexBasic,
0094                                                      std::vector<std::vector<KFParticle>>& selectedDaughtersBasic,
0095                                                      const std::vector<KFParticle>& daughterParticlesBasic,
0096                                                      const std::vector<int>& goodTrackIndexBasic,
0097                                                      const std::vector<KFParticle>& primaryVerticesBasic, PHCompositeNode* topNode)
0098 {
0099   std::vector<std::vector<int>> goodTracksThatMeet = findTwoProngs(daughterParticlesBasic, goodTrackIndexBasic, m_num_tracks);
0100   for (int p = 3; p < m_num_tracks + 1; ++p)
0101   {
0102     goodTracksThatMeet = findNProngs(daughterParticlesBasic, goodTrackIndexBasic, goodTracksThatMeet, m_num_tracks, p);
0103   }
0104 
0105   getCandidateDecay(selectedMotherBasic, selectedVertexBasic, selectedDaughtersBasic, daughterParticlesBasic,
0106                     goodTracksThatMeet, primaryVerticesBasic, 0, m_num_tracks, false, 0, true, topNode);
0107 }
0108 
0109 /*
0110  *  This function is used to build a more complicated decay with intermediate particles such as D's or J/psi's
0111  */
0112 void KFParticle_eventReconstruction::buildChain(std::vector<KFParticle>& selectedMotherAdv,
0113                                                 std::vector<KFParticle>& selectedVertexAdv,
0114                                                 std::vector<std::vector<KFParticle>>& selectedDaughtersAdv,
0115                                                 std::vector<std::vector<KFParticle>>& selectedIntermediatesAdv,
0116                                                 const std::vector<KFParticle>& daughterParticlesAdv,
0117                                                 const std::vector<int>& goodTrackIndexAdv,
0118                                                 const std::vector<KFParticle>& primaryVerticesAdv, PHCompositeNode* topNode)
0119 {
0120   int track_start = 0;
0121   int track_stop = m_num_tracks_from_intermediate[0];
0122   std::vector<KFParticle> goodCandidates;
0123   std::vector<KFParticle> goodVertex;
0124   std::vector<KFParticle> *goodDaughters = new std::vector<KFParticle>[m_num_tracks];
0125   std::vector<KFParticle> *goodIntermediates = new std::vector<KFParticle>[m_num_intermediate_states];
0126   std::vector<KFParticle> *potentialIntermediates = new std::vector<KFParticle>[m_num_intermediate_states];
0127   std::vector<std::vector<KFParticle>> *potentialDaughters = new std::vector<std::vector<KFParticle>>[m_num_intermediate_states];
0128   for (int i = 0; i < m_num_intermediate_states; ++i)
0129   {
0130     std::vector<KFParticle> vertices;
0131     std::vector<std::vector<int>> goodTracksThatMeet = findTwoProngs(daughterParticlesAdv, goodTrackIndexAdv, m_num_tracks_from_intermediate[i]);
0132     for (int p = 3; p <= m_num_tracks_from_intermediate[i]; ++p)
0133     {
0134       goodTracksThatMeet = findNProngs(daughterParticlesAdv,
0135                                        goodTrackIndexAdv,
0136                                        goodTracksThatMeet,
0137                                        m_num_tracks_from_intermediate[i], p);
0138     }
0139     getCandidateDecay(potentialIntermediates[i], vertices, potentialDaughters[i], daughterParticlesAdv,
0140                       goodTracksThatMeet, primaryVerticesAdv, track_start, track_stop, true, i, m_constrain_int_mass, topNode);
0141     track_start += track_stop;
0142     track_stop += m_num_tracks_from_intermediate[i + 1];
0143   }
0144   int num_tracks_used_by_intermediates = 0;
0145   for (int i = 0; i < m_num_intermediate_states; ++i)
0146   {
0147     num_tracks_used_by_intermediates += m_num_tracks_from_intermediate[i];
0148   }
0149 
0150   int num_remaining_tracks = m_num_tracks - num_tracks_used_by_intermediates;
0151   unsigned int num_pot_inter_a, num_pot_inter_b, num_pot_inter_c, num_pot_inter_d;  // Number of potential intermediates found
0152   num_pot_inter_a = potentialIntermediates[0].size();
0153   num_pot_inter_b = m_num_intermediate_states < 2 ? 1 : potentialIntermediates[1].size();  // Ensure the code inside the loop below is executed
0154   num_pot_inter_c = m_num_intermediate_states < 3 ? 1 : potentialIntermediates[2].size();
0155   num_pot_inter_d = m_num_intermediate_states < 4 ? 1 : potentialIntermediates[3].size();
0156 
0157   for (unsigned int a = 0; a < num_pot_inter_a; ++a)
0158   {
0159     for (unsigned int b = 0; b < num_pot_inter_b; ++b)
0160     {
0161       for (unsigned int c = 0; c < num_pot_inter_c; ++c)
0162       {
0163         for (unsigned int d = 0; d < num_pot_inter_d; ++d)
0164         {
0165           KFParticle candidate;
0166           bool isGood = false;
0167           const unsigned int matchIterators[4] = {a, b, c, d};
0168 
0169           int num_mother_decay_products = m_num_intermediate_states + num_remaining_tracks;
0170           assert(num_mother_decay_products > 0);
0171           KFParticle *motherDecayProducts = new KFParticle[num_mother_decay_products];
0172           std::vector<KFParticle> finalTracks = potentialDaughters[0][a];
0173 
0174           for (int i = 0; i < m_num_intermediate_states; ++i)
0175           {
0176             motherDecayProducts[i] = potentialIntermediates[i][matchIterators[i]];
0177           }
0178           for (int j = 1; j < m_num_intermediate_states; ++j)
0179           {
0180             finalTracks.insert(finalTracks.end(), potentialDaughters[j][matchIterators[j]].begin(), potentialDaughters[j][matchIterators[j]].end());
0181           }
0182 
0183           bool have_duplicate_track = false;
0184 
0185           for (size_t i = 0; i < finalTracks.size(); ++i)
0186           {
0187             for (size_t j = i + 1; j < finalTracks.size(); ++j)
0188             {
0189               if (finalTracks[i].Id() == finalTracks[j].Id())
0190               {
0191                 have_duplicate_track = true;
0192                 break;
0193               }
0194             }
0195             if (have_duplicate_track)
0196             {
0197               break;
0198             }
0199           }
0200 
0201           if (have_duplicate_track)
0202       {
0203         continue;
0204       }
0205 
0206           // If there are daughter tracks coming from the mother not an intermediate, need to ensure that the intermeditate decay tracks aren't used again
0207           std::vector<int> goodTrackIndexAdv_withoutIntermediates = goodTrackIndexAdv;
0208           for (int m = 0; m < num_tracks_used_by_intermediates; ++m)
0209           {
0210             int trackID_to_remove = finalTracks[m].Id();
0211             int trackElement_to_remove = -1;
0212 
0213             auto it = std::find_if(daughterParticlesAdv.begin(), daughterParticlesAdv.end(),
0214                                    [&trackID_to_remove](const KFParticle& obj)
0215                                    { return obj.Id() == trackID_to_remove; });
0216 
0217             if (it != daughterParticlesAdv.end())
0218             {
0219               trackElement_to_remove = std::distance(daughterParticlesAdv.begin(), it);
0220             }
0221 
0222             goodTrackIndexAdv_withoutIntermediates.erase(remove(goodTrackIndexAdv_withoutIntermediates.begin(),
0223                                                                 goodTrackIndexAdv_withoutIntermediates.end(), trackElement_to_remove),
0224                                                          goodTrackIndexAdv_withoutIntermediates.end());
0225           }
0226 
0227           float required_unique_vertexID = 0;
0228           for (int n = 0; n < m_num_intermediate_states; ++n)
0229           {
0230             required_unique_vertexID += m_intermediate_charge[n] * kfp_Tools_evtReco.getParticleMass(m_intermediate_name[n].c_str());
0231           }
0232 
0233           std::vector<std::vector<int>> uniqueCombinations;
0234           std::vector<std::vector<int>> listOfTracksToAppend;
0235 
0236           if (num_remaining_tracks != 0)
0237           {
0238             for (int i = num_tracks_used_by_intermediates; i < m_num_tracks; ++i)
0239             {
0240               required_unique_vertexID += m_daughter_charge[i] * kfp_Tools_evtReco.getParticleMass(m_daughter_name[i].c_str());
0241             }
0242 
0243             uniqueCombinations = findUniqueDaughterCombinations(num_tracks_used_by_intermediates, m_num_tracks);  // Unique comb of remaining trackIDs
0244 
0245             listOfTracksToAppend = appendTracksToIntermediates(motherDecayProducts, daughterParticlesAdv, goodTrackIndexAdv_withoutIntermediates, num_remaining_tracks);
0246 
0247             for (auto& uniqueCombination : uniqueCombinations)
0248             {
0249               for (const auto& element_of_intermediate : m_intermediate_name)
0250               {
0251                 uniqueCombination.insert(begin(uniqueCombination), kfp_Tools_evtReco.getParticleID(element_of_intermediate));
0252               }
0253             }
0254           }
0255           else
0256           {
0257             std::vector<int> m_intermediate_id;
0258             m_intermediate_id.clear();
0259             for (const auto& element_of_intermediate : m_intermediate_name)
0260             {
0261               m_intermediate_id.push_back(kfp_Tools_evtReco.getParticleID(element_of_intermediate));
0262             }
0263             uniqueCombinations.push_back(m_intermediate_id);
0264             listOfTracksToAppend.push_back({0});
0265           }
0266 
0267           for (auto& n_tracks : listOfTracksToAppend)
0268           {
0269             for (int n_trackID = 0; n_trackID < num_remaining_tracks; ++n_trackID)
0270             {
0271               int mDP_trackElem = m_num_intermediate_states + n_trackID;
0272               int dP_trackElem = n_tracks[n_trackID];
0273               motherDecayProducts[mDP_trackElem] = daughterParticlesAdv[dP_trackElem];
0274             }
0275 
0276             for (auto& uniqueCombination : uniqueCombinations)
0277             {
0278               for (const auto& i_pv : primaryVerticesAdv)
0279               {
0280                 std::tie(candidate, isGood) = getCombination(motherDecayProducts, &uniqueCombination[0], i_pv,
0281                                                              m_constrain_to_vertex, false, 0, num_mother_decay_products, m_constrain_int_mass, required_unique_vertexID, topNode);
0282                 if (isGood)
0283                 {
0284 
0285                   if (m_require_bunch_crossing_match)
0286                   {
0287                     KFParticle_truthAndDetTools toolSet;
0288                     std::vector<int> crossings;
0289                     for (int i = 0; i < num_tracks_used_by_intermediates; ++i)
0290                     {
0291                       SvtxTrack *thisTrack = toolSet.getTrack(finalTracks[i].Id(), m_dst_trackmap);
0292                       if (thisTrack)
0293                       {
0294                         crossings.push_back(thisTrack->get_crossing());
0295                       }
0296                     }
0297                 
0298                     for (int k = 0; k < num_remaining_tracks; ++k)
0299                     {
0300                       int trackArrayID = k + m_num_intermediate_states;
0301                       SvtxTrack *thisTrack = toolSet.getTrack(motherDecayProducts[trackArrayID].Id(), m_dst_trackmap);
0302                       if (thisTrack)
0303                       {
0304                         crossings.push_back(thisTrack->get_crossing());
0305                       }
0306                     }
0307                 
0308                     removeDuplicates(crossings);
0309                 
0310                     if (crossings.size() !=1)
0311                     {
0312                       continue;
0313                     }
0314                   }
0315 
0316                   goodCandidates.push_back(candidate);
0317                   if (m_constrain_to_vertex)
0318                   {
0319                     goodVertex.push_back(i_pv);
0320                   }
0321                   for (int k = 0; k < m_num_intermediate_states; ++k)
0322                   {
0323                     goodIntermediates[k].push_back(motherDecayProducts[k]);
0324                   }
0325                   for (int k = 0; k < num_tracks_used_by_intermediates; ++k)
0326                   {
0327                     goodDaughters[k].push_back(finalTracks[k]);
0328                   }
0329                   for (int k = 0; k < num_remaining_tracks; ++k)
0330                   {  // Need to deal with track mass and PID assignment for extra tracks
0331                     int trackArrayID = k + m_num_intermediate_states;
0332                     double slowTrackMass=-1;
0333                     int slowTrackPDG=0;
0334                     if ((Int_t) motherDecayProducts[trackArrayID].GetQ() != 0)
0335                     {
0336                       slowTrackMass = kfp_Tools_evtReco.getParticleMass((Int_t) motherDecayProducts[trackArrayID].GetQ() * uniqueCombination[trackArrayID]);
0337                       slowTrackPDG = (Int_t) motherDecayProducts[trackArrayID].GetQ() * uniqueCombination[trackArrayID];
0338                       // pi+ pdgid = 211, pi- pdgid = -211
0339                       // e+ pdgid = -11, e- pdgid = 11
0340                       // mu+ pdgid = -13, mu- pdgid = 13
0341                       // tau+ pdgid = -15, tau- pdgid = 15
0342                       if (abs(slowTrackPDG) == 11 || abs(slowTrackPDG) == 13 || abs(slowTrackPDG) == 15)
0343                       {
0344                         slowTrackPDG *= -1;
0345                       }
0346                     }
0347                     else if ((Int_t) motherDecayProducts[trackArrayID].GetQ() == 0)
0348                     {
0349                       slowTrackMass = kfp_Tools_evtReco.getParticleMass(uniqueCombination[trackArrayID]);
0350                       slowTrackPDG = uniqueCombination[trackArrayID];
0351                     }
0352                     KFParticle slowTrack;
0353                     slowTrack.Create(motherDecayProducts[trackArrayID].Parameters(),
0354                                      motherDecayProducts[trackArrayID].CovarianceMatrix(),
0355                                      (Int_t) motherDecayProducts[trackArrayID].GetQ(),
0356                                      slowTrackMass);
0357                     slowTrack.NDF() = motherDecayProducts[trackArrayID].GetNDF();
0358                     slowTrack.Chi2() = motherDecayProducts[trackArrayID].GetChi2();
0359                     slowTrack.SetId(motherDecayProducts[trackArrayID].Id());
0360                     slowTrack.SetPDG(slowTrackPDG);
0361                     goodDaughters[k + num_tracks_used_by_intermediates].push_back(slowTrack);
0362                   }
0363                 }
0364               }
0365             }
0366 
0367             if (goodCandidates.size() != 0)
0368             {
0369               int bestCombinationIndex = selectBestCombination(m_constrain_to_vertex, false, goodCandidates, goodVertex);
0370 
0371               selectedMotherAdv.push_back(goodCandidates[bestCombinationIndex]);
0372               if (m_constrain_to_vertex)
0373               {
0374                 selectedVertexAdv.push_back(goodVertex[bestCombinationIndex]);
0375               }
0376               std::vector<KFParticle> intermediates;
0377               intermediates.reserve(m_num_intermediate_states);
0378               for (int i = 0; i < m_num_intermediate_states; ++i)
0379               {
0380                 intermediates.push_back(goodIntermediates[i][bestCombinationIndex]);
0381               }
0382               selectedIntermediatesAdv.push_back(intermediates);
0383               std::vector<KFParticle> particles;
0384               particles.reserve(m_num_tracks);
0385               for (int i = 0; i < m_num_tracks; ++i)
0386               {
0387                 particles.push_back(goodDaughters[i][bestCombinationIndex]);
0388               }
0389               selectedDaughtersAdv.push_back(particles);
0390             }
0391             goodCandidates.clear();
0392             goodVertex.clear();
0393             for (int j = 0; j < m_num_intermediate_states; ++j)
0394             {
0395               goodIntermediates[j].clear();
0396             }
0397             for (int j = 0; j < m_num_tracks; ++j)
0398             {
0399               goodDaughters[j].clear();
0400             }
0401           }
0402       delete [] motherDecayProducts;
0403         }  // Close forth intermediate
0404       }    // Close third intermediate
0405     }      // Close second intermediate
0406   }        // Close first intermediate
0407   delete [] goodDaughters;
0408   delete [] goodIntermediates;
0409   delete [] potentialIntermediates;
0410 }
0411 
0412 void KFParticle_eventReconstruction::getCandidateDecay(std::vector<KFParticle>& selectedMotherCand,
0413                                                        std::vector<KFParticle>& selectedVertexCand,
0414                                                        std::vector<std::vector<KFParticle>>& selectedDaughtersCand,
0415                                                        const std::vector<KFParticle>& daughterParticlesCand,
0416                                                        const std::vector<std::vector<int>>& goodTracksThatMeetCand,
0417                                                        const std::vector<KFParticle>& primaryVerticesCand,
0418                                                        int n_track_start, int n_track_stop,
0419                                                        bool isIntermediate, int intermediateNumber, bool constrainMass, PHCompositeNode* topNode)
0420 {
0421   int nTracks = n_track_stop - n_track_start;
0422   std::vector<std::vector<int>> uniqueCombinations = findUniqueDaughterCombinations(n_track_start, n_track_stop);
0423   std::vector<KFParticle> goodCandidates, goodVertex;
0424   std::vector<KFParticle> *goodDaughters = new std::vector<KFParticle>[nTracks];
0425   KFParticle candidate;
0426   bool isGood;
0427   bool fixToPV = m_constrain_to_vertex && !isIntermediate;
0428 
0429   float required_unique_vertexID = 0;
0430   for (int i = n_track_start; i < n_track_stop; ++i)
0431   {
0432     required_unique_vertexID += m_daughter_charge[i] * kfp_Tools_evtReco.getParticleMass(m_daughter_name[i].c_str());
0433   }
0434 
0435   for (auto& i_comb : goodTracksThatMeetCand)  // Loop over all good track combinations
0436   {
0437     KFParticle *daughterTracks = new KFParticle[nTracks];
0438 
0439     for (int i_track = 0; i_track < nTracks; ++i_track)
0440     {
0441       daughterTracks[i_track] = daughterParticlesCand[i_comb[i_track]];
0442     }  // Build array of the good tracks in that combination
0443 
0444     for (auto& uniqueCombination : uniqueCombinations)  // Loop over unique track PID assignments
0445     {
0446       for (unsigned int i_pv = 0; i_pv < primaryVerticesCand.size(); ++i_pv)  // Loop over all PVs in the event
0447       {
0448         int* PDGIDofFirstParticleInCombination = &uniqueCombination[0];
0449         std::tie(candidate, isGood) = getCombination(daughterTracks, PDGIDofFirstParticleInCombination, primaryVerticesCand[i_pv], m_constrain_to_vertex,
0450                                                      isIntermediate, intermediateNumber, nTracks, constrainMass, required_unique_vertexID, topNode);
0451         if (isIntermediate && isGood)
0452         {
0453           float min_ip = 0;
0454           float min_ipchi2 = 0;
0455           float min_ip_xy = 0;
0456           float min_ipchi2_xy  = 0;
0457           calcMinIP(candidate, primaryVerticesCand, min_ip, min_ipchi2);
0458           calcMinIP(candidate, primaryVerticesCand, min_ip_xy , min_ipchi2_xy, false);
0459           if (!isInRange(m_intermediate_min_ip[intermediateNumber], min_ip, m_intermediate_max_ip[intermediateNumber])
0460            || !isInRange(m_intermediate_min_ipchi2[intermediateNumber], min_ipchi2, m_intermediate_max_ipchi2[intermediateNumber])
0461            || !isInRange(m_intermediate_min_ip_xy[intermediateNumber], min_ip_xy, m_intermediate_max_ip_xy[intermediateNumber])
0462            || !isInRange(m_intermediate_min_ipchi2_xy[intermediateNumber], min_ipchi2_xy, m_intermediate_max_ipchi2_xy[intermediateNumber]))
0463           {
0464             isGood = false;
0465           }
0466         }
0467 
0468         if (isGood)
0469         {
0470           goodCandidates.push_back(candidate);
0471           goodVertex.push_back(primaryVerticesCand[i_pv]);
0472           for (int i = 0; i < nTracks; ++i)
0473           {
0474             KFParticle intParticle;
0475             double intParticleMass=-1;
0476             int intParticlePDG=0;
0477             if ((Int_t) daughterTracks[i].GetQ() != 0)
0478             {
0479               intParticleMass = kfp_Tools_evtReco.getParticleMass((Int_t) daughterTracks[i].GetQ() * PDGIDofFirstParticleInCombination[i]);
0480               intParticlePDG = (Int_t) daughterTracks[i].GetQ() * PDGIDofFirstParticleInCombination[i];
0481               // pi+ pdgid = 211, pi- pdgid = -211
0482               // e+ pdgid = -11, e- pdgid = 11
0483               // mu+ pdgid = -13, mu- pdgid = 13
0484               // tau+ pdgid = -15, tau- pdgid = 15
0485               if (abs(intParticlePDG) == 11 || abs(intParticlePDG) == 13 || abs(intParticlePDG) == 15)
0486               {
0487                 intParticlePDG *= -1;
0488               }
0489             }
0490             else if ((Int_t) daughterTracks[i].GetQ() == 0)
0491             {
0492               intParticleMass = kfp_Tools_evtReco.getParticleMass(PDGIDofFirstParticleInCombination[i]);
0493               intParticlePDG = PDGIDofFirstParticleInCombination[i];
0494             }
0495             intParticle.Create(daughterTracks[i].Parameters(),
0496                                daughterTracks[i].CovarianceMatrix(),
0497                                (Int_t) daughterTracks[i].GetQ(),
0498                                intParticleMass);
0499             intParticle.NDF() = daughterTracks[i].GetNDF();
0500             intParticle.Chi2() = daughterTracks[i].GetChi2();
0501             intParticle.SetId(daughterTracks[i].Id());
0502             intParticle.SetPDG(intParticlePDG);
0503             goodDaughters[i].push_back(intParticle);
0504           }
0505         }
0506       }
0507     }
0508 
0509     if (goodCandidates.size() != 0)
0510     {
0511       int bestCombinationIndex = selectBestCombination(fixToPV, isIntermediate, goodCandidates, goodVertex);
0512 
0513       selectedMotherCand.push_back(goodCandidates[bestCombinationIndex]);
0514       if (fixToPV)
0515       {
0516         selectedVertexCand.push_back(goodVertex[bestCombinationIndex]);
0517       }
0518       std::vector<KFParticle> particles;
0519       particles.reserve(nTracks);
0520       for (int i = 0; i < nTracks; ++i)
0521       {
0522         particles.push_back(goodDaughters[i][bestCombinationIndex]);
0523       }
0524       selectedDaughtersCand.push_back(particles);
0525 
0526       goodCandidates.clear();
0527       goodVertex.clear();
0528       for (int j = 0; j < nTracks; ++j)
0529       {
0530         goodDaughters[j].clear();
0531       }
0532     }
0533     delete [] daughterTracks;
0534   }
0535   delete [] goodDaughters;
0536 }
0537 
0538 int KFParticle_eventReconstruction::selectBestCombination(bool PVconstraint, bool isAnInterMother,
0539                                                           std::vector<KFParticle> possibleCandidates,
0540                                                           std::vector<KFParticle> possibleVertex)
0541 {
0542   KFParticle smallestMassError = possibleCandidates[0];
0543   int bestCombinationIndex = 0;
0544   for (unsigned int i = 1; i < possibleCandidates.size(); ++i)
0545   {
0546     if (PVconstraint && !isAnInterMother)// && !m_require_track_and_vertex_match)
0547     {
0548       float current_IPchi2 = 0;
0549       float best_IPchi2 = 0;
0550    
0551       if (m_use_2D_matching_tools)
0552       {
0553         current_IPchi2 = possibleCandidates[i].GetDeviationFromVertexXY(possibleVertex[i]);
0554         best_IPchi2 = smallestMassError.GetDeviationFromVertexXY(possibleVertex[bestCombinationIndex]);
0555       }
0556       else
0557       {
0558         current_IPchi2 = possibleCandidates[i].GetDeviationFromVertex(possibleVertex[i]);
0559         best_IPchi2 = smallestMassError.GetDeviationFromVertex(possibleVertex[bestCombinationIndex]);
0560       }
0561 
0562       if (current_IPchi2 < best_IPchi2)
0563       {
0564         smallestMassError = possibleCandidates[i];
0565         bestCombinationIndex = i;
0566       }
0567     }
0568     else
0569     {
0570       if (possibleCandidates[i].GetErrMass() < smallestMassError.GetErrMass())
0571       {
0572         smallestMassError = possibleCandidates[i];
0573         bestCombinationIndex = i;
0574       }
0575     }
0576   }
0577   return bestCombinationIndex;
0578 }
0579 
0580 KFParticle KFParticle_eventReconstruction::createFakePV()
0581 {
0582   float f_vertexParameters[6] = {0};
0583 
0584   float f_vertexCovariance[21] = {0};
0585 
0586   KFParticle kfp_vertex;
0587   kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
0588   kfp_vertex.NDF() = 0;
0589   kfp_vertex.Chi2() = 0;
0590   kfp_vertex.SetId(0);
0591   return kfp_vertex;
0592 }
[ 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!