Fork me on GitHub

Changeset 70b9632 in git


Ignore:
Timestamp:
Aug 25, 2016, 2:04:58 PM (8 years ago)
Author:
Alexandre Mertens <alexandre.mertens@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
629e819
Parents:
7bb13cd (diff), 1408174 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

dev_01

Files:
15 added
28 edited

Legend:

Unmodified
Added
Removed

  • CHANGELOG

    r7bb13cd r70b9632  
     13.3.3:
     2- improved compatibility with ROOT >= 6.04
     3- removed test of praticle stability based on ROOT pdgtable (#347, #753, #821, #855)
     4- fixed bugs in DelphesCMSFWLite
     5- fixed bugs in PUPPI
     6- fixed compiler and linker flags for Pythia8
     7- added CMS Phase II cards
     8- added particle gun to DelphesPythia8
     9- added jet merging to DelphesPythia8 based on Pythia8 main32.cc
     10- added UseMiniCone option in Isolation module
     11- added RecoPuFilter module
     12- added back OldCalorimeter for CheckMate bwd compatibilty (#591)
     13- updated tracking and calorimeter resolution in CMS cards
     14
    1153.3.2:
    216- added pre-validated card for CMS upgrade Phase II studies

  • Makefile

    r7bb13cd r70b9632  
    250250        external/ExRootAnalysis/ExRootTreeBranch.h \
    251251        external/ExRootAnalysis/ExRootProgressBar.h
     252DelphesROOT$(ExeSuf): \
     253        tmp/readers/DelphesROOT.$(ObjSuf)
     254
     255tmp/readers/DelphesROOT.$(ObjSuf): \
     256        readers/DelphesROOT.cpp \
     257        modules/Delphes.h \
     258        classes/DelphesStream.h \
     259        classes/DelphesClasses.h \
     260        classes/DelphesFactory.h \
     261        external/ExRootAnalysis/ExRootTreeWriter.h \
     262        external/ExRootAnalysis/ExRootTreeReader.h \
     263        external/ExRootAnalysis/ExRootTreeBranch.h \
     264        external/ExRootAnalysis/ExRootProgressBar.h
    252265DelphesSTDHEP$(ExeSuf): \
    253266        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    … …  
    265278        DelphesHepMC$(ExeSuf) \
    266279        DelphesLHEF$(ExeSuf) \
     280        DelphesROOT$(ExeSuf) \
    267281        DelphesSTDHEP$(ExeSuf)
    268282
    … …  
    270284        tmp/readers/DelphesHepMC.$(ObjSuf) \
    271285        tmp/readers/DelphesLHEF.$(ObjSuf) \
     286        tmp/readers/DelphesROOT.$(ObjSuf) \
    272287        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    273288
    … …  
    384399        modules/EnergySmearing.h \
    385400        modules/MomentumSmearing.h \
     401        modules/TrackSmearing.h \
    386402        modules/ImpactParameterSmearing.h \
    387403        modules/TimeSmearing.h \
    … …  
    407423        modules/StatusPidFilter.h \
    408424        modules/PdgCodeFilter.h \
     425        modules/BeamSpotFilter.h \
     426        modules/RecoPuFilter.h \
    409427        modules/Cloner.h \
    410428        modules/Weighter.h \
    … …  
    412430        modules/JetFlavorAssociation.h \
    413431        modules/JetFakeParticle.h \
     432        modules/VertexSorter.h \
     433        modules/VertexFinder.h \
     434        modules/VertexFinderDA4D.h \
    414435        modules/ExampleModule.h
    415436tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    618639        classes/DelphesFactory.h \
    619640        classes/DelphesFormula.h
     641tmp/modules/BeamSpotFilter.$(ObjSuf): \
     642        modules/BeamSpotFilter.$(SrcSuf) \
     643        modules/BeamSpotFilter.h \
     644        classes/DelphesClasses.h \
     645        classes/DelphesFactory.h \
     646        classes/DelphesFormula.h \
     647        external/ExRootAnalysis/ExRootResult.h \
     648        external/ExRootAnalysis/ExRootFilter.h \
     649        external/ExRootAnalysis/ExRootClassifier.h
    620650tmp/modules/Calorimeter.$(ObjSuf): \
    621651        modules/Calorimeter.$(SrcSuf) \
    … …  
    850880        external/ExRootAnalysis/ExRootFilter.h \
    851881        external/ExRootAnalysis/ExRootClassifier.h
     882tmp/modules/RecoPuFilter.$(ObjSuf): \
     883        modules/RecoPuFilter.$(SrcSuf) \
     884        modules/RecoPuFilter.h \
     885        classes/DelphesClasses.h \
     886        classes/DelphesFactory.h \
     887        classes/DelphesFormula.h \
     888        external/ExRootAnalysis/ExRootResult.h \
     889        external/ExRootAnalysis/ExRootFilter.h \
     890        external/ExRootAnalysis/ExRootClassifier.h
    852891tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    853892        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    911950        modules/TrackPileUpSubtractor.$(SrcSuf) \
    912951        modules/TrackPileUpSubtractor.h \
     952        classes/DelphesClasses.h \
     953        classes/DelphesFactory.h \
     954        classes/DelphesFormula.h \
     955        external/ExRootAnalysis/ExRootResult.h \
     956        external/ExRootAnalysis/ExRootFilter.h \
     957        external/ExRootAnalysis/ExRootClassifier.h
     958tmp/modules/TrackSmearing.$(ObjSuf): \
     959        modules/TrackSmearing.$(SrcSuf) \
     960        modules/TrackSmearing.h \
    913961        classes/DelphesClasses.h \
    914962        classes/DelphesFactory.h \
    … …  
    933981        classes/DelphesFactory.h \
    934982        classes/DelphesFormula.h \
     983        external/ExRootAnalysis/ExRootResult.h \
     984        external/ExRootAnalysis/ExRootFilter.h \
     985        external/ExRootAnalysis/ExRootClassifier.h
     986tmp/modules/VertexFinder.$(ObjSuf): \
     987        modules/VertexFinder.$(SrcSuf) \
     988        modules/VertexFinder.h \
     989        classes/DelphesClasses.h \
     990        classes/DelphesFactory.h \
     991        classes/DelphesFormula.h \
     992        classes/DelphesPileUpReader.h \
     993        external/ExRootAnalysis/ExRootResult.h \
     994        external/ExRootAnalysis/ExRootFilter.h \
     995        external/ExRootAnalysis/ExRootClassifier.h
     996tmp/modules/VertexFinderDA4D.$(ObjSuf): \
     997        modules/VertexFinderDA4D.$(SrcSuf) \
     998        modules/VertexFinderDA4D.h \
     999        classes/DelphesClasses.h \
     1000        classes/DelphesFactory.h \
     1001        classes/DelphesFormula.h \
     1002        classes/DelphesPileUpReader.h \
     1003        external/ExRootAnalysis/ExRootResult.h \
     1004        external/ExRootAnalysis/ExRootFilter.h \
     1005        external/ExRootAnalysis/ExRootClassifier.h
     1006tmp/modules/VertexSorter.$(ObjSuf): \
     1007        modules/VertexSorter.$(SrcSuf) \
     1008        modules/VertexSorter.h \
     1009        classes/DelphesClasses.h \
     1010        classes/DelphesFactory.h \
     1011        classes/DelphesFormula.h \
     1012        classes/DelphesPileUpReader.h \
    9351013        external/ExRootAnalysis/ExRootResult.h \
    9361014        external/ExRootAnalysis/ExRootFilter.h \
    … …  
    9961074        tmp/modules/AngularSmearing.$(ObjSuf) \
    9971075        tmp/modules/BTagging.$(ObjSuf) \
     1076        tmp/modules/BeamSpotFilter.$(ObjSuf) \
    9981077        tmp/modules/Calorimeter.$(ObjSuf) \
    9991078        tmp/modules/Cloner.$(ObjSuf) \
    … …  
    10201099        tmp/modules/PileUpJetID.$(ObjSuf) \
    10211100        tmp/modules/PileUpMerger.$(ObjSuf) \
     1101        tmp/modules/RecoPuFilter.$(ObjSuf) \
    10221102        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    10231103        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    10281108        tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
    10291109        tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
     1110        tmp/modules/TrackSmearing.$(ObjSuf) \
    10301111        tmp/modules/TreeWriter.$(ObjSuf) \
    10311112        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
     1113        tmp/modules/VertexFinder.$(ObjSuf) \
     1114        tmp/modules/VertexFinderDA4D.$(ObjSuf) \
     1115        tmp/modules/VertexSorter.$(ObjSuf) \
    10321116        tmp/modules/Weighter.$(ObjSuf)
    10331117
    … …  
    16341718        tmp/external/tcl/tclVar.$(ObjSuf)
    16351719
     1720modules/VertexFinderDA4D.h: \
     1721        classes/DelphesModule.h \
     1722        classes/DelphesClasses.h
     1723        @touch $@
     1724
     1725modules/TrackSmearing.h: \
     1726        classes/DelphesModule.h
     1727        @touch $@
     1728
    16361729external/fastjet/ClusterSequence.hh: \
    16371730        external/fastjet/PseudoJet.hh \
    … …  
    18941987        @touch $@
    18951988
     1989modules/VertexSorter.h: \
     1990        classes/DelphesModule.h \
     1991        classes/DelphesClasses.h
     1992        @touch $@
     1993
    18961994modules/Delphes.h: \
    18971995        classes/DelphesModule.h
     1996        @touch $@
     1997
     1998modules/VertexFinder.h: \
     1999        classes/DelphesModule.h \
     2000        classes/DelphesClasses.h
    18982001        @touch $@
    18992002
    … …  
    19672070        @touch $@
    19682071
     2072modules/RecoPuFilter.h: \
     2073        classes/DelphesModule.h
     2074        @touch $@
     2075
    19692076modules/Hector.h: \
    19702077        classes/DelphesModule.h
    … …  
    20662173
    20672174modules/FastJetFinder.h: \
     2175        classes/DelphesModule.h
     2176        @touch $@
     2177
     2178modules/BeamSpotFilter.h: \
    20682179        classes/DelphesModule.h
    20692180        @touch $@

  • README

    r7bb13cd r70b9632  
    44Commands to get the code:
    55
    6    wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.3.2.tar.gz
     6   wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.3.3.tar.gz
    77
    8    tar -zxf Delphes-3.3.2.tar.gz
     8   tar -zxf Delphes-3.3.3.tar.gz
    99
    1010Commands to compile the code:
    1111
    12    cd Delphes-3.3.2
     12   cd Delphes-3.3.3
    1313
    1414   make

  • VERSION

    r7bb13cd r70b9632  
    1 3.3.2
     13.3.3

  • cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_0PU.tcl

    r7bb13cd r70b9632  
    3333  PhotonEnergySmearing
    3434  ElectronFilter
     35
    3536  TrackPileUpSubtractor
     37  RecoPuFilter
    3638
    3739  TowerMerger
    3840  NeutralEFlowMerger
    39   EFlowMergerAllTracks
     41
    4042  EFlowMerger
     43  EFlowMergerCHS
     44  Rho
    4145
    4246  LeptonFilterNoLep
    … …  
    4852
    4953  PhotonIsolation
     54  PhotonIsolationCHS
    5055  PhotonEfficiency
     56  PhotonEfficiencyCHS
    5157
    5258  ElectronIsolation
     59  ElectronIsolationCHS
     60
    5361  ElectronEfficiency
     62  ElectronEfficiencyCHS
    5463
    5564  MuonIsolation
     65  MuonIsolationCHS
     66
    5667  MuonLooseIdEfficiency
    5768  MuonTightIdEfficiency
     69
     70  MuonLooseIdEfficiencyCHS
     71  MuonTightIdEfficiencyCHS
    5872
    5973  NeutrinoFilter
    … …  
    6882  FastJetFinder
    6983  FastJetFinderAK8
     84  JetPileUpSubtractor
     85  JetPileUpSubtractorAK8
    7086  FastJetFinderPUPPI
    7187  FastJetFinderPUPPIAK8
    … …  
    121137
    122138  # maximum spread in the beam direction in m
    123   set ZVertexSpread 0
     139  set ZVertexSpread 0.25
    124140
    125141  # maximum spread in time in s
    126   set TVertexSpread 0
     142  set TVertexSpread 800E-12
    127143
    128144  # 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))}
    129   set VertexDistributionFormula 1
     145  set VertexDistributionFormula {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
    130146
    131147}
    … …  
    273289  source muonMomentumResolution.tcl
    274290}
    275 
    276 
    277 
    278291
    279292##############
    … …  
    527540}
    528541
     542########################
     543# Reco PU filter
     544########################
     545
     546module RecoPuFilter RecoPuFilter {
     547  set InputArray HCal/eflowTracks
     548  set OutputArray eflowTracks
     549}
    529550
    530551###################################################
    … …  
    539560}
    540561
    541 
    542562####################
    543563# Neutral eflow erger
    … …  
    551571}
    552572
    553 
    554 ####################
     573#####################
    555574# Energy flow merger
    556 ####################
     575#####################
    557576
    558577module Merger EFlowMerger {
    … …  
    564583}
    565584
    566 ##################################
    567 # Energy flow merger (all tracks)
    568 ##################################
    569 
    570 module Merger EFlowMergerAllTracks {
     585############################
     586# Energy flow merger no PU
     587############################
     588
     589module Merger EFlowMergerCHS {
    571590# add InputArray InputArray
    572   add InputArray TrackMerger/tracks
     591  add InputArray RecoPuFilter/eflowTracks
    573592  add InputArray PhotonEnergySmearing/eflowPhotons
    574593  add InputArray HCal/eflowNeutralHadrons
    … …  
    695714}
    696715
    697 
    698716#####################
    699717# MC truth jet finder
    … …  
    735753}
    736754
    737 ############
     755
     756#############
     757# Rho pile-up
     758#############
     759
     760module FastJetFinder Rho {
     761#  set InputArray Calorimeter/towers
     762  set InputArray EFlowMergerCHS/eflow
     763
     764  set ComputeRho true
     765  set RhoOutputArray rho
     766
     767  # 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
     768  set AreaAlgorithm 5
     769
     770  # jet algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
     771  set JetAlgorithm 4
     772  set ParameterR 0.4
     773  set GhostEtaMax 5.0
     774
     775  add RhoEtaRange -5.0 -4.0
     776  add RhoEtaRange -4.0 -1.5
     777  add RhoEtaRange -1.5 1.5
     778  add RhoEtaRange 1.5 4.0
     779  add RhoEtaRange 4.0 5.0
     780
     781  set JetPTMin 0.0
     782}
     783
     784
     785##############
    738786# Jet finder
    739 ############
     787##############
    740788
    741789module FastJetFinder FastJetFinder {
    742790#  set InputArray TowerMerger/towers
    743   set InputArray EFlowMerger/eflow
     791  set InputArray EFlowMergerCHS/eflow
    744792
    745793  set OutputArray jets
     794
     795  set AreaAlgorithm 5
    746796
    747797  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
    … …  
    755805module FastJetFinder FastJetFinderAK8 {
    756806#  set InputArray TowerMerger/towers
    757   set InputArray EFlowMerger/eflow
     807  set InputArray EFlowMergerCHS/eflow
    758808
    759809  set OutputArray jets
     810
     811  set AreaAlgorithm 5
    760812
    761813  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
    … …  
    784836}
    785837
     838###########################
     839# Jet Pile-Up Subtraction
     840###########################
     841
     842module JetPileUpSubtractor JetPileUpSubtractor {
     843  set JetInputArray FastJetFinder/jets
     844  set RhoInputArray Rho/rho
     845
     846  set OutputArray jets
     847
     848  set JetPTMin 15.0
     849}
     850
     851##############################
     852# Jet Pile-Up Subtraction AK8
     853##############################
     854
     855module JetPileUpSubtractor JetPileUpSubtractorAK8 {
     856  set JetInputArray FastJetFinderAK8/jets
     857  set RhoInputArray Rho/rho
     858
     859  set OutputArray jets
     860
     861  set JetPTMin 15.0
     862}
     863
    786864module FastJetFinder FastJetFinderPUPPI {
    787865#  set InputArray TowerMerger/towers
    … …  
    833911
    834912module EnergyScale JetEnergyScale {
    835   set InputArray FastJetFinder/jets
     913  set InputArray JetPileUpSubtractor/jets
    836914  set OutputArray jets
    837915
    … …  
    841919
    842920module EnergyScale JetEnergyScaleAK8 {
    843   set InputArray FastJetFinderAK8/jets
     921  set InputArray JetPileUpSubtractorAK8/jets
    844922  set OutputArray jets
    845923
    … …  
    907985}
    908986
     987
     988########################
     989# Photon isolation CHS #
     990########################
     991
     992module Isolation PhotonIsolationCHS {
     993
     994  # particle for which calculate the isolation
     995  set CandidateInputArray PhotonFilter/photons
     996
     997  # isolation collection
     998  set IsolationInputArray EFlowMerger/eflow
     999
     1000  # output array
     1001  set OutputArray photons
     1002
     1003  # isolation cone
     1004  set DeltaRMax 0.3
     1005
     1006  # minimum pT
     1007  set PTMin     1.0
     1008
     1009  # iso ratio to cut
     1010  set PTRatioMax 9999.
     1011
     1012}
    9091013
    9101014
    … …  
    9291033
    9301034
     1035#####################
     1036# Photon efficiency #
     1037#####################
     1038
     1039module Efficiency PhotonEfficiencyCHS {
     1040
     1041  ## input particles
     1042  set InputArray PhotonIsolationCHS/photons
     1043  ## output particles
     1044  set OutputArray photons
     1045  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
     1046  # efficiency formula for photons
     1047  set EfficiencyFormula {                      (pt <= 10.0) * (0.00) + \
     1048                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.9635) + \
     1049         (abs(eta) > 1.5 && abs(eta) <= 4.0) * (pt > 10.0)  * (0.9624) + \
     1050         (abs(eta) > 4.0)                                   * (0.00)}
     1051
     1052}
     1053
    9311054######################
    9321055# Electron isolation #
    … …  
    9491072}
    9501073
     1074
     1075##########################
     1076# Electron isolation CHS #
     1077##########################
     1078
     1079module Isolation ElectronIsolationCHS {
     1080
     1081  set CandidateInputArray ElectronFilter/electrons
     1082
     1083  # isolation collection
     1084  set IsolationInputArray EFlowMerger/eflow
     1085
     1086  set OutputArray electrons
     1087
     1088  set DeltaRMax 0.3
     1089  set PTMin 1.0
     1090  set PTRatioMax 9999.
     1091
     1092}
    9511093
    9521094
    … …  
    9951137}
    9961138
     1139###########################
     1140# Electron efficiency CHS #
     1141###########################
     1142
     1143module Efficiency ElectronEfficiencyCHS {
     1144
     1145  set InputArray ElectronIsolationCHS/electrons
     1146  set OutputArray electrons
     1147
     1148  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
     1149  # efficiency formula for electrons
     1150  set EfficiencyFormula {
     1151                                      (pt <= 4.0)  * (0.00) + \
     1152                         (abs(eta) <= 1.45 ) * (pt >  4.0 && pt <= 6.0)   * (0.50) + \
     1153                         (abs(eta) <= 1.45 ) * (pt >  6.0 && pt <= 8.0)   * (0.70) + \
     1154                         (abs(eta) <= 1.45 ) * (pt >  8.0 && pt <= 10.0)  * (0.85) + \
     1155                         (abs(eta) <= 1.45 ) * (pt > 10.0 && pt <= 30.0)  * (0.94) + \
     1156                         (abs(eta) <= 1.45 ) * (pt > 30.0 && pt <= 50.0)  * (0.97) + \
     1157                         (abs(eta) <= 1.45 ) * (pt > 50.0 && pt <= 70.0)  * (0.98) + \
     1158                         (abs(eta) <= 1.45 ) * (pt > 70.0 )  * (1.0) + \
     1159                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt >  4.0 && pt <= 10.0)   * (0.35) + \
     1160                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 10.0 && pt <= 30.0)   * (0.40) + \
     1161                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 30.0 && pt <= 70.0)   * (0.45) + \
     1162                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 70.0 )  * (0.55) + \
     1163                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt >  4.0 && pt <= 10.0)  * (0.75) + \
     1164                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 10.0 && pt <= 30.0)  * (0.85) + \
     1165                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 30.0 && pt <= 50.0)  * (0.95) + \
     1166                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 50.0 && pt <= 70.0)  * (0.95) + \
     1167                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 70.0 )  * (1.0) + \
     1168                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt >  4.0 && pt <= 10.0)  * (0.65) + \
     1169                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 10.0 && pt <= 30.0)  * (0.75) + \
     1170                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 30.0 && pt <= 50.0)  * (0.90) + \
     1171                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 50.0 && pt <= 70.0)  * (0.90) + \
     1172                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 70.0 )  * (0.90) + \
     1173                         (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 4.0 && pt <= 10.0) * (0.65) + \
     1174                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 10.0 && pt <= 30.0) * (0.75) + \
     1175                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 30.0 && pt <= 50.0) * (0.90) + \
     1176                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 50.0 && pt <= 70.0) * (0.90) + \
     1177                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 70.0 ) * (0.90) + \
     1178                                          (abs(eta) > 4.0) * (0.00)
     1179
     1180  }
     1181}
     1182
     1183
    9971184##################
    9981185# Muon isolation #
    … …  
    10131200}
    10141201
    1015 
    1016 ##################
    1017 # Muon Loose Id  #
    1018 ##################
     1202######################
     1203# Muon isolation CHS #
     1204######################
     1205
     1206module Isolation MuonIsolationCHS {
     1207  set CandidateInputArray MuonMomentumSmearing/muons
     1208
     1209  # isolation collection
     1210  set IsolationInputArray EFlowMerger/eflow
     1211
     1212  set OutputArray muons
     1213
     1214  set DeltaRMax 0.3
     1215  set PTMin 1.0
     1216  set PTRatioMax 9999.
     1217
     1218}
     1219
     1220
     1221#####################
     1222# Muon Loose Id     #
     1223#####################
    10191224
    10201225module Efficiency MuonLooseIdEfficiency {
    … …  
    10381243}
    10391244
     1245
     1246#####################
     1247# Muon Loose Id CHS #
     1248#####################
     1249
     1250module Efficiency MuonLooseIdEfficiencyCHS {
     1251    set InputArray MuonIsolationCHS/muons
     1252    set OutputArray muons
     1253    # tracking + LooseID efficiency formula for muons
     1254    source muonLooseId.tcl
     1255
     1256}
     1257
     1258
     1259######################
     1260# Muon Tight Id  CHS #
     1261######################
     1262
     1263module Efficiency MuonTightIdEfficiencyCHS {
     1264    set InputArray MuonIsolationCHS/muons
     1265    set OutputArray muons
     1266    # tracking + TightID efficiency formula for muons
     1267    source muonTightId.tcl
     1268}
    10401269
    10411270
    … …  
    12501479module StatusPidFilter GenParticleFilter {
    12511480
    1252     set InputArray  Delphes/allParticles
     1481    set InputArray Delphes/allParticles
    12531482    set OutputArray filteredParticles
    12541483    set PTMin 5.0
    … …  
    12781507  add Branch MuonLooseIdEfficiency/muons MuonLoose Muon
    12791508  add Branch MuonTightIdEfficiency/muons MuonTight Muon
     1509
     1510  add Branch PhotonEfficiencyCHS/photons PhotonCHS Photon
     1511  add Branch ElectronEfficiencyCHS/electrons ElectronCHS Electron
     1512  add Branch MuonLooseIdEfficiencyCHS/muons MuonLooseCHS Muon
     1513  add Branch MuonTightIdEfficiencyCHS/muons MuonTightCHS Muon
    12801514
    12811515  add Branch JetEnergyScale/jets Jet Jet
    … …  
    12891523  add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
    12901524  add Branch ScalarHT/energy ScalarHT ScalarHT
    1291 }
     1525
     1526}

  • cards/CMS_PhaseII/CMS_PhaseII_Substructure_PIX4022_200PU.tcl

    r7bb13cd r70b9632  
    3333  PhotonEnergySmearing
    3434  ElectronFilter
     35
    3536  TrackPileUpSubtractor
     37  RecoPuFilter
    3638
    3739  TowerMerger
    3840  NeutralEFlowMerger
    39   EFlowMergerAllTracks
     41
    4042  EFlowMerger
     43  EFlowMergerCHS
     44  Rho
    4145
    4246  LeptonFilterNoLep
    … …  
    4852
    4953  PhotonIsolation
     54  PhotonIsolationCHS
    5055  PhotonEfficiency
     56  PhotonEfficiencyCHS
    5157
    5258  ElectronIsolation
     59  ElectronIsolationCHS
     60
    5361  ElectronEfficiency
     62  ElectronEfficiencyCHS
    5463
    5564  MuonIsolation
     65  MuonIsolationCHS
     66
    5667  MuonLooseIdEfficiency
    5768  MuonTightIdEfficiency
     69
     70  MuonLooseIdEfficiencyCHS
     71  MuonTightIdEfficiencyCHS
    5872
    5973  NeutrinoFilter
    … …  
    6882  FastJetFinder
    6983  FastJetFinderAK8
     84  JetPileUpSubtractor
     85  JetPileUpSubtractorAK8
    7086  FastJetFinderPUPPI
    7187  FastJetFinderPUPPIAK8
    … …  
    126142  set TVertexSpread 800E-12
    127143
    128   # vertex smearing formula f(z,t) (z,t need to be respectively given in m,s)
     144  # 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))}
    129145  set VertexDistributionFormula {exp(-(t^2/160e-12^2/2))*exp(-(z^2/0.053^2/2))}
    130146
    … …  
    138154module ParticlePropagator ParticlePropagator {
    139155  set InputArray PileUpMerger/stableParticles
     156  #set InputArray Delphes/stableParticles
    140157
    141158  set OutputArray stableParticles
    … …  
    272289  source muonMomentumResolution.tcl
    273290}
    274 
    275 
    276 
    277291
    278292##############
    … …  
    526540}
    527541
     542########################
     543# Reco PU filter
     544########################
     545
     546module RecoPuFilter RecoPuFilter {
     547  set InputArray HCal/eflowTracks
     548  set OutputArray eflowTracks
     549}
    528550
    529551###################################################
    … …  
    538560}
    539561
    540 
    541562####################
    542563# Neutral eflow erger
    … …  
    550571}
    551572
    552 
    553 ####################
     573#####################
    554574# Energy flow merger
    555 ####################
     575#####################
    556576
    557577module Merger EFlowMerger {
    … …  
    563583}
    564584
    565 ##################################
    566 # Energy flow merger (all tracks)
    567 ##################################
    568 
    569 module Merger EFlowMergerAllTracks {
     585############################
     586# Energy flow merger no PU
     587############################
     588
     589module Merger EFlowMergerCHS {
    570590# add InputArray InputArray
    571   add InputArray TrackMerger/tracks
     591  add InputArray RecoPuFilter/eflowTracks
    572592  add InputArray PhotonEnergySmearing/eflowPhotons
    573593  add InputArray HCal/eflowNeutralHadrons
    … …  
    694714}
    695715
    696 
    697716#####################
    698717# MC truth jet finder
    … …  
    735754
    736755
    737 
    738 ############
     756#############
     757# Rho pile-up
     758#############
     759
     760module FastJetFinder Rho {
     761#  set InputArray Calorimeter/towers
     762  set InputArray EFlowMergerCHS/eflow
     763
     764  set ComputeRho true
     765  set RhoOutputArray rho
     766
     767  # 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
     768  set AreaAlgorithm 5
     769
     770  # jet algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
     771  set JetAlgorithm 4
     772  set ParameterR 0.4
     773  set GhostEtaMax 5.0
     774
     775  add RhoEtaRange -5.0 -4.0
     776  add RhoEtaRange -4.0 -1.5
     777  add RhoEtaRange -1.5 1.5
     778  add RhoEtaRange 1.5 4.0
     779  add RhoEtaRange 4.0 5.0
     780
     781  set JetPTMin 0.0
     782}
     783
     784
     785##############
    739786# Jet finder
    740 ############
     787##############
    741788
    742789module FastJetFinder FastJetFinder {
    743790#  set InputArray TowerMerger/towers
    744   set InputArray EFlowMerger/eflow
     791  set InputArray EFlowMergerCHS/eflow
    745792
    746793  set OutputArray jets
     794
     795  set AreaAlgorithm 5
    747796
    748797  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
    … …  
    756805module FastJetFinder FastJetFinderAK8 {
    757806#  set InputArray TowerMerger/towers
    758   set InputArray EFlowMerger/eflow
     807  set InputArray EFlowMergerCHS/eflow
    759808
    760809  set OutputArray jets
     810
     811  set AreaAlgorithm 5
    761812
    762813  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
    … …  
    785836}
    786837
     838###########################
     839# Jet Pile-Up Subtraction
     840###########################
     841
     842module JetPileUpSubtractor JetPileUpSubtractor {
     843  set JetInputArray FastJetFinder/jets
     844  set RhoInputArray Rho/rho
     845
     846  set OutputArray jets
     847
     848  set JetPTMin 15.0
     849}
     850
     851##############################
     852# Jet Pile-Up Subtraction AK8
     853##############################
     854
     855module JetPileUpSubtractor JetPileUpSubtractorAK8 {
     856  set JetInputArray FastJetFinderAK8/jets
     857  set RhoInputArray Rho/rho
     858
     859  set OutputArray jets
     860
     861  set JetPTMin 15.0
     862}
     863
    787864module FastJetFinder FastJetFinderPUPPI {
    788865#  set InputArray TowerMerger/towers
    … …  
    834911
    835912module EnergyScale JetEnergyScale {
    836   set InputArray FastJetFinder/jets
     913  set InputArray JetPileUpSubtractor/jets
    837914  set OutputArray jets
    838915
    … …  
    842919
    843920module EnergyScale JetEnergyScaleAK8 {
    844   set InputArray FastJetFinderAK8/jets
     921  set InputArray JetPileUpSubtractorAK8/jets
    845922  set OutputArray jets
    846923
    … …  
    908985}
    909986
     987
     988########################
     989# Photon isolation CHS #
     990########################
     991
     992module Isolation PhotonIsolationCHS {
     993
     994  # particle for which calculate the isolation
     995  set CandidateInputArray PhotonFilter/photons
     996
     997  # isolation collection
     998  set IsolationInputArray EFlowMerger/eflow
     999
     1000  # output array
     1001  set OutputArray photons
     1002
     1003  # isolation cone
     1004  set DeltaRMax 0.3
     1005
     1006  # minimum pT
     1007  set PTMin     1.0
     1008
     1009  # iso ratio to cut
     1010  set PTRatioMax 9999.
     1011
     1012}
    9101013
    9111014
    … …  
    9301033
    9311034
     1035#####################
     1036# Photon efficiency #
     1037#####################
     1038
     1039module Efficiency PhotonEfficiencyCHS {
     1040
     1041  ## input particles
     1042  set InputArray PhotonIsolationCHS/photons
     1043  ## output particles
     1044  set OutputArray photons
     1045  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
     1046  # efficiency formula for photons
     1047  set EfficiencyFormula {                      (pt <= 10.0) * (0.00) + \
     1048                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.9635) + \
     1049         (abs(eta) > 1.5 && abs(eta) <= 4.0) * (pt > 10.0)  * (0.9624) + \
     1050         (abs(eta) > 4.0)                                   * (0.00)}
     1051
     1052}
     1053
    9321054######################
    9331055# Electron isolation #
    … …  
    9501072}
    9511073
     1074
     1075##########################
     1076# Electron isolation CHS #
     1077##########################
     1078
     1079module Isolation ElectronIsolationCHS {
     1080
     1081  set CandidateInputArray ElectronFilter/electrons
     1082
     1083  # isolation collection
     1084  set IsolationInputArray EFlowMerger/eflow
     1085
     1086  set OutputArray electrons
     1087
     1088  set DeltaRMax 0.3
     1089  set PTMin 1.0
     1090  set PTRatioMax 9999.
     1091
     1092}
    9521093
    9531094
    … …  
    9961137}
    9971138
     1139###########################
     1140# Electron efficiency CHS #
     1141###########################
     1142
     1143module Efficiency ElectronEfficiencyCHS {
     1144
     1145  set InputArray ElectronIsolationCHS/electrons
     1146  set OutputArray electrons
     1147
     1148  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
     1149  # efficiency formula for electrons
     1150  set EfficiencyFormula {
     1151                                      (pt <= 4.0)  * (0.00) + \
     1152                         (abs(eta) <= 1.45 ) * (pt >  4.0 && pt <= 6.0)   * (0.50) + \
     1153                         (abs(eta) <= 1.45 ) * (pt >  6.0 && pt <= 8.0)   * (0.70) + \
     1154                         (abs(eta) <= 1.45 ) * (pt >  8.0 && pt <= 10.0)  * (0.85) + \
     1155                         (abs(eta) <= 1.45 ) * (pt > 10.0 && pt <= 30.0)  * (0.94) + \
     1156                         (abs(eta) <= 1.45 ) * (pt > 30.0 && pt <= 50.0)  * (0.97) + \
     1157                         (abs(eta) <= 1.45 ) * (pt > 50.0 && pt <= 70.0)  * (0.98) + \
     1158                         (abs(eta) <= 1.45 ) * (pt > 70.0 )  * (1.0) + \
     1159                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt >  4.0 && pt <= 10.0)   * (0.35) + \
     1160                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 10.0 && pt <= 30.0)   * (0.40) + \
     1161                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 30.0 && pt <= 70.0)   * (0.45) + \
     1162                         (abs(eta) > 1.45  && abs(eta) <= 1.55) * (pt > 70.0 )  * (0.55) + \
     1163                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt >  4.0 && pt <= 10.0)  * (0.75) + \
     1164                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 10.0 && pt <= 30.0)  * (0.85) + \
     1165                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 30.0 && pt <= 50.0)  * (0.95) + \
     1166                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 50.0 && pt <= 70.0)  * (0.95) + \
     1167                         (abs(eta) >= 1.55 && abs(eta) <= 2.0 ) * (pt > 70.0 )  * (1.0) + \
     1168                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt >  4.0 && pt <= 10.0)  * (0.65) + \
     1169                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 10.0 && pt <= 30.0)  * (0.75) + \
     1170                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 30.0 && pt <= 50.0)  * (0.90) + \
     1171                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 50.0 && pt <= 70.0)  * (0.90) + \
     1172                         (abs(eta) >= 2.0 && abs(eta) <= 2.5 ) * (pt > 70.0 )  * (0.90) + \
     1173                         (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 4.0 && pt <= 10.0) * (0.65) + \
     1174                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 10.0 && pt <= 30.0) * (0.75) + \
     1175                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 30.0 && pt <= 50.0) * (0.90) + \
     1176                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 50.0 && pt <= 70.0) * (0.90) + \
     1177                                          (abs(eta) > 2.5 && abs(eta) <= 4.0 ) * (pt > 70.0 ) * (0.90) + \
     1178                                          (abs(eta) > 4.0) * (0.00)
     1179
     1180  }
     1181}
     1182
     1183
    9981184##################
    9991185# Muon isolation #
    … …  
    10141200}
    10151201
    1016 
    1017 ##################
    1018 # Muon Loose Id  #
    1019 ##################
     1202######################
     1203# Muon isolation CHS #
     1204######################
     1205
     1206module Isolation MuonIsolationCHS {
     1207  set CandidateInputArray MuonMomentumSmearing/muons
     1208
     1209  # isolation collection
     1210  set IsolationInputArray EFlowMerger/eflow
     1211
     1212  set OutputArray muons
     1213
     1214  set DeltaRMax 0.3
     1215  set PTMin 1.0
     1216  set PTRatioMax 9999.
     1217
     1218}
     1219
     1220
     1221#####################
     1222# Muon Loose Id     #
     1223#####################
    10201224
    10211225module Efficiency MuonLooseIdEfficiency {
    … …  
    10391243}
    10401244
     1245
     1246#####################
     1247# Muon Loose Id CHS #
     1248#####################
     1249
     1250module Efficiency MuonLooseIdEfficiencyCHS {
     1251    set InputArray MuonIsolationCHS/muons
     1252    set OutputArray muons
     1253    # tracking + LooseID efficiency formula for muons
     1254    source muonLooseId.tcl
     1255
     1256}
     1257
     1258
     1259######################
     1260# Muon Tight Id  CHS #
     1261######################
     1262
     1263module Efficiency MuonTightIdEfficiencyCHS {
     1264    set InputArray MuonIsolationCHS/muons
     1265    set OutputArray muons
     1266    # tracking + TightID efficiency formula for muons
     1267    source muonTightId.tcl
     1268}
    10411269
    10421270
    … …  
    12511479module StatusPidFilter GenParticleFilter {
    12521480
    1253     set InputArray  Delphes/allParticles
     1481    set InputArray Delphes/allParticles
    12541482    set OutputArray filteredParticles
    12551483    set PTMin 5.0
    … …  
    12791507  add Branch MuonLooseIdEfficiency/muons MuonLoose Muon
    12801508  add Branch MuonTightIdEfficiency/muons MuonTight Muon
     1509
     1510  add Branch PhotonEfficiencyCHS/photons PhotonCHS Photon
     1511  add Branch ElectronEfficiencyCHS/electrons ElectronCHS Electron
     1512  add Branch MuonLooseIdEfficiencyCHS/muons MuonLooseCHS Muon
     1513  add Branch MuonTightIdEfficiencyCHS/muons MuonTightCHS Muon
    12811514
    12821515  add Branch JetEnergyScale/jets Jet Jet
    … …  
    12901523  add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
    12911524  add Branch ScalarHT/energy ScalarHT ScalarHT
    1292 }
     1525
     1526}

  • cards/converter_card.tcl

    r7bb13cd r70b9632  
    1313module TreeWriter TreeWriter {
    1414# add Branch InputArray BranchName BranchClass
    15   add Branch Delphes/allParticles Particle GenParticle
     15  add Branch Delphes/stableParticles Particle GenParticle
    1616}
    1717

  • cards/delphes_card_CMS.tcl

    r7bb13cd r70b9632  
    3535 
    3636  FastJetFinder
     37  FastCaloJetFinder
    3738
    3839  JetEnergyScale
    … …  
    210211  set HCalEnergyMin 1.0
    211212
    212   set ECalEnergySignificanceMin 1.0
    213   set HCalEnergySignificanceMin 1.0
     213  set ECalEnergySignificanceMin 2.0
     214  set HCalEnergySignificanceMin 2.0
    214215
    215216  set SmearTowerCenter true
    … …  
    497498}
    498499
     500################
     501# caloJet finder
     502################
     503
     504module FastJetFinder FastCaloJetFinder {
     505  set InputArray Calorimeter/towers
     506
     507  set OutputArray jets
     508
     509  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
     510  set JetAlgorithm 6
     511  set ParameterR 0.5
     512
     513  set JetPTMin 20.0
     514}
     515
    499516##################
    500517# Jet Energy Scale
    … …  
    610627
    611628  add Branch UniqueObjectFinder/jets Jet Jet
     629  add Branch FastCaloJetFinder/jets CaloJet Jet
    612630  add Branch UniqueObjectFinder/electrons Electron Electron
    613631  add Branch UniqueObjectFinder/photons Photon Photon

  • classes/DelphesClasses.cc

    r7bb13cd r70b9632  
    4141CompBase *Tower::fgCompare = CompE<Tower>::Instance();
    4242CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
     43CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
    4344CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
    4445
    … …  
    121122  Charge(0), Mass(0.0),
    122123  IsPU(0), IsRecoPU(0), IsConstituent(0), IsFromConversion(0),
     124  ClusterIndex(-1), ClusterNDF(0), ClusterSigma(0), SumPT2(0), BTVSumPT2(0), GenDeltaZ(0), GenSumPT2(0),
    123125  Flavor(0), FlavorAlgo(0), FlavorPhys(0),
    124126  BTag(0), BTagAlgo(0), BTagPhys(0),
    … …  
    127129  Momentum(0.0, 0.0, 0.0, 0.0),
    128130  Position(0.0, 0.0, 0.0, 0.0),
     131  PositionError(0.0, 0.0, 0.0, 0.0),
     132  InitialPosition(0.0, 0.0, 0.0, 0.0),
    129133  Area(0.0, 0.0, 0.0, 0.0),
    130   Dxy(0), SDxy(0), Xd(0), Yd(0), Zd(0),
     134  L(0),
     135  D0(0), ErrorD0(0),
     136  DZ(0), ErrorDZ(0),
     137  P(0),  ErrorP(0),
     138  PT(0), ErrorPT(0),
     139  CtgTheta(0), ErrorCtgTheta(0),
     140  Phi(0), ErrorPhi(0), 
     141  Xd(0), Yd(0), Zd(0),
    131142  TrackResolution(0),
    132143  NCharged(0),
    … …  
    245256  object.IsConstituent = IsConstituent;
    246257  object.IsFromConversion = IsFromConversion;
     258  object.ClusterIndex = ClusterIndex;
     259  object.ClusterNDF = ClusterNDF;
     260  object.ClusterSigma = ClusterSigma;
     261  object.SumPT2 = SumPT2;
     262  object.BTVSumPT2 = BTVSumPT2;
     263  object.GenDeltaZ = GenDeltaZ;
     264  object.GenSumPT2 = GenSumPT2;
    247265  object.Flavor = Flavor;
    248266  object.FlavorAlgo = FlavorAlgo;
    … …  
    262280  object.Momentum = Momentum;
    263281  object.Position = Position;
     282  object.InitialPosition = InitialPosition;
     283  object.PositionError = PositionError;
    264284  object.Area = Area;
    265   object.Dxy = Dxy;
    266   object.SDxy = SDxy;
     285  object.L = L;
     286  object.ErrorT = ErrorT;
     287  object.D0 = D0;
     288  object.ErrorD0 = ErrorD0;
     289  object.DZ = DZ;
     290  object.ErrorDZ = ErrorDZ;
     291  object.P = P;
     292  object.ErrorP = ErrorP;
     293  object.PT = PT;
     294  object.ErrorPT = ErrorPT;
     295  object.CtgTheta = CtgTheta ;
     296  object.ErrorCtgTheta = ErrorCtgTheta;
     297  object.Phi = Phi;
     298  object.ErrorPhi = ErrorPhi; 
    267299  object.Xd = Xd;
    268300  object.Yd = Yd;
    … …  
    282314  object.SumPtChargedPU = SumPtChargedPU;
    283315  object.SumPt = SumPt;
    284 
     316  object.ClusterIndex = ClusterIndex;
     317  object.ClusterNDF = ClusterNDF;
     318  object.ClusterSigma = ClusterSigma;
     319  object.SumPT2 = SumPT2;
     320 
    285321  object.FracPt[0] = FracPt[0];
    286322  object.FracPt[1] = FracPt[1];
    … …  
    363399  Momentum.SetXYZT(0.0, 0.0, 0.0, 0.0);
    364400  Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
     401  InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    365402  Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
    366   Dxy = 0.0;
    367   SDxy = 0.0;
     403  L = 0.0;
     404  ErrorT = 0.0;
     405  D0 = 0.0; 
     406  ErrorD0 = 0.0;
     407  DZ = 0.0;
     408  ErrorDZ = 0.0;
     409  P =0.0;
     410  ErrorP =0.0;
     411  PT = 0.0;
     412  ErrorPT = 0.0;
     413  CtgTheta = 0.0;
     414  ErrorCtgTheta = 0.0;
     415  Phi = 0.0;
     416  ErrorPhi = 0.0;
    368417  Xd = 0.0;
    369418  Yd = 0.0;
    … …  
    387436  SumPt = -999;
    388437
     438  ClusterIndex = -1;
     439  ClusterNDF = -99;
     440  ClusterSigma = 0.0;
     441  SumPT2 = 0.0;
     442  BTVSumPT2 = 0.0;
     443  GenDeltaZ = 0.0;
     444  GenSumPT2 = 0.0;
     445 
    389446  FracPt[0] = 0.0;
    390447  FracPt[1] = 0.0;

  • classes/DelphesClasses.h

    r7bb13cd r70b9632  
    147147  Float_t Pz; // particle momentum vector (z component) | hepevt.phep[number][2]
    148148
    149   Float_t PT; // particle transverse momentum
     149  Float_t D0;
     150  Float_t DZ;
     151  Float_t P;
     152  Float_t PT;
     153  Float_t CtgTheta;
     154  Float_t Phi;
    150155  Float_t Eta; // particle pseudorapidity
    151   Float_t Phi; // particle azimuthal angle
    152 
    153156  Float_t Rapidity; // particle rapidity
    154157
    … …  
    168171//---------------------------------------------------------------------------
    169172
    170 class Vertex: public TObject
     173class Vertex: public SortableObject
    171174{
    172175public:
    … …  
    176179  Float_t Z; // vertex position (z component)
    177180
    178   ClassDef(Vertex, 1)
     181  Double_t ErrorX;
     182  Double_t ErrorY;
     183  Double_t ErrorZ;
     184  Double_t ErrorT;
     185
     186  Int_t Index;
     187  Int_t NDF;
     188  Double_t Sigma;
     189  Double_t SumPT2;
     190  Double_t BTVSumPT2;
     191  Double_t GenDeltaZ;
     192  Double_t GenSumPT2;
     193
     194  TRefArray Constituents; // references to constituents
     195
     196  static CompBase *fgCompare; //!
     197  const CompBase *GetCompare() const { return fgCompare; }
     198
     199  ClassDef(Vertex, 3)
    179200};
    180201
    … …  
    397418  Int_t Charge; // track charge
    398419
    399   Float_t PT; // track transverse momentum
    400 
    401420  Float_t Eta; // track pseudorapidity
    402   Float_t Phi; // track azimuthal angle
    403421
    404422  Float_t EtaOuter; // track pseudorapidity at the tracker edge
    … …  
    415433  Float_t TOuter; // track position (z component) at the tracker edge
    416434
    417   Float_t Dxy;     // track signed transverse impact parameter
    418   Float_t SDxy;    // signed error on the track signed transverse impact parameter
     435  Float_t L; // track path length
     436  Float_t ErrorT; // error on the time measurement
     437
     438  Float_t D0;     // track signed transverse impact parameter
     439  Float_t ErrorD0;    // signed error on the track signed transverse impact parameter
     440
     441  Float_t DZ; // track transverse momentum
     442  Float_t ErrorDZ; // track transverse momentum error
     443
     444  Float_t P; // track transverse momentum
     445  Float_t ErrorP; // track transverse momentum error
     446
     447  Float_t PT; // track transverse momentum
     448  Float_t ErrorPT; // track transverse momentum error
     449
     450  Float_t CtgTheta; // track transverse momentum
     451  Float_t ErrorCtgTheta; // track transverse momentum error
     452
     453  Float_t Phi; // track azimuthal angle
     454  Float_t ErrorPhi; // track azimuthal angle
     455
    419456  Float_t Xd;      // X coordinate of point of closest approach to vertex
    420457  Float_t Yd;      // Y coordinate of point of closest approach to vertex
    … …  
    423460  TRef Particle; // reference to generated particle
    424461
     462  Int_t VertexIndex; // reference to vertex
     463
    425464  static CompBase *fgCompare; //!
    426465  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    526565  Float_t DeltaPhi;
    527566
    528   TLorentzVector Momentum, Position, Area;
    529 
    530   Float_t Dxy;
    531   Float_t SDxy;
     567  TLorentzVector Momentum, Position, InitialPosition, PositionError, Area;
     568
     569  Float_t L; // path length
     570  Float_t ErrorT; // path length
     571  Float_t D0;
     572  Float_t ErrorD0;
     573  Float_t DZ;
     574  Float_t ErrorDZ;
     575  Float_t P;
     576  Float_t ErrorP;
     577  Float_t PT;
     578  Float_t ErrorPT;
     579  Float_t CtgTheta;
     580  Float_t ErrorCtgTheta;
     581  Float_t Phi;
     582  Float_t ErrorPhi;
     583
    532584  Float_t Xd;
    533585  Float_t Yd;
    … …  
    535587
    536588  // tracking resolution
    537  
     589
    538590  Float_t TrackResolution;
    539591
    … …  
    562614  Float_t SumPt;
    563615
     616  // vertex variables
     617
     618  Int_t ClusterIndex;
     619  Int_t ClusterNDF;
     620  Double_t ClusterSigma;
     621  Double_t SumPT2;
     622  Double_t BTVSumPT2;
     623  Double_t GenDeltaZ;
     624  Double_t GenSumPT2;
     625
    564626  // N-subjettiness variables
    565627
    … …  
    595657  void SetFactory(DelphesFactory *factory) { fFactory = factory; }
    596658
    597   ClassDef(Candidate, 4)
     659  ClassDef(Candidate, 5)
    598660};
    599661

  • classes/SortableObject.h

    r7bb13cd r70b9632  
    156156      return -1;
    157157    else if(t1->ET < t2->ET)
     158      return 1;
     159    else
     160      return 0;
     161  }
     162};
     163
     164//---------------------------------------------------------------------------
     165
     166template <typename T>
     167class CompSumPT2: public CompBase
     168{
     169  CompSumPT2() {}
     170public:
     171  static CompSumPT2 *Instance()
     172  {
     173    static CompSumPT2 single;
     174    return &single;
     175  }
     176
     177  Int_t Compare(const TObject *obj1, const TObject *obj2) const
     178  {
     179    const T *t1 = static_cast<const T*>(obj1);
     180    const T *t2 = static_cast<const T*>(obj2);
     181    if(t1->SumPT2 > t2->SumPT2)
     182      return -1;
     183    else if(t1->SumPT2 < t2->SumPT2)
    158184      return 1;
    159185    else

  • doc/genMakefile.tcl

    r7bb13cd r70b9632  
    263263executableDeps {converters/*.cpp} {examples/*.cpp}
    264264
    265 executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp}
     265executableDeps {readers/DelphesHepMC.cpp} {readers/DelphesLHEF.cpp} {readers/DelphesSTDHEP.cpp} {readers/DelphesROOT.cpp}
    266266
    267267puts {ifeq ($(HAS_CMSSW),true)}

  • examples/Pythia8/configParticleGun.cmnd

    r7bb13cd r70b9632  
    33! 1) Settings used in the main program.
    44
    5 Main:numberOfEvents = 10000        ! number of events to generate
     5Main:numberOfEvents = 100000         ! number of events to generate
    66Main:timesAllowErrors = 3          ! how many aborts before run stops
    7 Main:spareFlag1 = on               ! true means particle gun
    8 Main:spareMode1 = ID               ! 1-5 - di-quark, 21 - di-gluon, 11 - single electron, 13 - single muon, 15 - single tau, 22 - single photon
    9 Main:spareParm1 = 10000            ! max pt
    10 Main:spareParm2 = 2.5              ! max eta
     7Main:spareFlag1 = on                ! true means particle gun
     8Main:spareMode1 = ID               ! 1-5 - di-quark, 21 - di-gluon, 11 - single electron, 13 - single muon, 22 - single photon
     9Main:spareParm1 = 10000           ! max pt
    1110
    1211! 2) Settings related to output in init(), next() and stat().

  • examples/Validation.cpp

    r7bb13cd r70b9632  
    2121#include <utility>
    2222#include <vector>
     23#include <typeinfo>
    2324
    2425#include "TROOT.h"
    … …  
    2829#include "TString.h"
    2930
     31#include "TH1.h"
    3032#include "TH2.h"
     33#include "TMath.h"
     34#include "TStyle.h"
     35#include "TGraph.h"
     36#include "TCanvas.h"
    3137#include "THStack.h"
    3238#include "TLegend.h"
    … …  
    3440#include "TClonesArray.h"
    3541#include "TLorentzVector.h"
     42#include "TGraphErrors.h"
     43#include "TMultiGraph.h"
    3644
    3745#include "classes/DelphesClasses.h"
    … …  
    4755//------------------------------------------------------------------------------
    4856
    49 // Here you can put your analysis macro
    50 
    51 #include "Validation.C"
     57double ptrangemin = 10;
     58double ptrangemax = 10000;
     59static const int Nbins = 20;
     60
     61int objStyle = 1;
     62int trackStyle = 7;
     63int towerStyle = 3;
     64
     65Color_t objColor = kBlack;
     66Color_t trackColor = kBlack;
     67Color_t towerColor = kBlack;
     68
     69double effLegXmin = 0.22;
     70double effLegXmax = 0.7;
     71double effLegYmin = 0.22;
     72double effLegYmax = 0.5;
     73
     74double resLegXmin = 0.62;
     75double resLegXmax = 0.9;
     76double resLegYmin = 0.52;
     77double resLegYmax = 0.85;
     78
     79double topLeftLegXmin = 0.22;
     80double topLeftLegXmax = 0.7;
     81double topLeftLegYmin = 0.52;
     82double topLeftLegYmax = 0.85;
     83
     84
     85struct resolPlot
     86{
     87  TH1 *cenResolHist;
     88  TH1 *fwdResolHist;
     89  double ptmin;
     90  double ptmax;
     91  double xmin;
     92  double xmax;
     93  TString obj;
     94
     95  resolPlot();
     96  resolPlot(double ptdown, double ptup, TString object);
     97  void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
     98  void print(){std::cout << ptmin << std::endl;}
     99};
     100
     101
     102resolPlot::resolPlot()
     103{
     104}
     105
     106resolPlot::resolPlot(double ptdown, double ptup, TString object)
     107{
     108  this->set(ptdown,ptup,object);
     109}
     110
     111void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
     112{
     113  ptmin = ptdown;
     114  ptmax = ptup;
     115  obj = object;
     116
     117  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);
     118  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);
     119}
     120
     121void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
     122{
     123  double width;
     124  double ptdown;
     125  double ptup;
     126  resolPlot ptemp;
     127
     128  for (int i = 0; i < Nbins; i++)
     129  {
     130    width = (ptmax - ptmin) / Nbins;
     131    ptdown = TMath::Power(10,ptmin + i * width );
     132    ptup = TMath::Power(10,ptmin + (i+1) * width );
     133    ptemp.set(ptdown, ptup, obj, xmin, xmax);
     134    histos->push_back(ptemp);
     135  }
     136}
     137
     138//------------------------------------------------------------------------------
     139
     140class ExRootResult;
     141class ExRootTreeReader;
     142
     143//------------------------------------------------------------------------------
     144
     145void BinLogX(TH1*h)
     146{
     147  TAxis *axis = h->GetXaxis();
     148  int bins = axis->GetNbins();
     149
     150  Axis_t from = axis->GetXmin();
     151  Axis_t to = axis->GetXmax();
     152  Axis_t width = (to - from) / bins;
     153  Axis_t *new_bins = new Axis_t[bins + 1];
     154
     155  for (int i = 0; i <= bins; i++)
     156  {
     157    new_bins[i] = TMath::Power(10, from + i * width);
     158  }
     159  axis->Set(bins, new_bins);
     160  delete new_bins;
     161}
     162
     163
     164//------------------------------------------------------------------------------
     165
     166template<typename T>
     167std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
     168{
     169
     170  cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
     171
     172  Long64_t allEntries = treeReader->GetEntries();
     173
     174  GenParticle *particle;
     175  T *recoObj;
     176
     177  TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
     178
     179  Float_t deltaR;
     180  Float_t pt, eta;
     181  Long64_t entry;
     182
     183  Int_t i, j;
     184
     185  TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
     186  TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
     187  TH1D *histGenPtfwd  = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
     188  TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
     189
     190
     191  BinLogX(histGenPtcen);
     192  BinLogX(histRecoPtcen);
     193  BinLogX(histGenPtfwd);
     194  BinLogX(histRecoPtfwd);
     195
     196  // Loop over all events
     197  for(entry = 0; entry < allEntries; ++entry)
     198  {
     199    // Load selected branches with data from specified event
     200    treeReader->ReadEntry(entry);
     201
     202    // Loop over all generated particle in event
     203    for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
     204    {
     205
     206      particle = (GenParticle*) branchParticle->At(i);
     207      genMomentum = particle->P4();
     208
     209      deltaR = 999;
     210
     211      if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
     212      {
     213
     214        // Loop over all reco object in event
     215        for(j = 0; j < branchReco->GetEntriesFast(); ++j)
     216        {
     217          recoObj = (T*)branchReco->At(j);
     218          recoMomentum = recoObj->P4();
     219          // this is simply to avoid warnings from initial state particle
     220          // having infite rapidity ...
     221        //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
     222
     223          // take the closest parton candidate
     224          if(TMath::Abs(pdgID) == 5)
     225          {
     226            Jet *jet = (Jet *)recoObj;
     227            if(jet->BTag != 1) continue;
     228          }
     229          if(TMath::Abs(pdgID) == 15)
     230          {
     231            Jet *jet = (Jet *)recoObj;
     232            if(jet->TauTag != 1) continue;
     233          }
     234          if(genMomentum.DeltaR(recoMomentum) < deltaR)
     235          {
     236            deltaR = genMomentum.DeltaR(recoMomentum);
     237            bestRecoMomentum = recoMomentum;
     238          }
     239        }
     240
     241        pt  = genMomentum.Pt();
     242        eta = genMomentum.Eta();
     243
     244        if (TMath::Abs(eta) < 1.5)
     245        {
     246          histGenPtcen->Fill(pt);
     247          if(deltaR < 0.3) { histRecoPtcen->Fill(pt); }
     248        }
     249        else if (TMath::Abs(eta) < 2.5)
     250        {
     251          histGenPtfwd->Fill(pt);
     252          if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); }
     253
     254        }
     255      }
     256    }
     257  }
     258
     259
     260  std::pair<TH1D*,TH1D*> histos;
     261
     262  histRecoPtcen->Divide(histGenPtcen);
     263  histRecoPtfwd->Divide(histGenPtfwd);
     264
     265  histos.first = histRecoPtcen;
     266  histos.second = histRecoPtfwd;
     267
     268  return histos;
     269}
     270
     271template<typename T>
     272void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
     273{
     274  Long64_t allEntries = treeReader->GetEntries();
     275
     276  cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
     277
     278  GenParticle *particle;
     279  T* recoObj;
     280
     281  TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
     282
     283  Float_t deltaR;
     284  Float_t pt, eta;
     285  Long64_t entry;
     286
     287  Int_t i, j, bin;
     288
     289  // Loop over all events
     290  for(entry = 0; entry < allEntries; ++entry)
     291  {
     292    // Load selected branches with data from specified event
     293    treeReader->ReadEntry(entry);
     294
     295    // Loop over all reconstructed jets in event
     296    for(i = 0; i < branchReco->GetEntriesFast(); ++i)
     297    {
     298      recoObj = (T*) branchReco->At(i);
     299      recoMomentum = recoObj->P4();
     300
     301      deltaR = 999;
     302
     303     // Loop over all hard partons in event
     304     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
     305     {
     306        particle = (GenParticle*) branchParticle->At(j);
     307        if (particle->PID == pdgID && particle->Status == 1)
     308        {
     309          genMomentum = particle->P4();
     310
     311          // this is simply to avoid warnings from initial state particle
     312          // having infite rapidity ...
     313          if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
     314
     315          // take the closest parton candidate
     316          if(genMomentum.DeltaR(recoMomentum) < deltaR)
     317          {
     318             deltaR = genMomentum.DeltaR(recoMomentum);
     319             bestGenMomentum = genMomentum;
     320          }
     321        }
     322      }
     323
     324      if(deltaR < 0.3)
     325      {
     326        pt  = bestGenMomentum.Pt();
     327        eta = TMath::Abs(bestGenMomentum.Eta());
     328
     329        for (bin = 0; bin < Nbins; bin++)
     330        {
     331          if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
     332          {
     333            if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
     334            else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
     335          }
     336        }
     337      }
     338    }
     339  }
     340}
     341void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
     342{
     343
     344  Long64_t allEntries = treeReader->GetEntries();
     345
     346  cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
     347
     348  Jet *jet, *genjet;
     349
     350  TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
     351
     352  Float_t deltaR;
     353  Float_t pt, eta;
     354  Long64_t entry;
     355
     356  Int_t i, j, bin;
     357
     358  // Loop over all events
     359  for(entry = 0; entry < allEntries; ++entry)
     360  {
     361    // Load selected branches with data from specified event
     362    treeReader->ReadEntry(entry);
     363
     364    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
     365
     366    // Loop over all reconstructed jets in event
     367    for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
     368    {
     369
     370      jet = (Jet*) branchJet->At(i);
     371      jetMomentum = jet->P4();
     372
     373      deltaR = 999;
     374
     375     // Loop over all hard partons in event
     376     for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
     377     {
     378        genjet = (Jet*) branchGenJet->At(j);
     379
     380        genJetMomentum = genjet->P4();
     381
     382        // this is simply to avoid warnings from initial state particle
     383        // having infite rapidity ...
     384        if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
     385
     386        // take the closest parton candidate
     387        if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
     388        {
     389           deltaR = genJetMomentum.DeltaR(jetMomentum);
     390           bestGenJetMomentum = genJetMomentum;
     391        }
     392      }
     393
     394      if(deltaR < 0.25)
     395      {
     396        pt  = genJetMomentum.Pt();
     397        eta = TMath::Abs(genJetMomentum.Eta());
     398
     399        for (bin = 0; bin < Nbins; bin++)
     400        {
     401            if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
     402            {
     403                histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
     404            }
     405            else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
     406            {
     407                histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
     408            }
     409
     410        }
     411      }
     412    }
     413  }
     414}
     415
     416void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
     417{
     418
     419  Long64_t allEntries = treeReader->GetEntries();
     420
     421  cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
     422
     423  MissingET *met;
     424  ScalarHT *scalarHT;
     425
     426  Long64_t entry;
     427
     428  Int_t bin;
     429  Double_t ht;
     430
     431  Jet *jet;
     432  TLorentzVector p1, p2;
     433
     434  // Loop over all events
     435  for(entry = 0; entry < allEntries; ++entry)
     436  {
     437    // Load selected branches with data from specified event
     438    treeReader->ReadEntry(entry);
     439
     440    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
     441
     442    if (branchJet->GetEntriesFast() > 1)
     443    {
     444
     445      jet = (Jet*) branchJet->At(0);
     446      p1 = jet->P4();
     447      jet = (Jet*) branchJet->At(1);
     448      p2 = jet->P4();
     449
     450      met = (MissingET*) branchMet->At(0);
     451      scalarHT = (ScalarHT*) branchScalarHT->At(0);
     452      ht = scalarHT->HT;
     453
     454      if(p1.Pt() < 0.75*ht/2) continue;
     455      if(p2.Pt() < 0.75*ht/2) continue;
     456
     457      for (bin = 0; bin < Nbins; bin++)
     458      {
     459        if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
     460        {
     461          histos->at(bin).cenResolHist->Fill(met->P4().Px());
     462          histos->at(bin).fwdResolHist->Fill(met->P4().Py());
     463        }
     464      }
     465    }
     466  }
     467}
     468
     469
     470std::pair<Double_t, Double_t> GausFit(TH1* hist)
     471{
     472  TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
     473  hist->Fit("f1","RQ");
     474  Double_t sig = f1->GetParameter(2);
     475  Double_t sigErr = f1->GetParError(2);
     476  delete f1;
     477  return make_pair (sig, sigErr);
     478}
     479
     480
     481TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
     482{
     483  Int_t bin;
     484  Int_t count = 0;
     485  TGraphErrors gr = TGraphErrors(Nbins/2);
     486  Double_t sig = 0;
     487  Double_t sigErr = 0;
     488  for (bin = 0; bin < Nbins; bin++)
     489  {
     490    if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
     491    {
     492      std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
     493      if (rms == true)
     494      {
     495        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
     496        gr.SetPointError(count,0, sigvalues.second); // to correct
     497      }
     498      else
     499      {
     500        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
     501        gr.SetPointError(count,0, sigvalues.second);
     502      }
     503      count++;
     504    }
     505
     506    else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10)
     507    {
     508      std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist);
     509      if (rms == true)
     510      {
     511        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
     512        gr.SetPointError(count,0, sigvalues.second); // to correct
     513      }
     514      else
     515      {
     516        gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
     517        gr.SetPointError(count,0, sigvalues.second);
     518      }
     519      count++;
     520    }
     521
     522  }
     523  return gr;
     524}
     525
     526
     527//------------------------------------------------------------------------------
     528
     529
     530// type 1 : object, 2 : track, 3 : tower
     531
     532void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
     533{
     534
     535  gr->SetLineWidth(2);
     536
     537  switch ( type )
     538  {
     539    case 1:
     540      gr->SetLineColor(objColor);
     541      gr->SetLineStyle(objStyle);
     542      std::cout << "Adding " << gr->GetName() << std::endl;
     543      mg->Add(gr);
     544      leg->AddEntry(gr,"Reco","l");
     545      break;
     546
     547    case 2:
     548      gr->SetLineColor(trackColor);
     549      gr->SetLineStyle(trackStyle);
     550      mg->Add(gr);
     551      leg->AddEntry(gr,"Track","l");
     552      break;
     553
     554    case 3:
     555      gr->SetLineColor(towerColor);
     556      gr->SetLineStyle(towerStyle);
     557      mg->Add(gr);
     558      leg->AddEntry(gr,"Tower","l");
     559      break;
     560
     561    case 0:
     562      gr->SetLineColor(objColor);
     563      gr->SetLineStyle(objStyle);
     564      mg->Add(gr);
     565      break;
     566
     567    default:
     568      std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
     569      break;
     570  }
     571}
     572
     573void addHist(TH1D *h, TLegend *leg, int type)
     574{
     575  h->SetLineWidth(2);
     576
     577  switch ( type )
     578  {
     579    case 1:
     580      h->SetLineColor(objColor);
     581      h->SetLineStyle(objStyle);
     582      leg->AddEntry(h,"Reco","l");
     583      break;
     584
     585    case 2:
     586      h->SetLineColor(trackColor);
     587      h->SetLineStyle(trackStyle);
     588      leg->AddEntry(h,"Track","l");
     589      break;
     590
     591    case 3:
     592      h->SetLineColor(towerColor);
     593      h->SetLineStyle(towerStyle);
     594      leg->AddEntry(h,"Tower","l");
     595      break;
     596
     597    case 0:
     598      h->SetLineColor(objColor);
     599      h->SetLineStyle(objStyle);
     600      break;
     601
     602    default:
     603      std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
     604      break;
     605  }
     606}
     607
     608void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
     609{
     610  mg->SetMinimum(0.);
     611  mg->SetMaximum(max);
     612  mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
     613  mg->GetYaxis()->SetTitle("resolution");
     614  mg->GetXaxis()->SetTitle("p_{T} [GeV]");
     615  mg->GetYaxis()->SetTitleSize(0.07);
     616  mg->GetXaxis()->SetTitleSize(0.07);
     617  mg->GetYaxis()->SetLabelSize(0.06);
     618  mg->GetXaxis()->SetLabelSize(0.06);
     619  mg->GetYaxis()->SetLabelOffset(0.03);
     620  mg->GetYaxis()->SetTitleOffset(1.4);
     621  mg->GetXaxis()->SetTitleOffset(1.4);
     622
     623  mg->GetYaxis()->SetNdivisions(505);
     624
     625  leg->SetBorderSize(0);
     626  leg->SetShadowColor(0);
     627  leg->SetFillColor(0);
     628  leg->SetFillStyle(0);
     629
     630  gStyle->SetOptTitle(0);
     631  gPad->SetLogx();
     632  gPad->SetBottomMargin(0.2);
     633  gPad->SetLeftMargin(0.2);
     634  gPad->Modified();
     635  gPad->Update();
     636
     637}
     638
     639void DrawAxis(TH1D *h, TLegend *leg, int type)
     640{
     641
     642  h->GetYaxis()->SetRangeUser(0,1.0);
     643  if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
     644  else h->GetXaxis()->SetTitle("#eta");
     645  h->GetYaxis()->SetTitle("efficiency");
     646  h->GetYaxis()->SetTitleSize(0.07);
     647  h->GetXaxis()->SetTitleSize(0.07);
     648  h->GetYaxis()->SetLabelSize(0.06);
     649  h->GetXaxis()->SetLabelSize(0.06);
     650  h->GetYaxis()->SetLabelOffset(0.03);
     651  h->GetYaxis()->SetTitleOffset(1.3);
     652  h->GetXaxis()->SetTitleOffset(1.4);
     653
     654  h->GetYaxis()->SetNdivisions(505);
     655
     656  leg->SetBorderSize(0);
     657  leg->SetShadowColor(0);
     658  leg->SetFillColor(0);
     659  leg->SetFillStyle(0);
     660
     661  gStyle->SetOptTitle(0);
     662  gStyle->SetOptStat(0);
     663  gPad->SetBottomMargin(0.2);
     664  gPad->SetLeftMargin(0.2);
     665
     666  gPad->Modified();
     667  gPad->Update();
     668
     669}
     670
     671
     672void Validation(const char *inputFileElectron, const char *inputFileMuon, const char *inputFilePhoton, const char *inputFileJet, const char *inputFileBJet, const char *inputFileTauJet, const char *outputFile)
     673{
     674  TChain *chainElectron = new TChain("Delphes");
     675  chainElectron->Add(inputFileElectron);
     676  ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
     677
     678  TChain *chainMuon = new TChain("Delphes");
     679  chainMuon->Add(inputFileMuon);
     680  ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
     681
     682  TChain *chainPhoton = new TChain("Delphes");
     683  chainPhoton->Add(inputFilePhoton);
     684  ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
     685
     686  TChain *chainJet = new TChain("Delphes");
     687  chainJet->Add(inputFileJet);
     688  ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
     689
     690  TChain *chainBJet = new TChain("Delphes");
     691  chainBJet->Add(inputFileBJet);
     692  ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
     693
     694  TChain *chainTauJet = new TChain("Delphes");
     695  chainTauJet->Add(inputFileTauJet);
     696  ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
     697
     698  TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
     699  TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
     700  TClonesArray *branchTowerElectron = treeReaderElectron->UseBranch("Tower");
     701  TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
     702
     703  TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
     704  TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
     705  TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
     706
     707  TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
     708  TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
     709  TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
     710
     711  TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
     712  TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet");
     713  TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
     714
     715  TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
     716  TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
     717
     718  TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
     719  TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
     720
     721  TClonesArray *branchScalarHT = treeReaderJet->UseBranch("ScalarHT");
     722  TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET");
     723
     724  ///////////////
     725  // Electrons //
     726  ///////////////
     727
     728  // Reconstruction efficiency
     729  TString elecs = "Electron";
     730  int elID = 11;
     731  std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticleElectron, "Electron", elID, treeReaderElectron);
     732
     733  // tracking reconstruction efficiency
     734  std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrackElectron, branchParticleElectron, "electronTrack", elID, treeReaderElectron);
     735
     736  // Tower reconstruction efficiency
     737  std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTowerElectron, branchParticleElectron, "electronTower", elID, treeReaderElectron);
     738
     739  // Electron Energy Resolution
     740  std::vector<resolPlot> plots_el;
     741  HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
     742  GetEres<Electron>(&plots_el, branchElectron, branchParticleElectron, elID, treeReaderElectron);
     743  TGraphErrors gr_el = EresGraph(&plots_el, true);
     744  TGraphErrors gr_elFwd = EresGraph(&plots_el, false);
     745  gr_el.SetName("Electron");
     746  gr_elFwd.SetName("ElectronFwd");
     747
     748  // Electron Track Energy Resolution
     749  std::vector<resolPlot> plots_eltrack;
     750  HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
     751  GetEres<Track>(&plots_eltrack, branchTrackElectron, branchParticleElectron, elID, treeReaderElectron);
     752  TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
     753  TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false);
     754  gr_eltrack.SetName("ElectronTracks");
     755  gr_eltrackFwd.SetName("ElectronTracksFwd");
     756
     757  // Electron Tower Energy Resolution
     758  std::vector<resolPlot> plots_eltower;
     759  HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
     760  GetEres<Tower>(&plots_eltower, branchTowerElectron, branchParticleElectron, elID, treeReaderElectron);
     761  TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
     762  TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false);
     763  gr_eltower.SetName("ElectronTower");
     764  gr_eltrackFwd.SetName("ElectronTracksFwd");
     765
     766  // Canvases
     767  TString elEff = "electronEff";
     768  TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600);
     769  C_el1->Divide(2);
     770  C_el1->cd(1);
     771  TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     772  leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
     773  leg_el1->AddEntry("","","");
     774
     775  gPad->SetLogx();
     776  histos_eltrack.first->Draw("][");
     777  addHist(histos_eltrack.first, leg_el1, 2);
     778  histos_el.first->Draw("same ][");
     779  addHist(histos_el.first, leg_el1, 1);
     780  DrawAxis(histos_eltrack.first, leg_el1, 0);
     781
     782  leg_el1->Draw();
     783
     784  C_el1->cd(2);
     785  TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     786  leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
     787  leg_el2->AddEntry("","","");
     788
     789  gPad->SetLogx();
     790  histos_eltrack.second->Draw("][");
     791  addHist(histos_eltrack.second, leg_el2, 2);
     792  histos_el.second->Draw("same ][");
     793  addHist(histos_el.second, leg_el2, 1);
     794
     795  DrawAxis(histos_eltrack.second, leg_el2, 0);
     796  leg_el2->Draw();
     797
     798  C_el1->cd(0);
     799
     800  TString elRes = "electronERes";
     801  TString elResFwd = "electronEResForward";
     802  TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600);
     803  C_el2->Divide(2);
     804  C_el2->cd(1);
     805  TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
     806  TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     807  leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
     808  leg_el->AddEntry("","","");
     809
     810  addGraph(mg_el, &gr_eltower, leg_el, 3);
     811  addGraph(mg_el, &gr_eltrack, leg_el, 2);
     812  addGraph(mg_el, &gr_el, leg_el, 1);
     813
     814  mg_el->Draw("ACX");
     815  leg_el->Draw();
     816
     817  DrawAxis(mg_el, leg_el, 0.1);
     818
     819  C_el2->cd(2);
     820  TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd);
     821  TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     822  leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
     823  leg_elFwd->AddEntry("","","");
     824
     825  addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3);
     826  addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2);
     827  addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1);
     828
     829  mg_elFwd->Draw("ACX");
     830  leg_elFwd->Draw();
     831
     832  DrawAxis(mg_elFwd, leg_elFwd, 0.2);
     833
     834  C_el1->Print("validation.pdf(","pdf");
     835  C_el2->Print("validation.pdf","pdf");
     836
     837  gDirectory->cd(0);
     838
     839  ///////////
     840  // Muons //
     841  ///////////
     842
     843  // Reconstruction efficiency
     844  int muID = 13;
     845  std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon);
     846
     847  // muon tracking reconstruction efficiency
     848  std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
     849
     850  // Muon Energy Resolution
     851  std::vector<resolPlot> plots_mu;
     852  HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
     853  GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
     854  TGraphErrors gr_mu = EresGraph(&plots_mu, true);
     855  TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
     856  gr_mu.SetName("Muon");
     857  gr_muFwd.SetName("MuonFwd");
     858
     859  // Muon Track Energy Resolution
     860  std::vector<resolPlot> plots_mutrack;
     861  HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
     862  GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
     863  TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
     864  TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
     865  gr_mutrackFwd.SetName("MuonTracksFwd");
     866
     867  // Canvas
     868
     869  TString muEff = "muonEff";
     870  TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
     871  C_mu1->Divide(2);
     872  C_mu1->cd(1);
     873  TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     874  leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
     875  leg_mu1->AddEntry("","","");
     876
     877
     878  gPad->SetLogx();
     879  histos_mutrack.first->Draw("][");
     880  addHist(histos_mutrack.first, leg_mu1, 2);
     881  histos_mu.first->Draw("same ][");
     882  addHist(histos_mu.first, leg_mu1, 1);
     883
     884  DrawAxis(histos_mutrack.first, leg_mu1, 0);
     885
     886  leg_mu1->Draw();
     887
     888  C_mu1->cd(2);
     889  TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     890  leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
     891  leg_mu2->AddEntry("","","");
     892
     893  gPad->SetLogx();
     894  histos_mutrack.second->Draw("][");
     895  addHist(histos_mutrack.second, leg_mu2, 2);
     896  histos_mu.second->Draw("same ][");
     897  addHist(histos_mu.second, leg_mu2, 1);
     898
     899  DrawAxis(histos_mutrack.second, leg_mu2, 1);
     900  leg_mu2->Draw();
     901
     902  TString muRes = "muonERes";
     903  TString muResFwd = "muonEResFwd";
     904
     905  TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
     906  C_mu->Divide(2);
     907  C_mu->cd(1);
     908  TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
     909  TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
     910  leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
     911  leg_mu->AddEntry("","","");
     912
     913  addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
     914  addGraph(mg_mu, &gr_mu, leg_mu, 1);
     915
     916  mg_mu->Draw("ACX");
     917  leg_mu->Draw();
     918
     919  DrawAxis(mg_mu, leg_mu, 0.3);
     920
     921  C_mu->cd(2);
     922  TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
     923  TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
     924  leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
     925  leg_muFwd->AddEntry("","","");
     926
     927  addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
     928  addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
     929
     930  mg_muFwd->Draw("ACX");
     931  leg_muFwd->Draw();
     932
     933  DrawAxis(mg_muFwd, leg_muFwd, 0.3);
     934
     935  //C_mu1->SaveAs(muEff+".eps");
     936  //C_mu->SaveAs(muRes+".eps");
     937
     938  C_mu1->Print("validation.pdf","pdf");
     939  C_mu->Print("validation.pdf","pdf");
     940
     941  gDirectory->cd(0);
     942
     943  /////////////
     944  // Photons //
     945  /////////////
     946
     947  // Reconstruction efficiency
     948  int phID = 22;
     949  std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
     950  std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
     951
     952  // Photon Energy Resolution
     953  std::vector<resolPlot> plots_ph;
     954  HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
     955  GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton);
     956  TGraphErrors gr_ph = EresGraph(&plots_ph, true);
     957  TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
     958  gr_ph.SetName("Photon");
     959  gr_phFwd.SetName("PhotonFwd");
     960
     961
     962  // Photon Tower Energy Resolution
     963  std::vector<resolPlot> plots_phtower;
     964  HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
     965  GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton);
     966  TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
     967  TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
     968  gr_phtower.SetName("PhotonTower");
     969  gr_phtowerFwd.SetName("PhotonTowerFwd");
     970
     971  // Canvas
     972
     973  TString phEff = "photonEff";
     974  TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
     975  C_ph1->Divide(2);
     976  C_ph1->cd(1);
     977  TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     978  leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
     979  leg_ph1->AddEntry("","","");
     980
     981
     982  gPad->SetLogx();
     983  histos_phtower.first->Draw("][");
     984  addHist(histos_phtower.first, leg_ph1, 3);
     985  histos_ph.first->Draw("same ][");
     986  addHist(histos_ph.first, leg_ph1, 1);
     987
     988  DrawAxis(histos_phtower.first, leg_ph1, 0);
     989  leg_ph1->Draw();
     990
     991  C_ph1->cd(2);
     992  TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
     993  leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
     994  leg_ph2->AddEntry("","","");
     995
     996
     997  gPad->SetLogx();
     998  histos_phtower.second->Draw("][");
     999  addHist(histos_phtower.second, leg_ph2, 3);
     1000  histos_ph.second->Draw("same ][");
     1001  addHist(histos_ph.second, leg_ph2, 1);
     1002
     1003  DrawAxis(histos_phtower.second, leg_ph2, 1);
     1004  leg_ph2->Draw();
     1005
     1006  C_ph1->SaveAs(phEff+".eps");
     1007
     1008  TString phRes = "phERes";
     1009  TString phResFwd = "phEResFwd";
     1010
     1011  TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
     1012  C_ph->Divide(2);
     1013  C_ph->cd(1);
     1014  TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
     1015  TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1016  leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
     1017  leg_ph->AddEntry("","","");
     1018
     1019  addGraph(mg_ph, &gr_phtower, leg_ph, 3);
     1020  addGraph(mg_ph, &gr_ph, leg_ph, 1);
     1021
     1022  mg_ph->Draw("ACX");
     1023  leg_ph->Draw();
     1024
     1025  DrawAxis(mg_ph, leg_ph, 0.3);
     1026
     1027  C_ph->cd(2);
     1028  TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
     1029  TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1030  leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
     1031  leg_phFwd->AddEntry("","","");
     1032
     1033  addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
     1034  addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
     1035
     1036  mg_phFwd->Draw("ACX");
     1037  leg_phFwd->Draw();
     1038
     1039  DrawAxis(mg_phFwd, leg_phFwd, 0.3);
     1040
     1041  C_ph->SaveAs(phRes+".eps");
     1042
     1043  C_ph1->Print("validation.pdf","pdf");
     1044  C_ph->Print("validation.pdf","pdf");
     1045
     1046  gDirectory->cd(0);
     1047
     1048  //////////
     1049  // Jets //
     1050  //////////
     1051
     1052  // BJets Reconstruction efficiency
     1053  int bID = 5;
     1054  std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet);
     1055
     1056  // TauJets Reconstruction efficiency
     1057  int tauID = 15;
     1058  std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet);
     1059
     1060  // PFJets Energy Resolution
     1061  std::vector<resolPlot> plots_pfjets;
     1062  HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
     1063  GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet);
     1064  TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
     1065  TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
     1066  gr_pfjets.SetName("pfJet");
     1067  gr_pfjetsFwd.SetName("pfJetFwd");
     1068
     1069  // CaloJets Energy Resolution
     1070  std::vector<resolPlot> plots_calojets;
     1071  HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
     1072  GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet);
     1073  TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
     1074  TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
     1075  gr_calojets.SetName("caloJet");
     1076  gr_calojetsFwd.SetName("caloJetFwd");
     1077
     1078  // MET Resolution vs HT
     1079  std::vector<resolPlot> plots_met;
     1080  HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
     1081  GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet);
     1082  TGraphErrors gr_met = EresGraph(&plots_met, true);
     1083  gr_calojets.SetName("MET");
     1084
     1085  // Canvas
     1086  TString btagEff = "btagEff";
     1087  TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
     1088  C_btag1->Divide(2);
     1089  C_btag1->cd(1);
     1090  TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1091  leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
     1092  leg_btag1->AddEntry("","","");
     1093
     1094  gPad->SetLogx();
     1095  histos_btag.first->Draw();
     1096  addHist(histos_btag.first, leg_btag1, 0);
     1097
     1098  DrawAxis(histos_btag.first, leg_btag1, 0);
     1099  leg_btag1->Draw();
     1100
     1101  C_btag1->cd(2);
     1102  TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
     1103  leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
     1104  leg_btag2->AddEntry("","","");
     1105
     1106  gPad->SetLogx();
     1107  histos_btag.second->Draw();
     1108  addHist(histos_btag.second, leg_btag2, 0);
     1109
     1110  DrawAxis(histos_btag.second, leg_btag2, 0);
     1111  leg_btag2->Draw();
     1112  C_btag1->cd(0);
     1113
     1114  TString tautagEff = "tautagEff";
     1115  TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
     1116  C_tautag1->Divide(2);
     1117  C_tautag1->cd(1);
     1118  TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1119  leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
     1120  leg_tautag1->AddEntry("","","");
     1121
     1122  gPad->SetLogx();
     1123  histos_tautag.first->Draw();
     1124  addHist(histos_tautag.first, leg_tautag1, 0);
     1125
     1126  DrawAxis(histos_tautag.first, leg_tautag1, 0);
     1127  leg_tautag1->Draw();
     1128
     1129  C_tautag1->cd(2);
     1130  TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
     1131  leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
     1132  leg_tautag2->AddEntry("","","");
     1133
     1134  gPad->SetLogx();
     1135  histos_tautag.second->Draw();
     1136  addHist(histos_tautag.second, leg_tautag2, 0);
     1137
     1138  DrawAxis(histos_tautag.second, leg_tautag2, 0);
     1139  leg_tautag2->Draw();
     1140  C_tautag1->cd(0);
     1141
     1142  TString jetRes = "jetERes";
     1143  TString jetResFwd = "jetEResFwd";
     1144  TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
     1145  C_jet->Divide(2);
     1146
     1147  C_jet->cd(1);
     1148  TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
     1149  TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1150  leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
     1151  leg_jet->AddEntry("","","");
     1152
     1153  addGraph(mg_jet, &gr_calojets, leg_jet, 3);
     1154  addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
     1155
     1156  mg_jet->Draw("ACX");
     1157  leg_jet->Draw();
     1158
     1159  DrawAxis(mg_jet, leg_jet, 0.25);
     1160
     1161  C_jet->cd(2);
     1162  TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
     1163  TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
     1164  leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
     1165  leg_jetFwd->AddEntry("","","");
     1166
     1167  addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
     1168  addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
     1169
     1170  mg_jetFwd->Draw("ACX");
     1171  leg_jetFwd->Draw();
     1172
     1173  DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
     1174
     1175  C_btag1->SaveAs(btagEff+".eps");
     1176  C_jet->SaveAs(jetRes+".eps");
     1177
     1178  TString metRes = "MetRes";
     1179  TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
     1180
     1181  TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
     1182  TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
     1183  leg_met->SetHeader("E_{T}^{miss}");
     1184  leg_met->AddEntry("","","");
     1185
     1186
     1187  addGraph(mg_met, &gr_met, leg_met, 0);
     1188
     1189  mg_met->Draw("ACX");
     1190  leg_met->Draw();
     1191
     1192  DrawAxis(mg_met, leg_met, 300);
     1193
     1194  mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
     1195  mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
     1196
     1197  C_met->SaveAs(metRes+".eps");
     1198
     1199  C_jet->Print("validation.pdf","pdf");
     1200  C_btag1->Print("validation.pdf","pdf");
     1201  C_tautag1->Print("validation.pdf","pdf");
     1202  C_met->Print("validation.pdf)","pdf");
     1203
     1204  TFile *fout = new TFile(outputFile,"recreate");
     1205
     1206  for (int bin = 0; bin < Nbins; bin++)
     1207  {
     1208    plots_el.at(bin).cenResolHist->Write();
     1209    plots_eltrack.at(bin).cenResolHist->Write();
     1210    plots_eltower.at(bin).cenResolHist->Write();
     1211    plots_el.at(bin).fwdResolHist->Write();
     1212    plots_eltrack.at(bin).fwdResolHist->Write();
     1213    plots_eltower.at(bin).fwdResolHist->Write();
     1214  }
     1215
     1216  histos_el.first->Write();
     1217  histos_el.second->Write();
     1218  histos_eltrack.first->Write();
     1219  histos_eltrack.second->Write();
     1220  histos_eltower.first->Write();
     1221  histos_eltower.second->Write();
     1222  C_el1->Write();
     1223  C_el2->Write();
     1224
     1225  histos_mu.first->Write();
     1226  histos_mu.second->Write();
     1227  histos_mutrack.first->Write();
     1228  histos_mutrack.second->Write();
     1229  C_mu1->Write();
     1230  C_mu->Write();
     1231
     1232  histos_ph.first->Write();
     1233  histos_ph.second->Write();
     1234  C_ph1->Write();
     1235  C_ph->Write();
     1236
     1237  for (int bin = 0; bin < Nbins; bin++)
     1238  {
     1239    plots_pfjets.at(bin).cenResolHist->Write();
     1240    plots_pfjets.at(bin).fwdResolHist->Write();
     1241    plots_calojets.at(bin).cenResolHist->Write();
     1242    plots_calojets.at(bin).fwdResolHist->Write();
     1243    plots_met.at(bin).cenResolHist->Write();
     1244  }
     1245  histos_btag.first->Write();
     1246  histos_btag.second->Write();
     1247  histos_tautag.first->Write();
     1248  histos_tautag.second->Write();
     1249  C_btag1->Write();
     1250  C_tautag1->Write();
     1251  C_jet->Write();
     1252  C_met->Write();
     1253
     1254  fout->Write();
     1255
     1256  cout << "** Exiting..." << endl;
     1257
     1258  delete treeReaderElectron;
     1259  delete treeReaderMuon;
     1260  delete treeReaderPhoton;
     1261  delete treeReaderJet;
     1262  delete treeReaderBJet;
     1263  delete treeReaderTauJet;
     1264  delete chainElectron;
     1265  delete chainMuon;
     1266  delete chainPhoton;
     1267  delete chainJet;
     1268  delete chainBJet;
     1269  delete chainTauJet;
     1270}
    521271
    531272//------------------------------------------------------------------------------
    … …  
    591278  if(argc != 3)
    601279  {
    61     cout << " Usage: " << appName << " input_file output_file" << endl;
    62     cout << " input_file  - input file in ROOT format ('Delphes' tree)," << endl;
     1280    cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
     1281    cout << " input_file_electron  - input file in ROOT format ('Delphes' tree)," << endl;
     1282    cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
     1283    cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
     1284    cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
     1285    cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
     1286    cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
    631287    cout << " output_file - output file in ROOT format" << endl;
    641288    return 1;
    … …  
    711295  TApplication app(appName, &appargc, appargv);
    721296
    73   TString inputFile(argv[1]);
    74   TString outputFile(argv[2]);
    75 
    76 
    77 //------------------------------------------------------------------------------
    78 // Here you call your macro's main function
    79 
    80   Validation(inputFile, outputFile);
    81 
    82 //------------------------------------------------------------------------------
    83 
    84 }
    85 
    86 
     1297  Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
     1298}
     1299
     1300

  • modules/Calorimeter.cc

    r7bb13cd r70b9632  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fECalResolutionFormula = new DelphesFormula;
    6362  fHCalResolutionFormula = new DelphesFormula;
    6463
    65   for(i = 0; i < 2; ++i)
    66   {
    67     fECalTowerTrackArray[i] = new TObjArray;
    68     fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
    69 
    70     fHCalTowerTrackArray[i] = new TObjArray;
    71     fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
    72   }
     64  fECalTowerTrackArray = new TObjArray;
     65  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
     66
     67  fHCalTowerTrackArray = new TObjArray;
     68  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
     69 
    7370}
    7471
    … …  
    7774Calorimeter::~Calorimeter()
    7875{
    79   Int_t i;
    80 
     76 
    8177  if(fECalResolutionFormula) delete fECalResolutionFormula;
    8278  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    8379
    84   for(i = 0; i < 2; ++i)
    85   {
    86     if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
    87     if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
    88 
    89     if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
    90     if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
    91   }
     80  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
     81  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
     82
     83  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
     84  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
     85 
    9286}
    9387
    … …  
    218212  Double_t ecalFraction, hcalFraction;
    219213  Double_t ecalEnergy, hcalEnergy;
    220   Double_t ecalSigma, hcalSigma;
    221214  Int_t pdgCode;
    222215
    … …  
    368361      fHCalTowerEnergy = 0.0;
    369362
    370       fECalTrackEnergy[0] = 0.0;
    371       fECalTrackEnergy[1] = 0.0;
    372 
    373       fHCalTrackEnergy[0] = 0.0;
    374       fHCalTrackEnergy[1] = 0.0;
    375 
     363      fECalTrackEnergy = 0.0;
     364      fHCalTrackEnergy = 0.0;
     365     
     366      fECalTrackSigma = 0.0;
     367      fHCalTrackSigma = 0.0;
     368     
    376369      fTowerTrackHits = 0;
    377370      fTowerPhotonHits = 0;
    378371
    379       fECalTowerTrackArray[0]->Clear();
    380       fECalTowerTrackArray[1]->Clear();
    381 
    382       fHCalTowerTrackArray[0]->Clear();
    383       fHCalTowerTrackArray[1]->Clear();
     372      fECalTowerTrackArray->Clear();
     373      fHCalTowerTrackArray->Clear();
     374   
    384375    }
    385376
    … …  
    406397      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    407398      {
    408         ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    409         if(ecalSigma/momentum.E() < track->TrackResolution)
    410         {
    411           fECalTrackEnergy[0] += ecalEnergy;
    412           fECalTowerTrackArray[0]->Add(track);
    413         }
    414         else
    415         {
    416           fECalTrackEnergy[1] += ecalEnergy;
    417           fECalTowerTrackArray[1]->Add(track);
    418         }
     399        fECalTrackEnergy += ecalEnergy;
     400        fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     401        fECalTowerTrackArray->Add(track);
    419402      }
     403     
    420404      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    421405      {
    422         hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    423         if(hcalSigma/momentum.E() < track->TrackResolution)
    424         {
    425           fHCalTrackEnergy[0] += hcalEnergy;
    426           fHCalTowerTrackArray[0]->Add(track);
    427         }
    428         else
    429         {
    430           fHCalTrackEnergy[1] += hcalEnergy;
    431           fHCalTowerTrackArray[1]->Add(track);
    432         }
     406        fHCalTrackEnergy += hcalEnergy;
     407        fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     408        fHCalTowerTrackArray->Add(track);
    433409      }
     410     
    434411      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    435412      {
    … …  
    476453  Double_t energy, pt, eta, phi;
    477454  Double_t ecalEnergy, hcalEnergy;
     455  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
     456 
    478457  Double_t ecalSigma, hcalSigma;
    479 
     458  Double_t ecalNeutralSigma, hcalNeutralSigma;
     459
     460  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
     461 
    480462  TLorentzVector momentum;
    481463  TFractionMap::iterator itFractionMap;
    … …  
    484466
    485467  if(!fTower) return;
    486 
    487468
    488469  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    554535    fTowerOutputArray->Add(fTower);
    555536  }
    556 
     537 
    557538  // fill energy flow candidates
    558 
    559   ecalEnergy -= fECalTrackEnergy[1];
    560   hcalEnergy -= fHCalTrackEnergy[1];
    561 
    562   fItECalTowerTrackArray[0]->Reset();
    563   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
    564   {
    565     mother = track;
    566     track = static_cast<Candidate*>(track->Clone());
    567     track->AddCandidate(mother);
    568 
    569     track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
    570 
    571     fEFlowTrackOutputArray->Add(track);
    572   }
    573 
    574   fItECalTowerTrackArray[1]->Reset();
    575   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
    576   {
    577     mother = track;
    578     track = static_cast<Candidate*>(track->Clone());
    579     track->AddCandidate(mother);
    580 
    581     fEFlowTrackOutputArray->Add(track);
    582   }
    583 
    584   fItHCalTowerTrackArray[0]->Reset();
    585   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
    586   {
    587     mother = track;
    588     track = static_cast<Candidate*>(track->Clone());
    589     track->AddCandidate(mother);
    590 
    591     track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
    592 
    593     fEFlowTrackOutputArray->Add(track);
    594   }
    595 
    596   fItHCalTowerTrackArray[1]->Reset();
    597   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
    598   {
    599     mother = track;
    600     track = static_cast<Candidate*>(track->Clone());
    601     track->AddCandidate(mother);
    602 
    603     fEFlowTrackOutputArray->Add(track);
    604   }
    605 
    606   if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
    607   if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0;
    608 
    609   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    610   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    611 
    612   if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
    613   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
    614 
    615   energy = ecalEnergy + hcalEnergy;
    616 
    617   if(ecalEnergy > 0.0)
     539  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
     540  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
     541
     542  //compute neutral excesses
     543  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
     544  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
     545 
     546  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
     547  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
     548 
     549   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
     550  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    618551  {
    619552    // create new photon tower
    620553    tower = static_cast<Candidate*>(fTower->Clone());
    621 
    622     pt = ecalEnergy / TMath::CosH(eta);
    623 
    624     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
    625     tower->Eem = ecalEnergy;
     554    pt =  ecalNeutralEnergy / TMath::CosH(eta);
     555   
     556    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     557    tower->Eem = ecalNeutralEnergy;
    626558    tower->Ehad = 0.0;
    627559    tower->PID = 22;
    628 
     560   
    629561    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron tower
     562   
     563    //clone tracks
     564    fItECalTowerTrackArray->Reset();
     565    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     566    {
     567      mother = track;
     568      track = static_cast<Candidate*>(track->Clone());
     569      track->AddCandidate(mother);
     570
     571      fEFlowTrackOutputArray->Add(track);
     572    }
     573 
     574  }
     575 
     576  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     577  else if(fECalTrackEnergy > 0.0)
     578  {
     579    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
     580    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
     581 
     582    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
     583    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     584
     585    //rescale tracks
     586    fItECalTowerTrackArray->Reset();
     587    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     588    { 
     589      mother = track;
     590      track = static_cast<Candidate*>(track->Clone());
     591      track->AddCandidate(mother);
     592
     593      track->Momentum *= rescaleFactor;
     594
     595      fEFlowTrackOutputArray->Add(track);
     596    }
     597  }
     598
     599
     600  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
     601  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
     602  {
     603    // create new photon tower
    634604    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     605    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     606   
     607    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     608    tower->Ehad = hcalNeutralEnergy;
    639609    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     610   
    642611    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     612   
     613    //clone tracks
     614    fItHCalTowerTrackArray->Reset();
     615    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     616    {
     617      mother = track;
     618      track = static_cast<Candidate*>(track->Clone());
     619      track->AddCandidate(mother);
     620
     621      fEFlowTrackOutputArray->Add(track);
     622    }
     623 
     624  }
     625 
     626  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     627  else if(fHCalTrackEnergy > 0.0)
     628  {
     629    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
     630    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
     631 
     632    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     633    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
     634
     635    //rescale tracks
     636    fItHCalTowerTrackArray->Reset();
     637    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     638    { 
     639      mother = track;
     640      track = static_cast<Candidate*>(track->Clone());
     641      track->AddCandidate(mother);
     642
     643      track->Momentum *= rescaleFactor;
     644
     645      fEFlowTrackOutputArray->Add(track);
     646    }
     647  }
     648 
     649 
    644650}
    645651

  • modules/Calorimeter.h

    r7bb13cd r70b9632  
    5858  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    5959  Double_t fECalTowerEnergy, fHCalTowerEnergy;
    60   Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
     60  Double_t fECalTrackEnergy, fHCalTrackEnergy;
    6161
    6262  Double_t fTimingEnergyMin;
    … …  
    7070  Double_t fECalEnergySignificanceMin;
    7171  Double_t fHCalEnergySignificanceMin;
     72
     73  Double_t fECalTrackSigma;
     74  Double_t fHCalTrackSigma;
    7275
    7376  Bool_t fSmearTowerCenter;
    … …  
    103106  TObjArray *fEFlowNeutralHadronOutputArray; //!
    104107
    105   TObjArray *fECalTowerTrackArray[2]; //!
    106   TIterator *fItECalTowerTrackArray[2]; //!
     108  TObjArray *fECalTowerTrackArray; //!
     109  TIterator *fItECalTowerTrackArray; //!
    107110
    108   TObjArray *fHCalTowerTrackArray[2]; //!
    109   TIterator *fItHCalTowerTrackArray[2]; //!
     111  TObjArray *fHCalTowerTrackArray; //!
     112  TIterator *fItHCalTowerTrackArray; //!
    110113
    111114  void FinalizeTower();

  • modules/ImpactParameterSmearing.cc

    r7bb13cd r70b9632  
    9696{
    9797  Candidate *candidate, *particle, *mother;
    98   Double_t xd, yd, zd, dxy, sx, sy, sz, ddxy;
     98  Double_t xd, yd, zd, d0, sx, sy, sz, dd0;
    9999  Double_t pt, eta, px, py, phi, e;
    100100
    … …  
    103103  {
    104104
    105     // take momentum before smearing (otherwise apply double smearing on dxy)
     105    // take momentum before smearing (otherwise apply double smearing on d0)
    106106    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    107107
    … …  
    131131
    132132    // calculate impact parameter (after-smearing)
    133     dxy = (xd*py - yd*px)/pt;
     133    d0 = (xd*py - yd*px)/pt;
    134134
    135     ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
     135    dd0 = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
    136136
    137137    // fill smeared values in candidate
    … …  
    143143    candidate->Zd = zd;
    144144
    145     candidate->Dxy = dxy;
    146     candidate->SDxy = ddxy;
     145    candidate->D0 = d0;
     146    candidate->ErrorD0 = dd0;
    147147
    148148    candidate->AddCandidate(mother);

  • modules/ModulesLinkDef.h

    r7bb13cd r70b9632  
    3535#include "modules/EnergySmearing.h"
    3636#include "modules/MomentumSmearing.h"
     37#include "modules/TrackSmearing.h"
    3738#include "modules/ImpactParameterSmearing.h"
    3839#include "modules/TimeSmearing.h"
    … …  
    5859#include "modules/StatusPidFilter.h"
    5960#include "modules/PdgCodeFilter.h"
     61#include "modules/BeamSpotFilter.h"
     62#include "modules/RecoPuFilter.h"
    6063#include "modules/Cloner.h"
    6164#include "modules/Weighter.h"
    … …  
    6366#include "modules/JetFlavorAssociation.h"
    6467#include "modules/JetFakeParticle.h"
     68#include "modules/VertexSorter.h"
     69#include "modules/VertexFinder.h"
     70#include "modules/VertexFinderDA4D.h"
    6571#include "modules/ExampleModule.h"
    6672
    … …  
    8086#pragma link C++ class EnergySmearing+;
    8187#pragma link C++ class MomentumSmearing+;
     88#pragma link C++ class TrackSmearing+;
    8289#pragma link C++ class ImpactParameterSmearing+;
    8390#pragma link C++ class TimeSmearing+;
    … …  
    103110#pragma link C++ class StatusPidFilter+;
    104111#pragma link C++ class PdgCodeFilter+;
     112#pragma link C++ class BeamSpotFilter+;
     113#pragma link C++ class RecoPuFilter+;
    105114#pragma link C++ class Cloner+;
    106115#pragma link C++ class Weighter+;
    … …  
    108117#pragma link C++ class JetFlavorAssociation+;
    109118#pragma link C++ class JetFakeParticle+;
     119#pragma link C++ class VertexSorter+;
     120#pragma link C++ class VertexFinder+;
     121#pragma link C++ class VertexFinderDA4D+;
    110122#pragma link C++ class ExampleModule+;
    111123

  • modules/ParticlePropagator.cc

    r7bb13cd r70b9632  
    6767}
    6868
     69
    6970//------------------------------------------------------------------------------
    7071
    … …  
    9192  fItInputArray = fInputArray->MakeIterator();
    9293
     94  // import beamspot
     95  try
     96  {
     97    fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
     98  }
     99  catch(runtime_error &e)
     100  {
     101    fBeamSpotInputArray = 0;
     102  }
    93103  // create output arrays
    94104
    … …  
    111121{
    112122  Candidate *candidate, *mother;
    113   TLorentzVector candidatePosition, candidateMomentum;
     123  TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
    114124  Double_t px, py, pz, pt, pt2, e, q;
    115125  Double_t x, y, z, t, r, phi;
    … …  
    120130  Double_t tmp, discr, discr2;
    121131  Double_t delta, gammam, omega, asinrho;
    122   Double_t rcu, rc2, dxy, xd, yd, zd;
     132  Double_t rcu, rc2, xd, yd, zd;
     133  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
     134  Double_t bsx, bsy, bsz;
    123135
    124136  const Double_t c_light = 2.99792458E8;
     137
     138  if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
     139    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
     140  else
     141  {
     142    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
     143    beamSpotPosition = beamSpotCandidate.Position;
     144  }
    125145
    126146  fItInputArray->Reset();
    … …  
    132152    y = candidatePosition.Y()*1.0E-3;
    133153    z = candidatePosition.Z()*1.0E-3;
     154
     155    bsx = beamSpotPosition.X()*1.0E-3;
     156    bsy = beamSpotPosition.Y()*1.0E-3;
     157    bsz = beamSpotPosition.Z()*1.0E-3;
     158
    134159    q = candidate->Charge;
    135160
    … …  
    182207      z_t = z + pz*t;
    183208
     209      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
     210
    184211      mother = candidate;
    185212      candidate = static_cast<Candidate*>(candidate->Clone());
    186213
     214      candidate->InitialPosition = candidatePosition;
    187215      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
     216      candidate->L = l*1.0E3;
    188217
    189218      candidate->Momentum = candidateMomentum;
    … …  
    239268      zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
    240269
    241       // calculate impact paramater
    242       dxy = (xd*py - yd*px)/pt;
     270      // use perigee momentum rather than original particle
     271      // momentum, since the orignal particle momentum isn't known
     272
     273      px = TMath::Sign(1.0,r) * pt * (-y_c / r_c);
     274      py = TMath::Sign(1.0,r) * pt * (x_c / r_c);
     275      etap = candidateMomentum.Eta();
     276      phip = TMath::ATan2(py, px);
     277
     278      candidateMomentum.SetPtEtaPhiE(pt, etap, phip, candidateMomentum.E());
     279
     280      // calculate additional track parameters (correct for beamspot position)
     281
     282      d0        = (  (x - bsx) * py - (y - bsy) * px) / pt;
     283      dz        = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
     284      p         = candidateMomentum.P();
     285      ctgTheta  = 1.0 / TMath::Tan (candidateMomentum.Theta ());
     286
    243287
    244288      // 3. time evaluation t = TMath::Min(t_r, t_z)
    … …  
    287331      r_t = TMath::Hypot(x_t, y_t);
    288332
     333
     334      // compute path length for an helix
     335
     336      alpha = pz*1.0E9 / c_light / gammam;
     337      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
     338
    289339      if(r_t > 0.0)
    290340      {
     341
     342        // store these variables before cloning
     343        candidate->D0 = d0*1.0E3;
     344        candidate->DZ = dz*1.0E3;
     345        candidate->P  = p;
     346        candidate->PT = pt;
     347        candidate->CtgTheta = ctgTheta;
     348        candidate->Phi = phip;
     349
    291350        mother = candidate;
    292351        candidate = static_cast<Candidate*>(candidate->Clone());
    293352
     353        candidate->InitialPosition = candidatePosition;
    294354        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
    295355
    296356        candidate->Momentum = candidateMomentum;
    297         candidate->Dxy = dxy*1.0E3;
    298         candidate->Xd = xd*1.0E3;
     357
     358            candidate->L  =  l*1.0E3;
     359
     360            candidate->Xd = xd*1.0E3;
    299361        candidate->Yd = yd*1.0E3;
    300362        candidate->Zd = zd*1.0E3;

  • modules/ParticlePropagator.h

    r7bb13cd r70b9632  
    3535class TClonesArray;
    3636class TIterator;
     37class TLorentzVector;
    3738
    3839class ParticlePropagator: public DelphesModule
    … …  
    5556
    5657  const TObjArray *fInputArray; //!
     58  const TObjArray *fBeamSpotInputArray; //!
    5759
    5860  TObjArray *fOutputArray; //!

  • modules/PhotonConversions.cc

    r7bb13cd r70b9632  
    5959  fItInputArray(0), fConversionMap(0), fDecayXsec(0)
    6060{
    61   fDecayXsec = new TF1;
     61  fDecayXsec = new TF1("decayXsec","1.0 - 4.0/3.0 * x * (1.0 - x)", 0.0, 1.0);
    6262  fConversionMap = new DelphesCylindricalFormula;
    6363}
    … …  
    8484
    8585  fConversionMap->Compile(GetString("ConversionMap", "0.0"));
    86 
    87 #if  ROOT_VERSION_CODE < ROOT_VERSION(6,04,00)
    88   fDecayXsec->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
    89 #else
    90   fDecayXsec->GetFormula()->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
    91 #endif
    92   fDecayXsec->SetRange(0.0, 1.0);
    9386
    9487  // import array with output from filter/classifier module

  • modules/PileUpMerger.cc

    r7bb13cd r70b9632  
    115115  TDatabasePDG *pdg = TDatabasePDG::Instance();
    116116  TParticlePDG *pdgParticle;
    117   Int_t pid;
     117  Int_t pid, nch, nvtx = -1;
    118118  Float_t x, y, z, t, vx, vy;
    119   Float_t px, py, pz, e;
    120   Double_t dz, dphi, dt;
     119  Float_t px, py, pz, e, pt;
     120  Double_t dz, dphi, dt, sumpt2, dz0, dt0;
    121121  Int_t numberOfEvents, event, numberOfParticles;
    122122  Long64_t allEntries, entry;
    … …  
    132132  fFunction->GetRandom2(dz, dt);
    133133
     134  dz0 = -1.0e6;
     135  dt0 = -1.0e6;
     136
    134137  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    135138  dz *= 1.0E3; // necessary in order to make z in mm
     139
     140  //cout<<dz<<","<<dt<<endl;
     141
    136142  vx = 0.0;
    137143  vy = 0.0;
     144
    138145  numberOfParticles = fInputArray->GetEntriesFast();
     146  nch = 0;
     147  sumpt2 = 0.0;
     148
     149  factory = GetFactory();
     150  vertex = factory->NewCandidate();
     151
    139152  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    140153  {
    … …  
    143156    z = candidate->Position.Z();
    144157    t = candidate->Position.T();
    145     candidate->Position.SetZ(z + dz);
    146     candidate->Position.SetT(t + dt);
     158    pt = candidate->Momentum.Pt();
     159
     160    // take postion and time from first stable particle
     161    if (dz0 < -999999.0)
     162      dz0 = z;
     163    if (dt0 < -999999.0)
     164      dt0 = t;
     165
     166    // cancel any possible offset in position and time the input file
     167    candidate->Position.SetZ(z - dz0 + dz);
     168    candidate->Position.SetT(t - dt0 + dt);
     169
     170    candidate->IsPU = 0;
     171
    147172    fParticleOutputArray->Add(candidate);
     173
     174    if(TMath::Abs(candidate->Charge) >  1.0E-9)
     175    {
     176      nch++;
     177      sumpt2 += pt*pt;
     178      vertex->AddCandidate(candidate);
     179    }
    148180  }
    149181
    150182  if(numberOfParticles > 0)
    151183  {
    152     vx /= numberOfParticles;
    153     vy /= numberOfParticles;
     184    vx /= sumpt2;
     185    vy /= sumpt2;
    154186  }
    155187
    156   factory = GetFactory();
    157 
    158   vertex = factory->NewCandidate();
     188  nvtx++;
    159189  vertex->Position.SetXYZT(vx, vy, dz, dt);
     190  vertex->ClusterIndex = nvtx;
     191  vertex->ClusterNDF = nch;
     192  vertex->SumPT2 = sumpt2;
     193  vertex->GenSumPT2 = sumpt2;
    160194  fVertexOutputArray->Add(vertex);
    161195
    … …  
    170204      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
    171205      break;
     206    case 2:
     207      numberOfEvents = fMeanPileUp;
     208      break;
    172209    default:
    173210      numberOfEvents = gRandom->Poisson(fMeanPileUp);
    … …  
    176213
    177214  allEntries = fReader->GetEntries();
     215
    178216
    179217  for(event = 0; event < numberOfEvents; ++event)
    … …  
    198236    vx = 0.0;
    199237    vy = 0.0;
     238
    200239    numberOfParticles = 0;
     240    sumpt2 = 0.0;
     241
     242    //factory = GetFactory();
     243    vertex = factory->NewCandidate();
     244
    201245    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    202246    {
    … …  
    215259      candidate->Momentum.SetPxPyPzE(px, py, pz, e);
    216260      candidate->Momentum.RotateZ(dphi);
     261      pt = candidate->Momentum.Pt();
    217262
    218263      x -= fInputBeamSpotX;
    … …  
    224269      vx += candidate->Position.X();
    225270      vy += candidate->Position.Y();
     271
    226272      ++numberOfParticles;
     273      if(TMath::Abs(candidate->Charge) >  1.0E-9)
     274      {
     275        nch++;
     276        sumpt2 += pt*pt;
     277        vertex->AddCandidate(candidate);
     278      }
    227279
    228280      fParticleOutputArray->Add(candidate);
    … …  
    235287    }
    236288
    237     vertex = factory->NewCandidate();
     289    nvtx++;
     290
    238291    vertex->Position.SetXYZT(vx, vy, dz, dt);
     292
     293    vertex->ClusterIndex = nvtx;
     294    vertex->ClusterNDF = nch;
     295    vertex->SumPT2 = sumpt2;
     296    vertex->GenSumPT2 = sumpt2;
     297
    239298    vertex->IsPU = 1;
    240299
    241300    fVertexOutputArray->Add(vertex);
     301
    242302  }
    243303}

  • modules/SimpleCalorimeter.cc

    r7bb13cd r70b9632  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fResolutionFormula = new DelphesFormula;
    63 
    64   for(i = 0; i < 2; ++i)
    65   {
    66     fTowerTrackArray[i] = new TObjArray;
    67     fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
    68   }
     62  fTowerTrackArray = new TObjArray;
     63  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
     64 
    6965}
    7066
    … …  
    7369SimpleCalorimeter::~SimpleCalorimeter()
    7470{
    75   Int_t i;
    76 
     71 
    7772  if(fResolutionFormula) delete fResolutionFormula;
    78 
    79   for(i = 0; i < 2; ++i)
    80   {
    81     if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
    82     if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
    83   }
     73  if(fTowerTrackArray) delete fTowerTrackArray;
     74  if(fItTowerTrackArray) delete fItTowerTrackArray;
     75 
    8476}
    8577
    … …  
    198190  Double_t fraction;
    199191  Double_t energy;
    200   Double_t sigma;
    201192  Int_t pdgCode;
    202193
    … …  
    340331      fTowerEnergy = 0.0;
    341332
    342       fTrackEnergy[0] = 0.0;
    343       fTrackEnergy[1] = 0.0;
     333      fTrackEnergy = 0.0;
     334      fTrackSigma = 0.0;
    344335
    345336      fTowerTime = 0.0;
    … …  
    351342      fTowerPhotonHits = 0;
    352343
    353       fTowerTrackArray[0]->Clear();
    354       fTowerTrackArray[1]->Clear();
    355     }
     344      fTowerTrackArray->Clear();
     345     }
    356346
    357347    // check for track hits
    … …  
    371361      if(fTrackFractions[number] > 1.0E-9)
    372362      {
    373         sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    374         if(sigma/momentum.E() < track->TrackResolution)
    375         {
    376           fTrackEnergy[0] += energy;
    377           fTowerTrackArray[0]->Add(track);
    378         }
    379         else
    380         {
    381           fTrackEnergy[1] += energy;
    382           fTowerTrackArray[1]->Add(track);
    383         }
     363             
     364       // compute total charged energy   
     365       fTrackEnergy += energy;
     366       fTrackSigma += ((track->TrackResolution)*momentum.E())*((track->TrackResolution)*momentum.E());
     367       
     368       fTowerTrackArray->Add(track);
     369     
    384370      }
     371       
    385372      else
    386373      {
    … …  
    403390    fTowerEnergy += energy;
    404391
    405     fTowerTime += TMath::Sqrt(energy)*position.T();
    406     fTowerTimeWeight += TMath::Sqrt(energy);
     392    fTowerTime += energy*position.T();
     393    fTowerTimeWeight += energy;
    407394
    408395    fTower->AddCandidate(particle);
    … …  
    418405{
    419406  Candidate *tower, *track, *mother;
    420   Double_t energy, pt, eta, phi;
    421   Double_t sigma;
     407  Double_t energy,neutralEnergy, pt, eta, phi;
     408  Double_t sigma, neutralSigma;
    422409  Double_t time;
     410   
     411  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    423412
    424413  TLorentzVector momentum;
    … …  
    436425
    437426  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
     427
    438428
    439429  if(fSmearTowerCenter)
    … …  
    464454  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    465455
    466   // fill e-flow candidates
    467 
    468   energy -= fTrackEnergy[1];
    469 
    470   fItTowerTrackArray[0]->Reset();
    471   while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
    472   {
    473     mother = track;
    474     track = static_cast<Candidate*>(track->Clone());
    475     track->AddCandidate(mother);
    476 
    477     track->Momentum *= energy/fTrackEnergy[0];
    478 
    479     fEFlowTrackOutputArray->Add(track);
    480   }
    481 
    482   fItTowerTrackArray[1]->Reset();
    483   while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
    484   {
    485     mother = track;
    486     track = static_cast<Candidate*>(track->Clone());
    487     track->AddCandidate(mother);
    488 
    489     fEFlowTrackOutputArray->Add(track);
    490   }
    491 
    492   if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
    493 
    494   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
    495   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    496 
    497   // save energy excess as an energy flow tower
    498   if(energy > 0.0)
     456 
     457  // e-flow candidates
     458
     459  //compute neutral excess
     460 
     461  fTrackSigma = TMath::Sqrt(fTrackSigma);
     462  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
     463 
     464  //compute sigma_trk total
     465  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
     466   
     467  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
     468  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
    499469  {
    500470    // create new photon tower
    501471    tower = static_cast<Candidate*>(fTower->Clone());
    502     pt = energy / TMath::CosH(eta);
    503 
    504     tower->Eem = (!fIsEcal) ? 0 : energy;
    505     tower->Ehad = (fIsEcal) ? 0 : energy;
     472    pt = neutralEnergy / TMath::CosH(eta);
     473
     474    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
     475    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     476    tower->PID = (fIsEcal) ? 22 : 0;
     477 
     478    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     479    fEFlowTowerOutputArray->Add(tower);
    506480   
    507     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     481    fItTowerTrackArray->Reset();
     482    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     483    {
     484      mother = track;
     485      track = static_cast<Candidate*>(track->Clone());
     486      track->AddCandidate(mother);
     487
     488      fEFlowTrackOutputArray->Add(track);
     489    }
     490  }
    508491   
    509     tower->PID = (fIsEcal) ? 22 : 0;
    510    
    511     fEFlowTowerOutputArray->Add(tower);
    512   }
    513 }
     492  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     493  else if (fTrackEnergy > 0.0)
     494  {
     495    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
     496    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
     497 
     498    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
     499    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
     500   
     501    fItTowerTrackArray->Reset();
     502    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     503    {
     504      mother = track;
     505      track = static_cast<Candidate*>(track->Clone());
     506      track->AddCandidate(mother);
     507
     508      track->Momentum *= rescaleFactor;
     509
     510      fEFlowTrackOutputArray->Add(track);
     511    }
     512  }
     513   
     514 }
    514515
    515516//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

    r7bb13cd r70b9632  
    5959  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    6060  Double_t fTowerEnergy;
    61   Double_t fTrackEnergy[2];
     61  Double_t fTrackEnergy;
    6262
    6363  Double_t fTowerTime;
    … …  
    7272
    7373  Double_t fEnergySignificanceMin;
     74
     75  Double_t fTrackSigma;
    7476
    7577  Bool_t fSmearTowerCenter;
    … …  
    102104  TObjArray *fEFlowTowerOutputArray; //!
    103105
    104   TObjArray *fTowerTrackArray[2]; //!
    105   TIterator *fItTowerTrackArray[2]; //!
     106  TObjArray *fTowerTrackArray; //!
     107  TIterator *fItTowerTrackArray; //!
    106108
    107109  void FinalizeTower();

  • modules/TimeSmearing.cc

    r7bb13cd r70b9632  
    9393{
    9494  Candidate *candidate, *mother;
    95   Double_t t;
     95  Double_t ti, tf_smeared, tf;
    9696  const Double_t c_light = 2.99792458E8;
    9797
    … …  
    9999  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    100100  {
    101     const TLorentzVector &candidatePosition = candidate->Position;
    102     t = candidatePosition.T()*1.0E-3/c_light;
     101    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
     102    const TLorentzVector &candidateFinalPosition = candidate->Position;
     103
     104    ti = candidateInitialPosition.T()*1.0E-3/c_light;
     105    tf = candidateFinalPosition.T()*1.0E-3/c_light;
    103106
    104107    // apply smearing formula
    105     t = gRandom->Gaus(t, fTimeResolution);
     108    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
     109    ti = ti + tf_smeared - tf;
    106110
    107111    mother = candidate;
    108112    candidate = static_cast<Candidate*>(candidate->Clone());
    109     candidate->Position.SetT(t*1.0E3*c_light);
     113    candidate->InitialPosition.SetT(ti*1.0E3*c_light);
     114    candidate->Position.SetT(tf*1.0E3*c_light);
     115
     116    candidate->ErrorT = fTimeResolution*1.0E3*c_light;
    110117
    111118    candidate->AddCandidate(mother);

  • modules/TrackCountingBTagging.cc

    r7bb13cd r70b9632  
    9696
    9797  Double_t jpx, jpy;
    98   Double_t dr, tpx, tpy, tpt;
    99   Double_t xd, yd, dxy, ddxy, ip, sip;
     98  Double_t dr, tpt;
     99  Double_t xd, yd, d0, dd0, ip, sip;
    100100
    101101  Int_t sign;
    … …  
    117117    {
    118118      const TLorentzVector &trkMomentum = track->Momentum;
    119 
     119     
    120120      dr = jetMomentum.DeltaR(trkMomentum);
    121 
    122121      tpt = trkMomentum.Pt();
    123       tpx = trkMomentum.Px();
    124       tpy = trkMomentum.Py();
    125 
    126122      xd = track->Xd;
    127123      yd = track->Yd;
    128       dxy = TMath::Hypot(xd, yd);
    129       ddxy = track->SDxy;
     124      d0 = TMath::Hypot(xd, yd);
     125      dd0 = track->ErrorD0;
    130126
    131127      if(tpt < fPtMin) continue;
    132128      if(dr > fDeltaR) continue;
    133       if(dxy > fIPmax) continue;
     129      if(d0 > fIPmax) continue;
    134130
    135131      sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
    136132
    137       ip = sign*dxy;
    138       sip = ip / TMath::Abs(ddxy);
     133      ip = sign*d0;
     134      sip = ip / TMath::Abs(dd0);
    139135
    140136      if(sip > fSigMin) count++;

  • modules/TreeWriter.cc

    r7bb13cd r70b9632  
    215215    entry->Pz = momentum.Pz();
    216216
     217    entry->D0            = candidate->D0;
     218    entry->DZ            = candidate->DZ;
     219    entry->P             = candidate->P;
     220    entry->PT            = candidate->PT;
     221    entry->CtgTheta      = candidate->CtgTheta;
     222    entry->Phi           = candidate->Phi;
     223
    217224    entry->Eta = eta;
    218225    entry->Phi = momentum.Phi();
    … …  
    233240{
    234241  TIter iterator(array);
    235   Candidate *candidate = 0;
     242  Candidate *candidate = 0, *constituent = 0;
    236243  Vertex *entry = 0;
    237244
    238245  const Double_t c_light = 2.99792458E8;
    239246
     247  Double_t x, y, z, t, xError, yError, zError, tError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
     248  UInt_t index, ndf;
     249
     250  CompBase *compare = Candidate::fgCompare;
     251  Candidate::fgCompare = CompSumPT2<Candidate>::Instance();
     252  array->Sort();
     253  Candidate::fgCompare = compare;
     254
    240255  // loop over all vertices
    241256  iterator.Reset();
    242257  while((candidate = static_cast<Candidate*>(iterator.Next())))
    243258  {
    244     const TLorentzVector &position = candidate->Position;
     259
     260    index = candidate->ClusterIndex;
     261    ndf = candidate->ClusterNDF;
     262    sigma = candidate->ClusterSigma;
     263    sumPT2 = candidate->SumPT2;
     264    btvSumPT2 = candidate->BTVSumPT2;
     265    genDeltaZ = candidate->GenDeltaZ;
     266    genSumPT2 = candidate->GenSumPT2;
     267
     268    x = candidate->Position.X();
     269    y = candidate->Position.Y();
     270    z = candidate->Position.Z();
     271    t = candidate->Position.T()*1.0E-3/c_light;
     272
     273    xError = candidate->PositionError.X ();
     274    yError = candidate->PositionError.Y ();
     275    zError = candidate->PositionError.Z ();
     276    tError = candidate->PositionError.T ()*1.0E-3/c_light;
    245277
    246278    entry = static_cast<Vertex*>(branch->NewEntry());
    247279
    248     entry->X = position.X();
    249     entry->Y = position.Y();
    250     entry->Z = position.Z();
    251     entry->T = position.T()*1.0E-3/c_light;
    252   }
    253 }
     280    entry->Index = index;
     281    entry->NDF = ndf;
     282    entry->Sigma = sigma;
     283    entry->SumPT2 = sumPT2;
     284    entry->BTVSumPT2 = btvSumPT2;
     285    entry->GenDeltaZ = genDeltaZ;
     286    entry->GenSumPT2 = genSumPT2;
     287
     288    entry->X = x;
     289    entry->Y = y;
     290    entry->Z = z;
     291    entry->T = t;
     292
     293    entry->ErrorX = xError;
     294    entry->ErrorY = yError;
     295    entry->ErrorZ = zError;
     296    entry->ErrorT = tError;
     297
     298
     299    TIter itConstituents(candidate->GetCandidates());
     300    itConstituents.Reset();
     301    entry->Constituents.Clear();
     302    while((constituent = static_cast<Candidate*>(itConstituents.Next())))
     303    {
     304      entry->Constituents.Add(constituent);
     305    }
     306
     307  }
     308}
     309
    254310
    255311//------------------------------------------------------------------------------
    … …  
    261317  Candidate *particle = 0;
    262318  Track *entry = 0;
    263   Double_t pt, signz, cosTheta, eta, rapidity;
     319  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
    264320  const Double_t c_light = 2.99792458E8;
    265321
    … …  
    292348    entry->TOuter = position.T()*1.0E-3/c_light;
    293349
    294     entry->Dxy = candidate->Dxy;
    295     entry->SDxy = candidate->SDxy ;
     350    entry->L = candidate->L;
     351
     352    entry->D0            = candidate->D0;
     353    entry->ErrorD0       = candidate->ErrorD0;
     354    entry->DZ            = candidate->DZ;
     355    entry->ErrorDZ       = candidate->ErrorDZ;
     356
     357    entry->ErrorP        = candidate->ErrorP;
     358    entry->ErrorPT       = candidate->ErrorPT;
     359    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
     360    entry->ErrorPhi      = candidate->ErrorPhi;
     361
    296362    entry->Xd = candidate->Xd;
    297363    entry->Yd = candidate->Yd;
    … …  
    301367
    302368    pt = momentum.Pt();
     369    p = momentum.P();
     370    phi = momentum.Phi();
     371    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1/TMath::Tan(momentum.Theta()) : 1e10;
     372
    303373    cosTheta = TMath::Abs(momentum.CosTheta());
    304374    signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
    … …  
    306376    rapidity = (cosTheta == 1.0 ? signz*999.9 : momentum.Rapidity());
    307377
     378    entry->PT  = pt;
    308379    entry->Eta = eta;
    309     entry->Phi = momentum.Phi();
    310     entry->PT = pt;
     380    entry->Phi = phi;
     381    entry->CtgTheta = ctgTheta;
    311382
    312383    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    … …  
    319390
    320391    entry->Particle = particle;
     392
     393    entry->VertexIndex = candidate->ClusterIndex;
     394
    321395  }
    322396}

  • readers/DelphesCMSFWLite.cpp

    r7bb13cd r70b9632  
    6565  ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchRwgt,
    6666  DelphesFactory *factory, TObjArray *allParticleOutputArray,
    67   TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray)
     67  TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray, Bool_t firstEvent)
    6868{
     69
    6970  fwlite::Handle< GenEventInfoProduct > handleGenEventInfo;
    70 
    7171  fwlite::Handle< LHEEventProduct > handleLHEEvent;
    72 
    7372  fwlite::Handle< vector< reco::GenParticle > > handleParticle;
     73
    7474  vector< reco::GenParticle >::const_iterator itParticle;
    7575
    … …  
    7878
    7979  handleGenEventInfo.getByLabel(event, "generator");
    80   handleLHEEvent.getByLabel(event, "externalLHEProducer");
    81   handleParticle.getByLabel(event, "genParticles");
     80
     81  if (!((handleLHEEvent.getBranchNameFor(event, "source")).empty()))
     82  {
     83    handleLHEEvent.getByLabel(event, "source");
     84  }
     85  else if (!((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty()))
     86  {
     87    handleLHEEvent.getByLabel(event, "externalLHEProducer");
     88  }
     89  else if (firstEvent)
     90  {
     91    std::cout<<"Wrong LHEEvent Label! Please, check the input file."<<std::endl;
     92  }
     93
     94  if (!((handleParticle.getBranchNameFor(event, "genParticles")).empty()))
     95  {
     96    handleParticle.getByLabel(event, "genParticles");
     97  }
     98  else if (!((handleParticle.getBranchNameFor(event, "prunedGenParticles")).empty()))
     99  {
     100    handleParticle.getByLabel(event, "prunedGenParticles");
     101  }
     102  else
     103  {
     104    std::cout<<"Wrong GenParticle Label! Please, check the input file."<<std::endl;
     105    exit(-1);
     106  }
     107
     108  Bool_t foundLHE = !((handleLHEEvent.getBranchNameFor(event, "source")).empty()) || !((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty());
    82109
    83110  HepMCEvent *element;
    … …  
    91118  Double_t px, py, pz, e, mass;
    92119  Double_t x, y, z;
    93 
    94   const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
    95   vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
    96120
    97121  element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
    … …  
    117141  element->ProcTime = 0.0;
    118142
    119   for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
    120   {
    121     weight = static_cast<Weight *>(branchRwgt->NewEntry());
    122     weight->Weight = itWeightsInfo->wgt;
     143
     144  if(foundLHE)
     145  {
     146    const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
     147    vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
     148   
     149    for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
     150    {
     151      weight = static_cast<Weight *>(branchRwgt->NewEntry());
     152      weight->Weight = itWeightsInfo->wgt;
     153    } 
    123154  }
    124155
    … …  
    207238  Int_t i;
    208239  Long64_t eventCounter, numberOfEvents;
     240  Bool_t firstEvent = kTRUE;
    209241
    210242  if(argc < 4)
    … …  
    240272
    241273    branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
    242     branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
     274    branchRwgt = treeWriter->NewBranch("Weight", Weight::Class());
    243275
    244276    confReader = new ExRootConfReader;
    … …  
    281313      modularDelphes->Clear();
    282314      treeWriter->Clear();
     315
    283316      for(event.toBegin(); !event.atEnd() && !interrupted; ++event)
    284317      {
    285318        ConvertInput(event, eventCounter, branchEvent, branchRwgt, factory,
    286           allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
     319          allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
    287320        modularDelphes->ProcessTask();
     321         
     322        firstEvent = kFALSE;
    288323
    289324        treeWriter->Fill();
Note: See TracChangeset for help on using the changeset viewer.