Index: classes/CMakeLists.txt
===================================================================
--- classes/CMakeLists.txt	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/CMakeLists.txt	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -9,11 +9,23 @@
 list(REMOVE_ITEM headers ${CMAKE_CURRENT_SOURCE_DIR}/ClassesLinkDef.h)
 
-DELPHES_GENERATE_DICTIONARY(ClassesDict 
+# the macro invocation for ROOT6 ensures that the headers are relocatable
+if (NOT ${ROOT_VERSION} VERSION_LESS "6.0.0")
+  DELPHES_GENERATE_DICTIONARY(ClassesDict 
+    classes/DelphesModule.h
+    classes/DelphesFactory.h
+    classes/SortableObject.h
+    classes/DelphesClasses.h
+    LINKDEF ClassesLinkDef.h
+  )
+else()
+  # for ROOT5 the above fails, keep the following as workaround
+  DELPHES_GENERATE_DICTIONARY(ClassesDict
   ${CMAKE_CURRENT_SOURCE_DIR}/DelphesModule.h
   ${CMAKE_CURRENT_SOURCE_DIR}/DelphesFactory.h
   ${CMAKE_CURRENT_SOURCE_DIR}/SortableObject.h
   ${CMAKE_CURRENT_SOURCE_DIR}/DelphesClasses.h
-  LINKDEF ClassesLinkDef.h
-)
+    LINKDEF ClassesLinkDef.h
+  )
+endif()
 
 add_library(classes OBJECT ${sources} ClassesDict.cxx)
Index: classes/ClassesLinkDef.h
===================================================================
--- classes/ClassesLinkDef.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/ClassesLinkDef.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -62,4 +62,5 @@
 #pragma link C++ class Track+;
 #pragma link C++ class Tower+;
+#pragma link C++ class ParticleFlowCandidate+;
 #pragma link C++ class HectorHit+;
 
Index: classes/DelphesClasses.cc
===================================================================
--- classes/DelphesClasses.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesClasses.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -28,5 +28,4 @@
 
 #include "classes/DelphesClasses.h"
-
 #include "classes/DelphesFactory.h"
 #include "classes/SortableObject.h"
@@ -40,4 +39,5 @@
 CompBase *Track::fgCompare = CompPT<Track>::Instance();
 CompBase *Tower::fgCompare = CompE<Tower>::Instance();
+CompBase *ParticleFlowCandidate::fgCompare = CompE<ParticleFlowCandidate>::Instance();
 CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
 CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
@@ -104,7 +104,48 @@
 {
   TLorentzVector vec;
-  vec.SetPtEtaPhiM(PT, Eta, Phi, 0.0);
-  return vec;
-}
+  vec.SetPtEtaPhiM(PT, Eta, Phi, Mass);
+  return vec;
+}
+
+//------------------------------------------------------------------------------
+
+TMatrixDSym Track::CovarianceMatrix() const
+{
+  TMatrixDSym Cv;
+  Cv.ResizeTo(5, 5);
+
+  // convert diagonal term to original units
+  Cv(0, 0)=TMath::Power(ErrorD0, 2.);
+  Cv(1, 1)=TMath::Power(ErrorPhi, 2.);
+  Cv(2, 2)=TMath::Power(ErrorC, 2.);
+  Cv(3, 3)=TMath::Power(ErrorDZ, 2.);
+  Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.);
+
+  // off diagonal terms
+  Cv(0, 1)=ErrorD0Phi;
+  Cv(0, 2)=ErrorD0C;
+  Cv(0, 3)=ErrorD0DZ;
+  Cv(0, 4)=ErrorD0CtgTheta;
+  Cv(1, 2)=ErrorPhiC;
+  Cv(1, 3)=ErrorPhiDZ;
+  Cv(1, 4)=ErrorPhiCtgTheta;
+  Cv(2, 3)=ErrorCDZ;
+  Cv(2, 4)=ErrorCCtgTheta;
+  Cv(3, 4)=ErrorDZCtgTheta;
+
+  Cv(1, 0)=Cv(0, 1);
+  Cv(2, 0)=Cv(0, 2);
+  Cv(3, 0)=Cv(0, 3);
+  Cv(4, 0)=Cv(0, 4);
+  Cv(2, 1)=Cv(1, 2);
+  Cv(3, 1)=Cv(1, 3);
+  Cv(4, 1)=Cv(1, 4);
+  Cv(3, 2)=Cv(2, 3);
+  Cv(4, 2)=Cv(2, 4);
+  Cv(4, 3)=Cv(3, 4);
+
+  return Cv;
+}
+
 
 //------------------------------------------------------------------------------
@@ -115,4 +156,53 @@
   vec.SetPtEtaPhiM(ET, Eta, Phi, 0.0);
   return vec;
+}
+
+//------------------------------------------------------------------------------
+
+TLorentzVector ParticleFlowCandidate::P4() const
+{
+  TLorentzVector vec;
+  vec.SetPtEtaPhiM(PT, Eta, Phi, Mass);
+  return vec;
+}
+
+//------------------------------------------------------------------------------
+
+TMatrixDSym ParticleFlowCandidate::CovarianceMatrix() const
+{
+  TMatrixDSym Cv;
+  Cv.ResizeTo(5, 5);
+
+  // convert diagonal term to original units
+  Cv(0, 0)=TMath::Power(ErrorD0, 2.);
+  Cv(1, 1)=TMath::Power(ErrorPhi, 2.);
+  Cv(2, 2)=TMath::Power(ErrorC, 2.);
+  Cv(3, 3)=TMath::Power(ErrorDZ, 2.);
+  Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.);
+
+  // off diagonal terms
+  Cv(0, 1)=ErrorD0Phi;
+  Cv(0, 2)=ErrorD0C;
+  Cv(0, 3)=ErrorD0DZ;
+  Cv(0, 4)=ErrorD0CtgTheta;
+  Cv(1, 2)=ErrorPhiC;
+  Cv(1, 3)=ErrorPhiDZ;
+  Cv(1, 4)=ErrorPhiCtgTheta;
+  Cv(2, 3)=ErrorCDZ;
+  Cv(2, 4)=ErrorCCtgTheta;
+  Cv(3, 4)=ErrorDZCtgTheta;
+
+  Cv(1, 0)=Cv(0, 1);
+  Cv(2, 0)=Cv(0, 2);
+  Cv(3, 0)=Cv(0, 3);
+  Cv(4, 0)=Cv(0, 4);
+  Cv(2, 1)=Cv(1, 2);
+  Cv(3, 1)=Cv(1, 3);
+  Cv(4, 1)=Cv(1, 4);
+  Cv(3, 2)=Cv(2, 3);
+  Cv(4, 2)=Cv(2, 4);
+  Cv(4, 3)=Cv(3, 4);
+
+  return Cv;
 }
 
@@ -125,5 +215,5 @@
   Flavor(0), FlavorAlgo(0), FlavorPhys(0),
   BTag(0), BTagAlgo(0), BTagPhys(0),
-  TauTag(0), TauWeight(0.0), Eem(0.0), Ehad(0.0),
+  TauTag(0), TauWeight(0.0), Eem(0.0), Ehad(0.0), Etrk(0.0),
   DeltaEta(0.0), DeltaPhi(0.0),
   Momentum(0.0, 0.0, 0.0, 0.0),
@@ -133,15 +223,21 @@
   PositionError(0.0, 0.0, 0.0, 0.0),
   Area(0.0, 0.0, 0.0, 0.0),
+  TrackCovariance(5),
   L(0),
   D0(0), ErrorD0(0),
   DZ(0), ErrorDZ(0),
   P(0), ErrorP(0),
+  C(0), ErrorC(0),
   PT(0), ErrorPT(0),
   CtgTheta(0), ErrorCtgTheta(0),
   Phi(0), ErrorPhi(0),
   Xd(0), Yd(0), Zd(0),
+  Nclusters(0.0),
+  dNdx(0.0),
   TrackResolution(0),
   NCharged(0),
   NNeutrals(0),
+  NeutralEnergyFraction(0),  // charged energy fraction
+  ChargedEnergyFraction(0),  // neutral energy fraction
   Beta(0),
   BetaStar(0),
@@ -163,4 +259,5 @@
   ExclYmerge45(0),
   ExclYmerge56(0),
+  ParticleDensity(0),
   fFactory(0),
   fArray(0)
@@ -284,4 +381,5 @@
   object.Eem = Eem;
   object.Ehad = Ehad;
+  object.Etrk = Etrk;
   object.Edges[0] = Edges[0];
   object.Edges[1] = Edges[1];
@@ -304,4 +402,6 @@
   object.P = P;
   object.ErrorP = ErrorP;
+  object.C = C;
+  object.ErrorC = ErrorC;
   object.PT = PT;
   object.ErrorPT = ErrorPT;
@@ -313,7 +413,11 @@
   object.Yd = Yd;
   object.Zd = Zd;
+  object.Nclusters = Nclusters;
+  object.dNdx = dNdx;
   object.TrackResolution = TrackResolution;
   object.NCharged = NCharged;
   object.NNeutrals = NNeutrals;
+  object.NeutralEnergyFraction = NeutralEnergyFraction;
+  object.ChargedEnergyFraction = ChargedEnergyFraction;
   object.Beta = Beta;
   object.BetaStar = BetaStar;
@@ -366,5 +470,5 @@
   object.SoftDroppedSubJet1 = SoftDroppedSubJet1;
   object.SoftDroppedSubJet2 = SoftDroppedSubJet2;
-
+  object.TrackCovariance = TrackCovariance;
   object.fFactory = fFactory;
   object.fArray = 0;
@@ -413,4 +517,5 @@
   Eem = 0.0;
   Ehad = 0.0;
+  Etrk = 0.0;
   Edges[0] = 0.0;
   Edges[1] = 0.0;
@@ -424,4 +529,5 @@
   DecayPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  TrackCovariance.Zero();
   L = 0.0;
   ErrorT = 0.0;
@@ -432,4 +538,6 @@
   P = 0.0;
   ErrorP = 0.0;
+  C = 0.0;
+  ErrorC = 0.0;
   PT = 0.0;
   ErrorPT = 0.0;
@@ -441,4 +549,6 @@
   Yd = 0.0;
   Zd = 0.0;
+  Nclusters = 0.0;
+  dNdx = 0.0;
   TrackResolution = 0.0;
   NCharged = 0;
Index: classes/DelphesClasses.h
===================================================================
--- classes/DelphesClasses.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesClasses.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -33,4 +33,5 @@
 
 #include "TLorentzVector.h"
+#include "TMatrixDSym.h"
 #include "TObject.h"
 #include "TRef.h"
@@ -150,10 +151,5 @@
   Float_t Eta; // particle pseudorapidity
   Float_t Phi; // particle azimuthal angle
-
   Float_t Rapidity; // particle rapidity
-  Float_t CtgTheta; // particle cotangent of theta
-
-  Float_t D0; // particle transverse impact parameter
-  Float_t DZ; // particle longitudinal impact parameter
 
   Float_t T; // particle vertex position (t component) | hepevt.vhep[number][3]
@@ -394,4 +390,8 @@
   Int_t NCharged; // number of charged constituents
   Int_t NNeutrals; // number of neutral constituents
+
+  Float_t NeutralEnergyFraction;  // charged energy fraction
+  Float_t ChargedEnergyFraction;  // neutral energy fraction
+
   Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
   Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents)
@@ -445,4 +445,6 @@
   Float_t Phi; // track azimuthal angle
   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
@@ -466,4 +468,6 @@
   Float_t D0; // track transverse impact parameter
   Float_t DZ; // track longitudinal impact parameter
+  Float_t Nclusters; // Number of ionization clusters
+  Float_t dNdx; // Number of ionization clusters
 
   Float_t ErrorP; // track momentum error
@@ -475,4 +479,17 @@
   Float_t ErrorD0; // track transverse impact parameter error
   Float_t ErrorDZ; // track longitudinal impact parameter error
+  Float_t ErrorC; // track curvature error
+
+  // 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;
 
   TRef Particle; // reference to generated particle
@@ -484,4 +501,5 @@
 
   TLorentzVector P4() const;
+  TMatrixDSym CovarianceMatrix() const;
 
   ClassDef(Track, 3)
@@ -504,4 +522,5 @@
   Float_t Eem; // calorimeter tower electromagnetic energy
   Float_t Ehad; // calorimeter tower hadronic energy
+  Float_t Etrk; // total charged energy hitting tower
 
   Float_t Edges[4]; // calorimeter tower edges
@@ -514,5 +533,91 @@
   TLorentzVector P4() const;
 
-  ClassDef(Tower, 2)
+  ClassDef(Tower, 3)
+};
+
+//---------------------------------------------------------------------------
+
+class ParticleFlowCandidate: public SortableObject
+{
+
+public:
+  Int_t PID; // HEP ID number
+
+  Int_t Charge; // track charge
+
+  Float_t E; // reconstructed energy [GeV]
+  Float_t P; // track momentum
+  Float_t PT; // track transverse momentum
+  Float_t Eta; // track pseudorapidity
+  Float_t Phi; // track azimuthal angle
+  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
+  Float_t PhiOuter; // track azimuthal angle at the tracker edge
+
+  Float_t T; // track vertex position (t component)
+  Float_t X; // track vertex position (x component)
+  Float_t Y; // track vertex position (y component)
+  Float_t Z; // track vertex position (z component)
+
+  Float_t TOuter; // track position (t component) at the tracker edge
+  Float_t XOuter; // track position (x component) at the tracker edge
+  Float_t YOuter; // track position (y component) at the tracker edge
+  Float_t ZOuter; // track position (z component) at the tracker edge
+
+  Float_t Xd; // X coordinate of point of closest approach to vertex
+  Float_t Yd; // Y coordinate of point of closest approach to vertex
+  Float_t Zd; // Z coordinate of point of closest approach to vertex
+
+  Float_t L; // track path length
+  Float_t D0; // track transverse impact parameter
+  Float_t DZ; // track longitudinal impact parameter
+  Float_t Nclusters; // Number of ionization clusters
+  Float_t dNdx; // Number of ionization clusters
+
+  Float_t ErrorP; // track momentum error
+  Float_t ErrorPT; // track transverse momentum error
+  Float_t ErrorPhi; // track azimuthal angle error
+  Float_t ErrorCtgTheta; // track cotangent of theta error
+
+  Float_t ErrorT; // time measurement error
+  Float_t ErrorD0; // track transverse impact parameter error
+  Float_t ErrorDZ; // track longitudinal impact parameter error
+  Float_t ErrorC; // track curvature error
+
+  // 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;
+
+  Int_t VertexIndex; // reference to vertex
+
+  static CompBase *fgCompare; //!
+  const CompBase *GetCompare() const { return fgCompare; }
+
+  TLorentzVector P4() const;
+  TMatrixDSym CovarianceMatrix() const;
+
+  Int_t NTimeHits; // number of hits contributing to time measurement
+
+  Float_t Eem; // calorimeter tower electromagnetic energy
+  Float_t Ehad; // calorimeter tower hadronic energy
+  Float_t Etrk; // total charged energy hitting tower
+
+  Float_t Edges[4]; // calorimeter tower edges
+
+  TRefArray Particles; // references to generated particles
+
+  ClassDef(ParticleFlowCandidate, 3)
+
 };
 
@@ -610,4 +715,5 @@
   Float_t Eem;
   Float_t Ehad;
+  Float_t Etrk;
 
   Float_t Edges[4];
@@ -618,9 +724,11 @@
 
   Float_t L; // path length
+  Float_t DZ;
+  Float_t ErrorDZ;
   Float_t ErrorT; // path length
   Float_t D0;
   Float_t ErrorD0;
-  Float_t DZ;
-  Float_t ErrorDZ;
+  Float_t C;
+  Float_t ErrorC;
   Float_t P;
   Float_t ErrorP;
@@ -632,4 +740,7 @@
   Float_t ErrorPhi;
 
+  Float_t Nclusters; // Number of ionization clusters
+  Float_t dNdx; // Number of ionization clusters per unit length
+
   Float_t Xd;
   Float_t Yd;
@@ -649,4 +760,7 @@
   Float_t PTD;
   Float_t FracPt[5];
+  Float_t NeutralEnergyFraction;  // charged energy fraction
+  Float_t ChargedEnergyFraction;  // neutral energy fraction
+
 
   // Timing information
@@ -663,4 +777,8 @@
   Float_t SumPtChargedPU;
   Float_t SumPt;
+
+  // ACTS compliant 6x6 track covariance (D0, phi, Curvature, dz, ctg(theta))
+
+  TMatrixDSym TrackCovariance;
 
   // vertex variables
@@ -698,4 +816,7 @@
   Double_t ExclYmerge56;
 
+  // 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: classes/DelphesFormula.cc
===================================================================
--- classes/DelphesFormula.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesFormula.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -18,4 +18,5 @@
 
 #include "classes/DelphesFormula.h"
+#include "classes/DelphesClasses.h"
 
 #include "TString.h"
@@ -63,4 +64,6 @@
   buffer.ReplaceAll("dz", "[1]");
   buffer.ReplaceAll("ctgTheta", "[2]");
+  buffer.ReplaceAll("radius", "[3]");
+  buffer.ReplaceAll("density", "[4]");
 
 #if ROOT_VERSION_CODE < ROOT_VERSION(6, 3, 0)
@@ -77,10 +80,18 @@
 //------------------------------------------------------------------------------
 
-Double_t DelphesFormula::Eval(Double_t pt, Double_t eta, Double_t phi,
-  Double_t energy, Double_t d0, Double_t dz,
-  Double_t ctgTheta)
+Double_t DelphesFormula::Eval(Double_t pt, Double_t eta, Double_t phi, Double_t energy, Candidate *candidate)
 {
+
+  Double_t d0 = 0., dz = 0., ctgTheta = 0., radius = 0., density = 0.;
+  if (candidate) {
+    d0 = candidate->D0;
+    dz = candidate->DZ;
+    ctgTheta = candidate->CtgTheta;
+    radius = candidate->Position.Pt();
+    density = candidate->ParticleDensity;
+  }
+    
   Double_t x[4] = {pt, eta, phi, energy};
-  Double_t params[3] = {d0, dz, ctgTheta};
+  Double_t params[5] = {d0, dz, ctgTheta, radius, density};
   return EvalPar(x, params);
 }
Index: classes/DelphesFormula.h
===================================================================
--- classes/DelphesFormula.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesFormula.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -22,4 +22,6 @@
 #include "TFormula.h"
 
+class Candidate;
+
 class DelphesFormula: public TFormula
 {
@@ -33,7 +35,5 @@
   Int_t Compile(const char *expression);
 
-  Double_t Eval(Double_t pt, Double_t eta = 0, Double_t phi = 0,
-    Double_t energy = 0, Double_t d0 = 0, Double_t dz = 0,
-    Double_t ctgTheta = 0);
+  Double_t Eval(Double_t pt, Double_t eta = 0, Double_t phi = 0, Double_t energy = 0, Candidate *candidate = nullptr);
 };
 
Index: classes/DelphesHepMC2Reader.cc
===================================================================
--- classes/DelphesHepMC2Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC2Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,489 @@
+/*
+ *  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 DelphesHepMC2Reader
+ *
+ *  Reads HepMC file
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesHepMC2Reader.h"
+
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+
+#include <map>
+#include <vector>
+
+#include <stdio.h>
+
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TObjArray.h"
+#include "TParticlePDG.h"
+#include "TStopwatch.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesStream.h"
+
+#include "ExRootAnalysis/ExRootTreeBranch.h"
+
+using namespace std;
+
+static const int kBufferSize = 16384;
+
+//---------------------------------------------------------------------------
+
+DelphesHepMC2Reader::DelphesHepMC2Reader() :
+  fInputFile(0), fBuffer(0), fPDG(0),
+  fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
+  fParticleCounter(0)
+{
+  fBuffer = new char[kBufferSize];
+
+  fPDG = TDatabasePDG::Instance();
+}
+
+//---------------------------------------------------------------------------
+
+DelphesHepMC2Reader::~DelphesHepMC2Reader()
+{
+  if(fBuffer) delete[] fBuffer;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::SetInputFile(FILE *inputFile)
+{
+  fInputFile = inputFile;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::Clear()
+{
+  fStateSize = 0;
+  fState.clear();
+  fWeightSize = 0;
+  fWeight.clear();
+  fMomentumCoefficient = 1.0;
+  fPositionCoefficient = 1.0;
+  fVertexCounter = -1;
+  fInCounter = -1;
+  fOutCounter = -1;
+  fMotherMap.clear();
+  fDaughterMap.clear();
+  fParticleCounter = 0;
+}
+
+//---------------------------------------------------------------------------
+
+bool DelphesHepMC2Reader::EventReady()
+{
+  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
+}
+
+//---------------------------------------------------------------------------
+
+bool DelphesHepMC2Reader::ReadBlock(DelphesFactory *factory,
+  TObjArray *allParticleOutputArray,
+  TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray)
+{
+  map<int, pair<int, int> >::iterator itMotherMap;
+  map<int, pair<int, int> >::iterator itDaughterMap;
+  char key, momentumUnit[4], positionUnit[3];
+  int i, rc, state;
+  double weight;
+
+  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
+
+  DelphesStream bufferStream(fBuffer + 1);
+
+  key = fBuffer[0];
+
+  if(key == 'E')
+  {
+    Clear();
+
+    rc = bufferStream.ReadInt(fEventNumber)
+      && bufferStream.ReadInt(fMPI)
+      && bufferStream.ReadDbl(fScale)
+      && bufferStream.ReadDbl(fAlphaQCD)
+      && bufferStream.ReadDbl(fAlphaQED)
+      && bufferStream.ReadInt(fProcessID)
+      && bufferStream.ReadInt(fSignalCode)
+      && bufferStream.ReadInt(fVertexCounter)
+      && bufferStream.ReadInt(fBeamCode[0])
+      && bufferStream.ReadInt(fBeamCode[1])
+      && bufferStream.ReadInt(fStateSize);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid event format" << endl;
+      return kFALSE;
+    }
+
+    for(i = 0; i < fStateSize; ++i)
+    {
+      rc = rc && bufferStream.ReadInt(state);
+      fState.push_back(state);
+    }
+
+    rc = rc && bufferStream.ReadInt(fWeightSize);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid event format" << endl;
+      return kFALSE;
+    }
+
+    for(i = 0; i < fWeightSize; ++i)
+    {
+      rc = rc && bufferStream.ReadDbl(weight);
+      fWeight.push_back(weight);
+    }
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid event format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'U')
+  {
+    rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
+
+    if(rc != 2)
+    {
+      cerr << "** ERROR: "
+           << "invalid units format" << endl;
+      return kFALSE;
+    }
+
+    if(strncmp(momentumUnit, "GEV", 3) == 0)
+    {
+      fMomentumCoefficient = 1.0;
+    }
+    else if(strncmp(momentumUnit, "MEV", 3) == 0)
+    {
+      fMomentumCoefficient = 0.001;
+    }
+
+    if(strncmp(positionUnit, "MM", 3) == 0)
+    {
+      fPositionCoefficient = 1.0;
+    }
+    else if(strncmp(positionUnit, "CM", 3) == 0)
+    {
+      fPositionCoefficient = 10.0;
+    }
+  }
+  else if(key == 'C')
+  {
+    rc = bufferStream.ReadDbl(fCrossSection)
+      && bufferStream.ReadDbl(fCrossSectionError);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid cross section format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'F')
+  {
+    rc = bufferStream.ReadInt(fID1)
+      && bufferStream.ReadInt(fID2)
+      && bufferStream.ReadDbl(fX1)
+      && bufferStream.ReadDbl(fX2)
+      && bufferStream.ReadDbl(fScalePDF)
+      && bufferStream.ReadDbl(fPDF1)
+      && bufferStream.ReadDbl(fPDF2);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid PDF format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'V' && fVertexCounter > 0)
+  {
+    rc = bufferStream.ReadInt(fOutVertexCode)
+      && bufferStream.ReadInt(fVertexID)
+      && bufferStream.ReadDbl(fX)
+      && bufferStream.ReadDbl(fY)
+      && bufferStream.ReadDbl(fZ)
+      && bufferStream.ReadDbl(fT)
+      && bufferStream.ReadInt(fInCounter)
+      && bufferStream.ReadInt(fOutCounter);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid vertex format" << endl;
+      return kFALSE;
+    }
+    --fVertexCounter;
+  }
+  else if(key == 'P' && fOutCounter > 0)
+  {
+    rc = bufferStream.ReadInt(fParticleCode)
+      && bufferStream.ReadInt(fPID)
+      && bufferStream.ReadDbl(fPx)
+      && bufferStream.ReadDbl(fPy)
+      && bufferStream.ReadDbl(fPz)
+      && bufferStream.ReadDbl(fE)
+      && bufferStream.ReadDbl(fMass)
+      && bufferStream.ReadInt(fStatus)
+      && bufferStream.ReadDbl(fTheta)
+      && bufferStream.ReadDbl(fPhi)
+      && bufferStream.ReadInt(fInVertexCode);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid particle format" << endl;
+      return kFALSE;
+    }
+
+    if(fInVertexCode < 0)
+    {
+      itMotherMap = fMotherMap.find(fInVertexCode);
+      if(itMotherMap == fMotherMap.end())
+      {
+        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
+      }
+      else
+      {
+        itMotherMap->second.second = fParticleCounter;
+      }
+    }
+
+    if(fInCounter <= 0)
+    {
+      itDaughterMap = fDaughterMap.find(fOutVertexCode);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
+      }
+      else
+      {
+        itDaughterMap->second.second = fParticleCounter;
+      }
+    }
+
+    AnalyzeParticle(factory, allParticleOutputArray,
+      stableParticleOutputArray, partonOutputArray);
+
+    if(fInCounter > 0)
+    {
+      --fInCounter;
+    }
+    else
+    {
+      --fOutCounter;
+    }
+
+    ++fParticleCounter;
+  }
+
+  if(EventReady())
+  {
+    FinalizeParticles(allParticleOutputArray);
+  }
+
+  return kTRUE;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
+  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
+{
+  HepMCEvent *element;
+
+  element = static_cast<HepMCEvent *>(branch->NewEntry());
+  element->Number = fEventNumber;
+
+  element->ProcessID = fProcessID;
+  element->MPI = fMPI;
+  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
+  element->CrossSection = fCrossSection;
+  element->CrossSectionError = fCrossSectionError;
+  element->Scale = fScale;
+  element->AlphaQED = fAlphaQED;
+  element->AlphaQCD = fAlphaQCD;
+
+  element->ID1 = fID1;
+  element->ID2 = fID2;
+  element->X1 = fX1;
+  element->X2 = fX2;
+  element->ScalePDF = fScalePDF;
+  element->PDF1 = fPDF1;
+  element->PDF2 = fPDF2;
+
+  element->ReadTime = readStopWatch->RealTime();
+  element->ProcTime = procStopWatch->RealTime();
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::AnalyzeWeight(ExRootTreeBranch *branch)
+{
+  Weight *element;
+  vector<double>::const_iterator itWeight;
+
+  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
+  {
+    element = static_cast<Weight *>(branch->NewEntry());
+
+    element->Weight = *itWeight;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::AnalyzeParticle(DelphesFactory *factory,
+  TObjArray *allParticleOutputArray,
+  TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray)
+{
+  Candidate *candidate;
+  TParticlePDG *pdgParticle;
+  int pdgCode;
+
+  candidate = factory->NewCandidate();
+
+  candidate->PID = fPID;
+  pdgCode = TMath::Abs(candidate->PID);
+
+  candidate->Status = fStatus;
+
+  pdgParticle = fPDG->GetParticle(fPID);
+  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
+  candidate->Mass = fMass;
+
+  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
+  if(fMomentumCoefficient != 1.0)
+  {
+    candidate->Momentum *= fMomentumCoefficient;
+  }
+
+  candidate->M2 = 1;
+  candidate->D2 = 1;
+  if(fInCounter > 0)
+  {
+    candidate->M1 = 1;
+    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  }
+  else
+  {
+    candidate->M1 = fOutVertexCode;
+    candidate->Position.SetXYZT(fX, fY, fZ, fT);
+    if(fPositionCoefficient != 1.0)
+    {
+      candidate->Position *= fPositionCoefficient;
+    }
+  }
+  if(fInVertexCode < 0)
+  {
+    candidate->D1 = fInVertexCode;
+  }
+  else
+  {
+    candidate->D1 = 1;
+  }
+
+  allParticleOutputArray->Add(candidate);
+
+  if(!pdgParticle) return;
+
+  if(fStatus == 1)
+  {
+    stableParticleOutputArray->Add(candidate);
+  }
+  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+  {
+    partonOutputArray->Add(candidate);
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC2Reader::FinalizeParticles(TObjArray *allParticleOutputArray)
+{
+  Candidate *candidate;
+  map<int, pair<int, int> >::iterator itMotherMap;
+  map<int, pair<int, int> >::iterator itDaughterMap;
+  int i;
+
+  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
+  {
+    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
+
+    if(candidate->M1 > 0)
+    {
+      candidate->M1 = -1;
+      candidate->M2 = -1;
+    }
+    else
+    {
+      itMotherMap = fMotherMap.find(candidate->M1);
+      if(itMotherMap == fMotherMap.end())
+      {
+        candidate->M1 = -1;
+        candidate->M2 = -1;
+      }
+      else
+      {
+        candidate->M1 = itMotherMap->second.first;
+        candidate->M2 = itMotherMap->second.second;
+      }
+    }
+    if(candidate->D1 > 0)
+    {
+      candidate->D1 = -1;
+      candidate->D2 = -1;
+    }
+    else
+    {
+      itDaughterMap = fDaughterMap.find(candidate->D1);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        candidate->D1 = -1;
+        candidate->D2 = -1;
+      }
+      else
+      {
+        candidate->D1 = itDaughterMap->second.first;
+        candidate->D2 = itDaughterMap->second.second;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
Index: classes/DelphesHepMC2Reader.h
===================================================================
--- classes/DelphesHepMC2Reader.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC2Reader.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,104 @@
+/*
+ *  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 DelphesHepMC2Reader_h
+#define DelphesHepMC2Reader_h
+
+/** \class DelphesHepMC2Reader
+ *
+ *  Reads HepMC file
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include <map>
+#include <vector>
+
+#include <stdio.h>
+
+class TObjArray;
+class TStopwatch;
+class TDatabasePDG;
+class ExRootTreeBranch;
+class DelphesFactory;
+
+class DelphesHepMC2Reader
+{
+public:
+  DelphesHepMC2Reader();
+  ~DelphesHepMC2Reader();
+
+  void SetInputFile(FILE *inputFile);
+
+  void Clear();
+  bool EventReady();
+
+  bool ReadBlock(DelphesFactory *factory,
+    TObjArray *allParticleOutputArray,
+    TObjArray *stableParticleOutputArray,
+    TObjArray *partonOutputArray);
+
+  void AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
+    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
+
+  void AnalyzeWeight(ExRootTreeBranch *branch);
+
+private:
+  void AnalyzeParticle(DelphesFactory *factory,
+    TObjArray *allParticleOutputArray,
+    TObjArray *stableParticleOutputArray,
+    TObjArray *partonOutputArray);
+
+  void FinalizeParticles(TObjArray *allParticleOutputArray);
+
+  FILE *fInputFile;
+
+  char *fBuffer;
+
+  TDatabasePDG *fPDG;
+
+  int fEventNumber, fMPI, fProcessID, fSignalCode, fVertexCounter, fBeamCode[2];
+  double fScale, fAlphaQCD, fAlphaQED;
+
+  double fMomentumCoefficient, fPositionCoefficient;
+
+  int fStateSize;
+  std::vector<int> fState;
+
+  int fWeightSize;
+  std::vector<double> fWeight;
+
+  double fCrossSection, fCrossSectionError;
+
+  int fID1, fID2;
+  double fX1, fX2, fScalePDF, fPDF1, fPDF2;
+
+  int fOutVertexCode, fVertexID, fInCounter, fOutCounter;
+  double fX, fY, fZ, fT;
+
+  int fParticleCode, fPID, fStatus, fInVertexCode;
+  double fPx, fPy, fPz, fE, fMass, fTheta, fPhi;
+
+  int fParticleCounter;
+
+  std::map<int, std::pair<int, int> > fMotherMap;
+  std::map<int, std::pair<int, int> > fDaughterMap;
+};
+
+#endif // DelphesHepMC2Reader_h
Index: classes/DelphesHepMC3Reader.cc
===================================================================
--- classes/DelphesHepMC3Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC3Reader.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,585 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2021  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 DelphesHepMC3Reader
+ *
+ *  Reads HepMC file
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesHepMC3Reader.h"
+
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+
+#include <map>
+#include <vector>
+
+#include <stdio.h>
+
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "TObjArray.h"
+#include "TParticlePDG.h"
+#include "TStopwatch.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesStream.h"
+
+#include "ExRootAnalysis/ExRootTreeBranch.h"
+
+using namespace std;
+
+static const int kBufferSize = 16384;
+
+//---------------------------------------------------------------------------
+
+DelphesHepMC3Reader::DelphesHepMC3Reader() :
+  fInputFile(0), fBuffer(0), fPDG(0),
+  fVertexCounter(-2), fParticleCounter(-1)
+{
+  fBuffer = new char[kBufferSize];
+
+  fPDG = TDatabasePDG::Instance();
+}
+
+//---------------------------------------------------------------------------
+
+DelphesHepMC3Reader::~DelphesHepMC3Reader()
+{
+  if(fBuffer) delete[] fBuffer;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::SetInputFile(FILE *inputFile)
+{
+  fInputFile = inputFile;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::Clear()
+{
+  fWeights.clear();
+  fMomentumCoefficient = 1.0;
+  fPositionCoefficient = 1.0;
+  fVertexCounter = -2;
+  fParticleCounter = -1;
+  fVertices.clear();
+  fParticles.clear();
+  fInVertexMap.clear();
+  fOutVertexMap.clear();
+  fMotherMap.clear();
+  fDaughterMap.clear();
+}
+
+//---------------------------------------------------------------------------
+
+bool DelphesHepMC3Reader::EventReady()
+{
+  return (fVertexCounter == -1) && (fParticleCounter == 0);
+}
+
+//---------------------------------------------------------------------------
+
+bool DelphesHepMC3Reader::ReadBlock(DelphesFactory *factory,
+  TObjArray *allParticleOutputArray,
+  TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray)
+{
+  map<int, pair<int, int> >::iterator itDaughterMap;
+  char key, momentumUnit[4], positionUnit[3];
+  int rc, code;
+  double weight;
+
+  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
+
+  DelphesStream bufferStream(fBuffer + 1);
+
+  key = fBuffer[0];
+
+  if(key == 'E')
+  {
+    Clear();
+
+    rc = bufferStream.ReadInt(fEventNumber)
+      && bufferStream.ReadInt(fVertexCounter)
+      && bufferStream.ReadInt(fParticleCounter);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid event format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'U')
+  {
+    rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
+
+    if(rc != 2)
+    {
+      cerr << "** ERROR: "
+           << "invalid units format" << endl;
+      return kFALSE;
+    }
+
+    if(strncmp(momentumUnit, "GEV", 3) == 0)
+    {
+      fMomentumCoefficient = 1.0;
+    }
+    else if(strncmp(momentumUnit, "MEV", 3) == 0)
+    {
+      fMomentumCoefficient = 0.001;
+    }
+
+    if(strncmp(positionUnit, "MM", 3) == 0)
+    {
+      fPositionCoefficient = 1.0;
+    }
+    else if(strncmp(positionUnit, "CM", 3) == 0)
+    {
+      fPositionCoefficient = 10.0;
+    }
+  }
+  else if(key == 'W')
+  {
+    while(bufferStream.ReadDbl(weight))
+    {
+      fWeights.push_back(weight);
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("mpi"))
+  {
+    rc = bufferStream.ReadInt(fMPI);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid MPI format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("signal_process_id"))
+  {
+    rc = bufferStream.ReadInt(fProcessID);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid process ID format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("event_scale"))
+  {
+    rc = bufferStream.ReadDbl(fScale);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid event scale format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("alphaQCD"))
+  {
+    rc = bufferStream.ReadDbl(fAlphaQCD);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid alphaQCD format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("alphaQED"))
+  {
+    rc = bufferStream.ReadDbl(fAlphaQED);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid alphaQED format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("GenCrossSection"))
+  {
+    rc = bufferStream.ReadDbl(fCrossSection)
+      && bufferStream.ReadDbl(fCrossSectionError);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid cross section format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'A' && bufferStream.FindStr("GenPdfInfo"))
+  {
+    rc = bufferStream.ReadInt(fID1)
+      && bufferStream.ReadInt(fID2)
+      && bufferStream.ReadDbl(fX1)
+      && bufferStream.ReadDbl(fX2)
+      && bufferStream.ReadDbl(fScalePDF)
+      && bufferStream.ReadDbl(fPDF1)
+      && bufferStream.ReadDbl(fPDF2);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid PDF format" << endl;
+      return kFALSE;
+    }
+  }
+  else if(key == 'V')
+  {
+    fParticles.clear();
+
+    fX = 0.0;
+    fY = 0.0;
+    fZ = 0.0;
+    fT = 0.0;
+
+    rc = bufferStream.ReadInt(fVertexCode)
+      && bufferStream.ReadInt(fVertexStatus);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid vertex format" << endl;
+      return kFALSE;
+    }
+
+    rc = bufferStream.FindChr('[');
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid vertex format" << endl;
+      return kFALSE;
+    }
+
+    while(bufferStream.ReadInt(code))
+    {
+      fParticles.push_back(code);
+      bufferStream.FindChr(',');
+    }
+
+    if(bufferStream.FindChr('@'))
+    {
+      rc = bufferStream.ReadDbl(fX)
+        && bufferStream.ReadDbl(fY)
+        && bufferStream.ReadDbl(fZ)
+        && bufferStream.ReadDbl(fT);
+
+      if(!rc)
+      {
+        cerr << "** ERROR: "
+             << "invalid vertex format" << endl;
+        return kFALSE;
+      }
+    }
+
+    AnalyzeVertex(factory, fVertexCode);
+  }
+  else if(key == 'P' && fParticleCounter > 0)
+  {
+    --fParticleCounter;
+
+    rc = bufferStream.ReadInt(fParticleCode)
+      && bufferStream.ReadInt(fOutVertexCode)
+      && bufferStream.ReadInt(fPID)
+      && bufferStream.ReadDbl(fPx)
+      && bufferStream.ReadDbl(fPy)
+      && bufferStream.ReadDbl(fPz)
+      && bufferStream.ReadDbl(fE)
+      && bufferStream.ReadDbl(fMass)
+      && bufferStream.ReadInt(fParticleStatus);
+
+    if(!rc)
+    {
+      cerr << "** ERROR: "
+           << "invalid particle format" << endl;
+      return kFALSE;
+    }
+
+    AnalyzeParticle(factory);
+  }
+
+  if(EventReady())
+  {
+    FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
+  }
+
+  return kTRUE;
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
+  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
+{
+  HepMCEvent *element;
+
+  element = static_cast<HepMCEvent *>(branch->NewEntry());
+  element->Number = fEventNumber;
+
+  element->ProcessID = fProcessID;
+  element->MPI = fMPI;
+  element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
+  element->CrossSection = fCrossSection;
+  element->CrossSectionError = fCrossSectionError;
+  element->Scale = fScale;
+  element->AlphaQED = fAlphaQED;
+  element->AlphaQCD = fAlphaQCD;
+
+  element->ID1 = fID1;
+  element->ID2 = fID2;
+  element->X1 = fX1;
+  element->X2 = fX2;
+  element->ScalePDF = fScalePDF;
+  element->PDF1 = fPDF1;
+  element->PDF2 = fPDF2;
+
+  element->ReadTime = readStopWatch->RealTime();
+  element->ProcTime = procStopWatch->RealTime();
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeWeight(ExRootTreeBranch *branch)
+{
+  Weight *element;
+  vector<double>::const_iterator itWeights;
+
+  for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
+  {
+    element = static_cast<Weight *>(branch->NewEntry());
+
+    element->Weight = *itWeights;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate)
+{
+  int index;
+  TLorentzVector *position;
+  TObjArray *array;
+  vector<int>::iterator itParticle;
+  map<int, int>::iterator itVertexMap;
+
+  itVertexMap = fOutVertexMap.find(code);
+  if(itVertexMap == fOutVertexMap.end())
+  {
+    --fVertexCounter;
+
+    index = fVertices.size();
+    fOutVertexMap[code] = index;
+    if(candidate && code > 0) fInVertexMap[code] = index;
+
+    position = factory->New<TLorentzVector>();
+    array = factory->NewArray();
+    position->SetXYZT(0.0, 0.0, 0.0, 0.0);
+    fVertices.push_back(make_pair(position, array));
+  }
+  else
+  {
+    index = itVertexMap->second;
+    position = fVertices[index].first;
+    array = fVertices[index].second;
+  }
+
+  if(candidate)
+  {
+    array->Add(candidate);
+  }
+  else
+  {
+    position->SetXYZT(fX, fY, fZ, fT);
+    for(itParticle = fParticles.begin(); itParticle != fParticles.end(); ++itParticle)
+    {
+      fInVertexMap[*itParticle] = index;
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory)
+{
+  Candidate *candidate;
+
+  candidate = factory->NewCandidate();
+
+  candidate->PID = fPID;
+
+  candidate->Status = fParticleStatus;
+
+  candidate->Mass = fMass;
+
+  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
+
+  candidate->D1 = fParticleCode;
+
+  AnalyzeVertex(factory, fOutVertexCode, candidate);
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray,
+  TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray)
+{
+  TLorentzVector *position;
+  TObjArray *array;
+  Candidate *candidate;
+  TParticlePDG *pdgParticle;
+  int pdgCode;
+  map<int, int >::iterator itVertexMap;
+  map<int, pair<int, int> >::iterator itMotherMap;
+  map<int, pair<int, int> >::iterator itDaughterMap;
+  int i, j, code, counter;
+
+  counter = 0;
+  for(i = 0; i < fVertices.size(); ++i)
+  {
+    position = fVertices[i].first;
+    array = fVertices[i].second;
+
+    for(j = 0; j < array->GetEntriesFast(); ++j)
+    {
+      candidate = static_cast<Candidate *>(array->At(j));
+
+      candidate->Position = *position;
+      if(fPositionCoefficient != 1.0)
+      {
+        candidate->Position *= fPositionCoefficient;
+      }
+
+      if(fMomentumCoefficient != 1.0)
+      {
+        candidate->Momentum *= fMomentumCoefficient;
+      }
+
+      candidate->M1 = i;
+
+      itDaughterMap = fDaughterMap.find(i);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        fDaughterMap[i] = make_pair(counter, counter);
+      }
+      else
+      {
+        itDaughterMap->second.second = counter;
+      }
+
+      code = candidate->D1;
+
+      itVertexMap = fInVertexMap.find(code);
+      if(itVertexMap == fInVertexMap.end())
+      {
+        candidate->D1 = -1;
+      }
+      else
+      {
+        code = itVertexMap->second;
+
+        candidate->D1 = code;
+
+        itMotherMap = fMotherMap.find(code);
+        if(itMotherMap == fMotherMap.end())
+        {
+          fMotherMap[code] = make_pair(counter, -1);
+        }
+        else
+        {
+          itMotherMap->second.second = counter;
+        }
+      }
+
+      allParticleOutputArray->Add(candidate);
+
+      ++counter;
+
+      pdgParticle = fPDG->GetParticle(candidate->PID);
+
+      candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
+
+      if(!pdgParticle) continue;
+
+      pdgCode = TMath::Abs(candidate->PID);
+
+      if(candidate->Status == 1)
+      {
+        stableParticleOutputArray->Add(candidate);
+      }
+      else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+      {
+        partonOutputArray->Add(candidate);
+      }
+    }
+  }
+
+  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
+  {
+    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
+
+    itMotherMap = fMotherMap.find(candidate->M1);
+    if(itMotherMap == fMotherMap.end())
+    {
+      candidate->M1 = -1;
+      candidate->M2 = -1;
+    }
+    else
+    {
+      candidate->M1 = itMotherMap->second.first;
+      candidate->M2 = itMotherMap->second.second;
+    }
+
+    if(candidate->D1 < 0)
+    {
+      candidate->D1 = -1;
+      candidate->D2 = -1;
+    }
+    else
+    {
+      itDaughterMap = fDaughterMap.find(candidate->D1);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        candidate->D1 = -1;
+        candidate->D2 = -1;
+      }
+      else
+      {
+        candidate->D1 = itDaughterMap->second.first;
+        candidate->D2 = itDaughterMap->second.second;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
Index: classes/DelphesHepMC3Reader.h
===================================================================
--- classes/DelphesHepMC3Reader.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
+++ classes/DelphesHepMC3Reader.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -0,0 +1,107 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2021  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 DelphesHepMC3Reader_h
+#define DelphesHepMC3Reader_h
+
+/** \class DelphesHepMC3Reader
+ *
+ *  Reads HepMC file
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include <map>
+#include <vector>
+
+#include <stdio.h>
+
+class TObjArray;
+class TStopwatch;
+class TDatabasePDG;
+class TLorentzVector;
+class ExRootTreeBranch;
+class DelphesFactory;
+class Candidate;
+
+class DelphesHepMC3Reader
+{
+public:
+  DelphesHepMC3Reader();
+  ~DelphesHepMC3Reader();
+
+  void SetInputFile(FILE *inputFile);
+
+  void Clear();
+  bool EventReady();
+
+  bool ReadBlock(DelphesFactory *factory,
+    TObjArray *allParticleOutputArray,
+    TObjArray *stableParticleOutputArray,
+    TObjArray *partonOutputArray);
+
+  void AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
+    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
+
+  void AnalyzeWeight(ExRootTreeBranch *branch);
+
+private:
+  void AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate = 0);
+
+  void AnalyzeParticle(DelphesFactory *factory);
+
+  void FinalizeParticles(TObjArray *allParticleOutputArray,
+    TObjArray *stableParticleOutputArray,
+    TObjArray *partonOutputArray);
+
+  FILE *fInputFile;
+
+  char *fBuffer;
+
+  TDatabasePDG *fPDG;
+
+  int fEventNumber, fMPI, fProcessID, fSignalCode, fVertexCounter, fParticleCounter;
+  double fScale, fAlphaQCD, fAlphaQED;
+
+  double fMomentumCoefficient, fPositionCoefficient;
+
+  std::vector<double> fWeights;
+
+  double fCrossSection, fCrossSectionError;
+
+  int fID1, fID2;
+  double fX1, fX2, fScalePDF, fPDF1, fPDF2;
+
+  int fVertexCode, fVertexStatus;
+  double fX, fY, fZ, fT;
+
+  int fParticleCode, fPID, fParticleStatus, fOutVertexCode;
+  double fPx, fPy, fPz, fE, fMass;
+
+  std::vector<std::pair<TLorentzVector *, TObjArray *> > fVertices;
+  std::vector<int> fParticles;
+
+  std::map<int, int> fInVertexMap;
+  std::map<int, int> fOutVertexMap;
+
+  std::map<int, std::pair<int, int> > fMotherMap;
+  std::map<int, std::pair<int, int> > fDaughterMap;
+};
+
+#endif // DelphesHepMC3Reader_h
Index: classes/DelphesHepMCReader.cc
===================================================================
--- classes/DelphesHepMCReader.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ 	(revision )
@@ -1,484 +1,0 @@
-/*
- *  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 DelphesHepMCReader
- *
- *  Reads HepMC file
- *
- *  \author P. Demin - UCL, Louvain-la-Neuve
- *
- */
-
-#include "classes/DelphesHepMCReader.h"
-
-#include <iostream>
-#include <sstream>
-#include <stdexcept>
-
-#include <map>
-#include <vector>
-
-#include <stdio.h>
-
-#include "TDatabasePDG.h"
-#include "TLorentzVector.h"
-#include "TObjArray.h"
-#include "TParticlePDG.h"
-#include "TStopwatch.h"
-
-#include "classes/DelphesClasses.h"
-#include "classes/DelphesFactory.h"
-#include "classes/DelphesStream.h"
-
-#include "ExRootAnalysis/ExRootTreeBranch.h"
-
-using namespace std;
-
-static const int kBufferSize = 16384;
-
-//---------------------------------------------------------------------------
-
-DelphesHepMCReader::DelphesHepMCReader() :
-  fInputFile(0), fBuffer(0), fPDG(0),
-  fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
-  fParticleCounter(0)
-{
-  fBuffer = new char[kBufferSize];
-
-  fPDG = TDatabasePDG::Instance();
-}
-
-//---------------------------------------------------------------------------
-
-DelphesHepMCReader::~DelphesHepMCReader()
-{
-  if(fBuffer) delete[] fBuffer;
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::SetInputFile(FILE *inputFile)
-{
-  fInputFile = inputFile;
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::Clear()
-{
-  fStateSize = 0;
-  fState.clear();
-  fWeightSize = 0;
-  fWeight.clear();
-  fMomentumCoefficient = 1.0;
-  fPositionCoefficient = 1.0;
-  fVertexCounter = -1;
-  fInCounter = -1;
-  fOutCounter = -1;
-  fMotherMap.clear();
-  fDaughterMap.clear();
-  fParticleCounter = 0;
-}
-
-//---------------------------------------------------------------------------
-
-bool DelphesHepMCReader::EventReady()
-{
-  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
-}
-
-//---------------------------------------------------------------------------
-
-bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
-  TObjArray *allParticleOutputArray,
-  TObjArray *stableParticleOutputArray,
-  TObjArray *partonOutputArray)
-{
-  map<int, pair<int, int> >::iterator itMotherMap;
-  map<int, pair<int, int> >::iterator itDaughterMap;
-  char key, momentumUnit[4], positionUnit[3];
-  int i, rc, state;
-  double weight;
-
-  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
-
-  DelphesStream bufferStream(fBuffer + 1);
-
-  key = fBuffer[0];
-
-  if(key == 'E')
-  {
-    Clear();
-
-    rc = bufferStream.ReadInt(fEventNumber)
-      && bufferStream.ReadInt(fMPI)
-      && bufferStream.ReadDbl(fScale)
-      && bufferStream.ReadDbl(fAlphaQCD)
-      && bufferStream.ReadDbl(fAlphaQED)
-      && bufferStream.ReadInt(fProcessID)
-      && bufferStream.ReadInt(fSignalCode)
-      && bufferStream.ReadInt(fVertexCounter)
-      && bufferStream.ReadInt(fBeamCode[0])
-      && bufferStream.ReadInt(fBeamCode[1])
-      && bufferStream.ReadInt(fStateSize);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid event format" << endl;
-      return kFALSE;
-    }
-
-    for(i = 0; i < fStateSize; ++i)
-    {
-      rc = rc && bufferStream.ReadInt(state);
-      fState.push_back(state);
-    }
-
-    rc = rc && bufferStream.ReadInt(fWeightSize);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid event format" << endl;
-      return kFALSE;
-    }
-
-    for(i = 0; i < fWeightSize; ++i)
-    {
-      rc = rc && bufferStream.ReadDbl(weight);
-      fWeight.push_back(weight);
-    }
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid event format" << endl;
-      return kFALSE;
-    }
-  }
-  else if(key == 'U')
-  {
-    rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
-
-    if(rc != 2)
-    {
-      cerr << "** ERROR: "
-           << "invalid units format" << endl;
-      return kFALSE;
-    }
-
-    if(strncmp(momentumUnit, "GEV", 3) == 0)
-    {
-      fMomentumCoefficient = 1.0;
-    }
-    else if(strncmp(momentumUnit, "MEV", 3) == 0)
-    {
-      fMomentumCoefficient = 0.001;
-    }
-
-    if(strncmp(positionUnit, "MM", 3) == 0)
-    {
-      fPositionCoefficient = 1.0;
-    }
-    else if(strncmp(positionUnit, "CM", 3) == 0)
-    {
-      fPositionCoefficient = 10.0;
-    }
-  }
-
-  else if(key == 'C')
-  {
-    rc = bufferStream.ReadDbl(fCrossSection)
-      && bufferStream.ReadDbl(fCrossSectionError);
-  }
-
-  else if(key == 'F')
-  {
-    rc = bufferStream.ReadInt(fID1)
-      && bufferStream.ReadInt(fID2)
-      && bufferStream.ReadDbl(fX1)
-      && bufferStream.ReadDbl(fX2)
-      && bufferStream.ReadDbl(fScalePDF)
-      && bufferStream.ReadDbl(fPDF1)
-      && bufferStream.ReadDbl(fPDF2);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid PDF format" << endl;
-      return kFALSE;
-    }
-  }
-  else if(key == 'V' && fVertexCounter > 0)
-  {
-    rc = bufferStream.ReadInt(fOutVertexCode)
-      && bufferStream.ReadInt(fVertexID)
-      && bufferStream.ReadDbl(fX)
-      && bufferStream.ReadDbl(fY)
-      && bufferStream.ReadDbl(fZ)
-      && bufferStream.ReadDbl(fT)
-      && bufferStream.ReadInt(fInCounter)
-      && bufferStream.ReadInt(fOutCounter);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid vertex format" << endl;
-      return kFALSE;
-    }
-    --fVertexCounter;
-  }
-  else if(key == 'P' && fOutCounter > 0)
-  {
-    rc = bufferStream.ReadInt(fParticleCode)
-      && bufferStream.ReadInt(fPID)
-      && bufferStream.ReadDbl(fPx)
-      && bufferStream.ReadDbl(fPy)
-      && bufferStream.ReadDbl(fPz)
-      && bufferStream.ReadDbl(fE)
-      && bufferStream.ReadDbl(fMass)
-      && bufferStream.ReadInt(fStatus)
-      && bufferStream.ReadDbl(fTheta)
-      && bufferStream.ReadDbl(fPhi)
-      && bufferStream.ReadInt(fInVertexCode);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: "
-           << "invalid particle format" << endl;
-      return kFALSE;
-    }
-
-    if(fInVertexCode < 0)
-    {
-      itMotherMap = fMotherMap.find(fInVertexCode);
-      if(itMotherMap == fMotherMap.end())
-      {
-        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
-      }
-      else
-      {
-        itMotherMap->second.second = fParticleCounter;
-      }
-    }
-
-    if(fInCounter <= 0)
-    {
-      itDaughterMap = fDaughterMap.find(fOutVertexCode);
-      if(itDaughterMap == fDaughterMap.end())
-      {
-        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
-      }
-      else
-      {
-        itDaughterMap->second.second = fParticleCounter;
-      }
-    }
-
-    AnalyzeParticle(factory, allParticleOutputArray,
-      stableParticleOutputArray, partonOutputArray);
-
-    if(fInCounter > 0)
-    {
-      --fInCounter;
-    }
-    else
-    {
-      --fOutCounter;
-    }
-
-    ++fParticleCounter;
-  }
-
-  if(EventReady())
-  {
-    FinalizeParticles(allParticleOutputArray);
-  }
-
-  return kTRUE;
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
-  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
-{
-  HepMCEvent *element;
-
-  element = static_cast<HepMCEvent *>(branch->NewEntry());
-  element->Number = fEventNumber;
-
-  element->ProcessID = fProcessID;
-  element->MPI = fMPI;
-  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
-  element->CrossSection = fCrossSection;
-  element->CrossSectionError = fCrossSectionError;
-  element->Scale = fScale;
-  element->AlphaQED = fAlphaQED;
-  element->AlphaQCD = fAlphaQCD;
-
-  element->ID1 = fID1;
-  element->ID2 = fID2;
-  element->X1 = fX1;
-  element->X2 = fX2;
-  element->ScalePDF = fScalePDF;
-  element->PDF1 = fPDF1;
-  element->PDF2 = fPDF2;
-
-  element->ReadTime = readStopWatch->RealTime();
-  element->ProcTime = procStopWatch->RealTime();
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::AnalyzeWeight(ExRootTreeBranch *branch)
-{
-  Weight *element;
-  vector<double>::const_iterator itWeight;
-
-  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
-  {
-    element = static_cast<Weight *>(branch->NewEntry());
-
-    element->Weight = *itWeight;
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
-  TObjArray *allParticleOutputArray,
-  TObjArray *stableParticleOutputArray,
-  TObjArray *partonOutputArray)
-{
-  Candidate *candidate;
-  TParticlePDG *pdgParticle;
-  int pdgCode;
-
-  candidate = factory->NewCandidate();
-
-  candidate->PID = fPID;
-  pdgCode = TMath::Abs(candidate->PID);
-
-  candidate->Status = fStatus;
-
-  pdgParticle = fPDG->GetParticle(fPID);
-  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
-  candidate->Mass = fMass;
-
-  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
-  if(fMomentumCoefficient != 1.0)
-  {
-    candidate->Momentum *= fMomentumCoefficient;
-  }
-
-  candidate->M2 = 1;
-  candidate->D2 = 1;
-  if(fInCounter > 0)
-  {
-    candidate->M1 = 1;
-    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
-  }
-  else
-  {
-    candidate->M1 = fOutVertexCode;
-    candidate->Position.SetXYZT(fX, fY, fZ, fT);
-    if(fPositionCoefficient != 1.0)
-    {
-      candidate->Position *= fPositionCoefficient;
-    }
-  }
-  if(fInVertexCode < 0)
-  {
-    candidate->D1 = fInVertexCode;
-  }
-  else
-  {
-    candidate->D1 = 1;
-  }
-
-  allParticleOutputArray->Add(candidate);
-
-  if(!pdgParticle) return;
-
-  if(fStatus == 1)
-  {
-    stableParticleOutputArray->Add(candidate);
-  }
-  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
-  {
-    partonOutputArray->Add(candidate);
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
-{
-  Candidate *candidate;
-  map<int, pair<int, int> >::iterator itMotherMap;
-  map<int, pair<int, int> >::iterator itDaughterMap;
-  int i;
-
-  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
-  {
-    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
-
-    if(candidate->M1 > 0)
-    {
-      candidate->M1 = -1;
-      candidate->M2 = -1;
-    }
-    else
-    {
-      itMotherMap = fMotherMap.find(candidate->M1);
-      if(itMotherMap == fMotherMap.end())
-      {
-        candidate->M1 = -1;
-        candidate->M2 = -1;
-      }
-      else
-      {
-        candidate->M1 = itMotherMap->second.first;
-        candidate->M2 = itMotherMap->second.second;
-      }
-    }
-    if(candidate->D1 > 0)
-    {
-      candidate->D1 = -1;
-      candidate->D2 = -1;
-    }
-    else
-    {
-      itDaughterMap = fDaughterMap.find(candidate->D1);
-      if(itDaughterMap == fDaughterMap.end())
-      {
-        candidate->D1 = -1;
-        candidate->D2 = -1;
-      }
-      else
-      {
-        candidate->D1 = itDaughterMap->second.first;
-        candidate->D2 = itDaughterMap->second.second;
-      }
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
Index: classes/DelphesHepMCReader.h
===================================================================
--- classes/DelphesHepMCReader.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ 	(revision )
@@ -1,104 +1,0 @@
-/*
- *  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 DelphesHepMCReader_h
-#define DelphesHepMCReader_h
-
-/** \class DelphesHepMCReader
- *
- *  Reads HepMC file
- *
- *  \author P. Demin - UCL, Louvain-la-Neuve
- *
- */
-
-#include <map>
-#include <vector>
-
-#include <stdio.h>
-
-class TObjArray;
-class TStopwatch;
-class TDatabasePDG;
-class ExRootTreeBranch;
-class DelphesFactory;
-
-class DelphesHepMCReader
-{
-public:
-  DelphesHepMCReader();
-  ~DelphesHepMCReader();
-
-  void SetInputFile(FILE *inputFile);
-
-  void Clear();
-  bool EventReady();
-
-  bool ReadBlock(DelphesFactory *factory,
-    TObjArray *allParticleOutputArray,
-    TObjArray *stableParticleOutputArray,
-    TObjArray *partonOutputArray);
-
-  void AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
-    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
-
-  void AnalyzeWeight(ExRootTreeBranch *branch);
-
-private:
-  void AnalyzeParticle(DelphesFactory *factory,
-    TObjArray *allParticleOutputArray,
-    TObjArray *stableParticleOutputArray,
-    TObjArray *partonOutputArray);
-
-  void FinalizeParticles(TObjArray *allParticleOutputArray);
-
-  FILE *fInputFile;
-
-  char *fBuffer;
-
-  TDatabasePDG *fPDG;
-
-  int fEventNumber, fMPI, fProcessID, fSignalCode, fVertexCounter, fBeamCode[2];
-  double fScale, fAlphaQCD, fAlphaQED;
-
-  double fMomentumCoefficient, fPositionCoefficient;
-
-  int fStateSize;
-  std::vector<int> fState;
-
-  int fWeightSize;
-  std::vector<double> fWeight;
-
-  double fCrossSection, fCrossSectionError;
-
-  int fID1, fID2;
-  double fX1, fX2, fScalePDF, fPDF1, fPDF2;
-
-  int fOutVertexCode, fVertexID, fInCounter, fOutCounter;
-  double fX, fY, fZ, fT;
-
-  int fParticleCode, fPID, fStatus, fInVertexCode;
-  double fPx, fPy, fPz, fE, fMass, fTheta, fPhi;
-
-  int fParticleCounter;
-
-  std::map<int, std::pair<int, int> > fMotherMap;
-  std::map<int, std::pair<int, int> > fDaughterMap;
-};
-
-#endif // DelphesHepMCReader_h
Index: classes/DelphesModule.cc
===================================================================
--- classes/DelphesModule.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesModule.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -130,4 +130,21 @@
 //------------------------------------------------------------------------------
 
+void DelphesModule::AddInfo(const char *name, Double_t value)
+{
+  stringstream message;
+  if(!fTreeWriter)
+  {
+    fTreeWriter = static_cast<ExRootTreeWriter *>(GetObject("TreeWriter", ExRootTreeWriter::Class()));
+    if(!fTreeWriter)
+    {
+      message << "can't access access tree writer";
+      throw runtime_error(message.str());
+    }
+  }
+  fTreeWriter->AddInfo(name, value);
+}
+
+//------------------------------------------------------------------------------
+
 ExRootResult *DelphesModule::GetPlots()
 {
Index: classes/DelphesModule.h
===================================================================
--- classes/DelphesModule.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesModule.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -55,4 +55,5 @@
 
   ExRootTreeBranch *NewBranch(const char *name, TClass *cl);
+  void AddInfo(const char *name, Double_t value);
 
   ExRootResult *GetPlots();
Index: classes/DelphesStream.cc
===================================================================
--- classes/DelphesStream.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesStream.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -106,2 +106,38 @@
 
 //------------------------------------------------------------------------------
+
+bool DelphesStream::FindChr(int value)
+{
+  char *position;
+  bool result = false;
+
+  position = strchr(fBuffer, value);
+
+  if(position)
+  {
+    result = true;
+    fBuffer = position + 1;
+  }
+
+  return result;
+}
+
+//------------------------------------------------------------------------------
+
+bool DelphesStream::FindStr(const char *value)
+{
+  char *position;
+  bool result = false;
+
+  position = strstr(fBuffer, value);
+
+  if(position)
+  {
+    result = true;
+    fBuffer = position + strlen(value);
+  }
+
+  return result;
+}
+
+//------------------------------------------------------------------------------
Index: classes/DelphesStream.h
===================================================================
--- classes/DelphesStream.h	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesStream.h	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -35,4 +35,6 @@
   bool ReadDbl(double &value);
   bool ReadInt(int &value);
+  bool FindChr(int value);
+  bool FindStr(const char *value);
 
 private:
Index: classes/DelphesTF2.cc
===================================================================
--- classes/DelphesTF2.cc	(revision a5af1df1cf2ca5a488fc610ab546c762b7cdb5f2)
+++ classes/DelphesTF2.cc	(revision d612deca94f89952da103f30d46e21c02039b8bc)
@@ -32,5 +32,7 @@
 {
 
-#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 3, 0)
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 22, 7)
+  fFormula = std::unique_ptr<TFormula>(new TFormula());
+#elif ROOT_VERSION_CODE >= ROOT_VERSION(6, 3, 0)
   fFormula = new TFormula();
 #endif
