TrackHistory/src/HistoryBase.cc

0001 #include <algorithm>
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "SimTracker/TrackHistory/interface/HistoryBase.h"
0006
0007 void HistoryBase::traceRecoGenHistory(reco::GenParticle const *genParticle) {
0008   // fill the trace of reco::GenParticle for correct reading of the History
0009   // flags etc. This is called from the TrackingParticle->genParticle_begin()
0010   // which is a reco::GenParticle
0011
0012   // Take onlt genParticles with a status() smaller than "depth_". Typically
0013   // depth is 2 and so you trace back untill you reach the hard partons see for
0014   // status():
0015   // https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookGenParticleCandidate#GenPCand
0016   if (genParticle->status() <= abs(depth_) && (genParticle->pdgId() < 88 || genParticle->pdgId() > 99)) {
0017     // If the particle is already in the history, it is looping and you should
0018     // stop
0019     if (recoGenParticleTrailHelper_.find(genParticle) != recoGenParticleTrailHelper_.end()) {
0020       return;
0021     }
0022     recoGenParticleTrail_.push_back(genParticle);
0023     recoGenParticleTrailHelper_.insert(genParticle);
0024     // Get the genParticle's mother and trace its history
0025     if (genParticle->mother() != nullptr) {
0026       traceRecoGenHistory((const reco::GenParticle *)genParticle->mother());
0027     }
0028   }
0029 }
0030
0031 void HistoryBase::traceGenHistory(HepMC::GenParticle const *genParticle) {
0032   // Third stop criteria: status abs(depth_) particles after the hadronization.
0033   // The after hadronization is done by detecting the pdg_id pythia code from 88
0034   // to 99
0035   if (genParticle->status() <= abs(depth_) && (genParticle->pdg_id() < 88 || genParticle->pdg_id() > 99)) {
0036     genParticleTrail_.push_back(genParticle);
0037     // Get the producer vertex and trace it history
0038     traceGenHistory(genParticle->production_vertex());
0039   }
0040 }
0041
0042 void HistoryBase::traceGenHistory(HepMC::GenVertex const *genVertex) {
0043   // Verify if has a vertex associated
0044   if (genVertex) {
0045     // Skip if already exist in the collection
0046     if (genVertexTrailHelper_.find(genVertex) != genVertexTrailHelper_.end())
0047       return;
0048     // Add vertex to the history
0049     genVertexTrail_.push_back(genVertex);
0050     genVertexTrailHelper_.insert(genVertex);
0051     // Verify if the vertex has incoming particles
0052     if (genVertex->particles_in_size())
0053       traceGenHistory(*(genVertex->particles_in_const_begin()));
0054   }
0055 }
0056
0057 bool HistoryBase::traceSimHistory(TrackingParticleRef const &trackingParticle, int depth) {
0058   // first stop condition: if the required depth is reached
0059   if (depth == depth_ && depth_ >= 0)
0060     return true;
0061
0062   // second stop condition: if a gen particle is associated to the TP
0063   if (!trackingParticle->genParticles().empty()) {
0064     LogDebug("TrackHistory") << "Particle " << trackingParticle->pdgId() << " has a GenParicle image." << std::endl;
0065
0066     traceRecoGenHistory(&(**(trackingParticle->genParticle_begin())));
0067   }
0068
0069   LogDebug("TrackHistory") << "No GenParticle image for " << trackingParticle->pdgId() << std::endl;
0070
0071   // if no association between TP and genParticles is found: get a reference to
0072   // the TP's parent vertex and trace its history
0073   return traceSimHistory(trackingParticle->parentVertex(), depth);
0074 }
0075
0076 bool HistoryBase::traceSimHistory(TrackingVertexRef const &trackingVertex, int depth) {
0077   // verify if the parent vertex exists
0078   if (trackingVertex.isNonnull()) {
0079     // save the vertex in the trail
0080     simVertexTrail_.push_back(trackingVertex);
0081
0082     if (!trackingVertex->sourceTracks().empty()) {
0083       LogDebug("TrackHistory") << "Moving on to the parent particle." << std::endl;
0084
0085       // select the original source in case of combined vertices
0086       bool flag = false;
0087       TrackingVertex::tp_iterator itd, its;
0088
0089       for (its = trackingVertex->sourceTracks_begin(); its != trackingVertex->sourceTracks_end(); its++) {
0090         for (itd = trackingVertex->daughterTracks_begin(); itd != trackingVertex->daughterTracks_end(); itd++)
0091           if (itd != its) {
0092             flag = true;
0093             break;
0094           }
0095         if (flag)
0096           break;
0097       }
0098
0099       if (!flag)
0100         return false;
0101
0102       // verify if the new particle is not in the trail (looping partiles)
0103       if (std::find(simParticleTrail_.begin(), simParticleTrail_.end(), *its) != simParticleTrail_.end()) {
0104         LogDebug("TrackHistory") << "WARNING: Looping track found." << std::endl;
0105         return false;
0106       }
0107
0108       // save particle in the trail
0109       simParticleTrail_.push_back(*its);
0110       return traceSimHistory(*its, --depth);
0111     } else if (!trackingVertex->genVertices().empty()) {
0112       // navigate over all the associated generated vertexes
0113       LogDebug("TrackHistory") << "Vertex has " << trackingVertex->genVertices().size() << "GenVertex image."
0114                                << std::endl;
0115       for (TrackingVertex::genv_iterator ivertex = trackingVertex->genVertices_begin();
0116            ivertex != trackingVertex->genVertices_end();
0117            ++ivertex)
0118         traceGenHistory(&(**(ivertex)));
0119       return true;
0120     } else {
0121       LogDebug("TrackHistory") << "WARNING: Source track for tracking vertex cannot be found." << std::endl;
0122     }
0123   } else {
0124     LogDebug("TrackHistory") << " WARNING: Vertex cannot be found.";
0125   }
0126
0127   return false;
0128 }