Fork me on GitHub

Changes in / [9343566:e33e6db] in git


Ignore:
Files:
27 deleted
47 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r9343566 re33e6db  
    257257        classes/DelphesClasses.h \
    258258        classes/DelphesFactory.h \
    259         classes/DelphesLHEFReader.h \
    260259        external/ExRootAnalysis/ExRootTreeWriter.h \
    261260        external/ExRootAnalysis/ExRootTreeBranch.h \
    … …  
    338337        modules/Weighter.h \
    339338        modules/Hector.h \
    340         modules/RunPUPPI.h \
    341         modules/JetFlavorAssociation.h \
    342339        modules/ExampleModule.h
    343340ModulesDict$(PcmSuf): \
    … …  
    524521tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
    525522        external/Hector/H_VerticalQuadrupole.$(SrcSuf)
    526 tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
    527         external/PUPPI/puppiCleanContainer.$(SrcSuf) \
    528         external/fastjet/Selector.hh
    529523tmp/modules/AngularSmearing.$(ObjSuf): \
    530524        modules/AngularSmearing.$(SrcSuf) \
    … …  
    658652        external/ExRootAnalysis/ExRootFilter.h \
    659653        external/ExRootAnalysis/ExRootClassifier.h
    660 tmp/modules/JetFlavorAssociation.$(ObjSuf): \
    661         modules/JetFlavorAssociation.$(SrcSuf) \
    662         modules/JetFlavorAssociation.h \
    663         classes/DelphesClasses.h \
    664         classes/DelphesFactory.h \
    665         classes/DelphesFormula.h \
    666         external/ExRootAnalysis/ExRootResult.h \
    667         external/ExRootAnalysis/ExRootFilter.h \
    668         external/ExRootAnalysis/ExRootClassifier.h
    669654tmp/modules/JetPileUpSubtractor.$(ObjSuf): \
    670655        modules/JetPileUpSubtractor.$(SrcSuf) \
    … …  
    745730        classes/DelphesClasses.h \
    746731        classes/DelphesFactory.h \
    747         classes/DelphesTF2.h \
     732        classes/DelphesFormula.h \
    748733        classes/DelphesPileUpReader.h \
    749734        external/ExRootAnalysis/ExRootResult.h \
    750735        external/ExRootAnalysis/ExRootFilter.h \
    751736        external/ExRootAnalysis/ExRootClassifier.h
    752 tmp/modules/RunPUPPI.$(ObjSuf): \
    753         modules/RunPUPPI.$(SrcSuf) \
    754         modules/RunPUPPI.h \
    755         external/PUPPI/puppiCleanContainer.hh \
    756         external/PUPPI/RecoObj.hh \
    757         external/PUPPI/puppiParticle.hh \
    758         external/PUPPI/puppiAlgoBin.hh \
    759         classes/DelphesClasses.h \
    760         classes/DelphesFactory.h \
    761         classes/DelphesFormula.h
    762737tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    763738        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    894869        tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
    895870        tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
    896         tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
    897871        tmp/modules/AngularSmearing.$(ObjSuf) \
    898872        tmp/modules/BTagging.$(ObjSuf) \
    … …  
    909883        tmp/modules/ImpactParameterSmearing.$(ObjSuf) \
    910884        tmp/modules/Isolation.$(ObjSuf) \
    911         tmp/modules/JetFlavorAssociation.$(ObjSuf) \
    912885        tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
    913886        tmp/modules/LeptonDressing.$(ObjSuf) \
    … …  
    918891        tmp/modules/PileUpJetID.$(ObjSuf) \
    919892        tmp/modules/PileUpMerger.$(ObjSuf) \
    920         tmp/modules/RunPUPPI.$(ObjSuf) \
    921893        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    922894        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    11081080tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf): \
    11091081        external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(SrcSuf)
    1110 tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf): \
    1111         external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(SrcSuf) \
    1112         external/fastjet/JetDefinition.hh \
    1113         external/fastjet/ClusterSequenceAreaBase.hh
    1114 tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf): \
    1115         external/fastjet/contribs/RecursiveTools/Recluster.$(SrcSuf)
    1116 tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf): \
    1117         external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(SrcSuf) \
    1118         external/fastjet/JetDefinition.hh \
    1119         external/fastjet/ClusterSequenceAreaBase.hh
    1120 tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf): \
    1121         external/fastjet/contribs/RecursiveTools/SoftDrop.$(SrcSuf)
    11221082tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf): \
    11231083        external/fastjet/contribs/SoftKiller/SoftKiller.$(SrcSuf)
    … …  
    12531213        external/fastjet/contribs/Nsubjettiness/Njettiness.hh \
    12541214        external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh \
    1255         external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh \
    1256         external/fastjet/tools/Filter.hh \
    1257         external/fastjet/tools/Pruner.hh \
    1258         external/fastjet/contribs/RecursiveTools/SoftDrop.hh
     1215        external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh
    12591216tmp/modules/FastJetGridMedianEstimator.$(ObjSuf): \
    12601217        modules/FastJetGridMedianEstimator.$(SrcSuf) \
    … …  
    13281285        tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf) \
    13291286        tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf) \
    1330         tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf) \
    1331         tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf) \
    1332         tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf) \
    1333         tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf) \
    13341287        tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf) \
    13351288        tmp/external/fastjet/plugins/ATLASCone/ATLASConePlugin.$(ObjSuf) \
    … …  
    13851338        display/Delphes3DGeometry.$(SrcSuf) \
    13861339        display/Delphes3DGeometry.h \
    1387         classes/DelphesClasses.h \
    1388         external/ExRootAnalysis/ExRootConfReader.h
     1340        external/ExRootAnalysis/ExRootConfReader.h \
     1341        classes/DelphesClasses.h
    13891342tmp/display/DelphesBranchElement.$(ObjSuf): \
    13901343        display/DelphesBranchElement.$(SrcSuf) \
    … …  
    13991352tmp/display/DelphesEventDisplay.$(ObjSuf): \
    14001353        display/DelphesEventDisplay.$(SrcSuf) \
     1354        external/ExRootAnalysis/ExRootConfReader.h \
     1355        external/ExRootAnalysis/ExRootTreeReader.h \
    14011356        display/DelphesCaloData.h \
    14021357        display/DelphesBranchElement.h \
    14031358        display/Delphes3DGeometry.h \
    14041359        display/DelphesEventDisplay.h \
    1405         display/DelphesDisplay.h \
    1406         display/Delphes3DGeometry.h \
    1407         display/DelphesHtmlSummary.h \
    1408         display/DelphesPlotSummary.h \
    1409         classes/DelphesClasses.h \
    1410         external/ExRootAnalysis/ExRootConfReader.h \
    1411         external/ExRootAnalysis/ExRootTreeReader.h
     1360        classes/DelphesClasses.h
    14121361tmp/display/DelphesHtmlSummary.$(ObjSuf): \
    14131362        display/DelphesHtmlSummary.$(SrcSuf) \
    … …  
    15921541        @touch $@
    15931542
     1543external/fastjet/Selector.hh: \
     1544        external/fastjet/PseudoJet.hh \
     1545        external/fastjet/RangeDefinition.hh
     1546        @touch $@
     1547
    15941548external/fastjet/internal/Dnn2piCylinder.hh: \
    15951549        external/fastjet/internal/DynamicNearestNeighbours.hh \
    15961550        external/fastjet/internal/DnnPlane.hh \
    15971551        external/fastjet/internal/numconsts.hh
    1598         @touch $@
    1599 
    1600 external/fastjet/Selector.hh: \
    1601         external/fastjet/PseudoJet.hh \
    1602         external/fastjet/RangeDefinition.hh
    16031552        @touch $@
    16041553
    … …  
    16901639        @touch $@
    16911640
    1692 modules/RunPUPPI.h: \
    1693         classes/DelphesModule.h
    1694         @touch $@
    1695 
    16961641modules/Cloner.h: \
    16971642        classes/DelphesModule.h
    … …  
    17301675        @touch $@
    17311676
     1677display/DelphesEventDisplay.h: \
     1678        external/ExRootAnalysis/ExRootTreeReader.h \
     1679        display/DelphesDisplay.h \
     1680        display/Delphes3DGeometry.h \
     1681        display/DelphesHtmlSummary.h \
     1682        display/DelphesPlotSummary.h
     1683        @touch $@
     1684
    17321685modules/TauTagging.h: \
    17331686        classes/DelphesModule.h \
    … …  
    17721725        @touch $@
    17731726
    1774 modules/JetFlavorAssociation.h: \
    1775         classes/DelphesModule.h \
    1776         classes/DelphesClasses.h
    1777         @touch $@
    1778 
    17791727modules/ParticlePropagator.h: \
    17801728        classes/DelphesModule.h
    … …  
    18091757        external/fastjet/AreaDefinition.hh \
    18101758        external/fastjet/ClusterSequenceAreaBase.hh
    1811         @touch $@
    1812 
    1813 external/PUPPI/puppiCleanContainer.hh: \
    1814         external/PUPPI/RecoObj.hh \
    1815         external/PUPPI/puppiParticle.hh \
    1816         external/PUPPI/puppiAlgoBin.hh \
    1817         external/fastjet/internal/base.hh \
    1818         external/fastjet/PseudoJet.hh
    18191759        @touch $@
    18201760

  • cards/delphes_card_ATLAS.tcl

    r9343566 re33e6db  
    295295
    296296module Efficiency PhotonEfficiency {
    297   set InputArray Calorimeter/eflowPhotons
     297  set InputArray Calorimeter/photons
    298298  set OutputArray photons
    299299

  • cards/delphes_card_ATLAS_PileUp.tcl

    r9343566 re33e6db  
    478478
    479479module Efficiency PhotonEfficiency {
    480   set InputArray Calorimeter/eflowPhotons
     480  set InputArray Calorimeter/photons
    481481  set OutputArray photons
    482482

  • cards/delphes_card_FCC_basic.tcl

    r9343566 re33e6db  
    226226  set EFlowTowerOutputArray eflowPhotons
    227227
    228   set IsEcal true
    229  
    230228  set EnergyMin 0.5
    231229  set EnergySignificanceMin 1.0
    … …  
    292290  set EFlowTowerOutputArray eflowNeutralHadrons
    293291
    294   set IsEcal false
    295  
    296292  set EnergyMin 1.0
    297293  set EnergySignificanceMin 1.0
    … …  
    435431
    436432module FastJetFinder FastJetFinder {
    437 #  set InputArray TowerMerger/towers
     433#  set InputArray Calorimeter/towers
    438434  set InputArray EFlowMerger/eflow
    439435

  • cards/delphes_card_LHCb.tcl

    r9343566 re33e6db  
    223223  set TowerOutputArray ecalTowers
    224224  set EFlowTowerOutputArray eflowPhotons
    225 
    226   set IsEcal true
    227225
    228226  set EnergyMin 0.0
    … …  
    305303  set TowerOutputArray hcalTowers
    306304  set EFlowTowerOutputArray eflowNeutralHadrons
    307 
    308   set IsEcal false
    309305
    310306  set EnergyMin 0.0

  • classes/ClassesLinkDef.h

    r9343566 re33e6db  
    4646#pragma link C++ class LHCOEvent+;
    4747#pragma link C++ class LHEFEvent+;
    48 #pragma link C++ class LHEFWeight+;
    4948#pragma link C++ class HepMCEvent+;
    5049#pragma link C++ class GenParticle+;

  • classes/DelphesClasses.cc

    r9343566 re33e6db  
    120120  PID(0), Status(0), M1(-1), M2(-1), D1(-1), D2(-1),
    121121  Charge(0), Mass(0.0),
    122   IsPU(0), IsRecoPU(0), IsConstituent(0),
    123   Flavor(0), FlavorAlgo(0), FlavorPhys(0),
    124   BTag(0), BTagAlgo(0), BTagPhys(0),
    125   TauTag(0), Eem(0.0), Ehad(0.0),
     122  IsPU(0), IsConstituent(0),
     123  BTag(0), TauTag(0), Eem(0.0), Ehad(0.0),
    126124  DeltaEta(0.0), DeltaPhi(0.0),
    127125  Momentum(0.0, 0.0, 0.0, 0.0),
    … …  
    131129  NCharged(0),
    132130  NNeutrals(0),
    133   NSubJetsTrimmed(0),
    134   NSubJetsPruned(0),
    135   NSubJetsSoftDropped(0),
    136131  Beta(0),
    137132  BetaStar(0),
    138133  MeanSqDeltaR(0),
    139134  PTD(0),
    140   Ntimes(-1),
    141   IsolationVar(-999),
    142   IsolationVarRhoCorr(-999),
    143   SumPtCharged(-999),
    144   SumPtNeutral(-999),
    145   SumPtChargedPU(-999),
    146   SumPt(-999),
    147135  fFactory(0),
    148136  fArray(0)
    149137{
    150   int i;
    151138  Edges[0] = 0.0;
    152139  Edges[1] = 0.0;
    … …  
    163150  Tau[3] = 0.0;
    164151  Tau[4] = 0.0;
    165   for(i = 0; i < 5; ++i)
    166   {
    167     TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    168     PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    169     SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    170   }
    171152}
    172153
    … …  
    243224  object.IsPU = IsPU;
    244225  object.IsConstituent = IsConstituent;
    245   object.Flavor = Flavor;
    246   object.FlavorAlgo = FlavorAlgo;
    247   object.FlavorPhys = FlavorPhys;
    248226  object.BTag = BTag;
    249   object.BTagAlgo = BTagAlgo;
    250   object.BTagPhys = BTagPhys;
    251227  object.TauTag = TauTag;
    252228  object.Eem = Eem;
    … …  
    273249  object.MeanSqDeltaR = MeanSqDeltaR;
    274250  object.PTD = PTD;
    275   object.Ntimes = Ntimes;
    276   object.IsolationVar = IsolationVar;
    277   object.IsolationVarRhoCorr = IsolationVarRhoCorr;
    278   object.SumPtCharged = SumPtCharged;
    279   object.SumPtNeutral = SumPtNeutral;
    280   object.SumPtChargedPU = SumPtChargedPU;
    281   object.SumPt = SumPt;
    282 
    283251  object.FracPt[0] = FracPt[0];
    284252  object.FracPt[1] = FracPt[1];
    … …  
    291259  object.Tau[3] = Tau[3];
    292260  object.Tau[4] = Tau[4];
    293  
    294   object.TrimmedP4[0] = TrimmedP4[0];
    295   object.TrimmedP4[1] = TrimmedP4[1];
    296   object.TrimmedP4[2] = TrimmedP4[2];
    297   object.TrimmedP4[3] = TrimmedP4[3];
    298   object.TrimmedP4[4] = TrimmedP4[4];
    299   object.PrunedP4[0] = PrunedP4[0];
    300   object.PrunedP4[1] = PrunedP4[1];
    301   object.PrunedP4[2] = PrunedP4[2];
    302   object.PrunedP4[3] = PrunedP4[3];
    303   object.PrunedP4[4] = PrunedP4[4];
    304   object.SoftDroppedP4[0] = SoftDroppedP4[0];
    305   object.SoftDroppedP4[1] = SoftDroppedP4[1];
    306   object.SoftDroppedP4[2] = SoftDroppedP4[2];
    307   object.SoftDroppedP4[3] = SoftDroppedP4[3];
    308   object.SoftDroppedP4[4] = SoftDroppedP4[4];
    309 
    310   object.NSubJetsTrimmed = NSubJetsTrimmed;
    311   object.NSubJetsPruned = NSubJetsPruned;
    312   object.NSubJetsSoftDropped = NSubJetsSoftDropped;
    313261
    314262  object.fFactory = fFactory;
    315263  object.fArray = 0;
    316 
    317   // Copy cluster timing info
    318   copy(Ecal_E_t.begin(),Ecal_E_t.end(),back_inserter(object.Ecal_E_t));
    319264
    320265  if(fArray && fArray->GetEntriesFast() > 0)
    … …  
    333278void Candidate::Clear(Option_t* option)
    334279{
    335   int i;
    336280  SetUniqueID(0);
    337281  ResetBit(kIsReferenced);
    … …  
    343287  IsPU = 0;
    344288  IsConstituent = 0;
    345   Flavor = 0;
    346   FlavorAlgo = 0;
    347   FlavorPhys = 0;
    348289  BTag = 0;
    349   BTagAlgo = 0;
    350   BTagPhys = 0;
    351290  TauTag = 0;
    352291  Eem = 0.0;
    … …  
    372311  MeanSqDeltaR = 0.0;
    373312  PTD = 0.0;
    374  
    375   Ntimes = 0;
    376   Ecal_E_t.clear();
    377  
    378   IsolationVar = -999;
    379   IsolationVarRhoCorr = -999;
    380   SumPtCharged = -999;
    381   SumPtNeutral = -999;
    382   SumPtChargedPU = -999;
    383   SumPt = -999;
    384  
    385313  FracPt[0] = 0.0;
    386314  FracPt[1] = 0.0;
    … …  
    393321  Tau[3] = 0.0;
    394322  Tau[4] = 0.0;
    395  
    396   for(i = 0; i < 5; ++i)
    397   {
    398     TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    399     PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    400     SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
    401   }
    402 
    403   NSubJetsTrimmed = 0;
    404   NSubJetsPruned = 0;
    405   NSubJetsSoftDropped = 0;
    406  
     323
    407324  fArray = 0;
    408325}

  • classes/DelphesClasses.h

    r9343566 re33e6db  
    8484//---------------------------------------------------------------------------
    8585
    86 class LHEFWeight: public TObject
    87 {
    88 public:
    89   Int_t ID; // weight ID
    90   Float_t Weight; // weight value
    91 
    92   ClassDef(LHEFWeight, 1)
    93 };
    94 
    95 //---------------------------------------------------------------------------
    96 
    9786class HepMCEvent: public Event
    9887{
    … …  
    242231  TRefArray Particles; // references to generated particles
    243232
    244   // Isolation variables
    245 
    246   Float_t IsolationVar;
    247   Float_t IsolationVarRhoCorr;
    248   Float_t SumPtCharged;
    249   Float_t SumPtNeutral;
    250   Float_t SumPtChargedPU;
    251   Float_t SumPt;
    252 
    253233  static CompBase *fgCompare; //!
    254234  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    277257  TRef Particle; // reference to generated particle
    278258
    279   // Isolation variables
    280 
    281   Float_t IsolationVar;
    282   Float_t IsolationVarRhoCorr;
    283   Float_t SumPtCharged;
    284   Float_t SumPtNeutral;
    285   Float_t SumPtChargedPU;
    286   Float_t SumPt;
    287 
    288259  static CompBase *fgCompare; //!
    289260  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    310281  TRef Particle; // reference to generated particle
    311282
    312    // Isolation variables
    313 
    314   Float_t IsolationVar;
    315   Float_t IsolationVarRhoCorr;
    316   Float_t SumPtCharged;
    317   Float_t SumPtNeutral;
    318   Float_t SumPtChargedPU;
    319   Float_t SumPt;
    320 
    321283  static CompBase *fgCompare; //!
    322284  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    344306  Float_t DeltaPhi;  // jet radius in azimuthal angle
    345307
    346   UInt_t Flavor;
    347   UInt_t FlavorAlgo;
    348   UInt_t FlavorPhys;
    349 
    350308  UInt_t BTag; // 0 or 1 for a jet that has been tagged as containing a heavy quark
    351   UInt_t BTagAlgo;
    352   UInt_t BTagPhys;
    353 
    354309  UInt_t TauTag; // 0 or 1 for a jet that has been tagged as a tau
    355310
    … …  
    358313  Float_t EhadOverEem; // ratio of the hadronic versus electromagnetic energy deposited in the calorimeter
    359314
    360   Int_t NCharged; // number of charged constituents
    361   Int_t NNeutrals; // number of neutral constituents
    362   Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
    363   Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents)
    364   Float_t MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
    365   Float_t PTD; // average pt between constituent and jet weighted by pt of constituent
    366   Float_t FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents)
    367 
    368   Float_t Tau[5]; // N-subjettiness
    369 
    370   TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
    371   TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
    372   TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
    373 
    374   Int_t NSubJetsTrimmed; // number of subjets trimmed
    375   Int_t NSubJetsPruned; // number of subjets pruned
    376   Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
     315  Int_t    NCharged; // number of charged constituents
     316  Int_t    NNeutrals; // number of neutral constituents
     317  Float_t  Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
     318  Float_t  BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents)
     319  Float_t  MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
     320  Float_t  PTD; // average pt between constituent and jet weighted by pt of constituent
     321  Float_t  FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents)
     322
     323  Float_t Tau1; // 1-subjettiness
     324  Float_t Tau2; // 2-subjettiness
     325  Float_t Tau3; // 3-subjettiness
     326  Float_t Tau4; // 4-subjettiness
     327  Float_t Tau5; // 5-subjettiness
    377328
    378329  TRefArray Constituents; // references to constituents
    … …  
    383334
    384335  TLorentzVector P4() const;
    385   TLorentzVector Area;
    386 
    387   ClassDef(Jet, 3)
     336
     337  ClassDef(Jet, 2)
    388338};
    389339
    … …  
    442392  Float_t E; // calorimeter tower energy
    443393
    444   Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
    445   Int_t   Ntimes; // number of hits contributing to time measurement
     394  Float_t T; //particle arrival time of flight
    446395
    447396  Float_t Eem; // calorimeter tower electromagnetic energy
    … …  
    503452
    504453  Int_t IsPU;
    505   Int_t IsRecoPU;
    506 
    507454  Int_t IsConstituent;
    508455
    509   UInt_t Flavor;
    510   UInt_t FlavorAlgo;
    511   UInt_t FlavorPhys;
    512 
    513456  UInt_t BTag;
    514   UInt_t BTagAlgo;
    515   UInt_t BTagPhys;
    516 
    517457  UInt_t TauTag;
    518458
    … …  
    542482  Float_t  FracPt[5];
    543483
    544   //Timing information
    545 
    546   Int_t    Ntimes;
    547   std::vector<std::pair<Float_t,Float_t> > Ecal_E_t;
    548 
    549   // Isolation variables
    550 
    551   Float_t IsolationVar;
    552   Float_t IsolationVarRhoCorr;
    553   Float_t SumPtCharged;
    554   Float_t SumPtNeutral;
    555   Float_t SumPtChargedPU;
    556   Float_t SumPt;
    557 
    558484  // N-subjettiness variables
    559485
    560486  Float_t Tau[5];
    561 
    562   // Other Substructure variables
    563 
    564   TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
    565   TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
    566   TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
    567 
    568   Int_t NSubJetsTrimmed; // number of subjets trimmed
    569   Int_t NSubJetsPruned; // number of subjets pruned
    570   Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
    571 
    572487
    573488  static CompBase *fgCompare; //!
    … …  
    589504  void SetFactory(DelphesFactory *factory) { fFactory = factory; }
    590505
    591   ClassDef(Candidate, 3)
     506  ClassDef(Candidate, 2)
    592507};
    593508

  • classes/DelphesFormula.cc

    r9343566 re33e6db  
    2323
    2424#include <stdexcept>
     25#include <string>
    2526
    2627using namespace std;
    … …  
    5051Int_t DelphesFormula::Compile(const char *expression)
    5152{
    52   TString buffer;
     53  string buffer;
    5354  const char *it;
    5455  for(it = expression; *it; ++it)
    5556  {
    5657    if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    57     buffer.Append(*it);
     58    buffer.push_back(*it);
    5859  }
    59   buffer.ReplaceAll("pt", "x");
    60   buffer.ReplaceAll("eta", "y");
    61   buffer.ReplaceAll("phi", "z");
    62   buffer.ReplaceAll("energy", "t");
    63   if(TFormula::Compile(buffer) != 0)
     60  if(TFormula::Compile(buffer.c_str()) != 0)
    6461  {
    6562    throw runtime_error("Invalid formula.");
    … …  
    7774
    7875//------------------------------------------------------------------------------
     76
     77Int_t DelphesFormula::DefinedVariable(TString &chaine, Int_t &action)
     78{
     79  action = kVariable;
     80  if(chaine == "pt")
     81  {
     82    if(fNdim < 1) fNdim = 1;
     83    return 0;
     84  }
     85  else if(chaine == "eta")
     86  {
     87    if(fNdim < 2) fNdim = 2;
     88    return 1;
     89  }
     90  else if(chaine == "phi")
     91  {
     92    if(fNdim < 3) fNdim = 3;
     93    return 2;
     94  }
     95  else if(chaine == "energy")
     96  {
     97    if(fNdim < 4) fNdim = 4;
     98    return 3;
     99  }
     100  return -1;
     101}
     102
     103//------------------------------------------------------------------------------

  • classes/DelphesFormula.h

    r9343566 re33e6db  
    3535
    3636  Double_t Eval(Double_t pt, Double_t eta = 0, Double_t phi = 0, Double_t energy = 0);
     37
     38  Int_t DefinedVariable(TString &variable, Int_t &action);
    3739};
    3840

  • classes/DelphesLHEFReader.cc

    r9343566 re33e6db  
    8282  fEventCounter = -1;
    8383  fParticleCounter = -1;
    84   fWeightList.clear();
     84  fRwgtList.clear();
    8585}
    8686
    … …  
    9999  TObjArray *partonOutputArray)
    100100{
    101   int rc, id;
     101  int rc;
    102102  char *pch;
    103103  double weight;
    … …  
    158158  else if(strstr(fBuffer, "<wgt"))
    159159  {
    160     pch = strpbrk(fBuffer, "\"'");
     160    pch = strstr(fBuffer, ">");
    161161    if(!pch)
    162162    {
    … …  
    165165    }
    166166
    167     DelphesStream idStream(pch + 1);
    168     rc = idStream.ReadInt(id);
    169 
    170     pch = strchr(fBuffer, '>');
    171     if(!pch)
     167    DelphesStream bufferStream(pch + 1);
     168    rc = bufferStream.ReadDbl(weight);
     169
     170    if(!rc)
    172171    {
    173172      cerr << "** ERROR: " << "invalid weight format" << endl;
    … …  
    175174    }
    176175
    177     DelphesStream weightStream(pch + 1);
    178     rc = weightStream.ReadDbl(weight);
    179 
    180     if(!rc)
    181     {
    182       cerr << "** ERROR: " << "invalid weight format" << endl;
    183       return kFALSE;
    184     }
    185 
    186     fWeightList.push_back(make_pair(id, weight));
     176    fRwgtList.push_back(weight);
    187177  }
    188178  else if(strstr(fBuffer, "</event>"))
    … …  
    216206//---------------------------------------------------------------------------
    217207
    218 void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch)
    219 {
    220   LHEFWeight *element;
    221   vector< pair< int, double > >::const_iterator itWeightList;
    222 
    223   for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList)
    224   {
    225     element = static_cast<LHEFWeight *>(branch->NewEntry());
    226 
    227     element->ID = itWeightList->first;
    228     element->Weight = itWeightList->second;
     208void DelphesLHEFReader::AnalyzeRwgt(ExRootTreeBranch *branch)
     209{
     210  Weight *element;
     211  vector<double>::const_iterator itRwgtList;
     212
     213  for(itRwgtList = fRwgtList.begin(); itRwgtList != fRwgtList.end(); ++itRwgtList)
     214  {
     215    element = static_cast<Weight *>(branch->NewEntry());
     216
     217    element->Weight = *itRwgtList;
    229218  }
    230219}

  • classes/DelphesLHEFReader.h

    r9343566 re33e6db  
    3131
    3232#include <vector>
    33 #include <utility>
    3433
    3534class TObjArray;
    … …  
    5958    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
    6059
    61   void AnalyzeWeight(ExRootTreeBranch *branch);
     60  void AnalyzeRwgt(ExRootTreeBranch *branch);
    6261
    6362private:
    … …  
    8483  double fPx, fPy, fPz, fE, fMass;
    8584 
    86   std::vector< std::pair< int, double > > fWeightList;
     85  std::vector<double> fRwgtList;
    8786};
    8887

  • classes/DelphesTF2.cc

    r9343566 re33e6db  
    1818
    1919#include "classes/DelphesTF2.h"
    20 
    21 #include "RVersion.h"
    2220#include "TString.h"
    23 
    2421#include <stdexcept>
     22#include <string>
    2523
    2624using namespace std;
    … …  
    3634
    3735DelphesTF2::DelphesTF2(const char *name, const char *expression) :
    38   TF2(name, expression)
     36  TF2(name,expression)
    3937{
    4038}
    … …  
    4846//------------------------------------------------------------------------------
    4947
    50 Int_t DelphesTF2::Compile(const char *expression)
     48Int_t DelphesTF2::DefinedVariable(TString &chaine, Int_t &action)
    5149{
    52   TString buffer;
    53   const char *it;
    54   for(it = expression; *it; ++it)
     50  action = kVariable;
     51  if(chaine == "z")
    5552  {
    56     if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    57     buffer.Append(*it);
     53    if(fNdim < 1) fNdim = 1;
     54    return 0;
    5855  }
    59   buffer.ReplaceAll("z", "x");
    60   buffer.ReplaceAll("t", "y");
    61 #if  ROOT_VERSION_CODE < ROOT_VERSION(6,04,00)
    62   if(TF2::Compile(buffer) != 0)
    63 #else
    64   if(TF2::GetFormula()->Compile(buffer) != 0)
    65 #endif
     56  else if(chaine == "t")
    6657  {
    67     throw runtime_error("Invalid formula.");
     58    if(fNdim < 2) fNdim = 2;
     59    return 1;
    6860  }
    69   return 0;
     61  return -1;
    7062}
    7163

  • classes/DelphesTF2.h

    r9343566 re33e6db  
    2121
    2222#include "TF2.h"
     23#include "TFormula.h"
     24
     25#include <string>
    2326
    2427class DelphesTF2: public TF2
    … …  
    3235  ~DelphesTF2();
    3336
    34   Int_t Compile(const char *expression);
     37  Int_t DefinedVariable(TString &variable, Int_t &action);
     38
    3539};
    3640
    3741#endif /* DelphesTF2_h */
     42

  • display/Delphes3DGeometry.cc

    r9343566 re33e6db  
    1717 */
    1818
     19#include "display/Delphes3DGeometry.h"
    1920#include <set>
    2021#include <map>
     22#include <string>
    2123#include <utility>
    2224#include <vector>
    … …  
    2426#include <sstream>
    2527#include <cassert>
    26 
    27 #include "TAxis.h"
    2828#include "TGeoManager.h"
    2929#include "TGeoVolume.h"
    … …  
    3535#include "TGeoCone.h"
    3636#include "TGeoArb8.h"
     37#include "external/ExRootAnalysis/ExRootConfReader.h"
     38#include "classes/DelphesClasses.h"
    3739#include "TF2.h"
    3840#include "TFormula.h"
    3941#include "TH1F.h"
    4042#include "TMath.h"
    41 #include "TString.h"
    42 
    43 #include "display/Delphes3DGeometry.h"
    44 
    45 #include "classes/DelphesClasses.h"
    46 #include "external/ExRootAnalysis/ExRootConfReader.h"
    4743
    4844using namespace std;
    … …  
    9793   tk_Bz_     = confReader->GetDouble("ParticlePropagator::Bz", 0.0);                           // tk_Bz
    9894   
    99    TString buffer;
     95   string buffer;
    10096   const char *it;
    10197 
    … …  
    107103   tkEffFormula.ReplaceAll("phi","0.");
    108104 
    109    buffer.Clear();
    110105   for(it = tkEffFormula.Data(); *it; ++it)
    111106   {
    112107     if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    113      buffer.Append(*it);
    114    }
    115 
    116    TF2* tkEffFunction = new TF2("tkEff",buffer,0,1000,-10,10);
     108     buffer.push_back(*it);
     109   }
     110
     111   TF2* tkEffFunction = new TF2("tkEff",buffer.c_str(),0,1000,-10,10);
    117112   TH1F etaHisto("eta","eta",100,5.,-5.);
    118113   Double_t pt,eta;
    … …  
    137132   muonEffFormula.ReplaceAll("phi","0.");
    138133   
    139    buffer.Clear();
     134   buffer.clear();
    140135   for(it = muonEffFormula.Data(); *it; ++it)
    141136   {
    142137     if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    143      buffer.Append(*it);
    144    }
    145 
    146    TF2* muEffFunction = new TF2("muEff",buffer,0,1000,-10,10);
     138     buffer.push_back(*it);
     139   }
     140
     141   TF2* muEffFunction = new TF2("muEff",buffer.c_str(),0,1000,-10,10);
    147142   TH1F etaHisto("eta2","eta2",100,5.,-5.);
    148143   Double_t pt,eta;

  • display/Delphes3DGeometry.h

    r9343566 re33e6db  
    11/*
    2   * Delphes: a framework for fast simulation of a generic collider experiment
    3   * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
     2 * Delphes: a framework for fast simulation of a generic collider experiment
     3 * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    44 *
    5   * This program is free software: you can redistribute it and/or modify
    6   * it under the terms of the GNU General Public License as published by
    7   * the Free Software Foundation, either version 3 of the License, or
    8   * (at your option) any later version.
     5 * This program is free software: you can redistribute it and/or modify
     6 * it under the terms of the GNU General Public License as published by
     7 * the Free Software Foundation, either version 3 of the License, or
     8 * (at your option) any later version.
    99 *
    10   * This program is distributed in the hope that it will be useful,
    11   * but WITHOUT ANY WARRANTY; without even the implied warranty of
    12   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    13   * GNU General Public License for more details.
     10 * This program is distributed in the hope that it will be useful,
     11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
     12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13 * GNU General Public License for more details.
    1414 *
    15   * You should have received a copy of the GNU General Public License
    16   * along with this program.  If not, see <http://www.gnu.org/licenses/>.
     15 * You should have received a copy of the GNU General Public License
     16 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
    1717 */
    1818
    … …  
    2222#include <set>
    2323#include <map>
     24#include <utility>
    2425#include <vector>
    25 
    26 #include "Rtypes.h"
    27 
    28 class TAxis;
    29 class TGeoManager;
    30 class TGeoVolume;
    31 class TGeoMedium;
     26#include <algorithm>
     27#include <sstream>
     28#include "TAxis.h"
     29#include "TGeoManager.h"
     30#include "TGeoVolume.h"
     31#include "TGeoMedium.h"
    3232
    3333// TODO: asymmetric detector
    … …  
    3838     ~Delphes3DGeometry() {}
    3939
    40      void readFile(const char *filename, const char *ParticlePropagator="ParticlePropagator",
    41                                          const char *TrackingEfficiency="ChargedHadronTrackingEfficiency",
    42                                          const char *MuonEfficiency="MuonEfficiency",
    43                                          const char *Calorimeters="Calorimeter");
     40     void readFile(const char* filename, const char* ParticlePropagator="ParticlePropagator",
     41                                         const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
     42                                         const char* MuonEfficiency="MuonEfficiency",
     43                                         const char* Calorimeters="Calorimeter");
    4444
    4545     void setContingency(Double_t contingency) { contingency_ = contingency; }
    … …  
    4848     void setMuonSystemThickness(Double_t thickness) { muonSystem_thickn_ = thickness; }
    4949
    50      TGeoVolume *getDetector(bool withTowers = true);
     50     TGeoVolume* getDetector(bool withTowers = true);
    5151
    5252     Double_t getTrackerRadius() const { return tk_radius_; }
    … …  
    5959   private:
    6060     std::pair<Double_t, Double_t> addTracker(TGeoVolume *top);
    61      std::pair<Double_t, Double_t> addCalorimeter(TGeoVolume *top, const char *name, Double_t innerBarrelRadius, Double_t innerBarrelLength, std::set< std::pair<Double_t, Int_t> >& caloBinning);
    62      std::pair<Double_t, Double_t> addMuonDets(TGeoVolume *top, const char *name, Double_t innerBarrelRadius, Double_t innerBarrelLength);
    63      void addCaloTowers(TGeoVolume *top, const char *name, Double_t innerBarrelRadius, Double_t innerBarrelLength, std::set< std::pair<Double_t, Int_t> >& caloBinning);
     61     std::pair<Double_t, Double_t> addCalorimeter(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, std::set< std::pair<Double_t, Int_t> >& caloBinning);
     62     std::pair<Double_t, Double_t> addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength);
     63     void addCaloTowers(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, std::set< std::pair<Double_t, Int_t> >& caloBinning);
    6464
    6565   private:
    … …  
    7272     TGeoMedium *mudetmed_;
    7373
    74      TAxis *etaAxis_;
    75      TAxis *phiAxis_;
     74     TAxis* etaAxis_;
     75     TAxis* phiAxis_;
    7676
    7777     Double_t contingency_;
    … …  
    9191     std::map<std::string, Double_t> muonSystem_etamax_;
    9292     std::map<std::string, std::set< std::pair<Double_t, Int_t> > > caloBinning_;
    93 
     93     
    9494};
    9595

  • display/DelphesEventDisplay.cc

    r9343566 re33e6db  
    11/*
    2   * Delphes: a framework for fast simulation of a generic collider experiment
    3   * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
     2 * Delphes: a framework for fast simulation of a generic collider experiment
     3 * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    44 *
    5   * This program is free software: you can redistribute it and/or modify
    6   * it under the terms of the GNU General Public License as published by
    7   * the Free Software Foundation, either version 3 of the License, or
    8   * (at your option) any later version.
     5 * This program is free software: you can redistribute it and/or modify
     6 * it under the terms of the GNU General Public License as published by
     7 * the Free Software Foundation, either version 3 of the License, or
     8 * (at your option) any later version.
    99 *
    10   * This program is distributed in the hope that it will be useful,
    11   * but WITHOUT ANY WARRANTY; without even the implied warranty of
    12   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    13   * GNU General Public License for more details.
     10 * This program is distributed in the hope that it will be useful,
     11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
     12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13 * GNU General Public License for more details.
    1414 *
    15   * You should have received a copy of the GNU General Public License
    16   * along with this program.  If not, see <http://www.gnu.org/licenses/>.
     15 * You should have received a copy of the GNU General Public License
     16 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
    1717 */
    1818
    … …  
    2121#include <utility>
    2222#include <algorithm>
    23 
    2423#include "TGeoManager.h"
    2524#include "TGeoVolume.h"
     25#include "external/ExRootAnalysis/ExRootConfReader.h"
     26#include "external/ExRootAnalysis/ExRootTreeReader.h"
     27#include "display/DelphesCaloData.h"
     28#include "display/DelphesBranchElement.h"
     29#include "display/Delphes3DGeometry.h"
     30#include "display/DelphesEventDisplay.h"
     31#include "classes/DelphesClasses.h"
    2632#include "TEveElement.h"
    2733#include "TEveJetCone.h"
    … …  
    4753#include "TCanvas.h"
    4854#include "TH1F.h"
    49 #include "TAxis.h"
    50 #include "TChain.h"
    51 #include "TGHtml.h"
    52 #include "TGStatusBar.h"
    53 
    54 #include "display/DelphesCaloData.h"
    55 #include "display/DelphesBranchElement.h"
    56 #include "display/Delphes3DGeometry.h"
    57 #include "display/DelphesEventDisplay.h"
    58 #include "display/DelphesDisplay.h"
    59 #include "display/Delphes3DGeometry.h"
    60 #include "display/DelphesHtmlSummary.h"
    61 #include "display/DelphesPlotSummary.h"
    62 
    63 #include "classes/DelphesClasses.h"
    64 #include "external/ExRootAnalysis/ExRootConfReader.h"
    65 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    6655
    6756DelphesEventDisplay::DelphesEventDisplay()
    … …  
    10998   TEveManager::Create(kTRUE, "IV");
    11099   fStatusBar_ = gEve->GetBrowser()->GetStatusBar();
    111    TGeoManager *geom = gGeoManager;
     100   TGeoManager* geom = gGeoManager;
    112101
    113102   // build the detector
    … …  
    119108   etaAxis_ = det3D.getCaloAxes().first;
    120109   phiAxis_ = det3D.getCaloAxes().second;
    121    TGeoVolume *top = det3D.getDetector(false);
     110   TGeoVolume* top = det3D.getDetector(false);
    122111   geom->SetTopVolume(top);
    123112   TEveElementList *geometry = new TEveElementList("Geometry");
    124    TObjArray *nodes = top->GetNodes();
     113   TObjArray* nodes = top->GetNodes();
    125114   TIter itNodes(nodes);
    126    TGeoNode *nodeobj;
    127    TEveGeoTopNode *node;
     115   TGeoNode* nodeobj;
     116   TEveGeoTopNode* node;
    128117   while((nodeobj = (TGeoNode*)itNodes.Next())) {
    129118     node = new TEveGeoTopNode(gGeoManager,nodeobj);
    … …  
    142131   // prepare data collections
    143132   readConfig(configFile, elements_);
    144    for(std::vector<DelphesBranchBase *>::iterator element = elements_.begin(); element<elements_.end(); ++element) {
    145      DelphesBranchElement<TEveTrackList> *item_v1 = dynamic_cast<DelphesBranchElement<TEveTrackList>*>(*element);
    146      DelphesBranchElement<TEveElementList> *item_v2 = dynamic_cast<DelphesBranchElement<TEveElementList>*>(*element);
     133   for(std::vector<DelphesBranchBase*>::iterator element = elements_.begin(); element<elements_.end(); ++element) {
     134     DelphesBranchElement<TEveTrackList>*   item_v1 = dynamic_cast<DelphesBranchElement<TEveTrackList>*>(*element);
     135     DelphesBranchElement<TEveElementList>* item_v2 = dynamic_cast<DelphesBranchElement<TEveElementList>*>(*element);
    147136     if(item_v1) gEve->AddElement(item_v1->GetContainer());
    148137     if(item_v2) gEve->AddElement(item_v2->GetContainer());
    … …  
    155144   delphesDisplay_->ImportGeomRhoZ(geometry);
    156145   // find the first calo data and use that to initialize the calo display
    157    for(std::vector<DelphesBranchBase *>::iterator data=elements_.begin();data<elements_.end();++data) {
     146   for(std::vector<DelphesBranchBase*>::iterator data=elements_.begin();data<elements_.end();++data) {
    158147     if(TString((*data)->GetType())=="Tower") { // we could also use GetClassName()=="DelphesCaloData"
    159        DelphesCaloData *container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>((*data))->GetContainer();
     148       DelphesCaloData* container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>((*data))->GetContainer();
    160149       assert(container);
    161150       TEveCalo3D *calo3d = new TEveCalo3D(container);
    … …  
    193182   ExRootConfParam branches = confReader->GetParam("TreeWriter::Branch");
    194183   Int_t nBranches = branches.GetSize()/3;
    195    DelphesBranchElement<TEveTrackList> *tlist;
    196    DelphesBranchElement<DelphesCaloData> *clist;
    197    DelphesBranchElement<TEveElementList> *elist;
     184   DelphesBranchElement<TEveTrackList>* tlist;
     185   DelphesBranchElement<DelphesCaloData>* clist;
     186   DelphesBranchElement<TEveElementList>* elist;
    198187   // first loop with all but tracks
    199188   for(Int_t b = 0; b<nBranches; ++b) {
    … …  
    284273
    285274   // update display
    286    TEveElement *top = (TEveElement*)gEve->GetCurrentEvent();
     275   TEveElement* top = (TEveElement*)gEve->GetCurrentEvent();
    287276   delphesDisplay_->DestroyEventRPhi();
    288277   delphesDisplay_->ImportEventRPhi(top);
    … …  
    367356
    368357   // add a tab on the left
    369    TEveBrowser *browser = gEve->GetBrowser();
     358   TEveBrowser* browser = gEve->GetBrowser();
    370359   browser->SetWindowName("Delphes Event Display");
    371360   browser->StartEmbedding(TRootBrowser::kLeft);
    372361
    373362   // set the main title
    374    TGMainFrame *frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
     363   TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
    375364   frmMain->SetWindowName("Delphes Event Display");
    376365   frmMain->SetCleanup(kDeepCleanup);
    … …  
    382371   if(!gSystem->OpenDirectory(icondir))
    383372     icondir = Form("%s/icons/", (const char*)gSystem->GetFromPipe("root-config --etcdir") );
    384    TGGroupFrame *vf = new TGGroupFrame(frmMain,"Event navigation",kVerticalFrame | kFitWidth );
     373   TGGroupFrame* vf = new TGGroupFrame(frmMain,"Event navigation",kVerticalFrame | kFitWidth );
    385374   {
    386      TGHorizontalFrame *hf = new TGHorizontalFrame(frmMain);
     375     TGHorizontalFrame* hf = new TGHorizontalFrame(frmMain);
    387376     {
    388         TGPictureButton *b = 0;
     377        TGPictureButton* b = 0;
    389378
    390379        b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif"));
    … …  
    392381        b->Connect("Clicked()", "DelphesEventDisplay", this, "Bck()");
    393382
    394         TGNumberEntry *numberEntry = new TGNumberEntry(hf,0,9,-1,TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, 0, treeReader_->GetEntries());
     383        TGNumberEntry* numberEntry = new TGNumberEntry(hf,0,9,-1,TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, 0, treeReader_->GetEntries());
    395384        hf->AddFrame(numberEntry, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY , 2, 0, 10, 10));
    396385        this->Connect("EventChanged(Int_t)","TGNumberEntry",numberEntry,"SetIntNumber(Long_t)");
    … …  
    405394     vf->AddFrame(hf, new TGLayoutHints(kLHintsExpandX , 2, 2, 2, 2));
    406395
    407      TGHProgressBar *progress = new TGHProgressBar(frmMain, TGProgressBar::kFancy, 100);
     396     TGHProgressBar* progress = new TGHProgressBar(frmMain, TGProgressBar::kFancy, 100);
    408397     progress->SetMax( treeReader_->GetEntries());
    409398     progress->ShowPosition(kTRUE, kFALSE, "Event %.0f");
    … …  
    429418   // the summary tab
    430419   htmlSummary_ = new DelphesHtmlSummary("Delphes Event Display Summary Table");
    431    TEveWindowSlot *slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
     420   TEveWindowSlot* slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
    432421   gHtml_ = new TGHtml(0, 100, 100);
    433422   TEveWindowFrame *wf = slot->MakeFrame(gHtml_);
    … …  
    437426   // plot tab
    438427   slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
    439    TEveWindowTab *tab = slot->MakeTab();
     428   TEveWindowTab* tab = slot->MakeTab();
    440429   tab->SetElementName("Summary plots");
    441430   tab->SetShowTitleBar(kFALSE);
    … …  
    446435}
    447436
    448 void DelphesEventDisplay::Fwd() { 
    449   if (event_id_ < treeReader_->GetEntries() - 2) {
    450      EventChanged(event_id_+1);
    451   } else {
    452      printf("Already at last event.\n");
    453   }
    454 }
    455 
    456 void DelphesEventDisplay::Bck() {
    457   if (event_id_ > 0) {
    458      EventChanged(event_id_-1);
    459   } else {
    460      printf("Already at first event.\n");
    461   }
    462 }
    463 
    464 void DelphesEventDisplay::PreSetEv(char *ev) {
    465   event_id_tmp_ = Int_t(atoi(ev));
    466 }
    467 
    468 void DelphesEventDisplay::GoTo() {
    469   if (event_id_tmp_>=0 && event_id_tmp_ < treeReader_->GetEntries()-1) {
    470     EventChanged(event_id_tmp_);
    471   } else {
    472     printf("Error: no such event.\n");
    473   }
    474 }
    475 
    476 void DelphesEventDisplay::InitSummaryPlots() {
    477   plotSummary_->FillSample(treeReader_, event_id_);
    478   plotSummary_->FillEvent();
    479   plotSummary_->Draw();
    480 }
    481 
    482 void DelphesEventDisplay::DisplayProgress(Int_t p) {
    483   fStatusBar_->SetText(Form("Processing... %d %%",p), 1);
    484   gSystem->ProcessEvents();
    485 }
     437

  • display/DelphesEventDisplay.h

    r9343566 re33e6db  
    11/*
    2   * Delphes: a framework for fast simulation of a generic collider experiment
    3   * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
     2 * Delphes: a framework for fast simulation of a generic collider experiment
     3 * Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    44 *
    5   * This program is free software: you can redistribute it and/or modify
    6   * it under the terms of the GNU General Public License as published by
    7   * the Free Software Foundation, either version 3 of the License, or
    8   * (at your option) any later version.
     5 * This program is free software: you can redistribute it and/or modify
     6 * it under the terms of the GNU General Public License as published by
     7 * the Free Software Foundation, either version 3 of the License, or
     8 * (at your option) any later version.
    99 *
    10   * This program is distributed in the hope that it will be useful,
    11   * but WITHOUT ANY WARRANTY; without even the implied warranty of
    12   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    13   * GNU General Public License for more details.
     10 * This program is distributed in the hope that it will be useful,
     11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
     12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13 * GNU General Public License for more details.
    1414 *
    15   * You should have received a copy of the GNU General Public License
    16   * along with this program.  If not, see <http://www.gnu.org/licenses/>.
     15 * You should have received a copy of the GNU General Public License
     16 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
    1717 */
    1818
    … …  
    2121
    2222#include <vector>
     23#include <iostream>
     24#include "external/ExRootAnalysis/ExRootTreeReader.h"
     25#include "display/DelphesDisplay.h"
     26#include "display/Delphes3DGeometry.h"
     27#include "display/DelphesHtmlSummary.h"
     28#include "display/DelphesPlotSummary.h"
     29#include "TSystem.h"
     30#include "TChain.h"
     31#include "TAxis.h"
     32#include "TGHtml.h"
     33#include "TClonesArray.h"
     34#include "TGStatusBar.h"
     35#include "TGNumberEntry.h"
     36#include "TGProgressBar.h"
     37#include <RQ_OBJECT.h>
    2338
    24 #include "Rtypes.h"
    25 #include "RQ_OBJECT.h"
    2639
    27 class TAxis;
    28 class TChain;
    29 class TGHtml;
    30 class TGStatusBar;
    31 class DelphesDisplay;
    32 class Delphes3DGeometry;
    33 class DelphesBranchBase;
    34 class DelphesHtmlSummary;
    35 class DelphesPlotSummary;
    36 class ExRootTreeReader;
    3740
    3841/*
    39   *assembly.C: sauvegarde as shape-extract -> implement in the geometry class (read/write)
    40   *histobrowser.C: intégration d'histogrammes dans le display (on pourrait avoir Pt, eta, phi pour les principales collections)
    41   *also from alice_esd: summary html table
    42   *
     42 * assembly.C: sauvegarde as shape-extract -> implement in the geometry class (read/write)
     43 * histobrowser.C: intégration d'histogrammes dans le display (on pourrait avoir Pt, eta, phi pour les principales collections)
     44 * also from alice_esd: summary html table
     45 *
    4346 */
    4447
    … …  
    5659    void make_gui();
    5760    void load_event();
    58     void readConfig(const char *configFile, std::vector<DelphesBranchBase *>& elements);
     61    void readConfig(const char *configFile, std::vector<DelphesBranchBase*>& elements);
    5962
    6063    // Configuration and global variables.
    … …  
    6467    Double_t tkRadius_, totRadius_, tkHalfLength_, muHalfLength_, bz_;
    6568    TAxis *etaAxis_, *phiAxis_;
    66     TChain *chain_;
    67     std::vector<DelphesBranchBase *> elements_;
     69    TChain* chain_;
     70    std::vector<DelphesBranchBase*> elements_;
    6871    DelphesDisplay *delphesDisplay_;
    6972    DelphesHtmlSummary *htmlSummary_;
    7073    TGHtml *gHtml_;
    7174    DelphesPlotSummary *plotSummary_;
    72     TGStatusBar *fStatusBar_;
     75    TGStatusBar* fStatusBar_;
    7376   
    7477    // gui controls
    7578  public:
    76      void Fwd();
     79     void Fwd() { 
     80        if (event_id_ < treeReader_->GetEntries() - 2) {
     81           EventChanged(event_id_+1);
     82        } else {
     83           printf("Already at last event.\n");
     84        }
     85     }
    7786
    78      void Bck();
     87     void Bck() {
     88        if (event_id_ > 0) {
     89           EventChanged(event_id_-1);
     90        } else {
     91           printf("Already at first event.\n");
     92        }
     93     }
    7994
    80     void PreSetEv(char *ev);
     95    void PreSetEv(char* ev) {
     96      event_id_tmp_ = Int_t(atoi(ev));
     97    }
    8198
    82     void GoTo();
     99    void GoTo() {
     100      if (event_id_tmp_>=0 && event_id_tmp_ < treeReader_->GetEntries()-1) {
     101        EventChanged(event_id_tmp_);
     102      } else {
     103        printf("Error: no such event.\n");
     104      }
     105    }
    83106
    84     void InitSummaryPlots();
     107    void InitSummaryPlots() {
     108      plotSummary_->FillSample(treeReader_, event_id_);
     109      plotSummary_->FillEvent();
     110      plotSummary_->Draw();
     111    }
    85112
    86     void DisplayProgress(Int_t p);
     113    void DisplayProgress(Int_t p) {
     114         fStatusBar_->SetText(Form("Processing... %d %%",p), 1);
     115         gSystem->ProcessEvents();
     116    }
    87117};
    88118

  • doc/genMakefile.tcl

    r9343566 re33e6db  
    283283dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h}
    284284
    285 sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} {external/PUPPI/*.cc}
     285sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc}
    286286
    287287sourceDeps {FASTJET} {modules/FastJet*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/fastjet/contribs/*/*.cc}

  • examples/EventDisplay.C

    r9343566 re33e6db  
    33 * root -l examples/EventDisplay.C'("cards/delphes_card_FCC_basic.tcl","delphes_output.root","ParticlePropagator","ChargedHadronTrackingEfficiency","MuonTrackingEfficiency","Ecal,Hcal")'
    44 */
    5 
    6 #ifdef __CLING__
    7 R__LOAD_LIBRARY(libEve)
    8 R__LOAD_LIBRARY(libDelphesDisplay)
    9 #include "display/DelphesEventDisplay.h"
    10 #include "display/Delphes3DGeometry.h"
    11 #endif
    125
    136void EventDisplay(const char *configfile = "delphes_card_CMS.tcl",
    … …  
    4033
    4134    // create the application
    42     DelphesEventDisplay *display = new DelphesEventDisplay(configfile, datafile, det3D);
     35    DelphesEventDisplay* display = new DelphesEventDisplay(configfile, datafile, det3D);
    4336  }
    4437}

  • examples/Example1.C

    r9343566 re33e6db  
    66root -l examples/Example1.C'("delphes_output.root")'
    77*/
    8 
    9 #ifdef __CLING__
    10 R__LOAD_LIBRARY(libDelphes)
    11 #include "classes/DelphesClasses.h"
    12 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    13 #endif
    148
    159//------------------------------------------------------------------------------

  • examples/Example2.C

    r9343566 re33e6db  
    88#include "TH1.h"
    99#include "TSystem.h"
    10 
    11 #ifdef __CLING__
    12 R__LOAD_LIBRARY(libDelphes)
    13 #include "classes/DelphesClasses.h"
    14 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    15 #include "external/ExRootAnalysis/ExRootResult.h"
    16 #endif
    1710
    1811//------------------------------------------------------------------------------

  • examples/Example3.C

    r9343566 re33e6db  
    66*/
    77
    8 #ifdef __CLING__
    9 R__LOAD_LIBRARY(libDelphes)
    10 #include "classes/DelphesClasses.h"
    11 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    12 #include "external/ExRootAnalysis/ExRootResult.h"
    13 #else
    14 class ExRootTreeReader;
    15 class ExRootResult;
    16 #endif
    17 
    188//------------------------------------------------------------------------------
    199
    … …  
    3525  TH1 *fJetDeltaPT;
    3626};
     27
     28//------------------------------------------------------------------------------
     29
     30class ExRootResult;
     31class ExRootTreeReader;
    3732
    3833//------------------------------------------------------------------------------

  • examples/Example4.C

    r9343566 re33e6db  
    11/*
     2
    23This macro shows how to compute jet energy scale.
    34root -l examples/Example4.C'("delphes_output.root", "plots.root")'
    45
    5 The output ROOT file contains the pT(MC)/pT(Reco) distributions for
    6 various pT(Reco) and |eta| bins. The peak value of such distribution is
    7 interpreted as the jet energy correction to be applied for that
    8 given pT(Reco), |eta| bin.
    9 
    10 This can be done by modifying the "ScaleFormula" input parameter to
    11 the JetEnergyScale module in the delphes_card_XXX.tcl
     6The output ROOT file contains the pT(MC)/pT(Reco) distributions for various pT(Reco) and |eta| bins.
     7The peak value of such distribution is interpreted as the jet energy correction to be applied for that given pT(Reco), |eta| bin.
     8
     9This can be done by modifying the "ScaleFormula" input parameter to the JetEnergyScale module in the delphes_card_XXX.tcl
     10
     11
    1212
    1313e.g  a smooth function:
    1414
    15   set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
     15
     16  set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0 ) }
     17
    1618
    1719or a binned function:
     20
    1821
    1922  set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0)  * (1.10) +
    … …  
    2427                    (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0)               * (1.00)}
    2528
    26 Be aware that a binned jet energy scale can produce "steps" in the corrected
    27 jet pt distribution ...
     29
     30Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ...
     31
     32
     33
    2834*/
    29 
    30 #ifdef __CLING__
    31 R__LOAD_LIBRARY(libDelphes)
    32 #include "classes/DelphesClasses.h"
    33 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    34 #include "external/ExRootAnalysis/ExRootResult.h"
    35 #else
    36 class ExRootTreeReader;
    37 class ExRootResult;
    38 #endif
    3935
    4036//------------------------------------------------------------------------------
    … …  
    165161
    166162  Jet *jet, *genjet;
    167   GenParticle *particle;
     163  GenParticle *part;
    168164  TObject *object;
    169165
    170   TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
    171 
    172   Float_t deltaR;
     166  TLorentzVector JetMom, GenJetMom, BestGenJetMom;
     167
     168  Float_t Dr;
    173169  Float_t pt, eta;
    174170  Long64_t entry;
    … …  
    181177    // Load selected branches with data from specified event
    182178    treeReader->ReadEntry(entry);
     179    //  cout<<"--  New event -- "<<endl;
    183180
    184181    if(entry%500 == 0) cout << "Event number: "<< entry <<endl;
    … …  
    189186
    190187      jet = (Jet*) branchJet->At(i);
    191       jetMomentum = jet->P4();
    192 
    193       plots->fJetPT->Fill(jetMomentum.Pt());
    194 
    195       deltaR = 999;
     188      JetMom = jet-> P4();
     189
     190      plots->fJetPT->Fill(JetMom.Pt());
     191
     192      Dr = 999;
    196193
    197194     // Loop over all hard partons in event
    198195     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
    199196     {
    200         particle = (GenParticle*) branchParticle->At(j);
    201 
    202         genJetMomentum = particle->P4();
    203 
    204         // this is simply to avoid warnings from initial state particle
    205         // having infite rapidity ...
    206         if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
    207 
    208         // take the closest parton candidate
    209         if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
     197
     198        part = (GenParticle*) branchParticle->At(j);
     199
     200        GenJetMom = part -> P4();
     201
     202        //this is simply to avoid warnings from initial state particle having infite rapidity ...
     203        if(GenJetMom.Px() == 0 && GenJetMom.Py() == 0) continue;
     204
     205        //take the closest parton candidate
     206        if( GenJetMom.DeltaR(JetMom) < Dr )
    210207        {
    211            deltaR = genJetMomentum.DeltaR(jetMomentum);
    212            bestGenJetMomentum = genJetMomentum;
     208           Dr = GenJetMom.DeltaR(JetMom);
     209           BestGenJetMom = GenJetMom;
    213210        }
     211
    214212      }
    215213
    216       if(deltaR < 0.3)
    217       {
    218         pt  = jetMomentum.Pt();
    219         eta = TMath::Abs(jetMomentum.Eta());
    220 
    221 
    222         if(pt > 20.0 && pt < 50.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_20_50_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    223         if(pt > 20.0 && pt < 50.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_20_50_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    224 
    225         if(pt > 50.0 && pt < 100.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_50_100_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    226         if(pt > 50.0 && pt < 100.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_50_100_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    227 
    228         if(pt > 100.0 && pt < 200.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_100_200_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    229         if(pt > 100.0 && pt < 200.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_100_200_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    230 
    231         if(pt > 200.0 && pt < 500.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_200_500_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    232         if(pt > 200.0 && pt < 500.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_200_500_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    233 
    234         if(pt > 500.0               && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_500_inf_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    235         if(pt > 500.0               && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_500_inf_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
    236 
    237       }
     214     if(Dr < 0.3)
     215     {
     216       pt  = JetMom.Pt();
     217       eta = TMath::Abs(JetMom.Eta());
     218
     219
     220       if( pt > 20.0 && pt < 50.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_20_50_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     221       if( pt > 20.0 && pt < 50.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_20_50_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     222
     223       if( pt > 50.0 && pt < 100.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_50_100_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     224       if( pt > 50.0 && pt < 100.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_50_100_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     225
     226       if( pt > 100.0 && pt < 200.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_100_200_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     227       if( pt > 100.0 && pt < 200.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_100_200_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     228
     229       if( pt > 200.0 && pt < 500.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_200_500_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     230       if( pt > 200.0 && pt < 500.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_200_500_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     231
     232       if( pt > 500.0               && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_500_inf_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     233       if( pt > 500.0               && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_500_inf_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
     234
     235
     236     }
     237
     238
    238239    }
    239240  }
    240241}
     242
    241243
    242244//------------------------------------------------------------------------------

  • modules/BTagging.cc

    r9343566 re33e6db  
    2222 *  Determines origin of jet,
    2323 *  applies b-tagging efficiency (miss identification rate) formulas
    24  *  and sets b-tagging flags
     24 *  and sets b-tagging flags 
    2525 *
    2626 *  \author P. Demin - UCL, Louvain-la-Neuve
    … …  
    4646#include "TLorentzVector.h"
    4747
    48 #include <algorithm>
     48#include <algorithm> 
    4949#include <stdexcept>
    5050#include <iostream>
    … …  
    6262
    6363  Int_t GetCategory(TObject *object);
    64 
     64 
    6565  Double_t fEtaMax, fPTMin;
    6666};
    … …  
    7575
    7676  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
    77 
     77 
    7878  pdgCode = TMath::Abs(parton->PID);
    7979  if(pdgCode != 21 && pdgCode > 5) return -1;
    … …  
    117117  param = GetParam("EfficiencyFormula");
    118118  size = param.GetSize();
    119 
     119 
    120120  fEfficiencyMap.clear();
    121121  for(i = 0; i < size/2; ++i)
    … …  
    143143
    144144  fFilter = new ExRootFilter(fPartonInputArray);
    145 
     145 
    146146  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
    147147  fItJetInputArray = fJetInputArray->MakeIterator();
    … …  
    170170void BTagging::Process()
    171171{
    172   Candidate *jet;
     172  Candidate *jet, *parton;
    173173  Double_t pt, eta, phi, e;
    174174  TObjArray *partonArray;
    175175  map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
    176176  DelphesFormula *formula;
     177  Int_t pdgCode, pdgCodeMax;
    177178
    178179  // select quark and gluons
    179180  fFilter->Reset();
    180181  partonArray = fFilter->GetSubArray(fClassifier, 0);
    181 
     182 
    182183  if(partonArray == 0) return;
    183184
    184185  TIter itPartonArray(partonArray);
    185 
     186 
    186187  // loop over all input jets
    187188  fItJetInputArray->Reset();
    … …  
    189190  {
    190191    const TLorentzVector &jetMomentum = jet->Momentum;
     192    pdgCodeMax = -1;
    191193    eta = jetMomentum.Eta();
    192194    phi = jetMomentum.Phi();
    193195    pt = jetMomentum.Pt();
    194196    e = jetMomentum.E();
     197   
     198    // loop over all input partons
     199    itPartonArray.Reset();
     200    while((parton = static_cast<Candidate*>(itPartonArray.Next())))
     201    {
     202      pdgCode = TMath::Abs(parton->PID);
     203      if(pdgCode == 21) pdgCode = 0;
     204      if(jetMomentum.DeltaR(parton->Momentum) <= fDeltaR)
     205      {
     206        if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
     207      }
     208    }
     209    if(pdgCodeMax == 0) pdgCodeMax = 21;
     210    if(pdgCodeMax == -1) pdgCodeMax = 0;
    195211
    196212    // find an efficency formula
    197     itEfficiencyMap = fEfficiencyMap.find(jet->Flavor);
     213    itEfficiencyMap = fEfficiencyMap.find(pdgCodeMax);
    198214    if(itEfficiencyMap == fEfficiencyMap.end())
    199215    {
    … …  
    204220    // apply an efficency formula
    205221    jet->BTag |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
    206 
    207     // find an efficency formula for algo flavor definition
    208     itEfficiencyMap = fEfficiencyMap.find(jet->FlavorAlgo);
    209     if(itEfficiencyMap == fEfficiencyMap.end())
    210     {
    211       itEfficiencyMap = fEfficiencyMap.find(0);
    212     }
    213     formula = itEfficiencyMap->second;
    214 
    215     // apply an efficency formula
    216     jet->BTagAlgo |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
    217 
    218     // find an efficency formula for phys flavor definition
    219     itEfficiencyMap = fEfficiencyMap.find(jet->FlavorPhys);
    220     if(itEfficiencyMap == fEfficiencyMap.end())
    221     {
    222       itEfficiencyMap = fEfficiencyMap.find(0);
    223     }
    224     formula = itEfficiencyMap->second;
    225 
    226     // apply an efficency formula
    227     jet->BTagPhys |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
    228   }
    229 }
    230 
    231 //------------------------------------------------------------------------------
     222  }
     223}
     224
     225//------------------------------------------------------------------------------

  • modules/Calorimeter.cc

    r9343566 re33e6db  
    150150*/
    151151
    152   // read min E value for timing measurement in ECAL
    153   fTimingEMin = GetDouble("TimingEMin",4.);
    154   // For timing
    155   // So far this flag needs to be false
    156   // Curved extrapolation not supported
    157   fElectronsFromTrack = false;
    158 
    159  
    160152  // read min E value for towers to be saved
    161153  fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
    … …  
    364356      fTrackHCalEnergy = 0.0;
    365357
     358      fTowerECalTime = 0.0;
     359      fTowerHCalTime = 0.0;
     360
     361      fTrackECalTime = 0.0;
     362      fTrackHCalTime = 0.0;
     363
     364      fTowerECalTimeWeight = 0.0;
     365      fTowerHCalTimeWeight = 0.0;
     366
    366367      fTowerTrackHits = 0;
    367368      fTowerPhotonHits = 0;
    … …  
    379380      position = track->Position;
    380381
     382
    381383      ecalEnergy = momentum.E() * fTrackECalFractions[number];
    382384      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
    … …  
    385387      fTrackHCalEnergy += hcalEnergy;
    386388
    387       bool dbg_scz = false;
    388       if (dbg_scz) {
    389         cout << "   Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T()
    390              << endl;
    391         Candidate *prt = static_cast<Candidate*>(track->GetCandidates()->Last());
    392         const TLorentzVector &ini = prt->Position;
    393 
    394         cout << "                and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T();
    395 
    396       }
    397      
    398       if (ecalEnergy > fTimingEMin && fTower) {
    399         if (fElectronsFromTrack) {
    400           //      cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
    401           //           << " PID=" << track->PID << endl;
    402           fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T()));
    403         } else {
    404           //      cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
    405           //           << " PID=" << track->PID << endl;
    406         }
    407       }
    408 
     389      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     390      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     391
     392      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
     393      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
    409394
    410395      fTowerTrackArray->Add(track);
    … …  
    412397      continue;
    413398    }
    414  
     399
    415400    // check for photon and electron hits in current tower
    416401    if(flags & 2) ++fTowerPhotonHits;
    … …  
    427412    fTowerHCalEnergy += hcalEnergy;
    428413
    429     if (ecalEnergy > fTimingEMin && fTower) {
    430       if (abs(particle->PID) != 11 || !fElectronsFromTrack) {
    431         //      cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
    432         //   << " PID=" << particle->PID << endl;
    433         fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T()));
    434       } else {
    435        
    436         // N.B. Only charged particle set to leave ecal energy is the electrons
    437         //      cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
    438         //           << " PID=" << particle->PID << endl;
    439        
    440       }
    441     }
     414    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     415    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     416
     417    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
     418    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
     419
    442420
    443421    fTower->AddCandidate(particle);
    … …  
    456434  Double_t ecalEnergy, hcalEnergy;
    457435  Double_t ecalSigma, hcalSigma;
    458  
     436  Double_t ecalTime, hcalTime, time;
     437
    459438  if(!fTower) return;
    460439
    … …  
    465444  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
    466445
     446  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
     447  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
     448
    467449  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    468450  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    … …  
    472454
    473455  energy = ecalEnergy + hcalEnergy;
    474    
     456  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
     457
    475458  if(fSmearTowerCenter)
    476459  {
    … …  
    486469  pt = energy / TMath::CosH(eta);
    487470
    488   // Time calculation for tower
    489   fTower->Ntimes = 0;
    490   Float_t tow_sumT = 0;
    491   Float_t tow_sumW = 0;
    492  
    493   for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++)
    494   {
    495     Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first);
    496     tow_sumT += w*fTower->Ecal_E_t[i].second;
    497     tow_sumW += w;
    498     fTower->Ntimes++;
    499   }
    500  
    501   if (tow_sumW > 0.) {
    502     fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW);
    503   } else {
    504     fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.);
    505   }
    506 
    507 
     471  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
    508472  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    509473  fTower->Eem = ecalEnergy;

  • modules/Calorimeter.h

    r9343566 re33e6db  
    6060  Double_t fTrackECalEnergy, fTrackHCalEnergy;
    6161
    62   Double_t fTimingEMin;
    63   Bool_t fElectronsFromTrack; // for timing
    64  
     62  Double_t fTowerECalTime, fTowerHCalTime;
     63  Double_t fTrackECalTime, fTrackHCalTime;
     64
     65  Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
     66  Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
     67
    6568  Int_t fTowerTrackHits, fTowerPhotonHits;
    6669

  • modules/FastJetFinder.cc

    r9343566 re33e6db  
    6565#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
    6666#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
    67 
    68 #include "fastjet/tools/Filter.hh"
    69 #include "fastjet/tools/Pruner.hh"
    70 #include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
    7167
    7268using namespace std;
    … …  
    124120  fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
    125121  fN = GetInt("N", 2);                  // used only if Njettiness is used as jet clustering algo (case 8)
    126 
    127   //-- Trimming parameters --
    128  
    129   fComputeTrimming = GetBool("ComputeTrimming", false);
    130   fRTrim = GetDouble("RTrim", 0.2);
    131   fPtFracTrim = GetDouble("PtFracTrim", 0.05);
    132  
    133 
    134   //-- Pruning parameters --
    135  
    136   fComputePruning = GetBool("ComputePruning", false);
    137   fZcutPrun = GetDouble("ZcutPrun", 0.1);
    138   fRcutPrun = GetDouble("RcutPrun", 0.5);
    139   fRPrun = GetDouble("RPrun", 0.8);
    140  
    141   //-- SoftDrop parameters --
    142  
    143   fComputeSoftDrop     = GetBool("ComputeSoftDrop", false);
    144   fBetaSoftDrop        = GetDouble("BetaSoftDrop", 0.0);
    145   fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
    146   fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
    147  
    148122
    149123  // ---  Jet Area Parameters ---
    … …  
    286260  PseudoJet jet, area;
    287261  ClusterSequence *sequence;
    288   vector< PseudoJet > inputList, outputList, subjets;
     262  vector< PseudoJet > inputList, outputList;
    289263  vector< PseudoJet >::iterator itInputList, itOutputList;
    290264  vector< TEstimatorStruct >::iterator itEstimators;
    … …  
    378352    candidate->DeltaEta = detaMax;
    379353    candidate->DeltaPhi = dphiMax;
    380        
    381     //------------------------------------
    382     // Trimming
    383     //------------------------------------
    384 
    385      if(fComputeTrimming)
    386     {
    387 
    388       fastjet::Filter    trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
    389       fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
    390      
    391       trimmed_jet = join(trimmed_jet.constituents());
    392      
    393       candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
    394        
    395       // four hardest subjets
    396       subjets.clear();
    397       subjets = trimmed_jet.pieces();
    398       subjets = sorted_by_pt(subjets);
    399      
    400       candidate->NSubJetsTrimmed = subjets.size();
    401 
    402       for (size_t i = 0; i < subjets.size() and i < 4; i++){
    403         if(subjets.at(i).pt() < 0) continue ;
    404         candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    405       }
    406     }
    407    
    408    
    409     //------------------------------------
    410     // Pruning
    411     //------------------------------------
    412    
    413    
    414     if(fComputePruning)
    415     {
    416 
    417       fastjet::Pruner    pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
    418       fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
    419 
    420       candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
    421          
    422       // four hardest subjet
    423       subjets.clear();
    424       subjets = pruned_jet.pieces();
    425       subjets = sorted_by_pt(subjets);
    426      
    427       candidate->NSubJetsPruned = subjets.size();
    428 
    429       for (size_t i = 0; i < subjets.size() and i < 4; i++){
    430         if(subjets.at(i).pt() < 0) continue ;
    431         candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    432       }
    433 
    434     }
    435      
    436     //------------------------------------
    437     // SoftDrop
    438     //------------------------------------
    439    
    440     if(fComputeSoftDrop)
    441     {
    442    
    443       contrib::SoftDrop  softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
    444       fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
    445      
    446       candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
    447        
    448       // four hardest subjet
    449      
    450       subjets.clear();
    451       subjets    = softdrop_jet.pieces();
    452       subjets    = sorted_by_pt(subjets);
    453       candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
    454 
    455       for (size_t i = 0; i < subjets.size()  and i < 4; i++){
    456         if(subjets.at(i).pt() < 0) continue ;
    457         candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    458       }
    459     }
    460  
     354
    461355    // --- compute N-subjettiness with N = 1,2,3,4,5 ----
    462356

  • modules/FastJetFinder.h

    r9343566 re33e6db  
    8383  Int_t fN ;
    8484
    85   //-- Trimming parameters --
    86  
    87   Bool_t fComputeTrimming;
    88   Double_t fRTrim;
    89   Double_t fPtFracTrim;
    90  
    91   //-- Pruning parameters --
    92 
    93   Bool_t fComputePruning;
    94   Double_t fZcutPrun;
    95   Double_t fRcutPrun;
    96   Double_t fRPrun;
    97 
    98   //-- SoftDrop parameters --
    99 
    100   Bool_t fComputeSoftDrop;
    101   Double_t fBetaSoftDrop;
    102   Double_t fSymmetryCutSoftDrop;
    103   Double_t fR0SoftDrop;
    104 
    10585  // --- FastJet Area method --------
    10686

  • modules/Isolation.cc

    r9343566 re33e6db  
    2525 *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
    2626 *
    27  *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
     27 *  \author P. Demin - UCL, Louvain-la-Neuve
    2828 *
    2929 */
    … …  
    109109  fUsePTSum = GetBool("UsePTSum", false);
    110110
    111   fVetoLeptons = GetBool("VetoLeptons", true); 
    112  
    113111  fClassifier->fPTMin = GetDouble("PTMin", 0.5);
    114 
    115112
    116113  // import input array(s)
    … …  
    156153  Candidate *candidate, *isolation, *object;
    157154  TObjArray *isolationArray;
    158   Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
     155  Double_t sum, ratio;
    159156  Int_t counter;
    160157  Double_t eta = 0.0;
    … …  
    183180
    184181    // loop over all input tracks
    185    
    186     sumNeutral = 0.0;
    187     sumCharged = 0.0;
    188     sumChargedPU = 0.0;
    189     sumAllParticles = 0.0;
    190    
     182    sum = 0.0;
    191183    counter = 0;
    192184    itIsolationArray.Reset();
    193    
    194185    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
    195186    {
    … …  
    197188
    198189      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
    199          candidate->GetUniqueID() != isolation->GetUniqueID() &&
    200          ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) )
     190         candidate->GetUniqueID() != isolation->GetUniqueID())
    201191      {
    202      
    203         sumAllParticles += isolationMomentum.Pt();
    204         if(isolation->Charge !=0)
    205         {
    206           sumCharged += isolationMomentum.Pt();
    207           if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt();
    208         }
    209  
    210         else sumNeutral += isolationMomentum.Pt();
    211  
     192        sum += isolationMomentum.Pt();
    212193        ++counter;
    213194      }
    … …  
    228209    }
    229210
    230        
    231      // correct sum for pile-up contamination
    232     sumDBeta = sumCharged + TMath::Max(sumNeutral-0.5*sumChargedPU,0.0);
    233     sumRhoCorr = sumCharged + TMath::Max(sumNeutral-TMath::Max(rho,0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(),0.0);
    234     ratioDBeta = sumDBeta/candidateMomentum.Pt();
    235     ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
    236    
    237     candidate->IsolationVar = ratioDBeta;
    238     candidate->IsolationVarRhoCorr = ratioRhoCorr;
    239     candidate->SumPtCharged = sumCharged;
    240     candidate->SumPtNeutral = sumNeutral;
    241     candidate->SumPtChargedPU = sumChargedPU;
    242     candidate->SumPt = sumAllParticles;
    243 
    244     if((fUsePTSum && sumDBeta > fPTSumMax) || ratioDBeta > fPTRatioMax) continue;
     211    // correct sum for pile-up contamination
     212    sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
     213
     214    ratio = sum/candidateMomentum.Pt();
     215    if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
     216
    245217    fOutputArray->Add(candidate);
    246218  }

  • modules/Isolation.h

    r9343566 re33e6db  
    5959  Bool_t fUsePTSum;
    6060
    61   Bool_t fVetoLeptons; 
    62  
    6361  IsolationClassifier *fClassifier; //!
    6462

  • modules/ModulesLinkDef.h

    r9343566 re33e6db  
    5858#include "modules/Weighter.h"
    5959#include "modules/Hector.h"
    60 #include "modules/RunPUPPI.h"
    61 #include "modules/JetFlavorAssociation.h"
    6260#include "modules/ExampleModule.h"
    6361
    … …  
    10098#pragma link C++ class Weighter+;
    10199#pragma link C++ class Hector+;
    102 #pragma link C++ class RunPUPPI+;
    103 #pragma link C++ class JetFlavorAssociation+;
    104100#pragma link C++ class ExampleModule+;
    105101

  • modules/PdgCodeFilter.cc

    r9343566 re33e6db  
    7474  fPTMin = GetDouble("PTMin", 0.0);
    7575
    76   fInvertPdg = GetBool("InvertPdg", false);
    77 
    78   fRequireStatus = GetBool("RequireStatus", false);
    79   fStatus = GetInt("Status", 1);
    80 
    8176  // import input array
    8277  fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
    … …  
    114109  Double_t pt;
    115110
    116   const Bool_t requireStatus = fRequireStatus;
    117   const Bool_t invertPdg = fInvertPdg;
    118   const int reqStatus = fStatus;
    119 
    120111  fItInputArray->Reset();
    121112  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    … …  
    125116    pt = candidateMomentum.Pt();
    126117
    127     if(pt < fPTMin) continue;
    128     if(requireStatus && (candidate->Status != reqStatus)) continue;
     118    pass = kTRUE;
    129119
    130     pass = kTRUE;
     120    if(pt < fPTMin) pass = kFALSE;
    131121    if(find(fPdgCodes.begin(), fPdgCodes.end(), pdgCode) != fPdgCodes.end()) pass = kFALSE;
    132122
    133     if (invertPdg) pass = !pass;
    134123    if(pass) fOutputArray->Add(candidate);
    135124  }

  • modules/PdgCodeFilter.h

    r9343566 re33e6db  
    5050
    5151  Double_t fPTMin; //!
    52   Bool_t fInvertPdg; //!
    53   Bool_t fRequireStatus; //!
    54   Int_t fStatus; //!
    5552
    5653  std::vector<Int_t> fPdgCodes;

  • modules/PileUpJetID.cc

    r9343566 re33e6db  
     1/*
     2 *  Delphes: a framework for fast simulation of a generic collider experiment
     3 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
     4 *
     5 *  This program is free software: you can redistribute it and/or modify
     6 *  it under the terms of the GNU General Public License as published by
     7 *  the Free Software Foundation, either version 3 of the License, or
     8 *  (at your option) any later version.
     9 *
     10 *  This program is distributed in the hope that it will be useful,
     11 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
     12 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13 *  GNU General Public License for more details.
     14 *
     15 *  You should have received a copy of the GNU General Public License
     16 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
     17 */
     18
    119
    220/** \class PileUpJetID
    321 *
    4  *  CMS PileUp Jet ID Variables
    5  *
    6  *  \author S. Zenz
     22 *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
     23 *
     24 *  \author S. Zenz, December 2013
    725 *
    826 */
    … …  
    2341#include "TRandom3.h"
    2442#include "TObjArray.h"
    25 //#include "TDatabasePDG.h"
     43#include "TDatabasePDG.h"
    2644#include "TLorentzVector.h"
    2745
    … …  
    3654
    3755PileUpJetID::PileUpJetID() :
    38   fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
     56  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
    3957{
    4058
    … …  
    5270void PileUpJetID::Init()
    5371{
    54 
    5572  fJetPTMin = GetDouble("JetPTMin", 20.0);
    5673  fParameterR = GetDouble("ParameterR", 0.5);
    5774  fUseConstituents = GetInt("UseConstituents", 0);
    5875
    59 
    60   /*
    61   Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
    62   Double_t fBetaMinBarrel; // |eta| < 2.5
    63   Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
    64   Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
    65   Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
    66   */
    67 
    68   fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
    69   fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
    70   fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
    71   fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
    72   fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
    73   fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
    74   fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
    75 
    76 
    77 
    78   cout << "  set MeanSqDeltaRMaxBarrel " << fMeanSqDeltaRMaxBarrel << endl;
    79   cout << "  set BetaMinBarrel " << fBetaMinBarrel << endl;
    80   cout << "  set MeanSqDeltaRMaxEndcap " << fMeanSqDeltaRMaxEndcap << endl;
    81   cout << "  set BetaMinEndcap " << fBetaMinEndcap << endl;
    82   cout << "  set MeanSqDeltaRMaxForward " << fMeanSqDeltaRMaxForward << endl;
    83 
    84 
    85 
    8676  fAverageEachTower = false; // for timing
    8777
    … …  
    9181  fItJetInputArray = fJetInputArray->MakeIterator();
    9282
    93 
    94   //  cout << "BeforE SCZ additions in init" << endl;
    95   //  cout << GetString("TrackInputArray", "ParticlePropagator/tracks") << endl;
    96   //  cout << GetString("EFlowTrackInputArray", "ParticlePropagator/tracks") << endl;
    97 
    98   fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
     83  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
    9984  fItTrackInputArray = fTrackInputArray->MakeIterator();
    10085
    101   fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
     86  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
    10287  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    10388
     89  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
     90  fItVertexInputArray = fVertexInputArray->MakeIterator();
     91
     92  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    10493
    10594  // create output array(s)
    10695
    10796  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
    108 
    109   fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
    110 
    111 
    112   //  cout << " end of INIT " << endl;
    113 
    11497}
    11598
    … …  
    118101void PileUpJetID::Finish()
    119102{
    120   //  cout << "In finish" << endl;
    121 
    122103  if(fItJetInputArray) delete fItJetInputArray;
    123104  if(fItTrackInputArray) delete fItTrackInputArray;
    124105  if(fItNeutralInputArray) delete fItNeutralInputArray;
    125 
     106  if(fItVertexInputArray) delete fItVertexInputArray;
    126107}
    127108
    … …  
    130111void PileUpJetID::Process()
    131112{
    132   //  cout << "start of process" << endl;
    133 
    134113  Candidate *candidate, *constituent;
    135114  TLorentzVector momentum, area;
    136 
    137   //  cout << "BeforE SCZ additions in process" << endl;
    138 
    139   // SCZ
    140   Candidate *trk;
     115  Int_t i, nc, nn;
     116  Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
     117  Double_t dr, pt, pt_ann[5];
     118  Double_t zvtx = 0.0;
     119
     120  Candidate *track;
     121
     122  // find z position of primary vertex
     123
     124  fItVertexInputArray->Reset();
     125  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
     126  {
     127    if(!candidate->IsPU)
     128    {
     129      zvtx = candidate->Position.Z();
     130      break;
     131    }
     132  }
    141133
    142134  // loop over all input candidates
    … …  
    147139    area = candidate->Area;
    148140
    149     float sumT0 = 0.;
    150     float sumT1 = 0.;
    151     float sumT10 = 0.;
    152     float sumT20 = 0.;
    153     float sumT30 = 0.;
    154     float sumT40 = 0.;
    155     float sumWeightsForT = 0.;
    156     candidate->Ntimes = 0;
    157 
    158     float sumpt = 0.;
    159     float sumptch = 0.;
    160     float sumptchpv = 0.;
    161     float sumptchpu = 0.;
    162     float sumdrsqptsq = 0.;
    163     float sumptsq = 0.;
    164     int nc = 0;
    165     int nn = 0;
    166     float pt_ann[5];
    167 
    168     for (int i = 0 ; i < 5 ; i++) {
    169       pt_ann[i] = 0.;
    170     }
    171 
    172     if (fUseConstituents) {
     141    sumpt = 0.0;
     142    sumptch = 0.0;
     143    sumptchpv = 0.0;
     144    sumptchpu = 0.0;
     145    sumdrsqptsq = 0.0;
     146    sumptsq = 0.0;
     147    nc = 0;
     148    nn = 0;
     149
     150    for(i = 0; i < 5; ++i)
     151    {
     152      pt_ann[i] = 0.0;
     153    }
     154
     155    if(fUseConstituents)
     156    {
    173157      TIter itConstituents(candidate->GetCandidates());
    174       while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
    175         float pt = constituent->Momentum.Pt();
    176         float dr = candidate->Momentum.DeltaR(constituent->Momentum);
    177         //      cout << " There exists a constituent with dr=" << dr << endl;
    178         sumpt += pt;
    179         sumdrsqptsq += dr*dr*pt*pt;
    180         sumptsq += pt*pt;
    181         if (constituent->Charge == 0) {
    182           nn++;
    183         } else {
    184           if (constituent->IsRecoPU) {
    185             sumptchpu += pt;
    186           } else {
    187             sumptchpv += pt;
    188           }
    189           sumptch += pt;
    190           nc++;
    191         }
    192         for (int i = 0 ; i < 5 ; i++) {
    193           if (dr > 0.1*i && dr < 0.1*(i+1)) {
    194             pt_ann[i] += pt;
    195           }
    196         }
    197         float tow_sumT = 0;
    198         float tow_sumW = 0;
    199         for (int i = 0 ; i < constituent->Ecal_E_t.size() ; i++) {
    200           float w = TMath::Sqrt(constituent->Ecal_E_t[i].first);
    201           if (fAverageEachTower) {
    202             tow_sumT += w*constituent->Ecal_E_t[i].second;
    203             tow_sumW += w;
    204           } else {
    205             sumT0 += w*constituent->Ecal_E_t[i].second;
    206             sumT1 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.001);
    207             sumT10 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.010);
    208             sumT20 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.020);
    209             sumT30 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.030);
    210             sumT40 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.040);
    211             sumWeightsForT += w;
    212             candidate->Ntimes++;
    213           }
    214         }
    215         if (fAverageEachTower && tow_sumW > 0.) {
    216           sumT0 += tow_sumT;
    217           sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
    218           sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
    219           sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
    220           sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
    221           sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
    222           sumWeightsForT += tow_sumW;
    223           candidate->Ntimes++;
    224         }
    225       }
    226     } else {
     158      while((constituent = static_cast<Candidate*>(itConstituents.Next())))
     159      {
     160        pt = constituent->Momentum.Pt();
     161        dr = candidate->Momentum.DeltaR(constituent->Momentum);
     162        sumpt += pt;
     163        sumdrsqptsq += dr*dr*pt*pt;
     164        sumptsq += pt*pt;
     165        if(constituent->Charge == 0)
     166        {
     167          // neutrals
     168          ++nn;
     169        }
     170        else
     171        {
     172          // charged
     173          if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
     174          {
     175            sumptchpu += pt;
     176          }
     177          else
     178          {
     179            sumptchpv += pt;
     180          }
     181          sumptch += pt;
     182          ++nc;
     183        }
     184        for(i = 0; i < 5; ++i)
     185        {
     186          if(dr > 0.1*i && dr < 0.1*(i + 1))
     187          {
     188            pt_ann[i] += pt;
     189          }
     190        }
     191      }
     192    }
     193    else
     194    {
    227195      // Not using constituents, using dr
    228196      fItTrackInputArray->Reset();
    229       while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
    230         if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
    231           float pt = trk->Momentum.Pt();
    232           sumpt += pt;
    233           sumptch += pt;
    234           if (trk->IsRecoPU) {
    235             sumptchpu += pt;
    236           } else {
    237             sumptchpv += pt;
    238           }
    239           float dr = candidate->Momentum.DeltaR(trk->Momentum);
    240           sumdrsqptsq += dr*dr*pt*pt;
    241           sumptsq += pt*pt;
    242           nc++;
    243           for (int i = 0 ; i < 5 ; i++) {
    244             if (dr > 0.1*i && dr < 0.1*(i+1)) {
    245               pt_ann[i] += pt;
    246             }
    247           }
    248         }
    249       }
     197      while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
     198      {
     199        if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
     200        {
     201          pt = track->Momentum.Pt();
     202          sumpt += pt;
     203          sumptch += pt;
     204          if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
     205          {
     206            sumptchpu += pt;
     207          }
     208          else
     209          {
     210            sumptchpv += pt;
     211          }
     212          dr = candidate->Momentum.DeltaR(track->Momentum);
     213          sumdrsqptsq += dr*dr*pt*pt;
     214          sumptsq += pt*pt;
     215          nc++;
     216          for(i = 0; i < 5; ++i)
     217          {
     218            if(dr > 0.1*i && dr < 0.1*(i + 1))
     219            {
     220              pt_ann[i] += pt;
     221            }
     222          }
     223        }
     224      }
     225
    250226      fItNeutralInputArray->Reset();
    251       while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
    252         if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
    253           float pt = constituent->Momentum.Pt();
    254           sumpt += pt;
    255           float dr = candidate->Momentum.DeltaR(constituent->Momentum);
    256           sumdrsqptsq += dr*dr*pt*pt;
    257           sumptsq += pt*pt;
    258           nn++;
    259           for (int i = 0 ; i < 5 ; i++) {
    260             if (dr > 0.1*i && dr < 0.1*(i+1)) {
    261             pt_ann[i] += pt;
    262             }
    263           }
    264         }
    265       }
    266     }
    267 
    268     if (sumptch > 0.) {
    269       candidate->Beta = sumptchpv/sumptch;
    270       candidate->BetaStar = sumptchpu/sumptch;
    271     } else {
    272       candidate->Beta = -999.;
    273       candidate->BetaStar = -999.;
    274     }
    275     if (sumptsq > 0.) {
     227      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
     228      {
     229        if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
     230        {
     231          pt = constituent->Momentum.Pt();
     232          sumpt += pt;
     233          dr = candidate->Momentum.DeltaR(constituent->Momentum);
     234          sumdrsqptsq += dr*dr*pt*pt;
     235          sumptsq += pt*pt;
     236          nn++;
     237          for(i = 0; i < 5; ++i)
     238          {
     239            if(dr > 0.1*i && dr < 0.1*(i + 1))
     240            {
     241              pt_ann[i] += pt;
     242            }
     243          }
     244        }
     245      }
     246    }
     247
     248    if(sumptch > 0.0)
     249    {
     250      candidate->Beta = sumptchpu/sumptch;
     251      candidate->BetaStar = sumptchpv/sumptch;
     252    }
     253    else
     254    {
     255      candidate->Beta = -999.0;
     256      candidate->BetaStar = -999.0;
     257    }
     258    if(sumptsq > 0.0)
     259    {
    276260      candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
    277     } else {
    278       candidate->MeanSqDeltaR = -999.;
     261    }
     262    else
     263    {
     264      candidate->MeanSqDeltaR = -999.0;
    279265    }
    280266    candidate->NCharged = nc;
    281267    candidate->NNeutrals = nn;
    282     if (sumpt > 0.) {
     268    if(sumpt > 0.0)
     269    {
    283270      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
    284       for (int i = 0 ; i < 5 ; i++) {
     271      for(i = 0; i < 5; ++i)
     272      {
    285273        candidate->FracPt[i] = pt_ann[i]/sumpt;
    286274      }
    287     } else {
    288       candidate->PTD = -999.;
    289       for (int i = 0 ; i < 5 ; i++) {
    290         candidate->FracPt[i] = -999.;
     275    }
     276    else
     277    {
     278      candidate->PTD = -999.0;
     279      for(i = 0; i < 5; ++i)
     280      {
     281        candidate->FracPt[i] = -999.0;
    291282      }
    292283    }
    293284
    294285    fOutputArray->Add(candidate);
    295 
    296     // New stuff
    297     /*
    298     fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
    299     fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
    300     fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
    301     fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
    302     fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
    303     */
    304 
    305     bool passId = false;
    306     if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
    307       if (fabs(candidate->Momentum.Eta())<1.5) {
    308         passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
    309       } else if (fabs(candidate->Momentum.Eta())<4.0) {
    310         passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
    311       } else {
    312         passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
    313       }
    314     }
    315 
    316     //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
    317     //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
    318 
    319     if (passId) {
    320       if (fUseConstituents) {
    321         TIter itConstituents(candidate->GetCandidates());
    322         while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
    323           if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
    324             fNeutralsInPassingJets->Add(constituent);
    325             //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
    326           }
    327         }
    328       } else { // use DeltaR
    329         fItNeutralInputArray->Reset();
    330         while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
    331           if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
    332             fNeutralsInPassingJets->Add(constituent);
    333             //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
    334           }
    335         }
    336       }
    337     }
    338 
    339 
    340286  }
    341287}
    342288
    343289//------------------------------------------------------------------------------
     290

  • modules/PileUpJetID.h

    r9343566 re33e6db  
     1/*
     2 *  Delphes: a framework for fast simulation of a generic collider experiment
     3 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
     4 *
     5 *  This program is free software: you can redistribute it and/or modify
     6 *  it under the terms of the GNU General Public License as published by
     7 *  the Free Software Foundation, either version 3 of the License, or
     8 *  (at your option) any later version.
     9 *
     10 *  This program is distributed in the hope that it will be useful,
     11 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
     12 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     13 *  GNU General Public License for more details.
     14 *
     15 *  You should have received a copy of the GNU General Public License
     16 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
     17 */
     18
    119#ifndef PileUpJetID_h
    220#define PileUpJetID_h
    … …  
    422/** \class PileUpJetID
    523 *
    6  *  CMS PileUp Jet ID Variables
     24 *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
    725 *
    8  *  \author S. Zenz
     26 *  \author S. Zenz, December 2013
    927 *
    1028 */
    … …  
    1634
    1735class TObjArray;
    18 class DelphesFormula;
    1936
    2037class PileUpJetID: public DelphesModule
    … …  
    3451  Double_t fParameterR;
    3552
    36   Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
    37   Double_t fBetaMinBarrel; // |eta| < 2.5
    38   Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
    39   Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
    40   Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
    41 
    42   Double_t fNeutralPTMin;
    43   Double_t fJetPTMinForNeutrals;
    44 
    45   /*
    46 JAY
    47 ---
    48 
    49 |Eta|<1.5
    50 
    51 meanSqDeltaR betaStar SigEff BgdEff
    52 0.13 0.92 96% 8%
    53 0.13 0.95 97% 16%
    54 0.13 0.97 98% 27%
    55 
    56 |Eta|>1.5
    57 
    58 meanSqDeltaR betaStar SigEff BgdEff
    59 0.14 0.91 95% 15%
    60 0.14 0.94 97% 19%
    61 0.14 0.97 98% 29%
    62 
    63 BRYAN
    64 -----
    65 
    66 Barrel (MeanSqDR, Beta, sig eff, bg eff):
    67 0.10, 0.08, 90%, 8%
    68 0.11, 0.12, 90%, 6%
    69 0.13, 0.16, 89%, 5%
    70 
    71 Endcap (MeanSqDR, Beta, sig eff, bg eff):
    72 0.07, 0.06, 89%, 4%
    73 0.08, 0.08, 92%, 6%
    74 0.09, 0.08, 95%, 10%
    75 0.10, 0.08, 97%, 13%
    76 
    77 SETH GUESSES FOR |eta| > 4.0
    78 ----------------------------
    79 
    80 MeanSqDeltaR
    81 0.07
    82 0.10
    83 0.14
    84 0.2
    85   */
    86 
    8753  // If set to true, may have weird results for PFCHS
    8854  // If set to false, uses everything within dR < fParameterR even if in other jets &c.
    8955  // Results should be very similar for PF
    90   Int_t fUseConstituents; 
     56  Int_t fUseConstituents;
    9157
    9258  Bool_t fAverageEachTower;
    … …  
    9662  const TObjArray *fJetInputArray; //!
    9763
    98   const TObjArray *fTrackInputArray; // SCZ
    99   const TObjArray *fNeutralInputArray;
     64  const TObjArray *fTrackInputArray; //!
     65  const TObjArray *fNeutralInputArray; //!
    10066
    101   TIterator *fItTrackInputArray; // SCZ
    102   TIterator *fItNeutralInputArray; // SCZ
     67  TIterator *fItTrackInputArray; //!
     68  TIterator *fItNeutralInputArray; //!
    10369
    10470  TObjArray *fOutputArray; //!
    105   TObjArray *fNeutralsInPassingJets; // SCZ
    10671
     72  TIterator *fItVertexInputArray; //!
     73  const TObjArray *fVertexInputArray; //!
    10774
    108   ClassDef(PileUpJetID, 2)
     75  Double_t fZVertexResolution;
     76
     77  ClassDef(PileUpJetID, 1)
    10978};
    11079
    11180#endif
     81

  • modules/PileUpMerger.cc

    r9343566 re33e6db  
    8080  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
    8181
    82   fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
    83   fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
    84   fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
    85   fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
    86 
    8782  // read vertex smearing formula
    8883
    … …  
    116111  TParticlePDG *pdgParticle;
    117112  Int_t pid;
    118   Float_t x, y, z, t, vx, vy;
     113  Float_t x, y, z, t;
    119114  Float_t px, py, pz, e;
    120115  Double_t dz, dphi, dt;
    121   Int_t numberOfEvents, event, numberOfParticles;
     116  Int_t numberOfEvents, event;
    122117  Long64_t allEntries, entry;
    123   Candidate *candidate, *vertex;
     118  Candidate *candidate, *vertexcandidate;
    124119  DelphesFactory *factory;
    125120
    … …  
    128123  fItInputArray->Reset();
    129124
    130   // --- Deal with primary vertex first  ------
     125  // --- Deal with Primary vertex first  ------
    131126
    132127  fFunction->GetRandom2(dz, dt);
    … …  
    134129  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    135130  dz *= 1.0E3; // necessary in order to make z in mm
    136   vx = 0.0;
    137   vy = 0.0;
    138   numberOfParticles = fInputArray->GetEntriesFast();
     131
    139132  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    140133  {
    141     vx += candidate->Position.X();
    142     vy += candidate->Position.Y();
    143134    z = candidate->Position.Z();
    144135    t = candidate->Position.T();
    … …  
    148139  }
    149140
    150   if(numberOfParticles > 0)
    151   {
    152     vx /= numberOfParticles;
    153     vy /= numberOfParticles;
    154   }
    155 
    156141  factory = GetFactory();
    157142
    158   vertex = factory->NewCandidate();
    159   vertex->Position.SetXYZT(vx, vy, dz, dt);
    160   fVertexOutputArray->Add(vertex);
     143  vertexcandidate = factory->NewCandidate();
     144  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
     145  fVertexOutputArray->Add(vertexcandidate);
    161146
    162147  // --- Then with pile-up vertices  ------
    … …  
    196181    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    197182
    198     vx = 0.0;
    199     vy = 0.0;
    200     numberOfParticles = 0;
     183    vertexcandidate = factory->NewCandidate();
     184    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
     185    vertexcandidate->IsPU = 1;
     186
     187    fVertexOutputArray->Add(vertexcandidate);
     188
    201189    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    202190    {
    … …  
    216204      candidate->Momentum.RotateZ(dphi);
    217205
    218       x -= fInputBeamSpotX;
    219       y -= fInputBeamSpotY;
    220206      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
    221207      candidate->Position.RotateZ(dphi);
    222       candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
    223 
    224       vx += candidate->Position.X();
    225       vy += candidate->Position.Y();
    226       ++numberOfParticles;
    227208
    228209      fParticleOutputArray->Add(candidate);
    229210    }
    230 
    231     if(numberOfParticles > 0)
    232     {
    233       vx /= numberOfParticles;
    234       vy /= numberOfParticles;
    235     }
    236 
    237     vertex = factory->NewCandidate();
    238     vertex->Position.SetXYZT(vx, vy, dz, dt);
    239     vertex->IsPU = 1;
    240 
    241     fVertexOutputArray->Add(vertex);
    242211  }
    243212}
    244213
    245214//------------------------------------------------------------------------------
     215

  • modules/PileUpMerger.h

    r9343566 re33e6db  
    5353  Double_t fTVertexSpread;
    5454
    55   Double_t fInputBeamSpotX;
    56   Double_t fInputBeamSpotY;
    57   Double_t fOutputBeamSpotX;
    58   Double_t fOutputBeamSpotY;
    59 
    6055  DelphesTF2 *fFunction; //!
    6156

  • modules/PileUpMergerPythia8.cc

    r9343566 re33e6db  
    2929#include "classes/DelphesClasses.h"
    3030#include "classes/DelphesFactory.h"
    31 #include "classes/DelphesTF2.h"
     31#include "classes/DelphesFormula.h"
    3232#include "classes/DelphesPileUpReader.h"
    3333
    … …  
    5656
    5757PileUpMergerPythia8::PileUpMergerPythia8() :
    58   fFunction(0), fPythia(0), fItInputArray(0)
     58  fPythia(0), fItInputArray(0)
    5959{
    60   fFunction = new DelphesTF2;
    6160}
    6261
    … …  
    6564PileUpMergerPythia8::~PileUpMergerPythia8()
    6665{
    67   delete fFunction;
    6866}
    6967
    … …  
    7472  const char *fileName;
    7573
    76   fPileUpDistribution = GetInt("PileUpDistribution", 0);
    77 
    7874  fMeanPileUp  = GetDouble("MeanPileUp", 10);
    79 
    80   fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
    81   fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
    82 
    83   fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
    84   fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
    85   fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
    86   fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
     75  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
    8776
    8877  fPTMin = GetDouble("PTMin", 0.0);
    89 
    90   fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
    91   fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
    9278
    9379  fileName = GetString("ConfigFile", "MinBias.cmnd");
    … …  
    10086
    10187  // create output arrays
    102   fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
    103   fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
     88  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
    10489}
    10590
    … …  
    118103  TParticlePDG *pdgParticle;
    119104  Int_t pid, status;
    120   Float_t x, y, z, t, vx, vy;
     105  Float_t x, y, z, t;
    121106  Float_t px, py, pz, e;
    122   Double_t dz, dphi, dt;
    123   Int_t numberOfEvents, event, numberOfParticles, i;
    124   Candidate *candidate, *vertex;
     107  Double_t dz, dphi;
     108  Int_t poisson, event, i;
     109  Candidate *candidate;
    125110  DelphesFactory *factory;
    126111
    127   const Double_t c_light = 2.99792458E8;
    128 
    129112  fItInputArray->Reset();
    130 
    131   // --- Deal with primary vertex first  ------
    132 
    133   fFunction->GetRandom2(dz, dt);
    134 
    135   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    136   dz *= 1.0E3; // necessary in order to make z in mm
    137   vx = 0.0;
    138   vy = 0.0;
    139   numberOfParticles = fInputArray->GetEntriesFast();
    140113  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    141114  {
    142     vx += candidate->Position.X();
    143     vy += candidate->Position.Y();
    144     z = candidate->Position.Z();
    145     t = candidate->Position.T();
    146     candidate->Position.SetZ(z + dz);
    147     candidate->Position.SetT(t + dt);
    148     fParticleOutputArray->Add(candidate);
    149   }
    150 
    151   if(numberOfParticles > 0)
    152   {
    153     vx /= numberOfParticles;
    154     vy /= numberOfParticles;
     115    fOutputArray->Add(candidate);
    155116  }
    156117
    157118  factory = GetFactory();
    158119
    159   vertex = factory->NewCandidate();
    160   vertex->Position.SetXYZT(vx, vy, dz, dt);
    161   fVertexOutputArray->Add(vertex);
     120  poisson = gRandom->Poisson(fMeanPileUp);
    162121
    163   // --- Then with pile-up vertices  ------
    164 
    165   switch(fPileUpDistribution)
    166   {
    167     case 0:
    168       numberOfEvents = gRandom->Poisson(fMeanPileUp);
    169       break;
    170     case 1:
    171       numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
    172       break;
    173     default:
    174       numberOfEvents = gRandom->Poisson(fMeanPileUp);
    175       break;
    176   }
    177 
    178   for(event = 0; event < numberOfEvents; ++event)
     122  for(event = 0; event < poisson; ++event)
    179123  {
    180124    while(!fPythia->next());
    181125
    182    // --- Pile-up vertex smearing
    183 
    184     fFunction->GetRandom2(dz, dt);
    185 
    186     dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    187     dz *= 1.0E3; // necessary in order to make z in mm
    188 
     126    dz = gRandom->Gaus(0.0, fZVertexSpread);
    189127    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    190128
    191     vx = 0.0;
    192     vy = 0.0;
    193     numberOfParticles = fPythia->event.size();
    194     for(i = 1; i < numberOfParticles; ++i)
     129    for(i = 1; i < fPythia->event.size(); ++i)
    195130    {
    196131      Pythia8::Particle &particle = fPythia->event[i];
    … …  
    208143      candidate->PID = pid;
    209144
    210       candidate->Status = 1;
     145      candidate->Status = status;
    211146
    212147      pdgParticle = pdg->GetParticle(pid);
    … …  
    219154      candidate->Momentum.RotateZ(dphi);
    220155
    221       x -= fInputBeamSpotX;
    222       y -= fInputBeamSpotY;
    223       candidate->Position.SetXYZT(x, y, z + dz, t + dt);
     156      candidate->Position.SetXYZT(x, y, z + dz, t);
    224157      candidate->Position.RotateZ(dphi);
    225       candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
    226158
    227       vx += candidate->Position.X();
    228       vy += candidate->Position.Y();
    229 
    230       fParticleOutputArray->Add(candidate);
     159      fOutputArray->Add(candidate);
    231160    }
    232 
    233     if(numberOfParticles > 0)
    234     {
    235       vx /= numberOfParticles;
    236       vy /= numberOfParticles;
    237     }
    238 
    239     vertex = factory->NewCandidate();
    240     vertex->Position.SetXYZT(vx, vy, dz, dt);
    241     vertex->IsPU = 1;
    242 
    243     fVertexOutputArray->Add(vertex);
    244161  }
    245162}
    246163
    247164//------------------------------------------------------------------------------
     165

  • modules/PileUpMergerPythia8.h

    r9343566 re33e6db  
    3131
    3232class TObjArray;
    33 class DelphesTF2;
    3433
    3534namespace Pythia8
    … …  
    5150private:
    5251
    53   Int_t fPileUpDistribution;
    5452  Double_t fMeanPileUp;
    55 
    5653  Double_t fZVertexSpread;
    57   Double_t fTVertexSpread;
    58 
    59   Double_t fInputBeamSpotX;
    60   Double_t fInputBeamSpotY;
    61   Double_t fOutputBeamSpotX;
    62   Double_t fOutputBeamSpotY;
    63 
    6454  Double_t fPTMin;
    65 
    66   DelphesTF2 *fFunction; //!
    6755
    6856  Pythia8::Pythia *fPythia; //!
    … …  
    7260  const TObjArray *fInputArray; //!
    7361
    74   TObjArray *fParticleOutputArray; //!
    75   TObjArray *fVertexOutputArray; //!
     62  TObjArray *fOutputArray; //!
    7663
    7764  ClassDef(PileUpMergerPythia8, 1)

  • modules/SimpleCalorimeter.cc

    r9343566 re33e6db  
    149149
    150150  fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
    151 
    152   // flag that says if current calo is Ecal of Hcal (will then fill correct values of Eem and Ehad)
    153   fIsEcal = GetBool("IsEcal", false);
    154151
    155152  // switch on or off the dithering of the center of calorimeter towers
    … …  
    428425  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    429426
    430   fTower->Eem = (!fIsEcal) ? 0 : energy;
    431   fTower->Ehad = (fIsEcal) ? 0 : energy;
    432 
    433427  fTower->Edges[0] = fTowerEdges[0];
    434428  fTower->Edges[1] = fTowerEdges[1];
    … …  
    453447    pt = energy / TMath::CosH(eta);
    454448
    455     tower->Eem = (!fIsEcal) ? 0 : energy;
    456     tower->Ehad = (fIsEcal) ? 0 : energy;
    457 
    458449    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    459450    fEFlowTowerOutputArray->Add(tower);

  • modules/SimpleCalorimeter.h

    r9343566 re33e6db  
    1 
    21/*
    32 *  Delphes: a framework for fast simulation of a generic collider experiment
    … …  
    7574  Bool_t fSmearTowerCenter;
    7675
    77   Bool_t fIsEcal; //!
    78 
    7976  TFractionMap fFractionMap; //!
    8077  TBinMap fBinMap; //!

  • modules/TrackPileUpSubtractor.cc

    r9343566 re33e6db  
    7474  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    7575
    76   fPTMin = GetDouble("PTMin", 0.);
    7776  // import arrays with output from other modules
    78    
     77
    7978  ExRootConfParam param = GetParam("InputArray");
    8079  Long_t i, size;
    … …  
    148147      // assume perfect pile-up subtraction for tracks outside fZVertexResolution
    149148     
    150       if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1;
    151       else
    152       {
    153          candidate->IsRecoPU = 0;
    154          if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
    155       }
     149      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
     150
     151      array->Add(candidate);
    156152    }
    157153  }

  • modules/TrackPileUpSubtractor.h

    r9343566 re33e6db  
    5050  Double_t fZVertexResolution;
    5151
    52   Double_t fPTMin;
    53 
    5452  std::map< TIterator *, TObjArray * > fInputMap; //!
    5553

  • modules/TreeWriter.cc

    r9343566 re33e6db  
    291291    entry->ZOuter = position.Z();
    292292    entry->TOuter = position.T()*1.0E-3/c_light;
    293 
     293   
    294294    entry->Dxy = candidate->Dxy;
    295295    entry->SDxy = candidate->SDxy ;
    … …  
    297297    entry->Yd = candidate->Yd;
    298298    entry->Zd = candidate->Zd;
    299 
     299   
    300300    const TLorentzVector &momentum = candidate->Momentum;
    301301
    … …  
    362362
    363363    entry->T = position.T()*1.0E-3/c_light;
    364     entry->Ntimes = candidate->Ntimes;
    365364
    366365    FillParticles(candidate, &entry->Particles);
    … …  
    404403    entry->T = position.T()*1.0E-3/c_light;
    405404
    406     // Isolation variables
    407 
    408     entry->IsolationVar = candidate->IsolationVar;
    409     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
    410     entry->SumPtCharged = candidate->SumPtCharged ;
    411     entry->SumPtNeutral = candidate->SumPtNeutral ;
    412     entry->SumPtChargedPU = candidate->SumPtChargedPU ;
    413     entry->SumPt = candidate->SumPt ;
    414 
    415405    entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
    416406
    … …  
    452442    entry->T = position.T()*1.0E-3/c_light;
    453443
    454     // Isolation variables
    455 
    456     entry->IsolationVar = candidate->IsolationVar;
    457     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
    458     entry->SumPtCharged = candidate->SumPtCharged ;
    459     entry->SumPtNeutral = candidate->SumPtNeutral ;
    460     entry->SumPtChargedPU = candidate->SumPtChargedPU ;
    461     entry->SumPt = candidate->SumPt ;
    462 
    463 
    464444    entry->Charge = candidate->Charge;
    465445
    … …  
    489469    const TLorentzVector &momentum = candidate->Momentum;
    490470    const TLorentzVector &position = candidate->Position;
     471
    491472
    492473    pt = momentum.Pt();
    … …  
    507488    entry->T = position.T()*1.0E-3/c_light;
    508489
    509     // Isolation variables
    510 
    511     entry->IsolationVar = candidate->IsolationVar;
    512     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
    513     entry->SumPtCharged = candidate->SumPtCharged ;
    514     entry->SumPtNeutral = candidate->SumPtNeutral ;
    515     entry->SumPtChargedPU = candidate->SumPtChargedPU ;
    516     entry->SumPt = candidate->SumPt ;
    517 
    518490    entry->Charge = candidate->Charge;
    519491
    … …  
    532504  Double_t ecalEnergy, hcalEnergy;
    533505  const Double_t c_light = 2.99792458E8;
    534   Int_t i;
    535506
    536507  array->Sort();
    … …  
    561532    entry->Mass = momentum.M();
    562533
    563     entry->Area = candidate->Area;
    564 
    565534    entry->DeltaEta = candidate->DeltaEta;
    566535    entry->DeltaPhi = candidate->DeltaPhi;
    567536
    568     entry->Flavor = candidate->Flavor;
    569     entry->FlavorAlgo = candidate->FlavorAlgo;
    570     entry->FlavorPhys = candidate->FlavorPhys;
    571 
    572537    entry->BTag = candidate->BTag;
    573 
    574     entry->BTagAlgo = candidate->BTagAlgo;
    575     entry->BTagPhys = candidate->BTagPhys;
    576 
    577538    entry->TauTag = candidate->TauTag;
    578539
    … …  
    600561    entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
    601562    entry->PTD = candidate->PTD;
    602 
    603     //--- Sub-structure variables ----
    604 
    605     entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
    606     entry->NSubJetsPruned = candidate->NSubJetsPruned;
    607     entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
    608 
    609     for(i = 0; i < 5; i++)
    610     {
    611       entry->FracPt[i] = candidate -> FracPt[i];
    612       entry->Tau[i] = candidate -> Tau[i];
    613       entry->TrimmedP4[i] = candidate -> TrimmedP4[i];
    614       entry->PrunedP4[i] = candidate -> PrunedP4[i];
    615       entry->SoftDroppedP4[i] = candidate -> SoftDroppedP4[i];
    616     }
    617 
     563    entry->FracPt[0] = candidate->FracPt[0];
     564    entry->FracPt[1] = candidate->FracPt[1];
     565    entry->FracPt[2] = candidate->FracPt[2];
     566    entry->FracPt[3] = candidate->FracPt[3];
     567    entry->FracPt[4] = candidate->FracPt[4];
     568   
     569    //--- N-subjettiness variables ----
     570   
     571    entry->Tau1 = candidate->Tau[0];
     572    entry->Tau2 = candidate->Tau[1];
     573    entry->Tau3 = candidate->Tau[2];
     574    entry->Tau4 = candidate->Tau[3];
     575    entry->Tau5 = candidate->Tau[4];
     576   
    618577    FillParticles(candidate, &entry->Particles);
    619578  }

  • readers/DelphesLHEF.cpp

    r9343566 re33e6db  
    6363  TStopwatch readStopWatch, procStopWatch;
    6464  ExRootTreeWriter *treeWriter = 0;
    65   ExRootTreeBranch *branchEvent = 0, *branchWeight = 0;
     65  ExRootTreeBranch *branchEvent = 0, *branchRwgt = 0;
    6666  ExRootConfReader *confReader = 0;
    6767  Delphes *modularDelphes = 0;
    … …  
    103103
    104104    branchEvent = treeWriter->NewBranch("Event", LHEFEvent::Class());
    105     branchWeight = treeWriter->NewBranch("Weight", Weight::Class());
     105    branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
    106106
    107107    confReader = new ExRootConfReader;
    … …  
    196196
    197197            reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch);
    198             reader->AnalyzeWeight(branchWeight);
     198            reader->AnalyzeRwgt(branchRwgt);
    199199
    200200            treeWriter->Fill();

  • readers/DelphesPythia8.cpp

    r9343566 re33e6db  
    2020#include <iostream>
    2121#include <sstream>
    22 #include <string>
    2322
    2423#include <signal.h>
    … …  
    3938#include "classes/DelphesClasses.h"
    4039#include "classes/DelphesFactory.h"
    41 #include "classes/DelphesLHEFReader.h"
    4240
    4341#include "ExRootAnalysis/ExRootTreeWriter.h"
    … …  
    151149  char appName[] = "DelphesPythia8";
    152150  stringstream message;
    153   FILE *inputFile = 0;
    154151  TFile *outputFile = 0;
    155152  TStopwatch readStopWatch, procStopWatch;
    156153  ExRootTreeWriter *treeWriter = 0;
    157154  ExRootTreeBranch *branchEvent = 0;
    158   ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0;
    159155  ExRootConfReader *confReader = 0;
    160156  Delphes *modularDelphes = 0;
    161157  DelphesFactory *factory = 0;
    162158  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
    163   TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0;
    164   DelphesLHEFReader *reader = 0;
    165159  Long64_t eventCounter, errorCounter;
    166160  Long64_t numberOfEvents, timesAllowErrors;
    … …  
    211205    partonOutputArray = modularDelphes->ExportArray("partons");
    212206
    213     // Initialize Pythia
     207    modularDelphes->InitTask();
     208
     209    // Initialize pythia
    214210    pythia = new Pythia8::Pythia;
    215211
    … …  
    225221    numberOfEvents = pythia->mode("Main:numberOfEvents");
    226222    timesAllowErrors = pythia->mode("Main:timesAllowErrors");
    227 
    228     inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
    229     if(inputFile)
    230     {
    231       reader = new DelphesLHEFReader;
    232       reader->SetInputFile(inputFile);
    233 
    234       branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
    235       branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
    236 
    237       allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
    238       stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
    239       partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
    240     }
    241 
    242     modularDelphes->InitTask();
    243223
    244224    pythia->init();
    … …  
    254234    for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
    255235    {
    256       while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
    257         stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
    258 
    259236      if(!pythia->next())
    260237      {
    261238        // If failure because reached end of file then exit event loop
    262         if(pythia->info.atEndOfFile())
     239        if (pythia->info.atEndOfFile())
    263240        {
    264241          cerr << "Aborted since reached end of Les Houches Event File" << endl;
    … …  
    267244
    268245        // First few failures write off as "acceptable" errors, then quit
    269         if(++errorCounter > timesAllowErrors)
    270         {
    271           cerr << "Event generation aborted prematurely, owing to error!" << endl;
    272           break;
    273         }
    274 
    275         modularDelphes->Clear();
    276         reader->Clear();
    277         continue;
     246        if (++errorCounter < timesAllowErrors) continue;
     247        cerr << "Event generation aborted prematurely, owing to error!" << endl;
     248        break;
    278249      }
    279250
    … …  
    287258      procStopWatch.Stop();
    288259
    289       if(reader)
    290       {
    291         reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch);
    292         reader->AnalyzeWeight(branchWeightLHEF);
    293       }
    294 
    295260      treeWriter->Fill();
    296261
    297262      treeWriter->Clear();
    298263      modularDelphes->Clear();
    299       if(reader) reader->Clear();
    300264
    301265      readStopWatch.Start();
    … …  
    313277    cout << "** Exiting..." << endl;
    314278
    315     delete reader;
    316279    delete pythia;
    317280    delete modularDelphes;
Note: See TracChangeset for help on using the changeset viewer.