Changes in / [d38348d:f53a4d2] in git

Files:

• Makefile (modified) (18 diffs)
• classes/ClassesLinkDef.h (modified) (1 diff)
• classes/DelphesClasses.cc (modified) (10 diffs)
• classes/DelphesClasses.h (modified) (11 diffs)
• classes/DelphesLHEFReader.cc (modified) (6 diffs)
• classes/DelphesLHEFReader.h (modified) (3 diffs)
• doc/genMakefile.tcl (modified) (1 diff)
• examples/Pythia8/configLHE.cmnd (deleted)
• examples/Pythia8/configNoLHE.cmnd (deleted)
• examples/Pythia8/events.lhe (deleted)
• external/PUPPI/RecoObj.hh (deleted)
• external/PUPPI/puppiAlgoBin.hh (deleted)
• external/PUPPI/puppiCleanContainer.cc (deleted)
• external/PUPPI/puppiCleanContainer.hh (deleted)
• external/PUPPI/puppiParticle.hh (deleted)
• external/fastjet/contribs/RecursiveTools/AUTHORS (deleted)
• external/fastjet/contribs/RecursiveTools/COPYING (deleted)
• external/fastjet/contribs/RecursiveTools/ChangeLog (deleted)
• external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.cc (deleted)
• external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.hh (deleted)
• external/fastjet/contribs/RecursiveTools/NEWS (deleted)
• external/fastjet/contribs/RecursiveTools/README (deleted)
• external/fastjet/contribs/RecursiveTools/Recluster.cc (deleted)
• external/fastjet/contribs/RecursiveTools/Recluster.hh (deleted)
• external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc (deleted)
• external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.hh (deleted)
• external/fastjet/contribs/RecursiveTools/SoftDrop.cc (deleted)
• external/fastjet/contribs/RecursiveTools/SoftDrop.hh (deleted)
• external/fastjet/contribs/RecursiveTools/TODO (deleted)
• external/fastjet/contribs/RecursiveTools/VERSION (deleted)
• modules/BTaggingCMS.cc (deleted)
• modules/BTaggingCMS.h (deleted)
• modules/Calorimeter.cc (modified) (10 diffs)
• modules/Calorimeter.h (modified) (1 diff)
• modules/FastJetFinder.cc (modified) (4 diffs)
• modules/FastJetFinder.h (modified) (1 diff)
• modules/Isolation.cc (modified) (6 diffs)
• modules/Isolation.h (modified) (1 diff)
• modules/JetFlavorAssociation.cc (deleted)
• modules/JetFlavorAssociation.h (deleted)
• modules/ModulesLinkDef.h (modified) (4 diffs)
• modules/PileUpJetID.cc (modified) (8 diffs)
• modules/PileUpJetID.h (modified) (5 diffs)
• modules/PileUpMerger.cc (modified) (7 diffs)
• modules/PileUpMerger.h (modified) (1 diff)
• modules/PileUpMergerPythia8.cc (modified) (8 diffs)
• modules/PileUpMergerPythia8.h (modified) (3 diffs)
• modules/RunPUPPI.cc (deleted)
• modules/RunPUPPI.h (deleted)
• modules/TrackPileUpSubtractor.cc (modified) (2 diffs)
• modules/TrackPileUpSubtractor.h (modified) (1 diff)
• modules/TreeWriter.cc (modified) (9 diffs)
• readers/DelphesLHEF.cpp (modified) (3 diffs)
• readers/DelphesPythia8.cpp (modified) (9 diffs)

Makefile

@@ -257 +257 @@
         classes/DelphesClasses.h \
         classes/DelphesFactory.h \
-        classes/DelphesLHEFReader.h \
         external/ExRootAnalysis/ExRootTreeWriter.h \
         external/ExRootAnalysis/ExRootTreeBranch.h \
@@ -322 +321 @@
         modules/UniqueObjectFinder.h \
         modules/TrackCountingBTagging.h \
-        modules/BTaggingCMS.h \
         modules/BTagging.h \
         modules/TauTagging.h \
@@ -339 +337 @@
         modules/Weighter.h \
         modules/Hector.h \
-        modules/RunPUPPI.h \
-        modules/JetFlavorAssociation.h \
         modules/ExampleModule.h
 ModulesDict$(PcmSuf): \
@@ -525 +521 @@
 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
         external/Hector/H_VerticalQuadrupole.$(SrcSuf)
-tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
-        external/PUPPI/puppiCleanContainer.$(SrcSuf) \
-        external/fastjet/Selector.hh
 tmp/modules/AngularSmearing.$(ObjSuf): \
         modules/AngularSmearing.$(SrcSuf) \
@@ -540 +533 @@
         modules/BTagging.$(SrcSuf) \
         modules/BTagging.h \
-        classes/DelphesClasses.h \
-        classes/DelphesFactory.h \
-        classes/DelphesFormula.h \
-        external/ExRootAnalysis/ExRootResult.h \
-        external/ExRootAnalysis/ExRootFilter.h \
-        external/ExRootAnalysis/ExRootClassifier.h
-tmp/modules/BTaggingCMS.$(ObjSuf): \
-        modules/BTaggingCMS.$(SrcSuf) \
-        modules/BTaggingCMS.h \
         classes/DelphesClasses.h \
         classes/DelphesFactory.h \
@@ -668 +652 @@
         external/ExRootAnalysis/ExRootFilter.h \
         external/ExRootAnalysis/ExRootClassifier.h
-tmp/modules/JetFlavorAssociation.$(ObjSuf): \
-        modules/JetFlavorAssociation.$(SrcSuf) \
-        modules/JetFlavorAssociation.h \
-        classes/DelphesClasses.h \
-        classes/DelphesFactory.h \
-        classes/DelphesFormula.h \
-        external/ExRootAnalysis/ExRootResult.h \
-        external/ExRootAnalysis/ExRootFilter.h \
-        external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/JetPileUpSubtractor.$(ObjSuf): \
         modules/JetPileUpSubtractor.$(SrcSuf) \
@@ -755 +730 @@
         classes/DelphesClasses.h \
         classes/DelphesFactory.h \
-        classes/DelphesTF2.h \
+        classes/DelphesFormula.h \
         classes/DelphesPileUpReader.h \
         external/ExRootAnalysis/ExRootResult.h \
         external/ExRootAnalysis/ExRootFilter.h \
         external/ExRootAnalysis/ExRootClassifier.h
-tmp/modules/RunPUPPI.$(ObjSuf): \
-        modules/RunPUPPI.$(SrcSuf) \
-        modules/RunPUPPI.h \
-        external/PUPPI/puppiCleanContainer.hh \
-        external/PUPPI/RecoObj.hh \
-        external/PUPPI/puppiParticle.hh \
-        external/PUPPI/puppiAlgoBin.hh \
-        classes/DelphesClasses.h \
-        classes/DelphesFactory.h \
-        classes/DelphesFormula.h
 tmp/modules/SimpleCalorimeter.$(ObjSuf): \
         modules/SimpleCalorimeter.$(SrcSuf) \
@@ -904 +869 @@
         tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
         tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
-        tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
         tmp/modules/AngularSmearing.$(ObjSuf) \
         tmp/modules/BTagging.$(ObjSuf) \
-        tmp/modules/BTaggingCMS.$(ObjSuf) \
         tmp/modules/Calorimeter.$(ObjSuf) \
         tmp/modules/Cloner.$(ObjSuf) \
@@ -920 +883 @@
         tmp/modules/ImpactParameterSmearing.$(ObjSuf) \
         tmp/modules/Isolation.$(ObjSuf) \
-        tmp/modules/JetFlavorAssociation.$(ObjSuf) \
         tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
         tmp/modules/LeptonDressing.$(ObjSuf) \
@@ -929 +891 @@
         tmp/modules/PileUpJetID.$(ObjSuf) \
         tmp/modules/PileUpMerger.$(ObjSuf) \
-        tmp/modules/RunPUPPI.$(ObjSuf) \
         tmp/modules/SimpleCalorimeter.$(ObjSuf) \
         tmp/modules/StatusPidFilter.$(ObjSuf) \
@@ -1119 +1080 @@
 tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf): \
         external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(SrcSuf)
-tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf): \
-        external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(SrcSuf) \
-        external/fastjet/JetDefinition.hh \
-        external/fastjet/ClusterSequenceAreaBase.hh
-tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf): \
-        external/fastjet/contribs/RecursiveTools/Recluster.$(SrcSuf)
-tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf): \
-        external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(SrcSuf) \
-        external/fastjet/JetDefinition.hh \
-        external/fastjet/ClusterSequenceAreaBase.hh
-tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf): \
-        external/fastjet/contribs/RecursiveTools/SoftDrop.$(SrcSuf)
 tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf): \
         external/fastjet/contribs/SoftKiller/SoftKiller.$(SrcSuf)
@@ -1264 +1213 @@
         external/fastjet/contribs/Nsubjettiness/Njettiness.hh \
         external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh \
-        external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh \
-        external/fastjet/tools/Filter.hh \
-        external/fastjet/tools/Pruner.hh \
-        external/fastjet/contribs/RecursiveTools/SoftDrop.hh
+        external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh
 tmp/modules/FastJetGridMedianEstimator.$(ObjSuf): \
         modules/FastJetGridMedianEstimator.$(SrcSuf) \
@@ -1339 +1285 @@
         tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf) \
         tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf) \
-        tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf) \
-        tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf) \
-        tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf) \
-        tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf) \
         tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf) \
         tmp/external/fastjet/plugins/ATLASCone/ATLASConePlugin.$(ObjSuf) \
@@ -1599 +1541 @@
         @touch $@
 
+external/fastjet/Selector.hh: \
+        external/fastjet/PseudoJet.hh \
+        external/fastjet/RangeDefinition.hh
+        @touch $@
+
 external/fastjet/internal/Dnn2piCylinder.hh: \
         external/fastjet/internal/DynamicNearestNeighbours.hh \
         external/fastjet/internal/DnnPlane.hh \
         external/fastjet/internal/numconsts.hh
-        @touch $@
-
-external/fastjet/Selector.hh: \
-        external/fastjet/PseudoJet.hh \
-        external/fastjet/RangeDefinition.hh
         @touch $@
 
@@ -1697 +1639 @@
         @touch $@
 
-modules/RunPUPPI.h: \
-        classes/DelphesModule.h
-        @touch $@
-
 modules/Cloner.h: \
         classes/DelphesModule.h
@@ -1773 +1711 @@
         @touch $@
 
-modules/BTaggingCMS.h: \
-        classes/DelphesModule.h \
-        classes/DelphesClasses.h
-        @touch $@
-
 modules/TrackCountingBTagging.h: \
         classes/DelphesModule.h
@@ -1790 +1723 @@
         external/fastjet/ClusterSequenceAreaBase.hh \
         external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh
-        @touch $@
-
-modules/JetFlavorAssociation.h: \
-        classes/DelphesModule.h \
-        classes/DelphesClasses.h
         @touch $@
 
@@ -1829 +1757 @@
         external/fastjet/AreaDefinition.hh \
         external/fastjet/ClusterSequenceAreaBase.hh
-        @touch $@
-
-external/PUPPI/puppiCleanContainer.hh: \
-        external/PUPPI/RecoObj.hh \
-        external/PUPPI/puppiParticle.hh \
-        external/PUPPI/puppiAlgoBin.hh \
-        external/fastjet/internal/base.hh \
-        external/fastjet/PseudoJet.hh
         @touch $@
 

classes/ClassesLinkDef.h

@@ -46 +46 @@
 #pragma link C++ class LHCOEvent+;
 #pragma link C++ class LHEFEvent+;
-#pragma link C++ class LHEFWeight+;
 #pragma link C++ class HepMCEvent+;
 #pragma link C++ class GenParticle+;

classes/DelphesClasses.cc

@@ -120 +120 @@
   PID(0), Status(0), M1(-1), M2(-1), D1(-1), D2(-1),
   Charge(0), Mass(0.0),
-  IsPU(0), IsRecoPU(0), IsConstituent(0),
-  BTag(0), BTagAlgo(0), BTagDefault(0), BTagPhysics(0), BTagNearest2(0), BTagNearest3(0), BTagHeaviest(0), BTagHighestPt(0),
-  FlavorAlgo(0), FlavorDefault(0), FlavorPhysics(0), FlavorNearest2(0), FlavorNearest3(0), FlavorHeaviest(0), FlavorHighestPt(0),
-  TauTag(0), Eem(0.0), Ehad(0.0),
+  IsPU(0), IsConstituent(0),
+  BTag(0), TauTag(0), Eem(0.0), Ehad(0.0),
   DeltaEta(0.0), DeltaPhi(0.0),
   Momentum(0.0, 0.0, 0.0, 0.0),
@@ -131 +129 @@
   NCharged(0),
   NNeutrals(0),
-  NSubJetsTrimmed(0),
-  NSubJetsPruned(0), 
-  NSubJetsSoftDropped(0), 
   Beta(0),
   BetaStar(0),
   MeanSqDeltaR(0),
   PTD(0),
-  Ntimes(-1),
-  IsolationVar(-999),
-  IsolationVarRhoCorr(-999),
-  SumPtCharged(-999),
-  SumPtNeutral(-999),
-  SumPtChargedPU(-999),
-  SumPt(-999),
   fFactory(0),
   fArray(0)
 {
-  int i;
   Edges[0] = 0.0;
   Edges[1] = 0.0;
@@ -163 +150 @@
   Tau[3] = 0.0;
   Tau[4] = 0.0;
-  for(i = 0; i < 5; ++i)
-  {
-    TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
-    PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
-    SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
-  }
 }
 
@@ -244 +225 @@
   object.IsConstituent = IsConstituent;
   object.BTag = BTag;
-  object.BTagAlgo = BTagAlgo;
-  object.BTagDefault = BTagDefault;
-  object.BTagPhysics = BTagPhysics;
-  object.BTagNearest2 = BTagNearest2;
-  object.BTagNearest3 = BTagNearest3;
-  object.BTagHeaviest = BTagHeaviest;
-  object.BTagHighestPt = BTagHighestPt;
-  object.FlavorAlgo = FlavorAlgo;
-  object.FlavorDefault = FlavorDefault;
-  object.FlavorPhysics = FlavorPhysics;
-  object.FlavorNearest2 = FlavorNearest2;
-  object.FlavorNearest3 = FlavorNearest3;
-  object.FlavorHeaviest = FlavorHeaviest;
-  object.FlavorHighestPt = FlavorHighestPt;
   object.TauTag = TauTag;
   object.Eem = Eem;
@@ -282 +249 @@
   object.MeanSqDeltaR = MeanSqDeltaR;
   object.PTD = PTD;
-  object.Ntimes = Ntimes; 
-  object.IsolationVar = IsolationVar;
-  object.IsolationVarRhoCorr = IsolationVarRhoCorr;
-  object.SumPtCharged = SumPtCharged;
-  object.SumPtNeutral = SumPtNeutral;
-  object.SumPtChargedPU = SumPtChargedPU;
-  object.SumPt = SumPt;
-
   object.FracPt[0] = FracPt[0];
   object.FracPt[1] = FracPt[1];
@@ -300 +259 @@
   object.Tau[3] = Tau[3];
   object.Tau[4] = Tau[4];
-  
-  object.TrimmedP4[0] = TrimmedP4[0];
-  object.TrimmedP4[1] = TrimmedP4[1];
-  object.TrimmedP4[2] = TrimmedP4[2];
-  object.TrimmedP4[3] = TrimmedP4[3];
-  object.TrimmedP4[4] = TrimmedP4[4];
-  object.PrunedP4[0] = PrunedP4[0];
-  object.PrunedP4[1] = PrunedP4[1];
-  object.PrunedP4[2] = PrunedP4[2];
-  object.PrunedP4[3] = PrunedP4[3];
-  object.PrunedP4[4] = PrunedP4[4];
-  object.SoftDroppedP4[0] = SoftDroppedP4[0];
-  object.SoftDroppedP4[1] = SoftDroppedP4[1];
-  object.SoftDroppedP4[2] = SoftDroppedP4[2];
-  object.SoftDroppedP4[3] = SoftDroppedP4[3];
-  object.SoftDroppedP4[4] = SoftDroppedP4[4];
-
-  object.NSubJetsTrimmed = NSubJetsTrimmed;
-  object.NSubJetsPruned = NSubJetsPruned;
-  object.NSubJetsSoftDropped = NSubJetsSoftDropped;
 
   object.fFactory = fFactory;
   object.fArray = 0;
-
-  // Copy cluster timing info
-  copy(Ecal_E_t.begin(),Ecal_E_t.end(),back_inserter(object.Ecal_E_t));
 
   if(fArray && fArray->GetEntriesFast() > 0)
@@ -342 +278 @@
 void Candidate::Clear(Option_t* option)
 {
-  int i;
   SetUniqueID(0);
   ResetBit(kIsReferenced);
@@ -353 +288 @@
   IsConstituent = 0;
   BTag = 0;
-  BTagAlgo = 0;
-  BTagDefault = 0;
-  BTagPhysics = 0;
-  BTagNearest2 = 0;
-  BTagNearest3 = 0;
-  BTagHeaviest = 0;
-  BTagHighestPt = 0;
-  FlavorAlgo = 0;
-  FlavorDefault = 0;
-  FlavorPhysics = 0;
-  FlavorNearest2 = 0;
-  FlavorNearest3 = 0;
-  FlavorHeaviest = 0;
-  FlavorHighestPt = 0;
   TauTag = 0;
   Eem = 0.0;
@@ -390 +311 @@
   MeanSqDeltaR = 0.0;
   PTD = 0.0;
-  
-  Ntimes = 0;
-  Ecal_E_t.clear();
-  
-  IsolationVar = -999;
-  IsolationVarRhoCorr = -999;
-  SumPtCharged = -999;
-  SumPtNeutral = -999;
-  SumPtChargedPU = -999;
-  SumPt = -999;
-  
   FracPt[0] = 0.0;
   FracPt[1] = 0.0;
@@ -411 +321 @@
   Tau[3] = 0.0;
   Tau[4] = 0.0;
-  
-  for(i = 0; i < 5; ++i)
-  {
-    TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
-    PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 
-    SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0);
-  }
-
-  NSubJetsTrimmed = 0;
-  NSubJetsPruned = 0;
-  NSubJetsSoftDropped = 0;
-  
+
   fArray = 0;
 }

classes/DelphesClasses.h

@@ -84 +84 @@
 //---------------------------------------------------------------------------
 
-class LHEFWeight: public TObject
-{
-public:
-  Int_t ID; // weight ID
-  Float_t Weight; // weight value
-
-  ClassDef(LHEFWeight, 1)
-};
-
-//---------------------------------------------------------------------------
-
 class HepMCEvent: public Event
 {
@@ -242 +231 @@
   TRefArray Particles; // references to generated particles
 
-  // Isolation variables
-
-  Float_t IsolationVar;
-  Float_t IsolationVarRhoCorr;
-  Float_t SumPtCharged;
-  Float_t SumPtNeutral;
-  Float_t SumPtChargedPU;
-  Float_t SumPt;
-
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
@@ -277 +257 @@
   TRef Particle; // reference to generated particle
 
-  // Isolation variables
-
-  Float_t IsolationVar;
-  Float_t IsolationVarRhoCorr;
-  Float_t SumPtCharged;
-  Float_t SumPtNeutral;
-  Float_t SumPtChargedPU;
-  Float_t SumPt;
-
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
@@ -310 +281 @@
   TRef Particle; // reference to generated particle
 
-   // Isolation variables
-
-  Float_t IsolationVar;
-  Float_t IsolationVarRhoCorr;
-  Float_t SumPtCharged;
-  Float_t SumPtNeutral;
-  Float_t SumPtChargedPU;
-  Float_t SumPt;
-
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
@@ -345 +307 @@
 
   UInt_t BTag; // 0 or 1 for a jet that has been tagged as containing a heavy quark
-
-  UInt_t BTagAlgo;
-  UInt_t BTagDefault;
-  UInt_t BTagPhysics;
-  UInt_t BTagNearest2;
-  UInt_t BTagNearest3;
-  UInt_t BTagHeaviest;
-  UInt_t BTagHighestPt;
-
-  UInt_t FlavorAlgo;
-  UInt_t FlavorDefault;
-  UInt_t FlavorPhysics;
-  UInt_t FlavorNearest2;
-  UInt_t FlavorNearest3;
-  UInt_t FlavorHeaviest;
-  UInt_t FlavorHighestPt;
-
   UInt_t TauTag; // 0 or 1 for a jet that has been tagged as a tau
 
@@ -368 +313 @@
   Float_t EhadOverEem; // ratio of the hadronic versus electromagnetic energy deposited in the calorimeter
 
-  Int_t NCharged; // number of charged constituents
-  Int_t NNeutrals; // number of neutral constituents
-  Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
-  Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents)
-  Float_t MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
-  Float_t PTD; // average pt between constituent and jet weighted by pt of constituent
-  Float_t FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents)
-
-  Float_t Tau[5]; // N-subjettiness
-
-  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
-  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
-  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
-
-  Int_t NSubJetsTrimmed; // number of subjets trimmed
-  Int_t NSubJetsPruned; // number of subjets pruned
-  Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
+  Int_t    NCharged; // number of charged constituents 
+  Int_t    NNeutrals; // number of neutral constituents 
+  Float_t  Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents) 
+  Float_t  BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents) 
+  Float_t  MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
+  Float_t  PTD; // average pt between constituent and jet weighted by pt of constituent
+  Float_t  FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents) 
+
+  Float_t Tau1; // 1-subjettiness
+  Float_t Tau2; // 2-subjettiness
+  Float_t Tau3; // 3-subjettiness
+  Float_t Tau4; // 4-subjettiness
+  Float_t Tau5; // 5-subjettiness
 
   TRefArray Constituents; // references to constituents
@@ -393 +334 @@
 
   TLorentzVector P4() const;
-  TLorentzVector Area;
-
-  ClassDef(Jet, 3)
+
+  ClassDef(Jet, 2)
 };
 
@@ -452 +392 @@
   Float_t E; // calorimeter tower energy
 
-  Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
-  Int_t   Ntimes; // number of hits contributing to time measurement
+  Float_t T; //particle arrival time of flight
 
   Float_t Eem; // calorimeter tower electromagnetic energy
@@ -513 +452 @@
 
   Int_t IsPU;
-  Int_t IsRecoPU;
-
   Int_t IsConstituent;
 
   UInt_t BTag;
-
-  UInt_t BTagAlgo;
-  UInt_t BTagDefault;
-  UInt_t BTagPhysics;
-  UInt_t BTagNearest2;
-  UInt_t BTagNearest3;
-  UInt_t BTagHeaviest;
-  UInt_t BTagHighestPt;
-
-  UInt_t FlavorAlgo;
-  UInt_t FlavorDefault;
-  UInt_t FlavorPhysics;
-  UInt_t FlavorNearest2;
-  UInt_t FlavorNearest3;
-  UInt_t FlavorHeaviest;
-  UInt_t FlavorHighestPt;
-
   UInt_t TauTag;
 
@@ -562 +482 @@
   Float_t  FracPt[5];
 
-  //Timing information
-
-  Int_t    Ntimes;
-  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t;
-
-  // Isolation variables
-
-  Float_t IsolationVar;
-  Float_t IsolationVarRhoCorr;
-  Float_t SumPtCharged;
-  Float_t SumPtNeutral;
-  Float_t SumPtChargedPU;
-  Float_t SumPt;
-
   // N-subjettiness variables
 
   Float_t Tau[5];
-
-  // Other Substructure variables
-
-  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
-  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
-  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
-
-  Int_t NSubJetsTrimmed; // number of subjets trimmed
-  Int_t NSubJetsPruned; // number of subjets pruned
-  Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
-
 
   static CompBase *fgCompare; //!
@@ -609 +504 @@
   void SetFactory(DelphesFactory *factory) { fFactory = factory; }
 
-  ClassDef(Candidate, 3)
+  ClassDef(Candidate, 2)
 };
 

classes/DelphesLHEFReader.cc

@@ -82 +82 @@
   fEventCounter = -1;
   fParticleCounter = -1;
-  fWeightList.clear();
+  fRwgtList.clear();
 }
 
@@ -99 +99 @@
   TObjArray *partonOutputArray)
 {
-  int rc, id;
+  int rc;
   char *pch;
   double weight;
@@ -158 +158 @@
   else if(strstr(fBuffer, "<wgt"))
   {
-    pch = strpbrk(fBuffer, "\"'");
+    pch = strstr(fBuffer, ">");
     if(!pch)
     {
@@ -165 +165 @@
     }
 
-    DelphesStream idStream(pch + 1);
-    rc = idStream.ReadInt(id);
-
-    pch = strchr(fBuffer, '>');
-    if(!pch)
+    DelphesStream bufferStream(pch + 1);
+    rc = bufferStream.ReadDbl(weight);
+
+    if(!rc)
     {
       cerr << "** ERROR: " << "invalid weight format" << endl;
@@ -175 +174 @@
     }
 
-    DelphesStream weightStream(pch + 1);
-    rc = weightStream.ReadDbl(weight);
-
-    if(!rc)
-    {
-      cerr << "** ERROR: " << "invalid weight format" << endl;
-      return kFALSE;
-    }
-
-    fWeightList.push_back(make_pair(id, weight));
+    fRwgtList.push_back(weight);
   }
   else if(strstr(fBuffer, "</event>"))
@@ -216 +206 @@
 //---------------------------------------------------------------------------
 
-void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch)
-{
-  LHEFWeight *element;
-  vector< pair< int, double > >::const_iterator itWeightList;
-
-  for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList)
-  {
-    element = static_cast<LHEFWeight *>(branch->NewEntry());
-
-    element->ID = itWeightList->first;
-    element->Weight = itWeightList->second;
+void DelphesLHEFReader::AnalyzeRwgt(ExRootTreeBranch *branch)
+{
+  Weight *element;
+  vector<double>::const_iterator itRwgtList;
+
+  for(itRwgtList = fRwgtList.begin(); itRwgtList != fRwgtList.end(); ++itRwgtList)
+  {
+    element = static_cast<Weight *>(branch->NewEntry());
+
+    element->Weight = *itRwgtList;
   }
 }

classes/DelphesLHEFReader.h

@@ -31 +31 @@
 
 #include <vector>
-#include <utility>
 
 class TObjArray;
@@ -59 +58 @@
     TStopwatch *readStopWatch, TStopwatch *procStopWatch);
 
-  void AnalyzeWeight(ExRootTreeBranch *branch);
+  void AnalyzeRwgt(ExRootTreeBranch *branch);
 
 private:
@@ -84 +83 @@
   double fPx, fPy, fPz, fE, fMass;
   
-  std::vector< std::pair< int, double > > fWeightList;
+  std::vector<double> fRwgtList;
 };
 

doc/genMakefile.tcl

@@ -283 +283 @@
 dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h}
 
-sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} {external/PUPPI/*.cc}
+sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc}
 
 sourceDeps {FASTJET} {modules/FastJet*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/fastjet/contribs/*/*.cc} 

modules/Calorimeter.cc

@@ -150 +150 @@
 */
 
-  // read min E value for timing measurement in ECAL
-  fTimingEMin = GetDouble("TimingEMin",4.);
-  // For timing
-  // So far this flag needs to be false
-  // Curved extrapolation not supported
-  fElectronsFromTrack = false;
-
-  
   // read min E value for towers to be saved
   fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
@@ -364 +356 @@
       fTrackHCalEnergy = 0.0;
 
+      fTowerECalTime = 0.0;
+      fTowerHCalTime = 0.0;
+
+      fTrackECalTime = 0.0;
+      fTrackHCalTime = 0.0;
+
+      fTowerECalTimeWeight = 0.0;
+      fTowerHCalTimeWeight = 0.0;
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -379 +380 @@
       position = track->Position;
 
+
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -385 +387 @@
       fTrackHCalEnergy += hcalEnergy;
 
-      bool dbg_scz = false;
-      if (dbg_scz) {
-        cout << "   Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T() 
-             << endl;
-        Candidate *prt = static_cast<Candidate*>(track->GetCandidates()->Last());
-        const TLorentzVector &ini = prt->Position;
-
-        cout << "                and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T();
-
-      }
-      
-      if (ecalEnergy > fTimingEMin && fTower) {
-        if (fElectronsFromTrack) {
-          //      cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 
-          //           << " PID=" << track->PID << endl;
-          fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T()));
-        } else {
-          //      cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 
-          //           << " PID=" << track->PID << endl;
-        }
-      }
-
+      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
+      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -412 +397 @@
       continue;
     }
-  
+
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
@@ -427 +412 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    if (ecalEnergy > fTimingEMin && fTower) {
-      if (abs(particle->PID) != 11 || !fElectronsFromTrack) {
-        //      cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 
-        //   << " PID=" << particle->PID << endl;
-        fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T()));
-      } else {
-        
-        // N.B. Only charged particle set to leave ecal energy is the electrons
-        //      cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 
-        //           << " PID=" << particle->PID << endl;
-        
-      }
-    }
+    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
+    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
+
 
     fTower->AddCandidate(particle);
@@ -456 +434 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
- 
+  Double_t ecalTime, hcalTime, time;
+
   if(!fTower) return;
 
@@ -465 +444 @@
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
 
+  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
+  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
+
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -472 +454 @@
 
   energy = ecalEnergy + hcalEnergy;
-    
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
+
   if(fSmearTowerCenter)
   {
@@ -486 +469 @@
   pt = energy / TMath::CosH(eta);
 
-  // Time calculation for tower
-  fTower->Ntimes = 0;
-  Float_t tow_sumT = 0;
-  Float_t tow_sumW = 0;
-  
-  for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++) 
-  {
-    Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first);
-    tow_sumT += w*fTower->Ecal_E_t[i].second;
-    tow_sumW += w;
-    fTower->Ntimes++;
-  }
-  
-  if (tow_sumW > 0.) {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW);
-  } else {
-    fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.);
-  }
-
-
+  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;

modules/Calorimeter.h

@@ -60 +60 @@
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
 
-  Double_t fTimingEMin;
-  Bool_t fElectronsFromTrack; // for timing
-  
+  Double_t fTowerECalTime, fTowerHCalTime;
+  Double_t fTrackECalTime, fTrackHCalTime;
+
+  Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
+  Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
+
   Int_t fTowerTrackHits, fTowerPhotonHits;
 

modules/FastJetFinder.cc

@@ -65 +65 @@
 #include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
 #include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
-
-#include "fastjet/tools/Filter.hh"
-#include "fastjet/tools/Pruner.hh"
-#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
 
 using namespace std;
@@ -124 +120 @@
   fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
   fN = GetInt("N", 2);                  // used only if Njettiness is used as jet clustering algo (case 8)
-
-  //-- Trimming parameters --
-  
-  fComputeTrimming = GetBool("ComputeTrimming", false);
-  fRTrim = GetDouble("RTrim", 0.2);
-  fPtFracTrim = GetDouble("PtFracTrim", 0.05);
-  
-
-  //-- Pruning parameters --
-  
-  fComputePruning = GetBool("ComputePruning", false);
-  fZcutPrun = GetDouble("ZcutPrun", 0.1);
-  fRcutPrun = GetDouble("RcutPrun", 0.5);
-  fRPrun = GetDouble("RPrun", 0.8);
- 
-  //-- SoftDrop parameters --
-  
-  fComputeSoftDrop     = GetBool("ComputeSoftDrop", false);
-  fBetaSoftDrop        = GetDouble("BetaSoftDrop", 0.0);
-  fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
-  fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
-  
 
   // ---  Jet Area Parameters ---
@@ -286 +260 @@
   PseudoJet jet, area;
   ClusterSequence *sequence;
-  vector< PseudoJet > inputList, outputList, subjets;
+  vector< PseudoJet > inputList, outputList;
   vector< PseudoJet >::iterator itInputList, itOutputList;
   vector< TEstimatorStruct >::iterator itEstimators;
@@ -378 +352 @@
     candidate->DeltaEta = detaMax;
     candidate->DeltaPhi = dphiMax;
-       
-    //------------------------------------
-    // Trimming
-    //------------------------------------
-
-     if(fComputeTrimming)
-    {
-
-      fastjet::Filter    trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
-      fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
-      
-      trimmed_jet = join(trimmed_jet.constituents());
-     
-      candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
-        
-      // four hardest subjets 
-      subjets.clear();
-      subjets = trimmed_jet.pieces();
-      subjets = sorted_by_pt(subjets);
-      
-      candidate->NSubJetsTrimmed = subjets.size();
-
-      for (size_t i = 0; i < subjets.size() and i < 4; i++){
-        if(subjets.at(i).pt() < 0) continue ; 
-        candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
-      }
-    }
-    
-    
-    //------------------------------------
-    // Pruning
-    //------------------------------------
-    
-    
-    if(fComputePruning)
-    {
-
-      fastjet::Pruner    pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
-      fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
-
-      candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
-         
-      // four hardest subjet 
-      subjets.clear();
-      subjets = pruned_jet.pieces();
-      subjets = sorted_by_pt(subjets);
-      
-      candidate->NSubJetsPruned = subjets.size();
-
-      for (size_t i = 0; i < subjets.size() and i < 4; i++){
-        if(subjets.at(i).pt() < 0) continue ; 
-        candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
-      }
-
-    } 
-     
-    //------------------------------------
-    // SoftDrop
-    //------------------------------------
-   
-    if(fComputeSoftDrop)
-    {
-    
-      contrib::SoftDrop  softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
-      fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
-      
-      candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
-        
-      // four hardest subjet 
-      
-      subjets.clear();
-      subjets    = softdrop_jet.pieces();
-      subjets    = sorted_by_pt(subjets);
-      candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
-
-      for (size_t i = 0; i < subjets.size()  and i < 4; i++){
-        if(subjets.at(i).pt() < 0) continue ; 
-        candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
-      }
-    }
-  
+
     // --- compute N-subjettiness with N = 1,2,3,4,5 ----
 

modules/FastJetFinder.h

@@ -83 +83 @@
   Int_t fN ;
 
-  //-- Trimming parameters --
-  
-  Bool_t fComputeTrimming;
-  Double_t fRTrim;
-  Double_t fPtFracTrim;
-  
-  //-- Pruning parameters --
-
-  Bool_t fComputePruning;
-  Double_t fZcutPrun;
-  Double_t fRcutPrun;
-  Double_t fRPrun;
-
-  //-- SoftDrop parameters --
-
-  Bool_t fComputeSoftDrop;
-  Double_t fBetaSoftDrop;
-  Double_t fSymmetryCutSoftDrop;
-  Double_t fR0SoftDrop;
-
   // --- FastJet Area method --------
 

modules/Isolation.cc

@@ -25 +25 @@
  *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
  *
- *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
+ *  \author P. Demin - UCL, Louvain-la-Neuve
  *
  */
@@ -109 +109 @@
   fUsePTSum = GetBool("UsePTSum", false);
 
-  fVetoLeptons = GetBool("VetoLeptons", true);  
-  
   fClassifier->fPTMin = GetDouble("PTMin", 0.5);
-
 
   // import input array(s)
@@ -156 +153 @@
   Candidate *candidate, *isolation, *object;
   TObjArray *isolationArray;
-  Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
+  Double_t sum, ratio;
   Int_t counter;
   Double_t eta = 0.0;
@@ -183 +180 @@
 
     // loop over all input tracks
-    
-    sumNeutral = 0.0;
-    sumCharged = 0.0;
-    sumChargedPU = 0.0;
-    sumAllParticles = 0.0;
-   
+    sum = 0.0;
     counter = 0;
     itIsolationArray.Reset();
-    
     while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
     {
@@ -197 +188 @@
 
       if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
-         candidate->GetUniqueID() != isolation->GetUniqueID() &&
-         ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) )
+         candidate->GetUniqueID() != isolation->GetUniqueID())
       {
-     
-        sumAllParticles += isolationMomentum.Pt();
-        if(isolation->Charge !=0) 
-        { 
-          sumCharged += isolationMomentum.Pt();
-          if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt();
-        } 
-  
-        else sumNeutral += isolationMomentum.Pt();
-  
+        sum += isolationMomentum.Pt();
         ++counter;
       }
@@ -228 +209 @@
     }
 
-        
-     // correct sum for pile-up contamination
-    sumDBeta = sumCharged + TMath::Max(sumNeutral-0.5*sumChargedPU,0.0);
-    sumRhoCorr = sumCharged + TMath::Max(sumNeutral-TMath::Max(rho,0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(),0.0);
-    ratioDBeta = sumDBeta/candidateMomentum.Pt();
-    ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
-    
-    candidate->IsolationVar = ratioDBeta;
-    candidate->IsolationVarRhoCorr = ratioRhoCorr;
-    candidate->SumPtCharged = sumCharged;
-    candidate->SumPtNeutral = sumNeutral;
-    candidate->SumPtChargedPU = sumChargedPU;
-    candidate->SumPt = sumAllParticles;
-
-    if((fUsePTSum && sumDBeta > fPTSumMax) || ratioDBeta > fPTRatioMax) continue;
+    // correct sum for pile-up contamination
+    sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
+
+    ratio = sum/candidateMomentum.Pt();
+    if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
+
     fOutputArray->Add(candidate);
   }

modules/Isolation.h

@@ -59 +59 @@
   Bool_t fUsePTSum;
 
-  Bool_t fVetoLeptons;  
- 
   IsolationClassifier *fClassifier; //!
 

modules/ModulesLinkDef.h

@@ -42 +42 @@
 #include "modules/UniqueObjectFinder.h"
 #include "modules/TrackCountingBTagging.h"
-#include "modules/BTaggingCMS.h"
 #include "modules/BTagging.h"
 #include "modules/TauTagging.h"
@@ -59 +58 @@
 #include "modules/Weighter.h"
 #include "modules/Hector.h"
-#include "modules/RunPUPPI.h"
-#include "modules/JetFlavorAssociation.h"
 #include "modules/ExampleModule.h"
 
@@ -85 +82 @@
 #pragma link C++ class UniqueObjectFinder+;
 #pragma link C++ class TrackCountingBTagging+;
-#pragma link C++ class BTaggingCMS+;
 #pragma link C++ class BTagging+;
 #pragma link C++ class TauTagging+;
@@ -102 +98 @@
 #pragma link C++ class Weighter+;
 #pragma link C++ class Hector+;
-#pragma link C++ class RunPUPPI+;
-#pragma link C++ class JetFlavorAssociation+;
 #pragma link C++ class ExampleModule+;
 

modules/PileUpJetID.cc

@@ +1 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
 
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables
- *
- *  \author S. Zenz
+ *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *
+ *  \author S. Zenz, December 2013
  *
  */
@@ -23 +41 @@
 #include "TRandom3.h"
 #include "TObjArray.h"
-//#include "TDatabasePDG.h"
+#include "TDatabasePDG.h"
 #include "TLorentzVector.h"
 
@@ -36 +54 @@
 
 PileUpJetID::PileUpJetID() :
-  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
+  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
 {
 
@@ -52 +70 @@
 void PileUpJetID::Init()
 {
-
   fJetPTMin = GetDouble("JetPTMin", 20.0);
   fParameterR = GetDouble("ParameterR", 0.5);
   fUseConstituents = GetInt("UseConstituents", 0);
 
-
-  /*
-  Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
-  Double_t fBetaMinBarrel; // |eta| < 2.5
-  Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
-  Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
-  Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
-  */
-
-  fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
-  fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
-  fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
-  fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
-  fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
-  fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
-  fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
-
-
-
-  cout << "  set MeanSqDeltaRMaxBarrel " << fMeanSqDeltaRMaxBarrel << endl;
-  cout << "  set BetaMinBarrel " << fBetaMinBarrel << endl;
-  cout << "  set MeanSqDeltaRMaxEndcap " << fMeanSqDeltaRMaxEndcap << endl;
-  cout << "  set BetaMinEndcap " << fBetaMinEndcap << endl;
-  cout << "  set MeanSqDeltaRMaxForward " << fMeanSqDeltaRMaxForward << endl;
-
-
-
   fAverageEachTower = false; // for timing
 
@@ -91 +81 @@
   fItJetInputArray = fJetInputArray->MakeIterator();
 
-
-  //  cout << "BeforE SCZ additions in init" << endl;
-  //  cout << GetString("TrackInputArray", "ParticlePropagator/tracks") << endl;
-  //  cout << GetString("EFlowTrackInputArray", "ParticlePropagator/tracks") << endl;
-
-  fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
   fItTrackInputArray = fTrackInputArray->MakeIterator();
 
-  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
+  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
   fItNeutralInputArray = fNeutralInputArray->MakeIterator();
 
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+
+  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
 
   // create output array(s)
 
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
-
-  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
-
-
-  //  cout << " end of INIT " << endl;
-
 }
 
@@ -118 +101 @@
 void PileUpJetID::Finish()
 {
-  //  cout << "In finish" << endl;
-
   if(fItJetInputArray) delete fItJetInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
   if(fItNeutralInputArray) delete fItNeutralInputArray;
-
+  if(fItVertexInputArray) delete fItVertexInputArray;
 }
 
@@ -130 +111 @@
 void PileUpJetID::Process()
 {
-  //  cout << "start of process" << endl;
-
   Candidate *candidate, *constituent;
   TLorentzVector momentum, area;
-
-  //  cout << "BeforE SCZ additions in process" << endl;
-
-  // SCZ
-  Candidate *trk;
+  Int_t i, nc, nn;
+  Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
+  Double_t dr, pt, pt_ann[5];
+  Double_t zvtx = 0.0;
+
+  Candidate *track;
+
+  // find z position of primary vertex
+
+  fItVertexInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
+  {
+    if(!candidate->IsPU)
+    {
+      zvtx = candidate->Position.Z();
+      break;
+    }
+  }
 
   // loop over all input candidates
@@ -147 +139 @@
     area = candidate->Area;
 
-    float sumT0 = 0.;
-    float sumT1 = 0.;
-    float sumT10 = 0.;
-    float sumT20 = 0.;
-    float sumT30 = 0.;
-    float sumT40 = 0.;
-    float sumWeightsForT = 0.;
-    candidate->Ntimes = 0;
-
-    float sumpt = 0.;
-    float sumptch = 0.;
-    float sumptchpv = 0.;
-    float sumptchpu = 0.;
-    float sumdrsqptsq = 0.;
-    float sumptsq = 0.;
-    int nc = 0;
-    int nn = 0;
-    float pt_ann[5];
-
-    for (int i = 0 ; i < 5 ; i++) {
-      pt_ann[i] = 0.;
-    }
-
-    if (fUseConstituents) {
+    sumpt = 0.0;
+    sumptch = 0.0;
+    sumptchpv = 0.0;
+    sumptchpu = 0.0;
+    sumdrsqptsq = 0.0;
+    sumptsq = 0.0;
+    nc = 0;
+    nn = 0;
+
+    for(i = 0; i < 5; ++i)
+    {
+      pt_ann[i] = 0.0;
+    }
+
+    if(fUseConstituents)
+    {
       TIter itConstituents(candidate->GetCandidates());
-      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
-        float pt = constituent->Momentum.Pt();
-        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
-        //      cout << " There exists a constituent with dr=" << dr << endl;
-        sumpt += pt;
-        sumdrsqptsq += dr*dr*pt*pt;
-        sumptsq += pt*pt;
-        if (constituent->Charge == 0) {
-          nn++;
-        } else {
-          if (constituent->IsRecoPU) {
-            sumptchpu += pt;
-          } else {
-            sumptchpv += pt;
-          }
-          sumptch += pt;
-          nc++;
-        }
-        for (int i = 0 ; i < 5 ; i++) {
-          if (dr > 0.1*i && dr < 0.1*(i+1)) {
-            pt_ann[i] += pt;
-          }
-        }
-        float tow_sumT = 0;
-        float tow_sumW = 0;
-        for (int i = 0 ; i < constituent->Ecal_E_t.size() ; i++) {
-          float w = TMath::Sqrt(constituent->Ecal_E_t[i].first);
-          if (fAverageEachTower) {
-            tow_sumT += w*constituent->Ecal_E_t[i].second;
-            tow_sumW += w;
-          } else {
-            sumT0 += w*constituent->Ecal_E_t[i].second;
-            sumT1 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.001);
-            sumT10 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.010);
-            sumT20 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.020);
-            sumT30 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.030);
-            sumT40 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.040);
-            sumWeightsForT += w;
-            candidate->Ntimes++;
-          }
-        }
-        if (fAverageEachTower && tow_sumW > 0.) {
-          sumT0 += tow_sumT;
-          sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
-          sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
-          sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
-          sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
-          sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
-          sumWeightsForT += tow_sumW;
-          candidate->Ntimes++;
-        }
-      }
-    } else {
+      while((constituent = static_cast<Candidate*>(itConstituents.Next())))
+      {
+        pt = constituent->Momentum.Pt();
+        dr = candidate->Momentum.DeltaR(constituent->Momentum);
+        sumpt += pt;
+        sumdrsqptsq += dr*dr*pt*pt;
+        sumptsq += pt*pt;
+        if(constituent->Charge == 0)
+        {
+          // neutrals
+          ++nn;
+        }
+        else
+        {
+          // charged
+          if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
+          {
+            sumptchpu += pt;
+          }
+          else
+          {
+            sumptchpv += pt;
+          }
+          sumptch += pt;
+          ++nc;
+        }
+        for(i = 0; i < 5; ++i)
+        {
+          if(dr > 0.1*i && dr < 0.1*(i + 1))
+          {
+            pt_ann[i] += pt;
+          }
+        }
+      }
+    }
+    else
+    {
       // Not using constituents, using dr
       fItTrackInputArray->Reset();
-      while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
-        if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
-          float pt = trk->Momentum.Pt();
-          sumpt += pt;
-          sumptch += pt;
-          if (trk->IsRecoPU) {
-            sumptchpu += pt;
-          } else {
-            sumptchpv += pt;
-          }
-          float dr = candidate->Momentum.DeltaR(trk->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nc++;
-          for (int i = 0 ; i < 5 ; i++) {
-            if (dr > 0.1*i && dr < 0.1*(i+1)) {
-              pt_ann[i] += pt;
-            }
-          }
-        }
-      }
+      while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+      {
+        if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
+        {
+          pt = track->Momentum.Pt();
+          sumpt += pt;
+          sumptch += pt;
+          if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
+          {
+            sumptchpu += pt;
+          }
+          else
+          {
+            sumptchpv += pt;
+          }
+          dr = candidate->Momentum.DeltaR(track->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nc++;
+          for(i = 0; i < 5; ++i)
+          {
+            if(dr > 0.1*i && dr < 0.1*(i + 1))
+            {
+              pt_ann[i] += pt;
+            }
+          }
+        }
+      }
+
       fItNeutralInputArray->Reset();
-      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
-        if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
-          float pt = constituent->Momentum.Pt();
-          sumpt += pt;
-          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nn++;
-          for (int i = 0 ; i < 5 ; i++) {
-            if (dr > 0.1*i && dr < 0.1*(i+1)) {
-            pt_ann[i] += pt;
-            }
-          }
-        }
-      }
-    }
-
-    if (sumptch > 0.) {
-      candidate->Beta = sumptchpv/sumptch;
-      candidate->BetaStar = sumptchpu/sumptch;
-    } else {
-      candidate->Beta = -999.;
-      candidate->BetaStar = -999.;
-    }
-    if (sumptsq > 0.) {
+      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
+      {
+        if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
+        {
+          pt = constituent->Momentum.Pt();
+          sumpt += pt;
+          dr = candidate->Momentum.DeltaR(constituent->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nn++;
+          for(i = 0; i < 5; ++i)
+          {
+            if(dr > 0.1*i && dr < 0.1*(i + 1))
+            {
+              pt_ann[i] += pt;
+            }
+          }
+        }
+      }
+    }
+
+    if(sumptch > 0.0)
+    {
+      candidate->Beta = sumptchpu/sumptch;
+      candidate->BetaStar = sumptchpv/sumptch;
+    }
+    else
+    {
+      candidate->Beta = -999.0;
+      candidate->BetaStar = -999.0;
+    }
+    if(sumptsq > 0.0)
+    {
       candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
-    } else {
-      candidate->MeanSqDeltaR = -999.;
+    }
+    else
+    {
+      candidate->MeanSqDeltaR = -999.0;
     }
     candidate->NCharged = nc;
     candidate->NNeutrals = nn;
-    if (sumpt > 0.) {
+    if(sumpt > 0.0)
+    {
       candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
-      for (int i = 0 ; i < 5 ; i++) {
+      for(i = 0; i < 5; ++i)
+      {
         candidate->FracPt[i] = pt_ann[i]/sumpt;
       }
-    } else {
-      candidate->PTD = -999.;
-      for (int i = 0 ; i < 5 ; i++) {
-        candidate->FracPt[i] = -999.;
+    }
+    else
+    {
+      candidate->PTD = -999.0;
+      for(i = 0; i < 5; ++i)
+      {
+        candidate->FracPt[i] = -999.0;
       }
     }
 
     fOutputArray->Add(candidate);
-
-    // New stuff
-    /*
-    fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
-    fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
-    fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
-    fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
-    fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
-    */
-
-    bool passId = false;
-    if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
-      if (fabs(candidate->Momentum.Eta())<1.5) {
-        passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
-      } else if (fabs(candidate->Momentum.Eta())<4.0) {
-        passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
-      } else {
-        passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
-      }
-    }
-
-    //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 
-    //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
-
-    if (passId) {
-      if (fUseConstituents) {
-        TIter itConstituents(candidate->GetCandidates());
-        while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
-          if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
-            fNeutralsInPassingJets->Add(constituent);
-            //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
-          }
-        }
-      } else { // use DeltaR
-        fItNeutralInputArray->Reset();
-        while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
-          if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
-            fNeutralsInPassingJets->Add(constituent);
-            //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
-          }
-        }
-      }
-    }
-
-
   }
 }
 
 //------------------------------------------------------------------------------
+

modules/PileUpJetID.h

@@ +1 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
 #ifndef PileUpJetID_h
 #define PileUpJetID_h
@@ -4 +22 @@
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables
+ *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
  *
- *  \author S. Zenz
+ *  \author S. Zenz, December 2013
  *
  */
@@ -16 +34 @@
 
 class TObjArray;
-class DelphesFormula;
 
 class PileUpJetID: public DelphesModule
@@ -34 +51 @@
   Double_t fParameterR;
 
-  Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
-  Double_t fBetaMinBarrel; // |eta| < 2.5
-  Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
-  Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
-  Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
-
-  Double_t fNeutralPTMin;
-  Double_t fJetPTMinForNeutrals;
-
-  /*
-JAY
----
-
-|Eta|<1.5
-
-meanSqDeltaR betaStar SigEff BgdEff
-0.13 0.92 96% 8%
-0.13 0.95 97% 16%
-0.13 0.97 98% 27%
-
-|Eta|>1.5
-
-meanSqDeltaR betaStar SigEff BgdEff
-0.14 0.91 95% 15%
-0.14 0.94 97% 19%
-0.14 0.97 98% 29%
-
-BRYAN
------
-
-Barrel (MeanSqDR, Beta, sig eff, bg eff): 
-0.10, 0.08, 90%, 8%
-0.11, 0.12, 90%, 6%
-0.13, 0.16, 89%, 5%
-
-Endcap (MeanSqDR, Beta, sig eff, bg eff):
-0.07, 0.06, 89%, 4%
-0.08, 0.08, 92%, 6%
-0.09, 0.08, 95%, 10%
-0.10, 0.08, 97%, 13%
-
-SETH GUESSES FOR |eta| > 4.0
-----------------------------
-
-MeanSqDeltaR 
-0.07 
-0.10 
-0.14 
-0.2
-  */
-
   // If set to true, may have weird results for PFCHS
   // If set to false, uses everything within dR < fParameterR even if in other jets &c.
   // Results should be very similar for PF
-  Int_t fUseConstituents; 
+  Int_t fUseConstituents;
 
   Bool_t fAverageEachTower;
@@ -96 +62 @@
   const TObjArray *fJetInputArray; //!
 
-  const TObjArray *fTrackInputArray; // SCZ
-  const TObjArray *fNeutralInputArray; 
+  const TObjArray *fTrackInputArray; //!
+  const TObjArray *fNeutralInputArray; //!
 
-  TIterator *fItTrackInputArray; // SCZ
-  TIterator *fItNeutralInputArray; // SCZ
+  TIterator *fItTrackInputArray; //!
+  TIterator *fItNeutralInputArray; //!
 
   TObjArray *fOutputArray; //!
-  TObjArray *fNeutralsInPassingJets; // SCZ
 
+  TIterator *fItVertexInputArray; //!
+  const TObjArray *fVertexInputArray; //!
 
-  ClassDef(PileUpJetID, 2)
+  Double_t fZVertexResolution;
+
+  ClassDef(PileUpJetID, 1)
 };
 
 #endif
+

modules/PileUpMerger.cc

@@ -80 +80 @@
   fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
 
-  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
-  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
-  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
-  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
-
   // read vertex smearing formula
 
@@ -116 +111 @@
   TParticlePDG *pdgParticle;
   Int_t pid;
-  Float_t x, y, z, t, vx, vy;
+  Float_t x, y, z, t;
   Float_t px, py, pz, e;
   Double_t dz, dphi, dt;
-  Int_t numberOfEvents, event, numberOfParticles;
+  Int_t numberOfEvents, event;
   Long64_t allEntries, entry;
-  Candidate *candidate, *vertex;
+  Candidate *candidate, *vertexcandidate;
   DelphesFactory *factory;
 
@@ -128 +123 @@
   fItInputArray->Reset();
 
-  // --- Deal with primary vertex first  ------
+  // --- Deal with Primary vertex first  ------
 
   fFunction->GetRandom2(dz, dt);
@@ -134 +129 @@
   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
   dz *= 1.0E3; // necessary in order to make z in mm
-  vx = 0.0;
-  vy = 0.0;
-  numberOfParticles = fInputArray->GetEntriesFast();
+
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
-    vx += candidate->Position.X();
-    vy += candidate->Position.Y();
     z = candidate->Position.Z();
     t = candidate->Position.T();
@@ -148 +139 @@
   }
 
-  if(numberOfParticles > 0)
-  {
-    vx /= numberOfParticles;
-    vy /= numberOfParticles;
-  }
-
   factory = GetFactory();
 
-  vertex = factory->NewCandidate();
-  vertex->Position.SetXYZT(vx, vy, dz, dt);
-  fVertexOutputArray->Add(vertex);
+  vertexcandidate = factory->NewCandidate();
+  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+  fVertexOutputArray->Add(vertexcandidate);
 
   // --- Then with pile-up vertices  ------
@@ -196 +181 @@
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    vx = 0.0;
-    vy = 0.0;
-    numberOfParticles = 0;
+    vertexcandidate = factory->NewCandidate();
+    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+    vertexcandidate->IsPU = 1;
+
+    fVertexOutputArray->Add(vertexcandidate);
+
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     {
@@ -216 +204 @@
       candidate->Momentum.RotateZ(dphi);
 
-      x -= fInputBeamSpotX;
-      y -= fInputBeamSpotY;
       candidate->Position.SetXYZT(x, y, z + dz, t + dt);
       candidate->Position.RotateZ(dphi);
-      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
-
-      vx += candidate->Position.X();
-      vy += candidate->Position.Y();
-      ++numberOfParticles;
 
       fParticleOutputArray->Add(candidate);
     }
-
-    if(numberOfParticles > 0)
-    {
-      vx /= numberOfParticles;
-      vy /= numberOfParticles;
-    }
-
-    vertex = factory->NewCandidate();
-    vertex->Position.SetXYZT(vx, vy, dz, dt);
-    vertex->IsPU = 1;
-
-    fVertexOutputArray->Add(vertex);
   }
 }
 
 //------------------------------------------------------------------------------
+

modules/PileUpMerger.h

@@ -53 +53 @@
   Double_t fTVertexSpread;
 
-  Double_t fInputBeamSpotX;
-  Double_t fInputBeamSpotY;
-  Double_t fOutputBeamSpotX;
-  Double_t fOutputBeamSpotY;
-
   DelphesTF2 *fFunction; //!
 

modules/PileUpMergerPythia8.cc

@@ -29 +29 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesTF2.h"
+#include "classes/DelphesFormula.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -56 +56 @@
 
 PileUpMergerPythia8::PileUpMergerPythia8() :
-  fFunction(0), fPythia(0), fItInputArray(0)
+  fPythia(0), fItInputArray(0)
 {
-  fFunction = new DelphesTF2;
 }
 
@@ -65 +64 @@
 PileUpMergerPythia8::~PileUpMergerPythia8()
 {
-  delete fFunction;
 }
 
@@ -74 +72 @@
   const char *fileName;
 
-  fPileUpDistribution = GetInt("PileUpDistribution", 0);
-
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-
-  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
-  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
-
-  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
-  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
-  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
-  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
+  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
 
   fPTMin = GetDouble("PTMin", 0.0);
-
-  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
-  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
 
   fileName = GetString("ConfigFile", "MinBias.cmnd");
@@ -100 +86 @@
 
   // create output arrays
-  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
-  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
 }
 
@@ -118 +103 @@
   TParticlePDG *pdgParticle;
   Int_t pid, status;
-  Float_t x, y, z, t, vx, vy;
+  Float_t x, y, z, t;
   Float_t px, py, pz, e;
-  Double_t dz, dphi, dt;
-  Int_t numberOfEvents, event, numberOfParticles, i;
-  Candidate *candidate, *vertex;
+  Double_t dz, dphi;
+  Int_t poisson, event, i;
+  Candidate *candidate;
   DelphesFactory *factory;
 
-  const Double_t c_light = 2.99792458E8;
-
   fItInputArray->Reset();
-
-  // --- Deal with primary vertex first  ------
-
-  fFunction->GetRandom2(dz, dt);
-
-  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
-  dz *= 1.0E3; // necessary in order to make z in mm
-  vx = 0.0;
-  vy = 0.0;
-  numberOfParticles = fInputArray->GetEntriesFast();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
-    vx += candidate->Position.X();
-    vy += candidate->Position.Y();
-    z = candidate->Position.Z();
-    t = candidate->Position.T();
-    candidate->Position.SetZ(z + dz);
-    candidate->Position.SetT(t + dt);
-    fParticleOutputArray->Add(candidate);
-  }
-
-  if(numberOfParticles > 0)
-  {
-    vx /= numberOfParticles;
-    vy /= numberOfParticles;
+    fOutputArray->Add(candidate);
   }
 
   factory = GetFactory();
 
-  vertex = factory->NewCandidate();
-  vertex->Position.SetXYZT(vx, vy, dz, dt);
-  fVertexOutputArray->Add(vertex);
+  poisson = gRandom->Poisson(fMeanPileUp);
 
-  // --- Then with pile-up vertices  ------
-
-  switch(fPileUpDistribution)
-  {
-    case 0:
-      numberOfEvents = gRandom->Poisson(fMeanPileUp);
-      break;
-    case 1:
-      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
-      break;
-    default:
-      numberOfEvents = gRandom->Poisson(fMeanPileUp);
-      break;
-  }
-
-  for(event = 0; event < numberOfEvents; ++event)
+  for(event = 0; event < poisson; ++event)
   {
     while(!fPythia->next());
 
-   // --- Pile-up vertex smearing
-
-    fFunction->GetRandom2(dz, dt);
-
-    dt *= c_light*1.0E3; // necessary in order to make t in mm/c
-    dz *= 1.0E3; // necessary in order to make z in mm
-
+    dz = gRandom->Gaus(0.0, fZVertexSpread);
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    vx = 0.0;
-    vy = 0.0;
-    numberOfParticles = fPythia->event.size();
-    for(i = 1; i < numberOfParticles; ++i)
+    for(i = 1; i < fPythia->event.size(); ++i)
     {
       Pythia8::Particle &particle = fPythia->event[i];
@@ -208 +143 @@
       candidate->PID = pid;
 
-      candidate->Status = 1;
+      candidate->Status = status;
 
       pdgParticle = pdg->GetParticle(pid);
@@ -219 +154 @@
       candidate->Momentum.RotateZ(dphi);
 
-      x -= fInputBeamSpotX;
-      y -= fInputBeamSpotY;
-      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
+      candidate->Position.SetXYZT(x, y, z + dz, t);
       candidate->Position.RotateZ(dphi);
-      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
 
-      vx += candidate->Position.X();
-      vy += candidate->Position.Y();
-
-      fParticleOutputArray->Add(candidate);
+      fOutputArray->Add(candidate);
     }
-
-    if(numberOfParticles > 0)
-    {
-      vx /= numberOfParticles;
-      vy /= numberOfParticles;
-    }
-
-    vertex = factory->NewCandidate();
-    vertex->Position.SetXYZT(vx, vy, dz, dt);
-    vertex->IsPU = 1;
-
-    fVertexOutputArray->Add(vertex);
   }
 }
 
 //------------------------------------------------------------------------------
+

modules/PileUpMergerPythia8.h

@@ -31 +31 @@
 
 class TObjArray;
-class DelphesTF2;
 
 namespace Pythia8
@@ -51 +50 @@
 private:
 
-  Int_t fPileUpDistribution;
   Double_t fMeanPileUp;
-
   Double_t fZVertexSpread;
-  Double_t fTVertexSpread;
-
-  Double_t fInputBeamSpotX;
-  Double_t fInputBeamSpotY;
-  Double_t fOutputBeamSpotX;
-  Double_t fOutputBeamSpotY;
-
   Double_t fPTMin;
-
-  DelphesTF2 *fFunction; //!
 
   Pythia8::Pythia *fPythia; //!
@@ -72 +60 @@
   const TObjArray *fInputArray; //!
 
-  TObjArray *fParticleOutputArray; //!
-  TObjArray *fVertexOutputArray; //!
+  TObjArray *fOutputArray; //!
 
   ClassDef(PileUpMergerPythia8, 1)

modules/TrackPileUpSubtractor.cc

@@ -74 +74 @@
   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
 
-  fPTMin = GetDouble("PTMin", 0.);
   // import arrays with output from other modules
-   
+
   ExRootConfParam param = GetParam("InputArray");
   Long_t i, size;
@@ -148 +147 @@
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
       
-      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1;
-      else 
-      {
-         candidate->IsRecoPU = 0;
-         if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
-      }
+      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
+
+      array->Add(candidate);
     }
   }

modules/TrackPileUpSubtractor.h

@@ -50 +50 @@
   Double_t fZVertexResolution;
 
-  Double_t fPTMin; 
-
   std::map< TIterator *, TObjArray * > fInputMap; //!
 

modules/TreeWriter.cc

@@ -291 +291 @@
     entry->ZOuter = position.Z();
     entry->TOuter = position.T()*1.0E-3/c_light;
-
+    
     entry->Dxy = candidate->Dxy;
     entry->SDxy = candidate->SDxy ;
@@ -297 +297 @@
     entry->Yd = candidate->Yd;
     entry->Zd = candidate->Zd;
-
+   
     const TLorentzVector &momentum = candidate->Momentum;
 
@@ -362 +362 @@
 
     entry->T = position.T()*1.0E-3/c_light;
-    entry->Ntimes = candidate->Ntimes;
 
     FillParticles(candidate, &entry->Particles);
@@ -404 +403 @@
     entry->T = position.T()*1.0E-3/c_light;
 
-    // Isolation variables
-
-    entry->IsolationVar = candidate->IsolationVar;
-    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
-    entry->SumPtCharged = candidate->SumPtCharged ;
-    entry->SumPtNeutral = candidate->SumPtNeutral ;
-    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
-    entry->SumPt = candidate->SumPt ;
-
     entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
 
@@ -452 +442 @@
     entry->T = position.T()*1.0E-3/c_light;
 
-    // Isolation variables
-
-    entry->IsolationVar = candidate->IsolationVar;
-    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
-    entry->SumPtCharged = candidate->SumPtCharged ;
-    entry->SumPtNeutral = candidate->SumPtNeutral ;
-    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
-    entry->SumPt = candidate->SumPt ;
-
-
     entry->Charge = candidate->Charge;
 
@@ -508 +488 @@
     entry->T = position.T()*1.0E-3/c_light;
 
-    // Isolation variables
-
-    entry->IsolationVar = candidate->IsolationVar;
-    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
-    entry->SumPtCharged = candidate->SumPtCharged ;
-    entry->SumPtNeutral = candidate->SumPtNeutral ;
-    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
-    entry->SumPt = candidate->SumPt ;
-
     entry->Charge = candidate->Charge;
 
@@ -533 +504 @@
   Double_t ecalEnergy, hcalEnergy;
   const Double_t c_light = 2.99792458E8;
-  Int_t i;
 
   array->Sort();
@@ -562 +532 @@
     entry->Mass = momentum.M();
 
-    entry->Area = candidate->Area;
-
     entry->DeltaEta = candidate->DeltaEta;
     entry->DeltaPhi = candidate->DeltaPhi;
 
     entry->BTag = candidate->BTag;
-
-    entry->BTagAlgo = candidate->BTagAlgo;
-    entry->BTagDefault = candidate->BTagDefault;
-    entry->BTagPhysics = candidate->BTagPhysics;
-    entry->BTagNearest2 = candidate->BTagNearest2;
-    entry->BTagNearest3 = candidate->BTagNearest3;
-    entry->BTagHeaviest = candidate->BTagHeaviest;
-    entry->BTagHighestPt = candidate->BTagHighestPt;
-
-    entry->FlavorAlgo = candidate->FlavorAlgo;
-    entry->FlavorDefault = candidate->FlavorDefault;
-    entry->FlavorPhysics = candidate->FlavorPhysics;
-    entry->FlavorNearest2 = candidate->FlavorNearest2;
-    entry->FlavorNearest3 = candidate->FlavorNearest3;
-    entry->FlavorHeaviest = candidate->FlavorHeaviest;
-    entry->FlavorHighestPt = candidate->FlavorHighestPt;
-
     entry->TauTag = candidate->TauTag;
 
@@ -610 +561 @@
     entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
     entry->PTD = candidate->PTD;
-
-    //--- Sub-structure variables ----
-
-    entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
-    entry->NSubJetsPruned = candidate->NSubJetsPruned;
-    entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
-
-    for(i = 0; i < 5; i++)
-    {
-      entry->FracPt[i] = candidate -> FracPt[i];
-      entry->Tau[i] = candidate -> Tau[i];
-      entry->TrimmedP4[i] = candidate -> TrimmedP4[i];
-      entry->PrunedP4[i] = candidate -> PrunedP4[i];
-      entry->SoftDroppedP4[i] = candidate -> SoftDroppedP4[i];
-    }
-
+    entry->FracPt[0] = candidate->FracPt[0];
+    entry->FracPt[1] = candidate->FracPt[1];
+    entry->FracPt[2] = candidate->FracPt[2];
+    entry->FracPt[3] = candidate->FracPt[3];
+    entry->FracPt[4] = candidate->FracPt[4];
+   
+    //--- N-subjettiness variables ----
+    
+    entry->Tau1 = candidate->Tau[0];
+    entry->Tau2 = candidate->Tau[1];
+    entry->Tau3 = candidate->Tau[2];
+    entry->Tau4 = candidate->Tau[3];
+    entry->Tau5 = candidate->Tau[4];
+   
     FillParticles(candidate, &entry->Particles);
   }

readers/DelphesLHEF.cpp

@@ -63 +63 @@
   TStopwatch readStopWatch, procStopWatch;
   ExRootTreeWriter *treeWriter = 0;
-  ExRootTreeBranch *branchEvent = 0, *branchWeight = 0;
+  ExRootTreeBranch *branchEvent = 0, *branchRwgt = 0;
   ExRootConfReader *confReader = 0;
   Delphes *modularDelphes = 0;
@@ -103 +103 @@
 
     branchEvent = treeWriter->NewBranch("Event", LHEFEvent::Class());
-    branchWeight = treeWriter->NewBranch("Weight", Weight::Class());
+    branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
 
     confReader = new ExRootConfReader;
@@ -196 +196 @@
 
             reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch);
-            reader->AnalyzeWeight(branchWeight);
+            reader->AnalyzeRwgt(branchRwgt);
 
             treeWriter->Fill();

readers/DelphesPythia8.cpp

@@ -20 +20 @@
 #include <iostream>
 #include <sstream>
-#include <string>
 
 #include <signal.h>
@@ -39 +38 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesLHEFReader.h"
 
 #include "ExRootAnalysis/ExRootTreeWriter.h"
@@ -151 +149 @@
   char appName[] = "DelphesPythia8";
   stringstream message;
-  FILE *inputFile = 0;
   TFile *outputFile = 0;
   TStopwatch readStopWatch, procStopWatch;
   ExRootTreeWriter *treeWriter = 0;
   ExRootTreeBranch *branchEvent = 0;
-  ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0;
   ExRootConfReader *confReader = 0;
   Delphes *modularDelphes = 0;
   DelphesFactory *factory = 0;
   TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
-  TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0;
-  DelphesLHEFReader *reader = 0;
   Long64_t eventCounter, errorCounter;
   Long64_t numberOfEvents, timesAllowErrors;
@@ -211 +205 @@
     partonOutputArray = modularDelphes->ExportArray("partons");
 
-    // Initialize Pythia
+    modularDelphes->InitTask();
+
+    // Initialize pythia
     pythia = new Pythia8::Pythia;
 
@@ -225 +221 @@
     numberOfEvents = pythia->mode("Main:numberOfEvents");
     timesAllowErrors = pythia->mode("Main:timesAllowErrors");
-
-    inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
-    if(inputFile)
-    {
-      reader = new DelphesLHEFReader;
-      reader->SetInputFile(inputFile);
-
-      branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
-      branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
-
-      allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
-      stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
-      partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
-    }
-
-    modularDelphes->InitTask();
 
     pythia->init();
@@ -254 +234 @@
     for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
     {
-      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
-        stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
-
       if(!pythia->next())
       {
         // If failure because reached end of file then exit event loop
-        if(pythia->info.atEndOfFile())
+        if (pythia->info.atEndOfFile())
         {
           cerr << "Aborted since reached end of Les Houches Event File" << endl;
@@ -267 +244 @@
 
         // First few failures write off as "acceptable" errors, then quit
-        if(++errorCounter > timesAllowErrors)
-        {
-          cerr << "Event generation aborted prematurely, owing to error!" << endl;
-          break;
-        }
-
-        modularDelphes->Clear();
-        reader->Clear();
-        continue;
+        if (++errorCounter < timesAllowErrors) continue;
+        cerr << "Event generation aborted prematurely, owing to error!" << endl;
+        break;
       }
 
@@ -287 +258 @@
       procStopWatch.Stop();
 
-      if(reader)
-      {
-        reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch);
-        reader->AnalyzeWeight(branchWeightLHEF);
-      }
-
       treeWriter->Fill();
 
       treeWriter->Clear();
       modularDelphes->Clear();
-      if(reader) reader->Clear();
 
       readStopWatch.Start();
@@ -313 +277 @@
     cout << "** Exiting..." << endl;
 
-    delete reader;
     delete pythia;
     delete modularDelphes;