Index: classes/DelphesCscClusterFormula.cc
===================================================================
--- classes/DelphesCscClusterFormula.cc	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
+++ classes/DelphesCscClusterFormula.cc	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
@@ -0,0 +1,83 @@
+/*
+ *  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/>.
+ */
+
+#include "classes/DelphesCscClusterFormula.h"
+#include "classes/DelphesClasses.h"
+
+#include "TString.h"
+
+#include <stdexcept>
+#include <iostream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+DelphesCscClusterFormula::DelphesCscClusterFormula() :
+  TFormula()
+{
+}
+
+//------------------------------------------------------------------------------
+
+DelphesCscClusterFormula::DelphesCscClusterFormula(const char *name, const char *expression) :
+  TFormula()
+{
+}
+
+//------------------------------------------------------------------------------
+
+DelphesCscClusterFormula::~DelphesCscClusterFormula()
+{
+}
+
+//------------------------------------------------------------------------------
+
+Int_t DelphesCscClusterFormula::Compile(const char *expression)
+{
+  TString buffer;
+  const char *it;
+  for(it = expression; *it; ++it)
+  {
+    if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\') continue;
+    buffer.Append(*it);
+  }
+  buffer.ReplaceAll("decayR", "x");
+  buffer.ReplaceAll("decayZ", "y");
+  buffer.ReplaceAll("Ehad", "z");
+  buffer.ReplaceAll("Eem", "t");
+
+#if ROOT_VERSION_CODE < ROOT_VERSION(6, 3, 0)
+  TFormula::SetMaxima(100000, 1000, 1000000);
+#endif
+  if(TFormula::Compile(buffer) != 0)
+  {
+    throw runtime_error("Invalid formula.");
+  }
+  return 0;
+}
+
+//------------------------------------------------------------------------------
+
+Double_t DelphesCscClusterFormula::Eval(Double_t decayR, Double_t decayZ, Double_t Ehad, Double_t Eem)
+{
+  Double_t x[4] = {decayR, decayZ, Ehad, Eem};
+  return EvalPar(x);
+}
+
+//------------------------------------------------------------------------------
Index: classes/DelphesCscClusterFormula.h
===================================================================
--- classes/DelphesCscClusterFormula.h	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
+++ classes/DelphesCscClusterFormula.h	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
@@ -0,0 +1,40 @@
+/*
+ *  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 DelphesCscClusterFormula_h
+#define DelphesCscClusterFormula_h
+
+#include "TFormula.h"
+
+class Candidate;
+
+class DelphesCscClusterFormula: public TFormula
+{
+public:
+  DelphesCscClusterFormula();
+
+  DelphesCscClusterFormula(const char *name, const char *expression);
+
+  ~DelphesCscClusterFormula();
+
+  Int_t Compile(const char *expression);
+
+  Double_t Eval(Double_t deltaR, Double_t deltaZ = 0, Double_t Ehad = 0, Double_t Eem = 0);
+};
+
+#endif /* DelphesCscClusterFormula_h */
Index: classes/DelphesHepMC2Reader.cc
===================================================================
--- classes/DelphesHepMC2Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC2Reader.cc	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
@@ -438,4 +438,5 @@
 {
   Candidate *candidate;
+  Candidate *candidateDaughter;
   map<int, pair<int, int> >::iterator itMotherMap;
   map<int, pair<int, int> >::iterator itDaughterMap;
@@ -446,4 +447,5 @@
     candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
 
+
     if(candidate->M1 > 0)
     {
@@ -477,4 +479,6 @@
         candidate->D1 = -1;
         candidate->D2 = -1;
+        const TLorentzVector &decayPosition = candidate->Position;
+        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
       }
       else
@@ -482,8 +486,12 @@
         candidate->D1 = itDaughterMap->second.first;
         candidate->D2 = itDaughterMap->second.second;
-      }
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
+        candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1));
+        const TLorentzVector &decayPosition = candidateDaughter->Position;
+        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
+
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
Index: classes/DelphesHepMC3Reader.cc
===================================================================
--- classes/DelphesHepMC3Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC3Reader.cc	(revision 29b722a3a3a43304999897845a8762c4d304ea09)
@@ -458,4 +458,5 @@
   TObjArray *array;
   Candidate *candidate;
+  Candidate *candidateDaughter;
   TParticlePDG *pdgParticle;
   int pdgCode;
@@ -573,4 +574,6 @@
         candidate->D1 = -1;
         candidate->D2 = -1;
+        const TLorentzVector &decayPosition = candidate->Position;
+        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
       }
       else
@@ -578,8 +581,11 @@
         candidate->D1 = itDaughterMap->second.first;
         candidate->D2 = itDaughterMap->second.second;
-      }
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
+        candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1));
+        const TLorentzVector &decayPosition = candidateDaughter->Position;
+        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
