Fork me on GitHub

Changes in / [e33e6db:9343566] in git


Ignore:
Files:
27 added
47 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    re33e6db r9343566  
    257257        classes/DelphesClasses.h \
    258258        classes/DelphesFactory.h \
     259        classes/DelphesLHEFReader.h \
    259260        external/ExRootAnalysis/ExRootTreeWriter.h \
    260261        external/ExRootAnalysis/ExRootTreeBranch.h \
    … …  
    337338        modules/Weighter.h \
    338339        modules/Hector.h \
     340        modules/RunPUPPI.h \
     341        modules/JetFlavorAssociation.h \
    339342        modules/ExampleModule.h
    340343ModulesDict$(PcmSuf): \
    … …  
    521524tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
    522525        external/Hector/H_VerticalQuadrupole.$(SrcSuf)
     526tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
     527        external/PUPPI/puppiCleanContainer.$(SrcSuf) \
     528        external/fastjet/Selector.hh
    523529tmp/modules/AngularSmearing.$(ObjSuf): \
    524530        modules/AngularSmearing.$(SrcSuf) \
    … …  
    652658        external/ExRootAnalysis/ExRootFilter.h \
    653659        external/ExRootAnalysis/ExRootClassifier.h
     660tmp/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
    654669tmp/modules/JetPileUpSubtractor.$(ObjSuf): \
    655670        modules/JetPileUpSubtractor.$(SrcSuf) \
    … …  
    730745        classes/DelphesClasses.h \
    731746        classes/DelphesFactory.h \
    732         classes/DelphesFormula.h \
     747        classes/DelphesTF2.h \
    733748        classes/DelphesPileUpReader.h \
    734749        external/ExRootAnalysis/ExRootResult.h \
    735750        external/ExRootAnalysis/ExRootFilter.h \
    736751        external/ExRootAnalysis/ExRootClassifier.h
     752tmp/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
    737762tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    738763        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    869894        tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
    870895        tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
     896        tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
    871897        tmp/modules/AngularSmearing.$(ObjSuf) \
    872898        tmp/modules/BTagging.$(ObjSuf) \
    … …  
    883909        tmp/modules/ImpactParameterSmearing.$(ObjSuf) \
    884910        tmp/modules/Isolation.$(ObjSuf) \
     911        tmp/modules/JetFlavorAssociation.$(ObjSuf) \
    885912        tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
    886913        tmp/modules/LeptonDressing.$(ObjSuf) \
    … …  
    891918        tmp/modules/PileUpJetID.$(ObjSuf) \
    892919        tmp/modules/PileUpMerger.$(ObjSuf) \
     920        tmp/modules/RunPUPPI.$(ObjSuf) \
    893921        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    894922        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    10801108tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf): \
    10811109        external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(SrcSuf)
     1110tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf): \
     1111        external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(SrcSuf) \
     1112        external/fastjet/JetDefinition.hh \
     1113        external/fastjet/ClusterSequenceAreaBase.hh
     1114tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf): \
     1115        external/fastjet/contribs/RecursiveTools/Recluster.$(SrcSuf)
     1116tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf): \
     1117        external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(SrcSuf) \
     1118        external/fastjet/JetDefinition.hh \
     1119        external/fastjet/ClusterSequenceAreaBase.hh
     1120tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf): \
     1121        external/fastjet/contribs/RecursiveTools/SoftDrop.$(SrcSuf)
    10821122tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf): \
    10831123        external/fastjet/contribs/SoftKiller/SoftKiller.$(SrcSuf)
    … …  
    12131253        external/fastjet/contribs/Nsubjettiness/Njettiness.hh \
    12141254        external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh \
    1215         external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.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
    12161259tmp/modules/FastJetGridMedianEstimator.$(ObjSuf): \
    12171260        modules/FastJetGridMedianEstimator.$(SrcSuf) \
    … …  
    12851328        tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf) \
    12861329        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) \
    12871334        tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf) \
    12881335        tmp/external/fastjet/plugins/ATLASCone/ATLASConePlugin.$(ObjSuf) \
    … …  
    13381385        display/Delphes3DGeometry.$(SrcSuf) \
    13391386        display/Delphes3DGeometry.h \
    1340         external/ExRootAnalysis/ExRootConfReader.h \
    1341         classes/DelphesClasses.h
     1387        classes/DelphesClasses.h \
     1388        external/ExRootAnalysis/ExRootConfReader.h
    13421389tmp/display/DelphesBranchElement.$(ObjSuf): \
    13431390        display/DelphesBranchElement.$(SrcSuf) \
    … …  
    13521399tmp/display/DelphesEventDisplay.$(ObjSuf): \
    13531400        display/DelphesEventDisplay.$(SrcSuf) \
    1354         external/ExRootAnalysis/ExRootConfReader.h \
    1355         external/ExRootAnalysis/ExRootTreeReader.h \
    13561401        display/DelphesCaloData.h \
    13571402        display/DelphesBranchElement.h \
    13581403        display/Delphes3DGeometry.h \
    13591404        display/DelphesEventDisplay.h \
    1360         classes/DelphesClasses.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
    13611412tmp/display/DelphesHtmlSummary.$(ObjSuf): \
    13621413        display/DelphesHtmlSummary.$(SrcSuf) \
    … …  
    15411592        @touch $@
    15421593
    1543 external/fastjet/Selector.hh: \
    1544         external/fastjet/PseudoJet.hh \
    1545         external/fastjet/RangeDefinition.hh
    1546         @touch $@
    1547 
    15481594external/fastjet/internal/Dnn2piCylinder.hh: \
    15491595        external/fastjet/internal/DynamicNearestNeighbours.hh \
    15501596        external/fastjet/internal/DnnPlane.hh \
    15511597        external/fastjet/internal/numconsts.hh
     1598        @touch $@
     1599
     1600external/fastjet/Selector.hh: \
     1601        external/fastjet/PseudoJet.hh \
     1602        external/fastjet/RangeDefinition.hh
    15521603        @touch $@
    15531604
    … …  
    16391690        @touch $@
    16401691
     1692modules/RunPUPPI.h: \
     1693        classes/DelphesModule.h
     1694        @touch $@
     1695
    16411696modules/Cloner.h: \
    16421697        classes/DelphesModule.h
    … …  
    16751730        @touch $@
    16761731
    1677 display/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 
    16851732modules/TauTagging.h: \
    16861733        classes/DelphesModule.h \
    … …  
    17251772        @touch $@
    17261773
     1774modules/JetFlavorAssociation.h: \
     1775        classes/DelphesModule.h \
     1776        classes/DelphesClasses.h
     1777        @touch $@
     1778
    17271779modules/ParticlePropagator.h: \
    17281780        classes/DelphesModule.h
    … …  
    17571809        external/fastjet/AreaDefinition.hh \
    17581810        external/fastjet/ClusterSequenceAreaBase.hh
     1811        @touch $@
     1812
     1813external/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
    17591819        @touch $@
    17601820

  • cards/delphes_card_ATLAS.tcl

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

  • cards/delphes_card_ATLAS_PileUp.tcl

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

  • cards/delphes_card_FCC_basic.tcl

    re33e6db r9343566  
    226226  set EFlowTowerOutputArray eflowPhotons
    227227
     228  set IsEcal true
     229 
    228230  set EnergyMin 0.5
    229231  set EnergySignificanceMin 1.0
    … …  
    290292  set EFlowTowerOutputArray eflowNeutralHadrons
    291293
     294  set IsEcal false
     295 
    292296  set EnergyMin 1.0
    293297  set EnergySignificanceMin 1.0
    … …  
    431435
    432436module FastJetFinder FastJetFinder {
    433 #  set InputArray Calorimeter/towers
     437#  set InputArray TowerMerger/towers
    434438  set InputArray EFlowMerger/eflow
    435439

  • cards/delphes_card_LHCb.tcl

    re33e6db r9343566  
    223223  set TowerOutputArray ecalTowers
    224224  set EFlowTowerOutputArray eflowPhotons
     225
     226  set IsEcal true
    225227
    226228  set EnergyMin 0.0
    … …  
    303305  set TowerOutputArray hcalTowers
    304306  set EFlowTowerOutputArray eflowNeutralHadrons
     307
     308  set IsEcal false
    305309
    306310  set EnergyMin 0.0

  • classes/ClassesLinkDef.h

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

  • classes/DelphesClasses.cc

    re33e6db r9343566  
    120120  PID(0), Status(0), M1(-1), M2(-1), D1(-1), D2(-1),
    121121  Charge(0), Mass(0.0),
    122   IsPU(0), IsConstituent(0),
    123   BTag(0), TauTag(0), Eem(0.0), Ehad(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),
    124126  DeltaEta(0.0), DeltaPhi(0.0),
    125127  Momentum(0.0, 0.0, 0.0, 0.0),
    … …  
    129131  NCharged(0),
    130132  NNeutrals(0),
     133  NSubJetsTrimmed(0),
     134  NSubJetsPruned(0),
     135  NSubJetsSoftDropped(0),
    131136  Beta(0),
    132137  BetaStar(0),
    133138  MeanSqDeltaR(0),
    134139  PTD(0),
     140  Ntimes(-1),
     141  IsolationVar(-999),
     142  IsolationVarRhoCorr(-999),
     143  SumPtCharged(-999),
     144  SumPtNeutral(-999),
     145  SumPtChargedPU(-999),
     146  SumPt(-999),
    135147  fFactory(0),
    136148  fArray(0)
    137149{
     150  int i;
    138151  Edges[0] = 0.0;
    139152  Edges[1] = 0.0;
    … …  
    150163  Tau[3] = 0.0;
    151164  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  }
    152171}
    153172
    … …  
    224243  object.IsPU = IsPU;
    225244  object.IsConstituent = IsConstituent;
     245  object.Flavor = Flavor;
     246  object.FlavorAlgo = FlavorAlgo;
     247  object.FlavorPhys = FlavorPhys;
    226248  object.BTag = BTag;
     249  object.BTagAlgo = BTagAlgo;
     250  object.BTagPhys = BTagPhys;
    227251  object.TauTag = TauTag;
    228252  object.Eem = Eem;
    … …  
    249273  object.MeanSqDeltaR = MeanSqDeltaR;
    250274  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
    251283  object.FracPt[0] = FracPt[0];
    252284  object.FracPt[1] = FracPt[1];
    … …  
    259291  object.Tau[3] = Tau[3];
    260292  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;
    261313
    262314  object.fFactory = fFactory;
    263315  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));
    264319
    265320  if(fArray && fArray->GetEntriesFast() > 0)
    … …  
    278333void Candidate::Clear(Option_t* option)
    279334{
     335  int i;
    280336  SetUniqueID(0);
    281337  ResetBit(kIsReferenced);
    … …  
    287343  IsPU = 0;
    288344  IsConstituent = 0;
     345  Flavor = 0;
     346  FlavorAlgo = 0;
     347  FlavorPhys = 0;
    289348  BTag = 0;
     349  BTagAlgo = 0;
     350  BTagPhys = 0;
    290351  TauTag = 0;
    291352  Eem = 0.0;
    … …  
    311372  MeanSqDeltaR = 0.0;
    312373  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 
    313385  FracPt[0] = 0.0;
    314386  FracPt[1] = 0.0;
    … …  
    321393  Tau[3] = 0.0;
    322394  Tau[4] = 0.0;
    323 
     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 
    324407  fArray = 0;
    325408}

  • classes/DelphesClasses.h

    re33e6db r9343566  
    8484//---------------------------------------------------------------------------
    8585
     86class LHEFWeight: public TObject
     87{
     88public:
     89  Int_t ID; // weight ID
     90  Float_t Weight; // weight value
     91
     92  ClassDef(LHEFWeight, 1)
     93};
     94
     95//---------------------------------------------------------------------------
     96
    8697class HepMCEvent: public Event
    8798{
    … …  
    231242  TRefArray Particles; // references to generated particles
    232243
     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
    233253  static CompBase *fgCompare; //!
    234254  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    257277  TRef Particle; // reference to generated particle
    258278
     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
    259288  static CompBase *fgCompare; //!
    260289  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    281310  TRef Particle; // reference to generated particle
    282311
     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
    283321  static CompBase *fgCompare; //!
    284322  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    306344  Float_t DeltaPhi;  // jet radius in azimuthal angle
    307345
     346  UInt_t Flavor;
     347  UInt_t FlavorAlgo;
     348  UInt_t FlavorPhys;
     349
    308350  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
    309354  UInt_t TauTag; // 0 or 1 for a jet that has been tagged as a tau
    310355
    … …  
    313358  Float_t EhadOverEem; // ratio of the hadronic versus electromagnetic energy deposited in the calorimeter
    314359
    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
     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
    328377
    329378  TRefArray Constituents; // references to constituents
    … …  
    334383
    335384  TLorentzVector P4() const;
    336 
    337   ClassDef(Jet, 2)
     385  TLorentzVector Area;
     386
     387  ClassDef(Jet, 3)
    338388};
    339389
    … …  
    392442  Float_t E; // calorimeter tower energy
    393443
    394   Float_t T; //particle arrival time of flight
     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
    395446
    396447  Float_t Eem; // calorimeter tower electromagnetic energy
    … …  
    452503
    453504  Int_t IsPU;
     505  Int_t IsRecoPU;
     506
    454507  Int_t IsConstituent;
    455508
     509  UInt_t Flavor;
     510  UInt_t FlavorAlgo;
     511  UInt_t FlavorPhys;
     512
    456513  UInt_t BTag;
     514  UInt_t BTagAlgo;
     515  UInt_t BTagPhys;
     516
    457517  UInt_t TauTag;
    458518
    … …  
    482542  Float_t  FracPt[5];
    483543
     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
    484558  // N-subjettiness variables
    485559
    486560  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
    487572
    488573  static CompBase *fgCompare; //!
    … …  
    504589  void SetFactory(DelphesFactory *factory) { fFactory = factory; }
    505590
    506   ClassDef(Candidate, 2)
     591  ClassDef(Candidate, 3)
    507592};
    508593

  • classes/DelphesFormula.cc

    re33e6db r9343566  
    2323
    2424#include <stdexcept>
    25 #include <string>
    2625
    2726using namespace std;
    … …  
    5150Int_t DelphesFormula::Compile(const char *expression)
    5251{
    53   string buffer;
     52  TString buffer;
    5453  const char *it;
    5554  for(it = expression; *it; ++it)
    5655  {
    5756    if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    58     buffer.push_back(*it);
     57    buffer.Append(*it);
    5958  }
    60   if(TFormula::Compile(buffer.c_str()) != 0)
     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)
    6164  {
    6265    throw runtime_error("Invalid formula.");
    … …  
    7477
    7578//------------------------------------------------------------------------------
    76 
    77 Int_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

    re33e6db r9343566  
    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);
    3937};
    4038

  • classes/DelphesLHEFReader.cc

    re33e6db r9343566  
    8282  fEventCounter = -1;
    8383  fParticleCounter = -1;
    84   fRwgtList.clear();
     84  fWeightList.clear();
    8585}
    8686
    … …  
    9999  TObjArray *partonOutputArray)
    100100{
    101   int rc;
     101  int rc, id;
    102102  char *pch;
    103103  double weight;
    … …  
    158158  else if(strstr(fBuffer, "<wgt"))
    159159  {
    160     pch = strstr(fBuffer, ">");
     160    pch = strpbrk(fBuffer, "\"'");
    161161    if(!pch)
    162162    {
    … …  
    165165    }
    166166
    167     DelphesStream bufferStream(pch + 1);
    168     rc = bufferStream.ReadDbl(weight);
     167    DelphesStream idStream(pch + 1);
     168    rc = idStream.ReadInt(id);
     169
     170    pch = strchr(fBuffer, '>');
     171    if(!pch)
     172    {
     173      cerr << "** ERROR: " << "invalid weight format" << endl;
     174      return kFALSE;
     175    }
     176
     177    DelphesStream weightStream(pch + 1);
     178    rc = weightStream.ReadDbl(weight);
    169179
    170180    if(!rc)
    … …  
    174184    }
    175185
    176     fRwgtList.push_back(weight);
     186    fWeightList.push_back(make_pair(id, weight));
    177187  }
    178188  else if(strstr(fBuffer, "</event>"))
    … …  
    206216//---------------------------------------------------------------------------
    207217
    208 void 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;
     218void 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;
    218229  }
    219230}

  • classes/DelphesLHEFReader.h

    re33e6db r9343566  
    3131
    3232#include <vector>
     33#include <utility>
    3334
    3435class TObjArray;
    … …  
    5859    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
    5960
    60   void AnalyzeRwgt(ExRootTreeBranch *branch);
     61  void AnalyzeWeight(ExRootTreeBranch *branch);
    6162
    6263private:
    … …  
    8384  double fPx, fPy, fPz, fE, fMass;
    8485 
    85   std::vector<double> fRwgtList;
     86  std::vector< std::pair< int, double > > fWeightList;
    8687};
    8788

  • classes/DelphesTF2.cc

    re33e6db r9343566  
    1818
    1919#include "classes/DelphesTF2.h"
     20
     21#include "RVersion.h"
    2022#include "TString.h"
     23
    2124#include <stdexcept>
    22 #include <string>
    2325
    2426using namespace std;
    … …  
    3436
    3537DelphesTF2::DelphesTF2(const char *name, const char *expression) :
    36   TF2(name,expression)
     38  TF2(name, expression)
    3739{
    3840}
    … …  
    4648//------------------------------------------------------------------------------
    4749
    48 Int_t DelphesTF2::DefinedVariable(TString &chaine, Int_t &action)
     50Int_t DelphesTF2::Compile(const char *expression)
    4951{
    50   action = kVariable;
    51   if(chaine == "z")
     52  TString buffer;
     53  const char *it;
     54  for(it = expression; *it; ++it)
    5255  {
    53     if(fNdim < 1) fNdim = 1;
    54     return 0;
     56    if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
     57    buffer.Append(*it);
    5558  }
    56   else if(chaine == "t")
     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
    5766  {
    58     if(fNdim < 2) fNdim = 2;
    59     return 1;
     67    throw runtime_error("Invalid formula.");
    6068  }
    61   return -1;
     69  return 0;
    6270}
    6371

  • classes/DelphesTF2.h

    re33e6db r9343566  
    2121
    2222#include "TF2.h"
    23 #include "TFormula.h"
    24 
    25 #include <string>
    2623
    2724class DelphesTF2: public TF2
    … …  
    3532  ~DelphesTF2();
    3633
    37   Int_t DefinedVariable(TString &variable, Int_t &action);
    38 
     34  Int_t Compile(const char *expression);
    3935};
    4036
    4137#endif /* DelphesTF2_h */
    42 

  • display/Delphes3DGeometry.cc

    re33e6db r9343566  
    1717 */
    1818
    19 #include "display/Delphes3DGeometry.h"
    2019#include <set>
    2120#include <map>
    22 #include <string>
    2321#include <utility>
    2422#include <vector>
    … …  
    2624#include <sstream>
    2725#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"
    3937#include "TF2.h"
    4038#include "TFormula.h"
    4139#include "TH1F.h"
    4240#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"
    4347
    4448using namespace std;
    … …  
    9397   tk_Bz_     = confReader->GetDouble("ParticlePropagator::Bz", 0.0);                           // tk_Bz
    9498   
    95    string buffer;
     99   TString buffer;
    96100   const char *it;
    97101 
    … …  
    103107   tkEffFormula.ReplaceAll("phi","0.");
    104108 
     109   buffer.Clear();
    105110   for(it = tkEffFormula.Data(); *it; ++it)
    106111   {
    107112     if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    108      buffer.push_back(*it);
    109    }
    110 
    111    TF2* tkEffFunction = new TF2("tkEff",buffer.c_str(),0,1000,-10,10);
     113     buffer.Append(*it);
     114   }
     115
     116   TF2* tkEffFunction = new TF2("tkEff",buffer,0,1000,-10,10);
    112117   TH1F etaHisto("eta","eta",100,5.,-5.);
    113118   Double_t pt,eta;
    … …  
    132137   muonEffFormula.ReplaceAll("phi","0.");
    133138   
    134    buffer.clear();
     139   buffer.Clear();
    135140   for(it = muonEffFormula.Data(); *it; ++it)
    136141   {
    137142     if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
    138      buffer.push_back(*it);
    139    }
    140 
    141    TF2* muEffFunction = new TF2("muEff",buffer.c_str(),0,1000,-10,10);
     143     buffer.Append(*it);
     144   }
     145
     146   TF2* muEffFunction = new TF2("muEff",buffer,0,1000,-10,10);
    142147   TH1F etaHisto("eta2","eta2",100,5.,-5.);
    143148   Double_t pt,eta;

  • display/Delphes3DGeometry.h

    re33e6db r9343566  
    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>
    2524#include <vector>
    26 #include <algorithm>
    27 #include <sstream>
    28 #include "TAxis.h"
    29 #include "TGeoManager.h"
    30 #include "TGeoVolume.h"
    31 #include "TGeoMedium.h"
     25
     26#include "Rtypes.h"
     27
     28class TAxis;
     29class TGeoManager;
     30class TGeoVolume;
     31class TGeoMedium;
    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

    re33e6db r9343566  
    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
    2324#include "TGeoManager.h"
    2425#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"
    3226#include "TEveElement.h"
    3327#include "TEveJetCone.h"
    … …  
    5347#include "TCanvas.h"
    5448#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"
    5566
    5667DelphesEventDisplay::DelphesEventDisplay()
    … …  
    98109   TEveManager::Create(kTRUE, "IV");
    99110   fStatusBar_ = gEve->GetBrowser()->GetStatusBar();
    100    TGeoManager* geom = gGeoManager;
     111   TGeoManager *geom = gGeoManager;
    101112
    102113   // build the detector
    … …  
    108119   etaAxis_ = det3D.getCaloAxes().first;
    109120   phiAxis_ = det3D.getCaloAxes().second;
    110    TGeoVolume* top = det3D.getDetector(false);
     121   TGeoVolume *top = det3D.getDetector(false);
    111122   geom->SetTopVolume(top);
    112123   TEveElementList *geometry = new TEveElementList("Geometry");
    113    TObjArray* nodes = top->GetNodes();
     124   TObjArray *nodes = top->GetNodes();
    114125   TIter itNodes(nodes);
    115    TGeoNode* nodeobj;
    116    TEveGeoTopNode* node;
     126   TGeoNode *nodeobj;
     127   TEveGeoTopNode *node;
    117128   while((nodeobj = (TGeoNode*)itNodes.Next())) {
    118129     node = new TEveGeoTopNode(gGeoManager,nodeobj);
    … …  
    131142   // prepare data collections
    132143   readConfig(configFile, elements_);
    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);
     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);
    136147     if(item_v1) gEve->AddElement(item_v1->GetContainer());
    137148     if(item_v2) gEve->AddElement(item_v2->GetContainer());
    … …  
    144155   delphesDisplay_->ImportGeomRhoZ(geometry);
    145156   // find the first calo data and use that to initialize the calo display
    146    for(std::vector<DelphesBranchBase*>::iterator data=elements_.begin();data<elements_.end();++data) {
     157   for(std::vector<DelphesBranchBase *>::iterator data=elements_.begin();data<elements_.end();++data) {
    147158     if(TString((*data)->GetType())=="Tower") { // we could also use GetClassName()=="DelphesCaloData"
    148        DelphesCaloData* container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>((*data))->GetContainer();
     159       DelphesCaloData *container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>((*data))->GetContainer();
    149160       assert(container);
    150161       TEveCalo3D *calo3d = new TEveCalo3D(container);
    … …  
    182193   ExRootConfParam branches = confReader->GetParam("TreeWriter::Branch");
    183194   Int_t nBranches = branches.GetSize()/3;
    184    DelphesBranchElement<TEveTrackList>* tlist;
    185    DelphesBranchElement<DelphesCaloData>* clist;
    186    DelphesBranchElement<TEveElementList>* elist;
     195   DelphesBranchElement<TEveTrackList> *tlist;
     196   DelphesBranchElement<DelphesCaloData> *clist;
     197   DelphesBranchElement<TEveElementList> *elist;
    187198   // first loop with all but tracks
    188199   for(Int_t b = 0; b<nBranches; ++b) {
    … …  
    273284
    274285   // update display
    275    TEveElement* top = (TEveElement*)gEve->GetCurrentEvent();
     286   TEveElement *top = (TEveElement*)gEve->GetCurrentEvent();
    276287   delphesDisplay_->DestroyEventRPhi();
    277288   delphesDisplay_->ImportEventRPhi(top);
    … …  
    356367
    357368   // add a tab on the left
    358    TEveBrowser* browser = gEve->GetBrowser();
     369   TEveBrowser *browser = gEve->GetBrowser();
    359370   browser->SetWindowName("Delphes Event Display");
    360371   browser->StartEmbedding(TRootBrowser::kLeft);
    361372
    362373   // set the main title
    363    TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
     374   TGMainFrame *frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
    364375   frmMain->SetWindowName("Delphes Event Display");
    365376   frmMain->SetCleanup(kDeepCleanup);
    … …  
    371382   if(!gSystem->OpenDirectory(icondir))
    372383     icondir = Form("%s/icons/", (const char*)gSystem->GetFromPipe("root-config --etcdir") );
    373    TGGroupFrame* vf = new TGGroupFrame(frmMain,"Event navigation",kVerticalFrame | kFitWidth );
     384   TGGroupFrame *vf = new TGGroupFrame(frmMain,"Event navigation",kVerticalFrame | kFitWidth );
    374385   {
    375      TGHorizontalFrame* hf = new TGHorizontalFrame(frmMain);
     386     TGHorizontalFrame *hf = new TGHorizontalFrame(frmMain);
    376387     {
    377         TGPictureButton* b = 0;
     388        TGPictureButton *b = 0;
    378389
    379390        b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif"));
    … …  
    381392        b->Connect("Clicked()", "DelphesEventDisplay", this, "Bck()");
    382393
    383         TGNumberEntry* numberEntry = new TGNumberEntry(hf,0,9,-1,TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, 0, treeReader_->GetEntries());
     394        TGNumberEntry *numberEntry = new TGNumberEntry(hf,0,9,-1,TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, 0, treeReader_->GetEntries());
    384395        hf->AddFrame(numberEntry, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY , 2, 0, 10, 10));
    385396        this->Connect("EventChanged(Int_t)","TGNumberEntry",numberEntry,"SetIntNumber(Long_t)");
    … …  
    394405     vf->AddFrame(hf, new TGLayoutHints(kLHintsExpandX , 2, 2, 2, 2));
    395406
    396      TGHProgressBar* progress = new TGHProgressBar(frmMain, TGProgressBar::kFancy, 100);
     407     TGHProgressBar *progress = new TGHProgressBar(frmMain, TGProgressBar::kFancy, 100);
    397408     progress->SetMax( treeReader_->GetEntries());
    398409     progress->ShowPosition(kTRUE, kFALSE, "Event %.0f");
    … …  
    418429   // the summary tab
    419430   htmlSummary_ = new DelphesHtmlSummary("Delphes Event Display Summary Table");
    420    TEveWindowSlot* slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
     431   TEveWindowSlot *slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
    421432   gHtml_ = new TGHtml(0, 100, 100);
    422433   TEveWindowFrame *wf = slot->MakeFrame(gHtml_);
    … …  
    426437   // plot tab
    427438   slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());
    428    TEveWindowTab* tab = slot->MakeTab();
     439   TEveWindowTab *tab = slot->MakeTab();
    429440   tab->SetElementName("Summary plots");
    430441   tab->SetShowTitleBar(kFALSE);
    … …  
    435446}
    436447
    437 
     448void 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
     456void DelphesEventDisplay::Bck() {
     457  if (event_id_ > 0) {
     458     EventChanged(event_id_-1);
     459  } else {
     460     printf("Already at first event.\n");
     461  }
     462}
     463
     464void DelphesEventDisplay::PreSetEv(char *ev) {
     465  event_id_tmp_ = Int_t(atoi(ev));
     466}
     467
     468void 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
     476void DelphesEventDisplay::InitSummaryPlots() {
     477  plotSummary_->FillSample(treeReader_, event_id_);
     478  plotSummary_->FillEvent();
     479  plotSummary_->Draw();
     480}
     481
     482void DelphesEventDisplay::DisplayProgress(Int_t p) {
     483  fStatusBar_->SetText(Form("Processing... %d %%",p), 1);
     484  gSystem->ProcessEvents();
     485}

  • display/DelphesEventDisplay.h

    re33e6db r9343566  
    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>
    3823
     24#include "Rtypes.h"
     25#include "RQ_OBJECT.h"
    3926
     27class TAxis;
     28class TChain;
     29class TGHtml;
     30class TGStatusBar;
     31class DelphesDisplay;
     32class Delphes3DGeometry;
     33class DelphesBranchBase;
     34class DelphesHtmlSummary;
     35class DelphesPlotSummary;
     36class ExRootTreeReader;
    4037
    4138/*
    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  *
     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  *
    4643 */
    4744
    … …  
    5956    void make_gui();
    6057    void load_event();
    61     void readConfig(const char *configFile, std::vector<DelphesBranchBase*>& elements);
     58    void readConfig(const char *configFile, std::vector<DelphesBranchBase *>& elements);
    6259
    6360    // Configuration and global variables.
    … …  
    6764    Double_t tkRadius_, totRadius_, tkHalfLength_, muHalfLength_, bz_;
    6865    TAxis *etaAxis_, *phiAxis_;
    69     TChain* chain_;
    70     std::vector<DelphesBranchBase*> elements_;
     66    TChain *chain_;
     67    std::vector<DelphesBranchBase *> elements_;
    7168    DelphesDisplay *delphesDisplay_;
    7269    DelphesHtmlSummary *htmlSummary_;
    7370    TGHtml *gHtml_;
    7471    DelphesPlotSummary *plotSummary_;
    75     TGStatusBar* fStatusBar_;
     72    TGStatusBar *fStatusBar_;
    7673   
    7774    // gui controls
    7875  public:
    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      }
     76     void Fwd();
    8677
    87      void Bck() {
    88         if (event_id_ > 0) {
    89            EventChanged(event_id_-1);
    90         } else {
    91            printf("Already at first event.\n");
    92         }
    93      }
     78     void Bck();
    9479
    95     void PreSetEv(char* ev) {
    96       event_id_tmp_ = Int_t(atoi(ev));
    97     }
     80    void PreSetEv(char *ev);
    9881
    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     }
     82    void GoTo();
    10683
    107     void InitSummaryPlots() {
    108       plotSummary_->FillSample(treeReader_, event_id_);
    109       plotSummary_->FillEvent();
    110       plotSummary_->Draw();
    111     }
     84    void InitSummaryPlots();
    11285
    113     void DisplayProgress(Int_t p) {
    114          fStatusBar_->SetText(Form("Processing... %d %%",p), 1);
    115          gSystem->ProcessEvents();
    116     }
     86    void DisplayProgress(Int_t p);
    11787};
    11888

  • doc/genMakefile.tcl

    re33e6db r9343566  
    283283dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h}
    284284
    285 sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc}
     285sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} {external/PUPPI/*.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

    re33e6db r9343566  
    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__
     7R__LOAD_LIBRARY(libEve)
     8R__LOAD_LIBRARY(libDelphesDisplay)
     9#include "display/DelphesEventDisplay.h"
     10#include "display/Delphes3DGeometry.h"
     11#endif
    512
    613void EventDisplay(const char *configfile = "delphes_card_CMS.tcl",
    … …  
    3340
    3441    // create the application
    35     DelphesEventDisplay* display = new DelphesEventDisplay(configfile, datafile, det3D);
     42    DelphesEventDisplay *display = new DelphesEventDisplay(configfile, datafile, det3D);
    3643  }
    3744}

  • examples/Example1.C

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

  • examples/Example2.C

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

  • examples/Example3.C

    re33e6db r9343566  
    66*/
    77
     8#ifdef __CLING__
     9R__LOAD_LIBRARY(libDelphes)
     10#include "classes/DelphesClasses.h"
     11#include "external/ExRootAnalysis/ExRootTreeReader.h"
     12#include "external/ExRootAnalysis/ExRootResult.h"
     13#else
     14class ExRootTreeReader;
     15class ExRootResult;
     16#endif
     17
    818//------------------------------------------------------------------------------
    919
    … …  
    2535  TH1 *fJetDeltaPT;
    2636};
    27 
    28 //------------------------------------------------------------------------------
    29 
    30 class ExRootResult;
    31 class ExRootTreeReader;
    3237
    3338//------------------------------------------------------------------------------

  • examples/Example4.C

    re33e6db r9343566  
    11/*
    2 
    32This macro shows how to compute jet energy scale.
    43root -l examples/Example4.C'("delphes_output.root", "plots.root")'
    54
    6 The output ROOT file contains the pT(MC)/pT(Reco) distributions for various pT(Reco) and |eta| bins.
    7 The peak value of such distribution is interpreted as the jet energy correction to be applied for that given pT(Reco), |eta| bin.
    8 
    9 This can be done by modifying the "ScaleFormula" input parameter to the JetEnergyScale module in the delphes_card_XXX.tcl
    10 
    11 
     5The output ROOT file contains the pT(MC)/pT(Reco) distributions for
     6various pT(Reco) and |eta| bins. The peak value of such distribution is
     7interpreted as the jet energy correction to be applied for that
     8given pT(Reco), |eta| bin.
     9
     10This can be done by modifying the "ScaleFormula" input parameter to
     11the JetEnergyScale module in the delphes_card_XXX.tcl
    1212
    1313e.g  a smooth function:
    1414
    15 
    16   set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0 ) }
    17 
     15  set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
    1816
    1917or a binned function:
    20 
    2118
    2219  set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0)  * (1.10) +
    … …  
    2724                    (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0)               * (1.00)}
    2825
    29 
    30 Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ...
    31 
    32 
    33 
     26Be aware that a binned jet energy scale can produce "steps" in the corrected
     27jet pt distribution ...
    3428*/
     29
     30#ifdef __CLING__
     31R__LOAD_LIBRARY(libDelphes)
     32#include "classes/DelphesClasses.h"
     33#include "external/ExRootAnalysis/ExRootTreeReader.h"
     34#include "external/ExRootAnalysis/ExRootResult.h"
     35#else
     36class ExRootTreeReader;
     37class ExRootResult;
     38#endif
    3539
    3640//------------------------------------------------------------------------------
    … …  
    161165
    162166  Jet *jet, *genjet;
    163   GenParticle *part;
     167  GenParticle *particle;
    164168  TObject *object;
    165169
    166   TLorentzVector JetMom, GenJetMom, BestGenJetMom;
    167 
    168   Float_t Dr;
     170  TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
     171
     172  Float_t deltaR;
    169173  Float_t pt, eta;
    170174  Long64_t entry;
    … …  
    177181    // Load selected branches with data from specified event
    178182    treeReader->ReadEntry(entry);
    179     //  cout<<"--  New event -- "<<endl;
    180183
    181184    if(entry%500 == 0) cout << "Event number: "<< entry <<endl;
    … …  
    186189
    187190      jet = (Jet*) branchJet->At(i);
    188       JetMom = jet-> P4();
    189 
    190       plots->fJetPT->Fill(JetMom.Pt());
    191 
    192       Dr = 999;
     191      jetMomentum = jet->P4();
     192
     193      plots->fJetPT->Fill(jetMomentum.Pt());
     194
     195      deltaR = 999;
    193196
    194197     // Loop over all hard partons in event
    195198     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
    196199     {
    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 )
     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)
    207210        {
    208            Dr = GenJetMom.DeltaR(JetMom);
    209            BestGenJetMom = GenJetMom;
     211           deltaR = genJetMomentum.DeltaR(jetMomentum);
     212           bestGenJetMomentum = genJetMomentum;
    210213        }
    211 
    212214      }
    213215
    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 
     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      }
    239238    }
    240239  }
    241240}
    242 
    243241
    244242//------------------------------------------------------------------------------

  • modules/BTagging.cc

    re33e6db r9343566  
    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, *parton;
     172  Candidate *jet;
    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;
    178177
    179178  // select quark and gluons
    180179  fFilter->Reset();
    181180  partonArray = fFilter->GetSubArray(fClassifier, 0);
    182  
     181
    183182  if(partonArray == 0) return;
    184183
    185184  TIter itPartonArray(partonArray);
    186  
     185
    187186  // loop over all input jets
    188187  fItJetInputArray->Reset();
    … …  
    190189  {
    191190    const TLorentzVector &jetMomentum = jet->Momentum;
    192     pdgCodeMax = -1;
    193191    eta = jetMomentum.Eta();
    194192    phi = jetMomentum.Phi();
    195193    pt = jetMomentum.Pt();
    196194    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;
    211195
    212196    // find an efficency formula
    213     itEfficiencyMap = fEfficiencyMap.find(pdgCodeMax);
     197    itEfficiencyMap = fEfficiencyMap.find(jet->Flavor);
    214198    if(itEfficiencyMap == fEfficiencyMap.end())
    215199    {
    … …  
    220204    // apply an efficency formula
    221205    jet->BTag |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
    222   }
    223 }
    224 
    225 //------------------------------------------------------------------------------
     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//------------------------------------------------------------------------------

  • modules/Calorimeter.cc

    re33e6db r9343566  
    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 
    152160  // read min E value for towers to be saved
    153161  fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
    … …  
    356364      fTrackHCalEnergy = 0.0;
    357365
    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 
    367366      fTowerTrackHits = 0;
    368367      fTowerPhotonHits = 0;
    … …  
    380379      position = track->Position;
    381380
    382 
    383381      ecalEnergy = momentum.E() * fTrackECalFractions[number];
    384382      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
    … …  
    387385      fTrackHCalEnergy += hcalEnergy;
    388386
    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);
     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
    394409
    395410      fTowerTrackArray->Add(track);
    … …  
    397412      continue;
    398413    }
    399 
     414 
    400415    // check for photon and electron hits in current tower
    401416    if(flags & 2) ++fTowerPhotonHits;
    … …  
    412427    fTowerHCalEnergy += hcalEnergy;
    413428
    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 
     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    }
    420442
    421443    fTower->AddCandidate(particle);
    … …  
    434456  Double_t ecalEnergy, hcalEnergy;
    435457  Double_t ecalSigma, hcalSigma;
    436   Double_t ecalTime, hcalTime, time;
    437 
     458 
    438459  if(!fTower) return;
    439460
    … …  
    444465  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
    445466
    446   ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
    447   hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
    448 
    449467  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    450468  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    … …  
    454472
    455473  energy = ecalEnergy + hcalEnergy;
    456   time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
    457 
     474   
    458475  if(fSmearTowerCenter)
    459476  {
    … …  
    469486  pt = energy / TMath::CosH(eta);
    470487
    471   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
     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
    472508  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    473509  fTower->Eem = ecalEnergy;

  • modules/Calorimeter.h

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

  • modules/FastJetFinder.cc

    re33e6db r9343566  
    6666#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
    6767
     68#include "fastjet/tools/Filter.hh"
     69#include "fastjet/tools/Pruner.hh"
     70#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
     71
    6872using namespace std;
    6973using namespace fastjet;
    … …  
    121125  fN = GetInt("N", 2);                  // used only if Njettiness is used as jet clustering algo (case 8)
    122126
     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 
     148
    123149  // ---  Jet Area Parameters ---
    124150  fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
    … …  
    260286  PseudoJet jet, area;
    261287  ClusterSequence *sequence;
    262   vector< PseudoJet > inputList, outputList;
     288  vector< PseudoJet > inputList, outputList, subjets;
    263289  vector< PseudoJet >::iterator itInputList, itOutputList;
    264290  vector< TEstimatorStruct >::iterator itEstimators;
    … …  
    352378    candidate->DeltaEta = detaMax;
    353379    candidate->DeltaPhi = dphiMax;
    354 
     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 
    355461    // --- compute N-subjettiness with N = 1,2,3,4,5 ----
    356462

  • modules/FastJetFinder.h

    re33e6db r9343566  
    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
    85105  // --- FastJet Area method --------
    86106

  • modules/Isolation.cc

    re33e6db r9343566  
    2525 *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
    2626 *
    27  *  \author P. Demin - UCL, Louvain-la-Neuve
     27 *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
    2828 *
    2929 */
    … …  
    109109  fUsePTSum = GetBool("UsePTSum", false);
    110110
     111  fVetoLeptons = GetBool("VetoLeptons", true); 
     112 
    111113  fClassifier->fPTMin = GetDouble("PTMin", 0.5);
     114
    112115
    113116  // import input array(s)
    … …  
    153156  Candidate *candidate, *isolation, *object;
    154157  TObjArray *isolationArray;
    155   Double_t sum, ratio;
     158  Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
    156159  Int_t counter;
    157160  Double_t eta = 0.0;
    … …  
    180183
    181184    // loop over all input tracks
    182     sum = 0.0;
     185   
     186    sumNeutral = 0.0;
     187    sumCharged = 0.0;
     188    sumChargedPU = 0.0;
     189    sumAllParticles = 0.0;
     190   
    183191    counter = 0;
    184192    itIsolationArray.Reset();
     193   
    185194    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
    186195    {
    … …  
    188197
    189198      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
    190          candidate->GetUniqueID() != isolation->GetUniqueID())
     199         candidate->GetUniqueID() != isolation->GetUniqueID() &&
     200         ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) )
    191201      {
    192         sum += isolationMomentum.Pt();
     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 
    193212        ++counter;
    194213      }
    … …  
    209228    }
    210229
    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 
     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;
    217245    fOutputArray->Add(candidate);
    218246  }

  • modules/Isolation.h

    re33e6db r9343566  
    5959  Bool_t fUsePTSum;
    6060
     61  Bool_t fVetoLeptons; 
     62 
    6163  IsolationClassifier *fClassifier; //!
    6264

  • modules/ModulesLinkDef.h

    re33e6db r9343566  
    5858#include "modules/Weighter.h"
    5959#include "modules/Hector.h"
     60#include "modules/RunPUPPI.h"
     61#include "modules/JetFlavorAssociation.h"
    6062#include "modules/ExampleModule.h"
    6163
    … …  
    98100#pragma link C++ class Weighter+;
    99101#pragma link C++ class Hector+;
     102#pragma link C++ class RunPUPPI+;
     103#pragma link C++ class JetFlavorAssociation+;
    100104#pragma link C++ class ExampleModule+;
    101105

  • modules/PdgCodeFilter.cc

    re33e6db r9343566  
    7474  fPTMin = GetDouble("PTMin", 0.0);
    7575
     76  fInvertPdg = GetBool("InvertPdg", false);
     77
     78  fRequireStatus = GetBool("RequireStatus", false);
     79  fStatus = GetInt("Status", 1);
     80
    7681  // import input array
    7782  fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
    … …  
    109114  Double_t pt;
    110115
     116  const Bool_t requireStatus = fRequireStatus;
     117  const Bool_t invertPdg = fInvertPdg;
     118  const int reqStatus = fStatus;
     119
    111120  fItInputArray->Reset();
    112121  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    … …  
    116125    pt = candidateMomentum.Pt();
    117126
     127    if(pt < fPTMin) continue;
     128    if(requireStatus && (candidate->Status != reqStatus)) continue;
     129
    118130    pass = kTRUE;
    119 
    120     if(pt < fPTMin) pass = kFALSE;
    121131    if(find(fPdgCodes.begin(), fPdgCodes.end(), pdgCode) != fPdgCodes.end()) pass = kFALSE;
    122132
     133    if (invertPdg) pass = !pass;
    123134    if(pass) fOutputArray->Add(candidate);
    124135  }

  • modules/PdgCodeFilter.h

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

  • modules/PileUpJetID.cc

    re33e6db r9343566  
    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 
    191
    202/** \class PileUpJetID
    213 *
    22  *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
     4 *  CMS PileUp Jet ID Variables
    235 *
    24  *  \author S. Zenz, December 2013
     6 *  \author S. Zenz
    257 *
    268 */
    … …  
    4123#include "TRandom3.h"
    4224#include "TObjArray.h"
    43 #include "TDatabasePDG.h"
     25//#include "TDatabasePDG.h"
    4426#include "TLorentzVector.h"
    4527
    … …  
    5436
    5537PileUpJetID::PileUpJetID() :
    56   fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
     38  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
    5739{
    5840
    … …  
    7052void PileUpJetID::Init()
    7153{
     54
    7255  fJetPTMin = GetDouble("JetPTMin", 20.0);
    7356  fParameterR = GetDouble("ParameterR", 0.5);
    7457  fUseConstituents = GetInt("UseConstituents", 0);
    7558
     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
    7686  fAverageEachTower = false; // for timing
    7787
    … …  
    8191  fItJetInputArray = fJetInputArray->MakeIterator();
    8292
    83   fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
     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"));
    8499  fItTrackInputArray = fTrackInputArray->MakeIterator();
    85100
    86   fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
     101  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
    87102  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    88103
    89   fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
    90   fItVertexInputArray = fVertexInputArray->MakeIterator();
    91 
    92   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    93104
    94105  // create output array(s)
    95106
    96107  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
     108
     109  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
     110
     111
     112  //  cout << " end of INIT " << endl;
     113
    97114}
    98115
    … …  
    101118void PileUpJetID::Finish()
    102119{
     120  //  cout << "In finish" << endl;
     121
    103122  if(fItJetInputArray) delete fItJetInputArray;
    104123  if(fItTrackInputArray) delete fItTrackInputArray;
    105124  if(fItNeutralInputArray) delete fItNeutralInputArray;
    106   if(fItVertexInputArray) delete fItVertexInputArray;
     125
    107126}
    108127
    … …  
    111130void PileUpJetID::Process()
    112131{
     132  //  cout << "start of process" << endl;
     133
    113134  Candidate *candidate, *constituent;
    114135  TLorentzVector momentum, area;
    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   }
     136
     137  //  cout << "BeforE SCZ additions in process" << endl;
     138
     139  // SCZ
     140  Candidate *trk;
    133141
    134142  // loop over all input candidates
    … …  
    139147    area = candidate->Area;
    140148
    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     {
     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) {
    157173      TIter itConstituents(candidate->GetCandidates());
    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     {
     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 {
    195227      // Not using constituents, using dr
    196228      fItTrackInputArray->Reset();
    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 
     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      }
    226250      fItNeutralInputArray->Reset();
    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     {
     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.) {
    260276      candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
    261     }
    262     else
    263     {
    264       candidate->MeanSqDeltaR = -999.0;
     277    } else {
     278      candidate->MeanSqDeltaR = -999.;
    265279    }
    266280    candidate->NCharged = nc;
    267281    candidate->NNeutrals = nn;
    268     if(sumpt > 0.0)
    269     {
     282    if (sumpt > 0.) {
    270283      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
    271       for(i = 0; i < 5; ++i)
    272       {
     284      for (int i = 0 ; i < 5 ; i++) {
    273285        candidate->FracPt[i] = pt_ann[i]/sumpt;
    274286      }
    275     }
    276     else
    277     {
    278       candidate->PTD = -999.0;
    279       for(i = 0; i < 5; ++i)
    280       {
    281         candidate->FracPt[i] = -999.0;
     287    } else {
     288      candidate->PTD = -999.;
     289      for (int i = 0 ; i < 5 ; i++) {
     290        candidate->FracPt[i] = -999.;
    282291      }
    283292    }
    284293
    285294    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
    286340  }
    287341}
    288342
    289343//------------------------------------------------------------------------------
    290 

  • modules/PileUpJetID.h

    re33e6db r9343566  
    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 
    191#ifndef PileUpJetID_h
    202#define PileUpJetID_h
    … …  
    224/** \class PileUpJetID
    235 *
    24  *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
     6 *  CMS PileUp Jet ID Variables
    257 *
    26  *  \author S. Zenz, December 2013
     8 *  \author S. Zenz
    279 *
    2810 */
    … …  
    3416
    3517class TObjArray;
     18class DelphesFormula;
    3619
    3720class PileUpJetID: public DelphesModule
    … …  
    5134  Double_t fParameterR;
    5235
     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  /*
     46JAY
     47---
     48
     49|Eta|<1.5
     50
     51meanSqDeltaR betaStar SigEff BgdEff
     520.13 0.92 96% 8%
     530.13 0.95 97% 16%
     540.13 0.97 98% 27%
     55
     56|Eta|>1.5
     57
     58meanSqDeltaR betaStar SigEff BgdEff
     590.14 0.91 95% 15%
     600.14 0.94 97% 19%
     610.14 0.97 98% 29%
     62
     63BRYAN
     64-----
     65
     66Barrel (MeanSqDR, Beta, sig eff, bg eff):
     670.10, 0.08, 90%, 8%
     680.11, 0.12, 90%, 6%
     690.13, 0.16, 89%, 5%
     70
     71Endcap (MeanSqDR, Beta, sig eff, bg eff):
     720.07, 0.06, 89%, 4%
     730.08, 0.08, 92%, 6%
     740.09, 0.08, 95%, 10%
     750.10, 0.08, 97%, 13%
     76
     77SETH GUESSES FOR |eta| > 4.0
     78----------------------------
     79
     80MeanSqDeltaR
     810.07
     820.10
     830.14
     840.2
     85  */
     86
    5387  // If set to true, may have weird results for PFCHS
    5488  // If set to false, uses everything within dR < fParameterR even if in other jets &c.
    5589  // Results should be very similar for PF
    56   Int_t fUseConstituents;
     90  Int_t fUseConstituents; 
    5791
    5892  Bool_t fAverageEachTower;
    … …  
    6296  const TObjArray *fJetInputArray; //!
    6397
    64   const TObjArray *fTrackInputArray; //!
    65   const TObjArray *fNeutralInputArray; //!
     98  const TObjArray *fTrackInputArray; // SCZ
     99  const TObjArray *fNeutralInputArray;
    66100
    67   TIterator *fItTrackInputArray; //!
    68   TIterator *fItNeutralInputArray; //!
     101  TIterator *fItTrackInputArray; // SCZ
     102  TIterator *fItNeutralInputArray; // SCZ
    69103
    70104  TObjArray *fOutputArray; //!
     105  TObjArray *fNeutralsInPassingJets; // SCZ
    71106
    72   TIterator *fItVertexInputArray; //!
    73   const TObjArray *fVertexInputArray; //!
    74107
    75   Double_t fZVertexResolution;
    76 
    77   ClassDef(PileUpJetID, 1)
     108  ClassDef(PileUpJetID, 2)
    78109};
    79110
    80111#endif
    81 

  • modules/PileUpMerger.cc

    re33e6db r9343566  
    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
    8287  // read vertex smearing formula
    8388
    … …  
    111116  TParticlePDG *pdgParticle;
    112117  Int_t pid;
    113   Float_t x, y, z, t;
     118  Float_t x, y, z, t, vx, vy;
    114119  Float_t px, py, pz, e;
    115120  Double_t dz, dphi, dt;
    116   Int_t numberOfEvents, event;
     121  Int_t numberOfEvents, event, numberOfParticles;
    117122  Long64_t allEntries, entry;
    118   Candidate *candidate, *vertexcandidate;
     123  Candidate *candidate, *vertex;
    119124  DelphesFactory *factory;
    120125
    … …  
    123128  fItInputArray->Reset();
    124129
    125   // --- Deal with Primary vertex first  ------
     130  // --- Deal with primary vertex first  ------
    126131
    127132  fFunction->GetRandom2(dz, dt);
    … …  
    129134  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    130135  dz *= 1.0E3; // necessary in order to make z in mm
    131 
     136  vx = 0.0;
     137  vy = 0.0;
     138  numberOfParticles = fInputArray->GetEntriesFast();
    132139  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    133140  {
     141    vx += candidate->Position.X();
     142    vy += candidate->Position.Y();
    134143    z = candidate->Position.Z();
    135144    t = candidate->Position.T();
    … …  
    139148  }
    140149
     150  if(numberOfParticles > 0)
     151  {
     152    vx /= numberOfParticles;
     153    vy /= numberOfParticles;
     154  }
     155
    141156  factory = GetFactory();
    142157
    143   vertexcandidate = factory->NewCandidate();
    144   vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
    145   fVertexOutputArray->Add(vertexcandidate);
     158  vertex = factory->NewCandidate();
     159  vertex->Position.SetXYZT(vx, vy, dz, dt);
     160  fVertexOutputArray->Add(vertex);
    146161
    147162  // --- Then with pile-up vertices  ------
    … …  
    181196    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    182197
    183     vertexcandidate = factory->NewCandidate();
    184     vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
    185     vertexcandidate->IsPU = 1;
    186 
    187     fVertexOutputArray->Add(vertexcandidate);
    188 
     198    vx = 0.0;
     199    vy = 0.0;
     200    numberOfParticles = 0;
    189201    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
    190202    {
    … …  
    204216      candidate->Momentum.RotateZ(dphi);
    205217
     218      x -= fInputBeamSpotX;
     219      y -= fInputBeamSpotY;
    206220      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
    207221      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;
    208227
    209228      fParticleOutputArray->Add(candidate);
    210229    }
    211   }
    212 }
    213 
    214 //------------------------------------------------------------------------------
    215 
     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);
     242  }
     243}
     244
     245//------------------------------------------------------------------------------

  • modules/PileUpMerger.h

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

  • modules/PileUpMergerPythia8.cc

    re33e6db r9343566  
    2929#include "classes/DelphesClasses.h"
    3030#include "classes/DelphesFactory.h"
    31 #include "classes/DelphesFormula.h"
     31#include "classes/DelphesTF2.h"
    3232#include "classes/DelphesPileUpReader.h"
    3333
    … …  
    5656
    5757PileUpMergerPythia8::PileUpMergerPythia8() :
    58   fPythia(0), fItInputArray(0)
    59 {
     58  fFunction(0), fPythia(0), fItInputArray(0)
     59{
     60  fFunction = new DelphesTF2;
    6061}
    6162
    … …  
    6465PileUpMergerPythia8::~PileUpMergerPythia8()
    6566{
     67  delete fFunction;
    6668}
    6769
    … …  
    7274  const char *fileName;
    7375
     76  fPileUpDistribution = GetInt("PileUpDistribution", 0);
     77
    7478  fMeanPileUp  = GetDouble("MeanPileUp", 10);
    75   fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
     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);
    7687
    7788  fPTMin = GetDouble("PTMin", 0.0);
     89
     90  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
     91  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
    7892
    7993  fileName = GetString("ConfigFile", "MinBias.cmnd");
    … …  
    86100
    87101  // create output arrays
    88   fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
     102  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
     103  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
    89104}
    90105
    … …  
    103118  TParticlePDG *pdgParticle;
    104119  Int_t pid, status;
    105   Float_t x, y, z, t;
     120  Float_t x, y, z, t, vx, vy;
    106121  Float_t px, py, pz, e;
    107   Double_t dz, dphi;
    108   Int_t poisson, event, i;
    109   Candidate *candidate;
     122  Double_t dz, dphi, dt;
     123  Int_t numberOfEvents, event, numberOfParticles, i;
     124  Candidate *candidate, *vertex;
    110125  DelphesFactory *factory;
    111126
     127  const Double_t c_light = 2.99792458E8;
     128
    112129  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();
    113140  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    114141  {
    115     fOutputArray->Add(candidate);
     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;
    116155  }
    117156
    118157  factory = GetFactory();
    119158
    120   poisson = gRandom->Poisson(fMeanPileUp);
    121 
    122   for(event = 0; event < poisson; ++event)
     159  vertex = factory->NewCandidate();
     160  vertex->Position.SetXYZT(vx, vy, dz, dt);
     161  fVertexOutputArray->Add(vertex);
     162
     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)
    123179  {
    124180    while(!fPythia->next());
    125181
    126     dz = gRandom->Gaus(0.0, fZVertexSpread);
     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
    127189    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    128190
    129     for(i = 1; i < fPythia->event.size(); ++i)
     191    vx = 0.0;
     192    vy = 0.0;
     193    numberOfParticles = fPythia->event.size();
     194    for(i = 1; i < numberOfParticles; ++i)
    130195    {
    131196      Pythia8::Particle &particle = fPythia->event[i];
    … …  
    143208      candidate->PID = pid;
    144209
    145       candidate->Status = status;
     210      candidate->Status = 1;
    146211
    147212      pdgParticle = pdg->GetParticle(pid);
    … …  
    154219      candidate->Momentum.RotateZ(dphi);
    155220
    156       candidate->Position.SetXYZT(x, y, z + dz, t);
     221      x -= fInputBeamSpotX;
     222      y -= fInputBeamSpotY;
     223      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
    157224      candidate->Position.RotateZ(dphi);
    158 
    159       fOutputArray->Add(candidate);
     225      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
     226
     227      vx += candidate->Position.X();
     228      vy += candidate->Position.Y();
     229
     230      fParticleOutputArray->Add(candidate);
    160231    }
    161   }
    162 }
    163 
    164 //------------------------------------------------------------------------------
    165 
     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);
     244  }
     245}
     246
     247//------------------------------------------------------------------------------

  • modules/PileUpMergerPythia8.h

    re33e6db r9343566  
    3131
    3232class TObjArray;
     33class DelphesTF2;
    3334
    3435namespace Pythia8
    … …  
    5051private:
    5152
     53  Int_t fPileUpDistribution;
    5254  Double_t fMeanPileUp;
     55
    5356  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
    5464  Double_t fPTMin;
     65
     66  DelphesTF2 *fFunction; //!
    5567
    5668  Pythia8::Pythia *fPythia; //!
    … …  
    6072  const TObjArray *fInputArray; //!
    6173
    62   TObjArray *fOutputArray; //!
     74  TObjArray *fParticleOutputArray; //!
     75  TObjArray *fVertexOutputArray; //!
    6376
    6477  ClassDef(PileUpMergerPythia8, 1)

  • modules/SimpleCalorimeter.cc

    re33e6db r9343566  
    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);
    151154
    152155  // switch on or off the dithering of the center of calorimeter towers
    … …  
    425428  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    426429
     430  fTower->Eem = (!fIsEcal) ? 0 : energy;
     431  fTower->Ehad = (fIsEcal) ? 0 : energy;
     432
    427433  fTower->Edges[0] = fTowerEdges[0];
    428434  fTower->Edges[1] = fTowerEdges[1];
    … …  
    447453    pt = energy / TMath::CosH(eta);
    448454
     455    tower->Eem = (!fIsEcal) ? 0 : energy;
     456    tower->Ehad = (fIsEcal) ? 0 : energy;
     457
    449458    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    450459    fEFlowTowerOutputArray->Add(tower);

  • modules/SimpleCalorimeter.h

    re33e6db r9343566  
     1
    12/*
    23 *  Delphes: a framework for fast simulation of a generic collider experiment
    … …  
    7475  Bool_t fSmearTowerCenter;
    7576
     77  Bool_t fIsEcal; //!
     78
    7679  TFractionMap fFractionMap; //!
    7780  TBinMap fBinMap; //!

  • modules/TrackPileUpSubtractor.cc

    re33e6db r9343566  
    7474  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    7575
     76  fPTMin = GetDouble("PTMin", 0.);
    7677  // import arrays with output from other modules
    77 
     78   
    7879  ExRootConfParam param = GetParam("InputArray");
    7980  Long_t i, size;
    … …  
    147148      // assume perfect pile-up subtraction for tracks outside fZVertexResolution
    148149     
    149       if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
    150 
    151       array->Add(candidate);
     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      }
    152156    }
    153157  }

  • modules/TrackPileUpSubtractor.h

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

  • modules/TreeWriter.cc

    re33e6db r9343566  
    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;
    364365
    365366    FillParticles(candidate, &entry->Particles);
    … …  
    403404    entry->T = position.T()*1.0E-3/c_light;
    404405
     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
    405415    entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
    406416
    … …  
    442452    entry->T = position.T()*1.0E-3/c_light;
    443453
     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
    444464    entry->Charge = candidate->Charge;
    445465
    … …  
    469489    const TLorentzVector &momentum = candidate->Momentum;
    470490    const TLorentzVector &position = candidate->Position;
    471 
    472491
    473492    pt = momentum.Pt();
    … …  
    488507    entry->T = position.T()*1.0E-3/c_light;
    489508
     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
    490518    entry->Charge = candidate->Charge;
    491519
    … …  
    504532  Double_t ecalEnergy, hcalEnergy;
    505533  const Double_t c_light = 2.99792458E8;
     534  Int_t i;
    506535
    507536  array->Sort();
    … …  
    532561    entry->Mass = momentum.M();
    533562
     563    entry->Area = candidate->Area;
     564
    534565    entry->DeltaEta = candidate->DeltaEta;
    535566    entry->DeltaPhi = candidate->DeltaPhi;
    536567
     568    entry->Flavor = candidate->Flavor;
     569    entry->FlavorAlgo = candidate->FlavorAlgo;
     570    entry->FlavorPhys = candidate->FlavorPhys;
     571
    537572    entry->BTag = candidate->BTag;
     573
     574    entry->BTagAlgo = candidate->BTagAlgo;
     575    entry->BTagPhys = candidate->BTagPhys;
     576
    538577    entry->TauTag = candidate->TauTag;
    539578
    … …  
    561600    entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
    562601    entry->PTD = candidate->PTD;
    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    
     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
    577618    FillParticles(candidate, &entry->Particles);
    578619  }

  • readers/DelphesLHEF.cpp

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

  • readers/DelphesPythia8.cpp

    re33e6db r9343566  
    2020#include <iostream>
    2121#include <sstream>
     22#include <string>
    2223
    2324#include <signal.h>
    … …  
    3839#include "classes/DelphesClasses.h"
    3940#include "classes/DelphesFactory.h"
     41#include "classes/DelphesLHEFReader.h"
    4042
    4143#include "ExRootAnalysis/ExRootTreeWriter.h"
    … …  
    149151  char appName[] = "DelphesPythia8";
    150152  stringstream message;
     153  FILE *inputFile = 0;
    151154  TFile *outputFile = 0;
    152155  TStopwatch readStopWatch, procStopWatch;
    153156  ExRootTreeWriter *treeWriter = 0;
    154157  ExRootTreeBranch *branchEvent = 0;
     158  ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0;
    155159  ExRootConfReader *confReader = 0;
    156160  Delphes *modularDelphes = 0;
    157161  DelphesFactory *factory = 0;
    158162  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
     163  TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0;
     164  DelphesLHEFReader *reader = 0;
    159165  Long64_t eventCounter, errorCounter;
    160166  Long64_t numberOfEvents, timesAllowErrors;
    … …  
    205211    partonOutputArray = modularDelphes->ExportArray("partons");
    206212
    207     modularDelphes->InitTask();
    208 
    209     // Initialize pythia
     213    // Initialize Pythia
    210214    pythia = new Pythia8::Pythia;
    211215
    … …  
    221225    numberOfEvents = pythia->mode("Main:numberOfEvents");
    222226    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();
    223243
    224244    pythia->init();
    … …  
    234254    for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
    235255    {
     256      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
     257        stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
     258
    236259      if(!pythia->next())
    237260      {
    238261        // If failure because reached end of file then exit event loop
    239         if (pythia->info.atEndOfFile())
     262        if(pythia->info.atEndOfFile())
    240263        {
    241264          cerr << "Aborted since reached end of Les Houches Event File" << endl;
    … …  
    244267
    245268        // First few failures write off as "acceptable" errors, then quit
    246         if (++errorCounter < timesAllowErrors) continue;
    247         cerr << "Event generation aborted prematurely, owing to error!" << endl;
    248         break;
     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;
    249278      }
    250279
    … …  
    258287      procStopWatch.Stop();
    259288
     289      if(reader)
     290      {
     291        reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch);
     292        reader->AnalyzeWeight(branchWeightLHEF);
     293      }
     294
    260295      treeWriter->Fill();
    261296
    262297      treeWriter->Clear();
    263298      modularDelphes->Clear();
     299      if(reader) reader->Clear();
    264300
    265301      readStopWatch.Start();
    … …  
    277313    cout << "** Exiting..." << endl;
    278314
     315    delete reader;
    279316    delete pythia;
    280317    delete modularDelphes;
Note: See TracChangeset for help on using the changeset viewer.