Index: Makefile
===================================================================
--- Makefile	(revision 95b4e9f53f115f744fd16d3ccbc8713b04dacea4)
+++ Makefile	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1669,6 +1669,5 @@
 
 modules/VertexFinderDA4D.h: \
-	classes/DelphesModule.h \
-	classes/DelphesClasses.h
+	classes/DelphesModule.h
 	@touch $@
 
@@ -1942,6 +1941,5 @@
 
 modules/VertexSorter.h: \
-	classes/DelphesModule.h \
-	classes/DelphesClasses.h
+	classes/DelphesModule.h
 	@touch $@
 
@@ -1951,6 +1949,5 @@
 
 modules/VertexFinder.h: \
-	classes/DelphesModule.h \
-	classes/DelphesClasses.h
+	classes/DelphesModule.h
 	@touch $@
 
Index: modules/VertexFinderDA4D.cc
===================================================================
--- modules/VertexFinderDA4D.cc	(revision 95b4e9f53f115f744fd16d3ccbc8713b04dacea4)
+++ modules/VertexFinderDA4D.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -42,11 +42,49 @@
 static const Double_t c_light   = 2.99792458e+8 * m/s;
 
-
-namespace {
-  bool recTrackLessZ1(const VertexFinderDA4D::track_t & tk1,
-                      const VertexFinderDA4D::track_t & tk2)
-  {
-    return tk1.z < tk2.z;
-  }
+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)
+  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;
+  double t;
+  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;
+};
+
+static bool split(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
+static double update1(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y);
+static double update2(double beta, std::vector<track_t> &tks, std::vector<vertex_t> &y, double &rho0, const double dzCutOff);
+static void dump(const double beta, const std::vector<vertex_t> & y, const std::vector<track_t> & tks);
+static bool merge(std::vector<vertex_t> &);
+static bool merge(std::vector<vertex_t> &, double &);
+static bool purge(std::vector<vertex_t> &, std::vector<track_t> & , double &, const double, const double);
+static void splitAll(std::vector<vertex_t> &y);
+static double beta0(const double betamax, std::vector<track_t> &tks, std::vector<vertex_t> &y, const double coolingFactor);
+static double Eik(const track_t &t, const vertex_t &k);
+
+static bool recTrackLessZ1(const track_t & tk1, const track_t & tk2)
+{
+  return tk1.z < tk2.z;
 }
 
@@ -214,5 +252,5 @@
 
      if (fVerbose){
-	   std::cout << "x,y,z";
+     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;
@@ -242,6 +280,5 @@
 //------------------------------------------------------------------------------
 
-
-void VertexFinderDA4D::clusterize(const TObjArray & tracks, TObjArray & clusters)
+void VertexFinderDA4D::clusterize(const TObjArray &tracks, TObjArray &clusters)
 {
   if(fVerbose) {
@@ -294,5 +331,4 @@
 }
 
-
 //------------------------------------------------------------------------------
 
@@ -303,5 +339,58 @@
   vector< Candidate* > clusters;
 
-  vector<track_t> tks = fill();
+  vector<track_t> tks;
+  track_t tr;
+  Double_t z, dz, t, l, dt, d0, d0error;
+
+  // 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);
+    }
+  }
 
   //print out input tracks
@@ -341,6 +430,6 @@
 
   // estimate first critical temperature
-  double beta=beta0(fBetaMax, tks, y);
-  niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
+  double beta=beta0(fBetaMax, tks, y, fCoolingFactor);
+  niter=0; while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
 
   // annealing loop, stop when T<Tmin  (i.e. beta>1/Tmin)
@@ -348,6 +437,6 @@
 
     if(fUseTc){
-      update(beta, tks,y);
-      while(merge(y,beta)){update(beta, tks,y);}
+      update1(beta, tks,y);
+      while(merge(y,beta)){update1(beta, tks,y);}
       split(beta, tks,y);
       beta=beta/fCoolingFactor;
@@ -358,5 +447,5 @@
 
    // make sure we are not too far from equilibrium before cooling further
-   niter=0; while((update(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
+   niter=0; while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){ }
 
   }
@@ -364,16 +453,16 @@
   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);}
+    update1(beta, tks,y);
+    while(merge(y,beta)){update1(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)){}
+      while((update1(beta, tks,y)>1.e-6)  && (niter++ < fMaxIterations)){}
       merge(y,beta);
-      update(beta, tks,y);
+      update1(beta, tks,y);
     }
   }else{
     // merge collapsed clusters
-    while(merge(y,beta)){update(beta, tks,y);}
+    while(merge(y,beta)){update1(beta, tks,y);}
     if(fVerbose ){ cout << "dump after 1st merging " << endl;  dump(beta,y,tks);}
   }
@@ -381,5 +470,5 @@
   // 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)){  }
+  niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-8)  && (niter++ < fMaxIterations)){  }
   if(fVerbose  ){ cout << "rho0=" << rho0 <<   " niter=" << niter <<  endl; dump(beta,y,tks);}
 
@@ -392,9 +481,9 @@
   // 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)){  }
+    while(purge(y,tks,rho0, beta, fDzCutOff)){
+      niter=0; while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
     }
     beta/=fCoolingFactor;
-    niter=0; while((update(beta, tks, y, rho0) > 1.e-6)  && (niter++ < fMaxIterations)){  }
+    niter=0; while((update2(beta, tks, y, rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
   }
 
@@ -402,5 +491,5 @@
 //   // 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)){  }
+//     niter=0; while((update2(beta, tks,y,rho0, fDzCutOff) > 1.e-6)  && (niter++ < fMaxIterations)){  }
 //   }
 
@@ -440,6 +529,6 @@
       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 ) ){
+  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;
@@ -482,80 +571,13 @@
 //------------------------------------------------------------------------------
 
-vector<VertexFinderDA4D::track_t>
-VertexFinderDA4D::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 VertexFinderDA4D::Eik(const track_t & t, const vertex_t &k ) const {
+static double Eik(const track_t & t, const vertex_t &k)
+{
   return std::pow(t.z-k.z,2.)/t.dz2 + std::pow(t.t - k.t,2.)/t.dt2;
 }
 
-//--------------------------------------------------------------------------------
-
-void VertexFinderDA4D::dump(const double beta, const vector<vertex_t> & y, const vector<track_t> & tks0)const{
-
+//------------------------------------------------------------------------------
+
+static void dump(const double beta, const vector<vertex_t> &y, const vector<track_t> &tks0)
+{
   // copy and sort for nicer printout
   vector<track_t> tks;
@@ -564,5 +586,5 @@
 
   cout << "-----DAClusterizerInZT::dump ----" << endl;
-  cout << " beta=" << beta << "   betamax= " << fBetaMax << endl;
+  cout << " beta=" << beta << endl;
   cout << "                                                               z= ";
   cout.precision(4);
@@ -587,48 +609,40 @@
   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 << "        ";
-	}
+  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 VertexFinderDA4D::update( double beta,
-                                  vector<track_t> & tks,
-                                  vector<vertex_t> & y ) const {
+}
+
+//------------------------------------------------------------------------------
+
+static double update1(double beta, vector<track_t> &tks, vector<vertex_t> &y)
+{
   //update weights and vertex positions
   // mass constrained annealing without noise
@@ -661,10 +675,10 @@
       // 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->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);
+  k->swE += w * Eik(tks[i],*k);
       }
     }else{
@@ -687,5 +701,5 @@
       k->Tc = 2.*k->swE/k->sw;
     }else{
-      if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
+      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl
       k->Tc=-1;
     }
@@ -700,8 +714,6 @@
 //------------------------------------------------------------------------------
 
-double VertexFinderDA4D::update( double beta,
-                                  vector<track_t> & tks,
-                                  vector<vertex_t> & y,
-                                  double & rho0 ) const {
+static double update2(double beta, vector<track_t> &tks, vector<vertex_t> &y, double &rho0, double dzCutOff)
+{
   // MVF style, no more vertex weights, update tracks weights and vertex positions, with noise
   // returns the squared sum of changes of vertex positions
@@ -720,5 +732,5 @@
 
     // update pik and Zi and Ti
-    double Zi = rho0*std::exp(-beta*(fDzCutOff*fDzCutOff));// cut-off (eventually add finite size in time)
+    double Zi = rho0*std::exp(-beta*(dzCutOff*dzCutOff));// 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++){
@@ -732,10 +744,10 @@
        // 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->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);
+  k->swE += w * Eik(tks[i],*k);
       }
     }
@@ -754,5 +766,5 @@
       k->Tc  = 2*k->swE/k->sw;
     }else{
-      if(fVerbose){cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;}
+      // cout << " a cluster melted away ?  pk=" << k->pk <<  " sumw=" << k->sw <<  endl;
       k->Tc = 0;
     }
@@ -766,5 +778,6 @@
 //------------------------------------------------------------------------------
 
-bool VertexFinderDA4D::merge(vector<vertex_t> & y)const{
+static bool merge(vector<vertex_t> &y)
+{
   // merge clusters that collapsed or never separated, return true if vertices were merged, false otherwise
 
@@ -794,5 +807,6 @@
 //------------------------------------------------------------------------------
 
-bool VertexFinderDA4D::merge(vector<vertex_t> & y, double & beta)const{
+static bool merge(vector<vertex_t> &y, double &beta)
+{
   // merge clusters that collapsed or never separated,
   // only merge if the estimated critical temperature of the merged vertex is below the current temperature
@@ -809,17 +823,17 @@
 
       if(Tc*beta<1){
-	if(rho>0){
-	  k->z = ( k->pk * k->z + (k+1)->z * (k+1)->pk)/rho;
+  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);
+  }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;
+  }
+  k->pk  = rho;
+  k->sw += (k+1)->sw;
+  k->swE = swE;
+  k->Tc  = Tc;
+  y.erase(k+1);
+  return true;
       }
     }
@@ -831,5 +845,6 @@
 //------------------------------------------------------------------------------
 
-bool VertexFinderDA4D::purge(vector<vertex_t> & y, vector<track_t> & tks, double & rho0, const double beta)const{
+static bool purge(vector<vertex_t> &y, vector<track_t> &tks, double & rho0, const double beta, const double dzCutOff)
+{
   // eliminate clusters with only one significant/unique track
   if(y.size()<2)  return false;
@@ -841,10 +856,10 @@
     int nUnique=0;
     double sump=0;
-    double pmax=k->pk/(k->pk+rho0*exp(-beta*fDzCutOff*fDzCutOff));
+    double pmax=k->pk/(k->pk+rho0*exp(-beta*dzCutOff*dzCutOff));
     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++; }
+  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++; }
       }
     }
@@ -857,6 +872,5 @@
 
   if(k0!=y.end()){
-    if(fVerbose){cout << "eliminating prototype at " << k0->z << "," << k0->t
-                      << " with sump=" << sumpmin << endl;}
+    //cout << "eliminating prototype at " << k0->z << "," << k0->t << " with sump=" << sumpmin << endl;
     //rho0+=k0->pk;
     y.erase(k0);
@@ -869,7 +883,6 @@
 //------------------------------------------------------------------------------
 
-double VertexFinderDA4D::beta0( double betamax,
-                                 vector<track_t> & tks,
-                                 vector<vertex_t> & y ) const {
+static double beta0(double betamax, vector<track_t> &tks, vector<vertex_t> &y, const double coolingFactor)
+{
 
   double T0=0;  // max Tc for beta=0
@@ -906,22 +919,21 @@
 
   if (T0>1./betamax){
-    return betamax/pow(fCoolingFactor, int(std::log(T0*betamax)/std::log(fCoolingFactor))-1 );
+    return betamax/pow(coolingFactor, int(std::log(T0*betamax)/std::log(coolingFactor))-1 );
   }else{
     // ensure at least one annealing step
-    return betamax/fCoolingFactor;
-  }
-}
-
-//------------------------------------------------------------------------------
-
-bool VertexFinderDA4D::split( double beta,
-                               vector<track_t> & tks,
-                               vector<vertex_t> & y) const {
+    return betamax/coolingFactor;
+  }
+}
+
+//------------------------------------------------------------------------------
+
+static bool split(double beta, vector<track_t> &tks, vector<vertex_t> &y)
+{
   // 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;
+  const double epsilon = 1e-3; // split all single vertices by 10 um
+  bool split = false;
 
   // avoid left-right biases by splitting highest Tc first
@@ -943,12 +955,12 @@
     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;
-	}
+  //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;
+  }
       }
     }
@@ -974,5 +986,5 @@
      // adjust remaining pointers
       for(unsigned int jc=ic; jc<critical.size(); jc++){
-	if (critical[jc].second>ik) {critical[jc].second++;}
+        if (critical[jc].second>ik) {critical[jc].second++;}
       }
     }
@@ -985,10 +997,11 @@
 //------------------------------------------------------------------------------
 
-void VertexFinderDA4D::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
+void splitAll(vector<vertex_t> &y)
+{
+
+
+  const double epsilon=1e-3;      // split all single vertices by 10 um
+  const double zsep=2*epsilon;    // split vertices that are isolated by at least zsep (vertices that haven't collapsed)
+  const double tsep=2*epsilon;    // check t as well
 
   vector<vertex_t> y1;
Index: modules/VertexFinderDA4D.h
===================================================================
--- modules/VertexFinderDA4D.h	(revision 95b4e9f53f115f744fd16d3ccbc8713b04dacea4)
+++ modules/VertexFinderDA4D.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -30,67 +30,6 @@
   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)
-    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;
-    double t;
-    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);
-
+  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:
Index: modules/VertexSorter.cc
===================================================================
--- modules/VertexSorter.cc	(revision 95b4e9f53f115f744fd16d3ccbc8713b04dacea4)
+++ modules/VertexSorter.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -115,6 +115,9 @@
   Candidate *candidate, *jetCandidate, *beamSpotCandidate;
   map<Int_t, UInt_t> clusterIDToIndex;
+  map<Int_t, UInt_t>::const_iterator itClusterIDToIndex;
   map<Int_t, Double_t> clusterIDToSumPT2;
+  map<Int_t, Double_t>::const_iterator itClusterIDToSumPT2;
   vector<pair<Int_t, Double_t> > sortedClusterIDs;
+  vector<pair<Int_t, Double_t> >::const_iterator itSortedClusterIDs;
 
   for (Int_t iCluster = 0; iCluster < fInputArray->GetEntries (); iCluster++)
@@ -161,6 +164,6 @@
         }
 
-      for (const auto &clusterID : clusterIDToSumPT2)
-        sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
+      for (itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
+        sortedClusterIDs.push_back (make_pair (itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
       sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
     }
@@ -193,14 +196,14 @@
           if (candidate->IsPU)
             continue;
-          for (const auto &clusterID : clusterIDToIndex)
+          for (itClusterIDToIndex = clusterIDToIndex.begin(); itClusterIDToIndex != clusterIDToIndex.end(); ++itClusterIDToIndex)
             {
-              if (candidate->ClusterIndex != clusterID.first)
+              if (candidate->ClusterIndex != itClusterIDToIndex->first)
                 continue;
-              clusterIDToSumPT2.at (clusterID.first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
+              clusterIDToSumPT2.at (itClusterIDToIndex->first) += candidate->Momentum.Pt () * candidate->Momentum.Pt ();
             }
         }
 
-      for (const auto &clusterID : clusterIDToSumPT2)
-        sortedClusterIDs.push_back (make_pair (clusterID.first, clusterID.second));
+      for (itClusterIDToSumPT2 = clusterIDToSumPT2.begin(); itClusterIDToSumPT2 != clusterIDToSumPT2.end(); ++itClusterIDToSumPT2)
+        sortedClusterIDs.push_back (make_pair (itClusterIDToSumPT2->first, itClusterIDToSumPT2->second));
       sort (sortedClusterIDs.begin (), sortedClusterIDs.end (), secondDescending);
     }
@@ -214,14 +217,13 @@
       throw 0;
     }
-
-  for (const auto &clusterID : sortedClusterIDs)
-    {
-      Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (clusterID.first));
+  for (itSortedClusterIDs = sortedClusterIDs.begin(); itSortedClusterIDs != sortedClusterIDs.end(); ++itSortedClusterIDs)
+    {
+      Candidate *cluster = (Candidate *) fInputArray->At (clusterIDToIndex.at (itSortedClusterIDs->first));
       if (fMethod == "BTV")
-        cluster->BTVSumPT2 = clusterID.second;
+        cluster->BTVSumPT2 = itSortedClusterIDs->second;
       else if (fMethod == "GenClosest")
-        cluster->GenDeltaZ = clusterID.second;
+        cluster->GenDeltaZ = itSortedClusterIDs->second;
       else if (fMethod == "GenBest")
-        cluster->GenSumPT2 = clusterID.second;
+        cluster->GenSumPT2 = itSortedClusterIDs->second;
       fOutputArray->Add (cluster);
     }
