Fork me on GitHub

Changes in / [a9c67b3a:3c46e17] in git


Ignore:
Files:
5 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    ra9c67b3a r3c46e17  
    10011001endif
    10021002
    1003 tmp/external/PUPPI/PuppiAlgo.$(ObjSuf): \
    1004         external/PUPPI/PuppiAlgo.$(SrcSuf) \
    1005         external/fastjet/internal/base.hh
    1006 tmp/external/PUPPI/PuppiContainer.$(ObjSuf): \
    1007         external/PUPPI/PuppiContainer.$(SrcSuf) \
    1008         external/fastjet/internal/base.hh \
     1003tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
     1004        external/PUPPI/puppiCleanContainer.$(SrcSuf) \
    10091005        external/fastjet/Selector.hh
    10101006tmp/external/fastjet/AreaDefinition.$(ObjSuf): \
    … …  
    13591355        modules/RunPUPPI.$(SrcSuf) \
    13601356        modules/RunPUPPI.h \
    1361         external/PUPPI/RecoObj2.hh \
    1362         external/PUPPI/AlgoObj.hh \
     1357        external/PUPPI/puppiCleanContainer.hh \
     1358        external/PUPPI/RecoObj.hh \
     1359        external/PUPPI/puppiParticle.hh \
     1360        external/PUPPI/puppiAlgoBin.hh \
    13631361        classes/DelphesClasses.h \
    13641362        classes/DelphesFactory.h \
    13651363        classes/DelphesFormula.h
    13661364FASTJET_OBJ +=  \
    1367         tmp/external/PUPPI/PuppiAlgo.$(ObjSuf) \
    1368         tmp/external/PUPPI/PuppiContainer.$(ObjSuf) \
     1365        tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
    13691366        tmp/external/fastjet/AreaDefinition.$(ObjSuf) \
    13701367        tmp/external/fastjet/BasicRandom.$(ObjSuf) \
    … …  
    17951792
    17961793modules/RunPUPPI.h: \
    1797         classes/DelphesModule.h \
    1798         external/PUPPI/PuppiContainer.hh
     1794        classes/DelphesModule.h
    17991795        @touch $@
    18001796
    … …  
    19281924        external/fastjet/AreaDefinition.hh \
    19291925        external/fastjet/ClusterSequenceAreaBase.hh
     1926        @touch $@
     1927
     1928external/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
    19301934        @touch $@
    19311935

  • cards/CMS_PhaseII/CMS_PhaseII_200PU.tcl

    ra9c67b3a r3c46e17  
    4040  EFlowMerger
    4141
    42   LeptonFilterNoLep
    43   LeptonFilterLep
    44   RunPUPPIBase
    4542  RunPUPPI
    46 
    4743
    4844  PhotonIsolation
    … …  
    5955
    6056  MissingET
    61   PuppiMissingET
    6257  GenMissingET
    6358  GenPileUpMissingET
    … …  
    9691
    9792  # pre-generated minbias input file
    98   set PileUpFile MinBias.pileup
     93  set PileUpFile ../eos/cms/store/group/upgrade/delphes/PhaseII/MinBias_100k.pileup
    9994
    10095  # average expected pile up
    … …  
    110105 
    111106  #set VertexDistributionFormula {exp(-(t^2/(2*(0.05/2.99792458E8*exp(-(z^2/(2*(0.05)^2))))^2)))}
    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)}
     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)}
    116113
    117114}
    … …  
    550547#########################################
    551548
    552 module 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 
    562 module 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 
    572 module RunPUPPI RunPUPPIBase {
     549module RunPUPPI RunPUPPI {
    573550  ## input information
    574   set TrackInputArray   LeptonFilterNoLep/eflowTracksNoLeptons
     551  set TrackInputArray   TrackMerger/tracks
    575552  set NeutralInputArray NeutralEFlowMerger/eflowTowers
    576553  set PVInputArray      PileUpMerger/vertices
    577554  set MinPuppiWeight    0.05
    578555  set UseExp            false
    579   set UseNoLep          false
    580556 
    581557  ## define puppi algorithm parameters (more than one for the same eta region is possible)                                                                                     
    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
     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
    595570
    596571  ## output name
    … …  
    600575}
    601576
    602 module Merger RunPUPPI {
    603   add InputArray RunPUPPIBase/PuppiParticles
    604   add InputArray LeptonFilterLep/eflowTracksLeptons
    605   set OutputArray PuppiParticles
    606 }
     577
    607578
    608579###################
    … …  
    614585#  add InputArray RunPUPPI/PuppiParticles
    615586  add InputArray EFlowMerger/eflow
    616   set MomentumOutputArray momentum
    617 }
    618 
    619 module Merger PuppiMissingET {
    620   #add InputArray InputArray
    621   add InputArray RunPUPPI/PuppiParticles
    622   #add InputArray EFlowMerger/eflow
    623587  set MomentumOutputArray momentum
    624588}
    … …  
    1007971 
    1008972  add Branch MissingET/momentum MissingET MissingET
    1009   add Branch PuppiMissingET/momentum PuppiMissingET MissingET
    1010973  add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
    1011974  add Branch ScalarHT/energy ScalarHT ScalarHT

  • cards/delphes_card_ATLAS.tcl

    ra9c67b3a r3c46e17  
    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
    … …  
    157157
    158158  # resolution formula for electrons
    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)}
     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)}
    162162}
    163163
    … …  
    171171
    172172  # set ResolutionFormula {resolution formula as a function of eta and pt}
     173
    173174  # resolution formula for muons
    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)}
     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)}
    177178}
    178179

  • cards/delphes_card_ATLAS_PileUp.tcl

    ra9c67b3a r3c46e17  
    182182
    183183  # resolution formula for charged hadrons
    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)}
     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)}
    188187}
    189188
    … …  
    199198
    200199  # resolution formula for electrons
    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)}
     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)}
    205203}
    206204
    … …  
    216214
    217215  # resolution formula for muons
    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)}
     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)}
    221219}
    222220

  • cards/delphes_card_CMS.tcl

    ra9c67b3a r3c46e17  
    142142  # resolution formula for charged hadrons
    143143  # 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)}
     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)}
    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.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)}
     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)}
    164164}
    165165
    … …  
    175175
    176176  # resolution formula for muons
    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)}
     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)}
    180180}
    181181

  • cards/delphes_card_CMS_PileUp.tcl

    ra9c67b3a r3c46e17  
    184184  # resolution formula for charged hadrons
    185185  # based on arXiv:1405.6569
    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)}
     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)}
    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.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)}
     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)}
    206206}
    207207
    … …  
    217217
    218218  # resolution formula for muons
    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)}
     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)}
    222222}
    223223

  • modules/RunPUPPI.cc

    ra9c67b3a r3c46e17  
    11#include "modules/RunPUPPI.h"
    22
    3 #include "PUPPI/RecoObj2.hh"
    4 #include "PUPPI/AlgoObj.hh"
    5 //#include "PUPPI/puppiParticle.hh"
    6 //#include "PUPPI/puppiAlgoBin.hh"
     3#include "PUPPI/puppiCleanContainer.hh"
     4#include "PUPPI/RecoObj.hh"
     5#include "PUPPI/puppiParticle.hh"
     6#include "PUPPI/puppiAlgoBin.hh"
    77
    88#include "classes/DelphesClasses.h"
    … …  
    3030
    3131void RunPUPPI::Init(){
     32
    3233  // input collection
    3334  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    … …  
    3738  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
    3839  fPVItInputArray      = fPVInputArray->MakeIterator();
    39   // puppi parameters                                 
    40   fApplyNoLep     = GetBool("UseNoLep", true);   
     40
     41
     42  // puppi parameters                                     
    4143  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    4244  fUseExp         = GetBool("UseExp", false);
    43   // read eta min ranges
     45
     46  // read eta min ranges                                                                                                                                                           
    4447  ExRootConfParam param = GetParam("EtaMinBin");
    4548  fEtaMinBin.clear();
    4649  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
    47   // read eta max ranges
     50
     51  // read eta max ranges                                                                                                                                                           
    4852  param = GetParam("EtaMaxBin");
    4953  fEtaMaxBin.clear();
    5054  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
    51   // read pt min value
     55
     56  // read pt min value                                                                                                                                                           
    5257  param = GetParam("PtMinBin");
    5358  fPtMinBin.clear();
    5459  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
     60
    5561  // read cone size                                                                                                                                                           
    5662  param = GetParam("ConeSizeBin");
    5763  fConeSizeBin.clear();
    5864  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
     65
    5966  // read RMS min pt                                                                                                                                             
    6067  param = GetParam("RMSPtMinBin");
    6168  fRMSPtMinBin.clear();
    6269  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
    63   // read RMS scale factor
     70
     71  // read RMS scale factor                                                                                                                                                           
    6472  param = GetParam("RMSScaleFactorBin");
    6573  fRMSScaleFactorBin.clear();
    6674  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
    67   // read neutral pt min cut
     75
     76  // read neutral pt min cut                                                                                                                                                           
    6877  param = GetParam("NeutralMinEBin");
    6978  fNeutralMinEBin.clear();
    7079  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
    71   // read neutral pt min slope
     80
     81  // read neutral pt min slope                                                                                                                                                           
    7282  param = GetParam("NeutralPtSlope");
    7383  fNeutralPtSlope.clear();
    7484  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
     85
    7586  // read apply chs                                                                                                                                                           
    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
     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                                                                                                                                                           
    8092  param = GetParam("UseCharged");
    8193  fUseCharged.clear();
    8294  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
    83   // read apply chs correction
     95
     96  // read apply chs correction                                                                                                                                                           
    8497  param = GetParam("ApplyLowPUCorr");
    8598  fApplyLowPUCorr.clear();
    8699  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
     100 
    87101  // read metric id                                                                                                                                                         
    88102  param = GetParam("MetricId");
    89103  fMetricId.clear();
    90104  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
    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());
     105
    95106  // create output array
    96107  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
    97108  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    98109  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);
    136110}
    137111
    … …  
    150124  TLorentzVector momentum;
    151125
    152   //DelphesFactory *factory = GetFactory();
     126  DelphesFactory *factory = GetFactory();
    153127
    154128  // loop over input objects
    155   fItTrackInputArray   ->Reset();
    156   fItNeutralInputArray ->Reset();
    157   fPVItInputArray      ->Reset();
     129  fItTrackInputArray->Reset();
     130  fItNeutralInputArray->Reset();
     131  fPVItInputArray->Reset();
    158132
    159133  std::vector<Candidate *> InputParticles;
    … …  
    168142  std::vector<RecoObj> puppiInputVector;
    169143  puppiInputVector.clear();
    170   int lNBad  = 0;
     144
    171145  // Loop on charge track candidate
    172146  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
     147
    173148      momentum = candidate->Momentum;
     149     
    174150      RecoObj curRecoObj;
    175151      curRecoObj.pt  = momentum.Pt();
    … …  
    178154      curRecoObj.m   = momentum.M(); 
    179155      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;
    182156      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
    183         lNBad++;
    184157        curRecoObj.id    = 2;
    185         curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
     158        curRecoObj.vtxId = candidate->IsPU;
    186159        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
    187160        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
    … …  
    210183  // Loop on neutral calo cells
    211184  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
     185
    212186      momentum = candidate->Momentum;
     187
    213188      RecoObj curRecoObj;
    214189      curRecoObj.pt  = momentum.Pt();
    … …  
    216191      curRecoObj.phi = momentum.Phi();
    217192      curRecoObj.m   = momentum.M();
    218       curRecoObj.charge = 0;
    219193      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
     194
    220195
    221196      if(candidate->Charge == 0){
    … …  
    235210      InputParticles.push_back(candidate);
    236211  }
     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
    237243  // Create PUPPI container
    238   fPuppi->initialize(puppiInputVector);
    239   fPuppi->puppiWeights();
    240   std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
     244  puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
     245  std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
    241246
    242247  // Loop on final particles

  • modules/RunPUPPI.h

    ra9c67b3a r3c46e17  
    33
    44#include "classes/DelphesModule.h"
    5 #include "PUPPI/PuppiContainer.hh"
    65#include <vector>
    76
    … …  
    3029  const TObjArray *fNeutralInputArray; //!
    3130  const TObjArray *fPVInputArray; //!                                                                                                                                                     
    32   PuppiContainer* fPuppi;
     31 
    3332  // puppi parameters
    34   bool fApplyNoLep;
    35   double fMinPuppiWeight;
     33  float fMinPuppiWeight;
    3634  bool fUseExp;
    3735 
    … …  
    4846  std::vector<bool>  fApplyLowPUCorr;
    4947  std::vector<int>   fMetricId;
    50   std::vector<int>   fCombId;
    5148
    5249  TObjArray *fOutputArray;

  • readers/DelphesCMSFWLite.cpp

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

  • readers/DelphesPythia8.cpp

    ra9c67b3a r3c46e17  
    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 
    154 void 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 
    174 void 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 
    203146
    204147//---------------------------------------------------------------------------
    … …  
    270213    // Initialize Pythia
    271214    pythia = new Pythia8::Pythia;
    272     //Pythia8::Event& event = pythia->event;
    273     //Pythia8::ParticleData& pdt = pythia->particleData;
    274215
    275216    if(pythia == NULL)
    … …  
    289230    timesAllowErrors = pythia->mode("Main:timesAllowErrors");
    290231
    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       }
     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");
    307244    }
    308245
    … …  
    323260      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
    324261        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       }
    334262
    335263      if(!pythia->next())
Note: See TracChangeset for help on using the changeset viewer.