Index: modules/HighMassVertexRecover.cc
===================================================================
--- modules/HighMassVertexRecover.cc	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
+++ modules/HighMassVertexRecover.cc	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -0,0 +1,267 @@
+/** \class HighMassVertexRecover
+ *
+ *  Try to assign a vertex also to tracks which have
+ *  not been clusterized changing the mass hypotesis.
+ *
+ *  \author Olmo Cerri, Caltech
+ *
+ */
+
+#include "modules/HighMassVertexRecover.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 "TObjArray.h"
+#include "TLorentzVector.h"
+
+#include <iostream>
+#include <sstream>
+
+
+using namespace std;
+namespace vtx_DAZT
+{
+  static const Double_t c_light = 2.99792458e+8; // [m/s]
+}
+using namespace vtx_DAZT;
+
+//------------------------------------------------------------------------------
+
+HighMassVertexRecover::HighMassVertexRecover()
+{
+}
+
+//------------------------------------------------------------------------------
+
+HighMassVertexRecover::~HighMassVertexRecover()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void HighMassVertexRecover::Init()
+{
+  // Get parameters from the card
+  fVerbose         = GetInt("Verbose", 0);
+
+  fSigmaCompatibility = GetDouble("SigmaCompatibility", 2.0);
+
+  ExRootConfParam param = GetParam("MassHypot");
+  UInt_t size = param.GetSize();
+  if(size > 0)
+  {
+    for(UInt_t i = 0; i < size; ++i)
+    {
+      Double_t mass = GetDouble(param[i].GetString(), 1.);
+      fMassList.push_back(mass);
+    }
+  }
+  else
+  {
+    // Ordered by production probability in the SM
+    fMassList.push_back(0.13957);   //pion (+/-)
+    fMassList.push_back(0.49367);   // k (+/-)
+    fMassList.push_back(0.93827);   // proton
+    fMassList.push_back(0.00051);   //electron
+    fMassList.push_back(0.10558);   //muon
+  }
+
+  if(fVerbose)
+  {
+    for(unsigned int i = 0; i < fMassList.size(); i++)
+    {
+      cout << "mass: " << fMassList[i] << endl;
+    }
+  }
+
+
+  // import input array
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "VertexFinderDAClusterizerZT/tracks"));
+  fItTrackInputArray = fTrackInputArray->MakeIterator();
+
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "VertexFinderDAClusterizerZT/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+
+  // create output array
+  fTrackOutputArray = ExportArray(GetString("TrackOutputArray", "tracks"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+}
+
+//------------------------------------------------------------------------------
+
+void HighMassVertexRecover::Finish()
+{
+  delete fItTrackInputArray;
+  delete fItVertexInputArray;
+  delete fItVertexOutputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void HighMassVertexRecover::Process()
+{
+  if(fVerbose > 1)
+  {
+    cout << endl << endl << endl << "------------------ EVENT ----------------------------------------------" << endl;
+  }
+  // Make a copy of input VerticesfItVertexInputArray->Reset();
+  Candidate * vertex;
+  fItVertexInputArray->Reset();
+  while((vertex = static_cast<Candidate*>(fItVertexInputArray->Next())))
+  {
+    fVertexOutputArray->Add(static_cast<Candidate*>(vertex->Clone()));
+  }
+  fItVertexOutputArray = fVertexOutputArray->MakeIterator();
+
+  fItTrackInputArray->Reset();
+  Candidate * track;
+  while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+  {
+    Candidate* mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    if(track->ClusterIndex == -1)
+    {
+      if(fVerbose > 1)
+      {
+        cout << endl << "PID: " << track->PID << endl;
+        cout << "Pt: " << track->Momentum.Pt() << endl;
+        cout << Form("Zd: %.3f +/- %.3f", track->Zd, track->ErrorDZ)  << endl;
+      }
+      vector<Candidate*> vtx;                       //vector on vertices compatible in z
+      // Fill it with onyl the vertexes compatible
+      fItVertexOutputArray->Reset();
+      Candidate * vertex;
+      while((vertex = static_cast<Candidate*>(fItVertexOutputArray->Next())))
+      {
+        double dz = vertex->Position.Z() - track->Zd;
+        dz /= TMath::Hypot(track->ErrorDZ, vertex->PositionError.Z());
+
+        Double_t sigd0 = fabs(track->D0)/track->ErrorD0;
+        if(fabs(dz) < fSigmaCompatibility && sigd0 < 3) //Should be done better...is assuming v_xy = 0 and the stat threatment is poor/wrong
+        {
+          UInt_t i = 0;
+          UInt_t nv = vtx.size();
+          while(i<nv)
+          {
+            if(vtx[i]->SumPT2 > vertex->SumPT2) i++;
+            else break;
+          }
+          vtx.insert(vtx.begin()+i, vertex);
+        }
+      }
+
+      if(fVerbose > 2)
+      {
+        cout << "Compatible vertexes: " << vtx.size() << endl;
+        for(unsigned int i = 0; i<vtx.size(); i++)
+        {
+          cout << Form("SumPT2: %.0f GeV^2      Z: %.3f +/- %.3f mm      T: %.3f +/- %.3f ps", vtx[i]->SumPT2, vtx[i]->Position.Z(), vtx[i]->PositionError.Z(), vtx[i]->Position.T(), vtx[i]->PositionError.T()) << endl;
+        }
+      }
+
+      if(vtx.size() > 0)
+      {
+        // Try at first to see if the mass hypothesys can work
+        // The vertexes are order by SumPT2 and the mass hypotesys for relevance
+        UInt_t i = 0, match = 0;
+        while(i<vtx.size() && match == 0)
+        {
+          UInt_t j = 0;
+          while(j<fMassList.size() && match == 0)
+          {
+            auto Tpair = ComputeCATime(track, fMassList[j]);
+            double dt = vtx[i]->Position.T()*1.E9/c_light - Tpair.first;  // [ps]
+            // dt /= Tpair.second;
+            dt /= TMath::Hypot(Tpair.second, vtx[i]->PositionError.T()*1.E9/c_light);
+
+            if(fabs(dt)<fSigmaCompatibility)
+            {
+              match = 1;
+              track->Mass = fMassList[j];
+              track->ClusterIndex = vtx[i]->ClusterIndex;
+              track->InitialPosition.SetT(vtx[i]->Position.T());
+              track->InitialPosition.SetZ(vtx[i]->Position.Z());
+              track->Td = Tpair.first * 1E-9 * c_light;
+
+              vtx[i]->SumPT2 += track->Momentum.Pt()*track->Momentum.Pt();
+              vtx[i]->SumPt += track->Momentum.Pt();
+              if(fVerbose>5)
+              {
+                cout << "SM VTX found at " << i << ". dt = " << dt << "sigma" << endl;
+                cout << "vtx=> SumPT2 " << vtx[i]->SumPT2 << " - DOF: " << vtx[i]->ClusterNDF << endl;
+                cout << "Mass: " << fMassList[j] << endl;
+              }
+            }
+
+            j++;
+          }
+
+          i++;
+        }
+
+        if(match == 0)
+        {
+          Double_t p = track->Momentum.Pt() * sqrt(1 + track->CtgTheta*track->CtgTheta);
+          // Only Z
+          // Double_t beta_z = vtx[0]->Position.Z() - track->Position.Z();
+          // beta_z /= vtx[0]->Position.T() - track->Position.T();
+          // Double_t e = track->Momentum.Pt() * track->CtgTheta / beta_z;
+
+          //Full path length
+          Double_t beta = track->L / (vtx[0]->Position.T() - track->Position.T());
+          Double_t e = p / beta;
+
+          track->Mass = sqrt(e*e - p*p);
+          track->ClusterIndex = vtx[0]->ClusterIndex;
+          track->InitialPosition.SetT(vtx[0]->Position.T());
+          track->InitialPosition.SetZ(vtx[0]->Position.Z());
+          track->Td = vtx[0]->Position.T();
+
+          vtx[0]->SumPT2 += track->Momentum.Pt()*track->Momentum.Pt();
+          vtx[0]->SumPt += track->Momentum.Pt();
+          if(fVerbose>5)
+          {
+            cout << "BSM VTX fitted:" << endl;
+            cout << "vtx=> SumPT2 " << vtx[0]->SumPT2 << " - DOF: " << vtx[0]->ClusterNDF << endl;
+            cout << "Mass: " << track->Mass << endl;
+          }
+        }
+
+      }
+    }
+
+    fTrackOutputArray->Add(track);
+  }
+}
+
+//------------------------------------------------------------------------------
+// Auxiliary function to compute the CA time and erro given the mass
+pair<Double_t, Double_t> HighMassVertexRecover::ComputeCATime(Candidate * tk, Double_t m)
+{
+  Double_t p = tk->Momentum.Pt() * sqrt(1 + tk->CtgTheta*tk->CtgTheta);
+  Double_t e = sqrt(p*p + m*m);
+
+  Double_t t = tk->Position.T()*1.E9/c_light; // from [mm] to [ps]
+  //Full path length
+  t -= tk->L*1E9/(c_light*p/e);
+
+  // Only Z
+  // Double_t bz = tk->Momentum.Pt() * tk->CtgTheta/e;
+  // t += (tk->Zd - tk->Position.Z())*1E9/(c_light*bz);
+
+  pair<Double_t, Double_t> out;
+  out.first = t;
+  out.second = tk->ErrorT*1.E9/c_light; // [ps]
+
+  return out;
+}
Index: modules/HighMassVertexRecover.h
===================================================================
--- modules/HighMassVertexRecover.h	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
+++ modules/HighMassVertexRecover.h	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -0,0 +1,68 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef HighMassVertexRecover_h
+#define HighMassVertexRecover_h
+
+/** \class HighMassVertexRecover
+ *
+ *  Performs transverse time smearing.
+ *
+ *  \author Olmo Cerri, Caltech
+ *
+ */
+
+ #include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+
+class TIterator;
+class TObjArray;
+
+class HighMassVertexRecover: public DelphesModule
+{
+public:
+
+  HighMassVertexRecover();
+  ~HighMassVertexRecover();
+
+  void Init();
+  void Process();
+  void Finish();
+
+  std::pair<Double_t, Double_t> ComputeCATime(Candidate * tk, Double_t m);
+
+private:
+
+  UInt_t fVerbose;
+  Double_t fSigmaCompatibility;
+  std::vector<Double_t> fMassList;
+
+  TIterator *fItTrackInputArray;
+  const TObjArray *fTrackInputArray;
+
+  TIterator *fItVertexInputArray;
+  const TObjArray *fVertexInputArray;
+
+  TObjArray *fTrackOutputArray;
+  TObjArray *fVertexOutputArray;
+  TIterator *fItVertexOutputArray;
+
+  ClassDef(HighMassVertexRecover, 1)
+};
+
+#endif
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision 2b5ff2cc431ea2bfb8117eeada677e4d16b9da58)
+++ modules/ModulesLinkDef.h	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -56,4 +56,5 @@
 #include "modules/JetPileUpSubtractor.h"
 #include "modules/TrackPileUpSubtractor.h"
+#include "modules/TrackTimingPileUpSubtractor.h"
 #include "modules/TaggingParticlesSkimmer.h"
 #include "modules/PileUpJetID.h"
@@ -72,4 +73,5 @@
 #include "modules/VertexFinder.h"
 #include "modules/VertexFinderDA4D.h"
+#include "modules/HighMassVertexRecover.h"
 #include "modules/DecayFilter.h"
 #include "modules/ExampleModule.h"
@@ -111,4 +113,5 @@
 #pragma link C++ class JetPileUpSubtractor+;
 #pragma link C++ class TrackPileUpSubtractor+;
+#pragma link C++ class TrackTimingPileUpSubtractor+;
 #pragma link C++ class TaggingParticlesSkimmer+;
 #pragma link C++ class PileUpJetID+;
@@ -127,4 +130,5 @@
 #pragma link C++ class VertexFinder+;
 #pragma link C++ class VertexFinderDA4D+;
+#pragma link C++ class HighMassVertexRecover+;
 #pragma link C++ class DecayFilter+;
 #pragma link C++ class ExampleModule+;
Index: modules/TrackPileUpSubtractor.cc
===================================================================
--- modules/TrackPileUpSubtractor.cc	(revision 2b5ff2cc431ea2bfb8117eeada677e4d16b9da58)
+++ modules/TrackPileUpSubtractor.cc	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -117,5 +117,5 @@
 void TrackPileUpSubtractor::Process()
 {
-  Candidate *candidate, *particle;
+  Candidate *candidate, *particle, *particleTest;
   map<TIterator *, TObjArray *>::iterator itInputMap;
   TIterator *iterator;
@@ -123,15 +123,37 @@
   Double_t z, zvtx = 0;
   Double_t pt, eta, phi, e;
-
+  Double_t sumPT2 = 0;
+  Double_t sumSquare = 0;
+  int counter = 0;
+  Double_t tempPTSquare = 0;
+  Double_t tempZVertex = 0; 
   // find z position of primary vertex
 
+  cout << " ---------- NEW EVENT --------- " << endl; 
+
   fItVertexInputArray->Reset();
+  cout << " NUMBER OF VERTICES : " << fVertexInputArray->GetEntriesFast() << endl;
   while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next())))
   {
-    if(!candidate->IsPU)
-    {
-      zvtx = candidate->Position.Z();
-      // break;
-    }
+      cout << " ---------- NEW VERTEX --------- " << candidate->IsPU << endl; 
+      tempZVertex = candidate->Position.Z();
+      tempPTSquare = candidate->SumPT2;
+      if(tempPTSquare > sumSquare)
+      {
+        sumSquare = tempPTSquare;
+        zvtx = tempZVertex;
+        cout << " Sum Square : " << sumSquare << " Z of Vertex : " << zvtx << "  Is PileUp : " << candidate->IsPU << endl;
+      } 
+      
+      /*
+      for(int i = 0; i < candidate->GetCandidates()->GetEntriesFast(); i++)
+      {
+        particleTest = static_cast<Candidate *>(candidate->GetCandidates()->At(i));
+        TLorentzVector &candidateMomentumNew = particleTest->Momentum;
+        if(candidateMomentumNew.Pt() > 1.0)
+          sumSquare += (candidateMomentumNew.Pt()*candidateMomentumNew.Pt());
+        cout << candidateMomentumNew.Pt() << "  " << candidate->SumPT2 << " counter : " << candidate->GetCandidates()->GetEntriesFast() << endl;
+      }  
+      */
   }
 
@@ -147,5 +169,5 @@
     {
       particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
-      const TLorentzVector &candidateMomentum = particle->Momentum;
+      const TLorentzVector &candidateMomentum = candidate->Momentum;
 
       eta = candidateMomentum.Eta();
@@ -154,12 +176,20 @@
       e = candidateMomentum.E();
 
-      z = particle->Position.Z();
+      z = candidate->Position.Z();
+      counter ++; 
+      //sum = pt*pt;
+      if(pt > 1.0)
+        sumPT2 += pt*pt;
 
       // apply pile-up subtraction
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
-
-      if(candidate->Charge != 0 && candidate->IsPU && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
+      // cout << particle->IsRecoPU << " ParticleSumSquare : " << sumPT2 << " VertexResult : " << sumSquare << " Added SumSquare : " << zvtx << endl;
+      //cout << pt << "  " << candidate->IsPU << "  " <<  TMath::Abs(z - zvtx) << "   " << sumPT2 << "   " << sumPT2 << endl;
+
+      
+      if(particle->Charge != 0 && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
       {
         candidate->IsRecoPU = 1;
+        //cout <<  TMath::Abs(z - zvtx) << "   " << fFormula->Eval(pt, eta, phi, e) * 1.0e3 << endl;
       }
       else
Index: modules/TrackTimingPileUpSubtractor.cc
===================================================================
--- modules/TrackTimingPileUpSubtractor.cc	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
+++ modules/TrackTimingPileUpSubtractor.cc	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -0,0 +1,190 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/** \class TrackPileUpSubtractor
+ *
+ *  Subtract pile-up contribution from tracks.
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/TrackTimingPileUpSubtractor.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootClassifier.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootResult.h"
+
+#include "TDatabasePDG.h"
+#include "TFormula.h"
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
+#include "TString.h"
+
+#include <algorithm>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+TrackTimingPileUpSubtractor::TrackTimingPileUpSubtractor() :
+  fFormula(0)
+{
+  fFormula = new DelphesFormula;
+}
+
+//------------------------------------------------------------------------------
+
+TrackTimingPileUpSubtractor::~TrackTimingPileUpSubtractor()
+{
+  if(fFormula) delete fFormula;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackTimingPileUpSubtractor::Init()
+{
+  // import input array
+
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+
+  // read resolution formula in m
+  fFormula->Compile(GetString("ZVertexResolution", "0.001"));
+
+  fPTMin = GetDouble("PTMin", 0.);
+
+  // import arrays with output from other modules
+
+  ExRootConfParam param = GetParam("InputArray");
+  Long_t i, size;
+  const TObjArray *array;
+  TIterator *iterator;
+
+  size = param.GetSize();
+  for(i = 0; i < size / 2; ++i)
+  {
+    array = ImportArray(param[i * 2].GetString());
+    iterator = array->MakeIterator();
+
+    fInputMap[iterator] = ExportArray(param[i * 2 + 1].GetString());
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void TrackTimingPileUpSubtractor::Finish()
+{
+  map<TIterator *, TObjArray *>::iterator itInputMap;
+  TIterator *iterator;
+
+  for(itInputMap = fInputMap.begin(); itInputMap != fInputMap.end(); ++itInputMap)
+  {
+    iterator = itInputMap->first;
+
+    if(iterator) delete iterator;
+  }
+
+  if(fItVertexInputArray) delete fItVertexInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackTimingPileUpSubtractor::Process()
+{
+  Candidate *candidate, *particle;
+  map<TIterator *, TObjArray *>::iterator itInputMap;
+  TIterator *iterator;
+  TObjArray *array;
+  Double_t z, zvtx = 0;
+  Double_t z_err, zvtx_err = 0;
+  Double_t t, tvtx = 0;
+  Double_t t_err, tvtx_err = 0;
+  Double_t sumPTSquare = 0;
+  Double_t tempPTSquare = 0;
+  Double_t pt, eta, phi, e;
+  Double_t distance = 0;
+
+  // find z position of primary vertex
+
+  fItVertexInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next())))
+  {
+    tempPTSquare = candidate->SumPT2;
+    if(tempPTSquare > sumPTSquare)
+    {
+      sumPTSquare = tempPTSquare;
+      zvtx = candidate->Position.Z();
+      zvtx_err = candidate->PositionError.Z();
+      tvtx = candidate->Position.T();
+      tvtx_err = candidate->PositionError.T();
+      cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;
+    } 
+  }
+
+  // loop over all input arrays
+  for(itInputMap = fInputMap.begin(); itInputMap != fInputMap.end(); ++itInputMap)
+  {
+    iterator = itInputMap->first;
+    array = itInputMap->second;
+
+    // loop over all candidates
+    iterator->Reset();
+    while((candidate = static_cast<Candidate *>(iterator->Next())))
+    {
+      particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+      const TLorentzVector &candidateMomentum = particle->Momentum;
+
+      eta = candidateMomentum.Eta();
+      pt = candidateMomentum.Pt();
+      phi = candidateMomentum.Phi();
+      e = candidateMomentum.E();
+
+      z = particle->Position.Z();
+      z_err = particle->PositionError.Z();
+      t = particle->Position.T();
+      t_err = particle->PositionError.T();
+
+      // apply pile-up subtraction
+      distance = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
+      //cout << " t : " << tvtx << "  t(particle)" << t << endl;
+      // here I calculated distance using Z and T of selected vertex (highest sum Pt square) and particles
+      // however z_err of vertices is gives 0 because of using CMS trackResolutionCMS.tcl (in that formula, there is limitation on |eta| < 2.5)
+      // thats why I used TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0
+
+      if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)
+      {
+        candidate->IsRecoPU = 1;
+      }
+      else
+      {
+        candidate->IsRecoPU = 0;
+        if(candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
+      }
+    }
+  }
+}
Index: modules/TrackTimingPileUpSubtractor.h
===================================================================
--- modules/TrackTimingPileUpSubtractor.h	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
+++ modules/TrackTimingPileUpSubtractor.h	(revision d4e5b6a29d976f833d389ea7af3115081d34d5d5)
@@ -0,0 +1,62 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef TrackTimingPileUpSubtractor_h
+#define TrackTimingPileUpSubtractor_h
+
+/** \class TrackPileUpSubtractor
+ *
+ *  Subtract pile-up contribution from tracks.
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+#include <map>
+
+class TIterator;
+class TObjArray;
+class DelphesFormula;
+
+class TrackTimingPileUpSubtractor: public DelphesModule
+{
+public:
+  TrackTimingPileUpSubtractor();
+  ~TrackTimingPileUpSubtractor();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+  DelphesFormula *fFormula; //!
+
+  Double_t fPTMin;
+
+  std::map<TIterator *, TObjArray *> fInputMap; //!
+
+  ClassDef(TrackTimingPileUpSubtractor, 1)
+
+    TIterator *fItVertexInputArray; //!
+
+  const TObjArray *fVertexInputArray; //!
+};
+
+#endif
