Fork me on GitHub

Changes in / [21eab4f:ec5e04b] in git


Ignore:
Files:
14 deleted
28 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r21eab4f rec5e04b  
    137137tmp/examples/Example1.$(ObjSuf): \
    138138        examples/Example1.cpp \
    139         classes/DelphesClasses.h \
    140         external/ExRootAnalysis/ExRootTreeReader.h \
    141         external/ExRootAnalysis/ExRootTreeWriter.h \
    142         external/ExRootAnalysis/ExRootTreeBranch.h \
    143         external/ExRootAnalysis/ExRootResult.h \
    144         external/ExRootAnalysis/ExRootUtilities.h
    145 Validation$(ExeSuf): \
    146         tmp/examples/Validation.$(ObjSuf)
    147 
    148 tmp/examples/Validation.$(ObjSuf): \
    149         examples/Validation.cpp \
    150139        classes/DelphesClasses.h \
    151140        external/ExRootAnalysis/ExRootTreeReader.h \
    … …  
    161150        root2pileup$(ExeSuf) \
    162151        stdhep2pileup$(ExeSuf) \
    163         Example1$(ExeSuf) \
    164         Validation$(ExeSuf)
     152        Example1$(ExeSuf)
    165153
    166154EXECUTABLE_OBJ +=  \
    … …  
    171159        tmp/converters/root2pileup.$(ObjSuf) \
    172160        tmp/converters/stdhep2pileup.$(ObjSuf) \
    173         tmp/examples/Example1.$(ObjSuf) \
    174         tmp/examples/Validation.$(ObjSuf)
     161        tmp/examples/Example1.$(ObjSuf)
    175162
    176163DelphesHepMC$(ExeSuf): \
    … …  
    196183        classes/DelphesLHEFReader.h \
    197184        external/ExRootAnalysis/ExRootTreeWriter.h \
    198         external/ExRootAnalysis/ExRootTreeBranch.h \
    199         external/ExRootAnalysis/ExRootProgressBar.h
    200 DelphesROOT$(ExeSuf): \
    201         tmp/readers/DelphesROOT.$(ObjSuf)
    202 
    203 tmp/readers/DelphesROOT.$(ObjSuf): \
    204         readers/DelphesROOT.cpp \
    205         modules/Delphes.h \
    206         classes/DelphesStream.h \
    207         classes/DelphesClasses.h \
    208         classes/DelphesFactory.h \
    209         external/ExRootAnalysis/ExRootTreeWriter.h \
    210         external/ExRootAnalysis/ExRootTreeReader.h \
    211185        external/ExRootAnalysis/ExRootTreeBranch.h \
    212186        external/ExRootAnalysis/ExRootProgressBar.h
    … …  
    226200        DelphesHepMC$(ExeSuf) \
    227201        DelphesLHEF$(ExeSuf) \
    228         DelphesROOT$(ExeSuf) \
    229202        DelphesSTDHEP$(ExeSuf)
    230203
    … …  
    232205        tmp/readers/DelphesHepMC.$(ObjSuf) \
    233206        tmp/readers/DelphesLHEF.$(ObjSuf) \
    234         tmp/readers/DelphesROOT.$(ObjSuf) \
    235207        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    236208
    … …  
    347319        modules/EnergySmearing.h \
    348320        modules/MomentumSmearing.h \
    349         modules/TrackSmearing.h \
    350321        modules/ImpactParameterSmearing.h \
    351322        modules/TimeSmearing.h \
    … …  
    371342        modules/StatusPidFilter.h \
    372343        modules/PdgCodeFilter.h \
    373         modules/BeamSpotFilter.h \
    374         modules/RecoPuFilter.h \
    375344        modules/Cloner.h \
    376345        modules/Weighter.h \
    … …  
    378347        modules/JetFlavorAssociation.h \
    379348        modules/JetFakeParticle.h \
    380         modules/VertexSorter.h \
    381         modules/VertexFinder.h \
    382         modules/VertexFinderDA4D.h \
    383349        modules/ExampleModule.h
    384350tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    587553        classes/DelphesFactory.h \
    588554        classes/DelphesFormula.h
    589 tmp/modules/BeamSpotFilter.$(ObjSuf): \
    590         modules/BeamSpotFilter.$(SrcSuf) \
    591         modules/BeamSpotFilter.h \
    592         classes/DelphesClasses.h \
    593         classes/DelphesFactory.h \
    594         classes/DelphesFormula.h \
    595         external/ExRootAnalysis/ExRootResult.h \
    596         external/ExRootAnalysis/ExRootFilter.h \
    597         external/ExRootAnalysis/ExRootClassifier.h
    598555tmp/modules/Calorimeter.$(ObjSuf): \
    599556        modules/Calorimeter.$(SrcSuf) \
    … …  
    828785        external/ExRootAnalysis/ExRootFilter.h \
    829786        external/ExRootAnalysis/ExRootClassifier.h
    830 tmp/modules/RecoPuFilter.$(ObjSuf): \
    831         modules/RecoPuFilter.$(SrcSuf) \
    832         modules/RecoPuFilter.h \
    833         classes/DelphesClasses.h \
    834         classes/DelphesFactory.h \
    835         classes/DelphesFormula.h \
    836         external/ExRootAnalysis/ExRootResult.h \
    837         external/ExRootAnalysis/ExRootFilter.h \
    838         external/ExRootAnalysis/ExRootClassifier.h
    839787tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    840788        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    898846        modules/TrackPileUpSubtractor.$(SrcSuf) \
    899847        modules/TrackPileUpSubtractor.h \
    900         classes/DelphesClasses.h \
    901         classes/DelphesFactory.h \
    902         classes/DelphesFormula.h \
    903         external/ExRootAnalysis/ExRootResult.h \
    904         external/ExRootAnalysis/ExRootFilter.h \
    905         external/ExRootAnalysis/ExRootClassifier.h
    906 tmp/modules/TrackSmearing.$(ObjSuf): \
    907         modules/TrackSmearing.$(SrcSuf) \
    908         modules/TrackSmearing.h \
    909848        classes/DelphesClasses.h \
    910849        classes/DelphesFactory.h \
    … …  
    929868        classes/DelphesFactory.h \
    930869        classes/DelphesFormula.h \
    931         external/ExRootAnalysis/ExRootResult.h \
    932         external/ExRootAnalysis/ExRootFilter.h \
    933         external/ExRootAnalysis/ExRootClassifier.h
    934 tmp/modules/VertexFinder.$(ObjSuf): \
    935         modules/VertexFinder.$(SrcSuf) \
    936         modules/VertexFinder.h \
    937         classes/DelphesClasses.h \
    938         classes/DelphesFactory.h \
    939         classes/DelphesFormula.h \
    940         classes/DelphesPileUpReader.h \
    941         external/ExRootAnalysis/ExRootResult.h \
    942         external/ExRootAnalysis/ExRootFilter.h \
    943         external/ExRootAnalysis/ExRootClassifier.h
    944 tmp/modules/VertexFinderDA4D.$(ObjSuf): \
    945         modules/VertexFinderDA4D.$(SrcSuf) \
    946         modules/VertexFinderDA4D.h \
    947         classes/DelphesClasses.h \
    948         classes/DelphesFactory.h \
    949         classes/DelphesFormula.h \
    950         classes/DelphesPileUpReader.h \
    951         external/ExRootAnalysis/ExRootResult.h \
    952         external/ExRootAnalysis/ExRootFilter.h \
    953         external/ExRootAnalysis/ExRootClassifier.h
    954 tmp/modules/VertexSorter.$(ObjSuf): \
    955         modules/VertexSorter.$(SrcSuf) \
    956         modules/VertexSorter.h \
    957         classes/DelphesClasses.h \
    958         classes/DelphesFactory.h \
    959         classes/DelphesFormula.h \
    960         classes/DelphesPileUpReader.h \
    961870        external/ExRootAnalysis/ExRootResult.h \
    962871        external/ExRootAnalysis/ExRootFilter.h \
    … …  
    1022931        tmp/modules/AngularSmearing.$(ObjSuf) \
    1023932        tmp/modules/BTagging.$(ObjSuf) \
    1024         tmp/modules/BeamSpotFilter.$(ObjSuf) \
    1025933        tmp/modules/Calorimeter.$(ObjSuf) \
    1026934        tmp/modules/Cloner.$(ObjSuf) \
    … …  
    1047955        tmp/modules/PileUpJetID.$(ObjSuf) \
    1048956        tmp/modules/PileUpMerger.$(ObjSuf) \
    1049         tmp/modules/RecoPuFilter.$(ObjSuf) \
    1050957        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    1051958        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    1056963        tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
    1057964        tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
    1058         tmp/modules/TrackSmearing.$(ObjSuf) \
    1059965        tmp/modules/TreeWriter.$(ObjSuf) \
    1060966        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
    1061         tmp/modules/VertexFinder.$(ObjSuf) \
    1062         tmp/modules/VertexFinderDA4D.$(ObjSuf) \
    1063         tmp/modules/VertexSorter.$(ObjSuf) \
    1064967        tmp/modules/Weighter.$(ObjSuf)
    1065968
    … …  
    16661569        tmp/external/tcl/tclVar.$(ObjSuf)
    16671570
    1668 modules/VertexFinderDA4D.h: \
    1669         classes/DelphesModule.h \
    1670         classes/DelphesClasses.h
    1671         @touch $@
    1672 
    1673 modules/TrackSmearing.h: \
    1674         classes/DelphesModule.h
    1675         @touch $@
    1676 
    16771571external/fastjet/ClusterSequence.hh: \
    16781572        external/fastjet/PseudoJet.hh \
    … …  
    19351829        @touch $@
    19361830
    1937 modules/VertexSorter.h: \
    1938         classes/DelphesModule.h \
    1939         classes/DelphesClasses.h
    1940         @touch $@
    1941 
    19421831modules/Delphes.h: \
    19431832        classes/DelphesModule.h
    1944         @touch $@
    1945 
    1946 modules/VertexFinder.h: \
    1947         classes/DelphesModule.h \
    1948         classes/DelphesClasses.h
    19491833        @touch $@
    19501834
    … …  
    20181902        @touch $@
    20191903
    2020 modules/RecoPuFilter.h: \
    2021         classes/DelphesModule.h
    2022         @touch $@
    2023 
    20241904modules/Hector.h: \
    20251905        classes/DelphesModule.h
    … …  
    21212001
    21222002modules/FastJetFinder.h: \
    2123         classes/DelphesModule.h
    2124         @touch $@
    2125 
    2126 modules/BeamSpotFilter.h: \
    21272003        classes/DelphesModule.h
    21282004        @touch $@

  • cards/converter_card.tcl

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

  • cards/delphes_card_ATLAS.tcl

    r21eab4f rec5e04b  
    141141
    142142  # resolution formula for charged hadrons
    143   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    144                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    145                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     143  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
     144                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
     145                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
    146146}
    147147

  • cards/delphes_card_ATLAS_PileUp.tcl

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

  • cards/delphes_card_CMS.tcl

    r21eab4f rec5e04b  
    1515
    1616  TrackMerger
    17  
    18   ECal
    19   HCal
    20  
    2117  Calorimeter
    2218  EFlowMerger
    … …  
    128124  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
    129125                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
    130                                            (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e3) * (0.99) +
    131                                            (abs(eta) <= 1.5) * (pt > 1.0e3 )               * (0.99 * exp(0.5 - pt*5.0e-4)) +
    132 
     126                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
    133127                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
    134                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e3) * (0.98) +
    135                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e3)                * (0.98 * exp(0.5 - pt*5.0e-4)) +
     128                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
    136129                         (abs(eta) > 2.5)                                                  * (0.00)}
    137130}
    … …  
    149142  # resolution formula for charged hadrons
    150143  # based on arXiv:1405.6569
    151   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    152                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    153                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     144  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
     145                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
     146                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
    154147}
    155148
    … …  
    199192}
    200193
    201 
    202 
    203194#############
    204 #   ECAL
     195# Calorimeter
    205196#############
    206197
    207 module SimpleCalorimeter ECal {
     198module Calorimeter Calorimeter {
    208199  set ParticleInputArray ParticlePropagator/stableParticles
    209200  set TrackInputArray TrackMerger/tracks
    210201
    211   set TowerOutputArray ecalTowers
     202  set TowerOutputArray towers
     203  set PhotonOutputArray photons
     204
    212205  set EFlowTrackOutputArray eflowTracks
    213   set EFlowTowerOutputArray eflowPhotons
    214 
    215   set IsEcal true
    216 
    217   set EnergyMin 0.5
    218   set EnergySignificanceMin 2.0
    219 
    220   set SmearTowerCenter true
    221 
    222   set pi [expr {acos(-1)}]
    223 
    224   # lists of the edges of each tower in eta and phi
    225   # each list starts with the lower edge of the first tower
    226   # the list ends with the higher edged of the last tower
    227 
    228   # assume 0.02 x 0.02 resolution in eta,phi in the barrel |eta| < 1.5
    229 
    230   set PhiBins {}
    231   for {set i -180} {$i <= 180} {incr i} {
    232     add PhiBins [expr {$i * $pi/180.0}]
    233   }
    234 
    235   # 0.02 unit in eta up to eta = 1.5 (barrel)
    236   for {set i -85} {$i <= 86} {incr i} {
    237     set eta [expr {$i * 0.0174}]
    238     add EtaPhiBins $eta $PhiBins
    239   }
    240 
    241   # assume 0.02 x 0.02 resolution in eta,phi in the endcaps 1.5 < |eta| < 3.0 (HGCAL- ECAL)
    242 
    243   set PhiBins {}
    244   for {set i -180} {$i <= 180} {incr i} {
    245     add PhiBins [expr {$i * $pi/180.0}]
    246   }
    247 
    248   # 0.02 unit in eta up to eta = 3
    249   for {set i 1} {$i <= 84} {incr i} {
    250     set eta [expr { -2.958 + $i * 0.0174}]
    251     add EtaPhiBins $eta $PhiBins
    252   }
    253 
    254   for {set i 1} {$i <= 84} {incr i} {
    255     set eta [expr { 1.4964 + $i * 0.0174}]
    256     add EtaPhiBins $eta $PhiBins
    257   }
    258 
    259   # take present CMS granularity for HF
    260 
    261   # 0.175 x (0.175 - 0.35) resolution in eta,phi in the HF 3.0 < |eta| < 5.0
    262   set PhiBins {}
    263   for {set i -18} {$i <= 18} {incr i} {
    264     add PhiBins [expr {$i * $pi/18.0}]
    265   }
    266 
    267   foreach eta {-5 -4.7 -4.525 -4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.958 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525 4.7 5} {
    268     add EtaPhiBins $eta $PhiBins
    269   }
    270 
    271 
    272   add EnergyFraction {0} {0.0}
    273   # energy fractions for e, gamma and pi0
    274   add EnergyFraction {11} {1.0}
    275   add EnergyFraction {22} {1.0}
    276   add EnergyFraction {111} {1.0}
    277   # energy fractions for muon, neutrinos and neutralinos
    278   add EnergyFraction {12} {0.0}
    279   add EnergyFraction {13} {0.0}
    280   add EnergyFraction {14} {0.0}
    281   add EnergyFraction {16} {0.0}
    282   add EnergyFraction {1000022} {0.0}
    283   add EnergyFraction {1000023} {0.0}
    284   add EnergyFraction {1000025} {0.0}
    285   add EnergyFraction {1000035} {0.0}
    286   add EnergyFraction {1000045} {0.0}
    287   # energy fractions for K0short and Lambda
    288   add EnergyFraction {310} {0.3}
    289   add EnergyFraction {3122} {0.3}
    290 
    291   # set ResolutionFormula {resolution formula as a function of eta and energy}
    292 
    293   # for the ECAL barrel (|eta| < 1.5), see hep-ex/1306.2016 and 1502.02701
    294 
    295   # set ECalResolutionFormula {resolution formula as a function of eta and energy}
    296   # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
    297   set ResolutionFormula {                      (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
    298                              (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
    299                              (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
    300 
    301 }
    302 
    303 
    304 #############
    305 #   HCAL
    306 #############
    307 
    308 module SimpleCalorimeter HCal {
    309   set ParticleInputArray ParticlePropagator/stableParticles
    310   set TrackInputArray ECal/eflowTracks
    311 
    312   set TowerOutputArray hcalTowers
    313   set EFlowTrackOutputArray eflowTracks
    314   set EFlowTowerOutputArray eflowNeutralHadrons
    315 
    316   set IsEcal false
    317 
    318   set EnergyMin 1.0
    319   set EnergySignificanceMin 2.0
     206  set EFlowPhotonOutputArray eflowPhotons
     207  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
     208
     209  set ECalEnergyMin 0.5
     210  set HCalEnergyMin 1.0
     211
     212  set ECalEnergySignificanceMin 1.0
     213  set HCalEnergySignificanceMin 1.0
    320214
    321215  set SmearTowerCenter true
    … …  
    355249
    356250  # default energy fractions {abs(PDG code)} {Fecal Fhcal}
    357   add EnergyFraction {0} {1.0}
     251  add EnergyFraction {0} {0.0 1.0}
    358252  # energy fractions for e, gamma and pi0
    359   add EnergyFraction {11} {0.0}
    360   add EnergyFraction {22} {0.0}
    361   add EnergyFraction {111} {0.0}
     253  add EnergyFraction {11} {1.0 0.0}
     254  add EnergyFraction {22} {1.0 0.0}
     255  add EnergyFraction {111} {1.0 0.0}
    362256  # energy fractions for muon, neutrinos and neutralinos
    363   add EnergyFraction {12} {0.0}
    364   add EnergyFraction {13} {0.0}
    365   add EnergyFraction {14} {0.0}
    366   add EnergyFraction {16} {0.0}
    367   add EnergyFraction {1000022} {0.0}
    368   add EnergyFraction {1000023} {0.0}
    369   add EnergyFraction {1000025} {0.0}
    370   add EnergyFraction {1000035} {0.0}
    371   add EnergyFraction {1000045} {0.0}
     257  add EnergyFraction {12} {0.0 0.0}
     258  add EnergyFraction {13} {0.0 0.0}
     259  add EnergyFraction {14} {0.0 0.0}
     260  add EnergyFraction {16} {0.0 0.0}
     261  add EnergyFraction {1000022} {0.0 0.0}
     262  add EnergyFraction {1000023} {0.0 0.0}
     263  add EnergyFraction {1000025} {0.0 0.0}
     264  add EnergyFraction {1000035} {0.0 0.0}
     265  add EnergyFraction {1000045} {0.0 0.0}
    372266  # energy fractions for K0short and Lambda
    373   add EnergyFraction {310} {0.7}
    374   add EnergyFraction {3122} {0.7}
     267  add EnergyFraction {310} {0.3 0.7}
     268  add EnergyFraction {3122} {0.3 0.7}
     269
     270  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
     271  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
     272  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
     273                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
     274                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
    375275
    376276  # set HCalResolutionFormula {resolution formula as a function of eta and energy}
    377   set ResolutionFormula {                      (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
     277  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
    378278                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
    379 
    380 }
    381 
    382 
    383 #################
    384 # Electron filter
    385 #################
    386 
    387 module PdgCodeFilter ElectronFilter {
    388   set InputArray HCal/eflowTracks
    389   set OutputArray electrons
    390   set Invert true
    391   add PdgCode {11}
    392   add PdgCode {-11}
    393 }
    394 
    395 ###################################################
    396 # Tower Merger (in case not using e-flow algorithm)
    397 ###################################################
    398 
    399 module Merger Calorimeter {
    400 # add InputArray InputArray
    401   add InputArray ECal/ecalTowers
    402   add InputArray HCal/hcalTowers
    403   set OutputArray towers
    404 }
    405 
    406 
     279}
    407280
    408281####################
    … …  
    412285module Merger EFlowMerger {
    413286# add InputArray InputArray
    414   add InputArray HCal/eflowTracks
    415   add InputArray ECal/eflowPhotons
    416   add InputArray HCal/eflowNeutralHadrons
     287  add InputArray Calorimeter/eflowTracks
     288  add InputArray Calorimeter/eflowPhotons
     289  add InputArray Calorimeter/eflowNeutralHadrons
    417290  set OutputArray eflow
    418291}
    … …  
    423296
    424297module Efficiency PhotonEfficiency {
    425   set InputArray ECal/eflowPhotons
     298  set InputArray Calorimeter/eflowPhotons
    426299  set OutputArray photons
    427300
    … …  
    452325}
    453326
     327#################
     328# Electron filter
     329#################
     330
     331module PdgCodeFilter ElectronFilter {
     332  set InputArray Calorimeter/eflowTracks
     333  set OutputArray electrons
     334  set Invert true
     335  add PdgCode {11}
     336  add PdgCode {-11}
     337}
    454338
    455339#####################
    … …  
    498382
    499383  # efficiency formula for muons
    500   set EfficiencyFormula {                                     (pt <= 10.0)                * (0.00) +
    501                                            (abs(eta) <= 1.5) * (pt > 10.0)                * (0.95) +
    502                          (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0)                * (0.95) +
     384  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
     385                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
     386                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
     387                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
     388                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
    503389                         (abs(eta) > 2.4)                                                 * (0.00)}
    504390}
    … …  
    716602  add Branch Calorimeter/towers Tower Tower
    717603
    718   add Branch HCal/eflowTracks EFlowTrack Track
    719   add Branch ECal/eflowPhotons EFlowPhoton Tower
    720   add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
     604  add Branch Calorimeter/eflowTracks EFlowTrack Track
     605  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
     606  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
    721607
    722608  add Branch GenJetFinder/jets GenJet Jet
    723609  add Branch GenMissingET/momentum GenMissingET MissingET
    724  
     610
    725611  add Branch UniqueObjectFinder/jets Jet Jet
    726612  add Branch UniqueObjectFinder/electrons Electron Electron

  • cards/delphes_card_CMS_NoFastJet.tcl

    r21eab4f rec5e04b  
    128128  # resolution formula for electrons
    129129  # based on arXiv:1405.6569
    130   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
    131                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
    132                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
     130  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
     131                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
     132                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
    133133}
    134134
    … …  
    144144
    145145  # resolution formula for muons
    146   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
    147                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
    148                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
    149 }
    150 
     146  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
     147                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
     148                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
     149}
    151150
    152151##############

  • cards/delphes_card_CMS_PileUp.tcl

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

  • classes/DelphesClasses.cc

    r21eab4f rec5e04b  
    4141CompBase *Tower::fgCompare = CompE<Tower>::Instance();
    4242CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
    43 CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
    4443CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
    4544
    … …  
    122121  Charge(0), Mass(0.0),
    123122  IsPU(0), IsRecoPU(0), IsConstituent(0), IsFromConversion(0),
    124   ClusterIndex(-1), ClusterNDF(0), ClusterSigma(0), SumPT2(0), BTVSumPT2(0), GenDeltaZ(0), GenSumPT2(0),
    125123  Flavor(0), FlavorAlgo(0), FlavorPhys(0),
    126124  BTag(0), BTagAlgo(0), BTagPhys(0),
    … …  
    129127  Momentum(0.0, 0.0, 0.0, 0.0),
    130128  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),
    133129  Area(0.0, 0.0, 0.0, 0.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),
     130  Dxy(0), SDxy(0), Xd(0), Yd(0), Zd(0),
    142131  TrackResolution(0),
    143132  NCharged(0),
    … …  
    256245  object.IsConstituent = IsConstituent;
    257246  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;
    265247  object.Flavor = Flavor;
    266248  object.FlavorAlgo = FlavorAlgo;
    … …  
    280262  object.Momentum = Momentum;
    281263  object.Position = Position;
    282   object.InitialPosition = InitialPosition;
    283   object.PositionError = PositionError;
    284264  object.Area = Area;
    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; 
     265  object.Dxy = Dxy;
     266  object.SDxy = SDxy;
    299267  object.Xd = Xd;
    300268  object.Yd = Yd;
    … …  
    314282  object.SumPtChargedPU = SumPtChargedPU;
    315283  object.SumPt = SumPt;
    316   object.ClusterIndex = ClusterIndex;
    317   object.ClusterNDF = ClusterNDF;
    318   object.ClusterSigma = ClusterSigma;
    319   object.SumPT2 = SumPT2;
    320  
     284
    321285  object.FracPt[0] = FracPt[0];
    322286  object.FracPt[1] = FracPt[1];
    … …  
    399363  Momentum.SetXYZT(0.0, 0.0, 0.0, 0.0);
    400364  Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
    401   InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    402365  Area.SetXYZT(0.0, 0.0, 0.0, 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;
     366  Dxy = 0.0;
     367  SDxy = 0.0;
    417368  Xd = 0.0;
    418369  Yd = 0.0;
    … …  
    436387  SumPt = -999;
    437388
    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  
    446389  FracPt[0] = 0.0;
    447390  FracPt[1] = 0.0;

  • classes/DelphesClasses.h

    r21eab4f rec5e04b  
    147147  Float_t Pz; // particle momentum vector (z component) | hepevt.phep[number][2]
    148148
    149   Float_t D0;
    150   Float_t DZ;
    151   Float_t P;
    152   Float_t PT;
    153   Float_t CtgTheta;
    154   Float_t Phi;
     149  Float_t PT; // particle transverse momentum
    155150  Float_t Eta; // particle pseudorapidity
     151  Float_t Phi; // particle azimuthal angle
     152
    156153  Float_t Rapidity; // particle rapidity
    157154
    … …  
    171168//---------------------------------------------------------------------------
    172169
    173 class Vertex: public SortableObject
     170class Vertex: public TObject
    174171{
    175172public:
    … …  
    179176  Float_t Z; // vertex position (z component)
    180177
    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)
     178  ClassDef(Vertex, 1)
    200179};
    201180
    … …  
    418397  Int_t Charge; // track charge
    419398
     399  Float_t PT; // track transverse momentum
     400
    420401  Float_t Eta; // track pseudorapidity
     402  Float_t Phi; // track azimuthal angle
    421403
    422404  Float_t EtaOuter; // track pseudorapidity at the tracker edge
    … …  
    433415  Float_t TOuter; // track position (z component) at the tracker edge
    434416
    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 
     417  Float_t Dxy;     // track signed transverse impact parameter
     418  Float_t SDxy;    // signed error on the track signed transverse impact parameter
    456419  Float_t Xd;      // X coordinate of point of closest approach to vertex
    457420  Float_t Yd;      // Y coordinate of point of closest approach to vertex
    … …  
    460423  TRef Particle; // reference to generated particle
    461424
    462   Int_t VertexIndex; // reference to vertex
    463 
    464425  static CompBase *fgCompare; //!
    465426  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    565526  Float_t DeltaPhi;
    566527
    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 
     528  TLorentzVector Momentum, Position, Area;
     529
     530  Float_t Dxy;
     531  Float_t SDxy;
    584532  Float_t Xd;
    585533  Float_t Yd;
    … …  
    587535
    588536  // tracking resolution
    589 
     537 
    590538  Float_t TrackResolution;
    591539
    … …  
    614562  Float_t SumPt;
    615563
    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 
    626564  // N-subjettiness variables
    627565
    … …  
    657595  void SetFactory(DelphesFactory *factory) { fFactory = factory; }
    658596
    659   ClassDef(Candidate, 5)
     597  ClassDef(Candidate, 4)
    660598};
    661599

  • classes/SortableObject.h

    r21eab4f rec5e04b  
    156156      return -1;
    157157    else if(t1->ET < t2->ET)
    158       return 1;
    159     else
    160       return 0;
    161   }
    162 };
    163 
    164 //---------------------------------------------------------------------------
    165 
    166 template <typename T>
    167 class CompSumPT2: public CompBase
    168 {
    169   CompSumPT2() {}
    170 public:
    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)
    184158      return 1;
    185159    else

  • doc/genMakefile.tcl

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

  • examples/Pythia8/configParticleGun.cmnd

    r21eab4f rec5e04b  
    33! 1) Settings used in the main program.
    44
    5 Main:numberOfEvents = 10000        ! number of events to generate
     5Main:numberOfEvents = 10000         ! 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 = 11               ! 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 = 11               ! 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().

  • modules/Calorimeter.cc

    r21eab4f rec5e04b  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60  
     60  Int_t i;
     61
    6162  fECalResolutionFormula = new DelphesFormula;
    6263  fHCalResolutionFormula = new DelphesFormula;
    6364
    64   fECalTowerTrackArray = new TObjArray;
    65   fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
    66 
    67   fHCalTowerTrackArray = new TObjArray;
    68   fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
    69  
     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  }
    7073}
    7174
    … …  
    7477Calorimeter::~Calorimeter()
    7578{
    76  
     79  Int_t i;
     80
    7781  if(fECalResolutionFormula) delete fECalResolutionFormula;
    7882  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    7983
    80   if(fECalTowerTrackArray) delete fECalTowerTrackArray;
    81   if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
    82 
    83   if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
    84   if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
    85  
     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  }
    8692}
    8793
    … …  
    213219  Double_t ecalEnergy, hcalEnergy;
    214220  Double_t ecalSigma, hcalSigma;
    215   Double_t energyGuess;
    216221  Int_t pdgCode;
    217222
    … …  
    363368      fHCalTowerEnergy = 0.0;
    364369
    365       fECalTrackEnergy = 0.0;
    366       fHCalTrackEnergy = 0.0;
    367      
    368       fECalTrackSigma = 0.0;
    369       fHCalTrackSigma = 0.0;
    370      
     370      fECalTrackEnergy[0] = 0.0;
     371      fECalTrackEnergy[1] = 0.0;
     372
     373      fHCalTrackEnergy[0] = 0.0;
     374      fHCalTrackEnergy[1] = 0.0;
     375
    371376      fTowerTrackHits = 0;
    372377      fTowerPhotonHits = 0;
    373378
    374       fECalTowerTrackArray->Clear();
    375       fHCalTowerTrackArray->Clear();
    376    
     379      fECalTowerTrackArray[0]->Clear();
     380      fECalTowerTrackArray[1]->Clear();
     381
     382      fHCalTowerTrackArray[0]->Clear();
     383      fHCalTowerTrackArray[1]->Clear();
    377384    }
    378385
    … …  
    399406      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    400407      {
    401         fECalTrackEnergy += ecalEnergy;
    402         ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());       
    403         if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;       
    404         else energyGuess = momentum.E();
    405 
    406         fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
    407         fECalTowerTrackArray->Add(track);
     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        }
    408419      }
    409      
    410420      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    411421      {
    412         fHCalTrackEnergy += hcalEnergy;
    413422        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    414         if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
    415         else energyGuess = momentum.E();
    416 
    417         fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
    418         fHCalTowerTrackArray->Add(track);
     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        }
    419433      }
    420      
    421434      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    422435      {
    … …  
    463476  Double_t energy, pt, eta, phi;
    464477  Double_t ecalEnergy, hcalEnergy;
    465   Double_t ecalNeutralEnergy, hcalNeutralEnergy;
    466  
    467478  Double_t ecalSigma, hcalSigma;
    468   Double_t ecalNeutralSigma, hcalNeutralSigma;
    469 
    470   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    471  
     479
    472480  TLorentzVector momentum;
    473481  TFractionMap::iterator itFractionMap;
    … …  
    476484
    477485  if(!fTower) return;
     486
    478487
    479488  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    545554    fTowerOutputArray->Add(fTower);
    546555  }
    547  
     556
    548557  // fill energy flow candidates
    549   fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
    550   fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
    551 
    552   //compute neutral excesses
    553   ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
    554   hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
    555  
    556   ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
    557   hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
    558  
    559    // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
    560   if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
     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)
    561618  {
    562619    // create new photon tower
    563620    tower = static_cast<Candidate*>(fTower->Clone());
    564     pt =  ecalNeutralEnergy / TMath::CosH(eta);
    565    
    566     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
    567     tower->Eem = ecalNeutralEnergy;
     621
     622    pt = ecalEnergy / TMath::CosH(eta);
     623
     624    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
     625    tower->Eem = ecalEnergy;
    568626    tower->Ehad = 0.0;
    569627    tower->PID = 22;
    570    
     628
    571629    fEFlowPhotonOutputArray->Add(tower);
    572    
    573     //clone tracks
    574     fItECalTowerTrackArray->Reset();
    575     while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
    576     {
    577       mother = track;
    578       track = static_cast<Candidate*>(track->Clone());
    579       track->AddCandidate(mother);
    580 
    581       fEFlowTrackOutputArray->Add(track);
    582     }
    583  
    584   }
    585  
    586   // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    587   else if(fECalTrackEnergy > 0.0)
    588   {
    589     weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
    590     weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
    591  
    592     bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
    593     rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
    594 
    595     //rescale tracks
    596     fItECalTowerTrackArray->Reset();
    597     while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
    598     { 
    599       mother = track;
    600       track = static_cast<Candidate*>(track->Clone());
    601       track->AddCandidate(mother);
    602 
    603       track->Momentum *= rescaleFactor;
    604 
    605       fEFlowTrackOutputArray->Add(track);
    606     }
    607   }
    608 
    609 
    610   // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
    611   if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
    612   {
    613     // create new photon tower
     630  }
     631  if(hcalEnergy > 0.0)
     632  {
     633    // create new neutral hadron tower
    614634    tower = static_cast<Candidate*>(fTower->Clone());
    615     pt =  hcalNeutralEnergy / TMath::CosH(eta);
    616    
    617     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
    618     tower->Ehad = hcalNeutralEnergy;
     635
     636    pt = hcalEnergy / TMath::CosH(eta);
     637
     638    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
    619639    tower->Eem = 0.0;
    620    
     640    tower->Ehad = hcalEnergy;
     641
    621642    fEFlowNeutralHadronOutputArray->Add(tower);
    622    
    623     //clone tracks
    624     fItHCalTowerTrackArray->Reset();
    625     while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
    626     {
    627       mother = track;
    628       track = static_cast<Candidate*>(track->Clone());
    629       track->AddCandidate(mother);
    630 
    631       fEFlowTrackOutputArray->Add(track);
    632     }
    633  
    634   }
    635  
    636   // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    637   else if(fHCalTrackEnergy > 0.0)
    638   {
    639     weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
    640     weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
    641  
    642     bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
    643     rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
    644 
    645     //rescale tracks
    646     fItHCalTowerTrackArray->Reset();
    647     while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
    648     { 
    649       mother = track;
    650       track = static_cast<Candidate*>(track->Clone());
    651       track->AddCandidate(mother);
    652 
    653       track->Momentum *= rescaleFactor;
    654 
    655       fEFlowTrackOutputArray->Add(track);
    656     }
    657   }
    658  
    659  
     643  }
    660644}
    661645

  • modules/Calorimeter.h

    r21eab4f rec5e04b  
    5858  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    5959  Double_t fECalTowerEnergy, fHCalTowerEnergy;
    60   Double_t fECalTrackEnergy, fHCalTrackEnergy;
     60  Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
    6161
    6262  Double_t fTimingEnergyMin;
    … …  
    7070  Double_t fECalEnergySignificanceMin;
    7171  Double_t fHCalEnergySignificanceMin;
    72 
    73   Double_t fECalTrackSigma;
    74   Double_t fHCalTrackSigma;
    7572
    7673  Bool_t fSmearTowerCenter;
    … …  
    106103  TObjArray *fEFlowNeutralHadronOutputArray; //!
    107104
    108   TObjArray *fECalTowerTrackArray; //!
    109   TIterator *fItECalTowerTrackArray; //!
     105  TObjArray *fECalTowerTrackArray[2]; //!
     106  TIterator *fItECalTowerTrackArray[2]; //!
    110107
    111   TObjArray *fHCalTowerTrackArray; //!
    112   TIterator *fItHCalTowerTrackArray; //!
     108  TObjArray *fHCalTowerTrackArray[2]; //!
     109  TIterator *fItHCalTowerTrackArray[2]; //!
    113110
    114111  void FinalizeTower();

  • modules/ImpactParameterSmearing.cc

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

  • modules/ModulesLinkDef.h

    r21eab4f rec5e04b  
    3535#include "modules/EnergySmearing.h"
    3636#include "modules/MomentumSmearing.h"
    37 #include "modules/TrackSmearing.h"
    3837#include "modules/ImpactParameterSmearing.h"
    3938#include "modules/TimeSmearing.h"
    … …  
    5958#include "modules/StatusPidFilter.h"
    6059#include "modules/PdgCodeFilter.h"
    61 #include "modules/BeamSpotFilter.h"
    6260#include "modules/RecoPuFilter.h"
    6361#include "modules/Cloner.h"
    … …  
    6664#include "modules/JetFlavorAssociation.h"
    6765#include "modules/JetFakeParticle.h"
    68 #include "modules/VertexSorter.h"
    69 #include "modules/VertexFinder.h"
    70 #include "modules/VertexFinderDA4D.h"
    7166#include "modules/ExampleModule.h"
    7267
    … …  
    8681#pragma link C++ class EnergySmearing+;
    8782#pragma link C++ class MomentumSmearing+;
    88 #pragma link C++ class TrackSmearing+;
    8983#pragma link C++ class ImpactParameterSmearing+;
    9084#pragma link C++ class TimeSmearing+;
    … …  
    110104#pragma link C++ class StatusPidFilter+;
    111105#pragma link C++ class PdgCodeFilter+;
    112 #pragma link C++ class BeamSpotFilter+;
    113106#pragma link C++ class RecoPuFilter+;
    114107#pragma link C++ class Cloner+;
    … …  
    117110#pragma link C++ class JetFlavorAssociation+;
    118111#pragma link C++ class JetFakeParticle+;
    119 #pragma link C++ class VertexSorter+;
    120 #pragma link C++ class VertexFinder+;
    121 #pragma link C++ class VertexFinderDA4D+;
    122112#pragma link C++ class ExampleModule+;
    123113

  • modules/MomentumSmearing.cc

    r21eab4f rec5e04b  
    9696{
    9797  Candidate *candidate, *mother;
    98   Double_t pt, eta, phi, e, res;
     98  Double_t pt, eta, phi, e;
    9999
    100100  fItInputArray->Reset();
    … …  
    107107    pt = candidateMomentum.Pt();
    108108    e = candidateMomentum.E();
    109     res = fFormula->Eval(pt, eta, phi, e);
    110  
     109
    111110    // apply smearing formula
    112     //pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     111    pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
    113112   
    114     res = ( res > 1.0 ) ? 1.0 : res;
    115 
    116     pt = LogNormal(pt, res * pt );
    117    
    118     //if(pt <= 0.0) continue;
     113    if(pt <= 0.0) continue;
    119114
    120115    mother = candidate;
    … …  
    123118    phi = candidateMomentum.Phi();
    124119    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt*TMath::CosH(eta));
    125     //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
    126     candidate->TrackResolution = res;
     120    candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
    127121    candidate->AddCandidate(mother);
    128122       
    … …  
    130124  }
    131125}
    132 //----------------------------------------------------------------
    133 
    134 Double_t MomentumSmearing::LogNormal(Double_t mean, Double_t sigma)
    135 {
    136   Double_t a, b;
    137 
    138   if(mean > 0.0)
    139   {
    140     b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
    141     a = TMath::Log(mean) - 0.5*b*b;
    142 
    143     return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
    144   }
    145   else
    146   {
    147     return 0.0;
    148   }
    149 }
    150 
    151126
    152127//------------------------------------------------------------------------------

  • modules/MomentumSmearing.h

    r21eab4f rec5e04b  
    4747private:
    4848
    49   Double_t LogNormal(Double_t mean, Double_t sigma);
    50 
    5149  DelphesFormula *fFormula; //!
    5250

  • modules/ParticlePropagator.cc

    r21eab4f rec5e04b  
    6666{
    6767}
    68 
    6968
    7069//------------------------------------------------------------------------------
    … …  
    9291  fItInputArray = fInputArray->MakeIterator();
    9392
    94   // import beamspot
    95   try
    96   {
    97     fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
    98   }
    99   catch(runtime_error &e)
    100   {
    101     fBeamSpotInputArray = 0;
    102   }
    10393  // create output arrays
    10494
    … …  
    121111{
    122112  Candidate *candidate, *mother;
    123   TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
     113  TLorentzVector candidatePosition, candidateMomentum;
    124114  Double_t px, py, pz, pt, pt2, e, q;
    125115  Double_t x, y, z, t, r, phi;
    … …  
    130120  Double_t tmp, discr, discr2;
    131121  Double_t delta, gammam, omega, asinrho;
    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;
     122  Double_t rcu, rc2, dxy, xd, yd, zd;
    135123
    136124  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   }
    145125
    146126  fItInputArray->Reset();
    … …  
    152132    y = candidatePosition.Y()*1.0E-3;
    153133    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 
    159134    q = candidate->Charge;
    160135
    … …  
    207182      z_t = z + pz*t;
    208183
    209       l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
    210 
    211184      mother = candidate;
    212185      candidate = static_cast<Candidate*>(candidate->Clone());
    213186
    214       candidate->InitialPosition = candidatePosition;
    215187      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;
    217188
    218189      candidate->Momentum = candidateMomentum;
    … …  
    268239      zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
    269240
    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 
     241      // calculate impact paramater
     242      dxy = (xd*py - yd*px)/pt;
    287243
    288244      // 3. time evaluation t = TMath::Min(t_r, t_z)
    … …  
    331287      r_t = TMath::Hypot(x_t, y_t);
    332288
    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 
    339289      if(r_t > 0.0)
    340290      {
    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 
    350291        mother = candidate;
    351292        candidate = static_cast<Candidate*>(candidate->Clone());
    352293
    353         candidate->InitialPosition = candidatePosition;
    354294        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
    355295
    356296        candidate->Momentum = candidateMomentum;
    357 
    358             candidate->L  =  l*1.0E3;
    359 
    360             candidate->Xd = xd*1.0E3;
     297        candidate->Dxy = dxy*1.0E3;
     298        candidate->Xd = xd*1.0E3;
    361299        candidate->Yd = yd*1.0E3;
    362300        candidate->Zd = zd*1.0E3;

  • modules/ParticlePropagator.h

    r21eab4f rec5e04b  
    3535class TClonesArray;
    3636class TIterator;
    37 class TLorentzVector;
    3837
    3938class ParticlePropagator: public DelphesModule
    … …  
    5655
    5756  const TObjArray *fInputArray; //!
    58   const TObjArray *fBeamSpotInputArray; //!
    5957
    6058  TObjArray *fOutputArray; //!

  • modules/PileUpMerger.cc

    r21eab4f rec5e04b  
    115115  TDatabasePDG *pdg = TDatabasePDG::Instance();
    116116  TParticlePDG *pdgParticle;
    117   Int_t pid, nch, nvtx = -1;
     117  Int_t pid;
    118118  Float_t x, y, z, t, vx, vy;
    119   Float_t px, py, pz, e, pt;
    120   Double_t dz, dphi, dt, sumpt2, dz0, dt0;
     119  Float_t px, py, pz, e;
     120  Double_t dz, dphi, dt;
    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 
    137134  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    138135  dz *= 1.0E3; // necessary in order to make z in mm
    139 
    140   //cout<<dz<<","<<dt<<endl;
    141 
    142136  vx = 0.0;
    143137  vy = 0.0;
    144 
    145138  numberOfParticles = fInputArray->GetEntriesFast();
    146   nch = 0;
    147   sumpt2 = 0.0;
    148 
    149   factory = GetFactory();
    150   vertex = factory->NewCandidate();
    151 
    152139  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    153140  {
    … …  
    156143    z = candidate->Position.Z();
    157144    t = candidate->Position.T();
    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 
     145    candidate->Position.SetZ(z + dz);
     146    candidate->Position.SetT(t + dt);
    172147    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     }
    180148  }
    181149
    182150  if(numberOfParticles > 0)
    183151  {
    184     vx /= sumpt2;
    185     vy /= sumpt2;
    186   }
    187 
    188   nvtx++;
     152    vx /= numberOfParticles;
     153    vy /= numberOfParticles;
     154  }
     155
     156  factory = GetFactory();
     157
     158  vertex = factory->NewCandidate();
    189159  vertex->Position.SetXYZT(vx, vy, dz, dt);
    190   vertex->ClusterIndex = nvtx;
    191   vertex->ClusterNDF = nch;
    192   vertex->SumPT2 = sumpt2;
    193   vertex->GenSumPT2 = sumpt2;
    194160  fVertexOutputArray->Add(vertex);
    195161
    … …  
    204170      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
    205171      break;
    206     case 2:
    207       numberOfEvents = fMeanPileUp;
    208       break;
    209172    default:
    210173      numberOfEvents = gRandom->Poisson(fMeanPileUp);
    … …  
    213176
    214177  allEntries = fReader->GetEntries();
    215 
    216178
    217179  for(event = 0; event < numberOfEvents; ++event)
    … …  
    236198    vx = 0.0;
    237199    vy = 0.0;
    238 
    239200    numberOfParticles = 0;
    240     sumpt2 = 0.0;
    241 
    242     //factory = GetFactory();
    243     vertex = factory->NewCandidate();
    244 
    245201    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    246202    {
    … …  
    259215      candidate->Momentum.SetPxPyPzE(px, py, pz, e);
    260216      candidate->Momentum.RotateZ(dphi);
    261       pt = candidate->Momentum.Pt();
    262217
    263218      x -= fInputBeamSpotX;
    … …  
    269224      vx += candidate->Position.X();
    270225      vy += candidate->Position.Y();
    271 
    272226      ++numberOfParticles;
    273       if(TMath::Abs(candidate->Charge) >  1.0E-9)
    274       {
    275         nch++;
    276         sumpt2 += pt*pt;
    277         vertex->AddCandidate(candidate);
    278       }
    279227
    280228      fParticleOutputArray->Add(candidate);
    … …  
    287235    }
    288236
    289     nvtx++;
    290 
     237    vertex = factory->NewCandidate();
    291238    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 
    298239    vertex->IsPU = 1;
    299240
    300241    fVertexOutputArray->Add(vertex);
    301 
    302   }
    303 }
    304 
    305 //------------------------------------------------------------------------------
     242  }
     243}
     244
     245//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.cc

    r21eab4f rec5e04b  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60  
     60  Int_t i;
     61
    6162  fResolutionFormula = new DelphesFormula;
    62   fTowerTrackArray = new TObjArray;
    63   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    64  
     63
     64  for(i = 0; i < 2; ++i)
     65  {
     66    fTowerTrackArray[i] = new TObjArray;
     67    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
     68  }
    6569}
    6670
    … …  
    6973SimpleCalorimeter::~SimpleCalorimeter()
    7074{
    71  
     75  Int_t i;
     76
    7277  if(fResolutionFormula) delete fResolutionFormula;
    73   if(fTowerTrackArray) delete fTowerTrackArray;
    74   if(fItTowerTrackArray) delete fItTowerTrackArray;
    75  
     78
     79  for(i = 0; i < 2; ++i)
     80  {
     81    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
     82    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
     83  }
    7684}
    7785
    … …  
    191199  Double_t energy;
    192200  Double_t sigma;
    193   Double_t energyGuess;
    194 
    195201  Int_t pdgCode;
    196202
    … …  
    334340      fTowerEnergy = 0.0;
    335341
    336       fTrackEnergy = 0.0;
    337       fTrackSigma = 0.0;
     342      fTrackEnergy[0] = 0.0;
     343      fTrackEnergy[1] = 0.0;
    338344
    339345      fTowerTime = 0.0;
    … …  
    345351      fTowerPhotonHits = 0;
    346352
    347       fTowerTrackArray->Clear();
    348      }
     353      fTowerTrackArray[0]->Clear();
     354      fTowerTrackArray[1]->Clear();
     355    }
    349356
    350357    // check for track hits
    … …  
    364371      if(fTrackFractions[number] > 1.0E-9)
    365372      {
    366              
    367        // compute total charged energy   
    368        fTrackEnergy += energy;
    369        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    370        if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy;
    371        else energyGuess = momentum.E();
    372 
    373        fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess);
    374        fTowerTrackArray->Add(track);
     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        }
    375384      }
    376        
    377385      else
    378386      {
    … …  
    395403    fTowerEnergy += energy;
    396404
    397     fTowerTime += energy*position.T();
    398     fTowerTimeWeight += energy;
     405    fTowerTime += TMath::Sqrt(energy)*position.T();
     406    fTowerTimeWeight += TMath::Sqrt(energy);
    399407
    400408    fTower->AddCandidate(particle);
    … …  
    410418{
    411419  Candidate *tower, *track, *mother;
    412   Double_t energy,neutralEnergy, pt, eta, phi;
    413   Double_t sigma, neutralSigma;
     420  Double_t energy, pt, eta, phi;
     421  Double_t sigma;
    414422  Double_t time;
    415    
    416   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    417423
    418424  TLorentzVector momentum;
    … …  
    430436
    431437  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    432 
    433438
    434439  if(fSmearTowerCenter)
    … …  
    459464  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    460465
    461  
    462   // e-flow candidates
    463 
    464   //compute neutral excess
    465  
    466   fTrackSigma = TMath::Sqrt(fTrackSigma);
    467   neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
    468  
    469   //compute sigma_trk total
    470   neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
    471    
    472   // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
    473   if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
     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)
    474499  {
    475500    // create new photon tower
    476501    tower = static_cast<Candidate*>(fTower->Clone());
    477     pt = neutralEnergy / TMath::CosH(eta);
    478 
    479     tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
    480     tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     502    pt = energy / TMath::CosH(eta);
     503
     504    tower->Eem = (!fIsEcal) ? 0 : energy;
     505    tower->Ehad = (fIsEcal) ? 0 : energy;
     506   
     507    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     508   
    481509    tower->PID = (fIsEcal) ? 22 : 0;
    482  
    483     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     510   
    484511    fEFlowTowerOutputArray->Add(tower);
    485    
    486     fItTowerTrackArray->Reset();
    487     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    488     {
    489       mother = track;
    490       track = static_cast<Candidate*>(track->Clone());
    491       track->AddCandidate(mother);
    492 
    493       fEFlowTrackOutputArray->Add(track);
    494     }
    495   }
    496    
    497   // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    498   else if (fTrackEnergy > 0.0)
    499   {
    500     weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
    501     weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
    502  
    503     bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
    504     rescaleFactor = bestEnergyEstimate/fTrackEnergy;
    505    
    506     fItTowerTrackArray->Reset();
    507     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    508     {
    509       mother = track;
    510       track = static_cast<Candidate*>(track->Clone());
    511       track->AddCandidate(mother);
    512 
    513       track->Momentum *= rescaleFactor;
    514 
    515       fEFlowTrackOutputArray->Add(track);
    516     }
    517   }
    518    
    519  }
     512  }
     513}
    520514
    521515//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

    r21eab4f rec5e04b  
    5959  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    6060  Double_t fTowerEnergy;
    61   Double_t fTrackEnergy;
     61  Double_t fTrackEnergy[2];
    6262
    6363  Double_t fTowerTime;
    … …  
    7272
    7373  Double_t fEnergySignificanceMin;
    74 
    75   Double_t fTrackSigma;
    7674
    7775  Bool_t fSmearTowerCenter;
    … …  
    104102  TObjArray *fEFlowTowerOutputArray; //!
    105103
    106   TObjArray *fTowerTrackArray; //!
    107   TIterator *fItTowerTrackArray; //!
     104  TObjArray *fTowerTrackArray[2]; //!
     105  TIterator *fItTowerTrackArray[2]; //!
    108106
    109107  void FinalizeTower();

  • modules/TimeSmearing.cc

    r21eab4f rec5e04b  
    9393{
    9494  Candidate *candidate, *mother;
    95   Double_t ti, tf_smeared, tf;
     95  Double_t t;
    9696  const Double_t c_light = 2.99792458E8;
    9797
    … …  
    9999  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    100100  {
    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;
     101    const TLorentzVector &candidatePosition = candidate->Position;
     102    t = candidatePosition.T()*1.0E-3/c_light;
    106103
    107104    // apply smearing formula
    108     tf_smeared = gRandom->Gaus(tf, fTimeResolution);
    109     ti = ti + tf_smeared - tf;
     105    t = gRandom->Gaus(t, fTimeResolution);
    110106
    111107    mother = candidate;
    112108    candidate = static_cast<Candidate*>(candidate->Clone());
    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;
     109    candidate->Position.SetT(t*1.0E3*c_light);
    117110
    118111    candidate->AddCandidate(mother);

  • modules/TrackCountingBTagging.cc

    r21eab4f rec5e04b  
    9696
    9797  Double_t jpx, jpy;
    98   Double_t dr, tpt;
    99   Double_t xd, yd, d0, dd0, ip, sip;
     98  Double_t dr, tpx, tpy, tpt;
     99  Double_t xd, yd, dxy, ddxy, ip, sip;
    100100
    101101  Int_t sign;
    … …  
    117117    {
    118118      const TLorentzVector &trkMomentum = track->Momentum;
    119      
     119
    120120      dr = jetMomentum.DeltaR(trkMomentum);
     121
    121122      tpt = trkMomentum.Pt();
     123      tpx = trkMomentum.Px();
     124      tpy = trkMomentum.Py();
     125
    122126      xd = track->Xd;
    123127      yd = track->Yd;
    124       d0 = TMath::Hypot(xd, yd);
    125       dd0 = track->ErrorD0;
     128      dxy = TMath::Hypot(xd, yd);
     129      ddxy = track->SDxy;
    126130
    127131      if(tpt < fPtMin) continue;
    128132      if(dr > fDeltaR) continue;
    129       if(d0 > fIPmax) continue;
     133      if(dxy > fIPmax) continue;
    130134
    131135      sign = (jpx*xd + jpy*yd > 0.0) ? 1 : -1;
    132136
    133       ip = sign*d0;
    134       sip = ip / TMath::Abs(dd0);
     137      ip = sign*dxy;
     138      sip = ip / TMath::Abs(ddxy);
    135139
    136140      if(sip > fSigMin) count++;

  • modules/TreeWriter.cc

    r21eab4f rec5e04b  
    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 
    224217    entry->Eta = eta;
    225218    entry->Phi = momentum.Phi();
    … …  
    240233{
    241234  TIter iterator(array);
    242   Candidate *candidate = 0, *constituent = 0;
     235  Candidate *candidate = 0;
    243236  Vertex *entry = 0;
    244237
    245238  const Double_t c_light = 2.99792458E8;
    246239
    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 
    255240  // loop over all vertices
    256241  iterator.Reset();
    257242  while((candidate = static_cast<Candidate*>(iterator.Next())))
    258243  {
    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;
     244    const TLorentzVector &position = candidate->Position;
    277245
    278246    entry = static_cast<Vertex*>(branch->NewEntry());
    279247
    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 
     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}
    310254
    311255//------------------------------------------------------------------------------
    … …  
    317261  Candidate *particle = 0;
    318262  Track *entry = 0;
    319   Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
     263  Double_t pt, signz, cosTheta, eta, rapidity;
    320264  const Double_t c_light = 2.99792458E8;
    321265
    … …  
    348292    entry->TOuter = position.T()*1.0E-3/c_light;
    349293
    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 
     294    entry->Dxy = candidate->Dxy;
     295    entry->SDxy = candidate->SDxy ;
    362296    entry->Xd = candidate->Xd;
    363297    entry->Yd = candidate->Yd;
    … …  
    367301
    368302    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 
    373303    cosTheta = TMath::Abs(momentum.CosTheta());
    374304    signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
    … …  
    376306    rapidity = (cosTheta == 1.0 ? signz*999.9 : momentum.Rapidity());
    377307
    378     entry->PT  = pt;
    379308    entry->Eta = eta;
    380     entry->Phi = phi;
    381     entry->CtgTheta = ctgTheta;
     309    entry->Phi = momentum.Phi();
     310    entry->PT = pt;
    382311
    383312    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    … …  
    390319
    391320    entry->Particle = particle;
    392 
    393     entry->VertexIndex = candidate->ClusterIndex;
    394 
    395321  }
    396322}

  • readers/DelphesCMSFWLite.cpp

    r21eab4f rec5e04b  
    272272
    273273    branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
    274     branchRwgt = treeWriter->NewBranch("Weight", Weight::Class());
     274    branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
    275275
    276276    confReader = new ExRootConfReader;

  • readers/DelphesPythia8.cpp

    r21eab4f rec5e04b  
    153153// from pythia8 example 21
    154154
    155 void fillParticle(int id, double pt_max, double eta_max,
    156   Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
    157 {
     155void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
     156  Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
     157
    158158  // Reset event record to allow for new event.
    159159  event.reset();
    160160
    161   // Generate uniform pt and eta.
    162   double pt, eta, phi, pp, ee, mm;
    163   pt = pow(10, 1.0 + (log10(pt_max) - 1.0) * rndm.flat());
    164   eta = (2.0 * rndm.flat() - 1.0) * eta_max;
    165   phi = 2.0 * M_PI * rndm.flat();
    166   pp = pt * cosh(eta);
    167   mm = pdt.mSel(id);
    168   ee = Pythia8::sqrtpos(pp*pp + mm*mm);
     161  // Angles uniform in solid angle.
     162  double cThe, sThe, phi, ee;
     163  cThe = 2. * rndm.flat() - 1.;
     164  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     165  phi = 2. * M_PI * rndm.flat();
     166  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     167  double mm = pdt.mSel(id);
     168  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
    169169
    170170  // Store the particle in the event record.
    171   event.append(id, 1, 0, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
     171  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     172
    172173}
    173174
    174 void fillPartons(int id, double pt_max, double eta_max,
    175   Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
    176 {
     175void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
     176  Pythia8::Rndm& rndm) {
    177177
    178178  // Reset event record to allow for new event.
    179179  event.reset();
    180180
     181  // Angles uniform in solid angle.
     182  double cThe, sThe, phi, ee;
     183
    181184  // Information on a q qbar system, to be hadronized.
    182185
    183   // Generate uniform pt and eta.
    184   double pt, eta, phi, pp, ee, mm;
    185   pt = pow(10, 1.0 + (log10(pt_max) - 1.0) * rndm.flat());
    186   eta = (2.0 * rndm.flat() - 1.0) * eta_max;
    187   phi = 2.0 * M_PI * rndm.flat();
    188   pp = pt * cosh(eta);
    189   mm = pdt.m0(id);
    190   ee = Pythia8::sqrtpos(pp*pp + mm*mm);
    191 
    192   if(id == 21)
    193   {
    194     event.append(21, 23, 101, 102, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee);
    195     event.append(21, 23, 102, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee);
     186  cThe = 2. * rndm.flat() - 1.;
     187  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     188  phi = 2. * M_PI * rndm.flat();
     189  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     190  double mm = pdt.m0(type);
     191  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     192  if (type == 21)
     193  {
     194    event.append( 21, 23, 101, 102, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee);
     195    event.append( 21, 23, 102, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee);
    196196  }
    197197  else
    198198  {
    199     event.append(id, 23, 101, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
    200     event.append(-id, 23, 0, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee, mm);
     199    event.append(  type, 23, 101,   0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     200    event.append( -type, 23,   0, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee, mm);
    201201  }
    202202}
    … …  
    338338      if (pythia->flag("Main:spareFlag1"))
    339339      {
    340         if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 15 || pythia->mode("Main:spareMode1") == 22)
     340        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 22)
    341341        {
    342           fillParticle(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
     342          fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
    343343        }
    344         else
    345         {
    346           fillPartons(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
    347         }
     344        else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
    348345      }
    349346
Note: See TracChangeset for help on using the changeset viewer.