Index: classes/DelphesClasses.h
===================================================================
--- classes/DelphesClasses.h	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ classes/DelphesClasses.h	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -385,5 +385,5 @@
 
   Float_t NeutralEnergyFraction;  // charged energy fraction
-  Float_t ChargedEnergyFraction;  // neutral energy fraction 
+  Float_t ChargedEnergyFraction;  // neutral energy fraction
 
   Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
@@ -439,4 +439,5 @@
   Float_t CtgTheta; // track cotangent of theta
   Float_t C; // track curvature inverse
+  Float_t Mass; // particle mass
 
   Float_t EtaOuter; // track pseudorapidity at the tracker edge
@@ -472,14 +473,14 @@
 
   // track covariance off-diagonal terms
-  Float_t ErrorD0Phi;	  
-  Float_t ErrorD0C;	  
-  Float_t ErrorD0DZ;	  
-  Float_t ErrorD0CtgTheta;   
-  Float_t ErrorPhiC;	  
-  Float_t ErrorPhiDZ;	  
-  Float_t ErrorPhiCtgTheta ; 
-  Float_t ErrorCDZ;	  
-  Float_t ErrorCCtgTheta;    
-  Float_t ErrorDZCtgTheta;   
+  Float_t ErrorD0Phi;
+  Float_t ErrorD0C;
+  Float_t ErrorD0DZ;
+  Float_t ErrorD0CtgTheta;
+  Float_t ErrorPhiC;
+  Float_t ErrorPhiDZ;
+  Float_t ErrorPhiCtgTheta ;
+  Float_t ErrorCDZ;
+  Float_t ErrorCCtgTheta;
+  Float_t ErrorDZCtgTheta;
 
   TRef Particle; // reference to generated particle
@@ -542,4 +543,5 @@
   Float_t CtgTheta; // track cotangent of theta
   Float_t C; // track curvature inverse
+  Float_t Mass; // particle mass
 
   Float_t EtaOuter; // track pseudorapidity at the tracker edge
@@ -575,14 +577,14 @@
 
   // track covariance off-diagonal terms
-  Float_t ErrorD0Phi;	  
-  Float_t ErrorD0C;	  
-  Float_t ErrorD0DZ;	  
-  Float_t ErrorD0CtgTheta;   
-  Float_t ErrorPhiC;	  
-  Float_t ErrorPhiDZ;	  
-  Float_t ErrorPhiCtgTheta ; 
-  Float_t ErrorCDZ;	  
-  Float_t ErrorCCtgTheta;    
-  Float_t ErrorDZCtgTheta;   
+  Float_t ErrorD0Phi;
+  Float_t ErrorD0C;
+  Float_t ErrorD0DZ;
+  Float_t ErrorD0CtgTheta;
+  Float_t ErrorPhiC;
+  Float_t ErrorPhiDZ;
+  Float_t ErrorPhiCtgTheta ;
+  Float_t ErrorCDZ;
+  Float_t ErrorCCtgTheta;
+  Float_t ErrorDZCtgTheta;
 
   Int_t VertexIndex; // reference to vertex
@@ -710,5 +712,5 @@
   Float_t FracPt[5];
   Float_t NeutralEnergyFraction;  // charged energy fraction
-  Float_t ChargedEnergyFraction;  // neutral energy fraction 
+  Float_t ChargedEnergyFraction;  // neutral energy fraction
 
 
@@ -767,5 +769,5 @@
   // event characteristics variables
   Double_t ParticleDensity; // particle multiplicity density in the proximity of the particle
-  
+
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
Index: modules/AngularSmearing.cc
===================================================================
--- modules/AngularSmearing.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/AngularSmearing.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -98,18 +98,17 @@
 {
   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, candidate));
     phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e, candidate));
@@ -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 b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/Calorimeter.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -559,5 +559,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 +646,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/DenseTrackFilter.cc
===================================================================
--- modules/DenseTrackFilter.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/DenseTrackFilter.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -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 b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/DualReadoutCalorimeter.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -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;
- 
+
 }
 
@@ -379,9 +379,9 @@
       fHCalTrackEnergy = 0.0;
       fTrackEnergy = 0.0;
-      
+
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
       fTrackSigma = 0.0;
-      
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -390,5 +390,5 @@
       fHCalTowerTrackArray->Clear();
       fTowerTrackArray->Clear();
-    
+
     }
 
@@ -414,40 +414,10 @@
       }
 
-      
-      /*
-      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)
@@ -455,13 +425,13 @@
         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
@@ -511,10 +481,10 @@
   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;
@@ -536,5 +506,5 @@
   //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)
   {
@@ -544,5 +514,5 @@
   }
 
-  // if hadronic fraction > 0, use HCAL resolution 
+  // if hadronic fraction > 0, use HCAL resolution
   else
   {
@@ -554,5 +524,5 @@
   energy = LogNormal(energy, sigma);
   //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
-  
+
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
 
@@ -626,5 +596,5 @@
 
   // fill energy flow candidates
-  
+
   fTrackSigma = TMath::Sqrt(fTrackSigma);
   neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
@@ -638,5 +608,5 @@
   if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
   {
-    
+
     //cout<<"significant neutral excess found:"<<endl;
     // create new photon tower
@@ -646,5 +616,5 @@
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
 
-    // if no hadronic energy, use ECAL resolution 
+    // if no hadronic energy, use ECAL resolution
     if (fHCalTowerEnergy <= fHCalEnergyMin)
     {
@@ -654,5 +624,5 @@
     }
 
-    // if hadronic fraction > 0, use HCAL resolution 
+    // if hadronic fraction > 0, use HCAL resolution
     else
     {
@@ -663,5 +633,5 @@
 
     fEFlowPhotonOutputArray->Add(tower);
-    
+
 
     //clone tracks
@@ -669,5 +639,4 @@
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
-      //cout<<"looping over tracks"<<endl;
       mother = track;
       track = static_cast<Candidate*>(track->Clone());
@@ -677,5 +646,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)
@@ -685,5 +654,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;
 
@@ -691,15 +660,13 @@
     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);
     }
   }
-  
+
 
 }
Index: modules/EnergySmearing.cc
===================================================================
--- modules/EnergySmearing.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/EnergySmearing.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -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/MomentumSmearing.cc
===================================================================
--- modules/MomentumSmearing.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/MomentumSmearing.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -95,5 +95,5 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e, res;
+  Double_t pt, eta, phi, e, m, res;
 
   fItInputArray->Reset();
@@ -106,4 +106,5 @@
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
+    m = candidateMomentum.M();
     res = fFormula->Eval(pt, eta, phi, e, candidate);
 
@@ -121,5 +122,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/SimpleCalorimeter.cc
===================================================================
--- modules/SimpleCalorimeter.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/SimpleCalorimeter.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -507,7 +507,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/TrackCovariance.cc
===================================================================
--- modules/TrackCovariance.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/TrackCovariance.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -1,3 +1,3 @@
-/*
+  /*
  *  Delphes: a framework for fast simulation of a generic collider experiment
  *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), Belgium
Index: modules/TrackSmearing.cc
===================================================================
--- modules/TrackSmearing.cc	(revision b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/TrackSmearing.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -158,5 +158,5 @@
   TLorentzVector beamSpotPosition;
   Candidate *candidate, *mother;
-  Double_t pt, eta, e, 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;
@@ -224,5 +224,6 @@
     eta = momentum.Eta();
     e = momentum.E();
-    
+    m = momentum.M();
+
     d0 = trueD0 = candidate->D0;
     dz = trueDZ = candidate->DZ;
@@ -332,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 b8a6aa396f5ef50a79bdd687a9c1bab9d6813b0b)
+++ modules/TreeWriter.cc	(revision fd4b326a0bceca9ef6bf17e978c87dedd995cf91)
@@ -320,5 +320,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;
 
@@ -387,4 +387,5 @@
     p = momentum.P();
     phi = momentum.Phi();
+    m = momentum.M();
     ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
 
@@ -400,4 +401,5 @@
     entry->CtgTheta = ctgTheta;
     entry->C = candidate->C;
+    entry->Mass = m;
 
     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
@@ -470,5 +472,5 @@
   Candidate *particle = 0;
   ParticleFlowCandidate *entry = 0;
-  Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
+  Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
   const Double_t c_light = 2.99792458E8;
 
@@ -541,4 +543,5 @@
     p = momentum.P();
     phi = momentum.Phi();
+    m = momentum.M();
     ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
 
@@ -550,4 +553,5 @@
     entry->CtgTheta = ctgTheta;
     entry->C = candidate->C;
+    entry->Mass = m;
 
     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
