Index: /trunk/Makefile
===================================================================
--- /trunk/Makefile	(revision 1347)
+++ /trunk/Makefile	(revision 1348)
@@ -327,4 +327,5 @@
 	modules/JetPileUpSubtractor.h \
 	modules/TrackPileUpSubtractor.h \
+	modules/PileUpJetID.h \
 	modules/ConstituentFilter.h \
 	modules/StatusPidFilter.h \
@@ -511,4 +512,13 @@
 	modules/EnergyScale.$(SrcSuf) \
 	modules/EnergyScale.h \
+	classes/DelphesClasses.h \
+	classes/DelphesFactory.h \
+	classes/DelphesFormula.h \
+	external/ExRootAnalysis/ExRootResult.h \
+	external/ExRootAnalysis/ExRootFilter.h \
+	external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/PileUpJetID.$(ObjSuf): \
+	modules/PileUpJetID.$(SrcSuf) \
+	modules/PileUpJetID.h \
 	classes/DelphesClasses.h \
 	classes/DelphesFactory.h \
@@ -948,4 +958,5 @@
 	tmp/modules/ExampleModule.$(ObjSuf) \
 	tmp/modules/EnergyScale.$(ObjSuf) \
+	tmp/modules/PileUpJetID.$(ObjSuf) \
 	tmp/modules/Cloner.$(ObjSuf) \
 	tmp/modules/TreeWriter.$(ObjSuf) \
@@ -1296,4 +1307,8 @@
 	@touch $@
 
+modules/PileUpJetID.h: \
+	classes/DelphesModule.h
+	@touch $@
+
 external/fastjet/version.hh: \
 	external/fastjet/config.h
Index: /trunk/classes/DelphesClasses.cc
===================================================================
--- /trunk/classes/DelphesClasses.cc	(revision 1347)
+++ /trunk/classes/DelphesClasses.cc	(revision 1348)
@@ -112,5 +112,11 @@
   Area(0.0, 0.0, 0.0, 0.0),
   fFactory(0),
-  fArray(0)
+  fArray(0),
+  NCharged(0),
+  NNeutrals(0),
+  Beta(0),
+  BetaStar(0),
+  MeanSqDeltaR(0),
+  PTD(0)
 {
   Edges[0] = 0.0;
@@ -118,4 +124,9 @@
   Edges[2] = 0.0;
   Edges[3] = 0.0;
+  FracPt[0] = 0.0;
+  FracPt[1] = 0.0;
+  FracPt[2] = 0.0;
+  FracPt[3] = 0.0;
+  FracPt[4] = 0.0;
 }
 
@@ -206,4 +217,16 @@
   object.Area = Area;
 
+  object.NCharged = NCharged;
+  object.NNeutrals = NNeutrals;
+  object.Beta = Beta;
+  object.BetaStar = BetaStar;
+  object.MeanSqDeltaR = MeanSqDeltaR;
+  object.PTD = PTD;
+  object.FracPt[0] = FracPt[0];
+  object.FracPt[1] = FracPt[1];
+  object.FracPt[2] = FracPt[2];
+  object.FracPt[3] = FracPt[3];
+  object.FracPt[4] = FracPt[4];
+
   object.fFactory = fFactory;
   object.fArray = 0;
@@ -246,4 +269,17 @@
   Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
   Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
+
+  NCharged = 0;
+  NNeutrals = 0;
+  Beta = 0.0;
+  BetaStar = 0.0;
+  MeanSqDeltaR = 0.0;
+  PTD = 0.0;
+  FracPt[0] = 0.0;
+  FracPt[1] = 0.0;
+  FracPt[2] = 0.0;
+  FracPt[3] = 0.0;
+  FracPt[4] = 0.0;
+
   fArray = 0;
 }
Index: /trunk/classes/DelphesClasses.h
===================================================================
--- /trunk/classes/DelphesClasses.h	(revision 1347)
+++ /trunk/classes/DelphesClasses.h	(revision 1348)
@@ -308,4 +308,16 @@
   TLorentzVector P4();
 
+  // -- PileUpJetID variables ---
+
+  Int_t    NCharged;
+  Int_t    NNeutrals;
+  Float_t  Beta;
+  Float_t  BetaStar;
+  Float_t  MeanSqDeltaR;
+  Float_t  PTD;
+  Float_t  FracPt[5];
+
+
+
   ClassDef(Jet, 2)
 };
@@ -421,4 +433,14 @@
   virtual void Clear(Option_t* option = "");
 
+// -- PileUpJetID variables ---
+
+  Int_t    NCharged;
+  Int_t    NNeutrals;
+  Float_t  Beta;
+  Float_t  BetaStar;
+  Float_t  MeanSqDeltaR;
+  Float_t  PTD;
+  Float_t  FracPt[5];
+
 private:
   DelphesFactory *fFactory; //!
Index: /trunk/modules/ModulesLinkDef.h
===================================================================
--- /trunk/modules/ModulesLinkDef.h	(revision 1347)
+++ /trunk/modules/ModulesLinkDef.h	(revision 1348)
@@ -32,4 +32,5 @@
 #include "modules/JetPileUpSubtractor.h"
 #include "modules/TrackPileUpSubtractor.h"
+#include "modules/PileUpJetID.h"
 #include "modules/ConstituentFilter.h"
 #include "modules/StatusPidFilter.h"
@@ -64,4 +65,5 @@
 #pragma link C++ class JetPileUpSubtractor+;
 #pragma link C++ class TrackPileUpSubtractor+;
+#pragma link C++ class PileUpJetID+;
 #pragma link C++ class ConstituentFilter+;
 #pragma link C++ class StatusPidFilter+;
Index: /trunk/modules/PileUpJetID.cc
===================================================================
--- /trunk/modules/PileUpJetID.cc	(revision 1348)
+++ /trunk/modules/PileUpJetID.cc	(revision 1348)
@@ -0,0 +1,218 @@
+/** \class PileUpJetID
+ *
+ *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *
+ *  \author S. Zenz, December 2013
+ *
+ *
+ */
+
+#include "modules/PileUpJetID.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+PileUpJetID::PileUpJetID() :
+  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
+{
+
+}
+
+//------------------------------------------------------------------------------
+
+PileUpJetID::~PileUpJetID()
+{
+
+}
+
+//------------------------------------------------------------------------------
+
+void PileUpJetID::Init()
+{
+  fJetPTMin = GetDouble("JetPTMin", 20.0);
+  fParameterR = GetDouble("ParameterR", 0.5);
+  fUseConstituents = GetInt("UseConstituents", 0);
+
+  fAverageEachTower = false; // for timing
+
+  // import input array(s)
+
+  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
+  fItJetInputArray = fJetInputArray->MakeIterator();
+
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
+  fItTrackInputArray = fTrackInputArray->MakeIterator();
+
+  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
+  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
+
+  // create output array(s)
+
+  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
+
+}
+
+//------------------------------------------------------------------------------
+
+void PileUpJetID::Finish()
+{
+
+  if(fItJetInputArray) delete fItJetInputArray;
+  if(fItTrackInputArray) delete fItTrackInputArray;
+  if(fItNeutralInputArray) delete fItNeutralInputArray;
+
+}
+
+//------------------------------------------------------------------------------
+
+void PileUpJetID::Process()
+{
+  Candidate *candidate, *constituent;
+  TLorentzVector momentum, area;
+
+  Candidate *trk;
+
+  // loop over all input candidates
+  fItJetInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
+  {
+    momentum = candidate->Momentum;
+    area = candidate->Area;
+
+    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) {
+      TIter itConstituents(candidate->GetCandidates());
+      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
+        float pt = constituent->Momentum.Pt();
+        float 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) {
+	    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;
+	  }
+	}
+      }
+    } 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->IsPU) {
+	    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;
+	    }
+	  }
+	}
+      }
+      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 = sumptchpu/sumptch;
+      candidate->BetaStar = sumptchpv/sumptch;
+    } else {
+      candidate->Beta = -999.;
+      candidate->BetaStar = -999.;
+    }
+    if (sumptsq > 0.) {
+      candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
+    } else {
+      candidate->MeanSqDeltaR = -999.;
+    }
+    candidate->NCharged = nc;
+    candidate->NNeutrals = nn;
+    if (sumpt > 0.) {
+      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
+      for (int 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.;
+      }
+    }
+
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
+
Index: /trunk/modules/PileUpJetID.h
===================================================================
--- /trunk/modules/PileUpJetID.h	(revision 1348)
+++ /trunk/modules/PileUpJetID.h	(revision 1348)
@@ -0,0 +1,61 @@
+#ifndef PileUpJetID_h
+#define PileUpJetID_h
+
+/** \class PileUpJetID
+ *
+ *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *
+ *  \author S. Zenz, December 2013
+ *  
+ *
+ */
+
+
+#include "classes/DelphesModule.h"
+
+#include <deque>
+
+class TObjArray;
+class DelphesFormula;
+
+class PileUpJetID: public DelphesModule
+{
+public:
+
+  PileUpJetID();
+  ~PileUpJetID();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  Double_t fJetPTMin;
+  Double_t fParameterR;
+
+  // 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; 
+
+  Bool_t fAverageEachTower;
+
+  TIterator *fItJetInputArray; //!
+
+  const TObjArray *fJetInputArray; //!
+
+  const TObjArray *fTrackInputArray; // SCZ
+  const TObjArray *fNeutralInputArray; 
+
+  TIterator *fItTrackInputArray; // SCZ
+  TIterator *fItNeutralInputArray; // SCZ
+
+  TObjArray *fOutputArray; //!
+
+
+  ClassDef(PileUpJetID, 1)
+};
+
+#endif
+
Index: /trunk/modules/TrackPileUpSubtractor.cc
===================================================================
--- /trunk/modules/TrackPileUpSubtractor.cc	(revision 1347)
+++ /trunk/modules/TrackPileUpSubtractor.cc	(revision 1348)
@@ -102,5 +102,5 @@
   TIterator *iterator;
   TObjArray *array;
-  Double_t z, zvtx;
+  Double_t z, zvtx=0;
 
   
Index: /trunk/modules/TreeWriter.cc
===================================================================
--- /trunk/modules/TreeWriter.cc	(revision 1347)
+++ /trunk/modules/TreeWriter.cc	(revision 1348)
@@ -533,4 +533,18 @@
 
     entry->EhadOverEem = ecalEnergy > 0.0 ? hcalEnergy/ecalEnergy : 999.9;
+  
+    //---   Pile-Up Jet ID variables ----
+
+    entry->NCharged = candidate->NCharged;
+    entry->NNeutrals = candidate->NNeutrals;
+    entry->Beta = candidate->Beta;
+    entry->BetaStar = candidate->BetaStar;
+    entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
+    entry->PTD = candidate->PTD;
+    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];
 
     FillParticles(candidate, &entry->Particles);
