Changes in / [7bb13cd:70b9632] in git

Files:

• CHANGELOG (modified) (1 diff)
• Makefile (modified) (17 diffs)
• README (modified) (1 diff)
• VERSION (modified) (1 diff)
• cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_0PU.tcl (modified) (24 diffs)
• cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_200PU.tcl (modified) (25 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.cpp (modified) (6 diffs)
• 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/PhotonConversions.cc (modified) (2 diffs)
• modules/PileUpMerger.cc (modified) (9 diffs)
• modules/RecoPuFilter.cc (added)
• modules/RecoPuFilter.h (added)
• 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) (7 diffs)
• readers/DelphesROOT.cpp (added)
• validation.sh (added)

CHANGELOG

@@ +1 @@
+3.3.3:
+- improved compatibility with ROOT >= 6.04
+- removed test of praticle stability based on ROOT pdgtable (#347, #753, #821, #855)
+- fixed bugs in DelphesCMSFWLite
+- fixed bugs in PUPPI
+- fixed compiler and linker flags for Pythia8
+- added CMS Phase II cards
+- added particle gun to DelphesPythia8
+- added jet merging to DelphesPythia8 based on Pythia8 main32.cc
+- added UseMiniCone option in Isolation module
+- added RecoPuFilter module
+- added back OldCalorimeter for CheckMate bwd compatibilty (#591)
+- updated tracking and calorimeter resolution in CMS cards
+
 3.3.2:
 - added pre-validated card for CMS upgrade Phase II studies

Makefile

@@ -250 +250 @@
         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
 DelphesSTDHEP$(ExeSuf): \
         tmp/readers/DelphesSTDHEP.$(ObjSuf)
@@ -265 +278 @@
         DelphesHepMC$(ExeSuf) \
         DelphesLHEF$(ExeSuf) \
+        DelphesROOT$(ExeSuf) \
         DelphesSTDHEP$(ExeSuf)
 
@@ -270 +284 @@
         tmp/readers/DelphesHepMC.$(ObjSuf) \
         tmp/readers/DelphesLHEF.$(ObjSuf) \
+        tmp/readers/DelphesROOT.$(ObjSuf) \
         tmp/readers/DelphesSTDHEP.$(ObjSuf)
 
@@ -384 +399 @@
         modules/EnergySmearing.h \
         modules/MomentumSmearing.h \
+        modules/TrackSmearing.h \
         modules/ImpactParameterSmearing.h \
         modules/TimeSmearing.h \
@@ -407 +423 @@
         modules/StatusPidFilter.h \
         modules/PdgCodeFilter.h \
+        modules/BeamSpotFilter.h \
+        modules/RecoPuFilter.h \
         modules/Cloner.h \
         modules/Weighter.h \
@@ -412 +430 @@
         modules/JetFlavorAssociation.h \
         modules/JetFakeParticle.h \
+        modules/VertexSorter.h \
+        modules/VertexFinder.h \
+        modules/VertexFinderDA4D.h \
         modules/ExampleModule.h
 tmp/modules/ModulesDict$(PcmSuf): \
@@ -618 +639 @@
         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) \
@@ -850 +880 @@
         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) \
@@ -911 +950 @@
         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 \
@@ -933 +981 @@
         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 \
@@ -996 +1074 @@
         tmp/modules/AngularSmearing.$(ObjSuf) \
         tmp/modules/BTagging.$(ObjSuf) \
+        tmp/modules/BeamSpotFilter.$(ObjSuf) \
         tmp/modules/Calorimeter.$(ObjSuf) \
         tmp/modules/Cloner.$(ObjSuf) \
@@ -1020 +1099 @@
         tmp/modules/PileUpJetID.$(ObjSuf) \
         tmp/modules/PileUpMerger.$(ObjSuf) \
+        tmp/modules/RecoPuFilter.$(ObjSuf) \
         tmp/modules/SimpleCalorimeter.$(ObjSuf) \
         tmp/modules/StatusPidFilter.$(ObjSuf) \
@@ -1028 +1108 @@
         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)
 
@@ -1634 +1718 @@
         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 \
@@ -1894 +1987 @@
         @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 $@
 
@@ -1967 +2070 @@
         @touch $@
 
+modules/RecoPuFilter.h: \
+        classes/DelphesModule.h
+        @touch $@
+
 modules/Hector.h: \
         classes/DelphesModule.h
@@ -2066 +2173 @@
 
 modules/FastJetFinder.h: \
+        classes/DelphesModule.h
+        @touch $@
+
+modules/BeamSpotFilter.h: \
         classes/DelphesModule.h
         @touch $@

README

@@ -4 +4 @@
 Commands to get the code:
 
-   wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.3.2.tar.gz
+   wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.3.3.tar.gz
 
-   tar -zxf Delphes-3.3.2.tar.gz
+   tar -zxf Delphes-3.3.3.tar.gz
 
 Commands to compile the code:
 
-   cd Delphes-3.3.2
+   cd Delphes-3.3.3
 
    make

VERSION

@@ -1 @@
-3.3.2
+3.3.3

cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_0PU.tcl

@@ -33 +33 @@
   PhotonEnergySmearing
   ElectronFilter
+
   TrackPileUpSubtractor
+  RecoPuFilter
 
   TowerMerger
   NeutralEFlowMerger
-  EFlowMergerAllTracks
+
   EFlowMerger
+  EFlowMergerCHS
+  Rho
 
   LeptonFilterNoLep
@@ -48 +52 @@
 
   PhotonIsolation
+  PhotonIsolationCHS
   PhotonEfficiency
+  PhotonEfficiencyCHS
 
   ElectronIsolation
+  ElectronIsolationCHS
+
   ElectronEfficiency
+  ElectronEfficiencyCHS
 
   MuonIsolation
+  MuonIsolationCHS
+
   MuonLooseIdEfficiency
   MuonTightIdEfficiency
+
+  MuonLooseIdEfficiencyCHS
+  MuonTightIdEfficiencyCHS
 
   NeutrinoFilter
@@ -68 +82 @@
   FastJetFinder
   FastJetFinderAK8
+  JetPileUpSubtractor
+  JetPileUpSubtractorAK8
   FastJetFinderPUPPI
   FastJetFinderPUPPIAK8
@@ -121 +137 @@
 
   # maximum spread in the beam direction in m
-  set ZVertexSpread 0
+  set ZVertexSpread 0.25
 
   # maximum spread in time in s
-  set TVertexSpread 0
+  set TVertexSpread 800E-12
 
   # vertex smearing formula f(z,t) (z,t need to be respectively given in m,s) - {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
-  set VertexDistributionFormula 1
+  set VertexDistributionFormula {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
 
 }
@@ -273 +289 @@
   source muonMomentumResolution.tcl
 }
-
-
-
 
 ##############
@@ -527 +540 @@
 }
 
+########################
+# Reco PU filter
+########################
+
+module RecoPuFilter RecoPuFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray eflowTracks
+}
 
 ###################################################
@@ -539 +560 @@
 }
 
-
 ####################
 # Neutral eflow erger
@@ -551 +571 @@
 }
 
-
-####################
+#####################
 # Energy flow merger
-####################
+#####################
 
 module Merger EFlowMerger {
@@ -564 +583 @@
 }
 
-##################################
-# Energy flow merger (all tracks)
-##################################
-
-module Merger EFlowMergerAllTracks {
+############################
+# Energy flow merger no PU
+############################
+
+module Merger EFlowMergerCHS {
 # add InputArray InputArray
-  add InputArray TrackMerger/tracks
+  add InputArray RecoPuFilter/eflowTracks
   add InputArray PhotonEnergySmearing/eflowPhotons
   add InputArray HCal/eflowNeutralHadrons
@@ -695 +714 @@
 }
 
-
 #####################
 # MC truth jet finder
@@ -735 +753 @@
 }
 
-############
+
+#############
+# Rho pile-up
+#############
+
+module FastJetFinder Rho {
+#  set InputArray Calorimeter/towers
+  set InputArray EFlowMergerCHS/eflow
+
+  set ComputeRho true
+  set RhoOutputArray rho
+
+  # area algorithm: 0 Do not compute area, 1 Active area explicit ghosts, 2 One ghost passive area, 3 Passive area, 4 Voronoi, 5 Active area
+  set AreaAlgorithm 5
+
+  # jet algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 4
+  set ParameterR 0.4
+  set GhostEtaMax 5.0
+
+  add RhoEtaRange -5.0 -4.0
+  add RhoEtaRange -4.0 -1.5
+  add RhoEtaRange -1.5 1.5
+  add RhoEtaRange 1.5 4.0
+  add RhoEtaRange 4.0 5.0
+
+  set JetPTMin 0.0
+}
+
+
+##############
 # Jet finder
-############
+##############
 
 module FastJetFinder FastJetFinder {
 #  set InputArray TowerMerger/towers
-  set InputArray EFlowMerger/eflow
+  set InputArray EFlowMergerCHS/eflow
 
   set OutputArray jets
+
+  set AreaAlgorithm 5
 
   # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
@@ -755 +805 @@
 module FastJetFinder FastJetFinderAK8 {
 #  set InputArray TowerMerger/towers
-  set InputArray EFlowMerger/eflow
+  set InputArray EFlowMergerCHS/eflow
 
   set OutputArray jets
+
+  set AreaAlgorithm 5
 
   # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
@@ -784 +836 @@
 }
 
+###########################
+# Jet Pile-Up Subtraction
+###########################
+
+module JetPileUpSubtractor JetPileUpSubtractor {
+  set JetInputArray FastJetFinder/jets
+  set RhoInputArray Rho/rho
+
+  set OutputArray jets
+
+  set JetPTMin 15.0
+}
+
+##############################
+# Jet Pile-Up Subtraction AK8
+##############################
+
+module JetPileUpSubtractor JetPileUpSubtractorAK8 {
+  set JetInputArray FastJetFinderAK8/jets
+  set RhoInputArray Rho/rho
+
+  set OutputArray jets
+
+  set JetPTMin 15.0
+}
+
 module FastJetFinder FastJetFinderPUPPI {
 #  set InputArray TowerMerger/towers
@@ -833 +911 @@
 
 module EnergyScale JetEnergyScale {
-  set InputArray FastJetFinder/jets
+  set InputArray JetPileUpSubtractor/jets
   set OutputArray jets
 
@@ -841 +919 @@
 
 module EnergyScale JetEnergyScaleAK8 {
-  set InputArray FastJetFinderAK8/jets
+  set InputArray JetPileUpSubtractorAK8/jets
   set OutputArray jets
 
@@ -907 +985 @@
 }
 
+
+########################
+# Photon isolation CHS #
+########################
+
+module Isolation PhotonIsolationCHS {
+
+  # particle for which calculate the isolation
+  set CandidateInputArray PhotonFilter/photons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  # output array
+  set OutputArray photons
+
+  # isolation cone
+  set DeltaRMax 0.3
+
+  # minimum pT
+  set PTMin     1.0
+
+  # iso ratio to cut
+  set PTRatioMax 9999.
+
+}
 
 
@@ -929 +1033 @@
 
 
+#####################
+# Photon efficiency #
+#####################
+
+module Efficiency PhotonEfficiencyCHS {
+
+  ## input particles
+  set InputArray PhotonIsolationCHS/photons
+  ## output particles
+  set OutputArray photons
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+  # efficiency formula for photons
+  set EfficiencyFormula {                      (pt <= 10.0) * (0.00) + \
+                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.9635) + \
+         (abs(eta) > 1.5 && abs(eta) <= 4.0) * (pt > 10.0)  * (0.9624) + \
+         (abs(eta) > 4.0)                                   * (0.00)}
+
+}
+
 ######################
 # Electron isolation #
@@ -949 +1072 @@
 }
 
+
+##########################
+# Electron isolation CHS #
+##########################
+
+module Isolation ElectronIsolationCHS {
+
+  set CandidateInputArray ElectronFilter/electrons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray electrons
+
+  set DeltaRMax 0.3
+  set PTMin 1.0
+  set PTRatioMax 9999.
+
+}
 
 
@@ -995 +1137 @@
 }
 
+###########################
+# Electron efficiency CHS #
+###########################
+
+module Efficiency ElectronEfficiencyCHS {
+
+  set InputArray ElectronIsolationCHS/electrons
+  set OutputArray electrons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+  # efficiency formula for electrons
+  set EfficiencyFormula {
+                                      (pt <= 4.0)  * (0.00) + \
+                         (abs(eta) <= 1.45 ) * (pt >  4.0 && pt <= 6.0)   * (0.50) + \
+                         (abs(eta) <= 1.45 ) * (pt >  6.0 && pt <= 8.0)   * (0.70) + \
+                         (abs(eta) <= 1.45 ) * (pt >  8.0 && pt <= 10.0)  * (0.85) + \
+                         (abs(eta) <= 1.45 ) * (pt > 10.0 && pt <= 30.0)  * (0.94) + \
+                         (abs(eta) <= 1.45 ) * (pt > 30.0 && pt <= 50.0)  * (0.97) + \
+                         (abs(eta) <= 1.45 ) * (pt > 50.0 && pt <= 70.0)  * (0.98) + \
+                         (abs(eta) <= 1.45 ) * (pt > 70.0 )  * (1.0) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt >  4.0 && pt <= 10.0)   * (0.35) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 10.0 && pt <= 30.0)   * (0.40) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 30.0 && pt <= 70.0)   * (0.45) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 70.0 )  * (0.55) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt >  4.0 && pt <= 10.0)  * (0.75) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 10.0 && pt <= 30.0)  * (0.85) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 30.0 && pt <= 50.0)  * (0.95) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 50.0 && pt <= 70.0)  * (0.95) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 70.0 )  * (1.0) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt >  4.0 && pt <= 10.0)  * (0.65) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 10.0 && pt <= 30.0)  * (0.75) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 30.0 && pt <= 50.0)  * (0.90) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 50.0 && pt <= 70.0)  * (0.90) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 70.0 )  * (0.90) + \
+                         (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 4.0 && pt <= 10.0) * (0.65) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 10.0 && pt <= 30.0) * (0.75) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 30.0 && pt <= 50.0) * (0.90) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 50.0 && pt <= 70.0) * (0.90) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 70.0 ) * (0.90) + \
+                                          (abs(eta) > 4.0) * (0.00)
+
+  }
+}
+
+
 ##################
 # Muon isolation #
@@ -1013 +1200 @@
 }
 
-
-##################
-# Muon Loose Id  #
-##################
+######################
+# Muon isolation CHS #
+######################
+
+module Isolation MuonIsolationCHS {
+  set CandidateInputArray MuonMomentumSmearing/muons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray muons
+
+  set DeltaRMax 0.3
+  set PTMin 1.0
+  set PTRatioMax 9999.
+
+}
+
+
+#####################
+# Muon Loose Id     #
+#####################
 
 module Efficiency MuonLooseIdEfficiency {
@@ -1038 +1243 @@
 }
 
+
+#####################
+# Muon Loose Id CHS #
+#####################
+
+module Efficiency MuonLooseIdEfficiencyCHS {
+    set InputArray MuonIsolationCHS/muons
+    set OutputArray muons
+    # tracking + LooseID efficiency formula for muons
+    source muonLooseId.tcl
+
+}
+
+
+######################
+# Muon Tight Id  CHS #
+######################
+
+module Efficiency MuonTightIdEfficiencyCHS {
+    set InputArray MuonIsolationCHS/muons
+    set OutputArray muons
+    # tracking + TightID efficiency formula for muons
+    source muonTightId.tcl
+}
 
 
@@ -1250 +1479 @@
 module StatusPidFilter GenParticleFilter {
 
-    set InputArray  Delphes/allParticles
+    set InputArray Delphes/allParticles
     set OutputArray filteredParticles
     set PTMin 5.0
@@ -1278 +1507 @@
   add Branch MuonLooseIdEfficiency/muons MuonLoose Muon
   add Branch MuonTightIdEfficiency/muons MuonTight Muon
+
+  add Branch PhotonEfficiencyCHS/photons PhotonCHS Photon
+  add Branch ElectronEfficiencyCHS/electrons ElectronCHS Electron
+  add Branch MuonLooseIdEfficiencyCHS/muons MuonLooseCHS Muon
+  add Branch MuonTightIdEfficiencyCHS/muons MuonTightCHS Muon
 
   add Branch JetEnergyScale/jets Jet Jet
@@ -1289 +1523 @@
   add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
-}
+
+}

cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_200PU.tcl

@@ -33 +33 @@
   PhotonEnergySmearing
   ElectronFilter
+
   TrackPileUpSubtractor
+  RecoPuFilter
 
   TowerMerger
   NeutralEFlowMerger
-  EFlowMergerAllTracks
+
   EFlowMerger
+  EFlowMergerCHS
+  Rho
 
   LeptonFilterNoLep
@@ -48 +52 @@
 
   PhotonIsolation
+  PhotonIsolationCHS
   PhotonEfficiency
+  PhotonEfficiencyCHS
 
   ElectronIsolation
+  ElectronIsolationCHS
+
   ElectronEfficiency
+  ElectronEfficiencyCHS
 
   MuonIsolation
+  MuonIsolationCHS
+
   MuonLooseIdEfficiency
   MuonTightIdEfficiency
+
+  MuonLooseIdEfficiencyCHS
+  MuonTightIdEfficiencyCHS
 
   NeutrinoFilter
@@ -68 +82 @@
   FastJetFinder
   FastJetFinderAK8
+  JetPileUpSubtractor
+  JetPileUpSubtractorAK8
   FastJetFinderPUPPI
   FastJetFinderPUPPIAK8
@@ -126 +142 @@
   set TVertexSpread 800E-12
 
-  # vertex smearing formula f(z,t) (z,t need to be respectively given in m,s)
+  # vertex smearing formula f(z,t) (z,t need to be respectively given in m,s) - {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
   set VertexDistributionFormula {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
 
@@ -138 +154 @@
 module ParticlePropagator ParticlePropagator {
   set InputArray PileUpMerger/stableParticles
+  #set InputArray Delphes/stableParticles
 
   set OutputArray stableParticles
@@ -272 +289 @@
   source muonMomentumResolution.tcl
 }
-
-
-
 
 ##############
@@ -526 +540 @@
 }
 
+########################
+# Reco PU filter
+########################
+
+module RecoPuFilter RecoPuFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray eflowTracks
+}
 
 ###################################################
@@ -538 +560 @@
 }
 
-
 ####################
 # Neutral eflow erger
@@ -550 +571 @@
 }
 
-
-####################
+#####################
 # Energy flow merger
-####################
+#####################
 
 module Merger EFlowMerger {
@@ -563 +583 @@
 }
 
-##################################
-# Energy flow merger (all tracks)
-##################################
-
-module Merger EFlowMergerAllTracks {
+############################
+# Energy flow merger no PU
+############################
+
+module Merger EFlowMergerCHS {
 # add InputArray InputArray
-  add InputArray TrackMerger/tracks
+  add InputArray RecoPuFilter/eflowTracks
   add InputArray PhotonEnergySmearing/eflowPhotons
   add InputArray HCal/eflowNeutralHadrons
@@ -694 +714 @@
 }
 
-
 #####################
 # MC truth jet finder
@@ -735 +754 @@
 
 
-
-############
+#############
+# Rho pile-up
+#############
+
+module FastJetFinder Rho {
+#  set InputArray Calorimeter/towers
+  set InputArray EFlowMergerCHS/eflow
+
+  set ComputeRho true
+  set RhoOutputArray rho
+
+  # area algorithm: 0 Do not compute area, 1 Active area explicit ghosts, 2 One ghost passive area, 3 Passive area, 4 Voronoi, 5 Active area
+  set AreaAlgorithm 5
+
+  # jet algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 4
+  set ParameterR 0.4
+  set GhostEtaMax 5.0
+
+  add RhoEtaRange -5.0 -4.0
+  add RhoEtaRange -4.0 -1.5
+  add RhoEtaRange -1.5 1.5
+  add RhoEtaRange 1.5 4.0
+  add RhoEtaRange 4.0 5.0
+
+  set JetPTMin 0.0
+}
+
+
+##############
 # Jet finder
-############
+##############
 
 module FastJetFinder FastJetFinder {
 #  set InputArray TowerMerger/towers
-  set InputArray EFlowMerger/eflow
+  set InputArray EFlowMergerCHS/eflow
 
   set OutputArray jets
+
+  set AreaAlgorithm 5
 
   # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
@@ -756 +805 @@
 module FastJetFinder FastJetFinderAK8 {
 #  set InputArray TowerMerger/towers
-  set InputArray EFlowMerger/eflow
+  set InputArray EFlowMergerCHS/eflow
 
   set OutputArray jets
+
+  set AreaAlgorithm 5
 
   # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
@@ -785 +836 @@
 }
 
+###########################
+# Jet Pile-Up Subtraction
+###########################
+
+module JetPileUpSubtractor JetPileUpSubtractor {
+  set JetInputArray FastJetFinder/jets
+  set RhoInputArray Rho/rho
+
+  set OutputArray jets
+
+  set JetPTMin 15.0
+}
+
+##############################
+# Jet Pile-Up Subtraction AK8
+##############################
+
+module JetPileUpSubtractor JetPileUpSubtractorAK8 {
+  set JetInputArray FastJetFinderAK8/jets
+  set RhoInputArray Rho/rho
+
+  set OutputArray jets
+
+  set JetPTMin 15.0
+}
+
 module FastJetFinder FastJetFinderPUPPI {
 #  set InputArray TowerMerger/towers
@@ -834 +911 @@
 
 module EnergyScale JetEnergyScale {
-  set InputArray FastJetFinder/jets
+  set InputArray JetPileUpSubtractor/jets
   set OutputArray jets
 
@@ -842 +919 @@
 
 module EnergyScale JetEnergyScaleAK8 {
-  set InputArray FastJetFinderAK8/jets
+  set InputArray JetPileUpSubtractorAK8/jets
   set OutputArray jets
 
@@ -908 +985 @@
 }
 
+
+########################
+# Photon isolation CHS #
+########################
+
+module Isolation PhotonIsolationCHS {
+
+  # particle for which calculate the isolation
+  set CandidateInputArray PhotonFilter/photons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  # output array
+  set OutputArray photons
+
+  # isolation cone
+  set DeltaRMax 0.3
+
+  # minimum pT
+  set PTMin     1.0
+
+  # iso ratio to cut
+  set PTRatioMax 9999.
+
+}
 
 
@@ -930 +1033 @@
 
 
+#####################
+# Photon efficiency #
+#####################
+
+module Efficiency PhotonEfficiencyCHS {
+
+  ## input particles
+  set InputArray PhotonIsolationCHS/photons
+  ## output particles
+  set OutputArray photons
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+  # efficiency formula for photons
+  set EfficiencyFormula {                      (pt <= 10.0) * (0.00) + \
+                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.9635) + \
+         (abs(eta) > 1.5 && abs(eta) <= 4.0) * (pt > 10.0)  * (0.9624) + \
+         (abs(eta) > 4.0)                                   * (0.00)}
+
+}
+
 ######################
 # Electron isolation #
@@ -950 +1072 @@
 }
 
+
+##########################
+# Electron isolation CHS #
+##########################
+
+module Isolation ElectronIsolationCHS {
+
+  set CandidateInputArray ElectronFilter/electrons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray electrons
+
+  set DeltaRMax 0.3
+  set PTMin 1.0
+  set PTRatioMax 9999.
+
+}
 
 
@@ -996 +1137 @@
 }
 
+###########################
+# Electron efficiency CHS #
+###########################
+
+module Efficiency ElectronEfficiencyCHS {
+
+  set InputArray ElectronIsolationCHS/electrons
+  set OutputArray electrons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+  # efficiency formula for electrons
+  set EfficiencyFormula {
+                                      (pt <= 4.0)  * (0.00) + \
+                         (abs(eta) <= 1.45 ) * (pt >  4.0 && pt <= 6.0)   * (0.50) + \
+                         (abs(eta) <= 1.45 ) * (pt >  6.0 && pt <= 8.0)   * (0.70) + \
+                         (abs(eta) <= 1.45 ) * (pt >  8.0 && pt <= 10.0)  * (0.85) + \
+                         (abs(eta) <= 1.45 ) * (pt > 10.0 && pt <= 30.0)  * (0.94) + \
+                         (abs(eta) <= 1.45 ) * (pt > 30.0 && pt <= 50.0)  * (0.97) + \
+                         (abs(eta) <= 1.45 ) * (pt > 50.0 && pt <= 70.0)  * (0.98) + \
+                         (abs(eta) <= 1.45 ) * (pt > 70.0 )  * (1.0) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt >  4.0 && pt <= 10.0)   * (0.35) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 10.0 && pt <= 30.0)   * (0.40) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 30.0 && pt <= 70.0)   * (0.45) + \
+                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 70.0 )  * (0.55) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt >  4.0 && pt <= 10.0)  * (0.75) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 10.0 && pt <= 30.0)  * (0.85) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 30.0 && pt <= 50.0)  * (0.95) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 50.0 && pt <= 70.0)  * (0.95) + \
+                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 70.0 )  * (1.0) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt >  4.0 && pt <= 10.0)  * (0.65) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 10.0 && pt <= 30.0)  * (0.75) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 30.0 && pt <= 50.0)  * (0.90) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 50.0 && pt <= 70.0)  * (0.90) + \
+                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 70.0 )  * (0.90) + \
+                         (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 4.0 && pt <= 10.0) * (0.65) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 10.0 && pt <= 30.0) * (0.75) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 30.0 && pt <= 50.0) * (0.90) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 50.0 && pt <= 70.0) * (0.90) + \
+                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 70.0 ) * (0.90) + \
+                                          (abs(eta) > 4.0) * (0.00)
+
+  }
+}
+
+
 ##################
 # Muon isolation #
@@ -1014 +1200 @@
 }
 
-
-##################
-# Muon Loose Id  #
-##################
+######################
+# Muon isolation CHS #
+######################
+
+module Isolation MuonIsolationCHS {
+  set CandidateInputArray MuonMomentumSmearing/muons
+
+  # isolation collection
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray muons
+
+  set DeltaRMax 0.3
+  set PTMin 1.0
+  set PTRatioMax 9999.
+
+}
+
+
+#####################
+# Muon Loose Id     #
+#####################
 
 module Efficiency MuonLooseIdEfficiency {
@@ -1039 +1243 @@
 }
 
+
+#####################
+# Muon Loose Id CHS #
+#####################
+
+module Efficiency MuonLooseIdEfficiencyCHS {
+    set InputArray MuonIsolationCHS/muons
+    set OutputArray muons
+    # tracking + LooseID efficiency formula for muons
+    source muonLooseId.tcl
+
+}
+
+
+######################
+# Muon Tight Id  CHS #
+######################
+
+module Efficiency MuonTightIdEfficiencyCHS {
+    set InputArray MuonIsolationCHS/muons
+    set OutputArray muons
+    # tracking + TightID efficiency formula for muons
+    source muonTightId.tcl
+}
 
 
@@ -1251 +1479 @@
 module StatusPidFilter GenParticleFilter {
 
-    set InputArray  Delphes/allParticles
+    set InputArray Delphes/allParticles
     set OutputArray filteredParticles
     set PTMin 5.0
@@ -1279 +1507 @@
   add Branch MuonLooseIdEfficiency/muons MuonLoose Muon
   add Branch MuonTightIdEfficiency/muons MuonTight Muon
+
+  add Branch PhotonEfficiencyCHS/photons PhotonCHS Photon
+  add Branch ElectronEfficiencyCHS/electrons ElectronCHS Electron
+  add Branch MuonLooseIdEfficiencyCHS/muons MuonLooseCHS Muon
+  add Branch MuonTightIdEfficiencyCHS/muons MuonTightCHS Muon
 
   add Branch JetEnergyScale/jets Jet Jet
@@ -1290 +1523 @@
   add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
-}
+
+}

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

examples/Validation.cpp

@@ -21 +21 @@
 #include <utility>
 #include <vector>
+#include <typeinfo>
 
 #include "TROOT.h"
@@ -28 +29 @@
 #include "TString.h"
 
+#include "TH1.h"
 #include "TH2.h"
+#include "TMath.h"
+#include "TStyle.h"
+#include "TGraph.h"
+#include "TCanvas.h"
 #include "THStack.h"
 #include "TLegend.h"
@@ -34 +40 @@
 #include "TClonesArray.h"
 #include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TMultiGraph.h"
 
 #include "classes/DelphesClasses.h"
@@ -47 +55 @@
 //------------------------------------------------------------------------------
 
-// Here you can put your analysis macro
-
-#include "Validation.C"
+double ptrangemin = 10;
+double ptrangemax = 10000;
+static const int Nbins = 20;
+
+int objStyle = 1;
+int trackStyle = 7;
+int towerStyle = 3;
+
+Color_t objColor = kBlack;
+Color_t trackColor = kBlack;
+Color_t towerColor = kBlack;
+
+double effLegXmin = 0.22;
+double effLegXmax = 0.7;
+double effLegYmin = 0.22;
+double effLegYmax = 0.5;
+
+double resLegXmin = 0.62;
+double resLegXmax = 0.9;
+double resLegYmin = 0.52;
+double resLegYmax = 0.85;
+
+double topLeftLegXmin = 0.22;
+double topLeftLegXmax = 0.7;
+double topLeftLegYmin = 0.52;
+double topLeftLegYmax = 0.85;
+
+
+struct resolPlot
+{
+  TH1 *cenResolHist;
+  TH1 *fwdResolHist;
+  double ptmin;
+  double ptmax;
+  double xmin;
+  double xmax;
+  TString obj;
+
+  resolPlot();
+  resolPlot(double ptdown, double ptup, TString object);
+  void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
+  void print(){std::cout << ptmin << std::endl;}
+};
+
+
+resolPlot::resolPlot()
+{
+}
+
+resolPlot::resolPlot(double ptdown, double ptup, TString object)
+{
+  this->set(ptdown,ptup,object);
+}
+
+void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
+{
+  ptmin = ptdown;
+  ptmax = ptup;
+  obj = object;
+
+  cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 200,  xmin, xmax);
+  fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 200,  xmin, xmax);
+}
+
+void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
+{
+  double width;
+  double ptdown;
+  double ptup;
+  resolPlot ptemp;
+
+  for (int i = 0; i < Nbins; i++)
+  {
+    width = (ptmax - ptmin) / Nbins;
+    ptdown = TMath::Power(10,ptmin + i * width );
+    ptup = TMath::Power(10,ptmin + (i+1) * width );
+    ptemp.set(ptdown, ptup, obj, xmin, xmax);
+    histos->push_back(ptemp);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+class ExRootResult;
+class ExRootTreeReader;
+
+//------------------------------------------------------------------------------
+
+void BinLogX(TH1*h)
+{
+  TAxis *axis = h->GetXaxis();
+  int bins = axis->GetNbins();
+
+  Axis_t from = axis->GetXmin();
+  Axis_t to = axis->GetXmax();
+  Axis_t width = (to - from) / bins;
+  Axis_t *new_bins = new Axis_t[bins + 1];
+
+  for (int i = 0; i <= bins; i++)
+  {
+    new_bins[i] = TMath::Power(10, from + i * width);
+  }
+  axis->Set(bins, new_bins);
+  delete new_bins;
+}
+
+
+//------------------------------------------------------------------------------
+
+template<typename T>
+std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
+{
+
+  cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
+
+  Long64_t allEntries = treeReader->GetEntries();
+
+  GenParticle *particle;
+  T *recoObj;
+
+  TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
+
+  Float_t deltaR;
+  Float_t pt, eta;
+  Long64_t entry;
+
+  Int_t i, j;
+
+  TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
+  TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
+  TH1D *histGenPtfwd  = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
+  TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
+
+
+  BinLogX(histGenPtcen);
+  BinLogX(histRecoPtcen);
+  BinLogX(histGenPtfwd);
+  BinLogX(histRecoPtfwd);
+
+  // Loop over all events
+  for(entry = 0; entry < allEntries; ++entry)
+  {
+    // Load selected branches with data from specified event
+    treeReader->ReadEntry(entry);
+
+    // Loop over all generated particle in event
+    for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
+    {
+
+      particle = (GenParticle*) branchParticle->At(i);
+      genMomentum = particle->P4();
+
+      deltaR = 999;
+
+      if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
+      {
+
+        // Loop over all reco object in event
+        for(j = 0; j < branchReco->GetEntriesFast(); ++j)
+        {
+          recoObj = (T*)branchReco->At(j);
+          recoMomentum = recoObj->P4();
+          // this is simply to avoid warnings from initial state particle
+          // having infite rapidity ...
+        //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
+
+          // take the closest parton candidate
+          if(TMath::Abs(pdgID) == 5)
+          {
+            Jet *jet = (Jet *)recoObj;
+            if(jet->BTag != 1) continue;
+          }
+          if(TMath::Abs(pdgID) == 15)
+          {
+            Jet *jet = (Jet *)recoObj;
+            if(jet->TauTag != 1) continue;
+          }
+          if(genMomentum.DeltaR(recoMomentum) < deltaR)
+          {
+            deltaR = genMomentum.DeltaR(recoMomentum);
+            bestRecoMomentum = recoMomentum;
+          }
+        }
+
+        pt  = genMomentum.Pt();
+        eta = genMomentum.Eta();
+
+        if (TMath::Abs(eta) < 1.5)
+        {
+          histGenPtcen->Fill(pt);
+          if(deltaR < 0.3) { histRecoPtcen->Fill(pt); }
+        }
+        else if (TMath::Abs(eta) < 2.5)
+        {
+          histGenPtfwd->Fill(pt);
+          if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); }
+
+        }
+      }
+    }
+  }
+
+
+  std::pair<TH1D*,TH1D*> histos;
+
+  histRecoPtcen->Divide(histGenPtcen);
+  histRecoPtfwd->Divide(histGenPtfwd);
+
+  histos.first = histRecoPtcen;
+  histos.second = histRecoPtfwd;
+
+  return histos;
+}
+
+template<typename T>
+void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
+{
+  Long64_t allEntries = treeReader->GetEntries();
+
+  cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
+
+  GenParticle *particle;
+  T* recoObj;
+
+  TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
+
+  Float_t deltaR;
+  Float_t pt, eta;
+  Long64_t entry;
+
+  Int_t i, j, bin;
+
+  // Loop over all events
+  for(entry = 0; entry < allEntries; ++entry)
+  {
+    // Load selected branches with data from specified event
+    treeReader->ReadEntry(entry);
+
+    // Loop over all reconstructed jets in event
+    for(i = 0; i < branchReco->GetEntriesFast(); ++i)
+    {
+      recoObj = (T*) branchReco->At(i);
+      recoMomentum = recoObj->P4();
+
+      deltaR = 999;
+
+     // Loop over all hard partons in event
+     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
+     {
+        particle = (GenParticle*) branchParticle->At(j);
+        if (particle->PID == pdgID && particle->Status == 1)
+        {
+          genMomentum = particle->P4();
+
+          // this is simply to avoid warnings from initial state particle
+          // having infite rapidity ...
+          if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
+
+          // take the closest parton candidate
+          if(genMomentum.DeltaR(recoMomentum) < deltaR)
+          {
+             deltaR = genMomentum.DeltaR(recoMomentum);
+             bestGenMomentum = genMomentum;
+          }
+        }
+      }
+
+      if(deltaR < 0.3)
+      {
+        pt  = bestGenMomentum.Pt();
+        eta = TMath::Abs(bestGenMomentum.Eta());
+
+        for (bin = 0; bin < Nbins; bin++)
+        {
+          if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
+          {
+            if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
+            else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
+          }
+        }
+      }
+    }
+  }
+}
+void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
+{
+
+  Long64_t allEntries = treeReader->GetEntries();
+
+  cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
+
+  Jet *jet, *genjet;
+
+  TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
+
+  Float_t deltaR;
+  Float_t pt, eta;
+  Long64_t entry;
+
+  Int_t i, j, bin;
+
+  // Loop over all events
+  for(entry = 0; entry < allEntries; ++entry)
+  {
+    // Load selected branches with data from specified event
+    treeReader->ReadEntry(entry);
+
+    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
+
+    // Loop over all reconstructed jets in event
+    for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
+    {
+
+      jet = (Jet*) branchJet->At(i);
+      jetMomentum = jet->P4();
+
+      deltaR = 999;
+
+     // Loop over all hard partons in event
+     for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
+     {
+        genjet = (Jet*) branchGenJet->At(j);
+
+        genJetMomentum = genjet->P4();
+
+        // this is simply to avoid warnings from initial state particle
+        // having infite rapidity ...
+        if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
+
+        // take the closest parton candidate
+        if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
+        {
+           deltaR = genJetMomentum.DeltaR(jetMomentum);
+           bestGenJetMomentum = genJetMomentum;
+        }
+      }
+
+      if(deltaR < 0.25)
+      {
+        pt  = genJetMomentum.Pt();
+        eta = TMath::Abs(genJetMomentum.Eta());
+
+        for (bin = 0; bin < Nbins; bin++)
+        {
+            if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
+            {
+                histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
+            }
+            else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
+            {
+                histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
+            }
+
+        }
+      }
+    }
+  }
+}
+
+void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
+{
+
+  Long64_t allEntries = treeReader->GetEntries();
+
+  cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
+
+  MissingET *met;
+  ScalarHT *scalarHT;
+
+  Long64_t entry;
+
+  Int_t bin;
+  Double_t ht;
+
+  Jet *jet;
+  TLorentzVector p1, p2;
+
+  // Loop over all events
+  for(entry = 0; entry < allEntries; ++entry)
+  {
+    // Load selected branches with data from specified event
+    treeReader->ReadEntry(entry);
+
+    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
+
+    if (branchJet->GetEntriesFast() > 1)
+    {
+
+      jet = (Jet*) branchJet->At(0);
+      p1 = jet->P4();
+      jet = (Jet*) branchJet->At(1);
+      p2 = jet->P4();
+
+      met = (MissingET*) branchMet->At(0);
+      scalarHT = (ScalarHT*) branchScalarHT->At(0);
+      ht = scalarHT->HT;
+
+      if(p1.Pt() < 0.75*ht/2) continue;
+      if(p2.Pt() < 0.75*ht/2) continue;
+
+      for (bin = 0; bin < Nbins; bin++)
+      {
+        if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
+        {
+          histos->at(bin).cenResolHist->Fill(met->P4().Px());
+          histos->at(bin).fwdResolHist->Fill(met->P4().Py());
+        }
+      }
+    }
+  }
+}
+
+
+std::pair<Double_t, Double_t> GausFit(TH1* hist)
+{
+  TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
+  hist->Fit("f1","RQ");
+  Double_t sig = f1->GetParameter(2);
+  Double_t sigErr = f1->GetParError(2);
+  delete f1;
+  return make_pair (sig, sigErr);
+}
+
+
+TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
+{
+  Int_t bin;
+  Int_t count = 0;
+  TGraphErrors gr = TGraphErrors(Nbins/2);
+  Double_t sig = 0;
+  Double_t sigErr = 0;
+  for (bin = 0; bin < Nbins; bin++)
+  {
+    if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
+    {
+      std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
+      if (rms == true)
+      {
+        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
+        gr.SetPointError(count,0, sigvalues.second); // to correct
+      }
+      else
+      {
+        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
+        gr.SetPointError(count,0, sigvalues.second);
+      }
+      count++;
+    }
+
+    else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10)
+    {
+      std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist);
+      if (rms == true)
+      {
+        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
+        gr.SetPointError(count,0, sigvalues.second); // to correct
+      }
+      else
+      {
+        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
+        gr.SetPointError(count,0, sigvalues.second);
+      }
+      count++;
+    }
+
+  }
+  return gr;
+}
+
+
+//------------------------------------------------------------------------------
+
+
+// type 1 : object, 2 : track, 3 : tower
+
+void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
+{
+
+  gr->SetLineWidth(2);
+
+  switch ( type )
+  {
+    case 1:
+      gr->SetLineColor(objColor);
+      gr->SetLineStyle(objStyle);
+      std::cout << "Adding " << gr->GetName() << std::endl;
+      mg->Add(gr);
+      leg->AddEntry(gr,"Reco","l");
+      break;
+
+    case 2:
+      gr->SetLineColor(trackColor);
+      gr->SetLineStyle(trackStyle);
+      mg->Add(gr);
+      leg->AddEntry(gr,"Track","l");
+      break;
+
+    case 3:
+      gr->SetLineColor(towerColor);
+      gr->SetLineStyle(towerStyle);
+      mg->Add(gr);
+      leg->AddEntry(gr,"Tower","l");
+      break;
+
+    case 0:
+      gr->SetLineColor(objColor);
+      gr->SetLineStyle(objStyle);
+      mg->Add(gr);
+      break;
+
+    default:
+      std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
+      break;
+  }
+}
+
+void addHist(TH1D *h, TLegend *leg, int type)
+{
+  h->SetLineWidth(2);
+
+  switch ( type )
+  {
+    case 1:
+      h->SetLineColor(objColor);
+      h->SetLineStyle(objStyle);
+      leg->AddEntry(h,"Reco","l");
+      break;
+
+    case 2:
+      h->SetLineColor(trackColor);
+      h->SetLineStyle(trackStyle);
+      leg->AddEntry(h,"Track","l");
+      break;
+
+    case 3:
+      h->SetLineColor(towerColor);
+      h->SetLineStyle(towerStyle);
+      leg->AddEntry(h,"Tower","l");
+      break;
+
+    case 0:
+      h->SetLineColor(objColor);
+      h->SetLineStyle(objStyle);
+      break;
+
+    default:
+      std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
+      break;
+  }
+}
+
+void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
+{
+  mg->SetMinimum(0.);
+  mg->SetMaximum(max);
+  mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
+  mg->GetYaxis()->SetTitle("resolution");
+  mg->GetXaxis()->SetTitle("p_{T} [GeV]");
+  mg->GetYaxis()->SetTitleSize(0.07);
+  mg->GetXaxis()->SetTitleSize(0.07);
+  mg->GetYaxis()->SetLabelSize(0.06);
+  mg->GetXaxis()->SetLabelSize(0.06);
+  mg->GetYaxis()->SetLabelOffset(0.03);
+  mg->GetYaxis()->SetTitleOffset(1.4);
+  mg->GetXaxis()->SetTitleOffset(1.4);
+
+  mg->GetYaxis()->SetNdivisions(505);
+
+  leg->SetBorderSize(0);
+  leg->SetShadowColor(0);
+  leg->SetFillColor(0);
+  leg->SetFillStyle(0);
+
+  gStyle->SetOptTitle(0);
+  gPad->SetLogx();
+  gPad->SetBottomMargin(0.2);
+  gPad->SetLeftMargin(0.2);
+  gPad->Modified();
+  gPad->Update();
+
+}
+
+void DrawAxis(TH1D *h, TLegend *leg, int type)
+{
+
+  h->GetYaxis()->SetRangeUser(0,1.0);
+  if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
+  else h->GetXaxis()->SetTitle("#eta");
+  h->GetYaxis()->SetTitle("efficiency");
+  h->GetYaxis()->SetTitleSize(0.07);
+  h->GetXaxis()->SetTitleSize(0.07);
+  h->GetYaxis()->SetLabelSize(0.06);
+  h->GetXaxis()->SetLabelSize(0.06);
+  h->GetYaxis()->SetLabelOffset(0.03);
+  h->GetYaxis()->SetTitleOffset(1.3);
+  h->GetXaxis()->SetTitleOffset(1.4);
+
+  h->GetYaxis()->SetNdivisions(505);
+
+  leg->SetBorderSize(0);
+  leg->SetShadowColor(0);
+  leg->SetFillColor(0);
+  leg->SetFillStyle(0);
+
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptStat(0);
+  gPad->SetBottomMargin(0.2);
+  gPad->SetLeftMargin(0.2);
+
+  gPad->Modified();
+  gPad->Update();
+
+}
+
+
+void Validation(const char *inputFileElectron, const char *inputFileMuon, const char *inputFilePhoton, const char *inputFileJet, const char *inputFileBJet, const char *inputFileTauJet, const char *outputFile)
+{
+  TChain *chainElectron = new TChain("Delphes");
+  chainElectron->Add(inputFileElectron);
+  ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
+
+  TChain *chainMuon = new TChain("Delphes");
+  chainMuon->Add(inputFileMuon);
+  ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
+
+  TChain *chainPhoton = new TChain("Delphes");
+  chainPhoton->Add(inputFilePhoton);
+  ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
+
+  TChain *chainJet = new TChain("Delphes");
+  chainJet->Add(inputFileJet);
+  ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
+
+  TChain *chainBJet = new TChain("Delphes");
+  chainBJet->Add(inputFileBJet);
+  ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
+
+  TChain *chainTauJet = new TChain("Delphes");
+  chainTauJet->Add(inputFileTauJet);
+  ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
+
+  TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
+  TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
+  TClonesArray *branchTowerElectron = treeReaderElectron->UseBranch("Tower");
+  TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
+
+  TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
+  TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
+  TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
+
+  TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
+  TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
+  TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
+
+  TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
+  TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet");
+  TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
+
+  TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
+  TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
+
+  TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
+  TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
+
+  TClonesArray *branchScalarHT = treeReaderJet->UseBranch("ScalarHT");
+  TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET");
+
+  ///////////////
+  // Electrons //
+  ///////////////
+
+  // Reconstruction efficiency
+  TString elecs = "Electron";
+  int elID = 11;
+  std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticleElectron, "Electron", elID, treeReaderElectron);
+
+  // tracking reconstruction efficiency
+  std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrackElectron, branchParticleElectron, "electronTrack", elID, treeReaderElectron);
+
+  // Tower reconstruction efficiency
+  std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTowerElectron, branchParticleElectron, "electronTower", elID, treeReaderElectron);
+
+  // Electron Energy Resolution
+  std::vector<resolPlot> plots_el;
+  HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
+  GetEres<Electron>(&plots_el, branchElectron, branchParticleElectron, elID, treeReaderElectron);
+  TGraphErrors gr_el = EresGraph(&plots_el, true);
+  TGraphErrors gr_elFwd = EresGraph(&plots_el, false);
+  gr_el.SetName("Electron");
+  gr_elFwd.SetName("ElectronFwd");
+
+  // Electron Track Energy Resolution
+  std::vector<resolPlot> plots_eltrack;
+  HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
+  GetEres<Track>(&plots_eltrack, branchTrackElectron, branchParticleElectron, elID, treeReaderElectron);
+  TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
+  TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false);
+  gr_eltrack.SetName("ElectronTracks");
+  gr_eltrackFwd.SetName("ElectronTracksFwd");
+
+  // Electron Tower Energy Resolution
+  std::vector<resolPlot> plots_eltower;
+  HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
+  GetEres<Tower>(&plots_eltower, branchTowerElectron, branchParticleElectron, elID, treeReaderElectron);
+  TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
+  TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false);
+  gr_eltower.SetName("ElectronTower");
+  gr_eltrackFwd.SetName("ElectronTracksFwd");
+
+  // Canvases
+  TString elEff = "electronEff";
+  TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600);
+  C_el1->Divide(2);
+  C_el1->cd(1);
+  TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
+  leg_el1->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_eltrack.first->Draw("][");
+  addHist(histos_eltrack.first, leg_el1, 2);
+  histos_el.first->Draw("same ][");
+  addHist(histos_el.first, leg_el1, 1);
+  DrawAxis(histos_eltrack.first, leg_el1, 0);
+
+  leg_el1->Draw();
+
+  C_el1->cd(2);
+  TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
+  leg_el2->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_eltrack.second->Draw("][");
+  addHist(histos_eltrack.second, leg_el2, 2);
+  histos_el.second->Draw("same ][");
+  addHist(histos_el.second, leg_el2, 1);
+
+  DrawAxis(histos_eltrack.second, leg_el2, 0);
+  leg_el2->Draw();
+
+  C_el1->cd(0);
+
+  TString elRes = "electronERes";
+  TString elResFwd = "electronEResForward";
+  TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600);
+  C_el2->Divide(2);
+  C_el2->cd(1);
+  TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
+  TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
+  leg_el->AddEntry("","","");
+
+  addGraph(mg_el, &gr_eltower, leg_el, 3);
+  addGraph(mg_el, &gr_eltrack, leg_el, 2);
+  addGraph(mg_el, &gr_el, leg_el, 1);
+
+  mg_el->Draw("ACX");
+  leg_el->Draw();
+
+  DrawAxis(mg_el, leg_el, 0.1);
+
+  C_el2->cd(2);
+  TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd);
+  TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
+  leg_elFwd->AddEntry("","","");
+
+  addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3);
+  addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2);
+  addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1);
+
+  mg_elFwd->Draw("ACX");
+  leg_elFwd->Draw();
+
+  DrawAxis(mg_elFwd, leg_elFwd, 0.2);
+
+  C_el1->Print("validation.pdf(","pdf");
+  C_el2->Print("validation.pdf","pdf");
+
+  gDirectory->cd(0);
+
+  ///////////
+  // Muons //
+  ///////////
+
+  // Reconstruction efficiency
+  int muID = 13;
+  std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon);
+
+  // muon tracking reconstruction efficiency
+  std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
+
+  // Muon Energy Resolution
+  std::vector<resolPlot> plots_mu;
+  HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
+  GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
+  TGraphErrors gr_mu = EresGraph(&plots_mu, true);
+  TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
+  gr_mu.SetName("Muon");
+  gr_muFwd.SetName("MuonFwd");
+
+  // Muon Track Energy Resolution
+  std::vector<resolPlot> plots_mutrack;
+  HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
+  GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
+  TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
+  TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
+  gr_mutrackFwd.SetName("MuonTracksFwd");
+
+  // Canvas
+
+  TString muEff = "muonEff";
+  TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
+  C_mu1->Divide(2);
+  C_mu1->cd(1);
+  TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
+  leg_mu1->AddEntry("","","");
+
+
+  gPad->SetLogx();
+  histos_mutrack.first->Draw("][");
+  addHist(histos_mutrack.first, leg_mu1, 2);
+  histos_mu.first->Draw("same ][");
+  addHist(histos_mu.first, leg_mu1, 1);
+
+  DrawAxis(histos_mutrack.first, leg_mu1, 0);
+
+  leg_mu1->Draw();
+
+  C_mu1->cd(2);
+  TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
+  leg_mu2->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_mutrack.second->Draw("][");
+  addHist(histos_mutrack.second, leg_mu2, 2);
+  histos_mu.second->Draw("same ][");
+  addHist(histos_mu.second, leg_mu2, 1);
+
+  DrawAxis(histos_mutrack.second, leg_mu2, 1);
+  leg_mu2->Draw();
+
+  TString muRes = "muonERes";
+  TString muResFwd = "muonEResFwd";
+
+  TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
+  C_mu->Divide(2);
+  C_mu->cd(1);
+  TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
+  TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
+  leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
+  leg_mu->AddEntry("","","");
+
+  addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
+  addGraph(mg_mu, &gr_mu, leg_mu, 1);
+
+  mg_mu->Draw("ACX");
+  leg_mu->Draw();
+
+  DrawAxis(mg_mu, leg_mu, 0.3);
+
+  C_mu->cd(2);
+  TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
+  TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
+  leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
+  leg_muFwd->AddEntry("","","");
+
+  addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
+  addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
+
+  mg_muFwd->Draw("ACX");
+  leg_muFwd->Draw();
+
+  DrawAxis(mg_muFwd, leg_muFwd, 0.3);
+
+  //C_mu1->SaveAs(muEff+".eps");
+  //C_mu->SaveAs(muRes+".eps");
+
+  C_mu1->Print("validation.pdf","pdf");
+  C_mu->Print("validation.pdf","pdf");
+
+  gDirectory->cd(0);
+
+  /////////////
+  // Photons //
+  /////////////
+
+  // Reconstruction efficiency
+  int phID = 22;
+  std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
+  std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
+
+  // Photon Energy Resolution
+  std::vector<resolPlot> plots_ph;
+  HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
+  GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton);
+  TGraphErrors gr_ph = EresGraph(&plots_ph, true);
+  TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
+  gr_ph.SetName("Photon");
+  gr_phFwd.SetName("PhotonFwd");
+
+
+  // Photon Tower Energy Resolution
+  std::vector<resolPlot> plots_phtower;
+  HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
+  GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton);
+  TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
+  TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
+  gr_phtower.SetName("PhotonTower");
+  gr_phtowerFwd.SetName("PhotonTowerFwd");
+
+  // Canvas
+
+  TString phEff = "photonEff";
+  TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
+  C_ph1->Divide(2);
+  C_ph1->cd(1);
+  TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
+  leg_ph1->AddEntry("","","");
+
+
+  gPad->SetLogx();
+  histos_phtower.first->Draw("][");
+  addHist(histos_phtower.first, leg_ph1, 3);
+  histos_ph.first->Draw("same ][");
+  addHist(histos_ph.first, leg_ph1, 1);
+
+  DrawAxis(histos_phtower.first, leg_ph1, 0);
+  leg_ph1->Draw();
+
+  C_ph1->cd(2);
+  TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
+  leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
+  leg_ph2->AddEntry("","","");
+
+
+  gPad->SetLogx();
+  histos_phtower.second->Draw("][");
+  addHist(histos_phtower.second, leg_ph2, 3);
+  histos_ph.second->Draw("same ][");
+  addHist(histos_ph.second, leg_ph2, 1);
+
+  DrawAxis(histos_phtower.second, leg_ph2, 1);
+  leg_ph2->Draw();
+
+  C_ph1->SaveAs(phEff+".eps");
+
+  TString phRes = "phERes";
+  TString phResFwd = "phEResFwd";
+
+  TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
+  C_ph->Divide(2);
+  C_ph->cd(1);
+  TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
+  TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
+  leg_ph->AddEntry("","","");
+
+  addGraph(mg_ph, &gr_phtower, leg_ph, 3);
+  addGraph(mg_ph, &gr_ph, leg_ph, 1);
+
+  mg_ph->Draw("ACX");
+  leg_ph->Draw();
+
+  DrawAxis(mg_ph, leg_ph, 0.3);
+
+  C_ph->cd(2);
+  TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
+  TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
+  leg_phFwd->AddEntry("","","");
+
+  addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
+  addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
+
+  mg_phFwd->Draw("ACX");
+  leg_phFwd->Draw();
+
+  DrawAxis(mg_phFwd, leg_phFwd, 0.3);
+
+  C_ph->SaveAs(phRes+".eps");
+
+  C_ph1->Print("validation.pdf","pdf");
+  C_ph->Print("validation.pdf","pdf");
+
+  gDirectory->cd(0);
+
+  //////////
+  // Jets //
+  //////////
+
+  // BJets Reconstruction efficiency
+  int bID = 5;
+  std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet);
+
+  // TauJets Reconstruction efficiency
+  int tauID = 15;
+  std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet);
+
+  // PFJets Energy Resolution
+  std::vector<resolPlot> plots_pfjets;
+  HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
+  GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet);
+  TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
+  TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
+  gr_pfjets.SetName("pfJet");
+  gr_pfjetsFwd.SetName("pfJetFwd");
+
+  // CaloJets Energy Resolution
+  std::vector<resolPlot> plots_calojets;
+  HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
+  GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet);
+  TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
+  TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
+  gr_calojets.SetName("caloJet");
+  gr_calojetsFwd.SetName("caloJetFwd");
+
+  // MET Resolution vs HT
+  std::vector<resolPlot> plots_met;
+  HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
+  GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet);
+  TGraphErrors gr_met = EresGraph(&plots_met, true);
+  gr_calojets.SetName("MET");
+
+  // Canvas
+  TString btagEff = "btagEff";
+  TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
+  C_btag1->Divide(2);
+  C_btag1->cd(1);
+  TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
+  leg_btag1->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_btag.first->Draw();
+  addHist(histos_btag.first, leg_btag1, 0);
+
+  DrawAxis(histos_btag.first, leg_btag1, 0);
+  leg_btag1->Draw();
+
+  C_btag1->cd(2);
+  TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
+  leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
+  leg_btag2->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_btag.second->Draw();
+  addHist(histos_btag.second, leg_btag2, 0);
+
+  DrawAxis(histos_btag.second, leg_btag2, 0);
+  leg_btag2->Draw();
+  C_btag1->cd(0);
+
+  TString tautagEff = "tautagEff";
+  TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
+  C_tautag1->Divide(2);
+  C_tautag1->cd(1);
+  TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
+  leg_tautag1->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_tautag.first->Draw();
+  addHist(histos_tautag.first, leg_tautag1, 0);
+
+  DrawAxis(histos_tautag.first, leg_tautag1, 0);
+  leg_tautag1->Draw();
+
+  C_tautag1->cd(2);
+  TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
+  leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
+  leg_tautag2->AddEntry("","","");
+
+  gPad->SetLogx();
+  histos_tautag.second->Draw();
+  addHist(histos_tautag.second, leg_tautag2, 0);
+
+  DrawAxis(histos_tautag.second, leg_tautag2, 0);
+  leg_tautag2->Draw();
+  C_tautag1->cd(0);
+
+  TString jetRes = "jetERes";
+  TString jetResFwd = "jetEResFwd";
+  TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
+  C_jet->Divide(2);
+
+  C_jet->cd(1);
+  TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
+  TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
+  leg_jet->AddEntry("","","");
+
+  addGraph(mg_jet, &gr_calojets, leg_jet, 3);
+  addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
+
+  mg_jet->Draw("ACX");
+  leg_jet->Draw();
+
+  DrawAxis(mg_jet, leg_jet, 0.25);
+
+  C_jet->cd(2);
+  TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
+  TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
+  leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
+  leg_jetFwd->AddEntry("","","");
+
+  addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
+  addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
+
+  mg_jetFwd->Draw("ACX");
+  leg_jetFwd->Draw();
+
+  DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
+
+  C_btag1->SaveAs(btagEff+".eps");
+  C_jet->SaveAs(jetRes+".eps");
+
+  TString metRes = "MetRes";
+  TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
+
+  TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
+  TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
+  leg_met->SetHeader("E_{T}^{miss}");
+  leg_met->AddEntry("","","");
+
+
+  addGraph(mg_met, &gr_met, leg_met, 0);
+
+  mg_met->Draw("ACX");
+  leg_met->Draw();
+
+  DrawAxis(mg_met, leg_met, 300);
+
+  mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
+  mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
+
+  C_met->SaveAs(metRes+".eps");
+
+  C_jet->Print("validation.pdf","pdf");
+  C_btag1->Print("validation.pdf","pdf");
+  C_tautag1->Print("validation.pdf","pdf");
+  C_met->Print("validation.pdf)","pdf");
+
+  TFile *fout = new TFile(outputFile,"recreate");
+
+  for (int bin = 0; bin < Nbins; bin++)
+  {
+    plots_el.at(bin).cenResolHist->Write();
+    plots_eltrack.at(bin).cenResolHist->Write();
+    plots_eltower.at(bin).cenResolHist->Write();
+    plots_el.at(bin).fwdResolHist->Write();
+    plots_eltrack.at(bin).fwdResolHist->Write();
+    plots_eltower.at(bin).fwdResolHist->Write();
+  }
+
+  histos_el.first->Write();
+  histos_el.second->Write();
+  histos_eltrack.first->Write();
+  histos_eltrack.second->Write();
+  histos_eltower.first->Write();
+  histos_eltower.second->Write();
+  C_el1->Write();
+  C_el2->Write();
+
+  histos_mu.first->Write();
+  histos_mu.second->Write();
+  histos_mutrack.first->Write();
+  histos_mutrack.second->Write();
+  C_mu1->Write();
+  C_mu->Write();
+
+  histos_ph.first->Write();
+  histos_ph.second->Write();
+  C_ph1->Write();
+  C_ph->Write();
+
+  for (int bin = 0; bin < Nbins; bin++)
+  {
+    plots_pfjets.at(bin).cenResolHist->Write();
+    plots_pfjets.at(bin).fwdResolHist->Write();
+    plots_calojets.at(bin).cenResolHist->Write();
+    plots_calojets.at(bin).fwdResolHist->Write();
+    plots_met.at(bin).cenResolHist->Write();
+  }
+  histos_btag.first->Write();
+  histos_btag.second->Write();
+  histos_tautag.first->Write();
+  histos_tautag.second->Write();
+  C_btag1->Write();
+  C_tautag1->Write();
+  C_jet->Write();
+  C_met->Write();
+
+  fout->Write();
+
+  cout << "** Exiting..." << endl;
+
+  delete treeReaderElectron;
+  delete treeReaderMuon;
+  delete treeReaderPhoton;
+  delete treeReaderJet;
+  delete treeReaderBJet;
+  delete treeReaderTauJet;
+  delete chainElectron;
+  delete chainMuon;
+  delete chainPhoton;
+  delete chainJet;
+  delete chainBJet;
+  delete chainTauJet;
+}
 
 //------------------------------------------------------------------------------
@@ -59 +1278 @@
   if(argc != 3)
   {
-    cout << " Usage: " << appName << " input_file output_file" << endl;
-    cout << " input_file  - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
+    cout << " input_file_electron  - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
+    cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
     cout << " output_file - output file in ROOT format" << endl;
     return 1;
@@ -71 +1295 @@
   TApplication app(appName, &appargc, appargv);
 
-  TString inputFile(argv[1]);
-  TString outputFile(argv[2]);
-
-
-//------------------------------------------------------------------------------
-// Here you call your macro's main function
-
-  Validation(inputFile, outputFile);
-
-//------------------------------------------------------------------------------
-
-}
-
-
+  Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
+}
+
+

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"
 #include "modules/Weighter.h"
@@ -63 +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"
 
@@ -80 +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+;
@@ -103 +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+;
 #pragma link C++ class Weighter+;
@@ -108 +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/PhotonConversions.cc

@@ -59 +59 @@
   fItInputArray(0), fConversionMap(0), fDecayXsec(0)
 {
-  fDecayXsec = new TF1;
+  fDecayXsec = new TF1("decayXsec","1.0 - 4.0/3.0 * x * (1.0 - x)", 0.0, 1.0);
   fConversionMap = new DelphesCylindricalFormula;
 }
@@ -84 +84 @@
 
   fConversionMap->Compile(GetString("ConversionMap", "0.0"));
-
-#if  ROOT_VERSION_CODE < ROOT_VERSION(6,04,00)
-  fDecayXsec->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
-#else
-  fDecayXsec->GetFormula()->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
-#endif
-  fDecayXsec->SetRange(0.0, 1.0);
 
   // import array with output from filter/classifier module

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

@@ -65 +65 @@
   ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchRwgt,
   DelphesFactory *factory, TObjArray *allParticleOutputArray,
-  TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray)
+  TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray, Bool_t firstEvent)
 {
+
   fwlite::Handle< GenEventInfoProduct > handleGenEventInfo;
-
   fwlite::Handle< LHEEventProduct > handleLHEEvent;
-
   fwlite::Handle< vector< reco::GenParticle > > handleParticle;
+
   vector< reco::GenParticle >::const_iterator itParticle;
 
@@ -78 +78 @@
 
   handleGenEventInfo.getByLabel(event, "generator");
-  handleLHEEvent.getByLabel(event, "externalLHEProducer");
-  handleParticle.getByLabel(event, "genParticles");
+
+  if (!((handleLHEEvent.getBranchNameFor(event, "source")).empty()))
+  { 
+    handleLHEEvent.getByLabel(event, "source");
+  }
+  else if (!((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty()))
+  {
+    handleLHEEvent.getByLabel(event, "externalLHEProducer");
+  }
+  else if (firstEvent)
+  {
+    std::cout<<"Wrong LHEEvent Label! Please, check the input file."<<std::endl;
+  }
+
+  if (!((handleParticle.getBranchNameFor(event, "genParticles")).empty()))
+  {
+    handleParticle.getByLabel(event, "genParticles");
+  }
+  else if (!((handleParticle.getBranchNameFor(event, "prunedGenParticles")).empty()))
+  {
+    handleParticle.getByLabel(event, "prunedGenParticles");
+  }
+  else
+  {
+    std::cout<<"Wrong GenParticle Label! Please, check the input file."<<std::endl;
+    exit(-1);
+  }
+
+  Bool_t foundLHE = !((handleLHEEvent.getBranchNameFor(event, "source")).empty()) || !((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty());
 
   HepMCEvent *element;
@@ -91 +118 @@
   Double_t px, py, pz, e, mass;
   Double_t x, y, z;
-
-  const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
-  vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
 
   element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
@@ -117 +141 @@
   element->ProcTime = 0.0;
 
-  for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
-  {
-    weight = static_cast<Weight *>(branchRwgt->NewEntry());
-    weight->Weight = itWeightsInfo->wgt;
+
+  if(foundLHE)
+  {
+    const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
+    vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
+    
+    for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
+    {
+      weight = static_cast<Weight *>(branchRwgt->NewEntry());
+      weight->Weight = itWeightsInfo->wgt;
+    }  
   }
 
@@ -207 +238 @@
   Int_t i;
   Long64_t eventCounter, numberOfEvents;
+  Bool_t firstEvent = kTRUE;
 
   if(argc < 4)
@@ -240 +272 @@
 
     branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
-    branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
+    branchRwgt = treeWriter->NewBranch("Weight", Weight::Class());
 
     confReader = new ExRootConfReader;
@@ -281 +313 @@
       modularDelphes->Clear();
       treeWriter->Clear();
+
       for(event.toBegin(); !event.atEnd() && !interrupted; ++event)
       {
         ConvertInput(event, eventCounter, branchEvent, branchRwgt, factory,
-          allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
+          allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
         modularDelphes->ProcessTask();
+          
+        firstEvent = kFALSE;
 
         treeWriter->Fill();