Index: modules/BeamSpotFilter.cc
===================================================================
--- modules/BeamSpotFilter.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/BeamSpotFilter.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,85 @@
+/** \class BeamSpotFilter
+ *
+ *  Extracts beam spot 
+ *
+ *  \author Michele Selvaggi
+ *
+ */
+
+
+#include "modules/BeamSpotFilter.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+BeamSpotFilter::BeamSpotFilter() :
+  fItInputArray(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+BeamSpotFilter::~BeamSpotFilter()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void BeamSpotFilter::Init()
+{
+ 
+  // import input array
+  fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array
+
+  fOutputArray = ExportArray(GetString("OutputArray", "filteredParticles"));
+}
+
+//------------------------------------------------------------------------------
+
+void BeamSpotFilter::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void BeamSpotFilter::Process()
+{
+  Candidate *candidate;
+  Bool_t passed = false;
+  
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())) && !passed)
+  {
+    if(candidate->IsPU == 0) passed = true;
+    fOutputArray->Add(candidate);
+  }
+  
+}
+
+ 
Index: modules/BeamSpotFilter.h
===================================================================
--- modules/BeamSpotFilter.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/BeamSpotFilter.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,43 @@
+//------------------------------------------------------------------------------
+
+#ifndef BeamSpotFilter_h
+#define BeamSpotFilter_h
+
+/** \class BeamSpotFilter
+ *
+ *  Extracts beam spot 
+ *
+ *  \author Michele Selvaggi
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+
+class BeamSpotFilter: public DelphesModule
+{
+public:
+
+  BeamSpotFilter();
+  ~BeamSpotFilter();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  Float_t fPassedOne;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(BeamSpotFilter, 1)
+};
+
+#endif
Index: modules/ImpactParameterSmearing.cc
===================================================================
--- modules/ImpactParameterSmearing.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/ImpactParameterSmearing.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -96,5 +96,5 @@
 {
   Candidate *candidate, *particle, *mother;
-  Double_t xd, yd, zd, dxy, sx, sy, sz, ddxy;
+  Double_t xd, yd, zd, d0, sx, sy, sz, dd0;
   Double_t pt, eta, px, py, phi, e;
 
@@ -103,5 +103,5 @@
   {
 
-    // take momentum before smearing (otherwise apply double smearing on dxy)
+    // take momentum before smearing (otherwise apply double smearing on d0)
     particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
 
@@ -131,7 +131,7 @@
 
     // calculate impact parameter (after-smearing)
-    dxy = (xd*py - yd*px)/pt;
+    d0 = (xd*py - yd*px)/pt;
 
-    ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
+    dd0 = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
 
     // fill smeared values in candidate
@@ -143,6 +143,6 @@
     candidate->Zd = zd;
 
-    candidate->Dxy = dxy;
-    candidate->SDxy = ddxy;
+    candidate->D0 = d0;
+    candidate->ErrorD0 = dd0;
 
     candidate->AddCandidate(mother);
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/ModulesLinkDef.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -35,4 +35,5 @@
 #include "modules/EnergySmearing.h"
 #include "modules/MomentumSmearing.h"
+#include "modules/TrackSmearing.h"
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
@@ -57,4 +58,5 @@
 #include "modules/StatusPidFilter.h"
 #include "modules/PdgCodeFilter.h"
+#include "modules/BeamSpotFilter.h"
 #include "modules/Cloner.h"
 #include "modules/Weighter.h"
@@ -62,4 +64,7 @@
 #include "modules/JetFlavorAssociation.h"
 #include "modules/JetFakeParticle.h"
+#include "modules/VertexSorter.h"
+#include "modules/VertexFinder.h"
+#include "modules/VertexFinder4D.h"
 #include "modules/ExampleModule.h"
 
@@ -79,4 +84,5 @@
 #pragma link C++ class EnergySmearing+;
 #pragma link C++ class MomentumSmearing+;
+#pragma link C++ class TrackSmearing+;
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
@@ -101,4 +107,5 @@
 #pragma link C++ class StatusPidFilter+;
 #pragma link C++ class PdgCodeFilter+;
+#pragma link C++ class BeamSpotFilter+;
 #pragma link C++ class Cloner+;
 #pragma link C++ class Weighter+;
@@ -106,4 +113,7 @@
 #pragma link C++ class JetFlavorAssociation+;
 #pragma link C++ class JetFakeParticle+;
+#pragma link C++ class VertexSorter+;
+#pragma link C++ class VertexFinder+;
+#pragma link C++ class VertexFinder4D+;
 #pragma link C++ class ExampleModule+;
 
Index: modules/ParticlePropagator.cc
===================================================================
--- modules/ParticlePropagator.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/ParticlePropagator.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -67,4 +67,5 @@
 }
 
+
 //------------------------------------------------------------------------------
 
@@ -91,4 +92,13 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // import beamspot
+  try
+  {
+    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+  }
+  catch(runtime_error &e)
+  {
+    fBeamSpotInputArray = 0;
+  }  
   // create output arrays
 
@@ -111,5 +121,5 @@
 {
   Candidate *candidate, *mother;
-  TLorentzVector candidatePosition, candidateMomentum;
+  TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
   Double_t px, py, pz, pt, pt2, e, q;
   Double_t x, y, z, t, r, phi;
@@ -120,8 +130,18 @@
   Double_t tmp, discr, discr2;
   Double_t delta, gammam, omega, asinrho;
-  Double_t rcu, rc2, dxy, xd, yd, zd;
-
+  Double_t rcu, rc2, xd, yd, zd;
+  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
+  Double_t bsx, bsy, bsz;
+  	 
   const Double_t c_light = 2.99792458E8;
-
+ 
+  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0) 
+    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  else
+  {
+    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
+    beamSpotPosition = beamSpotCandidate.Position;
+  }
+ 
   fItInputArray->Reset();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
@@ -132,4 +152,9 @@
     y = candidatePosition.Y()*1.0E-3;
     z = candidatePosition.Z()*1.0E-3;
+    
+    bsx = beamSpotPosition.X()*1.0E-3;
+    bsy = beamSpotPosition.Y()*1.0E-3;
+    bsz = beamSpotPosition.Z()*1.0E-3;
+   
     q = candidate->Charge;
 
@@ -181,10 +206,14 @@
       y_t = y + py*t;
       z_t = z + pz*t;
+      
+      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
 
       mother = candidate;
       candidate = static_cast<Candidate*>(candidate->Clone());
 
+      candidate->InitialPosition = candidatePosition;
       candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
-
+      candidate->L = l*1.0E3;
+   
       candidate->Momentum = candidateMomentum;
       candidate->AddCandidate(mother);
@@ -238,8 +267,22 @@
       yd = (rc2 > 0.0) ? yd / rc2 : -999;
       zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
-
-      // calculate impact paramater
-      dxy = (xd*py - yd*px)/pt;
-
+      
+      // use perigee momentum rather than original particle
+      // momentum, since the orignal particle momentum isn't known 
+      
+      px = TMath::Sign(1.0,r) * pt * (-y_c / r_c);
+      py = TMath::Sign(1.0,r) * pt * (x_c / r_c);
+      etap = candidateMomentum.Eta();
+      phip = TMath::ATan2(py, px);
+       
+      candidateMomentum.SetPtEtaPhiE(pt, etap, phip, candidateMomentum.E());
+      
+      // calculate additional track parameters (correct for beamspot position)
+       
+      d0        = (  (x - bsx) * py - (y - bsy) * px) / pt;
+      dz        = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
+      p         = candidateMomentum.P();
+      ctgTheta  = 1.0 / TMath::Tan (candidateMomentum.Theta ());
+          
       // 3. time evaluation t = TMath::Min(t_r, t_z)
       //    t_r : time to exit from the sides
@@ -287,14 +330,33 @@
       r_t = TMath::Hypot(x_t, y_t);
 
+      
+      // compute path length for an helix
+      
+      alpha = pz*1.0E9 / c_light / gammam;
+      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
+          
       if(r_t > 0.0)
       {
+        
+        // store these variables before cloning
+           
+        candidate->D0 = d0*1.0E3;
+        candidate->DZ = dz*1.0E3;
+        candidate->P  = p;
+        candidate->PT = pt;
+        candidate->CtgTheta = ctgTheta;
+        candidate->Phi = phip;
+      
         mother = candidate;
         candidate = static_cast<Candidate*>(candidate->Clone());
 
+        candidate->InitialPosition = candidatePosition;
         candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
 
         candidate->Momentum = candidateMomentum;
-        candidate->Dxy = dxy*1.0E3;
-        candidate->Xd = xd*1.0E3;
+        
+	    candidate->L  =  l*1.0E3;
+        
+	    candidate->Xd = xd*1.0E3;
         candidate->Yd = yd*1.0E3;
         candidate->Zd = zd*1.0E3;
Index: modules/ParticlePropagator.h
===================================================================
--- modules/ParticlePropagator.h	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/ParticlePropagator.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -35,4 +35,5 @@
 class TClonesArray;
 class TIterator;
+class TLorentzVector;
 
 class ParticlePropagator: public DelphesModule
@@ -55,4 +56,5 @@
 
   const TObjArray *fInputArray; //!
+  const TObjArray *fBeamSpotInputArray; //!
 
   TObjArray *fOutputArray; //!
Index: modules/PileUpMerger.cc
===================================================================
--- modules/PileUpMerger.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/PileUpMerger.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -115,8 +115,8 @@
   TDatabasePDG *pdg = TDatabasePDG::Instance();
   TParticlePDG *pdgParticle;
-  Int_t pid;
+  Int_t pid, nch, nvtx = -1;
   Float_t x, y, z, t, vx, vy;
-  Float_t px, py, pz, e;
-  Double_t dz, dphi, dt;
+  Float_t px, py, pz, e, pt;
+  Double_t dz, dphi, dt, sumpt2, dz0, dt0;
   Int_t numberOfEvents, event, numberOfParticles;
   Long64_t allEntries, entry;
@@ -132,9 +132,19 @@
   fFunction->GetRandom2(dz, dt);
 
+  dz0 = -1.0e6;
+  dt0 = -1.0e6;
+  
   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
   dz *= 1.0E3; // necessary in order to make z in mm
+  
+  //cout<<dz<<","<<dt<<endl;
+ 
   vx = 0.0;
   vy = 0.0;
+  
   numberOfParticles = fInputArray->GetEntriesFast();
+  nch = 0;
+  sumpt2 = 0.0; 
+   
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
@@ -143,19 +153,41 @@
     z = candidate->Position.Z();
     t = candidate->Position.T();
-    candidate->Position.SetZ(z + dz);
-    candidate->Position.SetT(t + dt);
+    pt = candidate->Momentum.Pt();
+       
+    // take postion and time from first stable particle
+    if (dz0 < -999999.0)
+      dz0 = z;
+    if (dt0 < -999999.0)
+      dt0 = t;
+
+    // cancel any possible offset in position and time the input file
+    candidate->Position.SetZ(z - dz0 + dz);
+    candidate->Position.SetT(t - dt0 + dt);
+    
     fParticleOutputArray->Add(candidate);
-  }
-
+ 
+    if(TMath::Abs(candidate->Charge) >  1.0E-9)
+    {
+      nch++;    
+      sumpt2 += pt*pt;
+    }  
+  }
+ 
   if(numberOfParticles > 0)
   {
-    vx /= numberOfParticles;
-    vy /= numberOfParticles;
-  }
-
+    vx /= sumpt2;
+    vy /= sumpt2;
+  }
+  
+  nvtx++;
   factory = GetFactory();
 
   vertex = factory->NewCandidate();
   vertex->Position.SetXYZT(vx, vy, dz, dt);
+  vertex->ClusterIndex = nvtx;
+  vertex->ClusterNDF = nch;
+  vertex->SumPT2 = sumpt2;
+  vertex->GenSumPT2 = sumpt2; 
+  
   fVertexOutputArray->Add(vertex);
 
@@ -170,4 +202,7 @@
       numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
       break;
+    case 2:
+      numberOfEvents = fMeanPileUp;
+      break;   
     default:
       numberOfEvents = gRandom->Poisson(fMeanPileUp);
@@ -198,5 +233,8 @@
     vx = 0.0;
     vy = 0.0;
+    
     numberOfParticles = 0;
+    sumpt2 = 0.0;
+    
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     {
@@ -224,6 +262,12 @@
       vx += candidate->Position.X();
       vy += candidate->Position.Y();
+      
       ++numberOfParticles;
-
+      if(TMath::Abs(candidate->Charge) >  1.0E-9)
+      {
+        nch++;    
+        sumpt2 += pt*pt;
+      } 
+       
       fParticleOutputArray->Add(candidate);
     }
@@ -234,12 +278,21 @@
       vy /= numberOfParticles;
     }
+    
+    nvtx++;
 
     vertex = factory->NewCandidate();
     vertex->Position.SetXYZT(vx, vy, dz, dt);
+    
+    vertex->ClusterIndex = nvtx; 
+    vertex->ClusterNDF = nch;
+    vertex->SumPT2 = sumpt2;
+    vertex->GenSumPT2 = sumpt2;
+    
     vertex->IsPU = 1;
 
     fVertexOutputArray->Add(vertex);
-  }
-}
-
-//------------------------------------------------------------------------------
+    
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/TimeSmearing.cc
===================================================================
--- modules/TimeSmearing.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/TimeSmearing.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -99,13 +99,14 @@
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
+    const TLorentzVector &candidatePosition = candidate->InitialPosition;
     t = candidatePosition.T()*1.0E-3/c_light;
 
     // apply smearing formula
     t = gRandom->Gaus(t, fTimeResolution);
-
+    
     mother = candidate;
     candidate = static_cast<Candidate*>(candidate->Clone());
-    candidate->Position.SetT(t*1.0E3*c_light);
+    candidate->InitialPosition.SetT(t*1.0E3*c_light);
+    candidate->ErrorT = fTimeResolution*1.0E3*c_light;
 
     candidate->AddCandidate(mother);
Index: modules/TrackCountingBTagging.cc
===================================================================
--- modules/TrackCountingBTagging.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/TrackCountingBTagging.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -96,6 +96,6 @@
 
   Double_t jpx, jpy;
-  Double_t dr, tpx, tpy, tpt;
-  Double_t xd, yd, dxy, ddxy, ip, sip;
+  Double_t dr, tpt;
+  Double_t xd, yd, d0, dd0, ip, sip;
 
   Int_t sign;
@@ -117,24 +117,20 @@
     {
       const TLorentzVector &trkMomentum = track->Momentum;
-
+      
       dr = jetMomentum.DeltaR(trkMomentum);
-
       tpt = trkMomentum.Pt();
-      tpx = trkMomentum.Px();
-      tpy = trkMomentum.Py();
-
       xd = track->Xd;
       yd = track->Yd;
-      dxy = TMath::Hypot(xd, yd);
-      ddxy = track->SDxy;
+      d0 = TMath::Hypot(xd, yd);
+      dd0 = track->ErrorD0;
 
       if(tpt < fPtMin) continue;
       if(dr > fDeltaR) continue;
-      if(dxy > fIPmax) continue;
+      if(d0 > fIPmax) continue;
 
       sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
 
-      ip = sign*dxy;
-      sip = ip / TMath::Abs(ddxy);
+      ip = sign*d0;
+      sip = ip / TMath::Abs(dd0);
 
       if(sip > fSigMin) count++;
Index: modules/TrackSmearing.cc
===================================================================
--- modules/TrackSmearing.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/TrackSmearing.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,380 @@
+/** \class TrackSmearing
+ *
+ *  Performs d0, dZ, p, Theta, Phi smearing of tracks.
+ *
+ *  \authors A. Hart, M. Selvaggi
+ *
+*/
+
+#include "modules/TrackSmearing.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TFile.h"
+#include "TProfile2D.h"
+
+#include <algorithm> 
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+TrackSmearing::TrackSmearing() :
+  fD0Formula(0), fDZFormula(0), fPFormula(0), fCtgThetaFormula(0), fPhiFormula(0), fItInputArray(0)
+{
+  fD0Formula = new DelphesFormula;
+  fDZFormula = new DelphesFormula;
+  fPFormula = new DelphesFormula;
+  fCtgThetaFormula = new DelphesFormula;
+  fPhiFormula = new DelphesFormula;
+}
+
+//------------------------------------------------------------------------------
+
+TrackSmearing::~TrackSmearing()
+{
+  if(fD0Formula) delete fD0Formula;
+  if(fDZFormula) delete fDZFormula;
+  if(fPFormula) delete fPFormula;
+  if(fCtgThetaFormula) delete fCtgThetaFormula;
+  if(fPhiFormula) delete fPhiFormula;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackSmearing::Init()
+{
+  // read resolution formula
+
+  // !!! IF WE WANT TO KEEP ROOT INPUT !!!
+  if (string (GetString("D0ResolutionFormula", "0.0")) != "0.0")
+    {
+      fD0Formula->Compile(GetString("D0ResolutionFormula", "0.0"));
+      fUseD0Formula = true;
+    }
+  else
+    {
+      fD0ResolutionFile = GetString("D0ResolutionFile", "errors.root");
+      fD0ResolutionHist = GetString("D0ResolutionHist", "d0");
+      fUseD0Formula = false;
+    }
+  if (string (GetString("DZResolutionFormula", "0.0")) != "0.0")
+    {
+      fDZFormula->Compile(GetString("DZResolutionFormula", "0.0"));
+      fUseDZFormula = true;
+    }
+  else
+    {
+      fDZResolutionFile = GetString("DZResolutionFile", "errors.root");
+      fDZResolutionHist = GetString("DZResolutionHist", "dz");
+      fUseDZFormula = false;
+    }
+  if (string (GetString("PResolutionFormula", "0.0")) != "0.0")
+    {
+      fPFormula->Compile(GetString("PResolutionFormula", "0.0"));
+      fUsePFormula = true;
+    }
+  else
+    {
+      fPResolutionFile = GetString("PResolutionFile", "errors.root");
+      fPResolutionHist = GetString("PResolutionHist", "p");
+      fUsePFormula = false;
+    }
+  if (string (GetString("CtgThetaResolutionFormula", "0.0")) != "0.0")
+    {
+      fCtgThetaFormula->Compile(GetString("CtgThetaResolutionFormula", "0.0"));
+      fUseCtgThetaFormula = true;
+    }
+  else
+    {
+      fCtgThetaResolutionFile = GetString("CtgThetaResolutionFile", "errors.root");
+      fCtgThetaResolutionHist = GetString("CtgThetaResolutionHist", "ctgTheta");
+      fUseCtgThetaFormula = false;
+    }
+  if (string (GetString("PhiResolutionFormula", "0.0")) != "0.0")
+    {
+      fPhiFormula->Compile(GetString("PhiResolutionFormula", "0.0"));
+      fUsePhiFormula = true;
+    }
+  else
+    {
+      fPhiResolutionFile = GetString("PhiResolutionFile", "errors.root");
+      fPhiResolutionHist = GetString("PhiResolutionHist", "phi");
+      fUsePhiFormula = false;
+    }
+
+  fApplyToPileUp = GetBool("ApplyToPileUp", true);
+
+  // import input array
+
+  fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles"));
+  fItInputArray = fInputArray->MakeIterator();
+ 
+  // import beamspot
+  try
+  {
+    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+  }
+  catch(runtime_error &e)
+  {
+    fBeamSpotInputArray = 0;
+  }  
+ 
+  // create output array
+
+  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
+}
+
+//------------------------------------------------------------------------------
+
+void TrackSmearing::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackSmearing::Process()
+{
+  Int_t iCandidate = 0;
+  TLorentzVector beamSpotPosition;
+  Candidate *candidate, *mother;
+  Double_t pt, eta, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;
+  Double_t x, y, z, t, px, py, pz, theta;
+  TProfile2D *d0ErrorHist = NULL,
+             *dzErrorHist = NULL,
+             *pErrorHist = NULL,
+             *ctgThetaErrorHist = NULL,
+             *phiErrorHist = NULL;
+
+  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0) 
+    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  else
+  {
+    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At (0));
+    beamSpotPosition = beamSpotCandidate.Position;
+  }
+
+  if (!fUseD0Formula)
+  {
+     TFile *fin = TFile::Open (fD0ResolutionFile.c_str ());
+     d0ErrorHist = (TProfile2D *) fin->Get (fD0ResolutionHist.c_str ());
+     d0ErrorHist->SetDirectory (0);
+     fin->Close ();
+  }
+  if (!fUseDZFormula)
+  {
+     TFile *fin = TFile::Open (fDZResolutionFile.c_str ());
+     dzErrorHist = (TProfile2D *) fin->Get (fDZResolutionHist.c_str ());
+     dzErrorHist->SetDirectory (0);
+     fin->Close ();
+  }
+  if (!fUsePFormula)
+  {
+     TFile *fin = TFile::Open (fPResolutionFile.c_str ());
+     pErrorHist = (TProfile2D *) fin->Get (fPResolutionHist.c_str ());
+     pErrorHist->SetDirectory (0);
+     fin->Close ();
+  }
+  if (!fUseCtgThetaFormula)
+  {
+     TFile *fin = TFile::Open (fCtgThetaResolutionFile.c_str ());
+     ctgThetaErrorHist = (TProfile2D *) fin->Get (fCtgThetaResolutionHist.c_str ());
+     ctgThetaErrorHist->SetDirectory (0);
+     fin->Close ();
+  }
+  if (!fUsePhiFormula)
+  {
+     TFile *fin = TFile::Open (fPhiResolutionFile.c_str ());
+     phiErrorHist = (TProfile2D *) fin->Get (fPhiResolutionHist.c_str ());
+     phiErrorHist->SetDirectory (0);
+     fin->Close ();
+  }
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+  {
+   
+    const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->InitialPosition;
+   
+    pt = momentum.Pt();
+    eta = momentum.Eta();
+
+    d0 = trueD0 = candidate->D0;
+    dz = trueDZ = candidate->DZ;
+     
+    p = trueP = candidate->P;
+    ctgTheta = trueCtgTheta = candidate->CtgTheta;
+    phi = truePhi = candidate->Phi;
+
+    if (fUseD0Formula)
+      d0Error = fD0Formula->Eval(pt, eta);
+    else
+    {
+       Int_t xbin, ybin;
+
+       xbin = pt < d0ErrorHist->GetXaxis ()->GetXmax () ? d0ErrorHist->GetXaxis ()->FindBin (pt)
+            : d0ErrorHist->GetXaxis ()->GetBinCenter (d0ErrorHist->GetXaxis ()->GetNbins ());
+       ybin = d0ErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
+       d0Error = d0ErrorHist->GetBinContent (xbin, ybin);
+       if (!d0Error)
+          d0Error = -1.0;
+    }
+    if (d0Error < 0.0)
+      continue;
+
+    if (fUseDZFormula)
+      dzError = fDZFormula->Eval(pt, eta);
+    else
+    {
+       Int_t xbin, ybin;
+
+       xbin = pt < dzErrorHist->GetXaxis ()->GetXmax () ? dzErrorHist->GetXaxis ()->FindBin (pt)
+            : dzErrorHist->GetXaxis ()->GetBinCenter (dzErrorHist->GetXaxis ()->GetNbins ());
+       ybin = dzErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
+       dzError = dzErrorHist->GetBinContent (xbin, ybin);
+       if (!dzError)
+          dzError = -1.0;
+    }
+    if (dzError < 0.0)
+      continue;
+
+    if (fUsePFormula)
+      pError = fPFormula->Eval(pt, eta) * p;
+    else
+    {
+       Int_t xbin, ybin;
+
+       xbin = pt < pErrorHist->GetXaxis ()->GetXmax () ? pErrorHist->GetXaxis ()->FindBin (pt)
+            : pErrorHist->GetXaxis ()->GetBinCenter (pErrorHist->GetXaxis ()->GetNbins ());
+       ybin = pErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
+       pError = pErrorHist->GetBinContent (xbin, ybin) * p;
+       if (!pError)
+          pError = -1.0;
+    }
+    if (pError < 0.0)
+      continue;
+
+    if (fUseCtgThetaFormula)
+      ctgThetaError = fCtgThetaFormula->Eval(pt, eta);
+    else
+    {
+       Int_t xbin, ybin;
+
+       xbin = pt < ctgThetaErrorHist->GetXaxis ()->GetXmax () ? ctgThetaErrorHist->GetXaxis ()->FindBin (pt)
+            : ctgThetaErrorHist->GetXaxis ()->GetBinCenter (ctgThetaErrorHist->GetXaxis ()->GetNbins ());
+       ybin = ctgThetaErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
+       ctgThetaError = ctgThetaErrorHist->GetBinContent (xbin, ybin);
+       if (!ctgThetaError)
+         ctgThetaError = -1.0;
+    }
+    if (ctgThetaError < 0.0)
+      continue;
+
+    if (fUsePhiFormula)
+      phiError = fPhiFormula->Eval(pt, eta);
+    else
+    {
+       Int_t xbin, ybin;
+
+       xbin = pt < phiErrorHist->GetXaxis ()->GetXmax () ? phiErrorHist->GetXaxis ()->FindBin (pt)
+            : phiErrorHist->GetXaxis ()->GetBinCenter (phiErrorHist->GetXaxis ()->GetNbins ());
+       ybin = phiErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta));
+       phiError = phiErrorHist->GetBinContent (xbin, ybin);
+       if (!phiError)
+          phiError = -1.0;
+    }
+    if (phiError < 0.0)
+      continue;
+
+    if (fApplyToPileUp || !candidate->IsPU)
+    {
+       d0 = gRandom->Gaus(d0, d0Error);
+       dz = gRandom->Gaus(dz, dzError);
+       p = gRandom->Gaus(p, pError);
+       ctgTheta = gRandom->Gaus(ctgTheta, ctgThetaError);
+       phi = gRandom->Gaus(phi, phiError);
+    }
+
+    if(p < 0.0) continue;
+    while (phi > TMath::Pi ()) phi -= TMath::TwoPi ();
+    while (phi <= -TMath::Pi ()) phi += TMath::TwoPi ();
+
+    mother = candidate;
+    candidate = static_cast<Candidate*>(candidate->Clone());
+    candidate->D0 = d0;
+    candidate->DZ = dz;
+    candidate->P = p;
+    candidate->CtgTheta = ctgTheta;
+    candidate->Phi = phi;
+
+    theta = TMath::ACos(ctgTheta / TMath::Sqrt (1.0 + ctgTheta * ctgTheta));
+    candidate->Momentum.SetPx (p * TMath::Cos (phi) * TMath::Sin (theta));
+    candidate->Momentum.SetPy (p * TMath::Sin (phi) * TMath::Sin (theta));
+    candidate->Momentum.SetPz (p * TMath::Cos (theta));
+    candidate->Momentum.SetE (candidate->Momentum.Pt () * TMath::CosH (eta));
+    candidate->PT = candidate->Momentum.Pt ();
+
+    x = position.X ();
+    y = position.Y ();
+    z = position.Z ();
+    t = position.T ();
+    px = candidate->Momentum.Px ();
+    py = candidate->Momentum.Py ();
+    pz = candidate->Momentum.Pz ();
+    pt = candidate->Momentum.Pt ();
+    
+    // -- solve for delta: d0' = ( (x+delta)*py' - (y+delta)*px' )/pt'
+    
+    candidate->InitialPosition.SetX (x + ((px * y - py * x + d0 * pt) / (py - px)));
+    candidate->InitialPosition.SetY (y + ((px * y - py * x + d0 * pt) / (py - px)));
+    x = candidate->InitialPosition.X ();
+    y = candidate->InitialPosition.Y ();
+    candidate->InitialPosition.SetZ (z + ((pz * (px * (x - beamSpotPosition.X ()) + py * (y - beamSpotPosition.Y ())) + pt * pt * (dz - z)) / (pt * pt)));
+  
+    candidate->InitialPosition.SetT(t);
+    
+    if (fApplyToPileUp || !candidate->IsPU)
+    {
+       candidate->ErrorD0 = d0Error;
+       candidate->ErrorDZ = dzError;
+       candidate->ErrorP = pError;
+       candidate->ErrorCtgTheta = ctgThetaError;
+       candidate->ErrorPhi = phiError;
+       candidate->ErrorPT = ptError (p, ctgTheta, pError, ctgThetaError);
+       candidate->TrackResolution = pError/p;
+    }
+    
+    candidate->AddCandidate(mother);
+    fOutputArray->Add(candidate);
+
+    iCandidate++;
+  }
+}
+
+Double_t TrackSmearing::ptError (const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta)
+{
+  Double_t a, b;
+  a = (p * p * ctgTheta * ctgTheta * dCtgTheta * dCtgTheta) / ((ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1));
+  b = (dP * dP) / (ctgTheta * ctgTheta + 1);
+  return sqrt (a + b);
+}
+
+//------------------------------------------------------------------------------
Index: modules/TrackSmearing.h
===================================================================
--- modules/TrackSmearing.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/TrackSmearing.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,72 @@
+#ifndef TrackSmearing_h
+#define TrackSmearing_h
+
+/** \class TrackSmearing
+ *
+ *  Performs d0, dZ, p, Theta, Phi smearing of tracks.
+ *
+ *
+ *
+ *  \author A. Hart, M. Selvaggi
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+class DelphesFormula;
+
+class TrackSmearing: public DelphesModule
+{
+public:
+
+  TrackSmearing();
+  ~TrackSmearing();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  Double_t ptError (const Double_t, const Double_t, const Double_t, const Double_t);
+
+  DelphesFormula *fD0Formula; //!
+  std::string fD0ResolutionFile;
+  std::string fD0ResolutionHist;
+  Bool_t fUseD0Formula;
+
+  DelphesFormula *fDZFormula; //!
+  std::string fDZResolutionFile;
+  std::string fDZResolutionHist;
+  Bool_t fUseDZFormula;
+
+  DelphesFormula *fPFormula; //!
+  std::string fPResolutionFile;
+  std::string fPResolutionHist;
+  Bool_t fUsePFormula;
+
+  DelphesFormula *fCtgThetaFormula; //!
+  std::string fCtgThetaResolutionFile;
+  std::string fCtgThetaResolutionHist;
+  Bool_t fUseCtgThetaFormula;
+
+  DelphesFormula *fPhiFormula; //!
+  std::string fPhiResolutionFile;
+  std::string fPhiResolutionHist;
+  Bool_t fUsePhiFormula;
+
+  Bool_t fApplyToPileUp;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+  const TObjArray *fBeamSpotInputArray; //!
+ 
+  TObjArray *fOutputArray; //!
+
+  ClassDef(TrackSmearing, 1)
+};
+
+#endif
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision e2339af03cbe76802d8e0c030c3ab7ca89f5c5a7)
+++ modules/TreeWriter.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -215,4 +215,11 @@
     entry->Pz = momentum.Pz();
 
+    entry->D0            = candidate->D0;
+    entry->DZ            = candidate->DZ;
+    entry->P             = candidate->P;
+    entry->PT            = candidate->PT;
+    entry->CtgTheta      = candidate->CtgTheta;
+    entry->Phi           = candidate->Phi;
+
     entry->Eta = eta;
     entry->Phi = momentum.Phi();
@@ -237,19 +244,53 @@
 
   const Double_t c_light = 2.99792458E8;
-
+ 
+  Double_t x, y, z, t, xError, yError, zError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
+  UInt_t index, ndf;
+ 
+  array->Sort();
+ 
   // loop over all vertices
   iterator.Reset();
   while((candidate = static_cast<Candidate*>(iterator.Next())))
   {
-    const TLorentzVector &position = candidate->Position;
+   
+    index = candidate->ClusterIndex;
+    ndf = candidate->ClusterNDF;
+    sigma = candidate->ClusterSigma;
+    sumPT2 = candidate->SumPT2;
+    btvSumPT2 = candidate->BTVSumPT2;
+    genDeltaZ = candidate->GenDeltaZ;
+    genSumPT2 = candidate->GenSumPT2;
+   
+    x = candidate->Position.X();
+    y = candidate->Position.Y();
+    z = candidate->Position.Z();
+    t = candidate->Position.T()*1.0E-3/c_light;
+  
+    xError = candidate->PositionError.X ();
+    yError = candidate->PositionError.Y ();
+    zError = candidate->PositionError.Z ();
 
     entry = static_cast<Vertex*>(branch->NewEntry());
 
-    entry->X = position.X();
-    entry->Y = position.Y();
-    entry->Z = position.Z();
-    entry->T = position.T()*1.0E-3/c_light;
-  }
-}
+    entry->Index = index;
+    entry->NDF = ndf;
+    entry->Sigma = sigma;
+    entry->SumPT2 = sumPT2;
+    entry->BTVSumPT2 = btvSumPT2;
+    entry->GenDeltaZ = genDeltaZ;
+    entry->GenSumPT2 = genSumPT2;
+    
+    entry->X = x;
+    entry->Y = y;
+    entry->Z = z;
+    entry->T = t;
+    
+    entry->ErrorX = xError;
+    entry->ErrorY = yError;
+    entry->ErrorZ = zError;
+  }
+}
+
 
 //------------------------------------------------------------------------------
@@ -291,7 +332,20 @@
     entry->ZOuter = position.Z();
     entry->TOuter = position.T()*1.0E-3/c_light;
-
-    entry->Dxy = candidate->Dxy;
-    entry->SDxy = candidate->SDxy ;
+ 
+    entry->L = candidate->L;
+  
+    entry->D0            = candidate->D0;
+    entry->ErrorD0       = candidate->ErrorD0;
+    entry->DZ            = candidate->DZ;
+    entry->ErrorDZ       = candidate->ErrorDZ;
+    entry->P             = candidate->P;
+    entry->ErrorP        = candidate->ErrorP;
+    entry->PT            = candidate->PT;
+    entry->ErrorPT       = candidate->ErrorPT;
+    entry->CtgTheta      = candidate->CtgTheta;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+    entry->Phi           = candidate->Phi;
+    entry->ErrorPhi      = candidate->ErrorPhi; 
+    
     entry->Xd = candidate->Xd;
     entry->Yd = candidate->Yd;
@@ -307,7 +361,5 @@
 
     entry->Eta = eta;
-    entry->Phi = momentum.Phi();
-    entry->PT = pt;
-
+    
     particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
     const TLorentzVector &initialPosition = particle->Position;
@@ -319,4 +371,7 @@
 
     entry->Particle = particle;
+
+    entry->VertexIndex = candidate->ClusterIndex;
+
   }
 }
Index: modules/VertexFinder.cc
===================================================================
--- modules/VertexFinder.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexFinder.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,376 @@
+/** \class VertexFinder
+ *
+ *  Cluster vertices from tracks
+ *
+ *  \authors A. Hart, M. Selvaggi
+ *
+ */
+
+
+#include "modules/VertexFinder.h"
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+#include "classes/DelphesPileUpReader.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TMatrixT.h"
+#include "TVector3.h"
+
+static const Double_t mm  = 1.;
+static const Double_t m = 1000.*mm;
+static const Double_t ns  = 1.;
+static const Double_t s = 1.e+9 *ns;
+static const Double_t c_light   = 2.99792458e+8 * m/s;
+
+//------------------------------------------------------------------------------
+
+VertexFinder::VertexFinder() :
+  fSigma(0), fMinPT(0), fMaxEta(0), fSeedMinPT(0), fMinNDF(0), fGrowSeeds(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+VertexFinder::~VertexFinder()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder::Init()
+{
+  fSigma = GetDouble("Sigma", 3.0);
+  fMinPT = GetDouble("MinPT", 0.1);
+  fMaxEta = GetDouble("MaxEta", 10.0);
+  fSeedMinPT = GetDouble("SeedMinPT", 5.0);
+  fMinNDF = GetInt("MinNDF", 4);
+  fGrowSeeds = GetInt("GrowSeeds", 1);
+
+  fInputArray = ImportArray(GetString("InputArray", "TrackSmearing/tracks"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+//
+Bool_t VertexFinder::secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
+{
+  return (pair0.second < pair1.second);
+}
+
+Bool_t VertexFinder::secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
+{
+  return (pair0.second > pair1.second);
+}
+
+void VertexFinder::Process()
+{
+  Candidate *candidate;
+
+////////////////////////////////////////////////////////////////////////////////
+// Clear the track and cluster maps before starting
+////////////////////////////////////////////////////////////////////////////////
+  trackIDToDouble.clear ();
+  trackIDToInt.clear ();
+  trackIDToBool.clear ();
+  clusterIDToDouble.clear ();
+  clusterIDToInt.clear ();
+  clusterIDToBool.clear ();
+  trackPT.clear ();
+  clusterSumPT2.clear ();
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// Create the initial cluster seeds
+////////////////////////////////////////////////////////////////////////////////
+  createSeeds ();
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// In order of descending seed pt, grow each cluster. If a cluster ends up with
+// fewer than MinNDF tracks, release the tracks for other clusters to claim.
+////////////////////////////////////////////////////////////////////////////////
+  sort (clusterSumPT2.begin (), clusterSumPT2.end (), secondDescending);
+  for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end (); cluster++)
+    {
+      // Skip the cluster if it no longer has any tracks
+      if (!clusterIDToInt.at (cluster->first).at ("ndf"))
+        continue;
+        
+      // Grow the cluster if GrowSeeds is true
+      if (fGrowSeeds)
+        growCluster (cluster->first);
+
+      // If the cluster still has fewer than MinNDF tracks, release the tracks;
+      // otherwise, mark the seed track as claimed
+      
+      if ((Int_t) clusterIDToInt.at (cluster->first).at ("ndf") < fMinNDF)
+        {
+          for (map<UInt_t, map<string, Int_t> >::iterator track = trackIDToInt.begin (); track != trackIDToInt.end (); track++)
+            {
+              if (track->second.at ("clusterIndex") != (Int_t) cluster->first)
+                continue;
+              track->second["clusterIndex"] = -1;
+              trackIDToBool[track->first]["claimed"] = false;
+            }
+        }
+      else
+        trackIDToBool[clusterIDToInt.at (cluster->first).at ("seed")]["claimed"] = true;
+    }
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// Add tracks to the output array after updating their ClusterIndex.
+////////////////////////////////////////////////////////////////////////////////
+  fItInputArray->Reset ();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+    {
+      if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta ()) > fMaxEta)
+        continue;
+      candidate->ClusterIndex = trackIDToInt.at (candidate->GetUniqueID ()).at ("clusterIndex");
+      fOutputArray->Add(candidate);
+    }
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// Add clusters with at least MinNDF tracks to the output array in order of
+// descending sum(pt**2).
+////////////////////////////////////////////////////////////////////////////////
+  clusterSumPT2.clear ();
+  for (map<UInt_t, map<string, Int_t> >::const_iterator cluster = clusterIDToInt.begin (); cluster != clusterIDToInt.end (); cluster++)
+  {
+    
+    if (cluster->second.at ("ndf") < fMinNDF)
+      continue;
+    clusterSumPT2.push_back (make_pair (cluster->first, clusterIDToDouble.at (cluster->first).at ("sumPT2")));
+  }
+  sort (clusterSumPT2.begin (), clusterSumPT2.end (), secondDescending);
+  
+  for (vector<pair<UInt_t, Double_t> >::const_iterator cluster = clusterSumPT2.begin (); cluster != clusterSumPT2.end (); cluster++)
+  {
+    DelphesFactory *factory = GetFactory();
+    candidate = factory->NewCandidate();
+
+    candidate->ClusterIndex = cluster->first;
+    candidate->ClusterNDF = clusterIDToInt.at (cluster->first).at ("ndf");
+    candidate->ClusterSigma = fSigma;
+    candidate->SumPT2 = cluster->second;
+    candidate->Position.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at ("z"), 0.0);
+    candidate->PositionError.SetXYZT(0.0, 0.0, clusterIDToDouble.at (cluster->first).at ("ez"), 0.0);
+
+    fVertexOutputArray->Add(candidate);
+  }
+////////////////////////////////////////////////////////////////////////////////
+}
+
+void
+VertexFinder::createSeeds ()
+{
+  Candidate *candidate;
+  UInt_t clusterIndex = 0, maxSeeds = 0;
+
+////////////////////////////////////////////////////////////////////////////////
+// Loop over all tracks, initializing some variables.
+////////////////////////////////////////////////////////////////////////////////
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+    {
+      if (candidate->Momentum.Pt () < fMinPT || fabs (candidate->Momentum.Eta ()) > fMaxEta)
+        continue;
+
+      trackIDToDouble[candidate->GetUniqueID ()]["pt"] = candidate->Momentum.Pt ();
+      trackIDToDouble[candidate->GetUniqueID ()]["ept"] = candidate->ErrorPT ? candidate->ErrorPT : 1.0e-15;;
+      trackIDToDouble[candidate->GetUniqueID ()]["eta"] = candidate->Momentum.Eta ();
+
+      trackIDToDouble[candidate->GetUniqueID ()]["z"] = candidate->DZ;
+      trackIDToDouble[candidate->GetUniqueID ()]["ez"] = candidate->ErrorDZ ? candidate->ErrorDZ : 1.0e-15;
+
+      trackIDToInt[candidate->GetUniqueID ()]["clusterIndex"] = -1;
+      trackIDToInt[candidate->GetUniqueID ()]["interactionIndex"] = candidate->IsPU;
+
+      trackIDToBool[candidate->GetUniqueID ()]["claimed"] = false;
+
+      trackPT.push_back (make_pair (candidate->GetUniqueID (), candidate->Momentum.Pt ()));
+    }
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// Sort tracks by pt and leave only the SeedMinPT highest pt ones in the
+// trackPT vector.
+////////////////////////////////////////////////////////////////////////////////
+  sort (trackPT.begin (), trackPT.end (), secondDescending);
+  for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end (); track++, maxSeeds++)
+    {
+      if (track->second < fSeedMinPT)
+        break;
+    }
+  // If there are no tracks with pt above MinSeedPT, create just one seed from
+  // the highest pt track.
+  if (!maxSeeds)
+    maxSeeds++;
+  if (trackPT.size () > maxSeeds)
+    {
+      trackPT.erase (trackPT.begin () + maxSeeds, trackPT.end ());
+    }
+////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
+// Create the seeds from the SeedMinPT highest pt tracks.
+////////////////////////////////////////////////////////////////////////////////
+  for (vector<pair<UInt_t, Double_t> >::const_iterator track = trackPT.begin (); track != trackPT.end (); track++, clusterIndex++)
+    {
+      addTrackToCluster (track->first, clusterIndex);
+      clusterSumPT2.push_back (make_pair (clusterIndex, track->second * track->second));
+    }
+////////////////////////////////////////////////////////////////////////////////
+}
+
+void
+VertexFinder::growCluster (const UInt_t clusterIndex)
+{
+  Bool_t done = false;
+  UInt_t nearestID;
+  Int_t oldClusterIndex;
+  Double_t nearestDistance;
+  vector<UInt_t> nearTracks;
+  nearTracks.clear ();
+
+////////////////////////////////////////////////////////////////////////////////
+// Grow the cluster until there are no more tracks within Sigma standard
+// deviations of the cluster.
+////////////////////////////////////////////////////////////////////////////////
+  while (!done)
+    {
+      done = true;
+      nearestID = 0;
+      nearestDistance = -1.0;
+
+      // These two loops are for finding the nearest track to the cluster. The
+      // first time, the ID of each track within 10*Sigma of the cluster is
+      // saved in the nearTracks vector; subsequently, to save time, only the
+      // tracks in this vector are checked.
+      if (!nearTracks.size ())
+        {
+        
+            for (map<UInt_t, map<string, Double_t> >::const_iterator track = trackIDToDouble.begin (); track != trackIDToDouble.end (); track++)
+            {
+              if (trackIDToBool.at (track->first).at ("claimed") || trackIDToInt.at (track->first).at ("clusterIndex") == (Int_t) clusterIndex)
+                continue;
+                
+              Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - track->second.at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), track->second.at ("ez"));
+              if (nearestDistance < 0.0 || distance < nearestDistance)
+                {
+                  nearestID = track->first;
+                  nearestDistance = distance;
+                }
+              if (distance < 10.0 * fSigma)
+                nearTracks.push_back (track->first);
+            }
+        }
+      
+      else
+        {
+          for (vector<UInt_t>::const_iterator track = nearTracks.begin (); track != nearTracks.end (); track++)
+            {
+              if (trackIDToBool.at (*track).at ("claimed") || trackIDToInt.at (*track).at ("clusterIndex") == (Int_t) clusterIndex)
+                continue;
+              Double_t distance = fabs (clusterIDToDouble.at (clusterIndex).at ("z") - trackIDToDouble.at (*track).at ("z")) / hypot (clusterIDToDouble.at (clusterIndex).at ("ez"), trackIDToDouble.at (*track).at ("ez"));
+              if (nearestDistance < 0.0 || distance < nearestDistance)
+                {
+                  nearestID = *track;
+                  nearestDistance = distance;
+                }
+            }
+        }
+      
+      // If no tracks within Sigma of the cluster were found, stop growing.
+      done = nearestDistance > fSigma || nearestDistance < 0.0;
+      if (done)
+        {
+          continue;
+        }
+
+      // Add the nearest track within Sigma to the cluster. If it already
+      // belonged to another cluster, remove it from that cluster first.
+      if (nearestDistance < fSigma)
+        {
+          oldClusterIndex = trackIDToInt.at (nearestID).at ("clusterIndex");
+          if (oldClusterIndex >= 0)
+            removeTrackFromCluster (nearestID, oldClusterIndex);
+
+          trackIDToBool[nearestID]["claimed"] = true;
+          addTrackToCluster (nearestID, clusterIndex);
+        }
+    }
+////////////////////////////////////////////////////////////////////////////////
+}
+
+Double_t
+VertexFinder::weight (const UInt_t trackID)
+{
+  return ((trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez"))) * (trackIDToDouble.at (trackID).at ("pt") / (trackIDToDouble.at (trackID).at ("ept") * trackIDToDouble.at (trackID).at ("ez"))));
+}
+
+void
+VertexFinder::removeTrackFromCluster (const UInt_t trackID, const UInt_t clusterID)
+{
+  Double_t wz = weight (trackID);
+
+  trackIDToInt[trackID]["clusterIndex"] = -1;
+  clusterIDToInt[clusterID]["ndf"]--;
+
+  clusterIDToDouble[clusterID]["sumZ"] -= wz * trackIDToDouble.at (trackID).at ("z");
+  clusterIDToDouble[clusterID]["errorSumZ"] -= wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at ("ez");
+  clusterIDToDouble[clusterID]["sumOfWeightsZ"] -= wz;
+  clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ");
+  clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ")));
+  clusterIDToDouble[clusterID]["sumPT2"] -= trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at ("pt");
+}
+
+void
+VertexFinder::addTrackToCluster (const UInt_t trackID, const UInt_t clusterID)
+{
+  Double_t wz = weight (trackID);
+
+  if (!clusterIDToInt.count (clusterID))
+    {
+      clusterIDToInt[clusterID]["ndf"] = 0;
+      clusterIDToInt[clusterID]["seed"] = trackID;
+      clusterIDToDouble[clusterID]["sumZ"] = 0.0;
+      clusterIDToDouble[clusterID]["errorSumZ"] = 0.0;
+      clusterIDToDouble[clusterID]["sumOfWeightsZ"] = 0.0;
+      clusterIDToDouble[clusterID]["sumPT2"] = 0.0;
+    }
+
+  trackIDToInt[trackID]["clusterIndex"] = clusterID;
+  clusterIDToInt[clusterID]["ndf"]++;
+
+  clusterIDToDouble[clusterID]["sumZ"] += wz * trackIDToDouble.at (trackID).at ("z");
+  clusterIDToDouble[clusterID]["errorSumZ"] += wz * trackIDToDouble.at (trackID).at ("ez") * trackIDToDouble.at (trackID).at ("ez");
+  clusterIDToDouble[clusterID]["sumOfWeightsZ"] += wz;
+  clusterIDToDouble[clusterID]["z"] = clusterIDToDouble.at (clusterID).at ("sumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ");
+  clusterIDToDouble[clusterID]["ez"] = sqrt ((1.0 / clusterIDToInt.at (clusterID).at ("ndf")) * (clusterIDToDouble.at (clusterID).at ("errorSumZ") / clusterIDToDouble.at (clusterID).at ("sumOfWeightsZ")));
+  clusterIDToDouble[clusterID]["sumPT2"] += trackIDToDouble.at (trackID).at ("pt") * trackIDToDouble.at (trackID).at ("pt");
+}
+
+//------------------------------------------------------------------------------
Index: modules/VertexFinder.h
===================================================================
--- modules/VertexFinder.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexFinder.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,76 @@
+#ifndef VertexFinder_h
+#define VertexFinder_h
+
+/** \class VertexFinder
+ *
+ *  Cluster vertices from tracks
+ *
+ *  \authors A. Hart, M. Selvaggi
+ *
+ */
+
+
+#include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+#include <utility>
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <vector>
+#include <unordered_map>
+
+using namespace std;
+
+class TObjArray;
+class Candidate;
+class TVector3;
+
+class VertexFinder: public DelphesModule
+{
+public:
+
+  VertexFinder();
+  ~VertexFinder();
+
+  void Init();
+  void Process();
+  void Finish();
+
+  static Bool_t secondDescending (pair<UInt_t, Double_t>, pair<UInt_t, Double_t>);
+  static Bool_t secondAscending (pair<UInt_t, Double_t>, pair<UInt_t, Double_t>);
+
+private:
+
+  void createSeeds ();
+  void growCluster (const UInt_t);
+  Double_t weight (const UInt_t);
+  void addTrackToCluster (const UInt_t, const UInt_t);
+  void removeTrackFromCluster (const UInt_t, const UInt_t);
+
+  Double_t fSigma;
+  Double_t fMinPT;
+  Double_t fMaxEta;
+  Double_t fSeedMinPT;
+  Int_t fMinNDF;
+  Int_t fGrowSeeds;
+
+  TObjArray *fInputArray;
+  TIterator *fItInputArray;
+
+  TObjArray *fOutputArray;
+  TObjArray *fVertexOutputArray;
+
+  map<UInt_t, map<string, Double_t> > trackIDToDouble;
+  map<UInt_t, map<string, Int_t> > trackIDToInt;
+  map<UInt_t, map<string, Bool_t> > trackIDToBool;
+
+  map<UInt_t, map<string, Double_t> > clusterIDToDouble;
+  map<UInt_t, map<string, Int_t> > clusterIDToInt;
+  map<UInt_t, map<string, Bool_t> > clusterIDToBool;
+  vector<pair<UInt_t, Double_t> > trackPT;
+  vector<pair<UInt_t, Double_t> > clusterSumPT2;
+
+  ClassDef(VertexFinder, 1)
+};
+
+#endif
Index: modules/VertexFinder4D.cc
===================================================================
--- modules/VertexFinder4D.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexFinder4D.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,1015 @@
+/** \class VertexFinder4D
+ *
+ *  Cluster vertices from tracks using deterministic annealing and timing information
+ *
+ *  \authors M. Selvaggi, L. Gray
+ *
+ */
+
+
+#include "modules/VertexFinder4D.h"
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+#include "classes/DelphesPileUpReader.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TMatrixT.h"
+#include "TVector3.h"
+
+static const Double_t mm  = 1.;
+static const Double_t m = 1000.*mm;
+static const Double_t ns  = 1.;
+static const Double_t s = 1.e+9 *ns;
+static const Double_t c_light   = 2.99792458e+8 * m/s;
+
+
+namespace {
+  bool recTrackLessZ1(const VertexFinder4D::track_t & tk1,
+                      const VertexFinder4D::track_t & tk2)
+  {
+    return tk1.z < tk2.z;
+  }
+}
+
+
+//------------------------------------------------------------------------------
+
+VertexFinder4D::VertexFinder4D() :
+  fVerbose(0), fMinPT(0), fVertexSpaceSize(0), fVertexTimeSize(0),
+  fUseTc(0), fBetaMax(0), fBetaStop(0), fCoolingFactor(0),
+  fMaxIterations(0), fDzCutOff(0), fD0CutOff(0), fDtCutOff(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+VertexFinder4D::~VertexFinder4D()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder4D::Init()
+{
+  
+  fVerbose = GetBool("Verbose", 1);
+  fMinPT = GetDouble("MinPT", 0.1);
+  fVertexSpaceSize = GetDouble("VertexSpaceSize", 0.5); //in mm
+  fVertexTimeSize = GetDouble("VertexTimeSize", 10E-12); //in s
+  fUseTc         = GetBool("UseTc", 1);
+  fBetaMax       = GetDouble("BetaMax ", 0.1);
+  fBetaStop      = GetDouble("BetaStop", 1.0);
+  fCoolingFactor = GetDouble("CoolingFactor", 0.8);
+  fMaxIterations = GetInt("MaxIterations", 100);
+  fDzCutOff      = GetDouble("DzCutOff", 40);  // Adaptive Fitter uses 30 mm but that appears to be a bit tight here sometimes
+  fD0CutOff      = GetDouble("D0CutOff", 30);
+  fDtCutOff      = GetDouble("DtCutOff", 100E-12);  // dummy
+  
+  // convert stuff in cm, ns
+  fVertexSpaceSize /= 10.0;  
+  fVertexTimeSize *= 1E9;  
+  fDzCutOff       /= 10.0;   // Adaptive Fitter uses 3.0 but that appears to be a bit tight here sometimes
+  fD0CutOff       /= 10.0; 
+  
+  fInputArray = ImportArray(GetString("InputArray", "TrackSmearing/tracks"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder4D::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder4D::Process()
+{
+  Candidate *candidate, *track;
+  TObjArray *ClusterArray;
+  ClusterArray = new TObjArray;
+  TIterator *ItClusterArray;
+  Int_t ivtx = 0;
+  
+  fInputArray->Sort();
+  
+  TLorentzVector pos, mom;
+  if (fVerbose)
+  {
+     cout<<" start processing vertices ..."<<endl;
+     cout<<" Found "<<fInputArray->GetEntriesFast()<<" input tracks"<<endl;
+     //loop over input tracks
+     fItInputArray->Reset();
+     while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+     {  
+        pos = candidate->InitialPosition;
+        mom = candidate->Momentum;
+       
+        cout<<"pt: "<<mom.Pt()<<", eta: "<<mom.Eta()<<", phi: "<<mom.Phi()<<", z: "<<candidate->DZ/10<<endl;
+     }
+  }   
+     
+  // clusterize tracks in Z
+  clusterize(*fInputArray, *ClusterArray);
+ 
+  if (fVerbose){std::cout <<  " clustering returned  "<< ClusterArray->GetEntriesFast() << " clusters  from " << fInputArray->GetEntriesFast() << " selected tracks" <<std::endl;}
+ 
+  //loop over vertex candidates
+  ItClusterArray = ClusterArray->MakeIterator();
+  ItClusterArray->Reset();
+  while((candidate = static_cast<Candidate*>(ItClusterArray->Next())))
+  {
+  
+     double meantime = 0.;
+     double expv_x2 = 0.;
+     double normw = 0.;      
+     double errtime = 0;
+     
+     double meanpos = 0.;
+     double meanerr2 = 0.;
+     double normpos = 0.;      
+     double errpos = 0.;
+     
+     double sumpt2 = 0.;
+    
+     int itr = 0;
+     
+     if(fVerbose)cout<<"this vertex has: "<<candidate->GetCandidates()->GetEntriesFast()<<" tracks"<<endl;
+     
+     // loop over tracks belonging to this vertex
+     TIter it1(candidate->GetCandidates());
+     it1.Reset();
+     
+     while((track = static_cast<Candidate*>(it1.Next())))
+     {
+        
+        itr++;
+        // TBC: the time is in ns for now TBC 
+        double t = track->InitialPosition.T()/c_light;
+        double dt = track->ErrorT/c_light;
+        const double time = t;
+        const double inverr = 1.0/dt;
+        meantime += time*inverr;
+        expv_x2  += time*time*inverr;
+        normw    += inverr;
+     
+        // compute error position TBC
+        const double pt = track->Momentum.Pt();
+        const double z = track->DZ/10.0;
+        const double err_pt = track->ErrorPT;
+        const double err_z = track->ErrorDZ;
+        
+        const double wi = (pt/(err_pt*err_z))*(pt/(err_pt*err_z));
+        meanpos += z*wi;
+        
+        meanerr2 += err_z*err_z*wi;
+        normpos += wi;
+        sumpt2 += pt*pt;
+      
+        // while we are here store cluster index in tracks
+        track->ClusterIndex = ivtx;  
+     }
+     
+     meantime = meantime/normw;
+     expv_x2 = expv_x2/normw;
+     errtime = TMath::Sqrt((expv_x2 - meantime*meantime)/itr); 
+
+     meanpos = meanpos/normpos;
+     meanerr2 = meanerr2/normpos;
+     errpos = TMath::Sqrt(meanerr2/itr);
+     
+     candidate->Position.SetXYZT(0.0, 0.0, meanpos*10.0 , meantime*c_light);
+     candidate->PositionError.SetXYZT(0.0, 0.0, errpos*10.0 , errtime*c_light);
+     candidate->SumPT2 = sumpt2;
+     candidate->ClusterNDF = itr;
+     candidate->ClusterIndex = ivtx;
+
+     fVertexOutputArray->Add(candidate);
+     
+     ivtx++;
+ 
+     if (fVerbose){
+	   std::cout << "x,y,z";
+       std::cout << ",t";
+       std::cout << "=" << candidate->Position.X()/10.0 <<" " << candidate->Position.Y()/10.0 << " " <<  candidate->Position.Z()/10.0;
+       std::cout << " " << candidate->Position.T()/c_light;
+       
+       std::cout << std::endl;
+       std::cout << "sumpt2 " << candidate->SumPT2<<endl;
+     }
+   }// end of cluster loop
+
+    
+    if(fVerbose){
+      std::cout << "PrimaryVertexProducerAlgorithm::vertices candidates =" << ClusterArray->GetEntriesFast() << std::endl;
+    }
+
+    //TBC maybe this can be done later  
+    // sort vertices by pt**2  vertex (aka signal vertex tagging)
+    /*if(pvs.size()>1){
+      sort(pvs.begin(), pvs.end(), VertexHigherPtSquared());
+    }
+     */
+        
+  delete ClusterArray; 
+
+}
+
+//------------------------------------------------------------------------------
+
+
+void VertexFinder4D::clusterize(const TObjArray & tracks, TObjArray & clusters)
+{
+  if(fVerbose) {
+    cout << "###################################################" << endl;
+    cout << "# VertexFinder4D::clusterize   nt="<<tracks.GetEntriesFast() << endl;
+    cout << "###################################################" << endl;
+  }
+ 
+  vector< Candidate* > pv = vertices();
+
+  if(fVerbose){ cout << "# VertexFinder4D::clusterize   pv.size="<<pv.size() << endl;  }
+  if (pv.size()==0){ return;  }
+
+  // convert into vector of candidates
+  //TObjArray *ClusterArray = pv.begin()->GetCandidates();
+  //Candidate *aCluster = static_cast<Candidate*>(&(pv.at(0)));
+   Candidate *aCluster = pv.at(0);
+   
+  // fill into clusters and merge
+   
+  
+  if( fVerbose ) { 
+      std::cout << '\t' << 0;
+      std::cout << ' ' << (*pv.begin())->Position.Z()/10.0 << ' ' << (*pv.begin())->Position.T()/c_light << std::endl; 
+    }
+
+  for(vector<Candidate*>::iterator k=pv.begin()+1; k!=pv.end(); k++){
+    if( fVerbose ) { 
+      std::cout << '\t' << std::distance(pv.begin(),k);
+      std::cout << ' ' << (*k)->Position.Z() << ' ' << (*k)->Position.T() << std::endl; 
+    }
+    
+    
+    // TBC - check units here
+    if ( std::abs((*k)->Position.Z() - (*(k-1))->Position.Z())/10.0 > (2*fVertexSpaceSize) ||
+         std::abs((*k)->Position.T() - (*(k-1))->Position.Z())/c_light > 2*0.010 ) {
+      // close a cluster
+      clusters.Add(aCluster);
+      //aCluster.clear();
+    }
+    //for(unsigned int i=0; i<k->GetCandidates().GetEntriesFast(); i++){ 
+      aCluster = *k; 
+    //}
+    
+  }
+  clusters.Add(aCluster);
+  
+  if(fVerbose) { std::cout << "# VertexFinder4D::clusterize clusters.size="<<clusters.GetEntriesFast() << std::endl; }
+ 
+}
+
+
+//------------------------------------------------------------------------------
+
+vector< Candidate* > VertexFinder4D::vertices() 
+{
+  Candidate *candidate;
+  UInt_t clusterIndex = 0;
+  vector< Candidate* > clusters;
+ 
+  vector<track_t> tks = fill();
+   
+  //print out input tracks
+   
+  if(fVerbose)
+  { 
+    std::cout<<" start processing vertices ..."<<std::endl;
+    std::cout<<" Found "<<tks.size()<<" input tracks"<<std::endl;
+    //loop over input tracks
+ 
+   
+   for(std::vector<track_t>::const_iterator it=tks.begin(); it!=tks.end(); it++){
+     double z = it->z;
+     double pt=it->pt;
+     double eta=it->eta;
+     double phi=it->phi;
+     double t = it->t;
+
+     std::cout<<"pt: "<<pt<<", eta: "<<eta<<", phi: "<<phi<<", z: "<<z<<", t: "<<t<<std::endl;
+   } 
+  }
+   
+  unsigned int nt=tks.size();
+  double rho0=0.0;  // start with no outlier rejection
+
+  if (tks.empty()) return clusters;
+
+  vector<vertex_t> y; // the vertex prototypes
+
+  // initialize:single vertex at infinite temperature
+  vertex_t vstart;
+  vstart.z=0.;
+  vstart.t=0.;
+  vstart.pk=1.;
+  y.push_back(vstart);
+  int niter=0;      // number of iterations
+ 
+  // estimate first critical temperature
+  double beta=beta0(fBetaMax, tks, y);
+  niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
+
+  // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
+  while(beta<fBetaMax){ 
+
+    if(fUseTc){
+      update(beta, tks,y);
+      while(merge(y,beta)){update(beta, tks,y);}
+      split(beta, tks,y);
+      beta=beta/fCoolingFactor;
+    }else{
+      beta=beta/fCoolingFactor;
+      splitAll(y);
+    }
+
+   // make sure we are not too far from equilibrium before cooling further
+   niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
+
+  }
+
+  if(fUseTc){
+    // last round of splitting, make sure no critical clusters are left
+    update(beta, tks,y);
+    while(merge(y,beta)){update(beta, tks,y);}
+    unsigned int ntry=0;
+    while( split(beta, tks,y) && (ntry++<10) ){
+      niter=0; 
+      while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){}
+      merge(y,beta);
+      update(beta, tks,y);
+    }
+  }else{
+    // merge collapsed clusters 
+    while(merge(y,beta)){update(beta, tks,y);}  
+    if(fVerbose ){ cout << "dump after 1st merging " << endl;  dump(beta,y,tks);}
+  }
+ 
+  // switch on outlier rejection
+  rho0=1./nt; for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ k->pk =1.; }  // democratic
+  niter=0; while((update(beta, tks,y,rho0) > 1.e-8)  && (niter++ < fMaxIterations)){  }
+  if(fVerbose  ){ cout << "rho0=" << rho0 <<   " niter=" << niter <<  endl; dump(beta,y,tks);}
+
+  
+  // merge again  (some cluster split by outliers collapse here)
+  while(merge(y)){}  
+  if(fVerbose  ){ cout << "dump after 2nd merging " << endl;  dump(beta,y,tks);}
+
+
+  // continue from freeze-out to Tstop (=1) without splitting, eliminate insignificant vertices
+  while(beta<=fBetaStop){
+    while(purge(y,tks,rho0, beta)){
+      niter=0; while((update(beta, tks, y, rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
+    } 
+    beta/=fCoolingFactor;
+    niter=0; while((update(beta, tks, y, rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
+  }
+
+
+//   // new, one last round of cleaning at T=Tstop
+//   while(purge(y,tks,rho0, beta)){
+//     niter=0; while((update(beta, tks,y,rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
+//   } 
+
+
+  if(fVerbose){
+   cout << "Final result, rho0=" << rho0 << endl;
+   dump(beta,y,tks);
+  }
+
+
+  // select significant tracks and use a TransientVertex as a container
+  //GlobalError dummyError; 
+  
+  // ensure correct normalization of probabilities, should make double assginment reasonably impossible
+  for(unsigned int i=0; i<nt; i++){  
+    tks[i].Z=rho0*exp(-beta*( fDzCutOff*fDzCutOff));
+    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
+      tks[i].Z += k->pk * exp(-beta*Eik(tks[i],*k));
+    }
+  }
+  
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
+    
+    DelphesFactory *factory = GetFactory();
+    candidate = factory->NewCandidate();
+    
+    //cout<<"new vertex"<<endl;
+    //GlobalPoint pos(0, 0, k->z);
+    double time = k->t;
+    double z = k->z;
+    //vector< reco::TransientTrack > vertexTracks;
+    //double max_track_time_err2 = 0;
+    double mean = 0.;
+    double expv_x2 = 0.;
+    double normw = 0.;
+    for(unsigned int i=0; i<nt; i++){
+      const double invdt = 1.0/std::sqrt(tks[i].dt2);
+      if(tks[i].Z>0){
+	double p = k->pk * exp(-beta*Eik(tks[i],*k)) / tks[i].Z;
+	if( (tks[i].pi>0) && ( p > 0.5 ) ){ 
+          //std::cout << "pushing back " << i << ' ' << tks[i].tt << std::endl;
+          //vertexTracks.push_back(*(tks[i].tt)); tks[i].Z=0; 
+          
+          candidate->AddCandidate(tks[i].tt); tks[i].Z=0; 
+         
+          mean     += tks[i].t*invdt*p;
+          expv_x2  += tks[i].t*tks[i].t*invdt*p;
+          normw    += invdt*p;
+        } // setting Z=0 excludes double assignment        
+      }
+    }
+    
+    mean = mean/normw;
+    expv_x2 = expv_x2/normw;
+    const double time_var = expv_x2 - mean*mean;    
+    const double crappy_error_guess = std::sqrt(time_var);
+    /*GlobalError dummyErrorWithTime(0,
+                                   0,0,
+                                   0,0,0,
+                                   0,0,0,crappy_error_guess);*/
+    //TransientVertex v(pos, time, dummyErrorWithTime, vertexTracks, 5);
+   
+    
+    candidate->ClusterIndex = clusterIndex++;;
+    candidate->Position.SetXYZT(0.0, 0.0, z*10.0 , time*c_light);
+    
+    // TBC - fill error later ... 
+    candidate->PositionError.SetXYZT(0.0, 0.0, 0.0 , crappy_error_guess*c_light);
+
+    clusterIndex++;   
+    clusters.push_back(candidate);
+  }
+  
+
+  return clusters;
+
+}
+
+//------------------------------------------------------------------------------
+
+vector<VertexFinder4D::track_t> 
+VertexFinder4D::fill()const{
+  
+  Candidate *candidate;
+  track_t tr;
+  Double_t z, dz, t, l, dt, d0, d0error;
+  
+  // prepare track data for clustering
+  vector< track_t > tks;
+  
+  // loop over input tracks
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+  {
+     //TBC everything in cm 
+     z = candidate->DZ/10;
+     tr.z = z;
+     dz = candidate->ErrorDZ/10;
+     tr.dz2 = dz*dz          // track error
+   // TBC: beamspot size induced error, take 0 for now.        
+   //   + (std::pow(beamspot.BeamWidthX()*cos(phi),2.)+std::pow(beamspot.BeamWidthY()*sin(phi),2.))/std::pow(tantheta,2.) // beam-width induced
+      + fVertexSpaceSize*fVertexSpaceSize;                     // intrinsic vertex size, safer for outliers and short lived decays
+    
+     // TBC: the time is in ns for now TBC 
+     //t = candidate->Position.T()/c_light;
+     t = candidate->InitialPosition.T()/c_light;
+     l = candidate->L/c_light;
+     
+     double pt = candidate->Momentum.Pt();
+     double eta = candidate->Momentum.Eta();
+     double phi = candidate->Momentum.Phi();
+
+     tr.pt = pt;
+     tr.eta = eta;
+     tr.phi = phi;
+     tr.t = t; // 
+     tr.dtz = 0.;
+     dt = candidate->ErrorT/c_light;
+     tr.dt2 = dt*dt + fVertexTimeSize*fVertexTimeSize;   // the ~injected~ timing error plus a small minimum vertex size in time
+     if (fD0CutOff>0){
+    
+     d0 = TMath::Abs(candidate->D0)/10.0;
+     d0error = candidate->ErrorD0/10.0;
+      
+     tr.pi=1./(1.+exp((d0*d0)/(d0error*d0error) - fD0CutOff*fD0CutOff));  // reduce weight for high ip tracks  
+     
+    }
+    else{
+      tr.pi=1.;
+    }
+    tr.tt=&(*candidate);    
+    tr.Z=1.;
+ 
+    // TBC now putting track selection here (> fPTMin)
+    if( tr.pi > 1e-3 && tr.pt > fMinPT) {
+      tks.push_back(tr);
+    }
+ 
+  }
+  
+  return tks;
+}
+
+
+
+
+//------------------------------------------------------------------------------
+
+double VertexFinder4D::Eik(const track_t & t, const vertex_t &k ) const {
+  return std::pow(t.z-k.z,2.)/t.dz2 + std::pow(t.t - k.t,2.)/t.dt2;
+}
+  
+//--------------------------------------------------------------------------------  
+
+void VertexFinder4D::dump(const double beta, const vector<vertex_t> & y, const vector<track_t> & tks0)const{
+
+  // copy and sort for nicer printout
+  vector<track_t> tks; 
+  for(vector<track_t>::const_iterator t=tks0.begin(); t!=tks0.end(); t++){tks.push_back(*t); }
+  std::stable_sort(tks.begin(), tks.end(), recTrackLessZ1);
+
+  cout << "-----DAClusterizerInZT::dump ----" << endl;
+  cout << " beta=" << beta << "   betamax= " << fBetaMax << endl;
+  cout << "                                                               z= ";
+  cout.precision(4);
+  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
+    cout  <<  setw(8) << fixed << k->z;
+  }
+  cout << endl << "                                                               t= ";
+  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
+    cout  <<  setw(8) << fixed << k->t;
+  }
+  cout << endl << "T=" << setw(15) << 1./beta <<"                                             Tc= ";
+  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
+    cout  <<  setw(8) << fixed << k->Tc ;
+  }
+ 
+  cout << endl << "                                                              pk=";
+  double sumpk=0;
+  for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
+    cout <<  setw(8) <<  setprecision(3) <<  fixed << k->pk;
+    sumpk+=k->pk;
+  }
+  cout  << endl;
+
+  if(fVerbose){
+    double E=0, F=0;
+    cout << endl;
+    cout << "----       z +/- dz        t +/- dt        ip +/-dip       pt    phi  eta    weights  ----" << endl;
+    cout.precision(4);
+    for(unsigned int i=0; i<tks.size(); i++){
+      if (tks[i].Z>0){	F-=log(tks[i].Z)/beta;}
+      double tz= tks[i].z;
+      double tt= tks[i].t;
+      cout <<  setw (3)<< i << ")" <<  setw (8) << fixed << setprecision(4)<<  tz << " +/-" <<  setw (6)<< sqrt(tks[i].dz2) 
+           << setw(8) << fixed << setprecision(4) << tt << " +/-" << setw(6) << std::sqrt(tks[i].dt2)  ;
+      
+      double sump=0.;
+      for(vector<vertex_t>::const_iterator k=y.begin(); k!=y.end(); k++){
+	if((tks[i].pi>0)&&(tks[i].Z>0)){
+	  //double p=pik(beta,tks[i],*k);
+	  double p=k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z; 
+	  if( p > 0.0001){
+	    cout <<  setw (8) <<  setprecision(3) << p;
+	  }else{
+	    cout << "    .   ";
+	  }
+	  E+=p*Eik(tks[i],*k);
+	  sump+=p;
+	}else{
+	    cout << "        ";
+	}
+      }
+      cout << endl;
+    }
+    cout << endl << "T=" << 1/beta  << " E=" << E << " n="<< y.size() << "  F= " << F <<  endl <<  "----------" << endl;
+  }
+}
+
+
+  
+  
+
+//------------------------------------------------------------------------------
+
+
+double VertexFinder4D::update( double beta,
+                                  vector<track_t> & tks,
+                                  vector<vertex_t> & y ) const {
+  //update weights and vertex positions
+  // mass constrained annealing without noise
+  // returns the squared sum of changes of vertex positions
+
+  unsigned int nt=tks.size();
+
+  //initialize sums
+  double sumpi=0;
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); ++k){
+    k->sw=0.; k->swz=0.; k->swt = 0.; k->se=0.;
+    k->swE=0.;  k->Tc=0.;
+  }
+
+
+  // loop over tracks
+  for(unsigned int i=0; i<nt; i++){
+
+    // update pik and Zi
+    double Zi = 0.;
+    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); ++k){
+      k->ei = std::exp(-beta*Eik(tks[i],*k));// cache exponential for one track at a time
+      Zi   += k->pk * k->ei;      
+    }
+    tks[i].Z=Zi;
+
+    // normalization for pk
+    if (tks[i].Z>0){
+      sumpi += tks[i].pi;
+      // accumulate weighted z and weights for vertex update
+      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); ++k){
+	k->se  += tks[i].pi* k->ei / Zi;
+	const double w = k->pk * tks[i].pi* k->ei / ( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
+	k->sw  += w;
+	k->swz += w * tks[i].z;
+        k->swt += w * tks[i].t;
+	k->swE += w * Eik(tks[i],*k);
+      }
+    }else{
+      sumpi += tks[i].pi;
+    }
+    
+
+  } // end of track loop
+
+
+  // now update z and pk
+  double delta=0;
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+    if ( k->sw > 0){
+      const double znew = k->swz/k->sw; 
+      const double tnew = k->swt/k->sw;
+      delta += std::pow(k->z-znew,2.) + std::pow(k->t-tnew,2.);
+      k->z  = znew;
+      k->t  = tnew;
+      k->Tc = 2.*k->swE/k->sw;
+    }else{
+      if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
+      k->Tc=-1;
+    }
+
+    k->pk = k->pk * k->se / sumpi;
+  }
+
+  // return how much the prototypes moved
+  return delta;
+}
+
+//------------------------------------------------------------------------------
+
+double VertexFinder4D::update( double beta,
+                                  vector<track_t> & tks,
+                                  vector<vertex_t> & y,
+                                  double & rho0 ) const {
+  // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise 
+  // returns the squared sum of changes of vertex positions
+
+  unsigned int nt=tks.size();
+
+  //initialize sums
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+    k->sw = 0.;   k->swz = 0.; k->swt = 0.; k->se = 0.;    
+    k->swE = 0.;  k->Tc=0.;
+  }
+
+
+  // loop over tracks
+  for(unsigned int i=0; i<nt; i++){
+
+    // update pik and Zi and Ti
+    double Zi = rho0*std::exp(-beta*(fDzCutOff*fDzCutOff));// cut-off (eventually add finite size in time)
+    //double Ti = 0.; // dt0*std::exp(-beta*fDtCutOff);
+    for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+      k->ei = std::exp(-beta*Eik(tks[i],*k));// cache exponential for one track at a time
+      Zi   += k->pk * k->ei;
+    }
+    tks[i].Z=Zi;
+    
+    // normalization
+    if (tks[i].Z>0){
+       // accumulate weighted z and weights for vertex update
+      for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+	k->se += tks[i].pi* k->ei / Zi;
+	double w = k->pk * tks[i].pi * k->ei /( Zi * ( tks[i].dz2 * tks[i].dt2 ) );
+	k->sw  += w;
+	k->swz += w * tks[i].z;
+        k->swt += w * tks[i].t;
+	k->swE += w * Eik(tks[i],*k);
+      }
+    }
+
+  } // end of track loop
+
+  // now update z
+  double delta=0;
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+    if ( k->sw > 0){
+      const double znew=k->swz/k->sw; 
+      const double tnew=k->swt/k->sw;
+      delta += std::pow(k->z-znew,2.) + std::pow(k->t-tnew,2.);
+      k->z   = znew;
+      k->t   = tnew;
+      k->Tc  = 2*k->swE/k->sw;
+    }else{
+      if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
+      k->Tc = 0;
+    }
+
+  }
+
+  // return how much the prototypes moved
+  return delta;
+}
+
+//------------------------------------------------------------------------------
+
+bool VertexFinder4D::merge(vector<vertex_t> & y)const{
+  // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
+  
+  if(y.size()<2)  return false;
+  
+  for(vector<vertex_t>::iterator k=y.begin(); (k+1)!=y.end(); k++){
+    if( std::abs( (k+1)->z - k->z ) < 1.e-3 &&
+        std::abs( (k+1)->t - k->t ) < 1.e-3    ){  // with fabs if only called after freeze-out (splitAll() at highter T)
+      double rho = k->pk + (k+1)->pk;
+      if(rho>0){
+        k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho;
+        k->t = ( k->pk * k->t + (k+1)->t * (k+1)->pk)/rho;
+      }else{
+        k->z = 0.5*(k->z + (k+1)->z);
+        k->t = 0.5*(k->t + (k+1)->t);
+      }
+      k->pk = rho;
+      
+      y.erase(k+1);
+      return true;  
+    }
+  }
+  
+  return false;
+}
+
+//------------------------------------------------------------------------------
+
+bool VertexFinder4D::merge(vector<vertex_t> & y, double & beta)const{
+  // merge clusters that collapsed or never separated, 
+  // only merge if the estimated critical temperature of the merged vertex is below the current temperature
+  // return true if vertices were merged, false otherwise
+  if(y.size()<2)  return false;
+
+  for(vector<vertex_t>::iterator k=y.begin(); (k+1)!=y.end(); k++){
+    if ( std::abs((k+1)->z - k->z) < 2.e-3 && 
+         std::abs((k+1)->t - k->t) < 2.e-3    ) { 
+      double rho=k->pk + (k+1)->pk;
+      double swE=k->swE+(k+1)->swE - k->pk * (k+1)->pk / rho * ( std::pow((k+1)->z - k->z,2.) + 
+                                                                 std::pow((k+1)->t - k->t,2.)   );
+      double Tc=2*swE/(k->sw+(k+1)->sw);
+      
+      if(Tc*beta<1){
+	if(rho>0){
+	  k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho;
+          k->t = ( k->pk * k->t + (k+1)->t * (k+1)->pk)/rho;
+	}else{
+	  k->z = 0.5*(k->z + (k+1)->z);
+          k->t = 0.5*(k->t + (k+1)->t);
+	}
+	k->pk  = rho;
+	k->sw += (k+1)->sw;
+	k->swE = swE;
+	k->Tc  = Tc;
+	y.erase(k+1);
+	return true; 
+      }
+    }
+  }
+
+  return false;
+}
+
+//------------------------------------------------------------------------------
+
+bool VertexFinder4D::purge(vector<vertex_t> & y, vector<track_t> & tks, double & rho0, const double beta)const{
+  // eliminate clusters with only one significant/unique track
+  if(y.size()<2)  return false;
+  
+  unsigned int nt=tks.size();
+  double sumpmin=nt;
+  vector<vertex_t>::iterator k0=y.end();
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){ 
+    int nUnique=0;
+    double sump=0;
+    double pmax=k->pk/(k->pk+rho0*exp(-beta*fDzCutOff*fDzCutOff));
+    for(unsigned int i=0; i<nt; i++){
+      if(tks[i].Z > 0){
+	double p = k->pk * std::exp(-beta*Eik(tks[i],*k)) / tks[i].Z ;
+	sump+=p;
+	if( (p > 0.9*pmax) && (tks[i].pi>0) ){ nUnique++; }
+      }
+    }
+
+    if((nUnique<2)&&(sump<sumpmin)){
+      sumpmin=sump;
+      k0=k;
+    }
+  }
+ 
+  if(k0!=y.end()){
+    if(fVerbose){cout << "eliminating prototype at " << k0->z << "," << k0->t 
+                      << " with sump=" << sumpmin << endl;}
+    //rho0+=k0->pk;
+    y.erase(k0);
+    return true;
+  }else{
+    return false;
+  }
+}
+
+//------------------------------------------------------------------------------
+
+double VertexFinder4D::beta0( double betamax,
+                                 vector<track_t> & tks,
+                                 vector<vertex_t> & y ) const {
+  
+  double T0=0;  // max Tc for beta=0
+  // estimate critical temperature from beta=0 (T=inf)
+  unsigned int nt=tks.size();
+
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+
+    // vertex fit at T=inf 
+    double sumwz=0.;
+    double sumwt=0.;
+    double sumw=0.;
+    for(unsigned int i=0; i<nt; i++){
+      double w = tks[i].pi/(tks[i].dz2 * tks[i].dt2);
+      sumwz += w*tks[i].z;
+      sumwt += w*tks[i].t;
+      sumw  += w;
+    }
+    k->z = sumwz/sumw;
+    k->t = sumwt/sumw;
+
+    // estimate Tcrit, eventually do this in the same loop
+    double a=0, b=0;
+    for(unsigned int i=0; i<nt; i++){
+      double dx = tks[i].z-(k->z);
+      double dt = tks[i].t-(k->t);
+      double w  = tks[i].pi/(tks[i].dz2 * tks[i].dt2);
+      a += w*(std::pow(dx,2.)/tks[i].dz2 + std::pow(dt,2.)/tks[i].dt2);
+      b += w;
+    }
+    double Tc= 2.*a/b;  // the critical temperature of this vertex
+    if(Tc>T0) T0=Tc;
+  }// vertex loop (normally there should be only one vertex at beta=0)
+  
+  if (T0>1./betamax){
+    return betamax/pow(fCoolingFactor, int(std::log(T0*betamax)/std::log(fCoolingFactor))-1 );
+  }else{
+    // ensure at least one annealing step
+    return betamax/fCoolingFactor;
+  }
+}
+
+//------------------------------------------------------------------------------
+
+bool VertexFinder4D::split( double beta,
+                               vector<track_t> & tks,
+                               vector<vertex_t> & y) const {
+  // split only critical vertices (Tc >~ T=1/beta   <==>   beta*Tc>~1)
+  // an update must have been made just before doing this (same beta, no merging)
+  // returns true if at least one cluster was split
+
+  constexpr double epsilon=1e-3;      // split all single vertices by 10 um 
+  bool split=false;
+  
+  // avoid left-right biases by splitting highest Tc first
+  
+  std::vector<std::pair<double, unsigned int> > critical;
+  for(unsigned int ik=0; ik<y.size(); ik++){
+    if (beta*y[ik].Tc > 1.){
+      critical.push_back( make_pair(y[ik].Tc, ik));
+    }
+  }
+  std::stable_sort(critical.begin(), critical.end(), std::greater<std::pair<double, unsigned int> >() );
+
+  for(unsigned int ic=0; ic<critical.size(); ic++){
+    unsigned int ik=critical[ic].second;
+    // estimate subcluster positions and weight
+    double p1=0, z1=0, t1=0, w1=0;
+    double p2=0, z2=0, t2=0, w2=0;
+    //double sumpi=0;
+    for(unsigned int i=0; i<tks.size(); i++){
+      if(tks[i].Z>0){
+	//sumpi+=tks[i].pi;
+	double p=y[ik].pk * exp(-beta*Eik(tks[i],y[ik])) / tks[i].Z*tks[i].pi;
+	double w=p/(tks[i].dz2 * tks[i].dt2);
+	if(tks[i].z < y[ik].z){
+	  p1+=p; z1+=w*tks[i].z; t1+=w*tks[i].t; w1+=w;
+	}else{
+	  p2+=p; z2+=w*tks[i].z; t2+=w*tks[i].t; w2+=w;
+	}
+      }
+    }
+    if(w1>0){  z1=z1/w1; t1=t1/w1;} else{ z1=y[ik].z-epsilon; t1=y[ik].t-epsilon; }
+    if(w2>0){  z2=z2/w2; t2=t2/w2;} else{ z2=y[ik].z+epsilon; t2=y[ik].t+epsilon;}
+
+    // reduce split size if there is not enough room
+    if( ( ik   > 0       ) && ( y[ik-1].z>=z1 ) ){ z1=0.5*(y[ik].z+y[ik-1].z); t1=0.5*(y[ik].t+y[ik-1].t); }
+    if( ( ik+1 < y.size()) && ( y[ik+1].z<=z2 ) ){ z2=0.5*(y[ik].z+y[ik+1].z); t2=0.5*(y[ik].t+y[ik+1].t); }
+
+    // split if the new subclusters are significantly separated
+    if( (z2-z1)>epsilon || std::abs(t2-t1) > epsilon){
+      split=true;
+      vertex_t vnew;
+      vnew.pk = p1*y[ik].pk/(p1+p2);
+      y[ik].pk= p2*y[ik].pk/(p1+p2);
+      vnew.z  = z1;
+      vnew.t  = t1;
+      y[ik].z = z2;
+      y[ik].t = t2;
+      y.insert(y.begin()+ik, vnew);
+
+     // adjust remaining pointers
+      for(unsigned int jc=ic; jc<critical.size(); jc++){
+	if (critical[jc].second>ik) {critical[jc].second++;}
+      }
+    }
+  }
+
+  //  stable_sort(y.begin(), y.end(), clusterLessZ);
+  return split;
+}
+
+//------------------------------------------------------------------------------
+
+void VertexFinder4D::splitAll( vector<vertex_t> & y ) const {
+
+
+  constexpr double epsilon=1e-3;      // split all single vertices by 10 um 
+  constexpr double zsep=2*epsilon;    // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
+  constexpr double tsep=2*epsilon;    // check t as well
+  
+  vector<vertex_t> y1;
+
+  for(vector<vertex_t>::iterator k=y.begin(); k!=y.end(); k++){
+    if ( ( (k==y.begin())|| (k-1)->z < k->z - zsep) && (((k+1)==y.end()  )|| (k+1)->z > k->z + zsep)) { 
+      // isolated prototype, split
+      vertex_t vnew;
+      vnew.z  = k->z - epsilon;
+      vnew.t  = k->t - epsilon;
+      (*k).z  = k->z + epsilon;
+      (*k).t  = k->t + epsilon;
+      vnew.pk= 0.5* (*k).pk;
+      (*k).pk= 0.5* (*k).pk;
+      y1.push_back(vnew);
+      y1.push_back(*k);
+
+    }else if( y1.empty() || (y1.back().z < k->z -zsep) || (y1.back().t < k->t - tsep) ){
+      y1.push_back(*k);
+    }else{
+      y1.back().z -= epsilon;
+      y1.back().t -= epsilon;
+      k->z += epsilon;
+      k->t += epsilon;
+      y1.push_back(*k);
+    }
+  }// vertex loop
+  
+  y=y1;
+}
+ 
+ 
+
Index: modules/VertexFinder4D.h
===================================================================
--- modules/VertexFinder4D.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexFinder4D.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,132 @@
+#ifndef VertexFinder4D_h
+#define VertexFinder4D_h
+
+/** \class VertexFinder4D
+ *
+ *  Cluster vertices from tracks using deterministic annealing and timing information
+ *
+ *  \authors M. Selvaggi, L. Gray
+ *
+ */
+
+
+
+#include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+#include <utility>
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <vector>
+#include <unordered_map>
+
+using namespace std;
+
+class TObjArray;
+class Candidate;
+class TVector3;
+
+class VertexFinder4D: public DelphesModule
+{
+public:
+
+  VertexFinder4D();
+  ~VertexFinder4D();
+
+  void Init();
+  void Process();
+  void Finish();
+
+  struct track_t{
+  double z;              // z-coordinate at point of closest approach to the beamline
+  double t;              // t-coordinate at point of closest approach to the beamline
+  double dz2;            // square of the error of z(pca)
+  double dtz;            // covariance of z-t
+  double dt2;            // square of the error of t(pca)
+  //const reco::TransientTrack* tt;  // a pointer to the Transient Track
+  Candidate* tt;  // a pointer to the Candidate Track
+  double Z;              // Z[i]   for DA clustering
+  double pi;             // track weight
+  double pt;
+  double eta;
+  double phi;
+
+};
+
+
+struct vertex_t{
+  double z;    //           z coordinate
+  double t;    //           t coordinate
+  double pk;   //           vertex weight for "constrained" clustering
+  // --- temporary numbers, used during update
+  double ei;
+  double sw;
+  double swz;
+  double swt;
+  double se;
+  // ---for Tc
+  double swE;
+  double Tc;
+};
+
+  
+
+  void clusterize(const TObjArray & tracks, TObjArray & clusters);
+
+  std::vector< Candidate* > vertices();
+
+  std::vector<track_t> fill() const;
+
+  bool split( double beta,
+             std::vector<track_t> & tks,
+             std::vector<vertex_t> & y) const;
+
+  double update( double beta,
+                std::vector<track_t> & tks,
+                std::vector<vertex_t> & y ) const;
+
+  double update(double beta,
+               std::vector<track_t> & tks,
+               std::vector<vertex_t> & y,
+               double & )const;
+
+  void dump(const double beta, const std::vector<vertex_t> & y, const std::vector<track_t> & tks) const;
+  bool merge(std::vector<vertex_t> &) const;
+  bool merge(std::vector<vertex_t> &,double & ) const;
+  bool purge(std::vector<vertex_t> &, std::vector<track_t> & , double &, const double ) const;
+
+  void splitAll( std::vector<vertex_t> & y ) const;
+
+  double beta0(const double betamax,
+	       std::vector<track_t> & tks,
+	       std::vector<vertex_t> & y )const;
+
+  double Eik(const track_t & t, const vertex_t & k)const;
+
+
+private:
+
+  Bool_t fVerbose;
+  Double_t fMinPT;
+  
+  Float_t fVertexSpaceSize;
+  Float_t fVertexTimeSize;
+  Bool_t fUseTc;
+  Float_t fBetaMax;
+  Float_t fBetaStop;
+  Double_t fCoolingFactor;
+  Int_t fMaxIterations;
+  Double_t fDzCutOff;
+  Double_t fD0CutOff;
+  Double_t fDtCutOff; // for when the beamspot has time
+
+  TObjArray *fInputArray;
+  TIterator *fItInputArray;
+
+  TObjArray *fOutputArray;
+  TObjArray *fVertexOutputArray;
+
+  ClassDef(VertexFinder4D, 1)
+};
+
+#endif
Index: modules/VertexSorter.cc
===================================================================
--- modules/VertexSorter.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexSorter.cc	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,217 @@
+/** \class VertexSorter
+ *
+ *
+ *  Sorts vertices according to different criteria
+ *
+ *  \authors A. Hart, M. Selvaggi
+ *
+ *
+*/
+
+#include "modules/VertexSorter.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+#include "classes/DelphesPileUpReader.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TMatrixT.h"
+#include "TVector3.h"
+
+static const Double_t mm  = 1.;
+static const Double_t m = 1000.*mm;
+static const Double_t ns  = 1.;
+static const Double_t s = 1.e+9 *ns;
+static const Double_t c_light   = 2.99792458e+8 * m/s;
+
+//------------------------------------------------------------------------------
+
+VertexSorter::VertexSorter() :
+  fInputArray(NULL), fTrackInputArray(NULL), fItTrackInputArray(NULL), fJetInputArray(NULL), fItJetInputArray(NULL), fOutputArray(NULL)
+{
+}
+
+//------------------------------------------------------------------------------
+
+VertexSorter::~VertexSorter()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void VertexSorter::Init()
+{
+  fInputArray = ImportArray(GetString("InputArray", "VertexFinder/vertices"));
+
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "VertexFinder/tracks"));
+  fItTrackInputArray = fTrackInputArray->MakeIterator();
+
+  if (string (GetString("JetInputArray", "")) != "")
+    {
+      fJetInputArray = ImportArray(GetString("JetInputArray", ""));
+      fItJetInputArray = fJetInputArray->MakeIterator();
+    }
+
+  // import beamspot
+  try
+  {
+    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+  }
+  catch(runtime_error &e)
+  {
+    fBeamSpotInputArray = 0;
+  }  
+ 
+  fOutputArray = ExportArray(GetString("OutputArray", "clusters"));
+
+  fMethod = GetString ("Method", "BTV");
+}
+
+//------------------------------------------------------------------------------
+
+void VertexSorter::Finish()
+{
+  if(fItTrackInputArray) delete fItTrackInputArray;
+  if(fItJetInputArray) delete fItJetInputArray;
+}
+
+//------------------------------------------------------------------------------
+//
+Bool_t VertexSorter::secondDescending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
+{
+  return (pair0.second > pair1.second);
+}
+Bool_t VertexSorter::secondAscending (pair<UInt_t, Double_t> pair0, pair<UInt_t, Double_t> pair1)
+{
+  return (pair0.second < pair1.second);
+}
+
+void VertexSorter::Process()
+{
+  Candidate *candidate, *jetCandidate, *beamSpotCandidate;
+  unordered_map<Int_t, UInt_t> clusterIDToIndex;
+  unordered_map<Int_t, Double_t> clusterIDToSumPT2;
+  vector<pair<Int_t, Double_t> > sortedClusterIDs;
+
+  for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
+    {
+      const Candidate &cluster = *((Candidate *) fInputArray->At (iCluster));
+      clusterIDToIndex[cluster.ClusterIndex] = iCluster;
+      clusterIDToSumPT2[cluster.ClusterIndex] = 0.0;
+    }
+
+  if (fMethod == "BTV")
+    {
+      if (!fJetInputArray)
+        {
+          cout << "BTV PV sorting selected, but no jet collection given!" << endl;
+          throw 0;
+        }
+
+      fItTrackInputArray->Reset();
+      while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next())))
+        {
+          if (candidate->Momentum.Pt () < 1.0)
+            continue;
+          if (candidate->ClusterIndex < 0)
+            continue;
+          TLorentzVector p (candidate->Momentum.Px (), candidate->Momentum.Py (), candidate->Momentum.Pz (), candidate->Momentum.E ());
+          Bool_t isInJet = false;
+
+          fItJetInputArray->Reset();
+          while((jetCandidate = static_cast<Candidate*>(fItJetInputArray->Next())))
+            {
+              if (jetCandidate->Momentum.Pt () < 30.0)
+                continue;
+              TLorentzVector q (jetCandidate->Momentum.Px (), jetCandidate->Momentum.Py (), jetCandidate->Momentum.Pz (), jetCandidate->Momentum.E ());
+
+              if (p.DeltaR (q) > 0.4)
+                continue;
+              isInJet = true;
+              break;
+            }
+          if (!isInJet)
+            continue;
+
+          clusterIDToSumPT2.at (candidate->ClusterIndex) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
+        }
+
+      for (const auto &clusterID : clusterIDToSumPT2)
+        sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
+      sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
+    }
+  else if (fMethod == "GenClosest")
+    {
+      if (!fBeamSpotInputArray)
+        {
+          cout << "GenClosest PV sorting selected, but no beamspot collection given!" << endl;
+          throw 0;
+        }
+      if (!fBeamSpotInputArray->GetSize ())
+        {
+          cout << "Beamspot collection is empty!" << endl;
+          throw 0;
+        }
+
+      beamSpotCandidate = (Candidate *) fBeamSpotInputArray->At (0);
+      for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
+        {
+          const Candidate &cluster = *((Candidate *) fInputArray->At (iCluster));
+          sortedClusterIDs.push_back (make_pair (cluster.ClusterIndex, fabs (cluster.Position.Z () - beamSpotCandidate->Position.Z ())));
+        }
+      sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondAscending);
+    }
+  else if (fMethod == "GenBest")
+    {
+      fItTrackInputArray->Reset();
+      while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next())))
+        {
+          if (candidate->IsPU)
+            continue;
+          for (const auto &clusterID : clusterIDToIndex)
+            {
+              if (candidate->ClusterIndex != clusterID.first)
+                continue;
+              clusterIDToSumPT2.at (clusterID.first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
+            }
+        }
+
+      for (const auto &clusterID : clusterIDToSumPT2)
+        sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
+      sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
+    }
+  else
+    {
+      cout << "\"" << fMethod << "\" is not a valid sorting method!" << endl;
+      cout << "Valid methods are:" << endl;
+      cout << "  BTV" << endl;
+      cout << "  GenClosest" << endl;
+      cout << "  GenBest" << endl;
+      throw 0;
+    }
+
+  for (const auto &clusterID : sortedClusterIDs)
+    {
+      Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (clusterID.first));
+      if (fMethod == "BTV")
+        cluster->BTVSumPT2 = clusterID.second;
+      else if (fMethod == "GenClosest")
+        cluster->GenDeltaZ = clusterID.second;
+      else if (fMethod == "GenBest")
+        cluster->GenSumPT2 = clusterID.second;
+      fOutputArray->Add (cluster);
+    }
+}
+
+//------------------------------------------------------------------------------
Index: modules/VertexSorter.h
===================================================================
--- modules/VertexSorter.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
+++ modules/VertexSorter.h	(revision bd0a1aaab80e7e2982ee11cf10e51500c72a75fa)
@@ -0,0 +1,65 @@
+#ifndef VertexSorter_h
+#define VertexSorter_h
+
+/** \class VertexSorter
+ *
+ *
+ *  Sorts vertices according to different criteria
+ *
+ *  \authors A. Hart, M. Selvaggi
+ *
+ *
+*/
+
+#include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+#include <utility>
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <vector>
+#include <string>
+#include <unordered_map>
+
+
+using namespace std;
+
+class TObjArray;
+class Candidate;
+class TVector3;
+
+class VertexSorter: public DelphesModule
+{
+public:
+
+  VertexSorter();
+  ~VertexSorter();
+
+  void Init();
+  void Process();
+  void Finish();
+
+  static Bool_t secondDescending (pair<UInt_t, Double_t>, pair<UInt_t, Double_t>);
+  static Bool_t secondAscending (pair<UInt_t, Double_t>, pair<UInt_t, Double_t>);
+
+private:
+
+  TObjArray *fInputArray;
+
+  TObjArray *fTrackInputArray;
+  TIterator *fItTrackInputArray;
+
+  TObjArray *fJetInputArray;
+  TIterator *fItJetInputArray;
+
+  TObjArray *fBeamSpotInputArray;
+  TIterator *fItBeamSpotInputArray;
+
+  TObjArray *fOutputArray;
+
+  string fMethod;
+
+  ClassDef(VertexSorter, 1)
+};
+
+#endif
