Index: modules/Calorimeter.cc
===================================================================
--- modules/Calorimeter.cc	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/Calorimeter.cc	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -330,4 +330,13 @@
       fTrackHCalEnergy = 0.0;
 
+      fTowerECalTime = 0.0;
+      fTowerHCalTime = 0.0;
+
+      fTrackECalTime = 0.0;
+      fTrackHCalTime = 0.0;
+
+      fTowerECalWeightTime = 0.0; 
+      fTowerHCalWeightTime = 0.0;
+     
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -343,5 +352,7 @@
       track = static_cast<Candidate*>(fTrackInputArray->At(number));
       momentum = track->Momentum;
-
+      position = track->Position;
+
+      
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -349,4 +360,10 @@
       fTrackECalEnergy += ecalEnergy;
       fTrackHCalEnergy += hcalEnergy;
+      
+      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+       
+      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
+      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -360,4 +377,5 @@
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
+    position = particle->Position;
 
     // fill current tower
@@ -367,4 +385,11 @@
     fTowerECalEnergy += ecalEnergy;
     fTowerHCalEnergy += hcalEnergy;
+
+    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
+    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
+    
 
     fTower->AddCandidate(particle);
@@ -383,4 +408,5 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
+  Double_t ecalTime, hcalTime, time;
 
   if(!fTower) return;
@@ -392,4 +418,5 @@
 
   ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
+  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
 
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
@@ -399,6 +426,8 @@
 
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
+  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
 
   energy = ecalEnergy + hcalEnergy;
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 
 
 //  eta = fTowerEta;
@@ -410,5 +439,6 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
+ // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
+  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
@@ -420,4 +450,5 @@
   fTower->Edges[3] = fTowerEdges[3];
 
+
   // fill calorimeter towers and photon candidates
   if(energy > 0.0)
Index: modules/Calorimeter.h
===================================================================
--- modules/Calorimeter.h	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/Calorimeter.h	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -45,4 +45,11 @@
   Double_t fTowerECalEnergy, fTowerHCalEnergy;
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
+  
+  Double_t fTowerECalTime, fTowerHCalTime;
+  Double_t fTrackECalTime, fTrackHCalTime;
+   
+  Double_t fTowerECalWeightTime, fTowerHCalWeightTime;
+  Double_t fTrackECalWeightTime, fTrackHCalWeightTime;
+  
   Int_t fTowerTrackHits, fTowerPhotonHits;
 
Index: modules/FastJetFinder.cc
===================================================================
--- modules/FastJetFinder.cc	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/FastJetFinder.cc	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -193,4 +193,5 @@
   TLorentzVector momentum;
   Double_t deta, dphi, detaMax, dphiMax;
+  Double_t time, weightTime, avTime;
   Int_t number;
   Double_t rho = 0;
@@ -259,4 +260,7 @@
     candidate = factory->NewCandidate();
 
+    time=0;
+    weightTime=0;
+
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
@@ -269,9 +273,15 @@
       if(deta > detaMax) detaMax = deta;
       if(dphi > dphiMax) dphiMax = dphi;
-
+      
+      time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
+      weightTime += TMath::Sqrt(constituent->Momentum.E());
+    
       candidate->AddCandidate(constituent);
     }
+   
+    avTime = time/weightTime;
 
     candidate->Momentum = momentum;
+    candidate->Position.SetT(avTime);
     candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
 
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/ModulesLinkDef.h	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -19,4 +19,5 @@
 #include "modules/EnergySmearing.h"
 #include "modules/MomentumSmearing.h"
+#include "modules/TimeSmearing.h"
 #include "modules/Calorimeter.h"
 #include "modules/Isolation.h"
@@ -50,4 +51,5 @@
 #pragma link C++ class EnergySmearing+;
 #pragma link C++ class MomentumSmearing+;
+#pragma link C++ class TimeSmearing+;
 #pragma link C++ class Calorimeter+;
 #pragma link C++ class Isolation+;
Index: modules/PileUpMerger.cc
===================================================================
--- modules/PileUpMerger.cc	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/PileUpMerger.cc	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -16,5 +16,5 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
+#include "classes/DelphesTF2.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -41,7 +41,9 @@
 
 PileUpMerger::PileUpMerger() :
-  fReader(0), fItInputArray(0)
-{
-}
+  fFunction(0), fReader(0), fItInputArray(0)
+{
+  fFunction = new DelphesTF2;
+}
+
 
 //------------------------------------------------------------------------------
@@ -49,4 +51,5 @@
 PileUpMerger::~PileUpMerger()
 {
+  delete fFunction;
 }
 
@@ -60,6 +63,13 @@
 
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
-
+  
+  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
+  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
+
+  // read vertex smearing formula
+	
+  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
+  fFunction->SetRange(-fZVertexSpread,-fTVertexSpread,fZVertexSpread,fTVertexSpread);
+ 
   fileName = GetString("PileUpFile", "MinBias.pileup");
   fReader = new DelphesPileUpReader(fileName);
@@ -90,18 +100,36 @@
   Float_t x, y, z, t;
   Float_t px, py, pz, e;
-  Double_t dz, dphi;
+  Double_t dz, dphi, dt;
   Int_t numberOfEvents, event;
   Long64_t allEntries, entry;
-  Candidate *candidate;
+  Candidate *candidate, *vertexcandidate;
   DelphesFactory *factory;
+  
+  const Double_t c_light = 2.99792458E8;
 
   fItInputArray->Reset();
+  
+    // --- Deal with Primary vertex first  ------
+  
+  fFunction->GetRandom2(dz,dt);
+ 
+  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
+  
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
+    candidate->Position.SetXYZT(x, y, z+dz, t+dt);
     fParticleOutputArray->Add(candidate);
   }
 
   factory = GetFactory();
-
+  
+  vertexcandidate = factory->NewCandidate();
+  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+  fVertexOutputArray->Add(vertexcandidate);
+ 
+  
+  // --- Then with pile-up vertices  ------
+ 
   switch(fPileUpDistribution)
   {
@@ -118,5 +146,5 @@
 
   allEntries = fReader->GetEntries();
-
+ 
   for(event = 0; event < numberOfEvents; ++event)
   {
@@ -128,11 +156,19 @@
 
     fReader->ReadEntry(entry);
-
-    dz = gRandom->Gaus(0.0, fZVertexSpread);
+ 
+   // --- Pile-up vertex smearing 
+    
+    fFunction->GetRandom2(dz,dt);
+ 
+    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
+ 
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    candidate = factory->NewCandidate();
-    candidate->Position.SetXYZT(0.0, 0.0, dz, 0.0);
-    fVertexOutputArray->Add(candidate);
+    vertexcandidate = factory->NewCandidate();
+    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+    vertexcandidate->IsPU = 1;
+    
+    fVertexOutputArray->Add(vertexcandidate);
 
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
@@ -149,11 +185,11 @@
 
       candidate->IsPU = 1;
-
+    
       candidate->Momentum.SetPxPyPzE(px, py, pz, e);
       candidate->Momentum.RotateZ(dphi);
 
-      candidate->Position.SetXYZT(x, y, z + dz, t);
+      candidate->Position.SetXYZT(x, y, z+dz, t+dt);
       candidate->Position.RotateZ(dphi);
-
+      
       fParticleOutputArray->Add(candidate);
     }
Index: modules/PileUpMerger.h
===================================================================
--- modules/PileUpMerger.h	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/PileUpMerger.h	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -19,4 +19,5 @@
 class TObjArray;
 class DelphesPileUpReader;
+class DelphesTF2;
 
 class PileUpMerger: public DelphesModule
@@ -35,8 +36,12 @@
   Int_t fPileUpDistribution;
   Double_t fMeanPileUp;
+  
   Double_t fZVertexSpread;
+  Double_t fTVertexSpread;
 
   DelphesPileUpReader *fReader;
-
+  
+  DelphesTF2 *fFunction; //!
+  
   TIterator *fItInputArray; //!
 
Index: modules/TrackPileUpSubtractor.cc
===================================================================
--- modules/TrackPileUpSubtractor.cc	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/TrackPileUpSubtractor.cc	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -53,4 +53,9 @@
 void TrackPileUpSubtractor::Init()
 {
+// import input array
+
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+  
   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
 
@@ -85,4 +90,6 @@
     if(iterator) delete iterator;
   }
+
+  if(fItVertexInputArray) delete fItVertexInputArray;
 }
 
@@ -95,5 +102,19 @@
   TIterator *iterator;
   TObjArray *array;
-  Double_t z;
+  Double_t z, zvtx;
+
+  
+  // find z position of primary vertex
+  
+  fItVertexInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
+  {
+    if(candidate->IsPU == 0)
+    {
+    zvtx = candidate->Position.Z();
+    break;
+    }
+  }
+
 
   // loop over all input arrays
@@ -112,5 +133,5 @@
       // apply pile-up subtraction
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
-      if(candidate->IsPU && TMath::Abs(z) > fZVertexResolution) continue;
+      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
 
       array->Add(candidate);
Index: modules/TrackPileUpSubtractor.h
===================================================================
--- modules/TrackPileUpSubtractor.h	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/TrackPileUpSubtractor.h	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -39,4 +39,9 @@
 
   ClassDef(TrackPileUpSubtractor, 1)
+
+  TIterator *fItVertexInputArray; //!
+
+  const TObjArray *fVertexInputArray; //!
+
 };
 
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision af88093d0731bc574671224b51d187598cdec1e0)
+++ modules/TreeWriter.cc	(revision 22dc7fdf210b2c79eae0ce91d9aa4751b6469493)
@@ -95,5 +95,5 @@
       continue;
     }
-
+    
     itClassMap = fClassMap.find(branchClass);
     if(itClassMap == fClassMap.end())
@@ -115,5 +115,4 @@
 void TreeWriter::Finish()
 {
-
 }
 
@@ -161,5 +160,7 @@
   GenParticle *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all particles
   iterator.Reset();
@@ -208,5 +209,5 @@
     entry->Y = position.Y();
     entry->Z = position.Z();
-    entry->T = position.T();
+    entry->T = position.T()*1.0E-3/c_light;
   }
 }
@@ -219,4 +220,6 @@
   Candidate *candidate = 0;
   Vertex *entry = 0;
+  
+  const Double_t c_light = 2.99792458E8;
 
   // loop over all vertices
@@ -231,5 +234,5 @@
     entry->Y = position.Y();
     entry->Z = position.Z();
-    entry->T = position.T();
+    entry->T = position.T()*1.0E-3/c_light;
   }
 }
@@ -244,5 +247,6 @@
   Track *entry = 0;
   Double_t pt, signz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all tracks
   iterator.Reset();
@@ -271,4 +275,5 @@
     entry->YOuter = position.Y();
     entry->ZOuter = position.Z();
+    entry->TOuter = position.T()*1.0E-3/c_light;
 
     const TLorentzVector &momentum = candidate->Momentum;
@@ -290,4 +295,5 @@
     entry->Y = initialPosition.Y();
     entry->Z = initialPosition.Z();
+    entry->T = initialPosition.T()*1.0E-3/c_light;
 
     entry->Particle = particle;
@@ -303,5 +309,6 @@
   Tower *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all towers
   iterator.Reset();
@@ -309,5 +316,6 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
-
+    const TLorentzVector &position = candidate->Position;
+    
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -331,5 +339,7 @@
     entry->Edges[2] = candidate->Edges[2];
     entry->Edges[3] = candidate->Edges[3];
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+    
     FillParticles(candidate, &entry->Particles);
   }
@@ -344,5 +354,6 @@
   Photon *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -353,4 +364,6 @@
     TIter it1(candidate->GetCandidates());
     const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+    
 
     pt = momentum.Pt();
@@ -366,5 +379,7 @@
     entry->PT = pt;
     entry->E = momentum.E();
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+   
     entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
 
@@ -381,5 +396,6 @@
   Electron *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -389,5 +405,6 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
-
+    const TLorentzVector &position = candidate->Position;
+    
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -401,5 +418,7 @@
     entry->Phi = momentum.Phi();
     entry->PT = pt;
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+    
     entry->Charge = candidate->Charge;
 
@@ -418,5 +437,7 @@
   Muon *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -426,4 +447,6 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+    
 
     pt = momentum.Pt();
@@ -442,4 +465,8 @@
     entry->PT = pt;
 
+    entry->T = position.T()*1.0E-3/c_light;
+    
+   // cout<<entry->PT<<","<<entry->T<<endl;
+
     entry->Charge = candidate->Charge;
 
@@ -457,5 +484,6 @@
   Double_t pt, signPz, cosTheta, eta, rapidity;
   Double_t ecalEnergy, hcalEnergy;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -465,6 +493,8 @@
   {
     TIter itConstituents(candidate->GetCandidates());
-    const TLorentzVector &momentum = candidate->Momentum;
-
+    
+    const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+   
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -479,4 +509,6 @@
     entry->PT = pt;
 
+    entry->T = position.T()*1.0E-3/c_light;
+    
     entry->Mass = momentum.M();
 
