Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-16 23:02:14

0001 /*
0002  *  VertexClassifier.C
0003  */
0004 
0005 #include <cmath>
0006 #include <cstdlib>
0007 #include <iostream>
0008 
0009 #include "HepPDT/ParticleID.hh"
0010 
0011 #include "SimTracker/TrackHistory/interface/VertexClassifier.h"
0012 
0013 #define update(a, b) \
0014   do {               \
0015     (a) = (a) | (b); \
0016   } while (0)
0017 
0018 VertexClassifier::VertexClassifier(edm::ParameterSet const &config, edm::ConsumesCollector collector)
0019     : VertexCategories(),
0020       tracer_(config, collector),
0021       hepMCLabel_(config.getUntrackedParameter<edm::InputTag>("hepMC")),
0022       particleDataTableToken_(collector.esConsumes()) {
0023   collector.consumes<edm::HepMCProduct>(hepMCLabel_);
0024   // Set the history depth after hadronization
0025   tracer_.depth(-2);
0026 
0027   // Set the minimum decay length for detecting long decays
0028   longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
0029 
0030   // Set the distance for clustering vertices
0031   vertexClusteringDistance_ = config.getUntrackedParameter<double>("vertexClusteringDistance");
0032 }
0033 
0034 void VertexClassifier::newEvent(edm::Event const &event, edm::EventSetup const &setup) {
0035   // Get the new event information for the tracer
0036   tracer_.newEvent(event, setup);
0037 
0038   // Get hepmc of the event
0039   event.getByLabel(hepMCLabel_, mcInformation_);
0040 
0041   // Get the partivle data table
0042   particleDataTable_ = setup.getHandle(particleDataTableToken_);
0043 
0044   // Create the list of primary vertices associated to the event
0045   genPrimaryVertices();
0046 }
0047 
0048 VertexClassifier const &VertexClassifier::evaluate(reco::VertexBaseRef const &vertex) {
0049   // Initializing the category vector
0050   reset();
0051 
0052   // Associate and evaluate the vertex history (check for fakes)
0053   if (tracer_.evaluate(vertex)) {
0054     // Get all the information related to the simulation details
0055     simulationInformation();
0056 
0057     // Get all the information related to decay process
0058     processesAtGenerator();
0059 
0060     // Get information about conversion and other interactions
0061     processesAtSimulation();
0062 
0063     // Get geometrical information about the vertices
0064     vertexInformation();
0065 
0066     // Check for unkown classification
0067     unknownVertex();
0068   } else
0069     flags_[Fake] = true;
0070 
0071   return *this;
0072 }
0073 
0074 VertexClassifier const &VertexClassifier::evaluate(TrackingVertexRef const &vertex) {
0075   // Initializing the category vector
0076   reset();
0077 
0078   // Trace the history for the given TP
0079   tracer_.evaluate(vertex);
0080 
0081   // Check for a reconstructed track
0082   if (tracer_.recoVertex().isNonnull())
0083     flags_[Reconstructed] = true;
0084   else
0085     flags_[Reconstructed] = false;
0086 
0087   // Get all the information related to the simulation details
0088   simulationInformation();
0089 
0090   // Get all the information related to decay process
0091   processesAtGenerator();
0092 
0093   // Get information about conversion and other interactions
0094   processesAtSimulation();
0095 
0096   // Get geometrical information about the vertices
0097   vertexInformation();
0098 
0099   // Check for unkown classification
0100   unknownVertex();
0101 
0102   return *this;
0103 }
0104 
0105 void VertexClassifier::simulationInformation() {
0106   // Get the event id for the initial TP.
0107   EncodedEventId eventId = tracer_.simVertex()->eventId();
0108   // Check for signal events
0109   flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
0110 }
0111 
0112 void VertexClassifier::processesAtGenerator() {
0113   // Get the generated vetices from track history
0114   VertexHistory::GenVertexTrail const &genVertexTrail = tracer_.genVertexTrail();
0115 
0116   // Loop over the generated vertices
0117   for (VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin(); ivertex != genVertexTrail.end();
0118        ++ivertex) {
0119     // Get the pointer to the vertex by removing the const-ness (no const methos
0120     // in HepMC::GenVertex)
0121     HepMC::GenVertex *vertex = const_cast<HepMC::GenVertex *>(*ivertex);
0122 
0123     // Loop over the sources looking for specific decays
0124     for (HepMC::GenVertex::particle_iterator iparent = vertex->particles_begin(HepMC::parents);
0125          iparent != vertex->particles_end(HepMC::parents);
0126          ++iparent) {
0127       // Collect the pdgid of the parent
0128       int pdgid = std::abs((*iparent)->pdg_id());
0129       // Get particle type
0130       HepPDT::ParticleID particleID(pdgid);
0131 
0132       // Check if the particle type is valid one
0133       if (particleID.isValid()) {
0134         // Get particle data
0135         ParticleData const *particleData = particleDataTable_->particle(particleID);
0136         // Check if the particle exist in the table
0137         if (particleData) {
0138           // Check if their life time is bigger than longLivedDecayLength_
0139           if (particleData->lifetime() > longLivedDecayLength_) {
0140             // Check for B, C weak decays and long lived decays
0141             update(flags_[BWeakDecay], particleID.hasBottom());
0142             update(flags_[CWeakDecay], particleID.hasCharm());
0143             update(flags_[LongLivedDecay], true);
0144           }
0145           // Check Tau, Ks and Lambda decay
0146           update(flags_[TauDecay], pdgid == 15);
0147           update(flags_[KsDecay], pdgid == 310);
0148           update(flags_[LambdaDecay], pdgid == 3122);
0149           update(flags_[JpsiDecay], pdgid == 443);
0150           update(flags_[XiDecay], pdgid == 3312);
0151           update(flags_[OmegaDecay], pdgid == 3334);
0152           update(flags_[SigmaPlusDecay], pdgid == 3222);
0153           update(flags_[SigmaMinusDecay], pdgid == 3112);
0154         }
0155       }
0156     }
0157   }
0158 }
0159 
0160 void VertexClassifier::processesAtSimulation() {
0161   VertexHistory::SimVertexTrail const &simVertexTrail = tracer_.simVertexTrail();
0162 
0163   for (VertexHistory::SimVertexTrail::const_iterator ivertex = simVertexTrail.begin(); ivertex != simVertexTrail.end();
0164        ++ivertex) {
0165     // pdgid of the real source parent vertex
0166     int pdgid = 0;
0167 
0168     // select the original source in case of combined vertices
0169     bool flag = false;
0170     TrackingVertex::tp_iterator itd, its;
0171 
0172     for (its = (*ivertex)->sourceTracks_begin(); its != (*ivertex)->sourceTracks_end(); ++its) {
0173       for (itd = (*ivertex)->daughterTracks_begin(); itd != (*ivertex)->daughterTracks_end(); ++itd)
0174         if (itd != its) {
0175           flag = true;
0176           break;
0177         }
0178       if (flag)
0179         break;
0180     }
0181     // Collect the pdgid of the original source track
0182     if (its != (*ivertex)->sourceTracks_end())
0183       pdgid = std::abs((*its)->pdgId());
0184     else
0185       pdgid = 0;
0186 
0187     // Geant4 process type is selected using first Geant4 vertex assigned to
0188     // the TrackingVertex
0189     unsigned int processG4 = 0;
0190 
0191     if ((*ivertex)->nG4Vertices() > 0) {
0192       processG4 = (*(*ivertex)->g4Vertices_begin()).processType();
0193     }
0194 
0195     unsigned int process = g4toCMSProcMap_.processId(processG4);
0196 
0197     // Flagging all the different processes
0198     update(flags_[KnownProcess], process != CMS::Undefined && process != CMS::Unknown && process != CMS::Primary);
0199 
0200     update(flags_[UndefinedProcess], process == CMS::Undefined);
0201     update(flags_[UnknownProcess], process == CMS::Unknown);
0202     update(flags_[PrimaryProcess], process == CMS::Primary);
0203     update(flags_[HadronicProcess], process == CMS::Hadronic);
0204     update(flags_[DecayProcess], process == CMS::Decay);
0205     update(flags_[ComptonProcess], process == CMS::Compton);
0206     update(flags_[AnnihilationProcess], process == CMS::Annihilation);
0207     update(flags_[EIoniProcess], process == CMS::EIoni);
0208     update(flags_[HIoniProcess], process == CMS::HIoni);
0209     update(flags_[MuIoniProcess], process == CMS::MuIoni);
0210     update(flags_[PhotonProcess], process == CMS::Photon);
0211     update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
0212     update(flags_[ConversionsProcess], process == CMS::Conversions);
0213     update(flags_[EBremProcess], process == CMS::EBrem);
0214     update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
0215     update(flags_[MuBremProcess], process == CMS::MuBrem);
0216     update(flags_[MuNuclProcess], process == CMS::MuNucl);
0217 
0218     // Loop over the simulated particles
0219     for (TrackingVertex::tp_iterator iparticle = (*ivertex)->daughterTracks_begin();
0220          iparticle != (*ivertex)->daughterTracks_end();
0221          ++iparticle) {
0222       if ((*iparticle)->numberOfTrackerLayers()) {
0223         // Special treatment for decays
0224         if (process == CMS::Decay) {
0225           // Get particle type
0226           HepPDT::ParticleID particleID(pdgid);
0227           // Check if the particle type is valid one
0228           if (particleID.isValid()) {
0229             // Get particle data
0230             ParticleData const *particleData = particleDataTable_->particle(particleID);
0231             // Check if the particle exist in the table
0232             if (particleData) {
0233               // Check if their life time is bigger than 1e-14
0234               if (particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_) {
0235                 // Check for B, C weak decays and long lived decays
0236                 update(flags_[BWeakDecay], particleID.hasBottom());
0237                 update(flags_[CWeakDecay], particleID.hasCharm());
0238                 update(flags_[LongLivedDecay], true);
0239               }
0240               // Check Tau, Ks and Lambda decay
0241               update(flags_[TauDecay], pdgid == 15);
0242               update(flags_[KsDecay], pdgid == 310);
0243               update(flags_[LambdaDecay], pdgid == 3122);
0244               update(flags_[JpsiDecay], pdgid == 443);
0245               update(flags_[XiDecay], pdgid == 3312);
0246               update(flags_[OmegaDecay], pdgid == 3334);
0247               update(flags_[SigmaPlusDecay], pdgid == 3222);
0248               update(flags_[SigmaMinusDecay], pdgid == 3112);
0249             }
0250           }
0251         }
0252       }
0253     }
0254   }
0255 }
0256 
0257 void VertexClassifier::vertexInformation() {
0258   // Helper class for clusterization
0259   typedef std::multimap<double, HepMC::ThreeVector> Clusters;
0260   typedef std::pair<double, HepMC::ThreeVector> ClusterPair;
0261 
0262   Clusters clusters;
0263 
0264   // Get the main primary vertex from the list
0265   GeneratedPrimaryVertex const &genpv = genpvs_.back();
0266 
0267   // Get the generated history of the tracks
0268   const VertexHistory::GenVertexTrail &genVertexTrail = tracer_.genVertexTrail();
0269 
0270   // Unit transformation from mm to cm
0271   double const mm = 0.1;
0272 
0273   // Loop over the generated vertexes
0274   for (VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin(); ivertex != genVertexTrail.end();
0275        ++ivertex) {
0276     // Check vertex exist
0277     if (*ivertex) {
0278       // Measure the distance2 respecto the primary vertex
0279       HepMC::ThreeVector p = (*ivertex)->point3d();
0280       double distance =
0281           sqrt(pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2));
0282 
0283       // If there is not any clusters add the first vertex.
0284       if (clusters.empty()) {
0285         clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0286         continue;
0287       }
0288 
0289       // Check if there is already a cluster in the given distance from primary
0290       // vertex
0291       Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
0292 
0293       if (icluster == clusters.upper_bound(distance + vertexClusteringDistance_)) {
0294         clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0295         continue;
0296       }
0297 
0298       bool cluster = false;
0299 
0300       // Looping over the vertex clusters of a given distance from primary
0301       // vertex
0302       for (; icluster != clusters.upper_bound(distance + vertexClusteringDistance_); ++icluster) {
0303         double difference = sqrt(pow(p.x() * mm - icluster->second.x(), 2) + pow(p.y() * mm - icluster->second.y(), 2) +
0304                                  pow(p.z() * mm - icluster->second.z(), 2));
0305 
0306         if (difference < vertexClusteringDistance_) {
0307           cluster = true;
0308           break;
0309         }
0310       }
0311 
0312       if (!cluster)
0313         clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0314     }
0315   }
0316 
0317   const VertexHistory::SimVertexTrail &simVertexTrail = tracer_.simVertexTrail();
0318 
0319   // Loop over the generated particles
0320   for (VertexHistory::SimVertexTrail::const_reverse_iterator ivertex = simVertexTrail.rbegin();
0321        ivertex != simVertexTrail.rend();
0322        ++ivertex) {
0323     // Look for those with production vertex
0324     TrackingVertex::LorentzVector p = (*ivertex)->position();
0325 
0326     double distance = sqrt(pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2));
0327 
0328     // If there is not any clusters add the first vertex.
0329     if (clusters.empty()) {
0330       clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0331       continue;
0332     }
0333 
0334     // Check if there is already a cluster in the given distance from primary
0335     // vertex
0336     Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
0337 
0338     if (icluster == clusters.upper_bound(distance + vertexClusteringDistance_)) {
0339       clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0340       continue;
0341     }
0342 
0343     bool cluster = false;
0344 
0345     // Looping over the vertex clusters of a given distance from primary vertex
0346     for (; icluster != clusters.upper_bound(distance + vertexClusteringDistance_); ++icluster) {
0347       double difference = sqrt(pow(p.x() - icluster->second.x(), 2) + pow(p.y() - icluster->second.y(), 2) +
0348                                pow(p.z() - icluster->second.z(), 2));
0349 
0350       if (difference < vertexClusteringDistance_) {
0351         cluster = true;
0352         break;
0353       }
0354     }
0355 
0356     if (!cluster)
0357       clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0358   }
0359 
0360   if (clusters.size() == 1)
0361     flags_[PrimaryVertex] = true;
0362   else if (clusters.size() == 2)
0363     flags_[SecondaryVertex] = true;
0364   else
0365     flags_[TertiaryVertex] = true;
0366 }
0367 
0368 bool VertexClassifier::isFinalstateParticle(const HepMC::GenParticle *p) {
0369   return !p->end_vertex() && p->status() == 1;
0370 }
0371 
0372 bool VertexClassifier::isCharged(const HepMC::GenParticle *p) {
0373   const ParticleData *part = particleDataTable_->particle(p->pdg_id());
0374   if (part)
0375     return part->charge() != 0;
0376   else {
0377     // the new/improved particle table doesn't know anti-particles
0378     return particleDataTable_->particle(-p->pdg_id()) != nullptr;
0379   }
0380 }
0381 
0382 void VertexClassifier::genPrimaryVertices() {
0383   genpvs_.clear();
0384 
0385   const HepMC::GenEvent *event = mcInformation_->GetEvent();
0386 
0387   if (event) {
0388     int idx = 0;
0389 
0390     // Loop over the different GenVertex
0391     for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
0392          ++ivertex) {
0393       bool hasParentVertex = false;
0394 
0395       // Loop over the parents looking to see if they are coming from a
0396       // production vertex
0397       for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0398            iparent != (*ivertex)->particles_end(HepMC::parents);
0399            ++iparent)
0400         if ((*iparent)->production_vertex()) {
0401           hasParentVertex = true;
0402           break;
0403         }
0404 
0405       // Reject those vertices with parent vertices
0406       if (hasParentVertex)
0407         continue;
0408 
0409       // Get the position of the vertex
0410       HepMC::FourVector pos = (*ivertex)->position();
0411 
0412       double const mm = 0.1;
0413 
0414       GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);
0415 
0416       std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
0417 
0418       // Search for a VERY close vertex in the list
0419       for (; ientry != genpvs_.end(); ++ientry) {
0420         double distance = sqrt(pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2));
0421         if (distance < vertexClusteringDistance_)
0422           break;
0423       }
0424 
0425       // Check if there is not a VERY close vertex and added to the list
0426       if (ientry == genpvs_.end())
0427         ientry = genpvs_.insert(ientry, pv);
0428 
0429       // Add the vertex barcodes to the new or existent vertices
0430       ientry->genVertex.push_back((*ivertex)->barcode());
0431 
0432       // Collect final state descendants
0433       for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
0434            idecendants != (*ivertex)->particles_end(HepMC::descendants);
0435            ++idecendants) {
0436         if (isFinalstateParticle(*idecendants))
0437           if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
0438               ientry->finalstateParticles.end()) {
0439             ientry->finalstateParticles.push_back((*idecendants)->barcode());
0440             HepMC::FourVector m = (*idecendants)->momentum();
0441 
0442             ientry->ptot.setPx(ientry->ptot.px() + m.px());
0443             ientry->ptot.setPy(ientry->ptot.py() + m.py());
0444             ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
0445             ientry->ptot.setE(ientry->ptot.e() + m.e());
0446             ientry->ptsq += m.perp() * m.perp();
0447 
0448             if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
0449               ientry->nGenTrk++;
0450           }
0451       }
0452       idx++;
0453     }
0454   }
0455 
0456   std::sort(genpvs_.begin(), genpvs_.end());
0457 }