CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/SimTracker/TrackHistory/src/TrackClassifier.cc

Go to the documentation of this file.
00001 
00002 #include <math.h>
00003 #include <cstdlib>
00004 #include <iostream>
00005 
00006 #include "HepPDT/ParticleID.hh"
00007 
00008 #include "SimTracker/TrackHistory/interface/TrackClassifier.h"
00009 
00010 #define update(a, b) do { (a) = (a) | (b); } while(0)
00011 
00012 TrackClassifier::TrackClassifier(edm::ParameterSet const & config) : TrackCategories(),
00013         hepMCLabel_( config.getUntrackedParameter<edm::InputTag>("hepMC") ),
00014         beamSpotLabel_( config.getUntrackedParameter<edm::InputTag>("beamSpot") ),
00015         tracer_(config),
00016         quality_(config)
00017 {
00018     // Set the history depth after hadronization
00019     tracer_.depth(-2);
00020 
00021     // Set the maximum d0pull for the bad category
00022     badPull_ = config.getUntrackedParameter<double>("badPull");
00023 
00024     // Set the minimum decay length for detecting long decays
00025     longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
00026 
00027     // Set the distance for clustering vertices
00028     float vertexClusteringDistance = config.getUntrackedParameter<double>("vertexClusteringDistance");
00029     vertexClusteringSqDistance_ = vertexClusteringDistance * vertexClusteringDistance;
00030 
00031     // Set the number of innermost layers to check for bad hits
00032     numberOfInnerLayers_ = config.getUntrackedParameter<unsigned int>("numberOfInnerLayers");
00033 
00034     // Set the minimum number of simhits in the tracker
00035     minTrackerSimHits_ = config.getUntrackedParameter<unsigned int>("minTrackerSimHits");
00036 }
00037 
00038 
00039 void TrackClassifier::newEvent ( edm::Event const & event, edm::EventSetup const & setup )
00040 {
00041     // Get the new event information for the tracer
00042     tracer_.newEvent(event, setup);
00043 
00044     // Get the new event information for the track quality analyser
00045     quality_.newEvent(event, setup);
00046 
00047     // Get hepmc of the event
00048     event.getByLabel(hepMCLabel_, mcInformation_);
00049 
00050     // Magnetic field
00051     setup.get<IdealMagneticFieldRecord>().get(magneticField_);
00052 
00053     // Get the partivle data table
00054     setup.getData(particleDataTable_);
00055 
00056     // get the beam spot
00057     event.getByLabel(beamSpotLabel_, beamSpot_);
00058 
00059     // Transient track builder
00060     setup.get<TransientTrackRecord>().get("TransientTrackBuilder", transientTrackBuilder_);
00061 
00062     // Create the list of primary vertices associated to the event
00063     genPrimaryVertices();
00064 
00065     //Retrieve tracker topology from geometry
00066     edm::ESHandle<TrackerTopology> tTopoHand;
00067     setup.get<IdealGeometryRecord>().get(tTopoHand);
00068     tTopo_=tTopoHand.product();
00069 }
00070 
00071 
00072 TrackClassifier const & TrackClassifier::evaluate (reco::TrackBaseRef const & track)
00073 {
00074     // Initializing the category vector
00075     reset();
00076 
00077     // Associate and evaluate the track history (check for fakes)
00078     if ( tracer_.evaluate(track) )
00079     {
00080         // Classify all the tracks by their association and reconstruction information
00081         reconstructionInformation(track);
00082 
00083         // Get all the information related to the simulation details
00084         simulationInformation();
00085 
00086         // Analyse the track reconstruction quality
00087         qualityInformation(track);
00088 
00089         // Get hadron flavor of the initial hadron
00090         hadronFlavor();
00091 
00092         // Get all the information related to decay process
00093         processesAtGenerator();
00094 
00095         // Get information about conversion and other interactions
00096         processesAtSimulation();
00097 
00098         // Get geometrical information about the vertices
00099         vertexInformation();
00100 
00101         // Check for unkown classification
00102         unknownTrack();
00103     }
00104     else
00105         flags_[Fake] = true;
00106 
00107     return *this;
00108 }
00109 
00110 
00111 TrackClassifier const & TrackClassifier::evaluate (TrackingParticleRef const & track)
00112 {
00113     // Initializing the category vector
00114     reset();
00115 
00116     // Trace the history for the given TP
00117     tracer_.evaluate(track);
00118 
00119     // Collect the associated reco track
00120     const reco::TrackBaseRef & recotrack = tracer_.recoTrack();
00121 
00122     // If there is a reco truck then evaluate the simulated history
00123     if ( recotrack.isNonnull() )
00124     {
00125         flags_[Reconstructed] = true;
00126         // Classify all the tracks by their association and reconstruction information
00127         reconstructionInformation(recotrack);
00128         // Analyse the track reconstruction quality
00129         qualityInformation(recotrack);
00130     }
00131     else
00132         flags_[Reconstructed] = false;
00133 
00134     // Get all the information related to the simulation details
00135     simulationInformation();
00136 
00137     // Get hadron flavor of the initial hadron
00138     hadronFlavor();
00139 
00140     // Get all the information related to decay process
00141     processesAtGenerator();
00142 
00143     // Get information about conversion and other interactions
00144     processesAtSimulation();
00145 
00146     // Get geometrical information about the vertices
00147     vertexInformation();
00148 
00149     // Check for unkown classification
00150     unknownTrack();
00151 
00152     return *this;
00153 }
00154 
00155 
00156 void TrackClassifier::reconstructionInformation(reco::TrackBaseRef const & track)
00157 {
00158     TrackingParticleRef tpr = tracer_.simParticle();
00159 
00160     // Compute tracking particle parameters at point of closest approach to the beamline
00161 
00162     const SimTrack * assocTrack = &(*tpr->g4Track_begin());
00163 
00164     FreeTrajectoryState ftsAtProduction(
00165         GlobalPoint(
00166             tpr->vertex().x(),
00167             tpr->vertex().y(),
00168             tpr->vertex().z()
00169         ),
00170         GlobalVector(
00171             assocTrack->momentum().x(),
00172             assocTrack->momentum().y(),
00173             assocTrack->momentum().z()
00174         ),
00175         TrackCharge(track->charge()),
00176         magneticField_.product()
00177     );
00178 
00179     try
00180     {
00181         TSCPBuilderNoMaterial tscpBuilder;
00182         TrajectoryStateClosestToPoint tsAtClosestApproach = tscpBuilder(
00183                     ftsAtProduction,
00184                     GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())
00185                 );
00186 
00187         GlobalVector v = tsAtClosestApproach.theState().position()
00188                          - GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0());
00189         GlobalVector p = tsAtClosestApproach.theState().momentum();
00190 
00191         // Simulated dxy
00192         double dxySim = -v.x()*sin(p.phi()) + v.y()*cos(p.phi());
00193 
00194         // Simulated dz
00195         double dzSim = v.z() - (v.x()*p.x() + v.y()*p.y())*p.z()/p.perp2();
00196 
00197         // Calculate the dxy pull
00198         double dxyPull = std::abs(
00199                              track->dxy( reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()) ) - dxySim
00200                          ) / track->dxyError();
00201 
00202         // Calculate the dx pull
00203         double dzPull = std::abs(
00204                             track->dz( reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()) ) - dzSim
00205                         ) / track->dzError();
00206 
00207         // Return true if d0Pull > badD0Pull sigmas
00208         flags_[Bad] = (dxyPull > badPull_ || dzPull > badPull_);
00209 
00210     }
00211     catch (cms::Exception exception)
00212     {
00213         flags_[Bad] = true;
00214     }
00215 }
00216 
00217 
00218 void TrackClassifier::simulationInformation()
00219 {
00220     // Get the event id for the initial TP.
00221     EncodedEventId eventId = tracer_.simParticle()->eventId();
00222     // Check for signal events
00223     flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
00224     // Check for muons
00225     flags_[Muon] = (abs(tracer_.simParticle()->pdgId()) == 13);
00226     // Check for the number of psimhit in tracker
00227     flags_[TrackerSimHits] = tracer_.simParticle()->numberOfTrackerLayers() >= (int)minTrackerSimHits_;
00228 }
00229 
00230 
00231 void TrackClassifier::qualityInformation(reco::TrackBaseRef const & track)
00232 {
00233     // run the hit-by-hit reconstruction quality analysis
00234   quality_.evaluate(tracer_.simParticleTrail(), track, tTopo_);
00235 
00236     unsigned int maxLayers = std::min(numberOfInnerLayers_, quality_.numberOfLayers());
00237 
00238     // check the innermost layers for bad hits
00239     for (unsigned int i = 0; i < maxLayers; i++)
00240     {
00241         const TrackQuality::Layer &layer = quality_.layer(i);
00242 
00243         // check all hits in that layer
00244         for (unsigned int j = 0; j < layer.hits.size(); j++)
00245         {
00246             const TrackQuality::Layer::Hit &hit = layer.hits[j];
00247 
00248             // In those cases the bad hit was used by track reconstruction
00249             if (hit.state == TrackQuality::Layer::Noise ||
00250                     hit.state == TrackQuality::Layer::Misassoc)
00251                 flags_[BadInnerHits] = true;
00252             else if (hit.state == TrackQuality::Layer::Shared)
00253                 flags_[SharedInnerHits] = true;
00254         }
00255     }
00256 }
00257 
00258 
00259 void TrackClassifier::hadronFlavor()
00260 {
00261     // Get the initial hadron
00262     const HepMC::GenParticle * particle = tracer_.genParticle();
00263 
00264     // Check for the initial hadron
00265     if (particle)
00266     {
00267         HepPDT::ParticleID pid(particle->pdg_id());
00268         flags_[Bottom] = pid.hasBottom();
00269         flags_[Charm] =  pid.hasCharm();
00270         flags_[Light] = !pid.hasCharm() && !pid.hasBottom();
00271     }
00272 }
00273 
00274 
00275 void TrackClassifier::processesAtGenerator()
00276 {
00277     // pdgid of the "in" particle to the production vertex
00278     int pdgid = 0;
00279 
00280     // Get the generated particles from track history
00281     TrackHistory::GenParticleTrail const & genParticleTrail = tracer_.genParticleTrail();
00282 
00283     // Loop over the generated particles
00284     for (TrackHistory::GenParticleTrail::const_iterator iparticle = genParticleTrail.begin(); iparticle != genParticleTrail.end(); ++iparticle)
00285     {
00286         // Get the source vertex for the particle
00287         HepMC::GenVertex * productionVertex = (*iparticle)->production_vertex();
00288 
00289         // Get the pointer to the vertex by removing the const-ness (no const methos in HepMC::GenVertex)
00290         // HepMC::GenVertex * vertex = const_cast<HepMC::GenVertex *>(*ivertex);
00291 
00292         // Check for a non-null pointer to the production vertex
00293         if (productionVertex)
00294         {
00295             // Only case track history will navegate (one in or source particle per vertex)
00296             if ( productionVertex->particles_in_size() == 1 )
00297             {
00298                 // Look at the pdgid of the first "in" particle to the vertex
00299                 pdgid = std::abs((*productionVertex->particles_in_const_begin())->pdg_id());
00300                 // Get particle type
00301                 HepPDT::ParticleID particleID(pdgid);
00302 
00303                 // Check if the particle type is valid one
00304                 if (particleID.isValid())
00305                 {
00306                     // Get particle data
00307                     ParticleData const * particleData = particleDataTable_->particle(particleID);
00308                     // Check if the particle exist in the table
00309                     if (particleData)
00310                     {
00311                         // Check if their life time is bigger than longLivedDecayLength_
00312                         if ( particleData->lifetime() > longLivedDecayLength_ )
00313                             update(flags_[LongLivedDecay], true);
00314                         // Check for B and C weak decays
00315                         update(flags_[BWeakDecay], particleID.hasBottom());
00316                         update(flags_[CWeakDecay], particleID.hasCharm());
00317                         // Check for B and C pure leptonic decay
00318                         int daughterId = abs((*iparticle)->pdg_id());
00319                         update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughterId == 13);
00320                         update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughterId == 13);
00321                     }
00322                     // Check Tau, Ks and Lambda decay
00323                     update(flags_[ChargePionDecay], pdgid == 211);
00324                     update(flags_[ChargeKaonDecay], pdgid == 321);
00325                     update(flags_[TauDecay], pdgid == 15);
00326                     update(flags_[KsDecay], pdgid == 310);
00327                     update(flags_[LambdaDecay], pdgid == 3122);
00328                     update(flags_[JpsiDecay], pdgid == 443);
00329                     update(flags_[XiDecay], pdgid == 3312);
00330                     update(flags_[SigmaPlusDecay], pdgid == 3222);
00331                     update(flags_[SigmaMinusDecay], pdgid == 3112);
00332                 }
00333             }
00334         }
00335     }
00336     // Decays in flight
00337     update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
00338     update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
00339     update(flags_[DecayOnFlightMuon], (flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]));
00340 }
00341 
00342 
00343 void TrackClassifier::processesAtSimulation()
00344 {
00345     TrackHistory::SimParticleTrail const & simParticleTrail = tracer_.simParticleTrail();
00346 
00347     // Loop over the simulated particles
00348     for (
00349         TrackHistory::SimParticleTrail::const_iterator iparticle = simParticleTrail.begin();
00350         iparticle != simParticleTrail.end();
00351         ++iparticle
00352     )
00353     {
00354         // pdgid of the real source parent vertex
00355         int pdgid = 0;
00356 
00357         // Get a reference to the TP's parent vertex
00358         TrackingVertexRef const & parentVertex = (*iparticle)->parentVertex();
00359 
00360         // Look for the original source track
00361         if ( parentVertex.isNonnull() )
00362         {
00363             // select the original source in case of combined vertices
00364             bool flag = false;
00365             TrackingVertex::tp_iterator itd, its;
00366 
00367             for (its = parentVertex->sourceTracks_begin(); its != parentVertex->sourceTracks_end(); ++its)
00368             {
00369                 for (itd = parentVertex->daughterTracks_begin(); itd != parentVertex->daughterTracks_end(); ++itd)
00370                     if (itd != its)
00371                     {
00372                         flag = true;
00373                         break;
00374                     }
00375                 if (flag)
00376                     break;
00377             }
00378 
00379             // Collect the pdgid of the original source track
00380             if ( its != parentVertex->sourceTracks_end() )
00381                 pdgid = std::abs((*its)->pdgId());
00382             else
00383                 pdgid = 0;
00384         }
00385 
00386         // Check for the number of psimhit if different from zero
00387         if ((*iparticle)->numberOfHits()) continue;
00388 
00389         // Collect the G4 process of the first psimhit (it should be the same for all of them)
00390 #warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
00391 #ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
00392         unsigned short process = (*iparticle)->pSimHit_begin()->processType();
00393 #else
00394         unsigned short process = 0;
00395 #endif
00396         // Flagging all the different processes
00397 
00398         update(
00399             flags_[KnownProcess],
00400             process != G4::Undefined &&
00401             process != G4::Unknown &&
00402             process != G4::Primary
00403         );
00404 
00405         update(flags_[UndefinedProcess], process == G4::Undefined);
00406         update(flags_[UnknownProcess], process == G4::Unknown);
00407         update(flags_[PrimaryProcess], process == G4::Primary);
00408         update(flags_[HadronicProcess], process == G4::Hadronic);
00409         update(flags_[DecayProcess], process == G4::Decay);
00410         update(flags_[ComptonProcess], process == G4::Compton);
00411         update(flags_[AnnihilationProcess], process == G4::Annihilation);
00412         update(flags_[EIoniProcess], process == G4::EIoni);
00413         update(flags_[HIoniProcess], process == G4::HIoni);
00414         update(flags_[MuIoniProcess], process == G4::MuIoni);
00415         update(flags_[PhotonProcess], process == G4::Photon);
00416         update(flags_[MuPairProdProcess], process == G4::MuPairProd);
00417         update(flags_[ConversionsProcess], process == G4::Conversions);
00418         update(flags_[EBremProcess], process == G4::EBrem);
00419         update(flags_[SynchrotronRadiationProcess], process == G4::SynchrotronRadiation);
00420         update(flags_[MuBremProcess], process == G4::MuBrem);
00421         update(flags_[MuNuclProcess], process == G4::MuNucl);
00422 
00423         // Get particle type
00424         HepPDT::ParticleID particleID(pdgid);
00425 
00426         // Check if the particle type is valid one
00427         if (particleID.isValid())
00428         {
00429             // Get particle data
00430             ParticleData const * particleData = particleDataTable_->particle(particleID);
00431             // Special treatment for decays
00432             if (process == G4::Decay)
00433             {
00434                 // Check if the particle exist in the table
00435                 if (particleData)
00436                 {
00437                     // Check if their life time is bigger than 1e-14
00438                     if ( particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_ )
00439                         update(flags_[LongLivedDecay], true);
00440 
00441                     // Check for B and C weak decays
00442                     update(flags_[BWeakDecay], particleID.hasBottom());
00443                     update(flags_[CWeakDecay], particleID.hasCharm());
00444 
00445                     // Check for B or C pure leptonic decays
00446                     int daughtId = abs((*iparticle)->pdgId());
00447                     update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughtId == 13);
00448                     update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughtId == 13);
00449                 }
00450                 // Check decays
00451                 update(flags_[ChargePionDecay], pdgid == 211);
00452                 update(flags_[ChargeKaonDecay], pdgid == 321);
00453                 update(flags_[TauDecay], pdgid == 15);
00454                 update(flags_[KsDecay], pdgid == 310);
00455                 update(flags_[LambdaDecay], pdgid == 3122);
00456                 update(flags_[JpsiDecay], pdgid == 443);
00457                 update(flags_[XiDecay], pdgid == 3312);
00458                 update(flags_[OmegaDecay], pdgid == 3334);
00459                 update(flags_[SigmaPlusDecay], pdgid == 3222);
00460                 update(flags_[SigmaMinusDecay], pdgid == 3112);
00461             }
00462         }
00463     }
00464     // Decays in flight
00465     update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
00466     update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
00467     update(flags_[DecayOnFlightMuon], flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]);
00468 }
00469 
00470 
00471 void TrackClassifier::vertexInformation()
00472 {
00473     // Get the main primary vertex from the list
00474     GeneratedPrimaryVertex const & genpv = genpvs_.back();
00475 
00476     // Get the generated history of the tracks
00477     TrackHistory::GenParticleTrail & genParticleTrail = const_cast<TrackHistory::GenParticleTrail &> (tracer_.genParticleTrail());
00478 
00479     // Vertex counter
00480     int counter = 0;
00481 
00482     // Unit transformation from mm to cm
00483     double const mm = 0.1;
00484 
00485     double oldX = genpv.x;
00486     double oldY = genpv.y;
00487     double oldZ = genpv.z;
00488 
00489     // Loop over the generated particles
00490     for (
00491         TrackHistory::GenParticleTrail::reverse_iterator iparticle = genParticleTrail.rbegin();
00492         iparticle != genParticleTrail.rend();
00493         ++iparticle
00494     )
00495     {
00496         // Look for those with production vertex
00497         HepMC::GenVertex * parent = (*iparticle)->production_vertex();
00498         if (parent)
00499         {
00500             HepMC::ThreeVector p = parent->point3d();
00501 
00502             double distance2   = pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2);
00503             double difference2 = pow(p.x() * mm - oldX, 2)    + pow(p.y() * mm - oldY, 2)    + pow(p.z() * mm - oldZ, 2);
00504 
00505             // std::cout << "Distance2 : " << distance2 << " (" << p.x() * mm << "," << p.y() * mm << "," << p.z() * mm << ")" << std::endl;
00506             // std::cout << "Difference2 : " << difference2 << std::endl;
00507 
00508             if ( difference2 > vertexClusteringSqDistance_ )
00509             {
00510                 if ( distance2 > vertexClusteringSqDistance_ ) counter++;
00511                 oldX = p.x() * mm;
00512                 oldY = p.y() * mm;
00513                 oldZ = p.z() * mm;
00514             }
00515         }
00516     }
00517 
00518     TrackHistory::SimParticleTrail & simParticleTrail = const_cast<TrackHistory::SimParticleTrail &> (tracer_.simParticleTrail());
00519 
00520     // Loop over the generated particles
00521     for (
00522         TrackHistory::SimParticleTrail::reverse_iterator iparticle = simParticleTrail.rbegin();
00523         iparticle != simParticleTrail.rend();
00524         ++iparticle
00525     )
00526     {
00527         // Look for those with production vertex
00528         TrackingParticle::Point p = (*iparticle)->vertex();
00529 
00530         double distance2   = pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2);
00531         double difference2 = pow(p.x() - oldX, 2)    + pow(p.y() - oldY, 2)    + pow(p.z() - oldZ, 2);
00532 
00533         // std::cout << "Distance2 : " << distance2 << " (" << p.x() << "," << p.y() << "," << p.z() << ")" << std::endl;
00534         // std::cout << "Difference2 : " << difference2 << std::endl;
00535 
00536         if ( difference2 > vertexClusteringSqDistance_ )
00537         {
00538             if ( distance2 > vertexClusteringSqDistance_ ) counter++;
00539             oldX = p.x();
00540             oldY = p.y();
00541             oldZ = p.z();
00542         }
00543     }
00544 
00545     if ( !counter )
00546         flags_[PrimaryVertex] = true;
00547     else if ( counter == 1 )
00548         flags_[SecondaryVertex] = true;
00549     else
00550         flags_[TertiaryVertex] = true;
00551 }
00552 
00553 
00554 bool TrackClassifier::isFinalstateParticle(const HepMC::GenParticle * p)
00555 {
00556     return !p->end_vertex() && p->status() == 1;
00557 }
00558 
00559 
00560 bool TrackClassifier::isCharged(const HepMC::GenParticle * p)
00561 {
00562     const ParticleData * part = particleDataTable_->particle( p->pdg_id() );
00563     if (part)
00564         return part->charge()!=0;
00565     else
00566     {
00567         // the new/improved particle table doesn't know anti-particles
00568         return  particleDataTable_->particle( -p->pdg_id() ) != 0;
00569     }
00570 }
00571 
00572 
00573 void TrackClassifier::genPrimaryVertices()
00574 {
00575     genpvs_.clear();
00576 
00577     const HepMC::GenEvent * event = mcInformation_->GetEvent();
00578 
00579     if (event)
00580     {
00581         int idx = 0;
00582 
00583         // Loop over the different GenVertex
00584         for ( HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end(); ++ivertex )
00585         {
00586             bool hasParentVertex = false;
00587 
00588             // Loop over the parents looking to see if they are coming from a production vertex
00589             for (
00590                 HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
00591                 iparent != (*ivertex)->particles_end(HepMC::parents);
00592                 ++iparent
00593             )
00594                 if ( (*iparent)->production_vertex() )
00595                 {
00596                     hasParentVertex = true;
00597                     break;
00598                 }
00599 
00600             // Reject those vertices with parent vertices
00601             if (hasParentVertex) continue;
00602 
00603             // Get the position of the vertex
00604             HepMC::FourVector pos = (*ivertex)->position();
00605 
00606             double const mm = 0.1;
00607 
00608             GeneratedPrimaryVertex pv(pos.x()*mm, pos.y()*mm, pos.z()*mm);
00609 
00610             std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
00611 
00612             // Search for a VERY close vertex in the list
00613             for (; ientry != genpvs_.end(); ++ientry)
00614             {
00615                 double distance2 = pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2);
00616                 if ( distance2 < vertexClusteringSqDistance_ )
00617                     break;
00618             }
00619 
00620             // Check if there is not a VERY close vertex and added to the list
00621             if (ientry == genpvs_.end())
00622                 ientry = genpvs_.insert(ientry,pv);
00623 
00624             // Add the vertex barcodes to the new or existent vertices
00625             ientry->genVertex.push_back((*ivertex)->barcode());
00626 
00627             // Collect final state descendants
00628             for (
00629                 HepMC::GenVertex::particle_iterator idecendants  = (*ivertex)->particles_begin(HepMC::descendants);
00630                 idecendants != (*ivertex)->particles_end(HepMC::descendants);
00631                 ++idecendants
00632             )
00633             {
00634                 if (isFinalstateParticle(*idecendants))
00635                     if ( find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) == ientry->finalstateParticles.end() )
00636                     {
00637                         ientry->finalstateParticles.push_back((*idecendants)->barcode());
00638                         HepMC::FourVector m = (*idecendants)->momentum();
00639 
00640                         ientry->ptot.setPx(ientry->ptot.px() + m.px());
00641                         ientry->ptot.setPy(ientry->ptot.py() + m.py());
00642                         ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
00643                         ientry->ptot.setE(ientry->ptot.e() + m.e());
00644                         ientry->ptsq += m.perp() * m.perp();
00645 
00646                         if ( m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants) ) ientry->nGenTrk++;
00647                     }
00648             }
00649             idx++;
00650         }
00651     }
00652 
00653     std::sort(genpvs_.begin(), genpvs_.end());
00654 }