Index: Makefile
===================================================================
--- Makefile	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ Makefile	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -316,4 +316,5 @@
 	modules/EnergySmearing.h \
 	modules/MomentumSmearing.h \
+	modules/TrackSmearing.h \
 	modules/ImpactParameterSmearing.h \
 	modules/TimeSmearing.h \
@@ -338,4 +339,5 @@
 	modules/StatusPidFilter.h \
 	modules/PdgCodeFilter.h \
+	modules/BeamSpotFilter.h \
 	modules/Cloner.h \
 	modules/Weighter.h \
@@ -546,4 +548,13 @@
 	classes/DelphesFactory.h \
 	classes/DelphesFormula.h
+tmp/modules/BeamSpotFilter.$(ObjSuf): \
+	modules/BeamSpotFilter.$(SrcSuf) \
+	modules/BeamSpotFilter.h \
+	classes/DelphesClasses.h \
+	classes/DelphesFactory.h \
+	classes/DelphesFormula.h \
+	external/ExRootAnalysis/ExRootResult.h \
+	external/ExRootAnalysis/ExRootFilter.h \
+	external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/Calorimeter.$(ObjSuf): \
 	modules/Calorimeter.$(SrcSuf) \
@@ -830,4 +841,13 @@
 	modules/TrackPileUpSubtractor.$(SrcSuf) \
 	modules/TrackPileUpSubtractor.h \
+	classes/DelphesClasses.h \
+	classes/DelphesFactory.h \
+	classes/DelphesFormula.h \
+	external/ExRootAnalysis/ExRootResult.h \
+	external/ExRootAnalysis/ExRootFilter.h \
+	external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/TrackSmearing.$(ObjSuf): \
+	modules/TrackSmearing.$(SrcSuf) \
+	modules/TrackSmearing.h \
 	classes/DelphesClasses.h \
 	classes/DelphesFactory.h \
@@ -915,4 +935,5 @@
 	tmp/modules/AngularSmearing.$(ObjSuf) \
 	tmp/modules/BTagging.$(ObjSuf) \
+	tmp/modules/BeamSpotFilter.$(ObjSuf) \
 	tmp/modules/Calorimeter.$(ObjSuf) \
 	tmp/modules/Cloner.$(ObjSuf) \
@@ -946,4 +967,5 @@
 	tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
 	tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
+	tmp/modules/TrackSmearing.$(ObjSuf) \
 	tmp/modules/TreeWriter.$(ObjSuf) \
 	tmp/modules/UniqueObjectFinder.$(ObjSuf) \
@@ -1545,4 +1567,8 @@
 	tmp/external/tcl/tclVar.$(ObjSuf)
 
+modules/TrackSmearing.h: \
+	classes/DelphesModule.h
+	@touch $@
+
 external/fastjet/ClusterSequence.hh: \
 	external/fastjet/PseudoJet.hh \
@@ -1977,4 +2003,8 @@
 
 modules/FastJetFinder.h: \
+	classes/DelphesModule.h
+	@touch $@
+
+modules/BeamSpotFilter.h: \
 	classes/DelphesModule.h
 	@touch $@
Index: classes/DelphesClasses.cc
===================================================================
--- classes/DelphesClasses.cc	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ classes/DelphesClasses.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -127,6 +127,14 @@
   Momentum(0.0, 0.0, 0.0, 0.0),
   Position(0.0, 0.0, 0.0, 0.0),
+  InitialPosition(0.0, 0.0, 0.0, 0.0), 
   Area(0.0, 0.0, 0.0, 0.0),
-  Dxy(0), SDxy(0), Xd(0), Yd(0), Zd(0),
+  L(0),
+  D0(0), ErrorD0(0), 
+  DZ(0), ErrorDZ(0), 
+  P(0),  ErrorP(0), 
+  PT(0), ErrorPT(0), 
+  CtgTheta(0), ErrorCtgTheta(0), 
+  Phi(0), ErrorPhi(0),  
+  Xd(0), Yd(0), Zd(0), 
   TrackResolution(0),
   NCharged(0),
@@ -262,7 +270,19 @@
   object.Momentum = Momentum;
   object.Position = Position;
+  object.InitialPosition = Position;
   object.Area = Area;
-  object.Dxy = Dxy;
-  object.SDxy = SDxy;
+  object.L = L;
+  object.D0 = D0;
+  object.ErrorD0 = ErrorD0;
+  object.DZ = DZ;
+  object.ErrorDZ = ErrorDZ;
+  object.P = P;
+  object.ErrorP = ErrorP;
+  object.PT = PT;
+  object.ErrorPT = ErrorPT;
+  object.CtgTheta = CtgTheta ;
+  object.ErrorCtgTheta = ErrorCtgTheta;
+  object.Phi = Phi;
+  object.ErrorPhi = ErrorPhi;  
   object.Xd = Xd;
   object.Yd = Yd;
@@ -363,7 +383,19 @@
   Momentum.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
-  Dxy = 0.0;
-  SDxy = 0.0;
+  L = 0.0;
+  D0 = 0.0;  
+  ErrorD0 = 0.0;
+  DZ = 0.0;
+  ErrorDZ = 0.0;
+  P =0.0;
+  ErrorP =0.0;
+  PT = 0.0;
+  ErrorPT = 0.0;
+  CtgTheta = 0.0;
+  ErrorCtgTheta = 0.0;
+  Phi = 0.0;
+  ErrorPhi = 0.0;
   Xd = 0.0;
   Yd = 0.0;
Index: classes/DelphesClasses.h
===================================================================
--- classes/DelphesClasses.h	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ classes/DelphesClasses.h	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -415,6 +415,6 @@
   Float_t TOuter; // track position (z component) at the tracker edge
 
-  Float_t Dxy;     // track signed transverse impact parameter
-  Float_t SDxy;    // signed error on the track signed transverse impact parameter
+  Float_t D0;     // track signed transverse impact parameter
+  Float_t ErrorD0;    // signed error on the track signed transverse impact parameter
   Float_t Xd;      // X coordinate of point of closest approach to vertex
   Float_t Yd;      // Y coordinate of point of closest approach to vertex
@@ -526,8 +526,20 @@
   Float_t DeltaPhi;
 
-  TLorentzVector Momentum, Position, Area;
-
-  Float_t Dxy;
-  Float_t SDxy;
+  TLorentzVector Momentum, Position, InitialPosition, Area;
+
+  Float_t L; // path length
+  Float_t D0;
+  Float_t ErrorD0;
+  Float_t DZ;
+  Float_t ErrorDZ;
+  Float_t P;
+  Float_t ErrorP;
+  Float_t PT;
+  Float_t ErrorPT;
+  Float_t CtgTheta;
+  Float_t ErrorCtgTheta;
+  Float_t Phi;
+  Float_t ErrorPhi;
+
   Float_t Xd;
   Float_t Yd;
Index: modules/BeamSpotFilter.cc
===================================================================
--- modules/BeamSpotFilter.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
+++ modules/BeamSpotFilter.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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 cb4656839b9fb1ef75a65135a9839a1ba4dae840)
+++ modules/BeamSpotFilter.h	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/ImpactParameterSmearing.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/ModulesLinkDef.h	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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"
@@ -79,4 +81,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 +104,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+;
Index: modules/ParticlePropagator.cc
===================================================================
--- modules/ParticlePropagator.cc	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/ParticlePropagator.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -66,4 +66,5 @@
 {
 }
+
 
 //------------------------------------------------------------------------------
@@ -91,4 +92,8 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // import beamspot
+ 
+  fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+
   // create output arrays
 
@@ -111,5 +116,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 +125,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->GetSize () || !fBeamSpotInputArray->At(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 +147,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,4 +201,6 @@
       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;
@@ -186,5 +208,6 @@
 
       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 +261,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 + (py - 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,4 +324,10 @@
       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)
       {
@@ -295,6 +338,14 @@
 
         candidate->Momentum = candidateMomentum;
-        candidate->Dxy = dxy*1.0E3;
-        candidate->Xd = xd*1.0E3;
+        
+	candidate->L =  l*1.0E3;
+        candidate->D0 = d0*1.0E3;
+        candidate->DZ = dz*1.0E3;
+        candidate->P = p;
+        candidate->PT = pt;
+        candidate->CtgTheta = ctgTheta;
+        candidate->Phi = phi;
+
+	candidate->Xd = xd*1.0E3;
         candidate->Yd = yd*1.0E3;
         candidate->Zd = zd*1.0E3;
Index: modules/ParticlePropagator.h
===================================================================
--- modules/ParticlePropagator.h	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/ParticlePropagator.h	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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/TrackCountingBTagging.cc
===================================================================
--- modules/TrackCountingBTagging.cc	(revision 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/TrackCountingBTagging.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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 cb4656839b9fb1ef75a65135a9839a1ba4dae840)
+++ modules/TrackSmearing.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -0,0 +1,374 @@
+/** \class TrackSmearing
+ *
+ *  Performs d0, dZ, p, Theta, Phi smearing of tracks.
+ *
+ *
+ *
+ *  \author 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();
+ 
+  fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+
+  // 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->GetSize () || !fBeamSpotInputArray->At(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;
+    }
+    
+    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 cb4656839b9fb1ef75a65135a9839a1ba4dae840)
+++ modules/TrackSmearing.h	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -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 3b99510ce5c93329c4bbbac4f7aff27ab7c2d3c3)
+++ modules/TreeWriter.cc	(revision cb4656839b9fb1ef75a65135a9839a1ba4dae840)
@@ -292,6 +292,6 @@
     entry->TOuter = position.T()*1.0E-3/c_light;
 
-    entry->Dxy = candidate->Dxy;
-    entry->SDxy = candidate->SDxy ;
+    entry->D0 = candidate->D0;
+    entry->ErrorD0 = candidate->ErrorD0 ;
     entry->Xd = candidate->Xd;
     entry->Yd = candidate->Yd;
