TrackHistory/src/VertexClassifier.cc

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     // Loop over the different GenVertex
0389     for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
0390          ++ivertex) {
0391       bool hasParentVertex = false;
0392
0393       // Loop over the parents looking to see if they are coming from a
0394       // production vertex
0395       for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0396            iparent != (*ivertex)->particles_end(HepMC::parents);
0397            ++iparent)
0398         if ((*iparent)->production_vertex()) {
0399           hasParentVertex = true;
0400           break;
0401         }
0402
0403       // Reject those vertices with parent vertices
0404       if (hasParentVertex)
0405         continue;
0406
0407       // Get the position of the vertex
0408       HepMC::FourVector pos = (*ivertex)->position();
0409
0410       double const mm = 0.1;
0411
0412       GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);
0413
0414       std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
0415
0416       // Search for a VERY close vertex in the list
0417       for (; ientry != genpvs_.end(); ++ientry) {
0418         double distance = sqrt(pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2));
0419         if (distance < vertexClusteringDistance_)
0420           break;
0421       }
0422
0423       // Check if there is not a VERY close vertex and added to the list
0424       if (ientry == genpvs_.end())
0425         ientry = genpvs_.insert(ientry, pv);
0426
0427       // Add the vertex barcodes to the new or existent vertices
0428       ientry->genVertex.push_back((*ivertex)->barcode());
0429
0430       // Collect final state descendants
0431       for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
0432            idecendants != (*ivertex)->particles_end(HepMC::descendants);
0433            ++idecendants) {
0434         if (isFinalstateParticle(*idecendants))
0435           if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
0436               ientry->finalstateParticles.end()) {
0437             ientry->finalstateParticles.push_back((*idecendants)->barcode());
0438             HepMC::FourVector m = (*idecendants)->momentum();
0439
0440             ientry->ptot.setPx(ientry->ptot.px() + m.px());
0441             ientry->ptot.setPy(ientry->ptot.py() + m.py());
0442             ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
0443             ientry->ptot.setE(ientry->ptot.e() + m.e());
0444             ientry->ptsq += m.perp() * m.perp();
0445
0446             if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
0447               ientry->nGenTrk++;
0448           }
0449       }
0450     }
0451   }
0452
0453   std::sort(genpvs_.begin(), genpvs_.end());
0454 }