Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/ModulesLinkDef.h	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -56,5 +56,4 @@
 #include "modules/PileUpMerger.h"
 #include "modules/JetPileUpSubtractor.h"
-#include "modules/TrackPileUpSubtractor.h"
 #include "modules/TrackTimingPileUpSubtractor.h"
 #include "modules/TaggingParticlesSkimmer.h"
@@ -114,5 +113,4 @@
 #pragma link C++ class PileUpMerger+;
 #pragma link C++ class JetPileUpSubtractor+;
-#pragma link C++ class TrackPileUpSubtractor+;
 #pragma link C++ class TrackTimingPileUpSubtractor+;
 #pragma link C++ class TaggingParticlesSkimmer+;
Index: modules/TimeSmearing.cc
===================================================================
--- modules/TimeSmearing.cc	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/TimeSmearing.cc	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -53,7 +53,6 @@
 
 TimeSmearing::TimeSmearing() :
-  fFormula(0), fItInputArray(0)
+  fItInputArray(0)
 {
-  fFormula = new DelphesFormula;
 }
 
@@ -61,6 +60,5 @@
 
 TimeSmearing::~TimeSmearing()
-{
-  if(fFormula) delete fFormula;
+{  
 }
 
@@ -71,5 +69,5 @@
   // read resolution formula
 
-  fFormula->Compile(GetString("TimeResolution", "1.0"));
+  fTimeResolution = GetDouble("TimeResolution", 1.);
 
   // import input array
@@ -77,5 +75,4 @@
   fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
   fItInputArray = fInputArray->MakeIterator();
-
   // create output array
 
@@ -95,9 +92,12 @@
 {
   Candidate *candidate, *mother;
-  Double_t ti, tf_smeared, tf, timeResolution;
-  Double_t pt, eta, phi, e, d0, dz, ctgTheta;
+  Double_t ti, tf_smeared, tf;
+  Double_t pt, eta, phi, e, p, l;
+  Double_t sigma_t, beta_particle;
 
 
   const Double_t c_light = 2.99792458E8;
+
+  cout << " STARTINNNGG ---------->" << endl;
 
   fItInputArray->Reset();
@@ -114,28 +114,29 @@
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    d0 = candidate->D0;
-    dz = candidate->DZ;
-    ctgTheta = candidate->CtgTheta;
+    p = candidateMomentum.P();
+    beta_particle = p/e;
+    l = candidate->L;
+    // apply smearing formula
 
-    // apply smearing formula
-    
-    timeResolution = fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta);
-    if(fabs(candidate->Position.Eta())<fEtaMax)
+    if(candidate->Charge != 0)
     { 
-      tf_smeared = tf + timeResolution*gRandom->Gaus(0, 1);
+      tf_smeared = tf + fTimeResolution*gRandom->Gaus(0, 1);
+      //mother = candidate;
+      //candidate = static_cast<Candidate*>(candidate->Clone());  // I am not sure that we need these lines !!!
+      //candidate->AddCandidate(mother);
+      candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);
+      candidate->Position.SetT(tf_smeared*1.0E3*c_light);
+      candidate->ErrorT = fTimeResolution*1.0E3*c_light;
+      fOutputArray->Add(candidate);
     }
-    else continue;
-
-    // double beta_particle = candidate->Momentum.P()/candidate->Momentum.E();
-    // ti = tf_smeared - candidate->Ld*1.0E-3/(c_light*beta_particle);
-
-    mother = candidate;
-    candidate = static_cast<Candidate*>(candidate->Clone());
-    candidate->AddCandidate(mother);
-    candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);
-    candidate->Position.SetT(tf_smeared*1.0E3*c_light);
-    candidate->ErrorT = timeResolution*1.0E3*c_light;
-
-    fOutputArray->Add(candidate);
+    else
+    {
+    	sigma_t = sqrt(pow(20,2) + pow(150,2)/e);
+    	ti = sigma_t - l*1.0E3/(c_light*beta_particle);   
+    	candidate->InitialPosition.SetT(ti);
+    	candidate->ErrorT = sigma_t*1.0E3*c_light;		// Do we need to sum with 100 like in upside ?
+    	fOutputArray->Add(candidate);
+    }
+   
   }
 }
Index: modules/TimeSmearing.h
===================================================================
--- modules/TimeSmearing.h	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/TimeSmearing.h	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -45,5 +45,5 @@
 
 private:
-  DelphesFormula *fFormula; //!
+  Double_t fTimeResolution; //!
   Double_t fEtaMax;
 
Index: modules/TimeSmearingNeutral.cc
===================================================================
--- modules/TimeSmearingNeutral.cc	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
+++ modules/TimeSmearingNeutral.cc	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -0,0 +1,130 @@
+/*
+ *  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 TimeSmearingNeutral
+ *
+ *  Performs transverse momentum resolution smearing.
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/TimeSmearingNeutral.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;
+
+//------------------------------------------------------------------------------
+
+TimeSmearingNeutral::TimeSmearingNeutral() :
+  fFormula(0), fItInputArray(0)
+{
+  fFormula = new DelphesFormula;
+}
+
+//------------------------------------------------------------------------------
+
+TimeSmearingNeutral::~TimeSmearingNeutral()
+{
+  if(fFormula) delete fFormula;
+}
+
+//------------------------------------------------------------------------------
+
+void TimeSmearingNeutral::Init()
+{
+  // read resolution formula
+
+  //fFormula->Compile(GetString("TimeResolution", "1.0"));
+  fTimeResolution = GetDouble("TimeResolution", 1.);
+  // import input array
+  fEtaMax = GetDouble("EtaMax", 6.);
+  fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array
+
+  fOutputArray = ExportArray(GetString("OutputArray", "muons"));
+}
+
+//------------------------------------------------------------------------------
+
+void TimeSmearingNeutral::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TimeSmearingNeutral::Process()
+{
+  Candidate *candidate, *mother;
+  Double_t ti, tf_smeared, tf, timeResolution;
+  Double_t pt, eta, phi, e, p, d0, dz, ctgTheta;
+  Double_t beta_particle;
+  Double_t l;
+  Double_t sigma_t;
+
+  const Double_t c_light = 2.99792458E8;
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+  {
+    ti = candidate->InitialPosition.T()*1.0E-3/c_light;
+    tf = candidate->Position.T()*1.0E-3/c_light;
+    e = candidate->Momentum.E();
+    p = candidate->Momentum.P();
+    beta_particle = p/e;
+    l = candidate->L;
+
+    if(candidate->Charge == 0)
+    {
+      sigma_t = sqrt(pow(20,2) + pow(150,2)/e)*1.0E-12;
+    }else
+      sigma_t = fTimeResolution*gRandom->Gaus(0, 1);;
+
+    ti = sigma_t - l*1.0E3/(c_light*beta_particle);
+
+
+    candidate->InitialPosition.SetT(ti);
+    candidate->ErrorT = sigma_t*1.0E3*c_light;
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/TimeSmearingNeutral.h
===================================================================
--- modules/TimeSmearingNeutral.h	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
+++ modules/TimeSmearingNeutral.h	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -0,0 +1,60 @@
+/*
+ *  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 TimeSmearingNeutral_h
+#define TimeSmearingNeutral_h
+
+/** \class TimeSmearingNeutral
+ *
+ *  Performs transverse time smearing.
+ *
+ *  \author Michele Selvaggi - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+class DelphesFormula;
+
+class TimeSmearingNeutral: public DelphesModule
+{
+public:
+  TimeSmearingNeutral();
+  ~TimeSmearingNeutral();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+  DelphesFormula *fFormula; //!
+  Double_t fEtaMax;
+  Double_t fTimeResolution;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(TimeSmearingNeutral, 1)
+};
+
+#endif
Index: modules/TrackTimingPileUpSubtractor.cc
===================================================================
--- modules/TrackTimingPileUpSubtractor.cc	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/TrackTimingPileUpSubtractor.cc	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -75,5 +75,6 @@
 
   // read resolution formula in m
-  fFormula->Compile(GetString("ZVertexResolution", "0.001"));
+  fChargedMinSignificance = GetDouble("ChargedMinSignificance", 3);
+  fNeutralMinSignificance = GetDouble("NeutralMinSignificance", 3);
 
   fPTMin = GetDouble("PTMin", 0.);
@@ -128,5 +129,5 @@
   Double_t tempPTSquare = 0;
   Double_t pt, eta, phi, e;
-  Double_t distance = 0;
+  Double_t distanceCharged, distanceNeutral = 0;
 
   // find z position of primary vertex
@@ -143,5 +144,4 @@
       tvtx = candidate->Position.T();
       tvtx_err = candidate->PositionError.T();
-      cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;
     } 
   }
@@ -167,18 +167,18 @@
       z = particle->Position.Z();
       z_err = particle->PositionError.Z();
-      t = particle->Position.T();
+      t = particle->InitialPosition.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
+      distanceCharged = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
+      distanceNeutral = pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
 
-      if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)
+      if(candidate->Charge != 0 && distanceCharged < fChargedMinSignificance)
       {
         candidate->IsRecoPU = 1;
       }
+      else if(candidate->Charge == 0 && distanceNeutral < fNeutralMinSignificance)
+      {
+        candidate->IsRecoPU = 1;
+      }  
       else
       {
Index: modules/TrackTimingPileUpSubtractor.h
===================================================================
--- modules/TrackTimingPileUpSubtractor.h	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/TrackTimingPileUpSubtractor.h	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -49,4 +49,7 @@
   DelphesFormula *fFormula; //!
 
+  Double_t fChargedMinSignificance;
+  Double_t fNeutralMinSignificance;
+
   Double_t fPTMin;
 
Index: modules/VertexFinderDA4D.cc
===================================================================
--- modules/VertexFinderDA4D.cc	(revision b4a3c552b76f4dfc936284ef81738bad3f98ac27)
+++ modules/VertexFinderDA4D.cc	(revision 9e1e2eb9d384d56155304b81cb07a99f1c2a020b)
@@ -425,4 +425,5 @@
     candidate->ClusterIndex = k;
     candidate->Position.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);
+    candidate->InitialPosition.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);    
     candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light);
     candidate->SumPT2 = 0;
