Fork me on GitHub

Changeset 28027d5 in git


Ignore:
Timestamp:
Jun 26, 2015, 3:13:55 PM (9 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
f3c4047
Parents:
f53a4d2 (diff), fe0273c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge pull request #15 from delphes/mergeCMS

merge with the CMS code

Files:
27 added
26 edited

Legend:

Unmodified
Added
Removed
  • Makefile

    rf53a4d2 r28027d5  
    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
  • classes/ClassesLinkDef.h

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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/DelphesLHEFReader.cc

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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
  • doc/genMakefile.tcl

    rf53a4d2 r28027d5  
    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}
  • modules/BTagging.cc

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    5959  Bool_t fUsePTSum;
    6060
     61  Bool_t fVetoLeptons; 
     62 
    6163  IsolationClassifier *fClassifier; //!
    6264
  • modules/ModulesLinkDef.h

    rf53a4d2 r28027d5  
    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/PileUpJetID.cc

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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/TrackPileUpSubtractor.cc

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    5050  Double_t fZVertexResolution;
    5151
     52  Double_t fPTMin;
     53
    5254  std::map< TIterator *, TObjArray * > fInputMap; //!
    5355
  • modules/TreeWriter.cc

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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

    rf53a4d2 r28027d5  
    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.