Changes in / [ec5e04b:7993cad] in git

Files:

• Makefile (modified) (11 diffs)
• cards/converter_card.tcl (modified) (1 diff)
• cards/delphes_card_CMS.tcl (modified) (4 diffs)
• cards/delphes_card_CMS_simple.tcl (added)
• classes/DelphesClasses.cc (modified) (8 diffs)
• classes/DelphesClasses.h (modified) (10 diffs)
• classes/SortableObject.h (modified) (1 diff)
• doc/genMakefile.tcl (modified) (1 diff)
• examples/Pythia8/configParticleGun.cmnd (modified) (1 diff)
• examples/Validation.C (added)
• examples/Validation.cpp (added)
• modules/BeamSpotFilter.cc (added)
• modules/BeamSpotFilter.h (added)
• modules/Calorimeter.cc (modified) (8 diffs)
• modules/Calorimeter.h (modified) (3 diffs)
• modules/ImpactParameterSmearing.cc (modified) (4 diffs)
• modules/ModulesLinkDef.h (modified) (6 diffs)
• modules/ParticlePropagator.cc (modified) (8 diffs)
• modules/ParticlePropagator.h (modified) (2 diffs)
• modules/PileUpMerger.cc (modified) (9 diffs)
• modules/SimpleCalorimeter.cc (modified) (10 diffs)
• modules/SimpleCalorimeter.h (modified) (3 diffs)
• modules/TimeSmearing.cc (modified) (2 diffs)
• modules/TrackCountingBTagging.cc (modified) (2 diffs)
• modules/TrackSmearing.cc (added)
• modules/TrackSmearing.h (added)
• modules/TreeWriter.cc (modified) (7 diffs)
• modules/VertexFinder.cc (added)
• modules/VertexFinder.h (added)
• modules/VertexFinderDA4D.cc (added)
• modules/VertexFinderDA4D.h (added)
• modules/VertexSorter.cc (added)
• modules/VertexSorter.h (added)
• readers/DelphesCMSFWLite.cpp (modified) (1 diff)
• readers/DelphesROOT.cpp (added)
• validation.sh (added)

Makefile

@@ -319 +319 @@
         modules/EnergySmearing.h \
         modules/MomentumSmearing.h \
+        modules/TrackSmearing.h \
         modules/ImpactParameterSmearing.h \
         modules/TimeSmearing.h \
@@ -342 +343 @@
         modules/StatusPidFilter.h \
         modules/PdgCodeFilter.h \
+        modules/BeamSpotFilter.h \
         modules/Cloner.h \
         modules/Weighter.h \
@@ -347 +349 @@
         modules/JetFlavorAssociation.h \
         modules/JetFakeParticle.h \
+        modules/VertexSorter.h \
+        modules/VertexFinder.h \
         modules/ExampleModule.h
 tmp/modules/ModulesDict$(PcmSuf): \
@@ -553 +557 @@
         classes/DelphesFactory.h \
         classes/DelphesFormula.h
+tmp/modules/BeamSpotFilter.$(ObjSuf): \
+        modules/BeamSpotFilter.$(SrcSuf) \
+        modules/BeamSpotFilter.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/Calorimeter.$(ObjSuf): \
         modules/Calorimeter.$(SrcSuf) \
@@ -852 +865 @@
         external/ExRootAnalysis/ExRootFilter.h \
         external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/TrackSmearing.$(ObjSuf): \
+        modules/TrackSmearing.$(SrcSuf) \
+        modules/TrackSmearing.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/TreeWriter.$(ObjSuf): \
         modules/TreeWriter.$(SrcSuf) \
@@ -868 +890 @@
         classes/DelphesFactory.h \
         classes/DelphesFormula.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/VertexFinder.$(ObjSuf): \
+        modules/VertexFinder.$(SrcSuf) \
+        modules/VertexFinder.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        classes/DelphesPileUpReader.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/VertexSorter.$(ObjSuf): \
+        modules/VertexSorter.$(SrcSuf) \
+        modules/VertexSorter.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        classes/DelphesPileUpReader.h \
         external/ExRootAnalysis/ExRootResult.h \
         external/ExRootAnalysis/ExRootFilter.h \
@@ -931 +973 @@
         tmp/modules/AngularSmearing.$(ObjSuf) \
         tmp/modules/BTagging.$(ObjSuf) \
+        tmp/modules/BeamSpotFilter.$(ObjSuf) \
         tmp/modules/Calorimeter.$(ObjSuf) \
         tmp/modules/Cloner.$(ObjSuf) \
@@ -963 +1006 @@
         tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
         tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
+        tmp/modules/TrackSmearing.$(ObjSuf) \
         tmp/modules/TreeWriter.$(ObjSuf) \
         tmp/modules/UniqueObjectFinder.$(ObjSuf) \
+        tmp/modules/VertexFinder.$(ObjSuf) \
+        tmp/modules/VertexSorter.$(ObjSuf) \
         tmp/modules/Weighter.$(ObjSuf)
 
@@ -1569 +1615 @@
         tmp/external/tcl/tclVar.$(ObjSuf)
 
+modules/TrackSmearing.h: \
+        classes/DelphesModule.h
+        @touch $@
+
 external/fastjet/ClusterSequence.hh: \
         external/fastjet/PseudoJet.hh \
@@ -1829 +1879 @@
         @touch $@
 
+modules/VertexSorter.h: \
+        classes/DelphesModule.h \
+        classes/DelphesClasses.h
+        @touch $@
+
 modules/Delphes.h: \
         classes/DelphesModule.h
+        @touch $@
+
+modules/VertexFinder.h: \
+        classes/DelphesModule.h \
+        classes/DelphesClasses.h
         @touch $@
 
@@ -2001 +2061 @@
 
 modules/FastJetFinder.h: \
+        classes/DelphesModule.h
+        @touch $@
+
+modules/BeamSpotFilter.h: \
         classes/DelphesModule.h
         @touch $@

cards/converter_card.tcl

@@ -13 +13 @@
 module TreeWriter TreeWriter {
 # add Branch InputArray BranchName BranchClass
-  add Branch Delphes/allParticles Particle GenParticle
+  add Branch Delphes/stableParticles Particle GenParticle
 }
 

cards/delphes_card_CMS.tcl

@@ -35 +35 @@
   
   FastJetFinder
+  FastCaloJetFinder
 
   JetEnergyScale
@@ -210 +211 @@
   set HCalEnergyMin 1.0
 
-  set ECalEnergySignificanceMin 1.0
-  set HCalEnergySignificanceMin 1.0
+  set ECalEnergySignificanceMin 2.0
+  set HCalEnergySignificanceMin 2.0
 
   set SmearTowerCenter true
@@ -497 +498 @@
 }
 
+################
+# caloJet finder
+################
+
+module FastJetFinder FastCaloJetFinder {
+  set InputArray Calorimeter/towers
+
+  set OutputArray jets
+
+  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 6
+  set ParameterR 0.5
+
+  set JetPTMin 20.0
+}
+
 ##################
 # Jet Energy Scale
@@ -610 +627 @@
 
   add Branch UniqueObjectFinder/jets Jet Jet
+  add Branch FastCaloJetFinder/jets CaloJet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon

classes/DelphesClasses.cc

@@ -41 +41 @@
 CompBase *Tower::fgCompare = CompE<Tower>::Instance();
 CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
+CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
 CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
 
@@ -121 +122 @@
   Charge(0), Mass(0.0),
   IsPU(0), IsRecoPU(0), IsConstituent(0), IsFromConversion(0),
+  ClusterIndex(-1), ClusterNDF(0), ClusterSigma(0), SumPT2(0), BTVSumPT2(0), GenDeltaZ(0), GenSumPT2(0),
   Flavor(0), FlavorAlgo(0), FlavorPhys(0),
   BTag(0), BTagAlgo(0), BTagPhys(0),
@@ -127 +129 @@
   Momentum(0.0, 0.0, 0.0, 0.0),
   Position(0.0, 0.0, 0.0, 0.0),
+  PositionError(0.0, 0.0, 0.0, 0.0),
+  InitialPosition(0.0, 0.0, 0.0, 0.0),
   Area(0.0, 0.0, 0.0, 0.0),
-  Dxy(0), SDxy(0), Xd(0), Yd(0), Zd(0),
+  L(0),
+  D0(0), ErrorD0(0), 
+  DZ(0), ErrorDZ(0), 
+  P(0),  ErrorP(0), 
+  PT(0), ErrorPT(0), 
+  CtgTheta(0), ErrorCtgTheta(0), 
+  Phi(0), ErrorPhi(0),  
+  Xd(0), Yd(0), Zd(0), 
   TrackResolution(0),
   NCharged(0),
@@ -245 +256 @@
   object.IsConstituent = IsConstituent;
   object.IsFromConversion = IsFromConversion;
+  object.ClusterIndex = ClusterIndex;
+  object.ClusterNDF = ClusterNDF;
+  object.ClusterSigma = ClusterSigma;
+  object.SumPT2 = SumPT2;
+  object.BTVSumPT2 = BTVSumPT2;
+  object.GenDeltaZ = GenDeltaZ;
+  object.GenSumPT2 = GenSumPT2;
   object.Flavor = Flavor;
   object.FlavorAlgo = FlavorAlgo;
@@ -262 +280 @@
   object.Momentum = Momentum;
   object.Position = Position;
+  object.InitialPosition = InitialPosition;
+  object.PositionError = PositionError;
   object.Area = Area;
-  object.Dxy = Dxy;
-  object.SDxy = SDxy;
+  object.L = L;
+  object.ErrorT = ErrorT;
+  object.D0 = D0;
+  object.ErrorD0 = ErrorD0;
+  object.DZ = DZ;
+  object.ErrorDZ = ErrorDZ;
+  object.P = P;
+  object.ErrorP = ErrorP;
+  object.PT = PT;
+  object.ErrorPT = ErrorPT;
+  object.CtgTheta = CtgTheta ;
+  object.ErrorCtgTheta = ErrorCtgTheta;
+  object.Phi = Phi;
+  object.ErrorPhi = ErrorPhi;  
   object.Xd = Xd;
   object.Yd = Yd;
@@ -282 +314 @@
   object.SumPtChargedPU = SumPtChargedPU;
   object.SumPt = SumPt;
-
+  object.ClusterIndex = ClusterIndex;
+  object.ClusterNDF = ClusterNDF;
+  object.ClusterSigma = ClusterSigma; 
+  object.SumPT2 = SumPT2;
+  
   object.FracPt[0] = FracPt[0];
   object.FracPt[1] = FracPt[1];
@@ -363 +399 @@
   Momentum.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
-  Dxy = 0.0;
-  SDxy = 0.0;
+  L = 0.0;
+  ErrorT = 0.0;
+  D0 = 0.0;  
+  ErrorD0 = 0.0;
+  DZ = 0.0;
+  ErrorDZ = 0.0;
+  P =0.0;
+  ErrorP =0.0;
+  PT = 0.0;
+  ErrorPT = 0.0;
+  CtgTheta = 0.0;
+  ErrorCtgTheta = 0.0;
+  Phi = 0.0;
+  ErrorPhi = 0.0;
   Xd = 0.0;
   Yd = 0.0;
@@ -387 +436 @@
   SumPt = -999;
 
+  ClusterIndex = -1;
+  ClusterNDF = -99;
+  ClusterSigma = 0.0; 
+  SumPT2 = 0.0;
+  BTVSumPT2 = 0.0;
+  GenDeltaZ = 0.0;
+  GenSumPT2 = 0.0; 
+  
   FracPt[0] = 0.0;
   FracPt[1] = 0.0;

classes/DelphesClasses.h

@@ -147 +147 @@
   Float_t Pz; // particle momentum vector (z component) | hepevt.phep[number][2]
 
-  Float_t PT; // particle transverse momentum
+  Float_t D0;
+  Float_t DZ;
+  Float_t P;
+  Float_t PT;
+  Float_t CtgTheta;
+  Float_t Phi;
   Float_t Eta; // particle pseudorapidity
-  Float_t Phi; // particle azimuthal angle
-
   Float_t Rapidity; // particle rapidity
 
@@ -168 +171 @@
 //---------------------------------------------------------------------------
 
-class Vertex: public TObject
+class Vertex: public SortableObject
 {
 public:
@@ -176 +179 @@
   Float_t Z; // vertex position (z component)
 
-  ClassDef(Vertex, 1)
+  Double_t ErrorX;
+  Double_t ErrorY;
+  Double_t ErrorZ;
+  Double_t ErrorT;
+
+  Int_t Index;
+  Int_t NDF;
+  Double_t Sigma;
+  Double_t SumPT2;
+  Double_t BTVSumPT2;
+  Double_t GenDeltaZ;
+  Double_t GenSumPT2;
+
+  TRefArray Constituents; // references to constituents
+
+  static CompBase *fgCompare; //!
+  const CompBase *GetCompare() const { return fgCompare; }
+
+  ClassDef(Vertex, 3)
 };
 
@@ -397 +418 @@
   Int_t Charge; // track charge
 
-  Float_t PT; // track transverse momentum
-
   Float_t Eta; // track pseudorapidity
-  Float_t Phi; // track azimuthal angle
 
   Float_t EtaOuter; // track pseudorapidity at the tracker edge
@@ -415 +433 @@
   Float_t TOuter; // track position (z component) at the tracker edge
 
-  Float_t Dxy;     // track signed transverse impact parameter
-  Float_t SDxy;    // signed error on the track signed transverse impact parameter
+  Float_t L; // track path length
+  Float_t ErrorT; // error on the time measurement
+
+  Float_t D0;     // track signed transverse impact parameter
+  Float_t ErrorD0;    // signed error on the track signed transverse impact parameter
+
+  Float_t DZ; // track transverse momentum
+  Float_t ErrorDZ; // track transverse momentum error
+
+  Float_t P; // track transverse momentum
+  Float_t ErrorP; // track transverse momentum error
+
+  Float_t PT; // track transverse momentum
+  Float_t ErrorPT; // track transverse momentum error
+
+  Float_t CtgTheta; // track transverse momentum
+  Float_t ErrorCtgTheta; // track transverse momentum error
+
+  Float_t Phi; // track azimuthal angle
+  Float_t ErrorPhi; // track azimuthal angle
+
   Float_t Xd;      // X coordinate of point of closest approach to vertex
   Float_t Yd;      // Y coordinate of point of closest approach to vertex
@@ -423 +460 @@
   TRef Particle; // reference to generated particle
 
+  Int_t VertexIndex; // reference to vertex
+
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
@@ -526 +565 @@
   Float_t DeltaPhi;
 
-  TLorentzVector Momentum, Position, Area;
-
-  Float_t Dxy;
-  Float_t SDxy;
+  TLorentzVector Momentum, Position, InitialPosition, PositionError, Area;
+
+  Float_t L; // path length
+  Float_t ErrorT; // path length
+  Float_t D0;
+  Float_t ErrorD0;
+  Float_t DZ;
+  Float_t ErrorDZ;
+  Float_t P;
+  Float_t ErrorP;
+  Float_t PT;
+  Float_t ErrorPT;
+  Float_t CtgTheta;
+  Float_t ErrorCtgTheta;
+  Float_t Phi;
+  Float_t ErrorPhi;
+
   Float_t Xd;
   Float_t Yd;
@@ -535 +587 @@
 
   // tracking resolution
-  
+
   Float_t TrackResolution;
 
@@ -562 +614 @@
   Float_t SumPt;
 
+  // vertex variables
+
+  Int_t ClusterIndex;
+  Int_t ClusterNDF;
+  Double_t ClusterSigma;
+  Double_t SumPT2;
+  Double_t BTVSumPT2;
+  Double_t GenDeltaZ;
+  Double_t GenSumPT2;
+
   // N-subjettiness variables
 
@@ -595 +657 @@
   void SetFactory(DelphesFactory *factory) { fFactory = factory; }
 
-  ClassDef(Candidate, 4)
+  ClassDef(Candidate, 5)
 };
 

classes/SortableObject.h

@@ -156 +156 @@
       return -1;
     else if(t1->ET < t2->ET)
+      return 1;
+    else
+      return 0;
+  }
+};
+
+//---------------------------------------------------------------------------
+
+template <typename T>
+class CompSumPT2: public CompBase
+{
+  CompSumPT2() {}
+public:
+  static CompSumPT2 *Instance()
+  {
+    static CompSumPT2 single;
+    return &single;
+  }
+
+  Int_t Compare(const TObject *obj1, const TObject *obj2) const
+  {
+    const T *t1 = static_cast<const T*>(obj1);
+    const T *t2 = static_cast<const T*>(obj2);
+    if(t1->SumPT2 > t2->SumPT2)
+      return -1;
+    else if(t1->SumPT2 < t2->SumPT2)
       return 1;
     else

doc/genMakefile.tcl

@@ -263 +263 @@
 executableDeps {converters/*.cpp} {examples/*.cpp}
 
-executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp}
+executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp} {readers/DelphesROOT.cpp}
 
 puts {ifeq ($(HAS_CMSSW),true)}

examples/Pythia8/configParticleGun.cmnd

@@ -3 +3 @@
 ! 1) Settings used in the main program.
 
-Main:numberOfEvents = 10000         ! number of events to generate
+Main:numberOfEvents = 100000         ! number of events to generate
 Main:timesAllowErrors = 3          ! how many aborts before run stops
 Main:spareFlag1 = on                ! true means particle gun
-Main:spareMode1 = 11               ! 1-5 - di-quark, 21 - di-gluon, 11 - single electron, 13 - single muon, 22 - single photon
+Main:spareMode1 = ID               ! 1-5 - di-quark, 21 - di-gluon, 11 - single electron, 13 - single muon, 22 - single photon
 Main:spareParm1 = 10000           ! max pt
 

modules/Calorimeter.cc

@@ -58 +58 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  Int_t i;
-
+  
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
 
-  for(i = 0; i < 2; ++i)
-  {
-    fECalTowerTrackArray[i] = new TObjArray;
-    fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
-
-    fHCalTowerTrackArray[i] = new TObjArray;
-    fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
-  }
+  fECalTowerTrackArray = new TObjArray;
+  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
+
+  fHCalTowerTrackArray = new TObjArray;
+  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
+  
 }
 
@@ -77 +74 @@
 Calorimeter::~Calorimeter()
 {
-  Int_t i;
-
+  
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
 
-  for(i = 0; i < 2; ++i)
-  {
-    if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
-    if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
-
-    if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
-    if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
-  }
+  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
+  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
+
+  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
+  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
+ 
 }
 
@@ -218 +212 @@
   Double_t ecalFraction, hcalFraction;
   Double_t ecalEnergy, hcalEnergy;
-  Double_t ecalSigma, hcalSigma;
   Int_t pdgCode;
 
@@ -368 +361 @@
       fHCalTowerEnergy = 0.0;
 
-      fECalTrackEnergy[0] = 0.0;
-      fECalTrackEnergy[1] = 0.0;
-
-      fHCalTrackEnergy[0] = 0.0;
-      fHCalTrackEnergy[1] = 0.0;
-
+      fECalTrackEnergy = 0.0;
+      fHCalTrackEnergy = 0.0;
+      
+      fECalTrackSigma = 0.0;
+      fHCalTrackSigma = 0.0;
+      
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
 
-      fECalTowerTrackArray[0]->Clear();
-      fECalTowerTrackArray[1]->Clear();
-
-      fHCalTowerTrackArray[0]->Clear();
-      fHCalTowerTrackArray[1]->Clear();
+      fECalTowerTrackArray->Clear();
+      fHCalTowerTrackArray->Clear();
+    
     }
 
@@ -406 +397 @@
       if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
-        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(ecalSigma/momentum.E() < track->TrackResolution)
-        {
-          fECalTrackEnergy[0] += ecalEnergy;
-          fECalTowerTrackArray[0]->Add(track);
-        }
-        else
-        {
-          fECalTrackEnergy[1] += ecalEnergy;
-          fECalTowerTrackArray[1]->Add(track);
-        }
+        fECalTrackEnergy += ecalEnergy;
+        fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
+        fECalTowerTrackArray->Add(track);
       }
+     
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
       {
-        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(hcalSigma/momentum.E() < track->TrackResolution)
-        {
-          fHCalTrackEnergy[0] += hcalEnergy;
-          fHCalTowerTrackArray[0]->Add(track);
-        }
-        else
-        {
-          fHCalTrackEnergy[1] += hcalEnergy;
-          fHCalTowerTrackArray[1]->Add(track);
-        }
+        fHCalTrackEnergy += hcalEnergy;
+        fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
+        fHCalTowerTrackArray->Add(track);
       }
+      
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
@@ -476 +453 @@
   Double_t energy, pt, eta, phi;
   Double_t ecalEnergy, hcalEnergy;
+  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
+  
   Double_t ecalSigma, hcalSigma;
-
+  Double_t ecalNeutralSigma, hcalNeutralSigma;
+
+  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
+  
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
@@ -484 +466 @@
 
   if(!fTower) return;
-
 
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
@@ -554 +535 @@
     fTowerOutputArray->Add(fTower);
   }
-
+  
   // fill energy flow candidates
-
-  ecalEnergy -= fECalTrackEnergy[1];
-  hcalEnergy -= fHCalTrackEnergy[1];
-
-  fItECalTowerTrackArray[0]->Reset();
-  while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  fItECalTowerTrackArray[1]->Reset();
-  while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  fItHCalTowerTrackArray[0]->Reset();
-  while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  fItHCalTowerTrackArray[1]->Reset();
-  while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
-  if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0;
-
-  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
-  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
-
-  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
-  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
-
-  energy = ecalEnergy + hcalEnergy;
-
-  if(ecalEnergy > 0.0)
+  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
+  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
+
+  //compute neutral excesses
+  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
+  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
+  
+  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
+  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
+  
+   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
+  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-
-    pt = ecalEnergy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
-    tower->Eem = ecalEnergy;
+    pt =  ecalNeutralEnergy / TMath::CosH(eta);
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
+    tower->Eem = ecalNeutralEnergy;
     tower->Ehad = 0.0;
     tower->PID = 22;
-
+    
     fEFlowPhotonOutputArray->Add(tower);
-  }
-  if(hcalEnergy > 0.0)
-  {
-    // create new neutral hadron tower
+   
+    //clone tracks
+    fItECalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  
+  }
+ 
+  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
+  else if(fECalTrackEnergy > 0.0)
+  {
+    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
+    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
+
+    //rescale tracks
+    fItECalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
+    {  
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+
+
+  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
+  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
+  {
+    // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-
-    pt = hcalEnergy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
+    pt =  hcalNeutralEnergy / TMath::CosH(eta);
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
+    tower->Ehad = hcalNeutralEnergy;
     tower->Eem = 0.0;
-    tower->Ehad = hcalEnergy;
-
+    
     fEFlowNeutralHadronOutputArray->Add(tower);
-  }
+   
+    //clone tracks
+    fItHCalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  
+  }
+ 
+  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
+  else if(fHCalTrackEnergy > 0.0)
+  {
+    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
+    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
+
+    //rescale tracks
+    fItHCalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
+    {  
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+  
+  
 }
 

modules/Calorimeter.h

@@ -58 +58 @@
   Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
   Double_t fECalTowerEnergy, fHCalTowerEnergy;
-  Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
+  Double_t fECalTrackEnergy, fHCalTrackEnergy;
 
   Double_t fTimingEnergyMin;
@@ -70 +70 @@
   Double_t fECalEnergySignificanceMin;
   Double_t fHCalEnergySignificanceMin;
+
+  Double_t fECalTrackSigma;
+  Double_t fHCalTrackSigma;
 
   Bool_t fSmearTowerCenter;
@@ -103 +106 @@
   TObjArray *fEFlowNeutralHadronOutputArray; //!
 
-  TObjArray *fECalTowerTrackArray[2]; //!
-  TIterator *fItECalTowerTrackArray[2]; //!
+  TObjArray *fECalTowerTrackArray; //!
+  TIterator *fItECalTowerTrackArray; //!
 
-  TObjArray *fHCalTowerTrackArray[2]; //!
-  TIterator *fItHCalTowerTrackArray[2]; //!
+  TObjArray *fHCalTowerTrackArray; //!
+  TIterator *fItHCalTowerTrackArray; //!
 
   void FinalizeTower();

modules/ImpactParameterSmearing.cc

@@ -96 +96 @@
 {
   Candidate *candidate, *particle, *mother;
-  Double_t xd, yd, zd, dxy, sx, sy, sz, ddxy;
+  Double_t xd, yd, zd, d0, sx, sy, sz, dd0;
   Double_t pt, eta, px, py, phi, e;
 
@@ -103 +103 @@
   {
 
-    // take momentum before smearing (otherwise apply double smearing on dxy)
+    // take momentum before smearing (otherwise apply double smearing on d0)
     particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
 
@@ -131 +131 @@
 
     // calculate impact parameter (after-smearing)
-    dxy = (xd*py - yd*px)/pt;
+    d0 = (xd*py - yd*px)/pt;
 
-    ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
+    dd0 = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
 
     // fill smeared values in candidate
@@ -143 +143 @@
     candidate->Zd = zd;
 
-    candidate->Dxy = dxy;
-    candidate->SDxy = ddxy;
+    candidate->D0 = d0;
+    candidate->ErrorD0 = dd0;
 
     candidate->AddCandidate(mother);

modules/ModulesLinkDef.h

@@ -35 +35 @@
 #include "modules/EnergySmearing.h"
 #include "modules/MomentumSmearing.h"
+#include "modules/TrackSmearing.h"
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
@@ -58 +59 @@
 #include "modules/StatusPidFilter.h"
 #include "modules/PdgCodeFilter.h"
+#include "modules/BeamSpotFilter.h"
 #include "modules/RecoPuFilter.h"
 #include "modules/Cloner.h"
@@ -64 +66 @@
 #include "modules/JetFlavorAssociation.h"
 #include "modules/JetFakeParticle.h"
+#include "modules/VertexSorter.h"
+#include "modules/VertexFinder.h"
+#include "modules/VertexFinderDA4D.h"
 #include "modules/ExampleModule.h"
 
@@ -81 +86 @@
 #pragma link C++ class EnergySmearing+;
 #pragma link C++ class MomentumSmearing+;
+#pragma link C++ class TrackSmearing+;
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
@@ -104 +110 @@
 #pragma link C++ class StatusPidFilter+;
 #pragma link C++ class PdgCodeFilter+;
+#pragma link C++ class BeamSpotFilter+;
 #pragma link C++ class RecoPuFilter+;
 #pragma link C++ class Cloner+;
@@ -110 +117 @@
 #pragma link C++ class JetFlavorAssociation+;
 #pragma link C++ class JetFakeParticle+;
+#pragma link C++ class VertexSorter+;
+#pragma link C++ class VertexFinder+;
+#pragma link C++ class VertexFinderDA4D+;
 #pragma link C++ class ExampleModule+;
 

modules/ParticlePropagator.cc

@@ -67 +67 @@
 }
 
+
 //------------------------------------------------------------------------------
 
@@ -91 +92 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // import beamspot
+  try
+  {
+    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
+  }
+  catch(runtime_error &e)
+  {
+    fBeamSpotInputArray = 0;
+  }
   // create output arrays
 
@@ -111 +121 @@
 {
   Candidate *candidate, *mother;
-  TLorentzVector candidatePosition, candidateMomentum;
+  TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
   Double_t px, py, pz, pt, pt2, e, q;
   Double_t x, y, z, t, r, phi;
@@ -120 +130 @@
   Double_t tmp, discr, discr2;
   Double_t delta, gammam, omega, asinrho;
-  Double_t rcu, rc2, dxy, xd, yd, zd;
+  Double_t rcu, rc2, xd, yd, zd;
+  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
+  Double_t bsx, bsy, bsz;
 
   const Double_t c_light = 2.99792458E8;
+
+  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
+    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  else
+  {
+    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
+    beamSpotPosition = beamSpotCandidate.Position;
+  }
 
   fItInputArray->Reset();
@@ -132 +152 @@
     y = candidatePosition.Y()*1.0E-3;
     z = candidatePosition.Z()*1.0E-3;
+
+    bsx = beamSpotPosition.X()*1.0E-3;
+    bsy = beamSpotPosition.Y()*1.0E-3;
+    bsz = beamSpotPosition.Z()*1.0E-3;
+
     q = candidate->Charge;
 
@@ -182 +207 @@
       z_t = z + pz*t;
 
+      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
+
       mother = candidate;
       candidate = static_cast<Candidate*>(candidate->Clone());
 
+      candidate->InitialPosition = candidatePosition;
       candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
+      candidate->L = l*1.0E3;
 
       candidate->Momentum = candidateMomentum;
@@ -239 +268 @@
       zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
 
-      // calculate impact paramater
-      dxy = (xd*py - yd*px)/pt;
+      // use perigee momentum rather than original particle
+      // momentum, since the orignal particle momentum isn't known
+
+      px = TMath::Sign(1.0,r) * pt * (-y_c / r_c);
+      py = TMath::Sign(1.0,r) * pt * (x_c / r_c);
+      etap = candidateMomentum.Eta();
+      phip = TMath::ATan2(py, px);
+
+      candidateMomentum.SetPtEtaPhiE(pt, etap, phip, candidateMomentum.E());
+
+      // calculate additional track parameters (correct for beamspot position)
+
+      d0        = (  (x - bsx) * py - (y - bsy) * px) / pt;
+      dz        = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
+      p         = candidateMomentum.P();
+      ctgTheta  = 1.0 / TMath::Tan (candidateMomentum.Theta ());
+
 
       // 3. time evaluation t = TMath::Min(t_r, t_z)
@@ -287 +331 @@
       r_t = TMath::Hypot(x_t, y_t);
 
+
+      // compute path length for an helix
+
+      alpha = pz*1.0E9 / c_light / gammam;
+      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
+
       if(r_t > 0.0)
       {
+
+        // store these variables before cloning
+        candidate->D0 = d0*1.0E3;
+        candidate->DZ = dz*1.0E3;
+        candidate->P  = p;
+        candidate->PT = pt;
+        candidate->CtgTheta = ctgTheta;
+        candidate->Phi = phip;
+
         mother = candidate;
         candidate = static_cast<Candidate*>(candidate->Clone());
 
+        candidate->InitialPosition = candidatePosition;
         candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
 
         candidate->Momentum = candidateMomentum;
-        candidate->Dxy = dxy*1.0E3;
-        candidate->Xd = xd*1.0E3;
+
+            candidate->L  =  l*1.0E3;
+
+            candidate->Xd = xd*1.0E3;
         candidate->Yd = yd*1.0E3;
         candidate->Zd = zd*1.0E3;

modules/ParticlePropagator.h

@@ -35 +35 @@
 class TClonesArray;
 class TIterator;
+class TLorentzVector;
 
 class ParticlePropagator: public DelphesModule
@@ -55 +56 @@
 
   const TObjArray *fInputArray; //!
+  const TObjArray *fBeamSpotInputArray; //!
 
   TObjArray *fOutputArray; //!

modules/PileUpMerger.cc

@@ -115 +115 @@
   TDatabasePDG *pdg = TDatabasePDG::Instance();
   TParticlePDG *pdgParticle;
-  Int_t pid;
+  Int_t pid, nch, nvtx = -1;
   Float_t x, y, z, t, vx, vy;
-  Float_t px, py, pz, e;
-  Double_t dz, dphi, dt;
+  Float_t px, py, pz, e, pt;
+  Double_t dz, dphi, dt, sumpt2, dz0, dt0;
   Int_t numberOfEvents, event, numberOfParticles;
   Long64_t allEntries, entry;
@@ -132 +132 @@
   fFunction->GetRandom2(dz, dt);
 
+  dz0 = -1.0e6;
+  dt0 = -1.0e6;
+
   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
   dz *= 1.0E3; // necessary in order to make z in mm
+
+  //cout<<dz<<","<<dt<<endl;
+
   vx = 0.0;
   vy = 0.0;
+
   numberOfParticles = fInputArray->GetEntriesFast();
+  nch = 0;
+  sumpt2 = 0.0;
+
+  factory = GetFactory();
+  vertex = factory->NewCandidate();
+
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
@@ -143 +156 @@
     z = candidate->Position.Z();
     t = candidate->Position.T();
-    candidate->Position.SetZ(z + dz);
-    candidate->Position.SetT(t + dt);
+    pt = candidate->Momentum.Pt();
+
+    // take postion and time from first stable particle
+    if (dz0 < -999999.0)
+      dz0 = z;
+    if (dt0 < -999999.0)
+      dt0 = t;
+
+    // cancel any possible offset in position and time the input file
+    candidate->Position.SetZ(z - dz0 + dz);
+    candidate->Position.SetT(t - dt0 + dt);
+
+    candidate->IsPU = 0;
+
     fParticleOutputArray->Add(candidate);
+
+    if(TMath::Abs(candidate->Charge) >  1.0E-9)
+    {
+      nch++;
+      sumpt2 += pt*pt;
+      vertex->AddCandidate(candidate);
+    }
   }
 
   if(numberOfParticles > 0)
   {
-    vx /= numberOfParticles;
-    vy /= numberOfParticles;
+    vx /= sumpt2;
+    vy /= sumpt2;
   }
 
-  factory = GetFactory();
-
-  vertex = factory->NewCandidate();
+  nvtx++;
   vertex->Position.SetXYZT(vx, vy, dz, dt);
+  vertex->ClusterIndex = nvtx;
+  vertex->ClusterNDF = nch;
+  vertex->SumPT2 = sumpt2;
+  vertex->GenSumPT2 = sumpt2;
   fVertexOutputArray->Add(vertex);
 
@@ -170 +204 @@
       numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
       break;
+    case 2:
+      numberOfEvents = fMeanPileUp;
+      break;
     default:
       numberOfEvents = gRandom->Poisson(fMeanPileUp);
@@ -176 +213 @@
 
   allEntries = fReader->GetEntries();
+
 
   for(event = 0; event < numberOfEvents; ++event)
@@ -198 +236 @@
     vx = 0.0;
     vy = 0.0;
+
     numberOfParticles = 0;
+    sumpt2 = 0.0;
+
+    //factory = GetFactory();
+    vertex = factory->NewCandidate();
+
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     {
@@ -215 +259 @@
       candidate->Momentum.SetPxPyPzE(px, py, pz, e);
       candidate->Momentum.RotateZ(dphi);
+      pt = candidate->Momentum.Pt();
 
       x -= fInputBeamSpotX;
@@ -224 +269 @@
       vx += candidate->Position.X();
       vy += candidate->Position.Y();
+
       ++numberOfParticles;
+      if(TMath::Abs(candidate->Charge) >  1.0E-9)
+      {
+        nch++;
+        sumpt2 += pt*pt;
+        vertex->AddCandidate(candidate);
+      }
 
       fParticleOutputArray->Add(candidate);
@@ -235 +287 @@
     }
 
-    vertex = factory->NewCandidate();
+    nvtx++;
+
     vertex->Position.SetXYZT(vx, vy, dz, dt);
+
+    vertex->ClusterIndex = nvtx;
+    vertex->ClusterNDF = nch;
+    vertex->SumPT2 = sumpt2;
+    vertex->GenSumPT2 = sumpt2;
+
     vertex->IsPU = 1;
 
     fVertexOutputArray->Add(vertex);
+
   }
 }

modules/SimpleCalorimeter.cc

@@ -58 +58 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  Int_t i;
-
+  
   fResolutionFormula = new DelphesFormula;
-
-  for(i = 0; i < 2; ++i)
-  {
-    fTowerTrackArray[i] = new TObjArray;
-    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
-  }
+  fTowerTrackArray = new TObjArray;
+  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
+ 
 }
 
@@ -73 +69 @@
 SimpleCalorimeter::~SimpleCalorimeter()
 {
-  Int_t i;
-
+ 
   if(fResolutionFormula) delete fResolutionFormula;
-
-  for(i = 0; i < 2; ++i)
-  {
-    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
-    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
-  }
+  if(fTowerTrackArray) delete fTowerTrackArray;
+  if(fItTowerTrackArray) delete fItTowerTrackArray;
+ 
 }
 
@@ -198 +190 @@
   Double_t fraction;
   Double_t energy;
-  Double_t sigma;
   Int_t pdgCode;
 
@@ -340 +331 @@
       fTowerEnergy = 0.0;
 
-      fTrackEnergy[0] = 0.0;
-      fTrackEnergy[1] = 0.0;
+      fTrackEnergy = 0.0;
+      fTrackSigma = 0.0;
 
       fTowerTime = 0.0;
@@ -351 +342 @@
       fTowerPhotonHits = 0;
 
-      fTowerTrackArray[0]->Clear();
-      fTowerTrackArray[1]->Clear();
-    }
+      fTowerTrackArray->Clear();
+     }
 
     // check for track hits
@@ -371 +361 @@
       if(fTrackFractions[number] > 1.0E-9)
       {
-        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(sigma/momentum.E() < track->TrackResolution)
-        {
-          fTrackEnergy[0] += energy;
-          fTowerTrackArray[0]->Add(track);
-        }
-        else
-        {
-          fTrackEnergy[1] += energy;
-          fTowerTrackArray[1]->Add(track);
-        }
+             
+       // compute total charged energy   
+       fTrackEnergy += energy;
+       fTrackSigma += ((track->TrackResolution)*momentum.E())*((track->TrackResolution)*momentum.E());
+       
+       fTowerTrackArray->Add(track);
+      
       }
+       
       else
       {
@@ -403 +390 @@
     fTowerEnergy += energy;
 
-    fTowerTime += TMath::Sqrt(energy)*position.T();
-    fTowerTimeWeight += TMath::Sqrt(energy);
+    fTowerTime += energy*position.T();
+    fTowerTimeWeight += energy;
 
     fTower->AddCandidate(particle);
@@ -418 +405 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy, pt, eta, phi;
-  Double_t sigma;
+  Double_t energy,neutralEnergy, pt, eta, phi;
+  Double_t sigma, neutralSigma;
   Double_t time;
+   
+  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
 
   TLorentzVector momentum;
@@ -436 +425 @@
 
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
 
   if(fSmearTowerCenter)
@@ -464 +454 @@
   if(energy > 0.0) fTowerOutputArray->Add(fTower);
 
-  // fill e-flow candidates
-
-  energy -= fTrackEnergy[1];
-
-  fItTowerTrackArray[0]->Reset();
-  while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    track->Momentum *= energy/fTrackEnergy[0];
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  fItTowerTrackArray[1]->Reset();
-  while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
-  {
-    mother = track;
-    track = static_cast<Candidate*>(track->Clone());
-    track->AddCandidate(mother);
-
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
-
-  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
-
-  // save energy excess as an energy flow tower
-  if(energy > 0.0)
+ 
+  // e-flow candidates
+
+  //compute neutral excess
+  
+  fTrackSigma = TMath::Sqrt(fTrackSigma);
+  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
+  
+  //compute sigma_trk total
+  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
+   
+  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
+  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt = energy / TMath::CosH(eta);
-
-    tower->Eem = (!fIsEcal) ? 0 : energy;
-    tower->Ehad = (fIsEcal) ? 0 : energy;
+    pt = neutralEnergy / TMath::CosH(eta);
+
+    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
+    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
+    tower->PID = (fIsEcal) ? 22 : 0;
+ 
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
+    fEFlowTowerOutputArray->Add(tower);
     
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    fItTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
     
-    tower->PID = (fIsEcal) ? 22 : 0;
-    
-    fEFlowTowerOutputArray->Add(tower);
-  }
-}
+  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
+  else if (fTrackEnergy > 0.0)
+  {
+    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
+    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
+   
+    fItTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+   
+ } 
 
 //------------------------------------------------------------------------------

modules/SimpleCalorimeter.h

@@ -59 +59 @@
   Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
   Double_t fTowerEnergy;
-  Double_t fTrackEnergy[2];
+  Double_t fTrackEnergy;
 
   Double_t fTowerTime;
@@ -72 +72 @@
 
   Double_t fEnergySignificanceMin;
+
+  Double_t fTrackSigma;
 
   Bool_t fSmearTowerCenter;
@@ -102 +104 @@
   TObjArray *fEFlowTowerOutputArray; //!
 
-  TObjArray *fTowerTrackArray[2]; //!
-  TIterator *fItTowerTrackArray[2]; //!
+  TObjArray *fTowerTrackArray; //!
+  TIterator *fItTowerTrackArray; //!
 
   void FinalizeTower();

modules/TimeSmearing.cc

@@ -93 +93 @@
 {
   Candidate *candidate, *mother;
-  Double_t t;
+  Double_t ti, tf_smeared, tf;
   const Double_t c_light = 2.99792458E8;
 
@@ -99 +99 @@
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
-    t = candidatePosition.T()*1.0E-3/c_light;
+    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
+    const TLorentzVector &candidateFinalPosition = candidate->Position;
+
+    ti = candidateInitialPosition.T()*1.0E-3/c_light;
+    tf = candidateFinalPosition.T()*1.0E-3/c_light;
 
     // apply smearing formula
-    t = gRandom->Gaus(t, fTimeResolution);
+    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
+    ti = ti + tf_smeared - tf;
 
     mother = candidate;
     candidate = static_cast<Candidate*>(candidate->Clone());
-    candidate->Position.SetT(t*1.0E3*c_light);
+    candidate->InitialPosition.SetT(ti*1.0E3*c_light);
+    candidate->Position.SetT(tf*1.0E3*c_light);
+
+    candidate->ErrorT = fTimeResolution*1.0E3*c_light;
 
     candidate->AddCandidate(mother);

modules/TrackCountingBTagging.cc

@@ -96 +96 @@
 
   Double_t jpx, jpy;
-  Double_t dr, tpx, tpy, tpt;
-  Double_t xd, yd, dxy, ddxy, ip, sip;
+  Double_t dr, tpt;
+  Double_t xd, yd, d0, dd0, ip, sip;
 
   Int_t sign;
@@ -117 +117 @@
     {
       const TLorentzVector &trkMomentum = track->Momentum;
-
+      
       dr = jetMomentum.DeltaR(trkMomentum);
-
       tpt = trkMomentum.Pt();
-      tpx = trkMomentum.Px();
-      tpy = trkMomentum.Py();
-
       xd = track->Xd;
       yd = track->Yd;
-      dxy = TMath::Hypot(xd, yd);
-      ddxy = track->SDxy;
+      d0 = TMath::Hypot(xd, yd);
+      dd0 = track->ErrorD0;
 
       if(tpt < fPtMin) continue;
       if(dr > fDeltaR) continue;
-      if(dxy > fIPmax) continue;
+      if(d0 > fIPmax) continue;
 
       sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
 
-      ip = sign*dxy;
-      sip = ip / TMath::Abs(ddxy);
+      ip = sign*d0;
+      sip = ip / TMath::Abs(dd0);
 
       if(sip > fSigMin) count++;

modules/TreeWriter.cc

@@ -215 +215 @@
     entry->Pz = momentum.Pz();
 
+    entry->D0            = candidate->D0;
+    entry->DZ            = candidate->DZ;
+    entry->P             = candidate->P;
+    entry->PT            = candidate->PT;
+    entry->CtgTheta      = candidate->CtgTheta;
+    entry->Phi           = candidate->Phi;
+
     entry->Eta = eta;
     entry->Phi = momentum.Phi();
@@ -233 +240 @@
 {
   TIter iterator(array);
-  Candidate *candidate = 0;
+  Candidate *candidate = 0, *constituent = 0;
   Vertex *entry = 0;
 
   const Double_t c_light = 2.99792458E8;
 
+  Double_t x, y, z, t, xError, yError, zError, tError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
+  UInt_t index, ndf;
+
+  CompBase *compare = Candidate::fgCompare;
+  Candidate::fgCompare = CompSumPT2<Candidate>::Instance();
+  array->Sort();
+  Candidate::fgCompare = compare;
+
   // loop over all vertices
   iterator.Reset();
   while((candidate = static_cast<Candidate*>(iterator.Next())))
   {
-    const TLorentzVector &position = candidate->Position;
+
+    index = candidate->ClusterIndex;
+    ndf = candidate->ClusterNDF;
+    sigma = candidate->ClusterSigma;
+    sumPT2 = candidate->SumPT2;
+    btvSumPT2 = candidate->BTVSumPT2;
+    genDeltaZ = candidate->GenDeltaZ;
+    genSumPT2 = candidate->GenSumPT2;
+
+    x = candidate->Position.X();
+    y = candidate->Position.Y();
+    z = candidate->Position.Z();
+    t = candidate->Position.T()*1.0E-3/c_light;
+
+    xError = candidate->PositionError.X ();
+    yError = candidate->PositionError.Y ();
+    zError = candidate->PositionError.Z ();
+    tError = candidate->PositionError.T ()*1.0E-3/c_light;
 
     entry = static_cast<Vertex*>(branch->NewEntry());
 
-    entry->X = position.X();
-    entry->Y = position.Y();
-    entry->Z = position.Z();
-    entry->T = position.T()*1.0E-3/c_light;
-  }
-}
+    entry->Index = index;
+    entry->NDF = ndf;
+    entry->Sigma = sigma;
+    entry->SumPT2 = sumPT2;
+    entry->BTVSumPT2 = btvSumPT2;
+    entry->GenDeltaZ = genDeltaZ;
+    entry->GenSumPT2 = genSumPT2;
+
+    entry->X = x;
+    entry->Y = y;
+    entry->Z = z;
+    entry->T = t;
+
+    entry->ErrorX = xError;
+    entry->ErrorY = yError;
+    entry->ErrorZ = zError;
+    entry->ErrorT = tError;
+
+
+    TIter itConstituents(candidate->GetCandidates());
+    itConstituents.Reset();
+    entry->Constituents.Clear();
+    while((constituent = static_cast<Candidate*>(itConstituents.Next())))
+    {
+      entry->Constituents.Add(constituent);
+    }
+
+  }
+}
+
 
 //------------------------------------------------------------------------------
@@ -261 +317 @@
   Candidate *particle = 0;
   Track *entry = 0;
-  Double_t pt, signz, cosTheta, eta, rapidity;
+  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
   const Double_t c_light = 2.99792458E8;
 
@@ -292 +348 @@
     entry->TOuter = position.T()*1.0E-3/c_light;
 
-    entry->Dxy = candidate->Dxy;
-    entry->SDxy = candidate->SDxy ;
+    entry->L = candidate->L;
+
+    entry->D0            = candidate->D0;
+    entry->ErrorD0       = candidate->ErrorD0;
+    entry->DZ            = candidate->DZ;
+    entry->ErrorDZ       = candidate->ErrorDZ;
+
+    entry->ErrorP        = candidate->ErrorP;
+    entry->ErrorPT       = candidate->ErrorPT;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+    entry->ErrorPhi      = candidate->ErrorPhi;
+
     entry->Xd = candidate->Xd;
     entry->Yd = candidate->Yd;
@@ -301 +367 @@
 
     pt = momentum.Pt();
+    p = momentum.P();
+    phi = momentum.Phi();
+    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1/TMath::Tan(momentum.Theta()) : 1e10;
+
     cosTheta = TMath::Abs(momentum.CosTheta());
     signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
@@ -306 +376 @@
     rapidity = (cosTheta == 1.0 ? signz*999.9 : momentum.Rapidity());
 
+    entry->PT  = pt;
     entry->Eta = eta;
-    entry->Phi = momentum.Phi();
-    entry->PT = pt;
+    entry->Phi = phi;
+    entry->CtgTheta = ctgTheta;
 
     particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
@@ -319 +390 @@
 
     entry->Particle = particle;
+
+    entry->VertexIndex = candidate->ClusterIndex;
+
   }
 }

readers/DelphesCMSFWLite.cpp

@@ -272 +272 @@
 
     branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
-    branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
+    branchRwgt = treeWriter->NewBranch("Weight", Weight::Class());
 
     confReader = new ExRootConfReader;