Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:34:58 

0001 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0002 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0003 #include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 
0006 #include <sstream>
0007 
0008 namespace spr {
0009 
0010   edm::SimTrackContainer::const_iterator matchedSimTrack(const edm::Event& iEvent,
0011                                                          edm::Handle<edm::SimTrackContainer>& SimTk,
0012                                                          edm::Handle<edm::SimVertexContainer>& SimVtx,
0013                                                          const reco::Track* pTrack,
0014                                                          TrackerHitAssociator& associate,
0015                                                          bool debug) {
0016     edm::SimTrackContainer::const_iterator itr = SimTk->end();
0017     ;
0018 
0019     //Get the vector of PsimHits associated to TrackerRecHits and select the
0020     //matching SimTrack on the basis of maximum occurance of trackIds
0021     std::vector<unsigned int> trkId, trkOcc;
0022     for (auto const& trkHit : pTrack->recHits()) {
0023       std::vector<PSimHit> matchedSimIds = associate.associateHit(*trkHit);
0024       for (unsigned int isim = 0; isim < matchedSimIds.size(); isim++) {
0025         unsigned tkId = matchedSimIds[isim].trackId();
0026         bool found = false;
0027         for (unsigned int j = 0; j < trkId.size(); j++) {
0028           if (tkId == trkId[j]) {
0029             trkOcc[j]++;
0030             found = true;
0031             break;
0032           }
0033         }
0034         if (!found) {
0035           trkId.push_back(tkId);
0036           trkOcc.push_back(1);
0037         }
0038       }
0039     }
0040 
0041     if (debug) {
0042       std::ostringstream st1;
0043       for (unsigned int isim = 0; isim < trkId.size(); isim++) {
0044         st1 << "\n trkId " << trkId[isim] << "  Occurance " << trkOcc[isim] << ", ";
0045       }
0046       edm::LogVerbatim("IsoTrack") << st1.str();
0047     }
0048     int matchedId = 0;
0049 
0050     unsigned int matchSimTrk = 0;
0051     if (!trkOcc.empty()) {
0052       unsigned int maxTrkOcc = 0, idxMax = 0;
0053       for (unsigned int j = 0; j < trkOcc.size(); j++) {
0054         if (trkOcc[j] > maxTrkOcc) {
0055           maxTrkOcc = trkOcc[j];
0056           idxMax = j;
0057         }
0058       }
0059       matchSimTrk = trkId[idxMax];
0060       for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0061         if (simTrkItr->trackId() == matchSimTrk) {
0062           matchedId = simTrkItr->type();
0063           if (debug)
0064             edm::LogVerbatim("IsoTrack") << "matched trackId (maximum occurance) " << matchSimTrk << " type "
0065                                          << matchedId;
0066           itr = simTrkItr;
0067           break;
0068         }
0069       }
0070     }
0071 
0072     if (matchedId == 0 && debug) {
0073       edm::LogVerbatim("IsoTrack") << "Could not find matched SimTrk and track history now ";
0074     }
0075     return itr;
0076   }
0077 
0078   //Returns a vector of TrackId originating from the matching SimTrack
0079   std::vector<int> matchedSimTrackId(const edm::Event& iEvent,
0080                                      edm::Handle<edm::SimTrackContainer>& SimTk,
0081                                      edm::Handle<edm::SimVertexContainer>& SimVtx,
0082                                      const reco::Track* pTrack,
0083                                      TrackerHitAssociator& associate,
0084                                      bool debug) {
0085     // get the matching SimTrack
0086     edm::SimTrackContainer::const_iterator trkInfo =
0087         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0088     unsigned int matchSimTrk = trkInfo->trackId();
0089     if (debug)
0090       edm::LogVerbatim("IsoTrack") << "matchedSimTrackId finds the SimTrk ID of the current track to be "
0091                                    << matchSimTrk;
0092     std::vector<int> matchTkid;
0093     if (trkInfo->type() != 0) {
0094       for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0095         if (validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, false))
0096           matchTkid.push_back(static_cast<int>(simTrkItr->trackId()));
0097       }
0098     }
0099     return matchTkid;
0100   }
0101 
0102   spr::simTkInfo matchedSimTrackInfo(unsigned int simTkId,
0103                                      edm::Handle<edm::SimTrackContainer>& SimTk,
0104                                      edm::Handle<edm::SimVertexContainer>& SimVtx,
0105                                      bool debug) {
0106     spr::simTkInfo info;
0107     for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0108       if (simTkId == simTrkItr->trackId()) {
0109         if (spr::validSimTrack(simTkId, simTrkItr, SimTk, SimVtx, debug)) {
0110           info.found = true;
0111           info.pdgId = simTrkItr->type();
0112           info.charge = simTrkItr->charge();
0113         } else {
0114           edm::SimTrackContainer::const_iterator parentItr = spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0115           if (debug) {
0116             if (parentItr != SimTk->end())
0117               edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId() << " "
0118                                            << parentItr->trackId() << ", " << parentItr->type();
0119             else
0120               edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId() << " not found";
0121           }
0122           if (parentItr != SimTk->end()) {
0123             info.found = true;
0124             info.pdgId = parentItr->type();
0125             info.charge = parentItr->charge();
0126           }
0127         }
0128         break;
0129       }
0130     }
0131     return info;
0132   }
0133 
0134   // Checks if this SimTrack=thisTrkItr originates from the one with trackId=simTkId
0135   bool validSimTrack(unsigned int simTkId,
0136                      edm::SimTrackContainer::const_iterator thisTrkItr,
0137                      edm::Handle<edm::SimTrackContainer>& SimTk,
0138                      edm::Handle<edm::SimVertexContainer>& SimVtx,
0139                      bool debug) {
0140     if (debug)
0141       edm::LogVerbatim("IsoTrack") << "Inside validSimTrack: trackId " << thisTrkItr->trackId() << " vtxIndex "
0142                                    << thisTrkItr->vertIndex() << " to be matched to " << simTkId;
0143 
0144     //This track originates from simTkId
0145     if (thisTrkItr->trackId() == simTkId)
0146       return true;
0147 
0148     //Otherwise trace back the history using SimTracks and SimVertices
0149     int vertIndex = thisTrkItr->vertIndex();
0150     if (vertIndex == -1 || vertIndex >= static_cast<int>(SimVtx->size()))
0151       return false;
0152 
0153     edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
0154     for (int iv = 0; iv < vertIndex; iv++)
0155       simVtxItr++;
0156     int parent = simVtxItr->parentIndex();
0157     if (debug)
0158       edm::LogVerbatim("IsoTrack") << "validSimTrack:: parent index " << parent << " ";
0159 
0160     if (parent < 0 && simVtxItr != SimVtx->begin()) {
0161       const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
0162       for (simVtxItr = SimVtx->begin(); simVtxItr != SimVtx->end(); ++simVtxItr) {
0163         if (simVtxItr->parentIndex() > 0) {
0164           const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
0165           double dist = pos2.P();
0166           if (dist < 0.001) {
0167             parent = simVtxItr->parentIndex();
0168             break;
0169           }
0170         }
0171       }
0172     }
0173 
0174     if (debug)
0175       edm::LogVerbatim("IsoTrack") << "final index " << parent;
0176     for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0177       if (static_cast<int>(simTrkItr->trackId()) == parent && simTrkItr != thisTrkItr)
0178         return validSimTrack(simTkId, simTrkItr, SimTk, SimVtx, debug);
0179     }
0180 
0181     return false;
0182   }
0183 
0184   //Returns the parent of a SimTrack
0185   edm::SimTrackContainer::const_iterator parentSimTrack(edm::SimTrackContainer::const_iterator thisTrkItr,
0186                                                         edm::Handle<edm::SimTrackContainer>& SimTk,
0187                                                         edm::Handle<edm::SimVertexContainer>& SimVtx,
0188                                                         bool debug) {
0189     edm::SimTrackContainer::const_iterator itr = SimTk->end();
0190 
0191     int vertIndex = thisTrkItr->vertIndex();
0192     if (debug)
0193       edm::LogVerbatim("IsoTrack") << "SimTrackParent " << thisTrkItr->trackId() << " Vertex " << vertIndex << " Type "
0194                                    << thisTrkItr->type() << " Charge " << static_cast<int>(thisTrkItr->charge());
0195     if (vertIndex == -1)
0196       return thisTrkItr;
0197     else if (vertIndex >= static_cast<int>(SimVtx->size()))
0198       return itr;
0199 
0200     edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
0201     for (int iv = 0; iv < vertIndex; iv++)
0202       simVtxItr++;
0203     int parent = simVtxItr->parentIndex();
0204 
0205     if (parent < 0 && simVtxItr != SimVtx->begin()) {
0206       const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
0207       for (simVtxItr = SimVtx->begin(); simVtxItr != SimVtx->end(); ++simVtxItr) {
0208         if (simVtxItr->parentIndex() > 0) {
0209           const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
0210           double dist = pos2.P();
0211           if (dist < 0.001) {
0212             parent = simVtxItr->parentIndex();
0213             break;
0214           }
0215         }
0216       }
0217     }
0218     for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0219       if (static_cast<int>(simTrkItr->trackId()) == parent && simTrkItr != thisTrkItr)
0220         return parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0221     }
0222 
0223     return thisTrkItr;
0224   }
0225 
0226 }  // namespace spr
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!