Index: modules/AngularSmearing.cc
===================================================================
--- modules/AngularSmearing.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/AngularSmearing.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -98,20 +98,19 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e;
+  Double_t pt, eta, phi, e, m;
 
   fItInputArray->Reset();
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
     const TLorentzVector &candidateMomentum = candidate->Momentum;
-    eta = candidatePosition.Eta();
-    phi = candidatePosition.Phi();
+    eta = candidateMomentum.Eta();
+    phi = candidateMomentum.Phi();
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
+    m = candidateMomentum.M();
 
     // apply smearing formula for eta,phi
-
-    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e));
-    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e));
+    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e, candidate));
+    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e, candidate));
 
     if(pt <= 0.0) continue;
@@ -119,7 +118,5 @@
     mother = candidate;
     candidate = static_cast<Candidate *>(candidate->Clone());
-    eta = candidateMomentum.Eta();
-    phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
     candidate->AddCandidate(mother);
 
Index: modules/Calorimeter.cc
===================================================================
--- modules/Calorimeter.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/Calorimeter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -231,8 +231,13 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
   while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -450,4 +455,5 @@
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
   }
 
@@ -461,5 +467,5 @@
 {
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi;
+  Double_t energy, pt, eta, phi, r;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy;
@@ -511,5 +517,5 @@
   for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
   {
-    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
+    weight = TMath::Power((fTower->ECalEnergyTimePairs[i].first),2);
     sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
     sumWeight += weight;
@@ -517,11 +523,17 @@
   }
 
+  // check whether barrel or endcap tower
+  if (fTower->Position.Perp() < fTowerRmax && TMath::Abs(eta) > 0.)
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  else
+    r = fTower->Position.Pt();
+
   if(sumWeight > 0.0)
   {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime / sumWeight);
+    fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime / sumWeight);
   }
   else
   {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
+    fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9);
   }
 
@@ -559,5 +571,5 @@
   if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
   {
-    // create new photon tower
+    // create new photon tower assuming null mass
     tower = static_cast<Candidate *>(fTower->Clone());
     pt = ecalNeutralEnergy / TMath::CosH(eta);
@@ -646,6 +658,6 @@
       track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
       track->Momentum *= rescaleFactor;
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
 
       fEFlowTrackOutputArray->Add(track);
Index: modules/Calorimeter.h
===================================================================
--- modules/Calorimeter.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/Calorimeter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -60,4 +60,5 @@
   Double_t fTimingEnergyMin;
   Bool_t fElectronsFromTrack;
+  Double_t fTowerRmax;
 
   Int_t fTowerTrackHits, fTowerPhotonHits;
Index: modules/ClusterCounting.cc
===================================================================
--- modules/ClusterCounting.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/ClusterCounting.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,139 @@
+   /*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), 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 ClusterCounting
+ *
+ *  Counts ionisation clusters of energy loss in drift chambers
+ *
+ *  \authors F. Bedeschi - INFN Pisa
+*            P. Demin - UCLouvain, Louvain-la-Neuve
+ *           M. Selvaggi - CERN
+ *
+ *
+ */
+
+#include "modules/ClusterCounting.h"
+#include "classes/DelphesClasses.h"
+#include "TrackCovariance/TrkUtil.h"
+
+#include "TLorentzVector.h"
+#include "TVectorD.h"
+#include "TMath.h"
+#include "TObjArray.h"
+
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+ClusterCounting::ClusterCounting() :
+  fTrackUtil(0)
+{
+  fTrackUtil = new TrkUtil();
+}
+
+//------------------------------------------------------------------------------
+
+ClusterCounting::~ClusterCounting()
+{
+  if(fTrackUtil) delete fTrackUtil;
+}
+
+//------------------------------------------------------------------------------
+
+void ClusterCounting::Init()
+{
+
+  // geometric acceptance
+  fRmin  = GetDouble("Rmin", 0.);
+  fRmax  = GetDouble("Rmax", 0.);
+  fZmin  = GetDouble("Zmin", 0.);
+  fZmax  = GetDouble("Zmax", 0.);
+
+  // magnetic field
+  fBz    = GetDouble("Bz", 0.);
+
+  // gas mix option: 0
+  // 0:  Helium 90 - Isobutane 10
+  // 1:  Helium 100
+  // 2:  Argon 50 - Ethane 50
+  // 3:  Argon 100
+  fGasOption = GetInt("GasOption", 0);
+
+  // initialize drift chamber geometry and gas mix
+  fTrackUtil->SetBfield(fBz);
+  fTrackUtil->SetDchBoundaries(fRmin, fRmax, fZmin, fZmax);
+  fTrackUtil->SetGasMix(fGasOption);
+
+  // import input array
+  fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+}
+
+//------------------------------------------------------------------------------
+
+void ClusterCounting::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void ClusterCounting::Process()
+{
+  Candidate *candidate, *mother, *particle;
+  Double_t mass, trackLength, Ncl;
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+
+    // converting to meters
+    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+
+    // converting to meters
+    const TLorentzVector &candidatePosition = particle->Position*1e-03;
+    const TLorentzVector &candidateMomentum = particle->Momentum;
+
+		TVectorD Par = TrkUtil::XPtoPar(candidatePosition.Vect(), candidateMomentum.Vect(), candidate->Charge, fBz);
+    mass = candidateMomentum.M();
+
+    trackLength = fTrackUtil->TrkLen(Par);
+
+    mother    = candidate;
+    candidate = static_cast<Candidate*>(candidate->Clone());
+
+    Ncl = 0.;
+    if (fTrackUtil->IonClusters(Ncl, mass, Par))
+    {
+      candidate->Nclusters = Ncl;
+      candidate->dNdx = (trackLength > 0.) ? Ncl/trackLength : -1;
+    }
+
+    candidate->AddCandidate(mother);
+
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/ClusterCounting.h
===================================================================
--- modules/ClusterCounting.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/ClusterCounting.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,69 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), 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 ClusterCounting_h
+#define ClusterCounting_h
+
+/** \class ClusterCounting
+ *
+ *  Counts ionisation clusters of energy loss in drift chambers
+ *
+ *  \authors F. Bedeschi - INFN
+ *           M. Selvaggi - CERN
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+
+class TrkUtil;
+
+class ClusterCounting: public DelphesModule
+{
+public:
+  ClusterCounting();
+  ~ClusterCounting();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  Double_t fRmin;
+	Double_t fRmax;
+	Double_t fZmin;
+	Double_t fZmax;
+  Double_t fBz;
+
+  Int_t fGasOption;
+
+  TrkUtil *fTrackUtil;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(ClusterCounting, 1)
+};
+
+#endif
Index: modules/DecayFilter.cc
===================================================================
--- modules/DecayFilter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/DecayFilter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,136 @@
+/*
+ *  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 DecayFilter
+ *
+ *  This module randomly generates decays along the particle trajectory length 
+ *  according to actual particle decay length, taking into account for the boost
+ *  and using ROOT TDatabasePDG as a source for the particle lifetime.
+ *
+ *  This module is to to be used after a PropagateParticle step or a similar module
+ *  that calculates and store a trajectory length.
+ *
+ *  Particles that decay are not added to the OutputArray.
+ *
+ *  \author R. Preghenella - INFN, Bologna
+ *
+ */
+
+#include "modules/DecayFilter.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 "TParticlePDG.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;
+
+//------------------------------------------------------------------------------
+
+DecayFilter::DecayFilter() :
+  fItInputArray(0)
+{}
+
+//------------------------------------------------------------------------------
+
+DecayFilter::~DecayFilter()
+{}
+
+//------------------------------------------------------------------------------
+
+void DecayFilter::Init()
+{
+
+  // import input array(s)
+
+  fInputArray = ImportArray(GetString("InputArray", "FastJetFinder/jets"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array(s)
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+}
+
+//------------------------------------------------------------------------------
+
+void DecayFilter::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void DecayFilter::Process()
+{
+  Candidate *candidate;
+  TDatabasePDG *pdgdb = TDatabasePDG::Instance();
+  const Double_t c = TMath::C(); // [m/s]
+  Double_t m, t, p, bgct, L, l;
+  Bool_t hasDecayed = kFALSE;
+  
+  // loop over all input candidates
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+    // get particle information from PDG
+    TParticlePDG *pdg = pdgdb->GetParticle(candidate->PID);
+    if (!pdg) { // don't know this particle
+      fOutputArray->Add(candidate);
+      continue;
+    }    
+    m = pdg->Mass();
+    t = pdg->Lifetime(); // [s]
+    if (t == 0.) { // does not decay
+      fOutputArray->Add(candidate);
+      continue;
+    }
+
+    // compute boosted decay length (beta gamma c tau)
+    p = candidate->P;
+    bgct = p / m * c * t; // [m]
+
+    // get full trajectory length and generate random decay length
+    L = candidate->L * 1.0E-3; // [m]
+    l = gRandom->Exp(bgct);
+
+    // if random decay happens before end of trajectory, reject track
+    if (l < L) continue;
+
+    // else particle did not decay within the trajectory
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/DecayFilter.h
===================================================================
--- modules/DecayFilter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/DecayFilter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,64 @@
+/*
+ *  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 DecayFilter_h
+#define DecayFilter_h
+
+/** \class DecayFilter
+ *
+ *  This module randomly generates decays along the particle trajectory length 
+ *  according to actual particle decay length, taking into account for the boost
+ *  and using ROOT TDatabasePDG as a source for the particle lifetime.
+ *
+ *  This module is to to be used after a PropagateParticle step or a similar module
+ *  that calculates and store a trajectory length.
+ *
+ *  Particles that decay are not added to the OutputArray.
+ *
+ *  \author R. Preghenella - INFN, Bologna
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+#include <deque>
+
+class TObjArray;
+
+class DecayFilter: public DelphesModule
+{
+public:
+  DecayFilter();
+  ~DecayFilter();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(DecayFilter, 1)
+};
+
+#endif
Index: modules/Delphes.cc
===================================================================
--- modules/Delphes.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/Delphes.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -61,6 +61,9 @@
   fFactory(0)
 {
-  TFolder *folder = new TFolder(name, "");
+  TFolder *folder;
+
   fFactory = new DelphesFactory("ObjectFactory");
+
+  folder = new TFolder(name, "");
 
   SetName(name);
@@ -80,4 +83,5 @@
   if(folder)
   {
+    gROOT->GetListOfBrowsables()->Remove(folder);
     folder->Clear();
     delete folder;
Index: modules/DenseTrackFilter.cc
===================================================================
--- modules/DenseTrackFilter.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/DenseTrackFilter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -237,5 +237,5 @@
 
   Candidate *candidate, *track;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, m;
   Int_t numberOfCandidates;
 
@@ -251,7 +251,9 @@
   eta = candidate->Momentum.Eta();
   phi = candidate->Momentum.Phi();
+  m = candidate->Momentum.M();
+
   eta = gRandom->Gaus(eta, fEtaPhiRes);
   phi = gRandom->Gaus(phi, fEtaPhiRes);
-  candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+  candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
   candidate->AddCandidate(track);
 
Index: modules/DualReadoutCalorimeter.cc
===================================================================
--- modules/DualReadoutCalorimeter.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/DualReadoutCalorimeter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -24,5 +24,5 @@
   // This implementation of dual calorimetry relies on several approximations:
   // - If hadronic energy is found in the tower the energy resolution then the full tower enrgy is smeared according to hadronic resolution (pessimistic for (e,n) or (pi+,gamma))
-  // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles. 
+  // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles.
   //   Now we assume that regarless of the number of particle hits per tower we can always distinguish e+ vs pi+, which is probably not true in the case (e+,n) vs (pi+,gamma) without longitudinal segmentation.
 
@@ -63,5 +63,5 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  
+
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
@@ -75,5 +75,5 @@
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
-  
+
 }
 
@@ -82,5 +82,5 @@
 DualReadoutCalorimeter::~DualReadoutCalorimeter()
 {
-  
+
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
@@ -94,5 +94,5 @@
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
- 
+
 }
 
@@ -172,7 +172,9 @@
   fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
   fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
+  fEnergyMin = GetDouble("EnergyMin", 0.0);
 
   fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
   fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
+  fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
 
   // switch on or off the dithering of the center of DualReadoutCalorimeter towers
@@ -245,8 +247,16 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
+
+  //cout<<"--------- new event ---------- "<<endl;
+
   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -377,16 +387,19 @@
       fHCalTrackEnergy = 0.0;
       fTrackEnergy = 0.0;
-      
+
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
       fTrackSigma = 0.0;
-      
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
+
+      fTowerTime = 0.0;
+      fTowerTimeWeight = 0.0;
 
       fECalTowerTrackArray->Clear();
       fHCalTowerTrackArray->Clear();
       fTowerTrackArray->Clear();
-    
+
     }
 
@@ -412,40 +425,9 @@
       }
 
-      
-      /*
-      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
-      {
-        fECalTrackEnergy += ecalEnergy;
-        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
-        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
-        else energyGuess = momentum.E();
-
-        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
-        fECalTowerTrackArray->Add(track);
-      }
-     
-      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
-      {
-        fHCalTrackEnergy += hcalEnergy;
-        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
-        else energyGuess = momentum.E();
-
-        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
-        fHCalTowerTrackArray->Add(track);
-      }
-      
-      // muons
-      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
-      {
-        fEFlowTrackOutputArray->Add(track);
-      }
-      */
-      
       // in Dual Readout we do not care if tracks are ECAL of HCAL
       if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
-      { 
+      {
         fTrackEnergy += energy;
-        // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)      
+        // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)
         sigma = 0.0;
         if(fHCalTrackFractions[number] > 0)
@@ -453,13 +435,12 @@
         else
           sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-          
-        if(sigma/momentum.E() < track->TrackResolution) 
-          energyGuess = ecalEnergy + hcalEnergy;     
+
+        if(sigma/momentum.E() < track->TrackResolution)
+          energyGuess = ecalEnergy + hcalEnergy;
         else
           energyGuess = momentum.E();
-             
+
         fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
         fTowerTrackArray->Add(track);
-      
       }
       else
@@ -478,4 +459,5 @@
     position = particle->Position;
 
+
     // fill current tower
     ecalEnergy = momentum.E() * fECalTowerFractions[number];
@@ -485,13 +467,10 @@
     fHCalTowerEnergy += hcalEnergy;
 
-    if(ecalEnergy > fTimingEnergyMin && fTower)
-    {
-      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
-      {
-        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
-      }
-    }
+    // assume combined timing measurements in ECAL/HCAL sections
+    fTowerTime += (ecalEnergy + hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
+    fTowerTimeWeight += ecalEnergy + hcalEnergy;
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
   }
 
@@ -506,13 +485,13 @@
 
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi;
+  Double_t energy, pt, eta, phi, r, time;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
-  
+
   Double_t ecalSigma, hcalSigma, sigma;
   Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
 
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
-  
+
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
@@ -522,35 +501,20 @@
   if(!fTower) return;
 
-
-  //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
-  //cout<<"----------- New tower ---------"<<endl;
-
-
-  // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution.
-  // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
-
-  // if no hadronic fraction at all, then use ECAL resolution
-
-  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl;
-
-  // if no hadronic energy, use ECAL resolution 
+  // if no hadronic energy, use ECAL resolution
   if (fHCalTowerEnergy <= fHCalEnergyMin)
   {
     energy = fECalTowerEnergy;
     sigma  = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
-  }
-
-  // if hadronic fraction > 0, use HCAL resolution 
+  }
+
+  // if hadronic fraction > 0, use HCAL resolution
   else
   {
     energy = fECalTowerEnergy + fHCalTowerEnergy;
     sigma  = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
   }
 
   energy = LogNormal(energy, sigma);
-  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
-  
+
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
 
@@ -562,4 +526,6 @@
   hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
 
+  time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
+
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -568,6 +534,4 @@
   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
-  //cout<<"Measured energy: "<<energy<<endl;
-
   if(fSmearTowerCenter)
   {
@@ -583,30 +547,23 @@
   pt = energy / TMath::CosH(eta);
 
-  // Time calculation for tower
-  fTower->NTimeHits = 0;
-  sumWeightedTime = 0.0;
-  sumWeight = 0.0;
-
-  for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
-  {
-    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
-    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
-    sumWeight += weight;
-    fTower->NTimeHits++;
-  }
-
-  if(sumWeight > 0.0)
-  {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
-  }
-  else
-  {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
-  }
+  // check whether barrel or endcap tower
+
+  // endcap
+  if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  }
+  // barrel
+  else {
+    r = fTower->Position.Pt();
+  }
+
+  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
+  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->L = fTower->Position.Vect().Mag();
 
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
   fTower->Ehad = hcalEnergy;
-
+  fTower->Etrk = fTrackEnergy;
   fTower->Edges[0] = fTowerEdges[0];
   fTower->Edges[1] = fTowerEdges[1];
@@ -624,29 +581,32 @@
 
   // fill energy flow candidates
-  
+
   fTrackSigma = TMath::Sqrt(fTrackSigma);
   neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
   neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
 
-  //cout<<"trackEnergy: "<<fTrackEnergy<<", trackSigma: "<<fTrackSigma<<", Ntracks: "<<fTowerTrackArray->GetEntries()<<endl;
-
-  //cout<<"neutralEnergy: "<<neutralEnergy<<", neutralSigma: "<<neutralSigma<<endl;
-
-   // For now, if neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack !!! -> Creating only photons !! EFlowNeutralHadron collection will be empy!!! TO BE FIXED
   if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
   {
-    
-    //cout<<"significant neutral excess found:"<<endl;
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  neutralEnergy / TMath::CosH(eta);
-    //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl;
+    pt = neutralEnergy / TMath::CosH(eta);
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
-    tower->Eem = neutralEnergy;
-    tower->Ehad = 0.0;
-    tower->PID = 22;
-
-    fEFlowPhotonOutputArray->Add(tower);
-
+
+    // if no hadronic energy, use ECAL resolution
+    if (fHCalTowerEnergy <= fHCalEnergyMin)
+    {
+      tower->Eem = neutralEnergy;
+      tower->Ehad = 0.0;
+      tower->PID = 22;
+      fEFlowPhotonOutputArray->Add(tower);
+    }
+    // if hadronic fraction > 0, use HCAL resolution
+    else
+    {
+      tower->Eem = 0;
+      tower->Ehad = neutralEnergy;
+      tower->PID = 130;
+      fEFlowNeutralHadronOutputArray->Add(tower);
+    }
 
     //clone tracks
@@ -654,5 +614,4 @@
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
-      //cout<<"looping over tracks"<<endl;
       mother = track;
       track = static_cast<Candidate*>(track->Clone());
@@ -662,5 +621,5 @@
   }
 
-  
+
   // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
   else if(fTrackEnergy > 0.0)
@@ -670,5 +629,5 @@
     weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
 
-    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
+    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
     rescaleFactor = bestEnergyEstimate/fTrackEnergy;
 
@@ -676,129 +635,14 @@
     fItTowerTrackArray->Reset();
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {  
+    {
       mother = track;
-      track = static_cast<Candidate*>(track->Clone());
+      track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
       fEFlowTrackOutputArray->Add(track);
     }
   }
-  
-
-  /*
-  // fill energy flow candidates
-  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
-  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
-
-  //compute neutral excesses
-  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
-  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
-  
-  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
-  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
-  
-   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
-  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
-  {
-    // create new photon tower
-    tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  ecalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
-    tower->Eem = ecalNeutralEnergy;
-    tower->Ehad = 0.0;
-    tower->PID = 22;
-    
-    fEFlowPhotonOutputArray->Add(tower);
-   
-    //clone tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  
-  }
- 
-  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
-  else if(fECalTrackEnergy > 0.0)
-  {
-    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
-    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
-  
-    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 
-    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
-
-    //rescale tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-
-
-  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
-  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
-  {
-    // create new photon tower
-    tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  hcalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
-    tower->Ehad = hcalNeutralEnergy;
-    tower->Eem = 0.0;
-    
-    fEFlowNeutralHadronOutputArray->Add(tower);
-   
-    //clone tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  
-  }
- 
-  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
-  else if(fHCalTrackEnergy > 0.0)
-  {
-    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
-    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
-  
-    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 
-    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
-
-    //rescale tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-
-  */
+
+
 }
 
Index: modules/DualReadoutCalorimeter.h
===================================================================
--- modules/DualReadoutCalorimeter.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/DualReadoutCalorimeter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -60,4 +60,5 @@
   Double_t fECalTrackEnergy, fHCalTrackEnergy;
   Double_t fTrackEnergy;
+  Double_t fTowerRmax;
 
   Double_t fTimingEnergyMin;
@@ -77,4 +78,7 @@
   Double_t fHCalTrackSigma;
   Double_t fTrackSigma;
+
+  Double_t fTowerTime;
+  Double_t fTowerTimeWeight;
 
   Bool_t fSmearTowerCenter;
Index: modules/Efficiency.cc
===================================================================
--- modules/Efficiency.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/Efficiency.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -78,4 +78,7 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // switch to compute efficiency based on momentum vector eta, phi
+  fUseMomentumVector = GetBool("UseMomentumVector", false);
+
   // create output array
 
@@ -95,5 +98,6 @@
 {
   Candidate *candidate;
-  Double_t pt, eta, phi, e, d0, dz, ctgTheta;
+  Double_t pt, eta, phi, e;
+
 
   fItInputArray->Reset();
@@ -104,12 +108,15 @@
     eta = candidatePosition.Eta();
     phi = candidatePosition.Phi();
+
+    if (fUseMomentumVector){
+      eta = candidateMomentum.Eta();
+      phi = candidateMomentum.Phi();
+    }
+
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    d0 = candidate->D0;
-    dz = candidate->DZ;
-    ctgTheta = candidate->CtgTheta;
 
     // apply an efficency formula
-    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta)) continue;
+    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, candidate)) continue;
 
     fOutputArray->Add(candidate);
Index: modules/Efficiency.h
===================================================================
--- modules/Efficiency.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/Efficiency.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -53,4 +53,6 @@
   TObjArray *fOutputArray; //!
 
+  Double_t fUseMomentumVector; //!
+
   ClassDef(Efficiency, 1)
 };
Index: modules/EnergySmearing.cc
===================================================================
--- modules/EnergySmearing.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/EnergySmearing.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -95,5 +95,5 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, energy, eta, phi;
+  Double_t pt, energy, eta, phi, m;
 
   fItInputArray->Reset();
@@ -107,4 +107,5 @@
     phi = candidatePosition.Phi();
     energy = candidateMomentum.E();
+    m = candidateMomentum.M();
 
     // apply smearing formula
@@ -117,5 +118,6 @@
     eta = candidateMomentum.Eta();
     phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(energy / TMath::CosH(eta), eta, phi, energy);
+    pt = (energy > m) ? TMath::Sqrt(energy*energy - m*m)/TMath::CosH(eta) : 0;
+    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     candidate->TrackResolution = fFormula->Eval(pt, eta, phi, energy) / candidateMomentum.E();
     candidate->AddCandidate(mother);
Index: modules/FastJetFinder.cc
===================================================================
--- modules/FastJetFinder.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/FastJetFinder.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -104,4 +104,5 @@
   fJetAlgorithm = GetInt("JetAlgorithm", 6);
   fParameterR = GetDouble("ParameterR", 0.5);
+  fParameterP = GetDouble("ParameterP", -1.0);
 
   fConeRadius = GetDouble("ConeRadius", 0.5);
@@ -247,4 +248,8 @@
     fDefinition = new JetDefinition(fValenciaPlugin);
     break;
+  case 10:
+    fDefinition = new JetDefinition(ee_genkt_algorithm,fParameterR,fParameterP);
+    break;
+
   }
 
@@ -316,4 +321,6 @@
   Double_t deta, dphi, detaMax, dphiMax;
   Double_t time, timeWeight;
+  Double_t neutralEnergyFraction, chargedEnergyFraction;
+
   Int_t number, ncharged, nneutrals;
   Int_t charge;
@@ -418,4 +425,7 @@
     nneutrals = 0;
 
+    neutralEnergyFraction =0.;
+    chargedEnergyFraction =0.;
+
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
@@ -432,8 +442,14 @@
 
       if(constituent->Charge == 0)
+      {
         nneutrals++;
+        neutralEnergyFraction += constituent->Momentum.E();
+      }
       else
+      {
         ncharged++;
-
+        chargedEnergyFraction += constituent->Momentum.E();
+      }
+      
       time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
       timeWeight += TMath::Sqrt(constituent->Momentum.E());
@@ -455,4 +471,7 @@
     candidate->NCharged = ncharged;
 
+    candidate->NeutralEnergyFraction = (momentum.E() > 0 ) ? neutralEnergyFraction/momentum.E() : 0.0;
+    candidate->ChargedEnergyFraction = (momentum.E() > 0 ) ? chargedEnergyFraction/momentum.E() : 0.0;
+
     //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
     candidate->ExclYmerge23 = excl_ymerge23;
@@ -470,7 +489,5 @@
       fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
       fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
-
-      trimmed_jet = join(trimmed_jet.constituents());
-
+      
       candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
 
Index: modules/FastJetFinder.h
===================================================================
--- modules/FastJetFinder.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/FastJetFinder.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -72,4 +72,5 @@
   Int_t fJetAlgorithm;
   Double_t fParameterR;
+  Double_t fParameterP;
 
   Double_t fJetPTMin;
Index: modules/IdentificationMap.cc
===================================================================
--- modules/IdentificationMap.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/IdentificationMap.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -170,4 +170,5 @@
       {
         // change PID of particle
+        candidate = static_cast<Candidate *>(candidate->Clone());
         if(pdgCodeOut != 0) candidate->PID = charge * pdgCodeOut;
         fOutputArray->Add(candidate);
Index: modules/JetFakeParticle.cc
===================================================================
--- modules/JetFakeParticle.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/JetFakeParticle.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -146,8 +146,7 @@
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
     const TLorentzVector &candidateMomentum = candidate->Momentum;
-    eta = candidatePosition.Eta();
-    phi = candidatePosition.Phi();
+    eta = candidateMomentum.Eta();
+    phi = candidateMomentum.Phi();
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/ModulesLinkDef.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -36,6 +36,9 @@
 #include "modules/MomentumSmearing.h"
 #include "modules/TrackSmearing.h"
+#include "modules/TrackCovariance.h"
+#include "modules/ClusterCounting.h"
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
+#include "modules/TimeOfFlight.h"
 #include "modules/SimpleCalorimeter.h"
 #include "modules/DenseTrackFilter.h"
@@ -72,4 +75,7 @@
 #include "modules/VertexFinder.h"
 #include "modules/VertexFinderDA4D.h"
+#include "modules/DecayFilter.h"
+#include "modules/ParticleDensity.h"
+#include "modules/TruthVertexFinder.h"
 #include "modules/ExampleModule.h"
 #include "modules/LLPFilter.h"
@@ -93,6 +99,9 @@
 #pragma link C++ class MomentumSmearing+;
 #pragma link C++ class TrackSmearing+;
+#pragma link C++ class TrackCovariance+;
+#pragma link C++ class ClusterCounting+;
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
+#pragma link C++ class TimeOfFlight+;
 #pragma link C++ class SimpleCalorimeter+;
 #pragma link C++ class DenseTrackFilter+;
@@ -129,4 +138,7 @@
 #pragma link C++ class VertexFinder+;
 #pragma link C++ class VertexFinderDA4D+;
+#pragma link C++ class DecayFilter+;
+#pragma link C++ class ParticleDensity+;
+#pragma link C++ class TruthVertexFinder+;
 #pragma link C++ class ExampleModule+;
 #pragma link C++ class LLPFilter+;
Index: modules/MomentumSmearing.cc
===================================================================
--- modules/MomentumSmearing.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/MomentumSmearing.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -78,4 +78,7 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // switch to compute momentum smearing based on momentum vector eta, phi
+  fUseMomentumVector = GetBool("UseMomentumVector", false);
+
   // create output array
 
@@ -95,5 +98,5 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e, res;
+  Double_t pt, eta, phi, e, m, res;
 
   fItInputArray->Reset();
@@ -104,7 +107,14 @@
     eta = candidatePosition.Eta();
     phi = candidatePosition.Phi();
+
+    if (fUseMomentumVector){
+      eta = candidateMomentum.Eta();
+      phi = candidateMomentum.Phi();
+    }
+
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    res = fFormula->Eval(pt, eta, phi, e);
+    m = candidateMomentum.M();
+    res = fFormula->Eval(pt, eta, phi, e, candidate);
 
     // apply smearing formula
@@ -121,5 +131,5 @@
     eta = candidateMomentum.Eta();
     phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
     //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     candidate->TrackResolution = res;
Index: modules/MomentumSmearing.h
===================================================================
--- modules/MomentumSmearing.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/MomentumSmearing.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -55,4 +55,6 @@
   TObjArray *fOutputArray; //!
 
+  Double_t fUseMomentumVector; //!
+
   ClassDef(MomentumSmearing, 1)
 };
Index: modules/ParticleDensity.cc
===================================================================
--- modules/ParticleDensity.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/ParticleDensity.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,144 @@
+/*
+ *  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 ParticleDensity
+ *
+ *  This module calculates the particle multiplicity density in eta-phi bins.
+ *  It then assigns the value to the candidates according to the candidate eta.
+ *
+ *  \author R. Preghenella - INFN, Bologna
+ *
+ */
+
+#include "modules/ParticleDensity.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 "TFormula.h"
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
+#include "TString.h"
+#include "TH2F.h"
+
+#include <algorithm>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+ParticleDensity::ParticleDensity() :
+  fItInputArray(0)
+{}
+
+//------------------------------------------------------------------------------
+
+ParticleDensity::~ParticleDensity()
+{}
+
+//------------------------------------------------------------------------------
+
+void ParticleDensity::Init()
+{
+
+  // import input array(s)
+
+  fInputArray = ImportArray(GetString("InputArray", "FastJetFinder/jets"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array(s)
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+
+  // create multiplicity histogram
+
+  ExRootConfParam paramEta = GetParam("EtaBins");
+  const Long_t sizeEta = paramEta.GetSize();
+  Int_t nbinsEta = sizeEta - 1;
+  Float_t binsEta[sizeEta];
+  for (Int_t i = 0; i < sizeEta; ++i) {
+    binsEta[i] = paramEta[i].GetDouble();
+  }
+  
+  ExRootConfParam paramPhi = GetParam("PhiBins");
+  const Long_t sizePhi = paramPhi.GetSize();
+  Int_t nbinsPhi = sizePhi - 1;
+  Float_t binsPhi[sizePhi];
+  for (Int_t i = 0; i < sizePhi; ++i) {
+    binsPhi[i] = paramPhi[i].GetDouble();
+  }
+  
+  fHisto = new TH2F("hParticleDensity", ";#eta;#varphi;d^{2}N/d#etad#varphi", nbinsEta, binsEta, nbinsPhi, binsPhi);
+
+  fUseMomentumVector = GetBool("UseMomentumVector", false);
+}
+
+//------------------------------------------------------------------------------
+
+void ParticleDensity::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+  if (fHisto) delete fHisto;
+}
+
+//------------------------------------------------------------------------------
+
+void ParticleDensity::Process()
+{
+  Candidate *candidate;
+  fHisto->Reset();
+  
+  // loop over all input candidates to fill histogram
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) {
+    if (fUseMomentumVector) fHisto->Fill(candidate->Momentum.Eta(), candidate->Momentum.Phi());
+    else                    fHisto->Fill(candidate->Position.Eta(), candidate->Position.Phi());
+  }
+
+  // normalise by bin width
+  fHisto->Scale(1., "width");
+  
+  // loop over all input candidates to assign multiplicity
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) {
+    Int_t ieta, iphi;
+    if (fUseMomentumVector) {
+      ieta = fHisto->GetXaxis()->FindBin(candidate->Momentum.Eta());
+      iphi = fHisto->GetYaxis()->FindBin(candidate->Momentum.Phi());
+    } else {
+      ieta = fHisto->GetXaxis()->FindBin(candidate->Position.Eta());
+      iphi = fHisto->GetYaxis()->FindBin(candidate->Position.Phi());
+    }
+    candidate->ParticleDensity = fHisto->GetBinContent(ieta, iphi);
+    fOutputArray->Add(candidate);
+  }
+
+  
+}
+
+//------------------------------------------------------------------------------
Index: modules/ParticleDensity.h
===================================================================
--- modules/ParticleDensity.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/ParticleDensity.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -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 ParticleDensity_h
+#define ParticleDensity_h
+
+/** \class ParticleDensity
+ *
+ *  This module calculates the particle multiplicity density in eta-phi bins.
+ *  It then assigns the value to the candidates according to the candidate eta.
+ *
+ *  \author R. Preghenella - INFN, Bologna
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+#include <deque>
+
+class TObjArray;
+class TH2F;
+
+class ParticleDensity: public DelphesModule
+{
+public:
+  ParticleDensity();
+  ~ParticleDensity();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  Bool_t fUseMomentumVector; // !
+  TH2F *fHisto; //!
+  
+  ClassDef(ParticleDensity, 1)
+};
+
+#endif
Index: modules/ParticlePropagator.cc
===================================================================
--- modules/ParticlePropagator.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/ParticlePropagator.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -125,19 +125,21 @@
   TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
   Double_t px, py, pz, pt, pt2, e, q;
-  Double_t x, y, z, t, r, phi;
-  Double_t x_c, y_c, r_c, phi_c, phi_0;
-  Double_t x_t, y_t, z_t, r_t;
-  Double_t t1, t2, t3, t4, t5, t6;
-  Double_t t_z, t_r, t_ra, t_rb;
-  Double_t tmp, discr, discr2;
-  Double_t delta, gammam, omega, asinrho;
-  Double_t rcu, rc2, xd, yd, zd;
-  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
+  Double_t x, y, z, t, r;
+  Double_t x_c, y_c, r_c, phi_0;
+  Double_t x_t, y_t, z_t, r_t, phi_t;
+  Double_t t_r, t_z;
+  Double_t tmp;
+  Double_t gammam, omega;
+  Double_t xd, yd, zd;
+  Double_t l, d0, dz, ctgTheta, alpha;
   Double_t bsx, bsy, bsz;
+  Double_t td, pio, phid, vz;
 
   const Double_t c_light = 2.99792458E8;
 
   if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
+  {
     beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  }
   else
   {
@@ -160,4 +162,5 @@
     particlePosition = particle->Position;
     particleMomentum = particle->Momentum;
+
     x = particlePosition.X() * 1.0E-3;
     y = particlePosition.Y() * 1.0E-3;
@@ -206,25 +209,9 @@
       // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
       tmp = px * y - py * x;
-      discr2 = pt2 * fRadius2 - tmp * tmp;
-
-      if(discr2 < 0.0)
-      {
-        // no solutions
-        continue;
-      }
-
-      tmp = px * x + py * y;
-      discr = TMath::Sqrt(discr2);
-      t1 = (-tmp + discr) / pt2;
-      t2 = (-tmp - discr) / pt2;
-      t = (t1 < 0.0) ? t2 : t1;
-
-      z_t = z + pz * t;
-      if(TMath::Abs(z_t) > fHalfLength)
-      {
-        t3 = (+fHalfLength - z) / pz;
-        t4 = (-fHalfLength - z) / pz;
-        t = (t3 < 0.0) ? t4 : t3;
-      }
+      t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
+
+      t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
+
+      t = TMath::Min(t_r, t_z);
 
       x_t = x + px * t;
@@ -245,4 +232,5 @@
 
       fOutputArray->Add(candidate);
+
       if(TMath::Abs(q) > 1.0E-9)
       {
@@ -267,12 +255,12 @@
     {
 
-      // 1.  initial transverse momentum p_{T0}: Part->pt
-      //     initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
-      //     relativistic gamma: gamma = E/mc^2; gammam = gamma * m
-      //     gyration frequency omega = q/(gamma m) fBz
-      //     helix radius r = p_{T0} / (omega gamma m)
+      // 1. initial transverse momentum p_{T0}: Part->pt
+      //    initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
+      //    relativistic gamma: gamma = E / mc^2; gammam = gamma * m
+      //    gyration frequency omega = q * Bz / (gammam)
+      //    helix radius r = p_{T0} / (omega * gammam)
 
       gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
-      omega = q * fBz / (gammam); // omega is here in [89875518/s]
+      omega = q * fBz / gammam; // omega is here in [89875518/s]
       r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
 
@@ -283,33 +271,32 @@
       y_c = y - r * TMath::Cos(phi_0);
       r_c = TMath::Hypot(x_c, y_c);
-      phi_c = TMath::ATan2(y_c, x_c);
-      phi = phi_c;
-      if(x_c < 0.0) phi += TMath::Pi();
-
-      rcu = TMath::Abs(r);
-      rc2 = r_c * r_c;
-
-      // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
-      xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
-      xd = (rc2 > 0.0) ? xd / rc2 : -999;
-      yd = y_c * (-rcu * r_c + rc2);
-      yd = (rc2 > 0.0) ? yd / rc2 : -999;
-      zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / 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 = particleMomentum.Eta();
-      phip = TMath::ATan2(py, px);
-
-      particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
+
+      // time of closest approach
+      td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
+
+      // remove all the modulo pi that might have come from the atan
+      pio = TMath::Abs(TMath::Pi() / omega);
+      while(TMath::Abs(td) > 0.5 * pio)
+      {
+        td -= TMath::Sign(1.0, td) * pio;
+      }
+
+      vz = pz * c_light / e;
+
+      // calculate coordinates of closest approach to z axis
+      phid = phi_0 - omega * td;
+      xd = x_c - r * TMath::Sin(phid);
+      yd = y_c + r * TMath::Cos(phid);
+      zd = z + vz * td;
+
+      // momentum at closest approach
+      px = pt * TMath::Cos(phid);
+      py = pt * TMath::Sin(phid);
+
+      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
 
       // calculate additional track parameters (correct for beamspot position)
-
-      d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
-      dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
-      p = particleMomentum.P();
+      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
+      dz = zd - bsz;
       ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
 
@@ -317,10 +304,5 @@
       //    t_r : time to exit from the sides
       //    t_z : time to exit from the front or the back
-      t_r = 0.0; // in [ns]
-      int sign_pz = (pz > 0.0) ? 1 : -1;
-      if(pz == 0.0)
-        t_z = 1.0E99;
-      else
-        t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
+      t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
 
       if(r_c + TMath::Abs(r) < fRadius)
@@ -331,42 +313,22 @@
       else
       {
-        asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
-        delta = phi_0 - phi;
-        if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
-        if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
-        t1 = (delta + asinrho) / omega;
-        t2 = (delta + TMath::Pi() - asinrho) / omega;
-        t3 = (delta + TMath::Pi() + asinrho) / omega;
-        t4 = (delta - asinrho) / omega;
-        t5 = (delta - TMath::Pi() - asinrho) / omega;
-        t6 = (delta - TMath::Pi() + asinrho) / omega;
-
-        if(t1 < 0.0) t1 = 1.0E99;
-        if(t2 < 0.0) t2 = 1.0E99;
-        if(t3 < 0.0) t3 = 1.0E99;
-        if(t4 < 0.0) t4 = 1.0E99;
-        if(t5 < 0.0) t5 = 1.0E99;
-        if(t6 < 0.0) t6 = 1.0E99;
-
-        t_ra = TMath::Min(t1, TMath::Min(t2, t3));
-        t_rb = TMath::Min(t4, TMath::Min(t5, t6));
-        t_r = TMath::Min(t_ra, t_rb);
+        alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
+        t_r = td + TMath::Abs(alpha / omega);
+
         t = TMath::Min(t_r, t_z);
       }
 
       // 4. position in terms of x(t), y(t), z(t)
-      x_t = x_c + r * TMath::Sin(omega * t - phi_0);
-      y_t = y_c + r * TMath::Cos(omega * t - phi_0);
-      z_t = z + pz * 1.0E9 / c_light / gammam * t;
+      phi_t = phi_0 - omega * t;
+      x_t = x_c - r * TMath::Sin(phi_t);
+      y_t = y_c + r * TMath::Cos(phi_t);
+      z_t = z + vz * t;
       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);
+      // lenght of the path from production to tracker
+      l = t * TMath::Hypot(vz, r * omega);
 
       if(r_t > 0.0)
       {
-
         // store these variables before cloning
         if(particle == candidate)
@@ -374,8 +336,8 @@
           particle->D0 = d0 * 1.0E3;
           particle->DZ = dz * 1.0E3;
-          particle->P = p;
+          particle->P = particleMomentum.P();
           particle->PT = pt;
           particle->CtgTheta = ctgTheta;
-          particle->Phi = phip;
+          particle->Phi = particleMomentum.Phi();
         }
 
Index: modules/SimpleCalorimeter.cc
===================================================================
--- modules/SimpleCalorimeter.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/SimpleCalorimeter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -208,8 +208,13 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
   while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -394,8 +399,10 @@
     fTowerEnergy += energy;
 
-    fTowerTime += energy * position.T();
-    fTowerTimeWeight += energy;
+    fTowerTime += energy * energy * position.T(); //sigma_t ~ 1/E
+    fTowerTimeWeight += energy * energy;
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
+
   }
 
@@ -409,5 +416,5 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy, neutralEnergy, pt, eta, phi;
+  Double_t energy, neutralEnergy, pt, eta, phi, r;
   Double_t sigma, neutralSigma;
   Double_t time;
@@ -443,9 +450,20 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  // endcap
+  if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  }
+  // barrel
+  else {
+    r = fTower->Position.Pt();
+  }
+
+  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->L = fTower->Position.Vect().Mag();
 
   fTower->Eem = (!fIsEcal) ? 0 : energy;
   fTower->Ehad = (fIsEcal) ? 0 : energy;
+  fTower->Etrk = fTrackEnergy;
 
   fTower->Edges[0] = fTowerEdges[0];
@@ -507,7 +525,5 @@
       track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
       fEFlowTrackOutputArray->Add(track);
     }
Index: modules/SimpleCalorimeter.h
===================================================================
--- modules/SimpleCalorimeter.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/SimpleCalorimeter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -62,4 +62,6 @@
   Double_t fTrackTime;
 
+  Double_t fTowerRmax;
+
   Double_t fTowerTimeWeight;
   Double_t fTrackTimeWeight;
Index: modules/TauTagging.cc
===================================================================
--- modules/TauTagging.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TauTagging.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -83,7 +83,7 @@
     daughter1 = static_cast<Candidate *>(fParticleInputArray->At(i));
     pdgCode = TMath::Abs(daughter1->PID);
-    if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15)
-      return -1;
-    else if(pdgCode == 24)
+    //if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15)
+    //  return -1;
+    if(pdgCode == 24)
     {
       if(daughter1->D1 < 0) return -1;
@@ -96,5 +96,4 @@
     }
   }
-
   return 0;
 }
@@ -190,5 +189,5 @@
 void TauTagging::Process()
 {
-  Candidate *jet, *tau, *daughter;
+  Candidate *jet, *tau, *daughter, *part;
   TLorentzVector tauMomentum;
   Double_t pt, eta, phi, e, eff;
@@ -204,6 +203,8 @@
   // loop over all input jets
   fItJetInputArray->Reset();
+
   while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
   {
+
     const TLorentzVector &jetMomentum = jet->Momentum;
     pdgCode = 0;
@@ -243,4 +244,31 @@
       }
     }
+    
+    // fake electrons and muons
+    
+    if (pdgCode == 0)
+    {
+     
+      Double_t drMin = fDeltaR;   
+      fItPartonInputArray->Reset();
+      while((part = static_cast<Candidate *>(fItPartonInputArray->Next())))
+      {
+        if(TMath::Abs(part->PID) == 11 || TMath::Abs(part->PID) == 13) 
+        {
+            tauMomentum = part->Momentum;
+            if (tauMomentum.Pt() < fClassifier->fPTMin) continue;
+            if (TMath::Abs(tauMomentum.Eta()) > fClassifier->fEtaMax) continue;
+
+            Double_t dr = jetMomentum.DeltaR(tauMomentum);
+            if( dr < drMin)
+            {
+               drMin = dr;
+               pdgCode = TMath::Abs(part->PID);
+               charge = part->Charge;
+            }  
+        }
+      }
+    }
+
     // find an efficency formula
     itEfficiencyMap = fEfficiencyMap.find(pdgCode);
Index: modules/TimeOfFlight.cc
===================================================================
--- modules/TimeOfFlight.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TimeOfFlight.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,223 @@
+/*
+ *  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 TimeOfFlight
+  *
+  *  Calculates Time-Of-Flight
+  *
+  *  \author Michele Selvaggi - CERN
+  *
+ */
+
+#include "modules/TimeOfFlight.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;
+//------------------------------------------------------------------------------
+
+TimeOfFlight::TimeOfFlight() :
+  fItInputArray(0), fItVertexInputArray(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+TimeOfFlight::~TimeOfFlight()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void TimeOfFlight::Init()
+{
+
+  // method to compute vertex time
+  fVertexTimeMode = GetInt("VertexTimeMode", 0);
+
+  // import track input array
+  fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // import vertex input array
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "TruthVertexFinder/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+
+  // create output array
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+}
+
+//------------------------------------------------------------------------------
+
+void TimeOfFlight::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+  if(fItVertexInputArray) delete fItVertexInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TimeOfFlight::Process()
+{
+  Candidate *candidate, *particle, *vertex, *constituent, *mother;
+  Double_t ti, t_truth, tf;
+  Double_t l, tof, beta, p,  mass;
+
+  const Double_t c_light = 2.99792458E8;
+
+  // first compute momenta of vertices based on reconstructed tracks
+  ComputeVertexMomenta();
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+
+    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+
+    const TLorentzVector &candidateInitialPosition = particle->Position;
+    const TLorentzVector &candidateInitialPositionSmeared = candidate->InitialPosition;
+    const TLorentzVector &candidateFinalPosition = candidate->Position;
+    const TLorentzVector &candidateMomentum = particle->Momentum;
+
+    // time at vertex from MC truth
+    t_truth = candidateInitialPosition.T() * 1.0E-3 / c_light;
+
+    // various options on how to calculate the vertex time
+    ti=0;
+    switch (fVertexTimeMode)
+	  {
+    	case 0:
+    	{
+        // assume ti from MC truth
+        // most aggressive, we are cheating and assume we can perfectly reconstruct time of primary and secondary vertices
+        ti = t_truth;
+        break;
+    	}
+    	case 1:
+    	{
+        // always assume t=0, most conservative assumption
+        // reasonable assumption for particles originating from PV, if beamSpot has small time spread compared to timing resolution
+        // probably bad assumption for particles from highly displaced vertices (i.e Ks)
+        ti=0;
+        break;
+    	}
+    	case 2:
+    	{
+        // same as 2 but attempt at estimate beta from vertex mass and momentum
+        beta = 1.;
+        fItVertexInputArray->Reset();
+        while((vertex = static_cast<Candidate *>(fItVertexInputArray->Next())))
+        {
+          TIter itGenParts(vertex->GetCandidates());
+          itGenParts.Reset();
+
+          while((constituent = static_cast<Candidate *>(itGenParts.Next())))
+          {
+            if (particle == constituent)
+            {
+              beta = vertex->Momentum.Beta();
+              break;
+            }
+          }
+        } // end vertex  loop
+
+        // track displacement to be possibily replaced by vertex fitted position
+        ti = candidateInitialPositionSmeared.Vect().Mag() * 1.0E-3 /(beta*c_light);
+      }
+      break;
+  	}
+
+    p = candidateMomentum.P();
+
+    // this quantity has already been smeared by another module
+    tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
+
+    // calculate time-of-flight
+    tof = tf - ti;
+    // path length of the full helix
+    l = candidate->L * 1.0E-3;
+
+    // particle velocity
+    beta = l/(c_light*tof);
+
+    // calculate particle mass (i.e particle ID)
+    mass = 0.;
+    if (beta<1) mass = p* TMath::Sqrt(1/(beta*beta) - 1);
+
+    mother    = candidate;
+    candidate = static_cast<Candidate*>(candidate->Clone());
+
+    // update time at vertex based on option
+    candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
+
+    // update particle mass based on TOF-based PID (commented for now, assume this calculation is done offline)
+    //candidate->Momentum.SetVectM(candidateMomentum.Vect(), mass);
+
+    candidate->AddCandidate(mother);
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void TimeOfFlight::ComputeVertexMomenta()
+{
+  Candidate *track, *constituent, *particle, *vertex;
+
+  fItVertexInputArray->Reset();
+  while((vertex = static_cast<Candidate *>(fItVertexInputArray->Next())))
+  {
+
+    TIter itGenParts(vertex->GetCandidates());
+    itGenParts.Reset();
+
+    while((constituent = static_cast<Candidate *>(itGenParts.Next())))
+    {
+      fItInputArray->Reset();
+      while((track = static_cast<Candidate *>(fItInputArray->Next())))
+      {
+        // get gen part that generated track
+        particle = static_cast<Candidate *>(track->GetCandidates()->At(0));
+        if (particle == constituent)
+        {
+          vertex->Momentum += track->Momentum;
+        }
+
+      } // end track loop
+    } // end vertex consitutent loop
+  } // end vertex  loop
+}
Index: modules/TimeOfFlight.h
===================================================================
--- modules/TimeOfFlight.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TimeOfFlight.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,61 @@
+/*
+ *  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 TimeOfFlight_h
+#define TimeOfFlight_h
+
+/** \class TimeOfFlight
+ *
+ *  Calculates Time-Of-Flight
+ *
+ *  \author Michele Selvaggi - CERN
+ *
+*/
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+
+class TimeOfFlight: public DelphesModule
+{
+public:
+  TimeOfFlight();
+  ~TimeOfFlight();
+
+  void Init();
+  void Process();
+  void Finish();
+  void ComputeVertexMomenta();
+
+private:
+
+  Int_t fVertexTimeMode;
+
+  TIterator *fItInputArray; //!
+  TIterator *fItVertexInputArray; //!
+
+  const TObjArray *fInputArray; //!
+  const TObjArray *fVertexInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(TimeOfFlight, 1)
+};
+
+#endif
Index: modules/TimeSmearing.cc
===================================================================
--- modules/TimeSmearing.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TimeSmearing.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -19,7 +19,7 @@
 /** \class TimeSmearing
  *
- *  Performs transverse momentum resolution smearing.
+ *  Performs time smearing.
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author M. Selvaggi - CERN
  *
  */
@@ -49,10 +49,10 @@
 
 using namespace std;
-
 //------------------------------------------------------------------------------
 
 TimeSmearing::TimeSmearing() :
-  fItInputArray(0)
+  fItInputArray(0), fResolutionFormula(0)
 {
+	fResolutionFormula = new DelphesFormula;
 }
 
@@ -61,4 +61,5 @@
 TimeSmearing::~TimeSmearing()
 {
+	if(fResolutionFormula) delete fResolutionFormula;
 }
 
@@ -69,13 +70,14 @@
   // read resolution formula
 
-  fTimeResolution = GetDouble("TimeResolution", 1.0E-10);
-  // import input array
+  // read time resolution formula in seconds
+  fResolutionFormula->Compile(GetString("TimeResolution", "30e-12"));
 
+  // import track input array
   fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
   fItInputArray = fInputArray->MakeIterator();
 
+
   // create output array
-
-  fOutputArray = ExportArray(GetString("OutputArray", "muons"));
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
 }
 
@@ -92,5 +94,8 @@
 {
   Candidate *candidate, *mother;
-  Double_t ti, tf_smeared, tf;
+  Double_t tf_smeared, tf;
+  Double_t eta, energy;
+  Double_t timeResolution;
+
   const Double_t c_light = 2.99792458E8;
 
@@ -98,27 +103,25 @@
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
+
     const TLorentzVector &candidateFinalPosition = candidate->Position;
+    const TLorentzVector &candidateMomentum = candidate->Momentum;
 
-    ti = candidateInitialPosition.T() * 1.0E-3 / c_light;
     tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
 
+    eta = candidateMomentum.Eta();
+    energy = candidateMomentum.E();
+
     // apply smearing formula
-    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
-    ti = ti + tf_smeared - tf;
-    tf = tf_smeared;
+    timeResolution = fResolutionFormula->Eval(0.0, eta, 0.0, energy);
+    tf_smeared = gRandom->Gaus(tf, timeResolution);
 
     mother = candidate;
     candidate = static_cast<Candidate *>(candidate->Clone());
-    candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
-    candidate->Position.SetT(tf * 1.0E3 * c_light);
 
-    candidate->ErrorT = fTimeResolution * 1.0E3 * c_light;
+    candidate->Position.SetT(tf_smeared * 1.0E3 * c_light);
+    candidate->ErrorT = timeResolution * 1.0E3 * c_light;
 
     candidate->AddCandidate(mother);
-
     fOutputArray->Add(candidate);
   }
 }
-
-//------------------------------------------------------------------------------
Index: modules/TimeSmearing.h
===================================================================
--- modules/TimeSmearing.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TimeSmearing.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -22,7 +22,7 @@
 /** \class TimeSmearing
  *
- *  Performs transverse time smearing.
+ *  Performs time smearing.
  *
- *  \author Michele Selvaggi - UCL, Louvain-la-Neuve
+ *  \author Michele Selvaggi - CERN
  *
  */
@@ -32,4 +32,5 @@
 class TIterator;
 class TObjArray;
+class DelphesFormula;
 
 class TimeSmearing: public DelphesModule
@@ -44,5 +45,7 @@
 
 private:
-  Double_t fTimeResolution;
+
+  DelphesFormula *fResolutionFormula;
+  Int_t fVertexTimeMode;
 
   TIterator *fItInputArray; //!
Index: modules/TrackCovariance.cc
===================================================================
--- modules/TrackCovariance.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TrackCovariance.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,188 @@
+  /*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), 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 TrackCovariance
+ *
+ *  Smears track parameters according to appropriate covariance matrix.
+ *
+ *  \authors F. Bedeschi - INFN Pisa
+*            P. Demin - UCLouvain, Louvain-la-Neuve
+ *           M. Selvaggi - CERN
+ *
+ *
+ */
+
+//FIXME add reference to Bedeschi-code
+//FIXME make sure about units of P, X
+//FIXME fix pt > 200 GeV issue and angle > 6.41
+
+#include "modules/TrackCovariance.h"
+
+#include "classes/DelphesClasses.h"
+
+#include "TrackCovariance/SolGeom.h"
+#include "TrackCovariance/SolGridCov.h"
+#include "TrackCovariance/ObsTrk.h"
+
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TObjArray.h"
+
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+TrackCovariance::TrackCovariance() :
+  fGeometry(0), fCovariance(0), fAcx(0), fItInputArray(0)
+{
+  fGeometry = new SolGeom();
+  fCovariance = new SolGridCov();
+}
+
+//------------------------------------------------------------------------------
+
+TrackCovariance::~TrackCovariance()
+{
+  if(fGeometry) delete fGeometry;
+  if(fCovariance) delete fCovariance;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Init()
+{
+  fBz = GetDouble("Bz", 0.0);
+  fGeometry->Read(GetString("DetectorGeometry", ""));
+  fNMinHits = GetInt("NMinHits", 6);
+
+  // load geometry
+  fCovariance->Calc(fGeometry);
+  fCovariance->SetMinHits(fNMinHits);
+  // load geometry
+  fAcx = fCovariance->AccPnt();
+
+  // import input array
+  fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Process()
+{
+  Candidate *candidate, *mother, *particle;
+  Double_t mass, p, pt, q, ct;
+  Double_t dd0, ddz, dphi, dct, dp, dpt, dC;
+  //
+  // Get cylindrical box for fast track simulation
+  //
+  Double_t Rin = fGeometry->GetRmin();
+  Double_t ZinPos = fGeometry->GetZminPos();
+  Double_t ZinNeg = fGeometry->GetZminNeg();
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+
+    // converting to meters
+    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+
+    // converting to meters
+    const TLorentzVector &candidatePosition = particle->Position*1e-03;
+    const TLorentzVector &candidateMomentum = particle->Momentum;
+
+    Bool_t inside = TrkUtil::IsInside(candidatePosition.Vect(), Rin, ZinNeg, ZinPos); // Check if in inner box
+    Bool_t Accept = kTRUE;
+    if(inside) Accept = fCovariance->IsAccepted(candidateMomentum.Vect());
+    else       Accept = fCovariance->IsAccepted(candidatePosition.Vect(),candidateMomentum.Vect(), fGeometry);
+    if(!Accept) continue;
+
+    mass = candidateMomentum.M();
+
+    ObsTrk track(candidatePosition.Vect(), candidateMomentum.Vect(), candidate->Charge, fCovariance, fGeometry);
+
+    mother    = candidate;
+    candidate = static_cast<Candidate*>(candidate->Clone());
+
+    candidate->Momentum.SetVectM(track.GetObsP(), mass);
+
+    // converting back to mm
+    candidate->InitialPosition.SetXYZT(track.GetObsX().X()*1e03,track.GetObsX().Y()*1e03,track.GetObsX().Z()*1e03,candidatePosition.T()*1e03);
+
+    // save full covariance 5x5 matrix internally (D0, phi, Curvature, dz, ctg(theta))
+    candidate->TrackCovariance = track.GetCov();
+
+    pt = candidate->Momentum.Pt();
+    p  = candidate->Momentum.P();
+    q  = track.GetObsQ();
+    ct = track.GetObsPar()[4];
+
+    candidate->Xd = track.GetObsX().X()*1e03;
+    candidate->Yd = track.GetObsX().Y()*1e03;
+    candidate->Zd = track.GetObsX().Z()*1e03;
+
+    candidate->D0       = track.GetObsPar()[0]*1e03;
+    candidate->Phi      = track.GetObsPar()[1];
+
+    // inverse of curvature
+    candidate->C        = track.GetObsPar()[2]*1e-03;
+    candidate->DZ       = track.GetObsPar()[3]*1e03;
+    candidate->CtgTheta = track.GetObsPar()[4];
+    candidate->P        = track.GetObsP().Mag();
+    candidate->PT       = pt;
+    candidate->Charge   = q;
+
+    dd0       = TMath::Sqrt(track.GetCov()(0, 0))*1e03;
+    ddz       = TMath::Sqrt(track.GetCov()(3, 3))*1e03;
+    dphi      = TMath::Sqrt(track.GetCov()(1, 1));
+    dct       = TMath::Sqrt(track.GetCov()(4, 4));
+    dpt       = 2 * TMath::Sqrt( track.GetCov()(2, 2))*pt*pt / (0.2998*fBz);
+    dp        = TMath::Sqrt((1.+ct*ct)*dpt*dpt + 4*pt*pt*ct*ct*dct*dct/(1.+ct*ct)/(1.+ct*ct));
+    dC        = TMath::Sqrt(track.GetCov()(2, 2))*1e-03;
+
+    candidate->ErrorD0 = dd0;
+    candidate->ErrorDZ = ddz;
+    candidate->ErrorP = dp;
+    candidate->ErrorC = dC;
+    candidate->ErrorCtgTheta = dct;
+    candidate->ErrorPhi = dphi;
+    candidate->ErrorPT = dpt;
+    //candidate->TrackResolution = dpt / pt;
+    candidate->TrackResolution = dp / p;
+
+    candidate->AddCandidate(mother);
+
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/TrackCovariance.h
===================================================================
--- modules/TrackCovariance.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TrackCovariance.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,68 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), 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 TrackCovariance_h
+#define TrackCovariance_h
+
+/** \class TrackCovariance
+ *
+ *  Smears track parameters according to appropriate covariance matrix.
+ *
+ *  \authors P. Demin - UCLouvain, Louvain-la-Neuve
+ *           M. Selvaggi - CERN
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+
+class SolGeom;
+class SolGridCov;
+class AcceptanceClx;
+
+class TrackCovariance: public DelphesModule
+{
+public:
+  TrackCovariance();
+  ~TrackCovariance();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+  Double_t fBz;
+  Int_t fNMinHits;
+
+  SolGeom *fGeometry;
+  SolGridCov *fCovariance;
+
+  AcceptanceClx *fAcx;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(TrackCovariance, 1)
+};
+
+#endif
Index: modules/TrackSmearing.cc
===================================================================
--- modules/TrackSmearing.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TrackSmearing.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -158,5 +158,5 @@
   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 pt, eta, e, m, 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;
   Double_t q, r;
@@ -223,4 +223,6 @@
     pt = momentum.Pt();
     eta = momentum.Eta();
+    e = momentum.E();
+    m = momentum.M();
 
     d0 = trueD0 = candidate->D0;
@@ -232,5 +234,5 @@
 
     if(fUseD0Formula)
-      d0Error = fD0Formula->Eval(pt, eta);
+      d0Error = fD0Formula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -247,5 +249,5 @@
 
     if(fUseDZFormula)
-      dzError = fDZFormula->Eval(pt, eta);
+      dzError = fDZFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -262,5 +264,5 @@
 
     if(fUsePFormula)
-      pError = fPFormula->Eval(pt, eta) * p;
+      pError = fPFormula->Eval(pt, eta, phi, e, candidate) * p;
     else
     {
@@ -277,5 +279,5 @@
 
     if(fUseCtgThetaFormula)
-      ctgThetaError = fCtgThetaFormula->Eval(pt, eta);
+      ctgThetaError = fCtgThetaFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -292,5 +294,5 @@
 
     if(fUsePhiFormula)
-      phiError = fPhiFormula->Eval(pt, eta);
+      phiError = fPhiFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -331,5 +333,5 @@
     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->Momentum.SetE(TMath::Sqrt(p*p + m*m));
     candidate->PT = candidate->Momentum.Pt();
 
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TreeWriter.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -72,4 +72,5 @@
   fClassMap[Track::Class()] = &TreeWriter::ProcessTracks;
   fClassMap[Tower::Class()] = &TreeWriter::ProcessTowers;
+  fClassMap[ParticleFlowCandidate::Class()] = &TreeWriter::ProcessParticleFlowCandidates;
   fClassMap[Photon::Class()] = &TreeWriter::ProcessPhotons;
   fClassMap[Electron::Class()] = &TreeWriter::ProcessElectrons;
@@ -123,4 +124,17 @@
     fBranchMap.insert(make_pair(branch, make_pair(itClassMap->second, array)));
   }
+
+  param = GetParam("Info");
+  TString infoName;
+  Double_t infoValue;
+
+  size = param.GetSize();
+  for(i = 0; i < size / 2; ++i)
+  {
+    infoName = param[i * 2].GetString();
+    infoValue = param[i * 2 + 1].GetDouble();
+
+    AddInfo(infoName, infoValue);
+  }
 }
 
@@ -138,4 +152,5 @@
   it1.Reset();
   array->Clear();
+
   while((candidate = static_cast<Candidate *>(it1.Next())))
   {
@@ -215,11 +230,4 @@
     entry->Pz = momentum.Pz();
 
-    entry->D0 = candidate->D0;
-    entry->DZ = candidate->DZ;
-    entry->P = momentum.P();
-    entry->PT = candidate->PT;
-    entry->CtgTheta = candidate->CtgTheta;
-    entry->Phi = candidate->Phi;
-
     entry->Eta = eta;
     entry->Phi = momentum.Phi();
@@ -322,5 +330,5 @@
   Candidate *particle = 0;
   Track *entry = 0;
-  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
+  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
   const Double_t c_light = 2.99792458E8;
 
@@ -356,12 +364,29 @@
 
     entry->D0 = candidate->D0;
-    entry->ErrorD0 = candidate->ErrorD0;
     entry->DZ = candidate->DZ;
-    entry->ErrorDZ = candidate->ErrorDZ;
+    entry->Nclusters = candidate->Nclusters;
+    entry->dNdx = candidate->dNdx;
 
     entry->ErrorP = candidate->ErrorP;
     entry->ErrorPT = candidate->ErrorPT;
+
+    // diagonal covariance matrix terms
+    entry->ErrorD0 = candidate->ErrorD0;
+    entry->ErrorC = candidate->ErrorC;
+    entry->ErrorPhi = candidate->ErrorPhi;
+    entry->ErrorDZ = candidate->ErrorDZ;
     entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
-    entry->ErrorPhi = candidate->ErrorPhi;
+
+    // add some offdiagonal covariance matrix elements
+    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
+    entry->ErrorD0C            = candidate->TrackCovariance(0,2);
+    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
+    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
+    entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
+    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
+    entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
+    entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
+    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
+    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
 
     entry->Xd = candidate->Xd;
@@ -374,4 +399,5 @@
     p = momentum.P();
     phi = momentum.Phi();
+    m = momentum.M();
     ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
 
@@ -386,7 +412,10 @@
     entry->Phi = phi;
     entry->CtgTheta = ctgTheta;
+    entry->C = candidate->C;
+    entry->Mass = m;
 
     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
-    const TLorentzVector &initialPosition = particle->Position;
+    //const TLorentzVector &initialPosition = particle->Position;
+    const TLorentzVector &initialPosition = candidate->InitialPosition;
 
     entry->X = initialPosition.X();
@@ -394,4 +423,5 @@
     entry->Z = initialPosition.Z();
     entry->T = initialPosition.T() * 1.0E-3 / c_light;
+    entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light;
 
     entry->Particle = particle;
@@ -435,4 +465,5 @@
     entry->Eem = candidate->Eem;
     entry->Ehad = candidate->Ehad;
+    entry->Etrk = candidate->Etrk;
     entry->Edges[0] = candidate->Edges[0];
     entry->Edges[1] = candidate->Edges[1];
@@ -444,4 +475,126 @@
 
     FillParticles(candidate, &entry->Particles);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void TreeWriter::ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array)
+{
+
+  TIter iterator(array);
+  Candidate *candidate = 0;
+  Candidate *particle = 0;
+  ParticleFlowCandidate *entry = 0;
+  Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
+  const Double_t c_light = 2.99792458E8;
+
+  // loop over all tracks
+  iterator.Reset();
+  while((candidate = static_cast<Candidate *>(iterator.Next())))
+  {
+    const TLorentzVector &position = candidate->Position;
+
+    cosTheta = TMath::Abs(position.CosTheta());
+    signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
+    eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
+    rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
+
+    entry = static_cast<ParticleFlowCandidate *>(branch->NewEntry());
+
+    entry->SetBit(kIsReferenced);
+    entry->SetUniqueID(candidate->GetUniqueID());
+
+    entry->PID = candidate->PID;
+
+    entry->Charge = candidate->Charge;
+
+    entry->EtaOuter = eta;
+    entry->PhiOuter = position.Phi();
+
+    entry->XOuter = position.X();
+    entry->YOuter = position.Y();
+    entry->ZOuter = position.Z();
+    entry->TOuter = position.T() * 1.0E-3 / c_light;
+
+    entry->L = candidate->L;
+
+    entry->D0 = candidate->D0;
+    entry->DZ = candidate->DZ;
+    entry->Nclusters = candidate->Nclusters;
+    entry->dNdx = candidate->dNdx;
+
+    entry->ErrorP = candidate->ErrorP;
+    entry->ErrorPT = candidate->ErrorPT;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+
+
+    // diagonal covariance matrix terms
+
+    entry->ErrorD0 = candidate->ErrorD0;
+    entry->ErrorC = candidate->ErrorC;
+    entry->ErrorPhi = candidate->ErrorPhi;
+    entry->ErrorDZ = candidate->ErrorDZ;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+
+    // add some offdiagonal covariance matrix elements
+    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
+    entry->ErrorD0C            = candidate->TrackCovariance(0,2);
+    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
+    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
+    entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
+    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
+    entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
+    entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
+    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
+    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
+
+    entry->Xd = candidate->Xd;
+    entry->Yd = candidate->Yd;
+    entry->Zd = candidate->Zd;
+
+    const TLorentzVector &momentum = candidate->Momentum;
+
+    e = momentum.E();
+    pt = momentum.Pt();
+    p = momentum.P();
+    phi = momentum.Phi();
+    m = momentum.M();
+    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
+
+    entry->E = e;
+    entry->P = p;
+    entry->PT = pt;
+    entry->Eta = eta;
+    entry->Phi = phi;
+    entry->CtgTheta = ctgTheta;
+    entry->C = candidate->C;
+    entry->Mass = m;
+
+    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+    //const TLorentzVector &initialPosition = particle->Position;
+    const TLorentzVector &initialPosition = candidate->InitialPosition;
+
+    entry->X = initialPosition.X();
+    entry->Y = initialPosition.Y();
+    entry->Z = initialPosition.Z();
+    entry->T = initialPosition.T() * 1.0E-3 / c_light;
+    entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light;
+
+    entry->VertexIndex = candidate->ClusterIndex;
+
+    entry->Eem = candidate->Eem;
+    entry->Ehad = candidate->Ehad;
+    entry->Etrk = candidate->Etrk;
+    entry->Edges[0] = candidate->Edges[0];
+    entry->Edges[1] = candidate->Edges[1];
+    entry->Edges[2] = candidate->Edges[2];
+    entry->Edges[3] = candidate->Edges[3];
+
+    //entry->T = position.T() * 1.0E-3 / c_light;
+    entry->NTimeHits = candidate->NTimeHits;
+
+    FillParticles(candidate, &entry->Particles);
+
   }
 }
@@ -687,4 +840,7 @@
     entry->NCharged = candidate->NCharged;
     entry->NNeutrals = candidate->NNeutrals;
+
+    entry->NeutralEnergyFraction = candidate->NeutralEnergyFraction;
+    entry->ChargedEnergyFraction = candidate->ChargedEnergyFraction;
     entry->Beta = candidate->Beta;
     entry->BetaStar = candidate->BetaStar;
Index: modules/TreeWriter.h
===================================================================
--- modules/TreeWriter.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ modules/TreeWriter.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -56,4 +56,5 @@
   void ProcessTracks(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessTowers(ExRootTreeBranch *branch, TObjArray *array);
+  void ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessPhotons(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array);
@@ -79,5 +80,5 @@
 #endif
 
-  ClassDef(TreeWriter, 1)
+  ClassDef(TreeWriter, 2)
 };
 
Index: modules/TruthVertexFinder.cc
===================================================================
--- modules/TruthVertexFinder.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TruthVertexFinder.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,158 @@
+/*
+ *  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 TruthVertexFinder
+ *
+ *  Merges particles from pile-up sample into event
+ *
+ *  \author M. Selvaggi - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/TruthVertexFinder.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesPileUpReader.h"
+#include "classes/DelphesTF2.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;
+
+//------------------------------------------------------------------------------
+
+TruthVertexFinder::TruthVertexFinder() :
+fItInputArray(0), fItOutputArray(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+TruthVertexFinder::~TruthVertexFinder()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void TruthVertexFinder::Init()
+{
+
+  fResolution = GetDouble("Resolution", 1E-06); // resolution in meters
+  // import input array
+  fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output arrays
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+  //fItOutputArray = fVertexOutputArray->MakeIterator();
+}
+
+//------------------------------------------------------------------------------
+
+void TruthVertexFinder::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+  if(fItOutputArray) delete fItOutputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TruthVertexFinder::Process()
+{
+  Int_t nvtx = -1;
+  Float_t pt;
+  Candidate *candidate, *vertex;
+  DelphesFactory *factory;
+
+
+  fItInputArray->Reset();
+
+  factory = GetFactory();
+  vertex = factory->NewCandidate();
+
+  TLorentzVector vertexPosition(0., 0., 0., 0.);
+
+  nvtx=0;
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+
+     const TLorentzVector &candidatePosition = candidate->Position;
+     const TLorentzVector &candidateMomentum = candidate->Momentum;
+
+     pt = candidateMomentum.Pt();
+     
+     // check whether vertex already included, if so add particle
+     Bool_t old_vertex=false;
+     fItOutputArray = fVertexOutputArray->MakeIterator();
+     fItOutputArray->Reset();
+     while((vertex = static_cast<Candidate *>(fItOutputArray->Next())))
+     {
+        const TLorentzVector &vertexPosition = vertex->Position;
+        // check whether spatial difference is < 1 um, in that case assume it is the same vertex
+        if ( TMath::Abs((candidatePosition.P() - vertexPosition.P())) < fResolution*1.E3)
+        {
+           old_vertex=true;
+           vertex->AddCandidate(candidate);
+           if (TMath::Abs(candidate->Charge) > 0)
+           {
+              vertex->ClusterNDF += 1;
+              vertex->GenSumPT2 += pt*pt;
+           }
+        }
+     }
+
+     // else fill new vertex
+     if (!old_vertex)
+     {
+        vertex = factory->NewCandidate();
+        vertex->Position = candidatePosition;
+        vertex->ClusterIndex = nvtx;
+
+        if (TMath::Abs(candidate->Charge) > 0)
+        {
+           vertex->ClusterNDF = 1;
+           vertex->GenSumPT2 = pt*pt;
+        }
+        else
+        {
+           vertex->ClusterNDF = 0;
+           vertex->GenSumPT2 = 0.;
+        }
+        fVertexOutputArray->Add(vertex);
+        nvtx++;
+      }
+   }
+}
+
+//------------------------------------------------------------------------------
Index: modules/TruthVertexFinder.h
===================================================================
--- modules/TruthVertexFinder.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ modules/TruthVertexFinder.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,57 @@
+/*
+ *  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 TruthVertexFinder_h
+#define TruthVertexFinder_h
+
+/** \class TruthVertexFinder
+ *
+ *  Produces list of MC truth vertices
+ *
+ *  \author M. Selvaggi - CERN
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TObjArray;
+
+class TruthVertexFinder: public DelphesModule
+{
+public:
+  TruthVertexFinder();
+  ~TruthVertexFinder();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  Double_t fResolution; //!
+
+  TIterator *fItInputArray; //!
+  TIterator *fItOutputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fVertexOutputArray; //!
+  ClassDef(TruthVertexFinder, 1)
+};
+
+#endif
