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

#include <VertexClassifier.h>

Inheritance diagram for VertexClassifier:

Classes
struct	G4

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)
	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
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 *	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 16 of file VertexClassifier.h.

Member Typedef Documentation

typedef VertexCategories VertexClassifier::Categories

Type to the associate category.

Definition at line 22 of file VertexClassifier.h.

Constructor & Destructor Documentation

VertexClassifier::VertexClassifier ( edm::ParameterSet const & pset )

Constructor by ParameterSet.

Definition at line 17 of file VertexClassifier.cc.

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

                                                                  : VertexCategories(), tracer_(config),
         hepMCLabel_( config.getUntrackedParameter<edm::InputTag>("hepMC") )
 {
     // 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 27 of file VertexClassifier.h.

27 {}

Member Function Documentation

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

Classify the RecoVertex in categories.

Definition at line 47 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 77 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 39 of file VertexClassifier.h.

References evaluate().

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

void VertexClassifier::genPrimaryVertices ( )

private

Definition at line 445 of file VertexClassifier.cc.

References abs, event(), spr::find(), genpvs_, isCharged(), isFinalstateParticle(), m, mcInformation_, parents, pos, funct::pow(), 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 45 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 432 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 426 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 31 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 119 of file VertexClassifier.cc.

References abs, VertexCategories::BWeakDecay, VertexCategories::CWeakDecay, VertexCategories::flags_, HistoryBase::genVertexTrail(), VertexCategories::JpsiDecay, VertexCategories::KsDecay, VertexCategories::LambdaDecay, VertexCategories::LongLivedDecay, longLivedDecayLength_, VertexCategories::OmegaDecay, parents, particleDataTable_, RecoTau_DiTaus_pt_20-420_cfg::ParticleID, 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 178 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;
 
         // Loop over the simulated particles
         for (
             TrackingVertex::tp_iterator iparticle = (*ivertex)->daughterTracks_begin();
             iparticle != (*ivertex)->daughterTracks_end();
             ++iparticle
         )
         {
 
             if ( (*iparticle)->matchedHit() )
             {
                 // Collect the G4 process of the first psimhit (it should be the same for all of them)
                 unsigned short process = (*iparticle)->pSimHit_begin()->processType();
 
                 // Flagging all the different processes
 
                 update(
                     flags_[KnownProcess],
                     process != G4::Undefined &&
                     process != G4::Unknown &&
                     process != G4::Primary
                 );
 
                 update(flags_[UndefinedProcess], process == G4::Undefined);
                 update(flags_[UnknownProcess], process == G4::Unknown);
                 update(flags_[PrimaryProcess], process == G4::Primary);
                 update(flags_[HadronicProcess], process == G4::Hadronic);
                 update(flags_[DecayProcess], process == G4::Decay);
                 update(flags_[ComptonProcess], process == G4::Compton);
                 update(flags_[AnnihilationProcess], process == G4::Annihilation);
                 update(flags_[EIoniProcess], process == G4::EIoni);
                 update(flags_[HIoniProcess], process == G4::HIoni);
                 update(flags_[MuIoniProcess], process == G4::MuIoni);
                 update(flags_[PhotonProcess], process == G4::Photon);
                 update(flags_[MuPairProdProcess], process == G4::MuPairProd);
                 update(flags_[ConversionsProcess], process == G4::Conversions);
                 update(flags_[EBremProcess], process == G4::EBrem);
                 update(flags_[SynchrotronRadiationProcess], process == G4::SynchrotronRadiation);
                 update(flags_[MuBremProcess], process == G4::MuBrem);
                 update(flags_[MuNuclProcess], process == G4::MuNucl);
 
                 // Special treatment for decays
                 if (process == G4::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 110 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 291 of file VertexClassifier.cc.

References VertexCategories::flags_, genpvs_, HistoryBase::genVertexTrail(), L1TEmulatorMonitor_cff::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
     VertexHistory::GenVertexTrail & genVertexTrail = const_cast<VertexHistory::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)) );
         }
     }
 
     VertexHistory::SimVertexTrail & simVertexTrail = const_cast<VertexHistory::SimVertexTrail &> (tracer_.simVertexTrail());
 
     // Loop over the generated particles
     for (
         VertexHistory::SimVertexTrail::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;
 }