Changes in / [ec5e04b:21eab4f] in git

Files:

• Makefile (modified) (20 diffs)
• cards/converter_card.tcl (modified) (1 diff)
• cards/delphes_card_ATLAS.tcl (modified) (1 diff)
• cards/delphes_card_ATLAS_PileUp.tcl (modified) (1 diff)
• cards/delphes_card_CMS.tcl (modified) (10 diffs)
• cards/delphes_card_CMS_NoFastJet.tcl (modified) (2 diffs)
• cards/delphes_card_CMS_PileUp.tcl (modified) (1 diff)
• cards/validation_card.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.cpp (added)
• examples/validation.sh (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/MomentumSmearing.cc (modified) (4 diffs)
• modules/MomentumSmearing.h (modified) (1 diff)
• 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/DelphesPythia8.cpp (modified) (2 diffs)
• readers/DelphesROOT.cpp (added)

Makefile

@@ -137 +137 @@
 tmp/examples/Example1.$(ObjSuf): \
         examples/Example1.cpp \
+        classes/DelphesClasses.h \
+        external/ExRootAnalysis/ExRootTreeReader.h \
+        external/ExRootAnalysis/ExRootTreeWriter.h \
+        external/ExRootAnalysis/ExRootTreeBranch.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootUtilities.h
+Validation$(ExeSuf): \
+        tmp/examples/Validation.$(ObjSuf)
+
+tmp/examples/Validation.$(ObjSuf): \
+        examples/Validation.cpp \
         classes/DelphesClasses.h \
         external/ExRootAnalysis/ExRootTreeReader.h \
@@ -150 +161 @@
         root2pileup$(ExeSuf) \
         stdhep2pileup$(ExeSuf) \
-        Example1$(ExeSuf)
+        Example1$(ExeSuf) \
+        Validation$(ExeSuf)
 
 EXECUTABLE_OBJ +=  \
@@ -159 +171 @@
         tmp/converters/root2pileup.$(ObjSuf) \
         tmp/converters/stdhep2pileup.$(ObjSuf) \
-        tmp/examples/Example1.$(ObjSuf)
+        tmp/examples/Example1.$(ObjSuf) \
+        tmp/examples/Validation.$(ObjSuf)
 
 DelphesHepMC$(ExeSuf): \
@@ -183 +196 @@
         classes/DelphesLHEFReader.h \
         external/ExRootAnalysis/ExRootTreeWriter.h \
+        external/ExRootAnalysis/ExRootTreeBranch.h \
+        external/ExRootAnalysis/ExRootProgressBar.h
+DelphesROOT$(ExeSuf): \
+        tmp/readers/DelphesROOT.$(ObjSuf)
+
+tmp/readers/DelphesROOT.$(ObjSuf): \
+        readers/DelphesROOT.cpp \
+        modules/Delphes.h \
+        classes/DelphesStream.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        external/ExRootAnalysis/ExRootTreeWriter.h \
+        external/ExRootAnalysis/ExRootTreeReader.h \
         external/ExRootAnalysis/ExRootTreeBranch.h \
         external/ExRootAnalysis/ExRootProgressBar.h
@@ -200 +226 @@
         DelphesHepMC$(ExeSuf) \
         DelphesLHEF$(ExeSuf) \
+        DelphesROOT$(ExeSuf) \
         DelphesSTDHEP$(ExeSuf)
 
@@ -205 +232 @@
         tmp/readers/DelphesHepMC.$(ObjSuf) \
         tmp/readers/DelphesLHEF.$(ObjSuf) \
+        tmp/readers/DelphesROOT.$(ObjSuf) \
         tmp/readers/DelphesSTDHEP.$(ObjSuf)
 
@@ -319 +347 @@
         modules/EnergySmearing.h \
         modules/MomentumSmearing.h \
+        modules/TrackSmearing.h \
         modules/ImpactParameterSmearing.h \
         modules/TimeSmearing.h \
@@ -342 +371 @@
         modules/StatusPidFilter.h \
         modules/PdgCodeFilter.h \
+        modules/BeamSpotFilter.h \
+        modules/RecoPuFilter.h \
         modules/Cloner.h \
         modules/Weighter.h \
@@ -347 +378 @@
         modules/JetFlavorAssociation.h \
         modules/JetFakeParticle.h \
+        modules/VertexSorter.h \
+        modules/VertexFinder.h \
+        modules/VertexFinderDA4D.h \
         modules/ExampleModule.h
 tmp/modules/ModulesDict$(PcmSuf): \
@@ -553 +587 @@
         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) \
@@ -785 +828 @@
         external/ExRootAnalysis/ExRootFilter.h \
         external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/RecoPuFilter.$(ObjSuf): \
+        modules/RecoPuFilter.$(SrcSuf) \
+        modules/RecoPuFilter.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/SimpleCalorimeter.$(ObjSuf): \
         modules/SimpleCalorimeter.$(SrcSuf) \
@@ -846 +898 @@
         modules/TrackPileUpSubtractor.$(SrcSuf) \
         modules/TrackPileUpSubtractor.h \
+        classes/DelphesClasses.h \
+        classes/DelphesFactory.h \
+        classes/DelphesFormula.h \
+        external/ExRootAnalysis/ExRootResult.h \
+        external/ExRootAnalysis/ExRootFilter.h \
+        external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/TrackSmearing.$(ObjSuf): \
+        modules/TrackSmearing.$(SrcSuf) \
+        modules/TrackSmearing.h \
         classes/DelphesClasses.h \
         classes/DelphesFactory.h \
@@ -868 +929 @@
         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/VertexFinderDA4D.$(ObjSuf): \
+        modules/VertexFinderDA4D.$(SrcSuf) \
+        modules/VertexFinderDA4D.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 +1022 @@
         tmp/modules/AngularSmearing.$(ObjSuf) \
         tmp/modules/BTagging.$(ObjSuf) \
+        tmp/modules/BeamSpotFilter.$(ObjSuf) \
         tmp/modules/Calorimeter.$(ObjSuf) \
         tmp/modules/Cloner.$(ObjSuf) \
@@ -955 +1047 @@
         tmp/modules/PileUpJetID.$(ObjSuf) \
         tmp/modules/PileUpMerger.$(ObjSuf) \
+        tmp/modules/RecoPuFilter.$(ObjSuf) \
         tmp/modules/SimpleCalorimeter.$(ObjSuf) \
         tmp/modules/StatusPidFilter.$(ObjSuf) \
@@ -963 +1056 @@
         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/VertexFinderDA4D.$(ObjSuf) \
+        tmp/modules/VertexSorter.$(ObjSuf) \
         tmp/modules/Weighter.$(ObjSuf)
 
@@ -1569 +1666 @@
         tmp/external/tcl/tclVar.$(ObjSuf)
 
+modules/VertexFinderDA4D.h: \
+        classes/DelphesModule.h \
+        classes/DelphesClasses.h
+        @touch $@
+
+modules/TrackSmearing.h: \
+        classes/DelphesModule.h
+        @touch $@
+
 external/fastjet/ClusterSequence.hh: \
         external/fastjet/PseudoJet.hh \
@@ -1829 +1935 @@
         @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 $@
 
@@ -1902 +2018 @@
         @touch $@
 
+modules/RecoPuFilter.h: \
+        classes/DelphesModule.h
+        @touch $@
+
 modules/Hector.h: \
         classes/DelphesModule.h
@@ -2001 +2121 @@
 
 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_ATLAS.tcl

@@ -141 +141 @@
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 

cards/delphes_card_ATLAS_PileUp.tcl

@@ -177 +177 @@
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 

cards/delphes_card_CMS.tcl

@@ -15 +15 @@
 
   TrackMerger
+ 
+  ECal
+  HCal
+ 
   Calorimeter
   EFlowMerger
@@ -124 +128 @@
   set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
                                            (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e3) * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3 )               * (0.99 * exp(0.5 - pt*5.0e-4)) +
+
                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e3) * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e3)                * (0.98 * exp(0.5 - pt*5.0e-4)) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
@@ -142 +149 @@
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 
@@ -192 +199 @@
 }
 
+
+
 #############
-# Calorimeter
+#   ECAL
 #############
 
-module Calorimeter Calorimeter {
+module SimpleCalorimeter ECal {
   set ParticleInputArray ParticlePropagator/stableParticles
   set TrackInputArray TrackMerger/tracks
 
-  set TowerOutputArray towers
-  set PhotonOutputArray photons
-
+  set TowerOutputArray ecalTowers
   set EFlowTrackOutputArray eflowTracks
-  set EFlowPhotonOutputArray eflowPhotons
-  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
-
-  set ECalEnergyMin 0.5
-  set HCalEnergyMin 1.0
-
-  set ECalEnergySignificanceMin 1.0
-  set HCalEnergySignificanceMin 1.0
+  set EFlowTowerOutputArray eflowPhotons
+
+  set IsEcal true
+
+  set EnergyMin 0.5
+  set EnergySignificanceMin 2.0
+
+  set SmearTowerCenter true
+
+  set pi [expr {acos(-1)}]
+
+  # lists of the edges of each tower in eta and phi
+  # each list starts with the lower edge of the first tower
+  # the list ends with the higher edged of the last tower
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the barrel |eta| < 1.5
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 1.5 (barrel)
+  for {set i -85} {$i <= 86} {incr i} {
+    set eta [expr {$i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the endcaps 1.5 < |eta| < 3.0 (HGCAL- ECAL)
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 3
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { -2.958 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { 1.4964 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # take present CMS granularity for HF
+
+  # 0.175 x (0.175 - 0.35) resolution in eta,phi in the HF 3.0 < |eta| < 5.0
+  set PhiBins {}
+  for {set i -18} {$i <= 18} {incr i} {
+    add PhiBins [expr {$i * $pi/18.0}]
+  }
+
+  foreach eta {-5 -4.7 -4.525 -4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.958 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525 4.7 5} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+
+  add EnergyFraction {0} {0.0}
+  # energy fractions for e, gamma and pi0
+  add EnergyFraction {11} {1.0}
+  add EnergyFraction {22} {1.0}
+  add EnergyFraction {111} {1.0}
+  # energy fractions for muon, neutrinos and neutralinos
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
+  # energy fractions for K0short and Lambda
+  add EnergyFraction {310} {0.3}
+  add EnergyFraction {3122} {0.3}
+
+  # set ResolutionFormula {resolution formula as a function of eta and energy}
+
+  # for the ECAL barrel (|eta| < 1.5), see hep-ex/1306.2016 and 1502.02701
+
+  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ResolutionFormula {                      (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+
+}
+
+
+#############
+#   HCAL
+#############
+
+module SimpleCalorimeter HCal {
+  set ParticleInputArray ParticlePropagator/stableParticles
+  set TrackInputArray ECal/eflowTracks
+
+  set TowerOutputArray hcalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowNeutralHadrons
+
+  set IsEcal false
+
+  set EnergyMin 1.0
+  set EnergySignificanceMin 2.0
 
   set SmearTowerCenter true
@@ -249 +355 @@
 
   # default energy fractions {abs(PDG code)} {Fecal Fhcal}
-  add EnergyFraction {0} {0.0 1.0}
+  add EnergyFraction {0} {1.0}
   # energy fractions for e, gamma and pi0
-  add EnergyFraction {11} {1.0 0.0}
-  add EnergyFraction {22} {1.0 0.0}
-  add EnergyFraction {111} {1.0 0.0}
+  add EnergyFraction {11} {0.0}
+  add EnergyFraction {22} {0.0}
+  add EnergyFraction {111} {0.0}
   # energy fractions for muon, neutrinos and neutralinos
-  add EnergyFraction {12} {0.0 0.0}
-  add EnergyFraction {13} {0.0 0.0}
-  add EnergyFraction {14} {0.0 0.0}
-  add EnergyFraction {16} {0.0 0.0}
-  add EnergyFraction {1000022} {0.0 0.0}
-  add EnergyFraction {1000023} {0.0 0.0}
-  add EnergyFraction {1000025} {0.0 0.0}
-  add EnergyFraction {1000035} {0.0 0.0}
-  add EnergyFraction {1000045} {0.0 0.0}
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
   # energy fractions for K0short and Lambda
-  add EnergyFraction {310} {0.3 0.7}
-  add EnergyFraction {3122} {0.3 0.7}
-
-  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
-  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  add EnergyFraction {310} {0.7}
+  add EnergyFraction {3122} {0.7}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
+  set ResolutionFormula {                      (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
-}
+
+}
+
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
+}
+
+###################################################
+# Tower Merger (in case not using e-flow algorithm)
+###################################################
+
+module Merger Calorimeter {
+# add InputArray InputArray
+  add InputArray ECal/ecalTowers
+  add InputArray HCal/hcalTowers
+  set OutputArray towers
+}
+
+
 
 ####################
@@ -285 +412 @@
 module Merger EFlowMerger {
 # add InputArray InputArray
-  add InputArray Calorimeter/eflowTracks
-  add InputArray Calorimeter/eflowPhotons
-  add InputArray Calorimeter/eflowNeutralHadrons
+  add InputArray HCal/eflowTracks
+  add InputArray ECal/eflowPhotons
+  add InputArray HCal/eflowNeutralHadrons
   set OutputArray eflow
 }
@@ -296 +423 @@
 
 module Efficiency PhotonEfficiency {
-  set InputArray Calorimeter/eflowPhotons
+  set InputArray ECal/eflowPhotons
   set OutputArray photons
 
@@ -325 +452 @@
 }
 
-#################
-# Electron filter
-#################
-
-module PdgCodeFilter ElectronFilter {
-  set InputArray Calorimeter/eflowTracks
-  set OutputArray electrons
-  set Invert true
-  add PdgCode {11}
-  add PdgCode {-11}
-}
 
 #####################
@@ -382 +498 @@
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
-                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+  set EfficiencyFormula {                                     (pt <= 10.0)                * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)                * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0)                * (0.95) +
                          (abs(eta) > 2.4)                                                 * (0.00)}
 }
@@ -602 +716 @@
   add Branch Calorimeter/towers Tower Tower
 
-  add Branch Calorimeter/eflowTracks EFlowTrack Track
-  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
-  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
+  add Branch HCal/eflowTracks EFlowTrack Track
+  add Branch ECal/eflowPhotons EFlowPhoton Tower
+  add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
   add Branch GenMissingET/momentum GenMissingET MissingET
-
+ 
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron

cards/delphes_card_CMS_NoFastJet.tcl

@@ -128 +128 @@
   # resolution formula for electrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
@@ -144 +144 @@
 
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
-}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
+}
+
 
 ##############

cards/delphes_card_CMS_PileUp.tcl

@@ -178 +178 @@
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 

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 = 10000        ! 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:spareParm1 = 10000           ! max pt
+Main:spareFlag1 = on               ! true means particle gun
+Main:spareMode1 = 11               ! 1-5 - di-quark, 21 - di-gluon, 11 - single electron, 13 - single muon, 15 - single tau, 22 - single photon
+Main:spareParm1 = 10000            ! max pt
+Main:spareParm2 = 2.5              ! max eta
 
 ! 2) Settings related to output in init(), next() and stat().

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;
+ 
 }
 
@@ -219 +213 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
+  Double_t energyGuess;
   Int_t pdgCode;
 
@@ -368 +363 @@
       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 +399 @@
       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;
+        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
+        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
+        else energyGuess = momentum.E();
+
+        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fECalTowerTrackArray->Add(track);
       }
+     
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
       {
+        fHCalTrackEnergy += hcalEnergy;
         hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(hcalSigma/momentum.E() < track->TrackResolution)
-        {
-          fHCalTrackEnergy[0] += hcalEnergy;
-          fHCalTowerTrackArray[0]->Add(track);
-        }
-        else
-        {
-          fHCalTrackEnergy[1] += hcalEnergy;
-          fHCalTowerTrackArray[1]->Add(track);
-        }
+        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
+        else energyGuess = momentum.E();
+
+        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fHCalTowerTrackArray->Add(track);
       }
+      
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
@@ -476 +463 @@
   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 +476 @@
 
   if(!fTower) return;
-
 
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
@@ -554 +545 @@
     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/MomentumSmearing.cc

@@ -96 +96 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e;
+  Double_t pt, eta, phi, e, res;
 
   fItInputArray->Reset();
@@ -107 +107 @@
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
+    res = fFormula->Eval(pt, eta, phi, e);
+ 
+    // apply smearing formula
+    //pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
+    
+    res = ( res > 1.0 ) ? 1.0 : res; 
 
-    // apply smearing formula
-    pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
+    pt = LogNormal(pt, res * pt );
     
-    if(pt <= 0.0) continue;
+    //if(pt <= 0.0) continue;
 
     mother = candidate;
@@ -118 +123 @@
     phi = candidateMomentum.Phi();
     candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt*TMath::CosH(eta));
-    candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
+    //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
+    candidate->TrackResolution = res;
     candidate->AddCandidate(mother);
         
@@ -124 +130 @@
   }
 }
+//----------------------------------------------------------------
+
+Double_t MomentumSmearing::LogNormal(Double_t mean, Double_t sigma)
+{
+  Double_t a, b;
+
+  if(mean > 0.0)
+  {
+    b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
+    a = TMath::Log(mean) - 0.5*b*b;
+
+    return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
+  }
+  else
+  {
+    return 0.0;
+  }
+}
+
 
 //------------------------------------------------------------------------------

modules/MomentumSmearing.h

@@ -47 +47 @@
 private:
 
+  Double_t LogNormal(Double_t mean, Double_t sigma);
+
   DelphesFormula *fFormula; //!
 

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;
+ 
 }
 
@@ -199 +191 @@
   Double_t energy;
   Double_t sigma;
+  Double_t energyGuess;
+
   Int_t pdgCode;
 
@@ -340 +334 @@
       fTowerEnergy = 0.0;
 
-      fTrackEnergy[0] = 0.0;
-      fTrackEnergy[1] = 0.0;
+      fTrackEnergy = 0.0;
+      fTrackSigma = 0.0;
 
       fTowerTime = 0.0;
@@ -351 +345 @@
       fTowerPhotonHits = 0;
 
-      fTowerTrackArray[0]->Clear();
-      fTowerTrackArray[1]->Clear();
-    }
+      fTowerTrackArray->Clear();
+     }
 
     // check for track hits
@@ -371 +364 @@
       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;
+       sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+       if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy;
+       else energyGuess = momentum.E();
+
+       fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess);
+       fTowerTrackArray->Add(track);
       }
+       
       else
       {
@@ -403 +395 @@
     fTowerEnergy += energy;
 
-    fTowerTime += TMath::Sqrt(energy)*position.T();
-    fTowerTimeWeight += TMath::Sqrt(energy);
+    fTowerTime += energy*position.T();
+    fTowerTimeWeight += energy;
 
     fTower->AddCandidate(particle);
@@ -418 +410 @@
 {
   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 +430 @@
 
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
 
   if(fSmearTowerCenter)
@@ -464 +459 @@
   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;

readers/DelphesPythia8.cpp

@@ -153 +153 @@
 // from pythia8 example 21
 
-void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
-  Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
-
+void fillParticle(int id, double pt_max, double eta_max,
+  Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
+{
   // Reset event record to allow for new event.
   event.reset();
 
-  // Angles uniform in solid angle.
-  double cThe, sThe, phi, ee;
-  cThe = 2. * rndm.flat() - 1.;
-  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
-  phi = 2. * M_PI * rndm.flat();
-  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
-  double mm = pdt.mSel(id);
-  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
+  // Generate uniform pt and eta.
+  double pt, eta, phi, pp, ee, mm;
+  pt = pow(10, 1.0 + (log10(pt_max) - 1.0) * rndm.flat());
+  eta = (2.0 * rndm.flat() - 1.0) * eta_max;
+  phi = 2.0 * M_PI * rndm.flat();
+  pp = pt * cosh(eta);
+  mm = pdt.mSel(id);
+  ee = Pythia8::sqrtpos(pp*pp + mm*mm);
 
   // Store the particle in the event record.
-  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
-
+  event.append(id, 1, 0, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
 }
 
-void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
-  Pythia8::Rndm& rndm) {
+void fillPartons(int id, double pt_max, double eta_max,
+  Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
+{
 
   // Reset event record to allow for new event.
   event.reset();
 
-  // Angles uniform in solid angle.
-  double cThe, sThe, phi, ee;
-
   // Information on a q qbar system, to be hadronized.
 
-  cThe = 2. * rndm.flat() - 1.;
-  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
-  phi = 2. * M_PI * rndm.flat();
-  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
-  double mm = pdt.m0(type);
-  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
-  if (type == 21)
-  {
-    event.append( 21, 23, 101, 102, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee);
-    event.append( 21, 23, 102, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee);
+  // Generate uniform pt and eta.
+  double pt, eta, phi, pp, ee, mm;
+  pt = pow(10, 1.0 + (log10(pt_max) - 1.0) * rndm.flat());
+  eta = (2.0 * rndm.flat() - 1.0) * eta_max;
+  phi = 2.0 * M_PI * rndm.flat();
+  pp = pt * cosh(eta);
+  mm = pdt.m0(id);
+  ee = Pythia8::sqrtpos(pp*pp + mm*mm);
+
+  if(id == 21)
+  {
+    event.append(21, 23, 101, 102, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee);
+    event.append(21, 23, 102, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee);
   }
   else
   {
-    event.append(  type, 23, 101,   0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
-    event.append( -type, 23,   0, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee, mm);
+    event.append(id, 23, 101, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
+    event.append(-id, 23, 0, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee, mm);
   }
 }
@@ -338 +338 @@
       if (pythia->flag("Main:spareFlag1"))
       {
-        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 22) 
+        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 15 || pythia->mode("Main:spareMode1") == 22) 
         { 
-          fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
+          fillParticle(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
         }
-        else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
+        else
+        {
+          fillPartons(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
+        }
       }