Back to home page

Project CMSSW displayed by LXR

 
 

    


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