Index: modules/Calorimeter.cc
===================================================================
--- modules/Calorimeter.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/Calorimeter.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -150,4 +150,12 @@
 */
 
+  // 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);
@@ -356,13 +364,4 @@
       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;
@@ -380,5 +379,4 @@
       position = track->Position;
 
-
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -387,9 +385,26 @@
       fTrackHCalEnergy += hcalEnergy;
 
-      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
-      fTrackHCalTimeWeight += TMath::Sqrt(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;
+	}
+      }
+
 
       fTowerTrackArray->Add(track);
@@ -397,5 +412,5 @@
       continue;
     }
-
+  
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
@@ -412,10 +427,17 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
-    fTowerHCalTimeWeight += TMath::Sqrt(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;
+	
+      }
+    }
 
     fTower->AddCandidate(particle);
@@ -434,6 +456,5 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
-  Double_t ecalTime, hcalTime, time;
-
+ 
   if(!fTower) return;
 
@@ -444,7 +465,4 @@
   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);
@@ -454,6 +472,5 @@
 
   energy = ecalEnergy + hcalEnergy;
-  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
-
+    
   if(fSmearTowerCenter)
   {
@@ -469,5 +486,24 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  // 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->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
Index: modules/Calorimeter.h
===================================================================
--- modules/Calorimeter.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/Calorimeter.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -60,10 +60,7 @@
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
 
-  Double_t fTowerECalTime, fTowerHCalTime;
-  Double_t fTrackECalTime, fTrackHCalTime;
-
-  Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
-  Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
-
+  Double_t fTimingEMin;
+  Bool_t fElectronsFromTrack; // for timing
+  
   Int_t fTowerTrackHits, fTowerPhotonHits;
 
Index: modules/FastJetFinder.cc
===================================================================
--- modules/FastJetFinder.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/FastJetFinder.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -66,4 +66,8 @@
 #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;
 using namespace fastjet;
@@ -121,4 +125,26 @@
   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 ---
   fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
@@ -260,5 +286,5 @@
   PseudoJet jet, area;
   ClusterSequence *sequence;
-  vector< PseudoJet > inputList, outputList;
+  vector< PseudoJet > inputList, outputList, subjets;
   vector< PseudoJet >::iterator itInputList, itOutputList;
   vector< TEstimatorStruct >::iterator itEstimators;
@@ -352,5 +378,85 @@
     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 ----
 
Index: modules/FastJetFinder.h
===================================================================
--- modules/FastJetFinder.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/FastJetFinder.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -83,4 +83,24 @@
   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 --------
 
Index: modules/Isolation.cc
===================================================================
--- modules/Isolation.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/Isolation.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -25,5 +25,5 @@
  *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
  *
  */
@@ -109,5 +109,8 @@
   fUsePTSum = GetBool("UsePTSum", false);
 
+  fVetoLeptons = GetBool("VetoLeptons", true);  
+  
   fClassifier->fPTMin = GetDouble("PTMin", 0.5);
+
 
   // import input array(s)
@@ -153,5 +156,5 @@
   Candidate *candidate, *isolation, *object;
   TObjArray *isolationArray;
-  Double_t sum, ratio;
+  Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
   Int_t counter;
   Double_t eta = 0.0;
@@ -180,7 +183,13 @@
 
     // loop over all input tracks
-    sum = 0.0;
+    
+    sumNeutral = 0.0;
+    sumCharged = 0.0;
+    sumChargedPU = 0.0;
+    sumAllParticles = 0.0;
+   
     counter = 0;
     itIsolationArray.Reset();
+    
     while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
     {
@@ -188,7 +197,17 @@
 
       if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
-         candidate->GetUniqueID() != isolation->GetUniqueID())
+         candidate->GetUniqueID() != isolation->GetUniqueID() &&
+	 ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) )
       {
-        sum += isolationMomentum.Pt();
+     
+        sumAllParticles += isolationMomentum.Pt();
+        if(isolation->Charge !=0) 
+	{ 
+	  sumCharged += isolationMomentum.Pt();
+          if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt();
+	} 
+  
+        else sumNeutral += isolationMomentum.Pt();
+  
         ++counter;
       }
@@ -209,10 +228,19 @@
     }
 
-    // correct sum for pile-up contamination
-    sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
-
-    ratio = sum/candidateMomentum.Pt();
-    if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
-
+        
+     // 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;
     fOutputArray->Add(candidate);
   }
Index: modules/Isolation.h
===================================================================
--- modules/Isolation.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/Isolation.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -59,4 +59,6 @@
   Bool_t fUsePTSum;
 
+  Bool_t fVetoLeptons;  
+ 
   IsolationClassifier *fClassifier; //!
 
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/ModulesLinkDef.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -58,4 +58,5 @@
 #include "modules/Weighter.h"
 #include "modules/Hector.h"
+#include "modules/RunPUPPI.h"
 #include "modules/ExampleModule.h"
 
@@ -98,4 +99,5 @@
 #pragma link C++ class Weighter+;
 #pragma link C++ class Hector+;
+#pragma link C++ class RunPUPPI+;
 #pragma link C++ class ExampleModule+;
 
Index: modules/PileUpJetID.cc
===================================================================
--- modules/PileUpJetID.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpJetID.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -1,26 +1,8 @@
-/*
- *  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, based on http://cds.cern.ch/record/1581583
+ *  CMS PileUp Jet ID Variables
  *
- *  \author S. Zenz, December 2013
+ *  \author S. Zenz
  *
  */
@@ -41,5 +23,5 @@
 #include "TRandom3.h"
 #include "TObjArray.h"
-#include "TDatabasePDG.h"
+//#include "TDatabasePDG.h"
 #include "TLorentzVector.h"
 
@@ -54,5 +36,5 @@
 
 PileUpJetID::PileUpJetID() :
-  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
+  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
 {
 
@@ -70,8 +52,36 @@
 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
 
@@ -81,18 +91,25 @@
   fItJetInputArray = fJetInputArray->MakeIterator();
 
-  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
+
+  //  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"));
   fItTrackInputArray = fTrackInputArray->MakeIterator();
 
-  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
+  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
   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;
+
 }
 
@@ -101,8 +118,10 @@
 void PileUpJetID::Finish()
 {
+  //  cout << "In finish" << endl;
+
   if(fItJetInputArray) delete fItJetInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
   if(fItNeutralInputArray) delete fItNeutralInputArray;
-  if(fItVertexInputArray) delete fItVertexInputArray;
+
 }
 
@@ -111,24 +130,13 @@
 void PileUpJetID::Process()
 {
+  //  cout << "start of process" << endl;
+
   Candidate *candidate, *constituent;
   TLorentzVector momentum, area;
-  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;
-    }
-  }
+
+  //  cout << "BeforE SCZ additions in process" << endl;
+
+  // SCZ
+  Candidate *trk;
 
   // loop over all input candidates
@@ -139,152 +147,197 @@
     area = candidate->Area;
 
-    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)
-    {
+    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) {
       TIter itConstituents(candidate->GetCandidates());
-      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
-    {
+      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 {
       // Not using constituents, using dr
       fItTrackInputArray->Reset();
-      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;
-            }
-          }
-        }
-      }
-
+      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;
+	    }
+	  }
+	}
+      }
       fItNeutralInputArray->Reset();
-      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)
-    {
+      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.) {
       candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
-    }
-    else
-    {
-      candidate->MeanSqDeltaR = -999.0;
+    } else {
+      candidate->MeanSqDeltaR = -999.;
     }
     candidate->NCharged = nc;
     candidate->NNeutrals = nn;
-    if(sumpt > 0.0)
-    {
+    if (sumpt > 0.) {
       candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
-      for(i = 0; i < 5; ++i)
-      {
+      for (int i = 0 ; i < 5 ; i++) {
         candidate->FracPt[i] = pt_ann[i]/sumpt;
       }
-    }
-    else
-    {
-      candidate->PTD = -999.0;
-      for(i = 0; i < 5; ++i)
-      {
-        candidate->FracPt[i] = -999.0;
+    } else {
+      candidate->PTD = -999.;
+      for (int i = 0 ; i < 5 ; i++) {
+        candidate->FracPt[i] = -999.;
       }
     }
 
     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;
+	  }
+	}
+      }
+    }
+
+
   }
 }
 
 //------------------------------------------------------------------------------
-
Index: modules/PileUpJetID.h
===================================================================
--- modules/PileUpJetID.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpJetID.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -1,20 +1,2 @@
-/*
- *  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
@@ -22,7 +4,7 @@
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *  CMS PileUp Jet ID Variables
  *
- *  \author S. Zenz, December 2013
+ *  \author S. Zenz
  *
  */
@@ -34,4 +16,5 @@
 
 class TObjArray;
+class DelphesFormula;
 
 class PileUpJetID: public DelphesModule
@@ -51,8 +34,59 @@
   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;
@@ -62,20 +96,16 @@
   const TObjArray *fJetInputArray; //!
 
-  const TObjArray *fTrackInputArray; //!
-  const TObjArray *fNeutralInputArray; //!
+  const TObjArray *fTrackInputArray; // SCZ
+  const TObjArray *fNeutralInputArray; 
 
-  TIterator *fItTrackInputArray; //!
-  TIterator *fItNeutralInputArray; //!
+  TIterator *fItTrackInputArray; // SCZ
+  TIterator *fItNeutralInputArray; // SCZ
 
   TObjArray *fOutputArray; //!
+  TObjArray *fNeutralsInPassingJets; // SCZ
 
-  TIterator *fItVertexInputArray; //!
-  const TObjArray *fVertexInputArray; //!
 
-  Double_t fZVertexResolution;
-
-  ClassDef(PileUpJetID, 1)
+  ClassDef(PileUpJetID, 2)
 };
 
 #endif
-
Index: modules/PileUpMerger.cc
===================================================================
--- modules/PileUpMerger.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpMerger.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -80,4 +80,9 @@
   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
 
@@ -111,10 +116,10 @@
   TParticlePDG *pdgParticle;
   Int_t pid;
-  Float_t x, y, z, t;
+  Float_t x, y, z, t, vx, vy;
   Float_t px, py, pz, e;
   Double_t dz, dphi, dt;
-  Int_t numberOfEvents, event;
+  Int_t numberOfEvents, event, numberOfParticles;
   Long64_t allEntries, entry;
-  Candidate *candidate, *vertexcandidate;
+  Candidate *candidate, *vertex;
   DelphesFactory *factory;
 
@@ -123,5 +128,5 @@
   fItInputArray->Reset();
 
-  // --- Deal with Primary vertex first  ------
+  // --- Deal with primary vertex first  ------
 
   fFunction->GetRandom2(dz, dt);
@@ -129,7 +134,11 @@
   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();
@@ -139,9 +148,15 @@
   }
 
+  if(numberOfParticles > 0)
+  {
+    vx /= numberOfParticles;
+    vy /= numberOfParticles;
+  }
+
   factory = GetFactory();
 
-  vertexcandidate = factory->NewCandidate();
-  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
-  fVertexOutputArray->Add(vertexcandidate);
+  vertex = factory->NewCandidate();
+  vertex->Position.SetXYZT(vx, vy, dz, dt);
+  fVertexOutputArray->Add(vertex);
 
   // --- Then with pile-up vertices  ------
@@ -181,10 +196,7 @@
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    vertexcandidate = factory->NewCandidate();
-    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
-    vertexcandidate->IsPU = 1;
-
-    fVertexOutputArray->Add(vertexcandidate);
-
+    vx = 0.0;
+    vy = 0.0;
+    numberOfParticles = 0;
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     {
@@ -204,12 +216,30 @@
       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);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/PileUpMerger.h
===================================================================
--- modules/PileUpMerger.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpMerger.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -53,4 +53,9 @@
   Double_t fTVertexSpread;
 
+  Double_t fInputBeamSpotX;
+  Double_t fInputBeamSpotY;
+  Double_t fOutputBeamSpotX;
+  Double_t fOutputBeamSpotY;
+
   DelphesTF2 *fFunction; //!
 
Index: modules/PileUpMergerPythia8.cc
===================================================================
--- modules/PileUpMergerPythia8.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpMergerPythia8.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -29,5 +29,5 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
+#include "classes/DelphesTF2.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -56,6 +56,7 @@
 
 PileUpMergerPythia8::PileUpMergerPythia8() :
-  fPythia(0), fItInputArray(0)
-{
+  fFunction(0), fPythia(0), fItInputArray(0)
+{
+  fFunction = new DelphesTF2;
 }
 
@@ -64,4 +65,5 @@
 PileUpMergerPythia8::~PileUpMergerPythia8()
 {
+  delete fFunction;
 }
 
@@ -72,8 +74,20 @@
   const char *fileName;
 
+  fPileUpDistribution = GetInt("PileUpDistribution", 0);
+
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
+
+  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);
 
   fPTMin = GetDouble("PTMin", 0.0);
+
+  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
+  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
 
   fileName = GetString("ConfigFile", "MinBias.cmnd");
@@ -86,5 +100,6 @@
 
   // create output arrays
-  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
+  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
 }
 
@@ -103,29 +118,79 @@
   TParticlePDG *pdgParticle;
   Int_t pid, status;
-  Float_t x, y, z, t;
+  Float_t x, y, z, t, vx, vy;
   Float_t px, py, pz, e;
-  Double_t dz, dphi;
-  Int_t poisson, event, i;
-  Candidate *candidate;
+  Double_t dz, dphi, dt;
+  Int_t numberOfEvents, event, numberOfParticles, i;
+  Candidate *candidate, *vertex;
   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())))
   {
-    fOutputArray->Add(candidate);
+    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;
   }
 
   factory = GetFactory();
 
-  poisson = gRandom->Poisson(fMeanPileUp);
-
-  for(event = 0; event < poisson; ++event)
+  vertex = factory->NewCandidate();
+  vertex->Position.SetXYZT(vx, vy, dz, dt);
+  fVertexOutputArray->Add(vertex);
+
+  // --- 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)
   {
     while(!fPythia->next());
 
-    dz = gRandom->Gaus(0.0, fZVertexSpread);
+   // --- 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
+
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    for(i = 1; i < fPythia->event.size(); ++i)
+    vx = 0.0;
+    vy = 0.0;
+    numberOfParticles = fPythia->event.size();
+    for(i = 1; i < numberOfParticles; ++i)
     {
       Pythia8::Particle &particle = fPythia->event[i];
@@ -143,5 +208,5 @@
       candidate->PID = pid;
 
-      candidate->Status = status;
+      candidate->Status = 1;
 
       pdgParticle = pdg->GetParticle(pid);
@@ -154,12 +219,29 @@
       candidate->Momentum.RotateZ(dphi);
 
-      candidate->Position.SetXYZT(x, y, z + dz, t);
+      x -= fInputBeamSpotX;
+      y -= fInputBeamSpotY;
+      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
       candidate->Position.RotateZ(dphi);
-
-      fOutputArray->Add(candidate);
+      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
+
+      vx += candidate->Position.X();
+      vy += candidate->Position.Y();
+
+      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);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/PileUpMergerPythia8.h
===================================================================
--- modules/PileUpMergerPythia8.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/PileUpMergerPythia8.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -31,4 +31,5 @@
 
 class TObjArray;
+class DelphesTF2;
 
 namespace Pythia8
@@ -50,7 +51,18 @@
 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; //!
@@ -60,5 +72,6 @@
   const TObjArray *fInputArray; //!
 
-  TObjArray *fOutputArray; //!
+  TObjArray *fParticleOutputArray; //!
+  TObjArray *fVertexOutputArray; //!
 
   ClassDef(PileUpMergerPythia8, 1)
Index: modules/RunPUPPI.cc
===================================================================
--- modules/RunPUPPI.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
+++ modules/RunPUPPI.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -0,0 +1,262 @@
+#include "modules/RunPUPPI.h"
+
+#include "PUPPI/puppiCleanContainer.hh"
+#include "PUPPI/RecoObj.hh"
+#include "PUPPI/puppiParticle.hh"
+#include "PUPPI/puppiAlgoBin.hh"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include <algorithm> 
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+RunPUPPI::RunPUPPI() :
+  fItTrackInputArray(0), 
+  fItNeutralInputArray(0)
+{}
+
+//------------------------------------------------------------------------------
+RunPUPPI::~RunPUPPI(){}
+
+//------------------------------------------------------------------------------
+
+void RunPUPPI::Init(){
+
+  // input collection
+  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
+  fItTrackInputArray   = fTrackInputArray->MakeIterator();
+  fNeutralInputArray   = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
+  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
+  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
+  fPVItInputArray      = fPVInputArray->MakeIterator();
+
+
+  // puppi parameters                                     
+  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
+  fUseExp         = GetBool("UseExp", false);
+
+  // read eta min ranges                                                                                                                                                           
+  ExRootConfParam param = GetParam("EtaMinBin");
+  fEtaMinBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
+
+  // read eta max ranges                                                                                                                                                           
+  param = GetParam("EtaMaxBin");
+  fEtaMaxBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
+
+  // read pt min value                                                                                                                                                           
+  param = GetParam("PtMinBin");
+  fPtMinBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
+
+  // read cone size                                                                                                                                                           
+  param = GetParam("ConeSizeBin");
+  fConeSizeBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
+
+  // read RMS min pt                                                                                                                                             
+  param = GetParam("RMSPtMinBin");
+  fRMSPtMinBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
+
+  // read RMS scale factor                                                                                                                                                           
+  param = GetParam("RMSScaleFactorBin");
+  fRMSScaleFactorBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
+
+  // read neutral pt min cut                                                                                                                                                           
+  param = GetParam("NeutralMinEBin");
+  fNeutralMinEBin.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
+
+  // read neutral pt min slope                                                                                                                                                           
+  param = GetParam("NeutralPtSlope");
+  fNeutralPtSlope.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
+
+  // read apply chs                                                                                                                                                           
+  param = GetParam("ApplyCHS");
+  fApplyCHS.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
+
+  // read use charged                                                                                                                                                           
+  param = GetParam("UseCharged");
+  fUseCharged.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
+
+  // read apply chs correction                                                                                                                                                           
+  param = GetParam("ApplyLowPUCorr");
+  fApplyLowPUCorr.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
+  
+  // read metric id                                                                                                                                                          
+  param = GetParam("MetricId");
+  fMetricId.clear();
+  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
+
+  // create output array
+  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
+  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
+  fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
+}
+
+//------------------------------------------------------------------------------
+
+void RunPUPPI::Finish(){
+  if(fItTrackInputArray)   delete fItTrackInputArray;
+  if(fItNeutralInputArray) delete fItNeutralInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void RunPUPPI::Process(){
+
+  Candidate *candidate, *particle;
+  TLorentzVector momentum;
+
+  DelphesFactory *factory = GetFactory();
+
+  // loop over input objects
+  fItTrackInputArray->Reset();
+  fItNeutralInputArray->Reset();
+  fPVItInputArray->Reset();
+
+  std::vector<Candidate> InputParticles;
+  InputParticles.clear();
+
+  // take the leading vertex 
+  float PVZ = 0.;
+  Candidate *pv = static_cast<Candidate*>(fPVItInputArray->Next());
+  if (pv) PVZ = pv->Position.Z();
+
+  // Fill input particles for puppi
+  std::vector<RecoObj> puppiInputVector;
+
+  // Loop on charge track candidate
+  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
+
+      momentum = candidate->Momentum;
+      
+      RecoObj curRecoObj;
+      curRecoObj.pt  = momentum.Pt();
+      curRecoObj.eta = momentum.Eta();
+      curRecoObj.phi = momentum.Phi();
+      curRecoObj.m   = momentum.M();  
+      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
+      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
+	curRecoObj.id    = 2;
+	curRecoObj.vtxId = candidate->IsPU;
+	if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
+	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
+	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
+        else curRecoObj.pfType = 1;
+        curRecoObj.dZ = particle->Position.Z()-PVZ;
+      } 
+      else if(!candidate->IsRecoPU and candidate->Charge !=0) {
+	curRecoObj.id    = 1;  // charge from LV
+        curRecoObj.vtxId = 1; // from PV
+	if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
+	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
+	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
+        else curRecoObj.pfType = 1;
+        curRecoObj.dZ = particle->Position.Z()-PVZ;
+      }
+      else {
+	std::cerr<<" RunPUPPI: problem with a charged track --> it has charge 0 "<<std::endl;
+        continue;
+      }
+
+      puppiInputVector.push_back(curRecoObj);
+      InputParticles.push_back(*candidate);
+  }
+
+  // Loop on neutral calo cells 
+  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
+
+      momentum = candidate->Momentum;
+
+      RecoObj curRecoObj;
+      curRecoObj.pt  = momentum.Pt();
+      curRecoObj.eta = momentum.Eta();
+      curRecoObj.phi = momentum.Phi();
+      curRecoObj.m   = momentum.M();
+      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
+
+
+      if(candidate->Charge == 0){
+	curRecoObj.id    = 0; // neutrals have id==0
+	curRecoObj.vtxId = 0; // neutrals have vtxId==0 
+	if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
+	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
+	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
+        else curRecoObj.pfType = 5;
+        curRecoObj.dZ = particle->Position.Z()-PVZ;
+      }
+      else{
+	std::cerr<<" RunPUPPI: problem with a neutrals cells --> it has charge !=0 "<<std::endl;
+        continue;
+      }
+      puppiInputVector.push_back(curRecoObj);
+      InputParticles.push_back(*candidate);
+  }
+
+  // Create algorithm list for puppi
+  std::vector<puppiAlgoBin> puppiAlgo;
+  if(puppiAlgo.empty()){
+   if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
+       and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size() 
+       and fEtaMinBin.size() == fApplyCHS.size()  and fEtaMinBin.size() == fUseCharged.size()
+       and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
+    std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
+    std::exit(EXIT_FAILURE);
+   } 
+
+   for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
+    puppiAlgoBin algoTmp ;
+    algoTmp.fEtaMin_ = fEtaMinBin.at(iAlgo);
+    algoTmp.fEtaMax_ = fEtaMaxBin.at(iAlgo);
+    algoTmp.fPtMin_  = fPtMinBin.at(iAlgo);
+    algoTmp.fConeSize_        = fConeSizeBin.at(iAlgo);
+    algoTmp.fRMSPtMin_        = fRMSPtMinBin.at(iAlgo);
+    algoTmp.fRMSScaleFactor_  = fRMSScaleFactorBin.at(iAlgo);
+    algoTmp.fNeutralMinE_     = fNeutralMinEBin.at(iAlgo);
+    algoTmp.fNeutralPtSlope_  = fNeutralPtSlope.at(iAlgo);
+    algoTmp.fApplyCHS_        = fApplyCHS.at(iAlgo);
+    algoTmp.fUseCharged_      = fUseCharged.at(iAlgo);
+    algoTmp.fApplyLowPUCorr_  = fApplyLowPUCorr.at(iAlgo);
+    algoTmp.fMetricId_        = fMetricId.at(iAlgo);
+    if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;    
+    puppiAlgo.push_back(algoTmp);     
+   }
+  }  
+
+  // Create PUPPI container
+  puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
+  std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
+
+  // Loop on final particles
+  for (std::vector<fastjet::PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
+    if(it->user_index() <= int(InputParticles.size())){      
+      candidate = factory->NewCandidate();
+      candidate = dynamic_cast<Candidate*>(InputParticles.at(it->user_index()).Clone());
+      candidate->Momentum.SetXYZT(it->px(),it->py(),it->pz(),it->e());
+      fOutputArray->Add(candidate);
+      if( puppiInputVector.at(it->user_index()).id == 1 or puppiInputVector.at(it->user_index()).id == 2) fOutputTrackArray->Add(candidate);
+      else if (puppiInputVector.at(it->user_index()).id == 0) fOutputNeutralArray->Add(candidate);
+    }
+    else{ 
+      std::cerr<<" particle not found in the input Array --> skip "<<std::endl;
+      continue;
+    }     
+  }
+
+}
Index: modules/RunPUPPI.h
===================================================================
--- modules/RunPUPPI.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
+++ modules/RunPUPPI.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -0,0 +1,56 @@
+#ifndef RunPUPPI_h
+#define RunPUPPI_h
+
+#include "classes/DelphesModule.h"
+#include <vector>
+
+class TObjArray;
+class TIterator;
+
+
+class RunPUPPI: public DelphesModule {
+
+ public:
+
+  RunPUPPI();
+  ~RunPUPPI();
+
+  void Init();
+  void Process();
+  void Finish();
+  
+ private:
+
+  TIterator *fItTrackInputArray;
+  TIterator *fItNeutralInputArray; //!
+  TIterator *fPVItInputArray; //!                                                                                                                                                      
+  
+  const TObjArray *fTrackInputArray;
+  const TObjArray *fNeutralInputArray; //!
+  const TObjArray *fPVInputArray; //!                                                                                                                                                     
+ 
+  // puppi parameters
+  float fMinPuppiWeight;
+  bool fUseExp;
+  
+  std::vector<float> fEtaMinBin ;
+  std::vector<float> fEtaMaxBin ;
+  std::vector<float> fPtMinBin ;
+  std::vector<float> fConeSizeBin ;
+  std::vector<float> fRMSPtMinBin ;
+  std::vector<float> fRMSScaleFactorBin ;
+  std::vector<float> fNeutralMinEBin;
+  std::vector<float> fNeutralPtSlope;
+  std::vector<bool>  fApplyCHS;
+  std::vector<bool>  fUseCharged;
+  std::vector<bool>  fApplyLowPUCorr;
+  std::vector<int>   fMetricId;
+
+  TObjArray *fOutputArray;
+  TObjArray *fOutputTrackArray;
+  TObjArray *fOutputNeutralArray;
+
+  ClassDef(RunPUPPI, 1)
+};
+
+#endif
Index: modules/TrackPileUpSubtractor.cc
===================================================================
--- modules/TrackPileUpSubtractor.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/TrackPileUpSubtractor.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -74,6 +74,7 @@
   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;
@@ -147,7 +148,10 @@
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
       
-      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
-
-      array->Add(candidate);
+      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1;
+      else 
+      {
+         candidate->IsRecoPU = 0;
+         if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
+      }
     }
   }
Index: modules/TrackPileUpSubtractor.h
===================================================================
--- modules/TrackPileUpSubtractor.h	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/TrackPileUpSubtractor.h	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -50,4 +50,6 @@
   Double_t fZVertexResolution;
 
+  Double_t fPTMin; 
+
   std::map< TIterator *, TObjArray * > fInputMap; //!
 
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision fc6300dd9fff8e3db0f7a612ff1e7368f74616d6)
+++ modules/TreeWriter.cc	(revision dd514ae16e0b84339c012c1e881b4cf794f46541)
@@ -362,4 +362,5 @@
 
     entry->T = position.T()*1.0E-3/c_light;
+    entry->Ntimes = candidate->Ntimes;    
 
     FillParticles(candidate, &entry->Particles);
@@ -403,4 +404,13 @@
     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;
 
@@ -442,4 +452,14 @@
     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;
 
@@ -487,4 +507,13 @@
 
     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;
@@ -504,4 +533,5 @@
   Double_t ecalEnergy, hcalEnergy;
   const Double_t c_light = 2.99792458E8;
+  Int_t i;
 
   array->Sort();
@@ -531,4 +561,6 @@
 
     entry->Mass = momentum.M();
+
+    entry->Area = candidate->Area;
 
     entry->DeltaEta = candidate->DeltaEta;
@@ -561,17 +593,19 @@
     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];
-   
-    //--- 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];
+    //--- 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];
+    }
    
     FillParticles(candidate, &entry->Particles);
