Get track history and classify it in function of their . More...

#include <VertexClassifier.h>

Inheritance diagram for VertexClassifier:

Classes
struct	GeneratedPrimaryVertex
	Auxiliary class holding simulated primary vertices. More...

Public Types
typedef VertexCategories	Categories
	Type to the associate category. More...

Public Types inherited from VertexCategories
enum	Category { Fake = 0, Reconstructed = Fake, SignalEvent, BWeakDecay, CWeakDecay, TauDecay, KsDecay, LambdaDecay, JpsiDecay, XiDecay, OmegaDecay, SigmaPlusDecay, SigmaMinusDecay, LongLivedDecay, KnownProcess, UndefinedProcess, UnknownProcess, PrimaryProcess, HadronicProcess, DecayProcess, ComptonProcess, AnnihilationProcess, EIoniProcess, HIoniProcess, MuIoniProcess, PhotonProcess, MuPairProdProcess, ConversionsProcess, EBremProcess, SynchrotronRadiationProcess, MuBremProcess, MuNuclProcess, PrimaryVertex, SecondaryVertex, TertiaryVertex, TierciaryVertex = TertiaryVertex, Unknown }
	Categories available to vertexes. More...

typedef std::vector< bool >	Flags
	Main types associated to the class. More...

Public Member Functions
VertexClassifier const &	evaluate (reco::VertexBaseRef const &)
	Classify the RecoVertex in categories. More...

VertexClassifier const &	evaluate (TrackingVertexRef const &)
	Classify the TrackingVertex in categories. More...

VertexClassifier const &	evaluate (reco::VertexRef const &vertex)
	Classify the RecoVertex in categories. More...

VertexHistory const &	history () const
	Returns a reference to the vertex history used in the classification. More...

virtual void	newEvent (edm::Event const &, edm::EventSetup const &)
	Pre-process event information (for accessing reconstraction information) More...

	VertexClassifier (edm::ParameterSet const &pset, edm::ConsumesCollector &&)
	Constructor by ParameterSet. More...

virtual	~VertexClassifier ()

Public Member Functions inherited from VertexCategories
const Flags &	flags () const
	Returns flags with the category descriptions. More...

bool	is (Category category) const
	Returns track flag for a given category. More...

	VertexCategories ()
	Void constructor. More...

Private Member Functions
void	genPrimaryVertices ()

bool	isCharged (const HepMC::GenParticle *)

bool	isFinalstateParticle (const HepMC::GenParticle *)

void	processesAtGenerator ()
	Get all the information related to decay process. More...

void	processesAtSimulation ()
	Get information about conversion and other interactions. More...

void	reconstructionInformation (reco::TrackBaseRef const &)
	Get reconstruction information. More...

void	simulationInformation ()
	Get all the information related to the simulation details. More...

void	vertexInformation ()
	Get geometrical information about the vertices. More...

Private Attributes
const G4toCMSLegacyProcTypeMap	g4toCMSProcMap_

std::vector < GeneratedPrimaryVertex >	genpvs_

const edm::InputTag	hepMCLabel_

double	longLivedDecayLength_

edm::Handle< edm::HepMCProduct >	mcInformation_

edm::ESHandle< ParticleDataTable >	particleDataTable_

VertexHistory	tracer_

double	vertexClusteringDistance_

Additional Inherited Members
Static Public Attributes inherited from VertexCategories
static const char *const	Names []
	Name of the different categories. More...

Protected Member Functions inherited from VertexCategories
void	reset ()
	Reset the categories flags. More...

void	unknownVertex ()

Protected Attributes inherited from VertexCategories
Flags	flags_
	Flag containers. More...

Detailed Description

Get track history and classify it in function of their .

Definition at line 18 of file VertexClassifier.h.

Member Typedef Documentation

typedef VertexCategories VertexClassifier::Categories

Type to the associate category.

Definition at line 24 of file VertexClassifier.h.

Constructor & Destructor Documentation

VertexClassifier::VertexClassifier	(	edm::ParameterSet const &	pset,
		edm::ConsumesCollector &&	collector
	)

Constructor by ParameterSet.

Definition at line 17 of file VertexClassifier.cc.

References HistoryBase::depth(), edm::ParameterSet::getUntrackedParameter(), hepMCLabel_, longLivedDecayLength_, tracer_, and vertexClusteringDistance_.

                                                                      : 
         VertexCategories(), 
         tracer_(config,std::move(collector)),
         hepMCLabel_( config.getUntrackedParameter<edm::InputTag>("hepMC") )
 {
     collector.consumes<edm::HepMCProduct>(hepMCLabel_);
     // Set the history depth after hadronization
     tracer_.depth(-2);
 
     // Set the minimum decay length for detecting long decays
     longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
 
     // Set the distance for clustering vertices
     vertexClusteringDistance_ = config.getUntrackedParameter<double>("vertexClusteringDistance");
 }

virtual VertexClassifier::~VertexClassifier ( )

inlinevirtual

Definition at line 30 of file VertexClassifier.h.

30 {}

Member Function Documentation

VertexClassifier const & VertexClassifier::evaluate ( reco::VertexBaseRef const & vertex )

Classify the RecoVertex in categories.

Definition at line 51 of file VertexClassifier.cc.

References VertexHistory::evaluate(), VertexCategories::Fake, VertexCategories::flags_, processesAtGenerator(), processesAtSimulation(), VertexCategories::reset(), simulationInformation(), tracer_, VertexCategories::unknownVertex(), and vertexInformation().

Referenced by VertexHistoryAnalyzer::analyze(), VertexClassifierByProxy< reco::SecondaryVertexTagInfoCollection >::evaluate(), and evaluate().

 {
     // Initializing the category vector
     reset();
 
     // Associate and evaluate the vertex history (check for fakes)
     if ( tracer_.evaluate(vertex) )
     {
         // Get all the information related to the simulation details
         simulationInformation();
 
         // Get all the information related to decay process
         processesAtGenerator();
 
         // Get information about conversion and other interactions
         processesAtSimulation();
 
         // Get geometrical information about the vertices
         vertexInformation();
 
         // Check for unkown classification
         unknownVertex();
     }
     else
         flags_[Fake] = true;
 
     return *this;
 }

VertexClassifier const & VertexClassifier::evaluate ( TrackingVertexRef const & vertex )

Classify the TrackingVertex in categories.

Definition at line 81 of file VertexClassifier.cc.

References VertexHistory::evaluate(), VertexCategories::flags_, edm::RefToBase< T >::isNonnull(), processesAtGenerator(), processesAtSimulation(), VertexCategories::Reconstructed, VertexHistory::recoVertex(), VertexCategories::reset(), simulationInformation(), tracer_, VertexCategories::unknownVertex(), and vertexInformation().

 {
     // Initializing the category vector
     reset();
 
     // Trace the history for the given TP
     tracer_.evaluate(vertex);
 
     // Check for a reconstructed track
     if ( tracer_.recoVertex().isNonnull() )
         flags_[Reconstructed] = true;
     else
         flags_[Reconstructed] = false;
 
     // Get all the information related to the simulation details
     simulationInformation();
 
     // Get all the information related to decay process
     processesAtGenerator();
 
     // Get information about conversion and other interactions
     processesAtSimulation();
 
     // Get geometrical information about the vertices
     vertexInformation();
 
     // Check for unkown classification
     unknownVertex();
 
     return *this;
 }

VertexClassifier const& VertexClassifier::evaluate ( reco::VertexRef const & vertex )

inline

Classify the RecoVertex in categories.

Definition at line 42 of file VertexClassifier.h.

References evaluate().

     {
         return evaluate( reco::VertexBaseRef(vertex) );
     }

void VertexClassifier::genPrimaryVertices ( )

private

Definition at line 457 of file VertexClassifier.cc.

References funct::abs(), HLT_25ns14e33_v1_cff::distance, event(), spr::find(), genpvs_, customizeTrackingMonitorSeedNumber::idx, isCharged(), isFinalstateParticle(), visualization-live-secondInstance_cfg::m, mcInformation_, parents, funct::pow(), MetAnalyzer::pv(), python.multivaluedict::sort(), mathSSE::sqrt(), and vertexClusteringDistance_.

Referenced by newEvent().

 {
     genpvs_.clear();
 
     const HepMC::GenEvent * event = mcInformation_->GetEvent();
 
     if (event)
     {
         int idx = 0;
 
         // Loop over the different GenVertex
         for ( HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end(); ++ivertex )
         {
             bool hasParentVertex = false;
 
             // Loop over the parents looking to see if they are coming from a production vertex
             for (
                 HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
                 iparent != (*ivertex)->particles_end(HepMC::parents);
                 ++iparent
             )
                 if ( (*iparent)->production_vertex() )
                 {
                     hasParentVertex = true;
                     break;
                 }
 
             // Reject those vertices with parent vertices
             if (hasParentVertex) continue;
 
             // Get the position of the vertex
             HepMC::FourVector pos = (*ivertex)->position();
 
             double const mm = 0.1;
 
             GeneratedPrimaryVertex pv(pos.x()*mm, pos.y()*mm, pos.z()*mm);
 
             std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
 
             // Search for a VERY close vertex in the list
             for (; ientry != genpvs_.end(); ++ientry)
             {
                 double distance = sqrt( pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2) );
                 if ( distance < vertexClusteringDistance_ )
                     break;
             }
 
             // Check if there is not a VERY close vertex and added to the list
             if (ientry == genpvs_.end())
                 ientry = genpvs_.insert(ientry,pv);
 
             // Add the vertex barcodes to the new or existent vertices
             ientry->genVertex.push_back((*ivertex)->barcode());
 
             // Collect final state descendants
             for (
                 HepMC::GenVertex::particle_iterator idecendants  = (*ivertex)->particles_begin(HepMC::descendants);
                 idecendants != (*ivertex)->particles_end(HepMC::descendants);
                 ++idecendants
             )
             {
                 if (isFinalstateParticle(*idecendants))
                     if ( find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) == ientry->finalstateParticles.end() )
                     {
                         ientry->finalstateParticles.push_back((*idecendants)->barcode());
                         HepMC::FourVector m = (*idecendants)->momentum();
 
                         ientry->ptot.setPx(ientry->ptot.px() + m.px());
                         ientry->ptot.setPy(ientry->ptot.py() + m.py());
                         ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
                         ientry->ptot.setE(ientry->ptot.e() + m.e());
                         ientry->ptsq += m.perp() * m.perp();
 
                         if ( m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants) ) ientry->nGenTrk++;
                     }
             }
             idx++;
         }
     }
 
     std::sort(genpvs_.begin(), genpvs_.end());
 }

VertexHistory const& VertexClassifier::history ( ) const

inline

Returns a reference to the vertex history used in the classification.

Definition at line 48 of file VertexClassifier.h.

References tracer_.

Referenced by VertexHistoryAnalyzer::analyze(), recoBSVTagInfoValidationAnalyzer::analyze(), and SVTagInfoValidationAnalyzer::analyze().

     {
         return tracer_;
     }

bool VertexClassifier::isCharged ( const HepMC::GenParticle * p )

private

Definition at line 444 of file VertexClassifier.cc.

References particleDataTable_.

Referenced by genPrimaryVertices().

 {
     const ParticleData * part = particleDataTable_->particle( p->pdg_id() );
     if (part)
         return part->charge()!=0;
     else
     {
         // the new/improved particle table doesn't know anti-particles
         return  particleDataTable_->particle( -p->pdg_id() ) != 0;
     }
 }

bool VertexClassifier::isFinalstateParticle ( const HepMC::GenParticle * p )

private

Definition at line 438 of file VertexClassifier.cc.

Referenced by genPrimaryVertices().

 {
     return !p->end_vertex() && p->status() == 1;
 }

void VertexClassifier::newEvent	(	edm::Event const &	event,
		edm::EventSetup const &	setup
	)

virtual

Pre-process event information (for accessing reconstraction information)

Reimplemented in VertexClassifierByProxy< Collection >, and VertexClassifierByProxy< reco::SecondaryVertexTagInfoCollection >.

Definition at line 35 of file VertexClassifier.cc.

References genPrimaryVertices(), edm::EventSetup::getData(), hepMCLabel_, mcInformation_, VertexHistory::newEvent(), particleDataTable_, and tracer_.

Referenced by VertexHistoryAnalyzer::analyze(), and VertexClassifierByProxy< reco::SecondaryVertexTagInfoCollection >::newEvent().

 {
     // Get the new event information for the tracer
     tracer_.newEvent(event, setup);
 
     // Get hepmc of the event
     event.getByLabel(hepMCLabel_, mcInformation_);
 
     // Get the partivle data table
     setup.getData(particleDataTable_);
 
     // Create the list of primary vertices associated to the event
     genPrimaryVertices();
 }

void VertexClassifier::processesAtGenerator ( )

private

Get all the information related to decay process.

Definition at line 123 of file VertexClassifier.cc.

References funct::abs(), VertexCategories::BWeakDecay, VertexCategories::CWeakDecay, VertexCategories::flags_, HistoryBase::genVertexTrail(), VertexCategories::JpsiDecay, VertexCategories::KsDecay, VertexCategories::LambdaDecay, VertexCategories::LongLivedDecay, longLivedDecayLength_, VertexCategories::OmegaDecay, parents, particleDataTable_, VertexCategories::SigmaMinusDecay, VertexCategories::SigmaPlusDecay, tracer_, update, and VertexCategories::XiDecay.

Referenced by evaluate().

 {
     // Get the generated vetices from track history
     VertexHistory::GenVertexTrail const & genVertexTrail = tracer_.genVertexTrail();
 
     // Loop over the generated vertices
     for (
         VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin();
         ivertex != genVertexTrail.end();
         ++ivertex
     )
     {
         // Get the pointer to the vertex by removing the const-ness (no const methos in HepMC::GenVertex)
         HepMC::GenVertex * vertex = const_cast<HepMC::GenVertex *>(*ivertex);
 
         // Loop over the sources looking for specific decays
         for (
             HepMC::GenVertex::particle_iterator iparent = vertex->particles_begin(HepMC::parents);
             iparent != vertex->particles_end(HepMC::parents);
             ++iparent
         )
         {
             // Collect the pdgid of the parent
             int pdgid = std::abs((*iparent)->pdg_id());
             // Get particle type
             HepPDT::ParticleID particleID(pdgid);
 
             // Check if the particle type is valid one
             if (particleID.isValid())
             {
                 // Get particle data
                 ParticleData const * particleData = particleDataTable_->particle(particleID);
                 // Check if the particle exist in the table
                 if (particleData)
                 {
                     // Check if their life time is bigger than longLivedDecayLength_
                     if ( particleData->lifetime() > longLivedDecayLength_ )
                     {
                         // Check for B, C weak decays and long lived decays
                         update(flags_[BWeakDecay], particleID.hasBottom());
                         update(flags_[CWeakDecay], particleID.hasCharm());
                         update(flags_[LongLivedDecay], true);
                     }
                     // Check Tau, Ks and Lambda decay
                     update(flags_[TauDecay], pdgid == 15);
                     update(flags_[KsDecay], pdgid == 310);
                     update(flags_[LambdaDecay], pdgid == 3122);
                     update(flags_[JpsiDecay], pdgid == 443);
                     update(flags_[XiDecay], pdgid == 3312);
                     update(flags_[OmegaDecay], pdgid == 3334);
                     update(flags_[SigmaPlusDecay], pdgid == 3222);
                     update(flags_[SigmaMinusDecay], pdgid == 3112);
                 }
             }
         }
     }
 }

void VertexClassifier::processesAtSimulation ( )

private

Get information about conversion and other interactions.

Definition at line 182 of file VertexClassifier.cc.

Referenced by evaluate().

 {
     VertexHistory::SimVertexTrail const & simVertexTrail = tracer_.simVertexTrail();
 
     for (
         VertexHistory::SimVertexTrail::const_iterator ivertex = simVertexTrail.begin();
         ivertex != simVertexTrail.end();
         ++ivertex
     )
     {
         // pdgid of the real source parent vertex
         int pdgid = 0;
 
         // select the original source in case of combined vertices
         bool flag = false;
         TrackingVertex::tp_iterator itd, its;
 
         for (its = (*ivertex)->sourceTracks_begin(); its != (*ivertex)->sourceTracks_end(); ++its)
         {
             for (itd = (*ivertex)->daughterTracks_begin(); itd != (*ivertex)->daughterTracks_end(); ++itd)
                 if (itd != its)
                 {
                     flag = true;
                     break;
                 }
             if (flag)
                 break;
         }
         // Collect the pdgid of the original source track
         if ( its != (*ivertex)->sourceTracks_end() )
             pdgid = std::abs((*its)->pdgId());
         else
             pdgid = 0;
 
     // Geant4 process type is selected using first Geant4 vertex assigned to 
         // the TrackingVertex
     unsigned int processG4 = 0;
     
         if((*ivertex)->nG4Vertices() > 0) {
       processG4 = (*(*ivertex)->g4Vertices_begin()).processType();
     }
     
     unsigned int process = g4toCMSProcMap_.processId(processG4);
     
     // Flagging all the different processes
     update(
            flags_[KnownProcess],
            process != CMS::Undefined &&
            process != CMS::Unknown &&
            process != CMS::Primary
            );
 
     update(flags_[UndefinedProcess], process == CMS::Undefined);
     update(flags_[UnknownProcess], process == CMS::Unknown);
     update(flags_[PrimaryProcess], process == CMS::Primary);
     update(flags_[HadronicProcess], process == CMS::Hadronic);
     update(flags_[DecayProcess], process == CMS::Decay);
     update(flags_[ComptonProcess], process == CMS::Compton);
     update(flags_[AnnihilationProcess], process == CMS::Annihilation);
     update(flags_[EIoniProcess], process == CMS::EIoni);
     update(flags_[HIoniProcess], process == CMS::HIoni);
     update(flags_[MuIoniProcess], process == CMS::MuIoni);
     update(flags_[PhotonProcess], process == CMS::Photon);
     update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
     update(flags_[ConversionsProcess], process == CMS::Conversions);
     update(flags_[EBremProcess], process == CMS::EBrem);
     update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
     update(flags_[MuBremProcess], process == CMS::MuBrem);
     update(flags_[MuNuclProcess], process == CMS::MuNucl);
 
 
         // Loop over the simulated particles
         for (
             TrackingVertex::tp_iterator iparticle = (*ivertex)->daughterTracks_begin();
             iparticle != (*ivertex)->daughterTracks_end();
             ++iparticle
         )
         {
 
             if ( (*iparticle)->numberOfTrackerLayers() )
             {
 
                 // Special treatment for decays
                 if (process == CMS::Decay)
                 {
                     // Get particle type
                     HepPDT::ParticleID particleID(pdgid);
                     // Check if the particle type is valid one
                     if (particleID.isValid())
                     {
                         // Get particle data
                         ParticleData const * particleData = particleDataTable_->particle(particleID);
                         // Check if the particle exist in the table
                         if (particleData)
                         {
                             // Check if their life time is bigger than 1e-14
                             if ( particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_ )
                             {
                                 // Check for B, C weak decays and long lived decays
                                 update(flags_[BWeakDecay], particleID.hasBottom());
                                 update(flags_[CWeakDecay], particleID.hasCharm());
                                 update(flags_[LongLivedDecay], true);
                             }
                             // Check Tau, Ks and Lambda decay
                             update(flags_[TauDecay], pdgid == 15);
                             update(flags_[KsDecay], pdgid == 310);
                             update(flags_[LambdaDecay], pdgid == 3122);
                             update(flags_[JpsiDecay], pdgid == 443);
                             update(flags_[XiDecay], pdgid == 3312);
                             update(flags_[OmegaDecay], pdgid == 3334);
                             update(flags_[SigmaPlusDecay], pdgid == 3222);
                             update(flags_[SigmaMinusDecay], pdgid == 3112);
                         }
                     }
                 }
             }
         }
     }
 }

void VertexClassifier::reconstructionInformation ( reco::TrackBaseRef const & )

private

Get reconstruction information.

void VertexClassifier::simulationInformation ( )

private

Get all the information related to the simulation details.

Definition at line 114 of file VertexClassifier.cc.

References EncodedEventId::bunchCrossing(), EncodedEventId::event(), VertexCategories::flags_, VertexCategories::SignalEvent, HistoryBase::simVertex(), and tracer_.

Referenced by evaluate().

 {
     // Get the event id for the initial TP.
     EncodedEventId eventId = tracer_.simVertex()->eventId();
     // Check for signal events
     flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
 }

void VertexClassifier::vertexInformation ( )

private

Get geometrical information about the vertices.

Definition at line 303 of file VertexClassifier.cc.

References HLT_25ns14e33_v1_cff::clusters, HLT_25ns14e33_v1_cff::distance, VertexCategories::flags_, genpvs_, HistoryBase::genVertexTrail(), AlCaHLTBitMon_ParallelJobs::p, funct::pow(), VertexCategories::PrimaryVertex, VertexCategories::SecondaryVertex, HistoryBase::simVertexTrail(), mathSSE::sqrt(), VertexCategories::TertiaryVertex, tracer_, vertexClusteringDistance_, VertexClassifier::GeneratedPrimaryVertex::x, VertexClassifier::GeneratedPrimaryVertex::y, and VertexClassifier::GeneratedPrimaryVertex::z.

Referenced by evaluate().

 {
     // Helper class for clusterization
     typedef std::multimap<double, HepMC::ThreeVector> Clusters;
     typedef std::pair<double, HepMC::ThreeVector> ClusterPair;
 
     Clusters clusters;
 
     // Get the main primary vertex from the list
     GeneratedPrimaryVertex const & genpv = genpvs_.back();
 
     // Get the generated history of the tracks
     const VertexHistory::GenVertexTrail & genVertexTrail = tracer_.genVertexTrail();
 
     // Unit transformation from mm to cm
     double const mm = 0.1;
 
     // Loop over the generated vertexes
     for (
         VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin();
         ivertex != genVertexTrail.end();
         ++ivertex
     )
     {
         // Check vertex exist
         if (*ivertex)
         {
             // Measure the distance2 respecto the primary vertex
             HepMC::ThreeVector p = (*ivertex)->point3d();
             double distance = sqrt( pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2) );
 
             // If there is not any clusters add the first vertex.
             if ( clusters.empty() )
             {
                 clusters.insert( ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)) );
                 continue;
             }
 
             // Check if there is already a cluster in the given distance from primary vertex
             Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
 
             if ( icluster == clusters.upper_bound(distance + vertexClusteringDistance_) )
             {
                 clusters.insert ( ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)) );
                 continue;
             }
 
             bool cluster = false;
 
             // Looping over the vertex clusters of a given distance from primary vertex
             for (;
                     icluster != clusters.upper_bound(distance + vertexClusteringDistance_);
                     ++icluster
                 )
             {
                 double difference = sqrt (
                                         pow(p.x() * mm - icluster->second.x(), 2) +
                                         pow(p.y() * mm - icluster->second.y(), 2) +
                                         pow(p.z() * mm - icluster->second.z(), 2)
                                     );
 
                 if ( difference < vertexClusteringDistance_ )
                 {
                     cluster = true;
                     break;
                 }
             }
 
             if (!cluster) clusters.insert ( ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)) );
         }
     }
 
     const VertexHistory::SimVertexTrail & simVertexTrail = tracer_.simVertexTrail();
 
     // Loop over the generated particles
     for (
         VertexHistory::SimVertexTrail::const_reverse_iterator ivertex = simVertexTrail.rbegin();
         ivertex != simVertexTrail.rend();
         ++ivertex
     )
     {
         // Look for those with production vertex
         TrackingVertex::LorentzVector p = (*ivertex)->position();
 
         double distance = sqrt( pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2) );
 
         // If there is not any clusters add the first vertex.
         if ( clusters.empty() )
         {
             clusters.insert( ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())) );
             continue;
         }
 
         // Check if there is already a cluster in the given distance from primary vertex
         Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
 
         if ( icluster == clusters.upper_bound(distance + vertexClusteringDistance_) )
         {
             clusters.insert ( ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())) );
             continue;
         }
 
         bool cluster = false;
 
         // Looping over the vertex clusters of a given distance from primary vertex
         for (;
                 icluster != clusters.upper_bound(distance + vertexClusteringDistance_);
                 ++icluster
             )
         {
             double difference = sqrt (
                                     pow(p.x() - icluster->second.x(), 2) +
                                     pow(p.y() - icluster->second.y(), 2) +
                                     pow(p.z() - icluster->second.z(), 2)
                                 );
 
             if ( difference < vertexClusteringDistance_ )
             {
                 cluster = true;
                 break;
             }
         }
 
         if (!cluster) clusters.insert ( ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())) );
     }
 
     if ( clusters.size() == 1 )
         flags_[PrimaryVertex] = true;
     else if ( clusters.size() == 2 )
         flags_[SecondaryVertex] = true;
     else
         flags_[TertiaryVertex] = true;
 }