Fork me on GitHub

Changes in / [3c46e17:a9c67b3a] in git


Ignore:
Files:
1 added
5 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r3c46e17 ra9c67b3a  
    10011001endif
    10021002
    1003 tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
    1004         external/PUPPI/puppiCleanContainer.$(SrcSuf) \
     1003tmp/external/PUPPI/PuppiAlgo.$(ObjSuf): \
     1004        external/PUPPI/PuppiAlgo.$(SrcSuf) \
     1005        external/fastjet/internal/base.hh
     1006tmp/external/PUPPI/PuppiContainer.$(ObjSuf): \
     1007        external/PUPPI/PuppiContainer.$(SrcSuf) \
     1008        external/fastjet/internal/base.hh \
    10051009        external/fastjet/Selector.hh
    10061010tmp/external/fastjet/AreaDefinition.$(ObjSuf): \
    … …  
    13551359        modules/RunPUPPI.$(SrcSuf) \
    13561360        modules/RunPUPPI.h \
    1357         external/PUPPI/puppiCleanContainer.hh \
    1358         external/PUPPI/RecoObj.hh \
    1359         external/PUPPI/puppiParticle.hh \
    1360         external/PUPPI/puppiAlgoBin.hh \
     1361        external/PUPPI/RecoObj2.hh \
     1362        external/PUPPI/AlgoObj.hh \
    13611363        classes/DelphesClasses.h \
    13621364        classes/DelphesFactory.h \
    13631365        classes/DelphesFormula.h
    13641366FASTJET_OBJ +=  \
    1365         tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
     1367        tmp/external/PUPPI/PuppiAlgo.$(ObjSuf) \
     1368        tmp/external/PUPPI/PuppiContainer.$(ObjSuf) \
    13661369        tmp/external/fastjet/AreaDefinition.$(ObjSuf) \
    13671370        tmp/external/fastjet/BasicRandom.$(ObjSuf) \
    … …  
    17921795
    17931796modules/RunPUPPI.h: \
    1794         classes/DelphesModule.h
     1797        classes/DelphesModule.h \
     1798        external/PUPPI/PuppiContainer.hh
    17951799        @touch $@
    17961800
    … …  
    19241928        external/fastjet/AreaDefinition.hh \
    19251929        external/fastjet/ClusterSequenceAreaBase.hh
    1926         @touch $@
    1927 
    1928 external/PUPPI/puppiCleanContainer.hh: \
    1929         external/PUPPI/RecoObj.hh \
    1930         external/PUPPI/puppiParticle.hh \
    1931         external/PUPPI/puppiAlgoBin.hh \
    1932         external/fastjet/internal/base.hh \
    1933         external/fastjet/PseudoJet.hh
    19341930        @touch $@
    19351931

  • cards/CMS_PhaseII/CMS_PhaseII_200PU.tcl

    r3c46e17 ra9c67b3a  
    4040  EFlowMerger
    4141
     42  LeptonFilterNoLep
     43  LeptonFilterLep
     44  RunPUPPIBase
    4245  RunPUPPI
     46
    4347
    4448  PhotonIsolation
    … …  
    5559
    5660  MissingET
     61  PuppiMissingET
    5762  GenMissingET
    5863  GenPileUpMissingET
    … …  
    9196
    9297  # pre-generated minbias input file
    93   set PileUpFile ../eos/cms/store/group/upgrade/delphes/PhaseII/MinBias_100k.pileup
     98  set PileUpFile MinBias.pileup
    9499
    95100  # average expected pile up
    … …  
    105110 
    106111  #set VertexDistributionFormula {exp(-(t^2/(2*(0.05/2.99792458E8*exp(-(z^2/(2*(0.05)^2))))^2)))}
    107  
    108  
    109   set VertexDistributionFormula { (abs(t) <= 1.0e-09) * (abs(z) <= 0.15) * (1.00) +
    110                                   (abs(t) >  1.0e-09) * (abs(z) <= 0.15) * (0.00) +
    111                                   (abs(t) <= 1.0e-09) * (abs(z) > 0.15)  * (0.00) +
    112                                   (abs(t) >  1.0e-09) * (abs(z) > 0.15)  * (0.00)}
     112  set VertexDistributionFormula { (abs(t) <= 1.6e-10) * (abs(z) <= 0.053) * (1.00) +
     113      (abs(t) >  1.6e-10) * (abs(z) <= 0.053) * (0.00) +
     114      (abs(t) <= 1.6e-10) * (abs(z) > 0.053)  * (0.00) +
     115      (abs(t) >  1.6e-10) * (abs(z) > 0.053)  * (0.00)}
    113116
    114117}
    … …  
    547550#########################################
    548551
    549 module RunPUPPI RunPUPPI {
     552module PdgCodeFilter LeptonFilterNoLep {
     553  set InputArray HCal/eflowTracks
     554  set OutputArray eflowTracksNoLeptons
     555  set Invert false
     556  add PdgCode {13}
     557  add PdgCode {-13}
     558  add PdgCode {11}
     559  add PdgCode {-11}
     560}
     561
     562module PdgCodeFilter LeptonFilterLep {
     563  set InputArray HCal/eflowTracks
     564  set OutputArray eflowTracksLeptons
     565  set Invert true
     566  add PdgCode {11}
     567  add PdgCode {-11}
     568  add PdgCode {13}
     569  add PdgCode {-13}
     570}
     571
     572module RunPUPPI RunPUPPIBase {
    550573  ## input information
    551   set TrackInputArray   TrackMerger/tracks
     574  set TrackInputArray   LeptonFilterNoLep/eflowTracksNoLeptons
    552575  set NeutralInputArray NeutralEFlowMerger/eflowTowers
    553576  set PVInputArray      PileUpMerger/vertices
    554577  set MinPuppiWeight    0.05
    555578  set UseExp            false
     579  set UseNoLep          false
    556580 
    557581  ## define puppi algorithm parameters (more than one for the same eta region is possible)                                                                                     
    558   add EtaMinBin           0.    2.5    2.5    3.0   3.0
    559   add EtaMaxBin           2.5   3.0    3.0    10.0  10.0
    560   add PtMinBin            0.    0.5    0.5    0.5   0.5
    561   add ConeSizeBin         0.25  0.25   0.25   0.25  0.25
    562   add RMSPtMinBin         0.1   0.5    0.5    0.5   0.5
    563   add RMSScaleFactorBin   1.0   1.0    1.0    1.0   1.0
    564   add NeutralMinEBin      0.2   1.0    1.0    1.5   1.5
    565   add NeutralPtSlope      0.02  0.02   0.02   0.02  0.02
    566   add ApplyCHS            true  true   true   true  true
    567   add UseCharged          true  false  false  false false
    568   add ApplyLowPUCorr      true  true   true   true  true
    569   add MetricId            5     5      1      5     1
     582  add EtaMinBin           0.0   1.5   4.0
     583  add EtaMaxBin           1.5   4.0   10.0
     584  add PtMinBin            0.0   0.0   0.0
     585  add ConeSizeBin         0.2   0.2   0.2
     586  add RMSPtMinBin         0.1   0.5   0.5
     587  add RMSScaleFactorBin   1.0   1.0   1.0
     588  add NeutralMinEBin      0.2   0.2   0.5
     589  add NeutralPtSlope      0.006 0.013 0.067
     590  add ApplyCHS            true  true  true
     591  add UseCharged          true  true  false
     592  add ApplyLowPUCorr      true  true  true
     593  add MetricId            5     5     5
     594  add CombId              0     0     0
    570595
    571596  ## output name
    … …  
    575600}
    576601
    577 
     602module Merger RunPUPPI {
     603  add InputArray RunPUPPIBase/PuppiParticles
     604  add InputArray LeptonFilterLep/eflowTracksLeptons
     605  set OutputArray PuppiParticles
     606}
    578607
    579608###################
    … …  
    585614#  add InputArray RunPUPPI/PuppiParticles
    586615  add InputArray EFlowMerger/eflow
     616  set MomentumOutputArray momentum
     617}
     618
     619module Merger PuppiMissingET {
     620  #add InputArray InputArray
     621  add InputArray RunPUPPI/PuppiParticles
     622  #add InputArray EFlowMerger/eflow
    587623  set MomentumOutputArray momentum
    588624}
    … …  
    9711007 
    9721008  add Branch MissingET/momentum MissingET MissingET
     1009  add Branch PuppiMissingET/momentum PuppiMissingET MissingET
    9731010  add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
    9741011  add Branch ScalarHT/energy ScalarHT ScalarHT

  • cards/delphes_card_ATLAS.tcl

    r3c46e17 ra9c67b3a  
    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
    … …  
    157157
    158158  # resolution formula for electrons
    159   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    160                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    161                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     159  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
     160                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
     161                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
    162162}
    163163
    … …  
    171171
    172172  # set ResolutionFormula {resolution formula as a function of eta and pt}
    173 
    174173  # resolution formula for muons
    175   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*2.0e-4^2) +
    176                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*3.0e-4^2) +
    177                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*6.0e-4^2)}
     174  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
     175                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
     176                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
    178177}
    179178

  • cards/delphes_card_ATLAS_PileUp.tcl

    r3c46e17 ra9c67b3a  
    182182
    183183  # resolution formula for charged hadrons
    184   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    185                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    186                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     184  # based on arXiv:1405.6569
     185  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
     186                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
     187                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
    187188}
    188189
    … …  
    198199
    199200  # resolution formula for electrons
    200   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    201                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    202                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     201  # based on arXiv:1405.6569
     202  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
     203                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
     204                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
    203205}
    204206
    … …  
    214216
    215217  # resolution formula for muons
    216   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*2.0e-4^2) +
    217                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*3.0e-4^2) +
    218                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
     218  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
     219                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
     220                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
    219221}
    220222

  • cards/delphes_card_CMS.tcl

    r3c46e17 ra9c67b3a  
    142142  # resolution formula for charged hadrons
    143143  # based on arXiv:1405.6569
    144   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    145                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    146                          (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)}
    147147}
    148148
    … …  
    159159  # resolution formula for electrons
    160160  # based on arXiv:1405.6569
    161   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    162                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    163                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     161  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
     162                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
     163                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
    164164}
    165165
    … …  
    175175
    176176  # resolution formula for muons
    177   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
    178                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
    179                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
     177  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
     178                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
     179                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
    180180}
    181181

  • cards/delphes_card_CMS_PileUp.tcl

    r3c46e17 ra9c67b3a  
    184184  # resolution formula for charged hadrons
    185185  # based on arXiv:1405.6569
    186   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    187                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    188                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     186  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
     187                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
     188                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
    189189}
    190190
    … …  
    201201  # resolution formula for electrons
    202202  # based on arXiv:1405.6569
    203   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
    204                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
    205                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
     203  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
     204                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
     205                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
    206206}
    207207
    … …  
    217217
    218218  # resolution formula for muons
    219   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
    220                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
    221                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
     219  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
     220                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
     221                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
    222222}
    223223

  • modules/RunPUPPI.cc

    r3c46e17 ra9c67b3a  
    11#include "modules/RunPUPPI.h"
    22
    3 #include "PUPPI/puppiCleanContainer.hh"
    4 #include "PUPPI/RecoObj.hh"
    5 #include "PUPPI/puppiParticle.hh"
    6 #include "PUPPI/puppiAlgoBin.hh"
     3#include "PUPPI/RecoObj2.hh"
     4#include "PUPPI/AlgoObj.hh"
     5//#include "PUPPI/puppiParticle.hh"
     6//#include "PUPPI/puppiAlgoBin.hh"
    77
    88#include "classes/DelphesClasses.h"
    … …  
    3030
    3131void RunPUPPI::Init(){
    32 
    3332  // input collection
    3433  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    … …  
    3837  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
    3938  fPVItInputArray      = fPVInputArray->MakeIterator();
    40 
    41 
    42   // puppi parameters                                     
     39  // puppi parameters                                 
     40  fApplyNoLep     = GetBool("UseNoLep", true);   
    4341  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    4442  fUseExp         = GetBool("UseExp", false);
    45 
    46   // read eta min ranges                                                                                                                                                           
     43  // read eta min ranges
    4744  ExRootConfParam param = GetParam("EtaMinBin");
    4845  fEtaMinBin.clear();
    4946  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
    50 
    51   // read eta max ranges                                                                                                                                                           
     47  // read eta max ranges
    5248  param = GetParam("EtaMaxBin");
    5349  fEtaMaxBin.clear();
    5450  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
    55 
    56   // read pt min value                                                                                                                                                           
     51  // read pt min value
    5752  param = GetParam("PtMinBin");
    5853  fPtMinBin.clear();
    5954  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
    60 
    6155  // read cone size                                                                                                                                                           
    6256  param = GetParam("ConeSizeBin");
    6357  fConeSizeBin.clear();
    6458  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
    65 
    6659  // read RMS min pt                                                                                                                                             
    6760  param = GetParam("RMSPtMinBin");
    6861  fRMSPtMinBin.clear();
    6962  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
    70 
    71   // read RMS scale factor                                                                                                                                                           
     63  // read RMS scale factor
    7264  param = GetParam("RMSScaleFactorBin");
    7365  fRMSScaleFactorBin.clear();
    7466  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
    75 
    76   // read neutral pt min cut                                                                                                                                                           
     67  // read neutral pt min cut
    7768  param = GetParam("NeutralMinEBin");
    7869  fNeutralMinEBin.clear();
    7970  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
    80 
    81   // read neutral pt min slope                                                                                                                                                           
     71  // read neutral pt min slope
    8272  param = GetParam("NeutralPtSlope");
    8373  fNeutralPtSlope.clear();
    8474  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
    85 
    8675  // read apply chs                                                                                                                                                           
    87   param = GetParam("ApplyCHS");
    88   fApplyCHS.clear();
    89   for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
    90 
    91   // read use charged                                                                                                                                                           
     76  //param = GetParam("ApplyCHS");
     77  //fApplyCHS.clear();
     78  //for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
     79  // read use charged
    9280  param = GetParam("UseCharged");
    9381  fUseCharged.clear();
    9482  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
    95 
    96   // read apply chs correction                                                                                                                                                           
     83  // read apply chs correction
    9784  param = GetParam("ApplyLowPUCorr");
    9885  fApplyLowPUCorr.clear();
    9986  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
    100  
    10187  // read metric id                                                                                                                                                         
    10288  param = GetParam("MetricId");
    10389  fMetricId.clear();
    10490  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
    105 
     91  // scheme for combining
     92  param = GetParam("CombId");
     93  fCombId.clear();
     94  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
    10695  // create output array
    10796  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
    10897  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    10998  fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
     99  // Create algorithm list for puppi
     100  std::vector<AlgoObj> puppiAlgo;
     101  if(puppiAlgo.empty()){
     102    if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
     103         and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size()
     104         and fEtaMinBin.size() == fUseCharged.size()
     105         and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
     106      std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
     107      std::exit(EXIT_FAILURE);
     108    }
     109  }
     110  for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
     111    AlgoObj algoTmp ;
     112    algoTmp.etaMin            = fEtaMinBin.at(iAlgo);
     113    algoTmp.etaMax            = fEtaMaxBin.at(iAlgo);
     114    algoTmp.ptMin             = fPtMinBin.at(iAlgo);
     115    algoTmp.minNeutralPt      = fNeutralMinEBin.at(iAlgo);
     116    algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
     117    //Eta Extrapolation stuff is missing
     118    //Loop through file requiring algos for same bins to be adjacent
     119    while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) {
     120      AlgoSubObj algoSubTmp;
     121      algoSubTmp.metricId          = fMetricId.at(iAlgo);
     122      algoSubTmp.useCharged        = fUseCharged.at(iAlgo);
     123      algoSubTmp.applyLowPUCorr    = fApplyLowPUCorr.at(iAlgo);
     124      algoSubTmp.combId            = fCombId.at(iAlgo);
     125      algoSubTmp.coneSize          = fConeSizeBin.at(iAlgo);
     126      algoSubTmp.rmsPtMin          = fRMSPtMinBin.at(iAlgo);
     127      algoSubTmp.rmsScaleFactor    = fRMSScaleFactorBin.at(iAlgo);
     128      algoTmp.subAlgos.push_back(algoSubTmp);
     129      iAlgo++;
     130    }
     131    iAlgo--;
     132    //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;   
     133    puppiAlgo.push_back(algoTmp);     
     134  }
     135  fPuppi  = new PuppiContainer(true,fUseExp,fMinPuppiWeight,puppiAlgo);
    110136}
    111137
    … …  
    124150  TLorentzVector momentum;
    125151
    126   DelphesFactory *factory = GetFactory();
     152  //DelphesFactory *factory = GetFactory();
    127153
    128154  // loop over input objects
    129   fItTrackInputArray->Reset();
    130   fItNeutralInputArray->Reset();
    131   fPVItInputArray->Reset();
     155  fItTrackInputArray   ->Reset();
     156  fItNeutralInputArray ->Reset();
     157  fPVItInputArray      ->Reset();
    132158
    133159  std::vector<Candidate *> InputParticles;
    … …  
    142168  std::vector<RecoObj> puppiInputVector;
    143169  puppiInputVector.clear();
    144 
     170  int lNBad  = 0;
    145171  // Loop on charge track candidate
    146172  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
    147 
    148173      momentum = candidate->Momentum;
    149      
    150174      RecoObj curRecoObj;
    151175      curRecoObj.pt  = momentum.Pt();
    … …  
    154178      curRecoObj.m   = momentum.M(); 
    155179      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
     180      //if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
     181      //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
    156182      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
     183        lNBad++;
    157184        curRecoObj.id    = 2;
    158         curRecoObj.vtxId = candidate->IsPU;
     185        curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
    159186        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
    160187        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
    … …  
    183210  // Loop on neutral calo cells
    184211  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
    185 
    186212      momentum = candidate->Momentum;
    187 
    188213      RecoObj curRecoObj;
    189214      curRecoObj.pt  = momentum.Pt();
    … …  
    191216      curRecoObj.phi = momentum.Phi();
    192217      curRecoObj.m   = momentum.M();
     218      curRecoObj.charge = 0;
    193219      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
    194 
    195220
    196221      if(candidate->Charge == 0){
    … …  
    210235      InputParticles.push_back(candidate);
    211236  }
    212 
    213   // Create algorithm list for puppi
    214   std::vector<puppiAlgoBin> puppiAlgo;
    215   if(puppiAlgo.empty()){
    216    if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
    217        and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size()
    218        and fEtaMinBin.size() == fApplyCHS.size()  and fEtaMinBin.size() == fUseCharged.size()
    219        and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
    220     std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
    221     std::exit(EXIT_FAILURE);
    222    }
    223 
    224    for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
    225     puppiAlgoBin algoTmp ;
    226     algoTmp.fEtaMin_ = fEtaMinBin.at(iAlgo);
    227     algoTmp.fEtaMax_ = fEtaMaxBin.at(iAlgo);
    228     algoTmp.fPtMin_  = fPtMinBin.at(iAlgo);
    229     algoTmp.fConeSize_        = fConeSizeBin.at(iAlgo);
    230     algoTmp.fRMSPtMin_        = fRMSPtMinBin.at(iAlgo);
    231     algoTmp.fRMSScaleFactor_  = fRMSScaleFactorBin.at(iAlgo);
    232     algoTmp.fNeutralMinE_     = fNeutralMinEBin.at(iAlgo);
    233     algoTmp.fNeutralPtSlope_  = fNeutralPtSlope.at(iAlgo);
    234     algoTmp.fApplyCHS_        = fApplyCHS.at(iAlgo);
    235     algoTmp.fUseCharged_      = fUseCharged.at(iAlgo);
    236     algoTmp.fApplyLowPUCorr_  = fApplyLowPUCorr.at(iAlgo);
    237     algoTmp.fMetricId_        = fMetricId.at(iAlgo);
    238     if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;   
    239     puppiAlgo.push_back(algoTmp);     
    240    }
    241   } 
    242 
    243237  // Create PUPPI container
    244   puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
    245   std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
     238  fPuppi->initialize(puppiInputVector);
     239  fPuppi->puppiWeights();
     240  std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
    246241
    247242  // Loop on final particles

  • modules/RunPUPPI.h

    r3c46e17 ra9c67b3a  
    33
    44#include "classes/DelphesModule.h"
     5#include "PUPPI/PuppiContainer.hh"
    56#include <vector>
    67
    … …  
    2930  const TObjArray *fNeutralInputArray; //!
    3031  const TObjArray *fPVInputArray; //!                                                                                                                                                     
    31  
     32  PuppiContainer* fPuppi;
    3233  // puppi parameters
    33   float fMinPuppiWeight;
     34  bool fApplyNoLep;
     35  double fMinPuppiWeight;
    3436  bool fUseExp;
    3537 
    … …  
    4648  std::vector<bool>  fApplyLowPUCorr;
    4749  std::vector<int>   fMetricId;
     50  std::vector<int>   fCombId;
    4851
    4952  TObjArray *fOutputArray;

  • readers/DelphesCMSFWLite.cpp

    r3c46e17 ra9c67b3a  
    164164    candidate->Momentum.SetPxPyPzE(px, py, pz, e);
    165165
    166     candidate->Position.SetXYZT(x, y, z, 0.0);
     166    candidate->Position.SetXYZT(x*10.0, y*10.0, z*10.0, 0.0);
    167167
    168168    allParticleOutputArray->Add(candidate);

  • readers/DelphesPythia8.cpp

    r3c46e17 ra9c67b3a  
    144144  interrupted = true;
    145145}
     146
     147
     148// Single-particle gun. The particle must be a colour singlet.
     149// Input: flavour, energy, direction (theta, phi).
     150// If theta < 0 then random choice over solid angle.
     151// Optional final argument to put particle at rest => E = m.
     152// from pythia8 example 21
     153
     154void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
     155  Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
     156
     157  // Reset event record to allow for new event.
     158  event.reset();
     159
     160  // Angles uniform in solid angle.
     161  double cThe, sThe, phi, ee;
     162  cThe = 2. * rndm.flat() - 1.;
     163  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     164  phi = 2. * M_PI * rndm.flat();
     165  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     166  double mm = pdt.mSel(id);
     167  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     168
     169  // Store the particle in the event record.
     170  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     171
     172}
     173
     174void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
     175  Pythia8::Rndm& rndm) {
     176
     177  // Reset event record to allow for new event.
     178  event.reset();
     179
     180  // Angles uniform in solid angle.
     181  double cThe, sThe, phi, ee;
     182
     183  // Information on a q qbar system, to be hadronized.
     184
     185  cThe = 2. * rndm.flat() - 1.;
     186  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     187  phi = 2. * M_PI * rndm.flat();
     188  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     189  double mm = pdt.m0(type);
     190  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     191  if (type == 21)
     192  {
     193    event.append( 21, 23, 101, 102, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee);
     194    event.append( 21, 23, 102, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee);
     195  }
     196  else
     197  {
     198    event.append(  type, 23, 101,   0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     199    event.append( -type, 23,   0, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee, mm);
     200  }
     201}
     202
    146203
    147204//---------------------------------------------------------------------------
    … …  
    213270    // Initialize Pythia
    214271    pythia = new Pythia8::Pythia;
     272    //Pythia8::Event& event = pythia->event;
     273    //Pythia8::ParticleData& pdt = pythia->particleData;
    215274
    216275    if(pythia == NULL)
    … …  
    230289    timesAllowErrors = pythia->mode("Main:timesAllowErrors");
    231290
    232     inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
    233     if(inputFile)
    234     {
    235       reader = new DelphesLHEFReader;
    236       reader->SetInputFile(inputFile);
    237 
    238       branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
    239       branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
    240 
    241       allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
    242       stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
    243       partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
     291    // Check if particle gun
     292    if (!pythia->flag("Main:spareFlag1"))
     293    {
     294      inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
     295      if(inputFile)
     296      {
     297        reader = new DelphesLHEFReader;
     298        reader->SetInputFile(inputFile);
     299
     300        branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
     301        branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
     302
     303        allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
     304        stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
     305        partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
     306      }
    244307    }
    245308
    … …  
    260323      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
    261324        stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
     325
     326      if (pythia->flag("Main:spareFlag1"))
     327      {
     328        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 22)
     329        {
     330          fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
     331        }
     332        else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
     333      }
    262334
    263335      if(!pythia->next())
Note: See TracChangeset for help on using the changeset viewer.