Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-02-05 23:51:52

0001 
0002 #include <cmath>
0003 #include <cstdlib>
0004 #include <iostream>
0005 
0006 #include "HepPDT/ParticleID.hh"
0007 
0008 #include "SimTracker/TrackHistory/interface/TrackClassifier.h"
0009 
0010 #define update(a, b)  \
0011   do {                \
0012     (a) = (a) || (b); \
0013   } while (0)
0014 
0015 TrackClassifier::TrackClassifier(edm::ParameterSet const &config, edm::ConsumesCollector &&collector)
0016     : TrackCategories(),
0017       hepMCLabel_(config.getUntrackedParameter<edm::InputTag>("hepMC")),
0018       beamSpotLabel_(config.getUntrackedParameter<edm::InputTag>("beamSpot")),
0019       tracer_(config, std::move(collector)),
0020       quality_(config, collector),
0021       magneticFieldToken_(collector.esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0022       particleDataTableToken_(collector.esConsumes()),
0023       transientTrackBuilderToken_(collector.esConsumes<TransientTrackBuilder, TransientTrackRecord>()),
0024       tTopoHandToken_(collector.esConsumes<TrackerTopology, TrackerTopologyRcd>()) {
0025   collector.consumes<edm::HepMCProduct>(hepMCLabel_);
0026   collector.consumes<reco::BeamSpot>(beamSpotLabel_);
0027 
0028   // Set the history depth after hadronization
0029   tracer_.depth(-2);
0030 
0031   // Set the maximum d0pull for the bad category
0032   badPull_ = config.getUntrackedParameter<double>("badPull");
0033 
0034   // Set the minimum decay length for detecting long decays
0035   longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
0036 
0037   // Set the distance for clustering vertices
0038   float vertexClusteringDistance = config.getUntrackedParameter<double>("vertexClusteringDistance");
0039   vertexClusteringSqDistance_ = vertexClusteringDistance * vertexClusteringDistance;
0040 
0041   // Set the number of innermost layers to check for bad hits
0042   numberOfInnerLayers_ = config.getUntrackedParameter<unsigned int>("numberOfInnerLayers");
0043 
0044   // Set the minimum number of simhits in the tracker
0045   minTrackerSimHits_ = config.getUntrackedParameter<unsigned int>("minTrackerSimHits");
0046 }
0047 
0048 void TrackClassifier::newEvent(edm::Event const &event, edm::EventSetup const &setup) {
0049   // Get the new event information for the tracer
0050   tracer_.newEvent(event, setup);
0051 
0052   // Get the new event information for the track quality analyser
0053   quality_.newEvent(event, setup);
0054 
0055   // Get hepmc of the event
0056   event.getByLabel(hepMCLabel_, mcInformation_);
0057 
0058   // Magnetic field
0059   magneticField_ = setup.getHandle(magneticFieldToken_);
0060 
0061   // Get the partivle data table
0062   particleDataTable_ = setup.getHandle(particleDataTableToken_);
0063 
0064   // get the beam spot
0065   event.getByLabel(beamSpotLabel_, beamSpot_);
0066 
0067   // Transient track builder
0068   transientTrackBuilder_ = setup.getHandle(transientTrackBuilderToken_);
0069 
0070   // Create the list of primary vertices associated to the event
0071   genPrimaryVertices();
0072 
0073   // Retrieve tracker topology from geometry
0074   tTopo_ = &setup.getData(tTopoHandToken_);
0075 }
0076 
0077 TrackClassifier const &TrackClassifier::evaluate(reco::TrackBaseRef const &track) {
0078   // Initializing the category vector
0079   reset();
0080 
0081   // Associate and evaluate the track history (check for fakes)
0082   if (tracer_.evaluate(track)) {
0083     // Classify all the tracks by their association and reconstruction
0084     // information
0085     reconstructionInformation(track);
0086 
0087     // Get all the information related to the simulation details
0088     simulationInformation();
0089 
0090     // Analyse the track reconstruction quality
0091     qualityInformation(track);
0092 
0093     // Get hadron flavor of the initial hadron
0094     hadronFlavor();
0095 
0096     // Get all the information related to decay process
0097     processesAtGenerator();
0098 
0099     // Get information about conversion and other interactions
0100     processesAtSimulation();
0101 
0102     // Get geometrical information about the vertices
0103     vertexInformation();
0104 
0105     // Check for unkown classification
0106     unknownTrack();
0107   } else
0108     flags_[Fake] = true;
0109 
0110   return *this;
0111 }
0112 
0113 void TrackClassifier::fillPSetDescription(edm::ParameterSetDescription &desc) {
0114   TrackHistory::fillPSetDescription(desc);
0115 
0116   edm::ParameterSetDescription trackQualityDesc;
0117   TrackQuality::fillPSetDescription(trackQualityDesc);
0118   desc.add<edm::ParameterSetDescription>("hitAssociator", trackQualityDesc);
0119 
0120   desc.addUntracked<edm::InputTag>("hepMC", edm::InputTag("generatorSmeared"));
0121   desc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0122   desc.addUntracked<double>("badPull", 3.0);
0123   desc.addUntracked<double>("longLivedDecayLength", 1e-14);
0124   desc.addUntracked<double>("vertexClusteringDistance", 0.0001);
0125   desc.addUntracked<unsigned int>("numberOfInnerLayers", 2);
0126   desc.addUntracked<unsigned int>("minTrackerSimHits", 3);
0127 }
0128 
0129 TrackClassifier const &TrackClassifier::evaluate(TrackingParticleRef const &track) {
0130   // Initializing the category vector
0131   reset();
0132 
0133   // Trace the history for the given TP
0134   tracer_.evaluate(track);
0135 
0136   // Collect the associated reco track
0137   const reco::TrackBaseRef &recotrack = tracer_.recoTrack();
0138 
0139   // If there is a reco truck then evaluate the simulated history
0140   if (recotrack.isNonnull()) {
0141     flags_[Reconstructed] = true;
0142     // Classify all the tracks by their association and reconstruction
0143     // information
0144     reconstructionInformation(recotrack);
0145     // Analyse the track reconstruction quality
0146     qualityInformation(recotrack);
0147   } else
0148     flags_[Reconstructed] = false;
0149 
0150   // Get all the information related to the simulation details
0151   simulationInformation();
0152 
0153   // Get hadron flavor of the initial hadron
0154   hadronFlavor();
0155 
0156   // Get all the information related to decay process
0157   processesAtGenerator();
0158 
0159   // Get information about conversion and other interactions
0160   processesAtSimulation();
0161 
0162   // Get geometrical information about the vertices
0163   vertexInformation();
0164 
0165   // Check for unkown classification
0166   unknownTrack();
0167 
0168   return *this;
0169 }
0170 
0171 void TrackClassifier::reconstructionInformation(reco::TrackBaseRef const &track) {
0172   TrackingParticleRef tpr = tracer_.simParticle();
0173 
0174   // Compute tracking particle parameters at point of closest approach to the
0175   // beamline
0176 
0177   const SimTrack *assocTrack = &(*tpr->g4Track_begin());
0178 
0179   FreeTrajectoryState ftsAtProduction(
0180       GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0181       GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
0182       TrackCharge(track->charge()),
0183       magneticField_.product());
0184 
0185   try {
0186     TSCPBuilderNoMaterial tscpBuilder;
0187     TrajectoryStateClosestToPoint tsAtClosestApproach =
0188         tscpBuilder(ftsAtProduction, GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()));
0189 
0190     GlobalVector v =
0191         tsAtClosestApproach.theState().position() - GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0());
0192     GlobalVector p = tsAtClosestApproach.theState().momentum();
0193 
0194     // Simulated dxy
0195     double dxySim = -v.x() * sin(p.phi()) + v.y() * cos(p.phi());
0196 
0197     // Simulated dz
0198     double dzSim = v.z() - (v.x() * p.x() + v.y() * p.y()) * p.z() / p.perp2();
0199 
0200     // Calculate the dxy pull
0201     double dxyPull =
0202         std::abs(track->dxy(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dxySim) /
0203         track->dxyError();
0204 
0205     // Calculate the dx pull
0206     double dzPull =
0207         std::abs(track->dz(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dzSim) /
0208         track->dzError();
0209 
0210     // Return true if d0Pull > badD0Pull sigmas
0211     flags_[Bad] = (dxyPull > badPull_ || dzPull > badPull_);
0212 
0213   } catch (cms::Exception const &) {
0214     flags_[Bad] = true;
0215   }
0216 }
0217 
0218 void TrackClassifier::simulationInformation() {
0219   // Get the event id for the initial TP.
0220   EncodedEventId eventId = tracer_.simParticle()->eventId();
0221   // Check for signal events
0222   flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
0223   // Check for muons
0224   flags_[Muon] = (abs(tracer_.simParticle()->pdgId()) == 13);
0225   // Check for the number of psimhit in tracker
0226   flags_[TrackerSimHits] = tracer_.simParticle()->numberOfTrackerLayers() >= (int)minTrackerSimHits_;
0227 }
0228 
0229 void TrackClassifier::qualityInformation(reco::TrackBaseRef const &track) {
0230   // run the hit-by-hit reconstruction quality analysis
0231   quality_.evaluate(tracer_.simParticleTrail(), track, tTopo_);
0232 
0233   unsigned int maxLayers = std::min(numberOfInnerLayers_, quality_.numberOfLayers());
0234 
0235   // check the innermost layers for bad hits
0236   for (unsigned int i = 0; i < maxLayers; i++) {
0237     const TrackQuality::Layer &layer = quality_.layer(i);
0238 
0239     // check all hits in that layer
0240     for (unsigned int j = 0; j < layer.hits.size(); j++) {
0241       const TrackQuality::Layer::Hit &hit = layer.hits[j];
0242 
0243       // In those cases the bad hit was used by track reconstruction
0244       if (hit.state == TrackQuality::Layer::Noise || hit.state == TrackQuality::Layer::Misassoc)
0245         flags_[BadInnerHits] = true;
0246       else if (hit.state == TrackQuality::Layer::Shared)
0247         flags_[SharedInnerHits] = true;
0248     }
0249   }
0250 }
0251 
0252 void TrackClassifier::hadronFlavor() {
0253   // Get the initial hadron from the recoGenParticleTrail
0254   const reco::GenParticle *particle = tracer_.recoGenParticle();
0255 
0256   // Check for the initial hadron
0257   if (particle) {
0258     HepPDT::ParticleID pid(particle->pdgId());
0259     flags_[Bottom] = pid.hasBottom();
0260     flags_[Charm] = pid.hasCharm();
0261     flags_[Light] = !pid.hasCharm() && !pid.hasBottom();
0262   }
0263 }
0264 
0265 void TrackClassifier::processesAtGenerator() {
0266   // pdgid of the "in" particle to the production vertex
0267   int pdgid = 0;
0268 
0269   // Get the generated particles from track history (reco::GenParticle in the
0270   // recoGenParticleTrail)
0271   TrackHistory::RecoGenParticleTrail const &recoGenParticleTrail = tracer_.recoGenParticleTrail();
0272 
0273   // Loop over the generated particles (reco::GenParticle in the
0274   // recoGenParticleTrail)
0275   for (TrackHistory::RecoGenParticleTrail::const_iterator iparticle = recoGenParticleTrail.begin();
0276        iparticle != recoGenParticleTrail.end();
0277        ++iparticle) {
0278     pdgid = std::abs((*iparticle)->pdgId());
0279     // Get particle type
0280     HepPDT::ParticleID particleID(pdgid);
0281 
0282     // Check if the particle type is valid one
0283     if (particleID.isValid()) {
0284       // Get particle data
0285       ParticleData const *particleData = particleDataTable_->particle(particleID);
0286       // Check if the particle exist in the table
0287       if (particleData) {
0288         // Check if their life time is bigger than longLivedDecayLength_
0289         if (particleData->lifetime() > longLivedDecayLength_)
0290           update(flags_[LongLivedDecay], true);
0291         // Check for B and C weak decays
0292         update(flags_[BWeakDecay], particleID.hasBottom());
0293         update(flags_[CWeakDecay], particleID.hasCharm());
0294         // Check for B and C pure leptonic decay
0295         std::set<int> daughterIds;
0296         size_t ndau = (*iparticle)->numberOfDaughters();
0297         for (size_t i = 0; i < ndau; ++i) {
0298           daughterIds.insert((*iparticle)->daughter(i)->pdgId());
0299         }
0300         update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && (daughterIds.find(13) != daughterIds.end()));
0301         update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && (daughterIds.find(13) != daughterIds.end()));
0302       }
0303       // Check Tau, Ks and Lambda decay
0304       update(flags_[ChargePionDecay], pdgid == 211);
0305       update(flags_[ChargeKaonDecay], pdgid == 321);
0306       update(flags_[TauDecay], pdgid == 15);
0307       update(flags_[KsDecay], pdgid == 310);
0308       update(flags_[LambdaDecay], pdgid == 3122);
0309       update(flags_[JpsiDecay], pdgid == 443);
0310       update(flags_[XiDecay], pdgid == 3312);
0311       update(flags_[SigmaPlusDecay], pdgid == 3222);
0312       update(flags_[SigmaMinusDecay], pdgid == 3112);
0313     }
0314   }
0315   // Decays in flight
0316   update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
0317   update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
0318   update(flags_[DecayOnFlightMuon], (flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]));
0319 }
0320 
0321 void TrackClassifier::processesAtSimulation() {
0322   TrackHistory::SimParticleTrail const &simParticleTrail = tracer_.simParticleTrail();
0323 
0324   // Loop over the simulated particles
0325   for (TrackHistory::SimParticleTrail::const_iterator iparticle = simParticleTrail.begin();
0326        iparticle != simParticleTrail.end();
0327        ++iparticle) {
0328     // pdgid of the real source parent vertex
0329     int pdgid = 0;
0330 
0331     // Get a reference to the TP's parent vertex
0332     TrackingVertexRef const &parentVertex = (*iparticle)->parentVertex();
0333 
0334     // Look for the original source track
0335     if (parentVertex.isNonnull()) {
0336       // select the original source in case of combined vertices
0337       bool flag = false;
0338       TrackingVertex::tp_iterator itd, its;
0339 
0340       for (its = parentVertex->sourceTracks_begin(); its != parentVertex->sourceTracks_end(); ++its) {
0341         for (itd = parentVertex->daughterTracks_begin(); itd != parentVertex->daughterTracks_end(); ++itd)
0342           if (itd != its) {
0343             flag = true;
0344             break;
0345           }
0346         if (flag)
0347           break;
0348       }
0349 
0350       // Collect the pdgid of the original source track
0351       if (its != parentVertex->sourceTracks_end())
0352         pdgid = std::abs((*its)->pdgId());
0353       else
0354         pdgid = 0;
0355     }
0356 
0357     unsigned int processG4 = 0;
0358 
0359     // Check existence of SimVerteces assigned
0360     if (parentVertex->nG4Vertices() > 0) {
0361       processG4 = (*(parentVertex->g4Vertices_begin())).processType();
0362     }
0363 
0364     unsigned int process = g4toCMSProcMap_.processId(processG4);
0365 
0366     // Flagging all the different processes
0367     update(flags_[KnownProcess], process != CMS::Undefined && process != CMS::Unknown && process != CMS::Primary);
0368 
0369     update(flags_[UndefinedProcess], process == CMS::Undefined);
0370     update(flags_[UnknownProcess], process == CMS::Unknown);
0371     update(flags_[PrimaryProcess], process == CMS::Primary);
0372     update(flags_[HadronicProcess], process == CMS::Hadronic);
0373     update(flags_[DecayProcess], process == CMS::Decay);
0374     update(flags_[ComptonProcess], process == CMS::Compton);
0375     update(flags_[AnnihilationProcess], process == CMS::Annihilation);
0376     update(flags_[EIoniProcess], process == CMS::EIoni);
0377     update(flags_[HIoniProcess], process == CMS::HIoni);
0378     update(flags_[MuIoniProcess], process == CMS::MuIoni);
0379     update(flags_[PhotonProcess], process == CMS::Photon);
0380     update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
0381     update(flags_[ConversionsProcess], process == CMS::Conversions);
0382     update(flags_[EBremProcess], process == CMS::EBrem);
0383     update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
0384     update(flags_[MuBremProcess], process == CMS::MuBrem);
0385     update(flags_[MuNuclProcess], process == CMS::MuNucl);
0386 
0387     // Get particle type
0388     HepPDT::ParticleID particleID(pdgid);
0389 
0390     // Check if the particle type is valid one
0391     if (particleID.isValid()) {
0392       // Get particle data
0393       ParticleData const *particleData = particleDataTable_->particle(particleID);
0394       // Special treatment for decays
0395       if (process == CMS::Decay) {
0396         // Check if the particle exist in the table
0397         if (particleData) {
0398           // Check if their life time is bigger than 1e-14
0399           if (particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_)
0400             update(flags_[LongLivedDecay], true);
0401 
0402           // Check for B and C weak decays
0403           update(flags_[BWeakDecay], particleID.hasBottom());
0404           update(flags_[CWeakDecay], particleID.hasCharm());
0405 
0406           // Check for B or C pure leptonic decays
0407           int daughtId = abs((*iparticle)->pdgId());
0408           update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughtId == 13);
0409           update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughtId == 13);
0410         }
0411         // Check decays
0412         update(flags_[ChargePionDecay], pdgid == 211);
0413         update(flags_[ChargeKaonDecay], pdgid == 321);
0414         update(flags_[TauDecay], pdgid == 15);
0415         update(flags_[KsDecay], pdgid == 310);
0416         update(flags_[LambdaDecay], pdgid == 3122);
0417         update(flags_[JpsiDecay], pdgid == 443);
0418         update(flags_[XiDecay], pdgid == 3312);
0419         update(flags_[OmegaDecay], pdgid == 3334);
0420         update(flags_[SigmaPlusDecay], pdgid == 3222);
0421         update(flags_[SigmaMinusDecay], pdgid == 3112);
0422       }
0423     }
0424   }
0425   // Decays in flight
0426   update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
0427   update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
0428   update(flags_[DecayOnFlightMuon], flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]);
0429 }
0430 
0431 void TrackClassifier::vertexInformation() {
0432   // Get the main primary vertex from the list
0433   GeneratedPrimaryVertex const &genpv = genpvs_.back();
0434 
0435   // Get the generated history of the tracks
0436   const TrackHistory::GenParticleTrail &genParticleTrail = tracer_.genParticleTrail();
0437 
0438   // Vertex counter
0439   int counter = 0;
0440 
0441   // Unit transformation from mm to cm
0442   double const mm = 0.1;
0443 
0444   double oldX = genpv.x;
0445   double oldY = genpv.y;
0446   double oldZ = genpv.z;
0447 
0448   // Loop over the generated particles
0449   for (TrackHistory::GenParticleTrail::const_reverse_iterator iparticle = genParticleTrail.rbegin();
0450        iparticle != genParticleTrail.rend();
0451        ++iparticle) {
0452     // Look for those with production vertex
0453     HepMC::GenVertex *parent = (*iparticle)->production_vertex();
0454     if (parent) {
0455       HepMC::ThreeVector p = parent->point3d();
0456 
0457       double distance2 = pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2);
0458       double difference2 = pow(p.x() * mm - oldX, 2) + pow(p.y() * mm - oldY, 2) + pow(p.z() * mm - oldZ, 2);
0459 
0460       // std::cout << "Distance2 : " << distance2 << " (" << p.x() * mm << ","
0461       // << p.y() * mm << "," << p.z() * mm << ")" << std::endl; std::cout <<
0462       // "Difference2 : " << difference2 << std::endl;
0463 
0464       if (difference2 > vertexClusteringSqDistance_) {
0465         if (distance2 > vertexClusteringSqDistance_)
0466           counter++;
0467         oldX = p.x() * mm;
0468         oldY = p.y() * mm;
0469         oldZ = p.z() * mm;
0470       }
0471     }
0472   }
0473 
0474   const TrackHistory::SimParticleTrail &simParticleTrail = tracer_.simParticleTrail();
0475 
0476   // Loop over the generated particles
0477   for (TrackHistory::SimParticleTrail::const_reverse_iterator iparticle = simParticleTrail.rbegin();
0478        iparticle != simParticleTrail.rend();
0479        ++iparticle) {
0480     // Look for those with production vertex
0481     TrackingParticle::Point p = (*iparticle)->vertex();
0482 
0483     double distance2 = pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2);
0484     double difference2 = pow(p.x() - oldX, 2) + pow(p.y() - oldY, 2) + pow(p.z() - oldZ, 2);
0485 
0486     // std::cout << "Distance2 : " << distance2 << " (" << p.x() << "," << p.y()
0487     // << "," << p.z() << ")" << std::endl; std::cout << "Difference2 : " <<
0488     // difference2 << std::endl;
0489 
0490     if (difference2 > vertexClusteringSqDistance_) {
0491       if (distance2 > vertexClusteringSqDistance_)
0492         counter++;
0493       oldX = p.x();
0494       oldY = p.y();
0495       oldZ = p.z();
0496     }
0497   }
0498 
0499   if (!counter)
0500     flags_[PrimaryVertex] = true;
0501   else if (counter == 1)
0502     flags_[SecondaryVertex] = true;
0503   else
0504     flags_[TertiaryVertex] = true;
0505 }
0506 
0507 bool TrackClassifier::isFinalstateParticle(const HepMC::GenParticle *p) { return !p->end_vertex() && p->status() == 1; }
0508 
0509 bool TrackClassifier::isCharged(const HepMC::GenParticle *p) {
0510   const ParticleData *part = particleDataTable_->particle(p->pdg_id());
0511   if (part)
0512     return part->charge() != 0;
0513   else {
0514     // the new/improved particle table doesn't know anti-particles
0515     return particleDataTable_->particle(-p->pdg_id()) != nullptr;
0516   }
0517 }
0518 
0519 void TrackClassifier::genPrimaryVertices() {
0520   genpvs_.clear();
0521 
0522   const HepMC::GenEvent *event = mcInformation_->GetEvent();
0523 
0524   if (event) {
0525     // Loop over the different GenVertex
0526     for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
0527          ++ivertex) {
0528       bool hasParentVertex = false;
0529 
0530       // Loop over the parents looking to see if they are coming from a
0531       // production vertex
0532       for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0533            iparent != (*ivertex)->particles_end(HepMC::parents);
0534            ++iparent)
0535         if ((*iparent)->production_vertex()) {
0536           hasParentVertex = true;
0537           break;
0538         }
0539 
0540       // Reject those vertices with parent vertices
0541       if (hasParentVertex)
0542         continue;
0543 
0544       // Get the position of the vertex
0545       HepMC::FourVector pos = (*ivertex)->position();
0546 
0547       double const mm = 0.1;
0548 
0549       GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);
0550 
0551       std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
0552 
0553       // Search for a VERY close vertex in the list
0554       for (; ientry != genpvs_.end(); ++ientry) {
0555         double distance2 = pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2);
0556         if (distance2 < vertexClusteringSqDistance_)
0557           break;
0558       }
0559 
0560       // Check if there is not a VERY close vertex and added to the list
0561       if (ientry == genpvs_.end())
0562         ientry = genpvs_.insert(ientry, pv);
0563 
0564       // Add the vertex barcodes to the new or existent vertices
0565       ientry->genVertex.push_back((*ivertex)->barcode());
0566 
0567       // Collect final state descendants
0568       for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
0569            idecendants != (*ivertex)->particles_end(HepMC::descendants);
0570            ++idecendants) {
0571         if (isFinalstateParticle(*idecendants))
0572           if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
0573               ientry->finalstateParticles.end()) {
0574             ientry->finalstateParticles.push_back((*idecendants)->barcode());
0575             HepMC::FourVector m = (*idecendants)->momentum();
0576 
0577             ientry->ptot.setPx(ientry->ptot.px() + m.px());
0578             ientry->ptot.setPy(ientry->ptot.py() + m.py());
0579             ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
0580             ientry->ptot.setE(ientry->ptot.e() + m.e());
0581             ientry->ptsq += m.perp() * m.perp();
0582 
0583             if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
0584               ientry->nGenTrk++;
0585           }
0586       }
0587     }
0588   }
0589 
0590   std::sort(genpvs_.begin(), genpvs_.end());
0591 }