Fork me on GitHub

Changeset 21eab4f in git


Ignore:
Timestamp:
Aug 26, 2016, 5:17:14 PM (8 years ago)
Author:
GitHub <noreply@…>
Parents:
ec5e04b (diff), 94cacb6 (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.
git-author:
AlexandreMertens <alexandre.mertens@…> (08/26/16 17:17:14)
git-committer:
GitHub <noreply@…> (08/26/16 17:17:14)
Message:

Merge 94cacb6ab322f96b75a06369dfd9b7389d0d6802 into ec5e04b014990ea45a8708f4f726b6b214cdb82b

Files:
14 added
28 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    rec5e04b r21eab4f  
    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
     145Validation$(ExeSuf): \
     146        tmp/examples/Validation.$(ObjSuf)
     147
     148tmp/examples/Validation.$(ObjSuf): \
     149        examples/Validation.cpp \
    139150        classes/DelphesClasses.h \
    140151        external/ExRootAnalysis/ExRootTreeReader.h \
    … …  
    150161        root2pileup$(ExeSuf) \
    151162        stdhep2pileup$(ExeSuf) \
    152         Example1$(ExeSuf)
     163        Example1$(ExeSuf) \
     164        Validation$(ExeSuf)
    153165
    154166EXECUTABLE_OBJ +=  \
    … …  
    159171        tmp/converters/root2pileup.$(ObjSuf) \
    160172        tmp/converters/stdhep2pileup.$(ObjSuf) \
    161         tmp/examples/Example1.$(ObjSuf)
     173        tmp/examples/Example1.$(ObjSuf) \
     174        tmp/examples/Validation.$(ObjSuf)
    162175
    163176DelphesHepMC$(ExeSuf): \
    … …  
    183196        classes/DelphesLHEFReader.h \
    184197        external/ExRootAnalysis/ExRootTreeWriter.h \
     198        external/ExRootAnalysis/ExRootTreeBranch.h \
     199        external/ExRootAnalysis/ExRootProgressBar.h
     200DelphesROOT$(ExeSuf): \
     201        tmp/readers/DelphesROOT.$(ObjSuf)
     202
     203tmp/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 \
    185211        external/ExRootAnalysis/ExRootTreeBranch.h \
    186212        external/ExRootAnalysis/ExRootProgressBar.h
    … …  
    200226        DelphesHepMC$(ExeSuf) \
    201227        DelphesLHEF$(ExeSuf) \
     228        DelphesROOT$(ExeSuf) \
    202229        DelphesSTDHEP$(ExeSuf)
    203230
    … …  
    205232        tmp/readers/DelphesHepMC.$(ObjSuf) \
    206233        tmp/readers/DelphesLHEF.$(ObjSuf) \
     234        tmp/readers/DelphesROOT.$(ObjSuf) \
    207235        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    208236
    … …  
    319347        modules/EnergySmearing.h \
    320348        modules/MomentumSmearing.h \
     349        modules/TrackSmearing.h \
    321350        modules/ImpactParameterSmearing.h \
    322351        modules/TimeSmearing.h \
    … …  
    342371        modules/StatusPidFilter.h \
    343372        modules/PdgCodeFilter.h \
     373        modules/BeamSpotFilter.h \
     374        modules/RecoPuFilter.h \
    344375        modules/Cloner.h \
    345376        modules/Weighter.h \
    … …  
    347378        modules/JetFlavorAssociation.h \
    348379        modules/JetFakeParticle.h \
     380        modules/VertexSorter.h \
     381        modules/VertexFinder.h \
     382        modules/VertexFinderDA4D.h \
    349383        modules/ExampleModule.h
    350384tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    553587        classes/DelphesFactory.h \
    554588        classes/DelphesFormula.h
     589tmp/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
    555598tmp/modules/Calorimeter.$(ObjSuf): \
    556599        modules/Calorimeter.$(SrcSuf) \
    … …  
    785828        external/ExRootAnalysis/ExRootFilter.h \
    786829        external/ExRootAnalysis/ExRootClassifier.h
     830tmp/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
    787839tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    788840        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    846898        modules/TrackPileUpSubtractor.$(SrcSuf) \
    847899        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
     906tmp/modules/TrackSmearing.$(ObjSuf): \
     907        modules/TrackSmearing.$(SrcSuf) \
     908        modules/TrackSmearing.h \
    848909        classes/DelphesClasses.h \
    849910        classes/DelphesFactory.h \
    … …  
    868929        classes/DelphesFactory.h \
    869930        classes/DelphesFormula.h \
     931        external/ExRootAnalysis/ExRootResult.h \
     932        external/ExRootAnalysis/ExRootFilter.h \
     933        external/ExRootAnalysis/ExRootClassifier.h
     934tmp/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
     944tmp/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
     954tmp/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 \
    870961        external/ExRootAnalysis/ExRootResult.h \
    871962        external/ExRootAnalysis/ExRootFilter.h \
    … …  
    9311022        tmp/modules/AngularSmearing.$(ObjSuf) \
    9321023        tmp/modules/BTagging.$(ObjSuf) \
     1024        tmp/modules/BeamSpotFilter.$(ObjSuf) \
    9331025        tmp/modules/Calorimeter.$(ObjSuf) \
    9341026        tmp/modules/Cloner.$(ObjSuf) \
    … …  
    9551047        tmp/modules/PileUpJetID.$(ObjSuf) \
    9561048        tmp/modules/PileUpMerger.$(ObjSuf) \
     1049        tmp/modules/RecoPuFilter.$(ObjSuf) \
    9571050        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    9581051        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    9631056        tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
    9641057        tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
     1058        tmp/modules/TrackSmearing.$(ObjSuf) \
    9651059        tmp/modules/TreeWriter.$(ObjSuf) \
    9661060        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
     1061        tmp/modules/VertexFinder.$(ObjSuf) \
     1062        tmp/modules/VertexFinderDA4D.$(ObjSuf) \
     1063        tmp/modules/VertexSorter.$(ObjSuf) \
    9671064        tmp/modules/Weighter.$(ObjSuf)
    9681065
    … …  
    15691666        tmp/external/tcl/tclVar.$(ObjSuf)
    15701667
     1668modules/VertexFinderDA4D.h: \
     1669        classes/DelphesModule.h \
     1670        classes/DelphesClasses.h
     1671        @touch $@
     1672
     1673modules/TrackSmearing.h: \
     1674        classes/DelphesModule.h
     1675        @touch $@
     1676
    15711677external/fastjet/ClusterSequence.hh: \
    15721678        external/fastjet/PseudoJet.hh \
    … …  
    18291935        @touch $@
    18301936
     1937modules/VertexSorter.h: \
     1938        classes/DelphesModule.h \
     1939        classes/DelphesClasses.h
     1940        @touch $@
     1941
    18311942modules/Delphes.h: \
    18321943        classes/DelphesModule.h
     1944        @touch $@
     1945
     1946modules/VertexFinder.h: \
     1947        classes/DelphesModule.h \
     1948        classes/DelphesClasses.h
    18331949        @touch $@
    18341950
    … …  
    19022018        @touch $@
    19032019
     2020modules/RecoPuFilter.h: \
     2021        classes/DelphesModule.h
     2022        @touch $@
     2023
    19042024modules/Hector.h: \
    19052025        classes/DelphesModule.h
    … …  
    20012121
    20022122modules/FastJetFinder.h: \
     2123        classes/DelphesModule.h
     2124        @touch $@
     2125
     2126modules/BeamSpotFilter.h: \
    20032127        classes/DelphesModule.h
    20042128        @touch $@

  • cards/converter_card.tcl

    rec5e04b r21eab4f  
    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_ATLAS.tcl

    rec5e04b r21eab4f  
    141141
    142142  # resolution formula for charged hadrons
    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)}
     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)}
    146146}
    147147

  • cards/delphes_card_ATLAS_PileUp.tcl

    rec5e04b r21eab4f  
    177177  # resolution formula for charged hadrons
    178178  # based on arXiv:1405.6569
    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)}
     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)}
    182182}
    183183

  • cards/delphes_card_CMS.tcl

    rec5e04b r21eab4f  
    1515
    1616  TrackMerger
     17 
     18  ECal
     19  HCal
     20 
    1721  Calorimeter
    1822  EFlowMerger
    … …  
    124128  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
    125129                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
    126                                            (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
     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
    127133                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
    128                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
     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)) +
    129136                         (abs(eta) > 2.5)                                                  * (0.00)}
    130137}
    … …  
    142149  # resolution formula for charged hadrons
    143150  # based on arXiv:1405.6569
    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)}
     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)}
    147154}
    148155
    … …  
    192199}
    193200
     201
     202
    194203#############
    195 # Calorimeter
     204#   ECAL
    196205#############
    197206
    198 module Calorimeter Calorimeter {
     207module SimpleCalorimeter ECal {
    199208  set ParticleInputArray ParticlePropagator/stableParticles
    200209  set TrackInputArray TrackMerger/tracks
    201210
    202   set TowerOutputArray towers
    203   set PhotonOutputArray photons
    204 
     211  set TowerOutputArray ecalTowers
    205212  set EFlowTrackOutputArray eflowTracks
    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
     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
     308module 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
    214320
    215321  set SmearTowerCenter true
    … …  
    249355
    250356  # default energy fractions {abs(PDG code)} {Fecal Fhcal}
    251   add EnergyFraction {0} {0.0 1.0}
     357  add EnergyFraction {0} {1.0}
    252358  # energy fractions for e, gamma and pi0
    253   add EnergyFraction {11} {1.0 0.0}
    254   add EnergyFraction {22} {1.0 0.0}
    255   add EnergyFraction {111} {1.0 0.0}
     359  add EnergyFraction {11} {0.0}
     360  add EnergyFraction {22} {0.0}
     361  add EnergyFraction {111} {0.0}
    256362  # energy fractions for muon, neutrinos and neutralinos
    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}
     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}
    266372  # energy fractions for K0short and Lambda
    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)}
     373  add EnergyFraction {310} {0.7}
     374  add EnergyFraction {3122} {0.7}
    275375
    276376  # set HCalResolutionFormula {resolution formula as a function of eta and energy}
    277   set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
     377  set ResolutionFormula {                      (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
    278378                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
    279 }
     379
     380}
     381
     382
     383#################
     384# Electron filter
     385#################
     386
     387module 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
     399module Merger Calorimeter {
     400# add InputArray InputArray
     401  add InputArray ECal/ecalTowers
     402  add InputArray HCal/hcalTowers
     403  set OutputArray towers
     404}
     405
     406
    280407
    281408####################
    … …  
    285412module Merger EFlowMerger {
    286413# add InputArray InputArray
    287   add InputArray Calorimeter/eflowTracks
    288   add InputArray Calorimeter/eflowPhotons
    289   add InputArray Calorimeter/eflowNeutralHadrons
     414  add InputArray HCal/eflowTracks
     415  add InputArray ECal/eflowPhotons
     416  add InputArray HCal/eflowNeutralHadrons
    290417  set OutputArray eflow
    291418}
    … …  
    296423
    297424module Efficiency PhotonEfficiency {
    298   set InputArray Calorimeter/eflowPhotons
     425  set InputArray ECal/eflowPhotons
    299426  set OutputArray photons
    300427
    … …  
    325452}
    326453
    327 #################
    328 # Electron filter
    329 #################
    330 
    331 module PdgCodeFilter ElectronFilter {
    332   set InputArray Calorimeter/eflowTracks
    333   set OutputArray electrons
    334   set Invert true
    335   add PdgCode {11}
    336   add PdgCode {-11}
    337 }
    338454
    339455#####################
    … …  
    382498
    383499  # efficiency formula for muons
    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)) +
     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) +
    389503                         (abs(eta) > 2.4)                                                 * (0.00)}
    390504}
    … …  
    602716  add Branch Calorimeter/towers Tower Tower
    603717
    604   add Branch Calorimeter/eflowTracks EFlowTrack Track
    605   add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
    606   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
     718  add Branch HCal/eflowTracks EFlowTrack Track
     719  add Branch ECal/eflowPhotons EFlowPhoton Tower
     720  add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
    607721
    608722  add Branch GenJetFinder/jets GenJet Jet
    609723  add Branch GenMissingET/momentum GenMissingET MissingET
    610 
     724 
    611725  add Branch UniqueObjectFinder/jets Jet Jet
    612726  add Branch UniqueObjectFinder/electrons Electron Electron

  • cards/delphes_card_CMS_NoFastJet.tcl

    rec5e04b r21eab4f  
    128128  # resolution formula for electrons
    129129  # based on arXiv:1405.6569
    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)}
     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)}
    133133}
    134134
    … …  
    144144
    145145  # resolution formula for muons
    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 }
     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
    150151
    151152##############

  • cards/delphes_card_CMS_PileUp.tcl

    rec5e04b r21eab4f  
    178178  # resolution formula for charged hadrons
    179179  # based on arXiv:1405.6569
    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)}
     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)}
    183183}
    184184

  • classes/DelphesClasses.cc

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

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

  • modules/Calorimeter.cc

    rec5e04b r21eab4f  
    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
    … …  
    219213  Double_t ecalEnergy, hcalEnergy;
    220214  Double_t ecalSigma, hcalSigma;
     215  Double_t energyGuess;
    221216  Int_t pdgCode;
    222217
    … …  
    368363      fHCalTowerEnergy = 0.0;
    369364
    370       fECalTrackEnergy[0] = 0.0;
    371       fECalTrackEnergy[1] = 0.0;
    372 
    373       fHCalTrackEnergy[0] = 0.0;
    374       fHCalTrackEnergy[1] = 0.0;
    375 
     365      fECalTrackEnergy = 0.0;
     366      fHCalTrackEnergy = 0.0;
     367     
     368      fECalTrackSigma = 0.0;
     369      fHCalTrackSigma = 0.0;
     370     
    376371      fTowerTrackHits = 0;
    377372      fTowerPhotonHits = 0;
    378373
    379       fECalTowerTrackArray[0]->Clear();
    380       fECalTowerTrackArray[1]->Clear();
    381 
    382       fHCalTowerTrackArray[0]->Clear();
    383       fHCalTowerTrackArray[1]->Clear();
     374      fECalTowerTrackArray->Clear();
     375      fHCalTowerTrackArray->Clear();
     376   
    384377    }
    385378
    … …  
    406399      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    407400      {
    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         }
     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);
    419408      }
     409     
    420410      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    421411      {
     412        fHCalTrackEnergy += hcalEnergy;
    422413        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         }
     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);
    433419      }
     420     
    434421      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    435422      {
    … …  
    476463  Double_t energy, pt, eta, phi;
    477464  Double_t ecalEnergy, hcalEnergy;
     465  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
     466 
    478467  Double_t ecalSigma, hcalSigma;
    479 
     468  Double_t ecalNeutralSigma, hcalNeutralSigma;
     469
     470  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
     471 
    480472  TLorentzVector momentum;
    481473  TFractionMap::iterator itFractionMap;
    … …  
    484476
    485477  if(!fTower) return;
    486 
    487478
    488479  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
    … …  
    554545    fTowerOutputArray->Add(fTower);
    555546  }
    556 
     547 
    557548  // 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)
     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)
    618561  {
    619562    // create new photon tower
    620563    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;
     564    pt =  ecalNeutralEnergy / TMath::CosH(eta);
     565   
     566    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     567    tower->Eem = ecalNeutralEnergy;
    626568    tower->Ehad = 0.0;
    627569    tower->PID = 22;
    628 
     570   
    629571    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron 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
    634614    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     615    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     616   
     617    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     618    tower->Ehad = hcalNeutralEnergy;
    639619    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     620   
    642621    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     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 
    644660}
    645661

  • modules/Calorimeter.h

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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"
    6062#include "modules/RecoPuFilter.h"
    6163#include "modules/Cloner.h"
    … …  
    6466#include "modules/JetFlavorAssociation.h"
    6567#include "modules/JetFakeParticle.h"
     68#include "modules/VertexSorter.h"
     69#include "modules/VertexFinder.h"
     70#include "modules/VertexFinderDA4D.h"
    6671#include "modules/ExampleModule.h"
    6772
    … …  
    8186#pragma link C++ class EnergySmearing+;
    8287#pragma link C++ class MomentumSmearing+;
     88#pragma link C++ class TrackSmearing+;
    8389#pragma link C++ class ImpactParameterSmearing+;
    8490#pragma link C++ class TimeSmearing+;
    … …  
    104110#pragma link C++ class StatusPidFilter+;
    105111#pragma link C++ class PdgCodeFilter+;
     112#pragma link C++ class BeamSpotFilter+;
    106113#pragma link C++ class RecoPuFilter+;
    107114#pragma link C++ class Cloner+;
    … …  
    110117#pragma link C++ class JetFlavorAssociation+;
    111118#pragma link C++ class JetFakeParticle+;
     119#pragma link C++ class VertexSorter+;
     120#pragma link C++ class VertexFinder+;
     121#pragma link C++ class VertexFinderDA4D+;
    112122#pragma link C++ class ExampleModule+;
    113123

  • modules/MomentumSmearing.cc

    rec5e04b r21eab4f  
    9696{
    9797  Candidate *candidate, *mother;
    98   Double_t pt, eta, phi, e;
     98  Double_t pt, eta, phi, e, res;
    9999
    100100  fItInputArray->Reset();
    … …  
    107107    pt = candidateMomentum.Pt();
    108108    e = candidateMomentum.E();
     109    res = fFormula->Eval(pt, eta, phi, e);
     110 
     111    // apply smearing formula
     112    //pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     113   
     114    res = ( res > 1.0 ) ? 1.0 : res;
    109115
    110     // apply smearing formula
    111     pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     116    pt = LogNormal(pt, res * pt );
    112117   
    113     if(pt <= 0.0) continue;
     118    //if(pt <= 0.0) continue;
    114119
    115120    mother = candidate;
    … …  
    118123    phi = candidateMomentum.Phi();
    119124    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt*TMath::CosH(eta));
    120     candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     125    //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     126    candidate->TrackResolution = res;
    121127    candidate->AddCandidate(mother);
    122128       
    … …  
    124130  }
    125131}
     132//----------------------------------------------------------------
     133
     134Double_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
    126151
    127152//------------------------------------------------------------------------------

  • modules/MomentumSmearing.h

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

  • modules/ParticlePropagator.cc

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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/PileUpMerger.cc

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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
    … …  
    199191  Double_t energy;
    200192  Double_t sigma;
     193  Double_t energyGuess;
     194
    201195  Int_t pdgCode;
    202196
    … …  
    340334      fTowerEnergy = 0.0;
    341335
    342       fTrackEnergy[0] = 0.0;
    343       fTrackEnergy[1] = 0.0;
     336      fTrackEnergy = 0.0;
     337      fTrackSigma = 0.0;
    344338
    345339      fTowerTime = 0.0;
    … …  
    351345      fTowerPhotonHits = 0;
    352346
    353       fTowerTrackArray[0]->Clear();
    354       fTowerTrackArray[1]->Clear();
    355     }
     347      fTowerTrackArray->Clear();
     348     }
    356349
    357350    // check for track hits
    … …  
    371364      if(fTrackFractions[number] > 1.0E-9)
    372365      {
    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         }
     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);
    384375      }
     376       
    385377      else
    386378      {
    … …  
    403395    fTowerEnergy += energy;
    404396
    405     fTowerTime += TMath::Sqrt(energy)*position.T();
    406     fTowerTimeWeight += TMath::Sqrt(energy);
     397    fTowerTime += energy*position.T();
     398    fTowerTimeWeight += energy;
    407399
    408400    fTower->AddCandidate(particle);
    … …  
    418410{
    419411  Candidate *tower, *track, *mother;
    420   Double_t energy, pt, eta, phi;
    421   Double_t sigma;
     412  Double_t energy,neutralEnergy, pt, eta, phi;
     413  Double_t sigma, neutralSigma;
    422414  Double_t time;
     415   
     416  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    423417
    424418  TLorentzVector momentum;
    … …  
    436430
    437431  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
     432
    438433
    439434  if(fSmearTowerCenter)
    … …  
    464459  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    465460
    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)
     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)
    499474  {
    500475    // create new photon tower
    501476    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;
     477    pt = neutralEnergy / TMath::CosH(eta);
     478
     479    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
     480    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     481    tower->PID = (fIsEcal) ? 22 : 0;
     482 
     483    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     484    fEFlowTowerOutputArray->Add(tower);
    506485   
    507     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     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  }
    508496   
    509     tower->PID = (fIsEcal) ? 22 : 0;
    510    
    511     fEFlowTowerOutputArray->Add(tower);
    512   }
    513 }
     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 }
    514520
    515521//------------------------------------------------------------------------------

  • modules/SimpleCalorimeter.h

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

    rec5e04b r21eab4f  
    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

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

  • readers/DelphesPythia8.cpp

    rec5e04b r21eab4f  
    153153// from pythia8 example 21
    154154
    155 void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
    156   Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
    157 
     155void fillParticle(int id, double pt_max, double eta_max,
     156  Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
     157{
    158158  // Reset event record to allow for new event.
    159159  event.reset();
    160160
    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);
     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);
    169169
    170170  // Store the particle in the event record.
    171   event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
    172 
     171  event.append(id, 1, 0, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
    173172}
    174173
    175 void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
    176   Pythia8::Rndm& rndm) {
     174void fillPartons(int id, double pt_max, double eta_max,
     175  Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
     176{
    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 
    184181  // Information on a q qbar system, to be hadronized.
    185182
    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);
     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);
    196196  }
    197197  else
    198198  {
    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);
     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);
    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") == 22)
     340        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 15 || pythia->mode("Main:spareMode1") == 22)
    341341        {
    342           fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
     342          fillParticle(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
    343343        }
    344         else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
     344        else
     345        {
     346          fillPartons(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
     347        }
    345348      }
    346349
Note: See TracChangeset for help on using the changeset viewer.