Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:45 

0001 //
0002 //
0003 #include "DataFormats/Common/interface/Ref.h"
0004 
0005 #include "TrackAssociatorByHitsImpl.h"
0006 
0007 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 
0011 //reco track
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0015 //TrackingParticle
0016 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0018 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
0019 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0020 //##---new stuff
0021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0025 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0026 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0027 
0028 using namespace reco;
0029 using namespace std;
0030 
0031 /* Constructor */
0032 
0033 TrackAssociatorByHitsImpl::TrackAssociatorByHitsImpl(edm::EDProductGetter const& productGetter,
0034                                                      std::unique_ptr<TrackerHitAssociator> iAssociate,
0035                                                      TrackerTopology const* iTopo,
0036                                                      SimHitTPAssociationList const* iSimHitsTPAssoc,
0037                                                      SimToRecoDenomType iSimToRecoDenominator,
0038                                                      double iQuality_SimToReco,
0039                                                      double iPurity_SimToReco,
0040                                                      double iCut_RecoToSim,
0041                                                      bool iUsePixels,
0042                                                      bool iUseGrouped,
0043                                                      bool iUseSplitting,
0044                                                      bool iThreeHitTracksAreSpecial,
0045                                                      bool iAbsoluteNumberOfHits)
0046     : productGetter_(&productGetter),
0047       associate(std::move(iAssociate)),
0048       tTopo(iTopo),
0049       simHitsTPAssoc(iSimHitsTPAssoc),
0050       SimToRecoDenominator(iSimToRecoDenominator),
0051       quality_SimToReco(iQuality_SimToReco),
0052       purity_SimToReco(iPurity_SimToReco),
0053       cut_RecoToSim(iCut_RecoToSim),
0054       UsePixels(iUsePixels),
0055       UseGrouped(iUseGrouped),
0056       UseSplitting(iUseSplitting),
0057       ThreeHitTracksAreSpecial(iThreeHitTracksAreSpecial),
0058       AbsoluteNumberOfHits(iAbsoluteNumberOfHits) {}
0059 
0060 /*
0061 TrackAssociatorByHitsImpl::TrackAssociatorByHitsImpl (const edm::ParameterSet& conf) :  
0062   conf_(conf),
0063   AbsoluteNumberOfHits(conf_.getParameter<bool>("AbsoluteNumberOfHits")),
0064   SimToRecoDenominator(denomnone),
0065   quality_SimToReco(conf_.getParameter<double>("Quality_SimToReco")),
0066   purity_SimToReco(conf_.getParameter<double>("Purity_SimToReco")),
0067   cut_RecoToSim(conf_.getParameter<double>("Cut_RecoToSim")),
0068   UsePixels(conf_.getParameter<bool>("UsePixels")),
0069   UseGrouped(conf_.getParameter<bool>("UseGrouped")),
0070   UseSplitting(conf_.getParameter<bool>("UseSplitting")),
0071   ThreeHitTracksAreSpecial(conf_.getParameter<bool>("ThreeHitTracksAreSpecial")),
0072   _simHitTpMapTag(conf_.getParameter<edm::InputTag>("simHitTpMapTag"))
0073 {
0074   std::string tmp = conf_.getParameter<string>("SimToRecoDenominator");
0075   if (tmp=="sim") {
0076     SimToRecoDenominator = denomsim;
0077   } else if (tmp == "reco") {
0078     SimToRecoDenominator = denomreco;
0079   } 
0080 
0081   if (SimToRecoDenominator == denomnone) {
0082     throw cms::Exception("TrackAssociatorByHitsImpl") << "SimToRecoDenominator not specified as sim or reco";
0083   }
0084 
0085 }
0086 */
0087 
0088 //
0089 //---member functions
0090 //
0091 
0092 RecoToSimCollection TrackAssociatorByHitsImpl::associateRecoToSim(
0093     const edm::RefToBaseVector<reco::Track>& tC,
0094     const edm::RefVector<TrackingParticleCollection>& TPCollectionH) const {
0095   //edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateRecoToSim - #tracks="<<tC.size()<<" #TPs="<<TPCollectionH.size();
0096   int nshared = 0;
0097   float quality = 0;  //fraction or absolute number of shared hits
0098   std::vector<SimHitIdpr> SimTrackIds;
0099   std::vector<SimHitIdpr> matchedIds;
0100   RecoToSimCollection outputCollection(productGetter_);
0101 
0102   //dereference the edm::Refs only once
0103   std::vector<TrackingParticle const*> tPC;
0104   tPC.reserve(tPC.size());
0105   for (auto const& ref : TPCollectionH) {
0106     tPC.push_back(&(*ref));
0107   }
0108 
0109   //get the ID of the recotrack  by hits
0110   int tindex = 0;
0111   for (edm::RefToBaseVector<reco::Track>::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0112     matchedIds.clear();
0113     int ri = 0;  //valid rechits
0114     //LogTrace("TrackAssociator") << "\nNEW TRACK - track number " << tindex <<" with pt =" << (*track)->pt() << " # valid=" << (*track)->found();
0115     getMatchedIds<trackingRecHit_iterator>(
0116         matchedIds, SimTrackIds, ri, (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get());
0117 
0118     //LogTrace("TrackAssociator") << "MATCHED IDS LIST BEGIN" ;
0119     //for(size_t j=0; j<matchedIds.size(); j++){
0120     //  LogTrace("TrackAssociator") << "matchedIds[j].first=" << matchedIds[j].first;
0121     //}
0122     //LogTrace("TrackAssociator") << "MATCHED IDS LIST END" ;
0123     //LogTrace("TrackAssociator") << "#matched ids=" << matchedIds.size() << " #tps=" << tPC.size();
0124 
0125     //save id for the track
0126     std::vector<SimHitIdpr> idcachev;
0127     if (!matchedIds.empty()) {
0128       int tpindex = 0;
0129       for (auto t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0130         //int nsimhit = (*t)->trackPSimHit(DetId::Tracker).size();
0131         //LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << (*t)->pdgId() << " with number of PSimHits: "  << nsimhit;
0132         idcachev.clear();
0133         nshared = getShared(matchedIds, idcachev, **t);
0134 
0135         //if electron subtract double counting
0136         if (abs((*t)->pdgId()) == 11 && ((*t)->g4Track_end() - (*t)->g4Track_begin()) > 1) {
0137           nshared -= getDoubleCount<trackingRecHit_iterator>(
0138               (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get(), **t);
0139         }
0140 
0141         if (AbsoluteNumberOfHits)
0142           quality = static_cast<double>(nshared);
0143         else if (ri != 0)
0144           quality = (static_cast<double>(nshared) / static_cast<double>(ri));
0145         else
0146           quality = 0;
0147         //cut on the fraction
0148         //float purity = 1.0*nshared/ri;
0149         if (quality > cut_RecoToSim && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0150           //if a track has just 3 hits we require that all 3 hits are shared
0151           outputCollection.insert(tC[tindex], std::make_pair(TPCollectionH[tpindex], quality));
0152           //          LogTrace("TrackAssociator") << "reco::Track number " << tindex  << " with #hits=" << ri <<" pt=" << (*track)->pt()
0153           //                                      << " associated to TP (pdgId, nb segments, p) = "
0154           //                                      << (*t).pdgId() << " " << (*t).g4Tracks().size()
0155           //                                      << " " << (*t).momentum() << " #hits=" << nsimhit
0156           //                                      << " TP index=" << tpindex<< " with quality =" << quality;
0157         } else {
0158           //          LogTrace("TrackAssociator") <<"reco::Track number " << tindex << " with #hits=" << ri
0159           //                                      << " NOT associated with quality =" << quality;
0160         }
0161       }  //TP loop
0162     }
0163   }
0164   //LogTrace("TrackAssociator") << "% of Assoc Tracks=" << ((double)outputCollection.size())/((double)tC.size());
0165   outputCollection.post_insert();
0166   return outputCollection;
0167 }
0168 
0169 SimToRecoCollection TrackAssociatorByHitsImpl::associateSimToReco(
0170     const edm::RefToBaseVector<reco::Track>& tC,
0171     const edm::RefVector<TrackingParticleCollection>& TPCollectionH) const {
0172   //edm::ESHandle<TrackerTopology> tTopoHand;
0173   //setup->get<IdealGeometryRecord>().get(tTopoHand);
0174 
0175   //  edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateSimToReco - #tracks="<<tC.size()<<" #TPs="<<TPCollectionH.size();
0176   float quality = 0;  //fraction or absolute number of shared hits
0177   int nshared = 0;
0178   std::vector<SimHitIdpr> SimTrackIds;
0179   std::vector<SimHitIdpr> matchedIds;
0180   SimToRecoCollection outputCollection(productGetter_);
0181 
0182   //dereferene the edm::Refs only once
0183   std::vector<TrackingParticle const*> tPC;
0184   tPC.reserve(tPC.size());
0185   for (auto const& ref : TPCollectionH) {
0186     tPC.push_back(&(*ref));
0187   }
0188 
0189   //for (TrackingParticleCollection::const_iterator t = tPC.begin(); t != tPC.end(); ++t) {
0190   //  LogTrace("TrackAssociator") << "NEW TP DUMP";
0191   //  for (TrackingParticle::g4t_iterator g4T = t -> g4Track_begin();g4T !=  t -> g4Track_end(); ++g4T) {
0192   //    LogTrace("TrackAssociator") << "(*g4T).trackId()=" <<(*g4T).trackId() ;
0193   //  }
0194   //}
0195 
0196   //cdj edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0197   //warning: make sure the TP collection used in the map is the same used in the associator!
0198   //e->getByLabel(_simHitTpMapTag,simHitsTPAssoc);
0199 
0200   //get the ID of the recotrack  by hits
0201   int tindex = 0;
0202   for (edm::RefToBaseVector<reco::Track>::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0203     //LogTrace("TrackAssociator") << "\nNEW TRACK - hits of track number " << tindex <<" with pt =" << (*track)->pt() << " # valid=" << (*track)->found();
0204     int ri = 0;  //valid rechits
0205     getMatchedIds<trackingRecHit_iterator>(
0206         matchedIds, SimTrackIds, ri, (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get());
0207 
0208     //save id for the track
0209     std::vector<SimHitIdpr> idcachev;
0210     if (!matchedIds.empty()) {
0211       int tpindex = 0;
0212       for (auto t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0213         idcachev.clear();
0214         float totsimhit = 0;
0215 
0216         //int nsimhit = trackerPSimHit.size();
0217         std::vector<PSimHit> tphits;
0218         //LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << (*t)->pdgId() << " with number of PSimHits: "  << nsimhit;
0219 
0220         nshared = getShared(matchedIds, idcachev, **t);
0221 
0222         //for(std::vector<PSimHit>::const_iterator TPhit = (*t)->trackerPSimHit_begin(); TPhit != (*t)->trackerPSimHit_end(); TPhit++){
0223         //  unsigned int detid = TPhit->detUnitId();
0224         //  DetId detId = DetId(TPhit->detUnitId());
0225         //  LogTrace("TrackAssociator") <<  " hit trackId= " << TPhit->trackId() << " det ID = " << detid
0226         //                              << " SUBDET = " << detId.subdetId() << " layer = " << LayerFromDetid(detId);
0227         //}
0228 
0229         if (nshared != 0) {  //do not waste time recounting when it is not needed!!!!
0230 
0231           //count the TP simhit
0232           //LogTrace("TrackAssociator") << "recounting of tp hits";
0233 
0234           std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0235               TPCollectionH[tpindex], TrackPSimHitRef());  //SimHit is dummy: for simHitTPAssociationListGreater
0236                                                            // sorting only the cluster is needed
0237           auto range = std::equal_range(simHitsTPAssoc->begin(),
0238                                         simHitsTPAssoc->end(),
0239                                         clusterTPpairWithDummyTP,
0240                                         SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0241           for (auto ip = range.first; ip != range.second; ++ip) {
0242             TrackPSimHitRef TPhit = ip->second;
0243             DetId dId = DetId(TPhit->detUnitId());
0244 
0245             unsigned int subdetId = static_cast<unsigned int>(dId.subdetId());
0246             if (!UsePixels && (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap))
0247               continue;
0248 
0249             //unsigned int dRawId = dId.rawId();
0250             SiStripDetId* stripDetId = nullptr;
0251             if (subdetId == SiStripDetId::TIB || subdetId == SiStripDetId::TOB || subdetId == SiStripDetId::TID ||
0252                 subdetId == SiStripDetId::TEC)
0253               stripDetId = new SiStripDetId(dId);
0254             //LogTrace("TrackAssociator") << "consider hit SUBDET = " << subdetId
0255             //                            << " layer = " << LayerFromDetid(dId)
0256             //                            << " id = " << dId.rawId();
0257             bool newhit = true;
0258             for (std::vector<PSimHit>::const_iterator TPhitOK = tphits.begin(); TPhitOK != tphits.end(); TPhitOK++) {
0259               DetId dIdOK = DetId(TPhitOK->detUnitId());
0260               //unsigned int dRawIdOK = dIdOK.rawId();
0261               //LogTrace("TrackAssociator") << "\t\tcompare with SUBDET = " << dIdOK.subdetId()
0262               //                            << " layer = " << LayerFromDetid(dIdOK)
0263               //                            << " id = " << dIdOK.rawId();
0264               //no grouped, no splitting
0265               if (!UseGrouped && !UseSplitting)
0266                 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId())
0267                   newhit = false;
0268               //no grouped, splitting
0269               if (!UseGrouped && UseSplitting)
0270                 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId() &&
0271                     (stripDetId == nullptr || stripDetId->partnerDetId() != dIdOK.rawId()))
0272                   newhit = false;
0273               //grouped, no splitting
0274               if (UseGrouped && !UseSplitting)
0275                 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId() &&
0276                     stripDetId != nullptr && stripDetId->partnerDetId() == dIdOK.rawId())
0277                   newhit = false;
0278               //grouped, splitting
0279               if (UseGrouped && UseSplitting)
0280                 newhit = true;
0281             }
0282             if (newhit) {
0283               //LogTrace("TrackAssociator") << "\t\tok";
0284               tphits.push_back(*TPhit);
0285             }
0286             //else LogTrace("TrackAssociator") << "\t\tno";
0287             delete stripDetId;
0288           }
0289           totsimhit = tphits.size();
0290         }
0291 
0292         if (AbsoluteNumberOfHits)
0293           quality = static_cast<double>(nshared);
0294         else if (SimToRecoDenominator == denomsim && totsimhit != 0)
0295           quality = ((double)nshared) / ((double)totsimhit);
0296         else if (SimToRecoDenominator == denomreco && ri != 0)
0297           quality = ((double)nshared) / ((double)ri);
0298         else
0299           quality = 0;
0300         //LogTrace("TrackAssociator") << "Final count: nhit(TP) = " << nsimhit << " re-counted = " << totsimhit
0301         //<< " nshared = " << nshared << " nrechit = " << ri;
0302 
0303         float purity = 1.0 * nshared / ri;
0304         if (quality > quality_SimToReco && !(ThreeHitTracksAreSpecial && totsimhit == 3 && nshared < 3) &&
0305             (AbsoluteNumberOfHits || (purity > purity_SimToReco))) {
0306           //if a track has just 3 hits we require that all 3 hits are shared
0307           outputCollection.insert(TPCollectionH[tpindex], std::make_pair(tC[tindex], quality));
0308           //          LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0309           //                                      << " re-counted = "  << totsimhit << " nshared = " << nshared
0310           //                                      << " associated to track number " << tindex << " with pt=" << (*track)->pt()
0311           //                                      << " with hit quality =" << quality ;
0312         } else {
0313           //          LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0314           //                                      << " re-counted = "  << totsimhit << " nshared = " << nshared
0315           //                                      << " NOT associated with quality =" << quality;
0316         }
0317       }
0318     }
0319   }
0320   //LogTrace("TrackAssociator") << "% of Assoc TPs=" << ((double)outputCollection.size())/((double)TPCollectionH.size());
0321   outputCollection.post_insert();
0322   return outputCollection;
0323 }
0324 
0325 RecoToSimCollectionSeed TrackAssociatorByHitsImpl::associateRecoToSim(
0326     const edm::Handle<edm::View<TrajectorySeed> >& seedCollectionH,
0327     const edm::Handle<TrackingParticleCollection>& TPCollectionH) const {
0328   edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateRecoToSim - #seeds="
0329                                       << seedCollectionH->size() << " #TPs=" << TPCollectionH->size();
0330   int nshared = 0;
0331   float quality = 0;  //fraction or absolute number of shared hits
0332   std::vector<SimHitIdpr> SimTrackIds;
0333   std::vector<SimHitIdpr> matchedIds;
0334   RecoToSimCollectionSeed outputCollection(productGetter_);
0335 
0336   const TrackingParticleCollection& tPC = *(TPCollectionH.product());
0337 
0338   const edm::View<TrajectorySeed> sC = *(seedCollectionH.product());
0339 
0340   //get the ID of the recotrack  by hits
0341   int tindex = 0;
0342   for (edm::View<TrajectorySeed>::const_iterator seed = sC.begin(); seed != sC.end(); seed++, tindex++) {
0343     matchedIds.clear();
0344     int ri = 0;  //valid rechits
0345     int nsimhit = seed->nHits();
0346     LogTrace("TrackAssociator") << "\nNEW SEED - seed number " << tindex << " # valid=" << nsimhit;
0347     getMatchedIds<edm::OwnVector<TrackingRecHit>::const_iterator>(
0348         matchedIds, SimTrackIds, ri, seed->recHits().begin(), seed->recHits().end(), associate.get());
0349 
0350     //save id for the track
0351     std::vector<SimHitIdpr> idcachev;
0352     if (!matchedIds.empty()) {
0353       int tpindex = 0;
0354       for (TrackingParticleCollection::const_iterator t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0355         LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << t->pdgId()
0356                                     << " with number of PSimHits: " << nsimhit;
0357         idcachev.clear();
0358         nshared = getShared(matchedIds, idcachev, *t);
0359 
0360         //if electron subtract double counting
0361         if (abs(t->pdgId()) == 11 && (t->g4Track_end() - t->g4Track_begin()) > 1) {
0362           nshared -= getDoubleCount<edm::OwnVector<TrackingRecHit>::const_iterator>(
0363               seed->recHits().begin(), seed->recHits().end(), associate.get(), *t);
0364         }
0365 
0366         if (AbsoluteNumberOfHits)
0367           quality = static_cast<double>(nshared);
0368         else if (ri != 0)
0369           quality = (static_cast<double>(nshared) / static_cast<double>(ri));
0370         else
0371           quality = 0;
0372         //cut on the fraction
0373         if (quality > cut_RecoToSim && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0374           //if a track has just 3 hits we require that all 3 hits are shared
0375           outputCollection.insert(
0376               edm::RefToBase<TrajectorySeed>(seedCollectionH, tindex),
0377               std::make_pair(edm::Ref<TrackingParticleCollection>(TPCollectionH, tpindex), quality));
0378           LogTrace("TrackAssociator") << "Seed number " << tindex << " with #hits=" << ri
0379                                       << "associated to TP (pdgId, nb segments, p) = " << (*t).pdgId() << " "
0380                                       << (*t).g4Tracks().size() << " " << (*t).momentum() << " number " << tpindex
0381                                       << " with quality =" << quality;
0382         } else {
0383           LogTrace("TrackAssociator") << "Seed number " << tindex << " with #hits=" << ri
0384                                       << " NOT associated with quality =" << quality;
0385         }
0386       }  //TP loop
0387     }
0388   }
0389   LogTrace("TrackAssociator") << "% of Assoc Seeds="
0390                               << ((double)outputCollection.size()) / ((double)seedCollectionH->size());
0391 
0392   outputCollection.post_insert();
0393   return outputCollection;
0394 }
0395 
0396 SimToRecoCollectionSeed TrackAssociatorByHitsImpl::associateSimToReco(
0397     const edm::Handle<edm::View<TrajectorySeed> >& seedCollectionH,
0398     const edm::Handle<TrackingParticleCollection>& TPCollectionH) const {
0399   edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateSimToReco - #seeds="
0400                                       << seedCollectionH->size() << " #TPs=" << TPCollectionH->size();
0401   float quality = 0;  //fraction or absolute number of shared hits
0402   int nshared = 0;
0403   std::vector<SimHitIdpr> SimTrackIds;
0404   std::vector<SimHitIdpr> matchedIds;
0405   SimToRecoCollectionSeed outputCollection(productGetter_);
0406 
0407   const TrackingParticleCollection& tPC = *TPCollectionH.product();
0408 
0409   const edm::View<TrajectorySeed> sC = *(seedCollectionH.product());
0410 
0411   //get the ID of the recotrack  by hits
0412   int tindex = 0;
0413   for (edm::View<TrajectorySeed>::const_iterator seed = sC.begin(); seed != sC.end(); seed++, tindex++) {
0414     int ri = 0;  //valid rechits
0415     LogTrace("TrackAssociator") << "\nNEW SEED - seed number " << tindex << " # valid=" << seed->nHits();
0416     getMatchedIds<edm::OwnVector<TrackingRecHit>::const_iterator>(
0417         matchedIds, SimTrackIds, ri, seed->recHits().begin(), seed->recHits().end(), associate.get());
0418 
0419     //save id for the track
0420     std::vector<SimHitIdpr> idcachev;
0421     if (!matchedIds.empty()) {
0422       int tpindex = 0;
0423       for (TrackingParticleCollection::const_iterator t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0424         idcachev.clear();
0425         int nsimhit = t->numberOfTrackerHits();
0426         LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << t->pdgId()
0427                                     << " with number of PSimHits: " << nsimhit;
0428         nshared = getShared(matchedIds, idcachev, *t);
0429 
0430         if (AbsoluteNumberOfHits)
0431           quality = static_cast<double>(nshared);
0432         else if (ri != 0)
0433           quality = ((double)nshared) / ((double)ri);
0434         else
0435           quality = 0;
0436         //LogTrace("TrackAssociator") << "Final count: nhit(TP) = " << nsimhit
0437         //<< " nshared = " << nshared
0438         //<< " nrechit = " << ri;
0439         if (quality > quality_SimToReco && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0440           outputCollection.insert(edm::Ref<TrackingParticleCollection>(TPCollectionH, tpindex),
0441                                   std::make_pair(edm::RefToBase<TrajectorySeed>(seedCollectionH, tindex), quality));
0442           LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0443                                       << " associated to seed number " << tindex << " with #hits=" << ri
0444                                       << " with hit quality =" << quality;
0445         } else {
0446           LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0447                                       << " NOT associated with quality =" << quality;
0448         }
0449       }
0450     }
0451   }
0452   LogTrace("TrackAssociator") << "% of Assoc TPs="
0453                               << ((double)outputCollection.size()) / ((double)TPCollectionH->size());
0454   outputCollection.post_insert();
0455   return outputCollection;
0456 }
0457 
0458 template <typename iter>
0459 void TrackAssociatorByHitsImpl::getMatchedIds(std::vector<SimHitIdpr>& matchedIds,
0460                                               std::vector<SimHitIdpr>& SimTrackIds,
0461                                               int& ri,
0462                                               iter begin,
0463                                               iter end,
0464                                               TrackerHitAssociator* associate) const {
0465   matchedIds.clear();
0466   ri = 0;  //valid rechits
0467   for (iter it = begin; it != end; it++) {
0468     const TrackingRecHit* hit = getHitPtr(it);
0469     if (hit->isValid()) {
0470       ri++;
0471       uint32_t t_detID = hit->geographicalId().rawId();
0472       SimTrackIds.clear();
0473       associate->associateHitId(*hit, SimTrackIds);
0474       //save all the id of matched simtracks
0475       if (!SimTrackIds.empty()) {
0476         for (size_t j = 0; j < SimTrackIds.size(); j++) {
0477           LogTrace("TrackAssociator") << " hit # " << ri << " valid=" << hit->isValid() << " det id = " << t_detID
0478                                       << " SimId " << SimTrackIds[j].first << " evt=" << SimTrackIds[j].second.event()
0479                                       << " bc=" << SimTrackIds[j].second.bunchCrossing();
0480           matchedIds.push_back(SimTrackIds[j]);
0481         }
0482       }
0483       ////debugging....****
0484       //****to be used when the crossing frame is in the event and with flag TrackAssociatorByHitsImplESProducer.associateRecoTracks = false
0485       //std::vector<PSimHit> tmpSimHits = associate->associateHit(*getHitPtr(it));
0486       ////cout << "tmpSimHits size=" << tmpSimHits.size() << endl;
0487       //for(size_t j=0; j<tmpSimHits.size(); j++) {
0488       //  LogTrace("TrackAssociator") << " hit # " << ri << " valid=" << getHitPtr(it)->isValid()
0489       //                              << " det id = " << t_detID << " SimId " << SimTrackIds[j].first
0490       //                              << " evt=" << SimTrackIds[j].second.event()
0491       //                              << " bc=" << SimTrackIds[j].second.bunchCrossing()
0492       //                              << " process=" << tmpSimHits[j].processType()
0493       //                              << " eloss=" << tmpSimHits[j].energyLoss()
0494       //                              << " part type=" << tmpSimHits[j].particleType()
0495       //                              << " pabs=" << tmpSimHits[j].pabs()
0496       //                              << " trId=" << tmpSimHits[j].trackId();
0497       //}
0498       ////******************
0499     } else {
0500       LogTrace("TrackAssociator") << "\t\t Invalid Hit On " << hit->geographicalId().rawId();
0501     }
0502   }  //trackingRecHit loop
0503 }
0504 
0505 int TrackAssociatorByHitsImpl::getShared(std::vector<SimHitIdpr>& matchedIds,
0506                                          std::vector<SimHitIdpr>& idcachev,
0507                                          TrackingParticle const& t) const {
0508   int nshared = 0;
0509   if (t.numberOfHits() == 0)
0510     return nshared;  //should use trackerPSimHit but is not const
0511 
0512   for (size_t j = 0; j < matchedIds.size(); j++) {
0513     //LogTrace("TrackAssociator") << "now matchedId=" << matchedIds[j].first;
0514     if (find(idcachev.begin(), idcachev.end(), matchedIds[j]) == idcachev.end()) {
0515       //only the first time we see this ID
0516       idcachev.push_back(matchedIds[j]);
0517 
0518       for (TrackingParticle::g4t_iterator g4T = t.g4Track_begin(); g4T != t.g4Track_end(); ++g4T) {
0519         //  LogTrace("TrackAssociator") << " TP   (ID, Ev, BC) = " << (*g4T).trackId()
0520         //                                  << ", " << t.eventId().event() << ", "<< t.eventId().bunchCrossing()
0521         //                                  << " Match(ID, Ev, BC) = " <<  matchedIds[j].first
0522         //                                  << ", " << matchedIds[j].second.event() << ", "
0523         //                                  << matchedIds[j].second.bunchCrossing() ;
0524         //<< "\t G4  Track Momentum " << (*g4T).momentum()
0525         //<< " \t reco Track Momentum " << track->momentum();
0526         if ((*g4T).trackId() == matchedIds[j].first && t.eventId() == matchedIds[j].second) {
0527           int countedhits = std::count(matchedIds.begin(), matchedIds.end(), matchedIds[j]);
0528           nshared += countedhits;
0529 
0530           //        LogTrace("TrackAssociator") << "hits shared by this segment : " << countedhits;
0531           //        LogTrace("TrackAssociator") << "hits shared so far : " << nshared;
0532         }
0533       }  //g4Tracks loop
0534     }
0535   }
0536   return nshared;
0537 }
0538 
0539 template <typename iter>
0540 int TrackAssociatorByHitsImpl::getDoubleCount(iter begin,
0541                                               iter end,
0542                                               TrackerHitAssociator* associate,
0543                                               TrackingParticle const& t) const {
0544   int doublecount = 0;
0545   std::vector<SimHitIdpr> SimTrackIdsDC;
0546   //  cout<<begin-end<<endl;
0547   for (iter it = begin; it != end; it++) {
0548     int idcount = 0;
0549     SimTrackIdsDC.clear();
0550     associate->associateHitId(*getHitPtr(it), SimTrackIdsDC);
0551     //    cout<<SimTrackIdsDC.size()<<endl;
0552     if (SimTrackIdsDC.size() > 1) {
0553       //     cout<<(t.g4Track_end()-t.g4Track_begin())<<endl;
0554       for (TrackingParticle::g4t_iterator g4T = t.g4Track_begin(); g4T != t.g4Track_end(); ++g4T) {
0555         if (find(SimTrackIdsDC.begin(),
0556                  SimTrackIdsDC.end(),
0557                  SimHitIdpr((*g4T).trackId(), SimTrackIdsDC.begin()->second)) != SimTrackIdsDC.end()) {
0558           idcount++;
0559         }
0560       }
0561     }
0562     if (idcount > 1)
0563       doublecount += (idcount - 1);
0564   }
0565   return doublecount;
0566 }
[ 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!