Index: Makefile
===================================================================
--- Makefile	(revision 7b45ff5e3840b1a7e6b83f7376cdb74e900f182c)
+++ Makefile	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -955,10 +955,6 @@
 endif
 
-tmp/external/PUPPI/PuppiAlgo.$(ObjSuf): \
-	external/PUPPI/PuppiAlgo.$(SrcSuf) \
-	external/fastjet/internal/base.hh
-tmp/external/PUPPI/PuppiContainer.$(ObjSuf): \
-	external/PUPPI/PuppiContainer.$(SrcSuf) \
-	external/fastjet/internal/base.hh \
+tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \
+	external/PUPPI/puppiCleanContainer.$(SrcSuf) \
 	external/fastjet/Selector.hh
 tmp/external/fastjet/AreaDefinition.$(ObjSuf): \
@@ -1313,12 +1309,13 @@
 	modules/RunPUPPI.$(SrcSuf) \
 	modules/RunPUPPI.h \
-	external/PUPPI/RecoObj2.hh \
-	external/PUPPI/AlgoObj.hh \
+	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
 FASTJET_OBJ +=  \
-	tmp/external/PUPPI/PuppiAlgo.$(ObjSuf) \
-	tmp/external/PUPPI/PuppiContainer.$(ObjSuf) \
+	tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \
 	tmp/external/fastjet/AreaDefinition.$(ObjSuf) \
 	tmp/external/fastjet/BasicRandom.$(ObjSuf) \
@@ -1745,6 +1742,5 @@
 
 modules/RunPUPPI.h: \
-	classes/DelphesModule.h \
-	external/PUPPI/PuppiContainer.hh
+	classes/DelphesModule.h
 	@touch $@
 
@@ -1868,4 +1864,12 @@
 	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 $@
 
Index: cards/CMS_PhaseII/CMS_PhaseII_200PU.tcl
===================================================================
--- cards/CMS_PhaseII/CMS_PhaseII_200PU.tcl	(revision 7b45ff5e3840b1a7e6b83f7376cdb74e900f182c)
+++ cards/CMS_PhaseII/CMS_PhaseII_200PU.tcl	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -40,9 +40,5 @@
   EFlowMerger
 
-  LeptonFilterNoLep
-  LeptonFilterLep
-  RunPUPPIBase
   RunPUPPI
-
 
   PhotonIsolation
@@ -59,5 +55,4 @@
 
   MissingET
-  PuppiMissingET
   GenMissingET
   GenPileUpMissingET
@@ -96,5 +91,5 @@
 
   # pre-generated minbias input file
-  set PileUpFile MinBias.pileup 
+  set PileUpFile ../eos/cms/store/group/upgrade/delphes/PhaseII/MinBias_100k.pileup 
 
   # average expected pile up
@@ -110,8 +105,10 @@
  
   #set VertexDistributionFormula {exp(-(t^2/(2*(0.05/2.99792458E8*exp(-(z^2/(2*(0.05)^2))))^2)))}
-  set VertexDistributionFormula { (abs(t) <= 1.6e-10) * (abs(z) <= 0.053) * (1.00) +
-      (abs(t) >  1.6e-10) * (abs(z) <= 0.053) * (0.00) +
-      (abs(t) <= 1.6e-10) * (abs(z) > 0.053)  * (0.00) +
-      (abs(t) >  1.6e-10) * (abs(z) > 0.053)  * (0.00)}
+  
+  
+  set VertexDistributionFormula { (abs(t) <= 1.0e-09) * (abs(z) <= 0.15) * (1.00) +
+                                  (abs(t) >  1.0e-09) * (abs(z) <= 0.15) * (0.00) +
+   				  (abs(t) <= 1.0e-09) * (abs(z) > 0.15)  * (0.00) +
+   				  (abs(t) >  1.0e-09) * (abs(z) > 0.15)  * (0.00)}
 
 }
@@ -550,47 +547,25 @@
 #########################################
 
-module PdgCodeFilter LeptonFilterNoLep {
-  set InputArray HCal/eflowTracks
-  set OutputArray eflowTracksNoLeptons
-  set Invert false
-  add PdgCode {13}
-  add PdgCode {-13}
-  add PdgCode {11}
-  add PdgCode {-11}
-}
-
-module PdgCodeFilter LeptonFilterLep {
-  set InputArray HCal/eflowTracks
-  set OutputArray eflowTracksLeptons
-  set Invert true
-  add PdgCode {11}
-  add PdgCode {-11}
-  add PdgCode {13}
-  add PdgCode {-13}
-}
-
-module RunPUPPI RunPUPPIBase {
+module RunPUPPI RunPUPPI {
   ## input information
-  set TrackInputArray   LeptonFilterNoLep/eflowTracksNoLeptons
+  set TrackInputArray   TrackMerger/tracks
   set NeutralInputArray NeutralEFlowMerger/eflowTowers
   set PVInputArray      PileUpMerger/vertices
   set MinPuppiWeight    0.05
   set UseExp            false
-  set UseNoLep          false
  
   ## define puppi algorithm parameters (more than one for the same eta region is possible)                                                                                      
-  add EtaMinBin           0.0   1.5   4.0
-  add EtaMaxBin           1.5   4.0   10.0
-  add PtMinBin            0.0   0.0   0.0
-  add ConeSizeBin         0.2   0.2   0.2
-  add RMSPtMinBin         0.1   0.5   0.5
-  add RMSScaleFactorBin   1.0   1.0   1.0
-  add NeutralMinEBin      0.2   0.2   0.5
-  add NeutralPtSlope      0.006 0.013 0.067
-  add ApplyCHS            true  true  true
-  add UseCharged          true  true  false
-  add ApplyLowPUCorr      true  true  true 
-  add MetricId            5     5     5
-  add CombId              0     0     0
+  add EtaMinBin           0.    2.5    2.5    3.0   3.0
+  add EtaMaxBin           2.5   3.0    3.0    10.0  10.0
+  add PtMinBin            0.    0.5    0.5    0.5   0.5
+  add ConeSizeBin         0.25  0.25   0.25   0.25  0.25
+  add RMSPtMinBin         0.1   0.5    0.5    0.5   0.5
+  add RMSScaleFactorBin   1.0   1.0    1.0    1.0   1.0
+  add NeutralMinEBin      0.2   1.0    1.0    1.5   1.5
+  add NeutralPtSlope      0.02  0.02   0.02   0.02  0.02
+  add ApplyCHS            true  true   true   true  true
+  add UseCharged          true  false  false  false false
+  add ApplyLowPUCorr      true  true   true   true  true
+  add MetricId            5     5      1      5     1
 
   ## output name
@@ -600,9 +575,5 @@
 } 
 
-module Merger RunPUPPI {
-  add InputArray RunPUPPIBase/PuppiParticles
-  add InputArray LeptonFilterLep/eflowTracksLeptons
-  set OutputArray PuppiParticles
-}
+
 
 ###################
@@ -614,11 +585,4 @@
 #  add InputArray RunPUPPI/PuppiParticles
   add InputArray EFlowMerger/eflow
-  set MomentumOutputArray momentum
-}
-
-module Merger PuppiMissingET {
-  #add InputArray InputArray
-  add InputArray RunPUPPI/PuppiParticles
-  #add InputArray EFlowMerger/eflow
   set MomentumOutputArray momentum
 }
@@ -1007,5 +971,4 @@
   
   add Branch MissingET/momentum MissingET MissingET
-  add Branch PuppiMissingET/momentum PuppiMissingET MissingET
   add Branch GenPileUpMissingET/momentum GenPileUpMissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
Index: external/PUPPI/RecoObj.hh
===================================================================
--- external/PUPPI/RecoObj.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
+++ external/PUPPI/RecoObj.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -0,0 +1,33 @@
+#ifndef BACONANA_DATAFORMATS_RECOOBJ_HH
+#define BACONANA_DATAFORMATS_RECOOBJ_HH
+
+class RecoObj 
+{
+public:
+      RecoObj():
+	pt(0), eta(0), phi(0), m(0),
+	pfType(-1),vtxId(-1),
+	trkChi2(0),vtxChi2(0),
+	id(0),time(0),depth(0)
+    {}
+    ~RecoObj(){}
+    
+    float         pt, eta, phi, m;  // kinematics
+    int           pfType;
+    int           vtxId;               // Vertex Id from Vertex Collection
+    float         trkChi2;             // Track Chi2
+    float         vtxChi2;             // Vertex Chi2
+    int           id;
+    float         time,depth;    // Usefule Info
+    float         expProb;
+    float         expChi2PU;
+    float         expChi2;
+    float         dZ;
+    float         d0;
+};
+#endif
+
+// Only need 4-vector and ID
+// id = 0:   neutral
+// id = 1; charged from PV
+// id = 2; charged not from PV
Index: external/PUPPI/puppiAlgoBin.hh
===================================================================
--- external/PUPPI/puppiAlgoBin.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
+++ external/PUPPI/puppiAlgoBin.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -0,0 +1,125 @@
+#ifndef PUPPIALGOBIN_HH
+#define PUPPIALGOBIN_HH
+
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+//...................... class to identify the puppi Algorithm parameters to be used                                                                                                       
+class puppiAlgoBin {
+
+ public :
+
+  puppiAlgoBin(){
+
+    fEtaMin_   = 0.;
+    fEtaMax_   = 0.;
+    fPtMin_    = 0.;
+    fConeSize_ = 0.;
+    fRMSPtMin_ = 0.;
+    fRMSScaleFactor_ = 1.0;
+    fNeutralMinE_    = 0.;
+    fNeutralPtSlope_ = 0.;
+    fApplyCHS_   = true;
+    fUseCharged_ = false ;
+    fApplyLowPUCorr_ = false;
+    fMetricId_ = 0;
+    fRMS_      = 0.;
+    fMean_     = 0.;
+    fMedian_   = 0.;
+  };
+
+  puppiAlgoBin(const float & etaMin,         const float & etaMax, const float & ptMin, const float & coneSize, const float & RMSPtMin,
+               const float & RMSScaleFactor, const float & neutralMinE, const float & neutralPtSlope,	                                     
+               const bool  & applyCHS,       const bool & useCharged, const bool & applyLowPUCorr, const int & metricId):
+    fEtaMin_(etaMin),
+    fEtaMax_(etaMax),
+    fPtMin_(ptMin),
+    fConeSize_(coneSize),
+    fRMSPtMin_(RMSPtMin),
+    fRMSScaleFactor_(RMSScaleFactor),
+    fNeutralMinE_(neutralMinE),
+    fNeutralPtSlope_(neutralPtSlope),
+    fApplyCHS_(applyCHS),
+    fUseCharged_(useCharged),
+    fApplyLowPUCorr_(applyLowPUCorr),
+    fMetricId_(metricId)
+  {
+    fRMS_    = 0.;
+    fMean_   = 0.;
+    fMedian_ = 0.;
+
+  };
+
+  virtual ~puppiAlgoBin(){};
+
+  bool operator == (const puppiAlgoBin & algo2) const {
+    if( fEtaMin_ == algo2.fEtaMin_ and 
+        fEtaMax_ == algo2.fEtaMax_ and 
+        fApplyCHS_ == algo2.fApplyCHS_ and 
+        fUseCharged_ == algo2.fUseCharged_ and 
+        fApplyLowPUCorr_ == algo2.fApplyLowPUCorr_ and 
+        fMetricId_ == algo2.fMetricId_ and 
+        fPtMin_ == algo2.fPtMin_ and 
+        fConeSize_ == algo2.fConeSize_ and 
+        fRMSPtMin_ == algo2.fRMSPtMin_ and 
+        fRMSScaleFactor_ == algo2.fRMSScaleFactor_ and 
+        fNeutralMinE_ == algo2.fNeutralMinE_ and 
+        fNeutralPtSlope_ == algo2.fNeutralPtSlope_ ) return true;
+    else return false;
+  };
+
+  void setPuppiParticles(const std::vector<puppiParticle> & puppiParticles){ // set the particles used by the current algorithm                                             
+    fPuppiParticlesPU_.clear();
+    fPuppiParticlesPV_.clear();
+    fPuppiParticlesNULL_.clear();
+  
+    for(size_t iPart = 0; iPart < puppiParticles.size(); iPart++){ // loop on the puppiParticles                            
+      if(puppiParticles.at(iPart).fPval_ == -999){ // default PVal
+        fPuppiParticlesNULL_.push_back(puppiParticles.at(iPart));        
+        continue ; 
+      }
+                                                 
+      if(puppiParticles.at(iPart).fPt_ <  fRMSPtMin_){  // should be discarded in the RMS computation
+        fPuppiParticlesNULL_.push_back(puppiParticles.at(iPart)); 
+        continue ; 
+      }
+     
+      if(puppiParticles.at(iPart).fPt_ < fPtMin_){ // under the chosen pt cut
+        fPuppiParticlesNULL_.push_back(puppiParticles.at(iPart)); 
+        continue ; 
+      }
+
+      // zero are neutral particle, 1 LV charged, 2PU      
+      if(fUseCharged_ && (fabs(puppiParticles.at(iPart).fParticleId_) >=0 && fabs(puppiParticles.at(iPart).fParticleId_) <=1)) fPuppiParticlesPV_.push_back(puppiParticles.at(iPart));
+      if(fUseCharged_ && fabs(puppiParticles.at(iPart).fParticleId_) < 2) continue;
+      fPuppiParticlesPU_.push_back(puppiParticles.at(iPart));
+    }
+  };
+
+  float fEtaMin_;
+  float fEtaMax_;
+  float fPtMin_ ;
+  float fConeSize_;
+  float fRMSPtMin_;
+  float fRMSScaleFactor_;
+  float fNeutralMinE_;
+  float fNeutralPtSlope_;
+
+  bool  fApplyCHS_;
+  bool  fUseCharged_;
+  bool  fApplyLowPUCorr_;
+
+  int   fMetricId_;
+
+  float fRMS_;
+  float fMean_;
+  float fMedian_;
+
+  std::vector<puppiParticle> fPuppiParticlesPU_;
+  std::vector<puppiParticle> fPuppiParticlesPV_;
+  std::vector<puppiParticle> fPuppiParticlesNULL_;
+
+};
+
+#endif
Index: external/PUPPI/puppiCleanContainer.cc
===================================================================
--- external/PUPPI/puppiCleanContainer.cc	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
+++ external/PUPPI/puppiCleanContainer.cc	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -0,0 +1,471 @@
+#include "puppiCleanContainer.hh"
+#include "Math/SpecFuncMathCore.h"
+#include "Math/ProbFunc.h"
+#include "TH2F.h"
+#include "fastjet/Selector.hh"
+
+#include <algorithm>
+
+#include "TMath.h"
+#include "Math/QuantFuncMathCore.h"
+#include "Math/SpecFuncMathCore.h"
+#include "Math/ProbFunc.h"
+
+using namespace std;
+
+// ------------- Constructor
+puppiCleanContainer::puppiCleanContainer(std::vector<RecoObj> inParticles, 
+                                         std::vector<puppiAlgoBin> puppiAlgo,
+                                         float minPuppiWeight,
+                                         bool  fUseExp){
+
+    // take the input particles
+    fRecoParticles_.clear();
+    fRecoParticles_ = inParticles;
+
+    // puppi algo 
+    puppiAlgo_.clear();
+    puppiAlgo_ = puppiAlgo; 
+
+    // min puppi weight
+    fMinPuppiWeight_ = minPuppiWeight;
+
+    //Clear everything
+    fPFParticles_.clear();
+    fPFchsParticles_.clear();
+    fChargedPV_.clear();
+    fChargedNoPV_.clear();
+    fPuppiWeights_.clear();
+
+    fNPV_    = 1 ;
+    fPVFrac_ = 0.;
+    fUseExp_ = fUseExp;
+
+    //Link to the RecoObjects --> loop on the input particles
+    for (unsigned int i = 0; i < fRecoParticles_.size(); i++){
+        fastjet::PseudoJet curPseudoJet;
+        curPseudoJet.reset_PtYPhiM (fRecoParticles_[i].pt,fRecoParticles_[i].eta,fRecoParticles_[i].phi,fRecoParticles_[i].m);
+        curPseudoJet.set_user_index(fRecoParticles_[i].id);  
+        // fill vector of pseudojets for internal references
+        fPFParticles_.push_back(curPseudoJet);
+        if(fRecoParticles_[i].id <= 1) fPFchsParticles_.push_back(curPseudoJet);    //Remove Charged particles associated to other vertex
+        if(fRecoParticles_[i].id == 1) fChargedPV_.push_back(curPseudoJet);         //Take Charged particles associated to PV
+        if(fRecoParticles_[i].id == 2) fChargedNoPV_.push_back(curPseudoJet);
+        if(fRecoParticles_[i].id >= 0) fPVFrac_++ ;
+	if(fNPV_ < fRecoParticles_[i].vtxId) fNPV_ = fRecoParticles_[i].vtxId;
+
+    }
+
+    fPVFrac_ = double(fChargedPV_.size())/fPVFrac_;
+}
+
+// ------------- De-Constructor
+puppiCleanContainer::~puppiCleanContainer(){}
+
+// main function to compute puppi Event
+std::vector<fastjet::PseudoJet> puppiCleanContainer::puppiEvent(){
+
+    // output particles
+    std::vector<fastjet::PseudoJet> particles;
+    particles.clear();
+
+    std::vector<int> pPupId ; 
+    std::vector<puppiParticle> partTmp ; // temp puppi particle vector; make a clone of the same particle for all the algo in which it is contained
+  
+    // calculate puppi metric, RMS and mean value for all the algorithms
+    for(size_t iPuppiAlgo = 0; iPuppiAlgo < puppiAlgo_.size(); iPuppiAlgo++){
+      getRMSAvg(iPuppiAlgo,fPFParticles_,fChargedPV_); // give all the particles in the event and the charged one
+    }
+  
+    int npart = 0;  
+
+    // Loop on all the incoming particles
+    for(size_t iPart = 0; iPart < fPFParticles_.size(); iPart++) {
+
+      float pWeight = 1; // default weight
+      pPupId.clear();
+      pPupId = getPuppiId(fPFParticles_[iPart].pt(),fPFParticles_[iPart].eta(),puppiAlgo_); // take into account only algo eta
+
+      //////////////////////////////////////////      
+      // acceptance check of the puppi algorithm
+      //////////////////////////////////////////      
+
+      if(pPupId.empty()) { // out acceptance... no algorithm found
+        fPuppiWeights_.push_back(pWeight); // take the particle as it is
+        fastjet::PseudoJet curjet(pWeight*fPFParticles_[iPart].px(),pWeight*fPFParticles_[iPart].py(),pWeight*fPFParticles_[iPart].pz(),pWeight*fPFParticles_[iPart].e());    
+        curjet.set_user_index(fPFParticles_[iPart].user_index());                                                                                                      
+        particles.push_back(curjet); // fill the output collection                                                                                                            
+	continue; //go to the next particle
+      }
+      
+      ///////////////      
+      //  PT check //
+      ///////////////      
+      for(size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){ // loop on all the available algo for that region       
+	if(fPFParticles_.at(iPart).pt() < puppiAlgo_.at(pPupId.at(iAlgo)).fPtMin_){ // low momentum particles should be cut by puppi method
+	pWeight = 0; // if this particle is under the pT threshold of one the algorithm, put the weight as zero      
+	break;
+	}
+      }
+
+      if(pWeight == 0){ 
+        fPuppiWeights_.push_back(0); // puppi weight is zero
+        continue;
+      }
+ 
+      /////////////////////////////////////      
+      // fill the p-values for Z-vertex //
+      /////////////////////////////////////      
+
+      double pChi2 = 0;   
+      if(fUseExp_){ // use vertex-z resolution
+       //Compute an Experimental Puppi Weight with delta Z info (very simple example)
+       if(iPart <= fRecoParticles_.size() and fRecoParticles_[iPart].id == fPFParticles_.at(iPart).user_index()){
+	pChi2 = getChi2FromdZ(fRecoParticles_[iPart].dZ); // get the probability fiven the dZ of the particle wrt the leading vertex
+        if(fRecoParticles_[iPart].pfType > 3) pChi2 = 0; // not use this info for neutrals
+       }
+      }
+      
+
+      /////////////////////////////////////      
+      // found  the particle in all the algorithm
+      /////////////////////////////////////      
+      partTmp.clear();
+      for(size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){ // loop on all the algo found
+
+       int found = 0;  //  found index
+       if(fabs(fPFParticles_[iPart].user_index()) <= 1 and puppiAlgo_.at(pPupId.at(iAlgo)).fUseCharged_){ // charged or neutral from PV
+	 for(size_t puppiIt = 0 ; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.size(); puppiIt++){ // Loop on PV particles
+	   if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.at(puppiIt).fPosition_ == int(iPart)){
+	    partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.at(puppiIt));  // take the puppi particle
+            found = 1 ;
+            break;
+	  }
+	}
+       }
+       else if ((fabs(fPFParticles_[iPart].user_index()) <= 1 and !puppiAlgo_.at(pPupId.at(iAlgo)).fUseCharged_) or fabs(fPFParticles_[iPart].user_index()) >= 2){
+	for(size_t puppiIt = 0; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.size(); puppiIt++){
+          if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.at(puppiIt).fPosition_ == int(iPart)){
+            partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.at(puppiIt));
+            found = 1;
+            break;
+	  }
+	}
+       }
+      
+      /////////////////////////////////////      
+       // means that is inside the NULL vector for some reasons
+      /////////////////////////////////////      
+
+       if(found == 0){
+	for(size_t puppiIt = 0; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.size(); puppiIt++){
+          if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.at(puppiIt).fPosition_ == int(iPart)){
+            partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.at(puppiIt));
+            found = 1 ;
+            break;
+	  }
+	}
+       }
+      }
+
+      if(partTmp.size() != pPupId.size()){ // not found the particle in one of the algorithms 
+        pWeight = 1 ;
+	fPuppiWeights_.push_back(pWeight);
+        fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());    
+        curjet.set_user_index(fPFParticles_[iPart].user_index());                                                                                                      
+        particles.push_back(curjet); // by default is one, so 4V is not chaged                                                                                                
+        continue;
+      }
+      
+      /////////////////////////////////////      
+      //Check the Pval  
+      /////////////////////////////////////      
+
+      bool badPVal = false ;
+      for(size_t iPuppi = 0; iPuppi < partTmp.size() ; iPuppi++){
+	if(partTmp.at(iPuppi).fPval_ == -999){ // if the default is found as PVal, leave the particle as it is in the output
+	 pWeight = 1 ;
+         badPVal = true ;
+	}
+      }
+
+      if(badPVal){
+ 	 fPuppiWeights_.push_back(pWeight);
+         fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());    
+         curjet.set_user_index(fPFParticles_[iPart].user_index());                                                                                                      
+         particles.push_back(curjet);                                                                                                                                                
+         continue;
+	
+      }
+
+      
+      // compute combining the weight for all the algorithm
+      pWeight = compute(pChi2,partTmp,puppiAlgo_,pPupId);      
+      
+      //Basic Weight Checks
+      if( std::isinf(pWeight) || std::isnan(pWeight)){
+	//std::cerr << "====> Weight is nan : pt " << fPFParticles_[iPart].pt() << " -- eta : " << fPFParticles_[iPart].eta() << " -- id : " << fPFParticles_[iPart].user_index() << std::endl;
+	pWeight = 1; // set the default to avoid problems
+      }
+
+      //Basic Cuts
+      if(pWeight < fMinPuppiWeight_) pWeight = 0; //==> Elminate the low Weight stuff
+      
+      //threshold cut on the neutral Pt
+      for(size_t iPuppi = 0; iPuppi < pPupId.size(); iPuppi++){
+        if(fPFParticles_[iPart].user_index() == 1 && puppiAlgo_.at(pPupId.at(iPuppi)).fApplyCHS_ ) pWeight = 1; // charged from LV 
+        if(fPFParticles_[iPart].user_index() == 2 && puppiAlgo_.at(pPupId.at(iPuppi)).fApplyCHS_ ) pWeight = 0; // charged from PU
+	if(pWeight*fPFParticles_[iPart].pt() < getNeutralPtCut(puppiAlgo_.at(pPupId.at(iPuppi)).fNeutralMinE_,puppiAlgo_.at(pPupId.at(iPuppi)).fNeutralPtSlope_,fNPV_) && fPFParticles_[iPart].user_index() == 0 ) // if don't pass one of the algo neutral pt condition the particle is cut
+       pWeight = 0; 
+      }
+
+      fPuppiWeights_.push_back(pWeight); // push back the weight
+
+      //Now get rid of the thrown out weights for the particle collection
+      if(pWeight == 0) continue; // if zero don't fill the particle in the output
+      npart++;
+
+      //Produce
+      fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());           
+      curjet.set_user_index(iPart);                                                                                                                            
+      particles.push_back(curjet);
+      
+    }
+    
+    return particles;
+      
+}
+
+// compute puppi metric, RMS and median for PU particle for each algo
+void puppiCleanContainer::getRMSAvg(const int & iPuppiAlgo, std::vector<fastjet::PseudoJet> & particlesAll, std::vector<fastjet::PseudoJet> &chargedPV) { 
+
+  std::vector<puppiParticle> puppiParticles; // puppi particles to be set for a specific algo
+  puppiParticles.clear();
+
+  // Loop on all the particles of the event  
+    
+  for(size_t iPart = 0; iPart < particlesAll.size(); iPart++ ) { 
+
+    float pVal    = -999;
+    bool  pPupId  = isGoodPuppiId(particlesAll[iPart].pt(),particlesAll[iPart].eta(),puppiAlgo_.at(iPuppiAlgo)); // get the puppi id algo asaf of eta and phi of the particle
+    // does not exsist and algorithm for this particle, store -999 as pVal
+    if(pPupId == false) continue;
+    // apply CHS in puppi metric computation -> use only LV hadrons to compute the metric for each particle
+    if(puppiAlgo_.at(iPuppiAlgo).fUseCharged_)  
+       pVal = goodVar(particlesAll[iPart], chargedPV,    puppiAlgo_.at(iPuppiAlgo).fMetricId_,puppiAlgo_.at(iPuppiAlgo).fConeSize_);
+    else if(!puppiAlgo_.at(iPuppiAlgo).fUseCharged_) 
+       pVal = goodVar(particlesAll[iPart], particlesAll, puppiAlgo_.at(iPuppiAlgo).fMetricId_,puppiAlgo_.at(iPuppiAlgo).fConeSize_);
+
+    // fill the value
+    if(std::isnan(pVal) || std::isinf(pVal)) std::cout << "====>  Value is Nan " << pVal << " == " << particlesAll[iPart].pt() << " -- " << particlesAll[iPart].eta() << std::endl;
+    if(std::isnan(pVal) || std::isinf(pVal)) continue;
+    
+    puppiParticles.push_back(puppiParticle(particlesAll.at(iPart).pt(),particlesAll.at(iPart).eta(),pVal,particlesAll.at(iPart).user_index(),iPart));
+  }
+  
+  // set the puppi particles for the algorithm
+  puppiAlgo_.at(iPuppiAlgo).setPuppiParticles(puppiParticles);
+  // compute RMS, median and mean value  
+  computeMedRMS(iPuppiAlgo);
+  
+}
+
+
+float puppiCleanContainer::goodVar(const fastjet::PseudoJet & particle, const std::vector<fastjet::PseudoJet> & particleAll, const int & pPupId, const float & coneSize) {
+  float lPup = 0;
+  lPup = var_within_R(pPupId,particleAll,particle,coneSize);
+  return lPup;
+}
+
+float puppiCleanContainer::var_within_R(const int & pPupId, const vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, const float & R){
+
+  if(pPupId == -1) return 1;
+  fastjet::Selector sel = fastjet::SelectorCircle(R);
+  sel.set_reference(centre);
+  std::vector<fastjet::PseudoJet> near_particles = sel(particles);
+  float var = 0;
+
+  for(size_t iPart = 0; iPart < near_particles.size(); iPart++){
+
+    double pDEta = near_particles[iPart].eta()-centre.eta();
+    double pDPhi = fabs(near_particles[iPart].phi()-centre.phi());
+    if(pDPhi > 2.*3.14159265-pDPhi) pDPhi = 2.*3.14159265-pDPhi;
+    double pDR = sqrt(pDEta*pDEta+pDPhi*pDPhi);
+
+    if(pDR < 0.0001) continue;
+    if(pDR == 0)    continue;
+
+    if(pPupId == 0) var += (near_particles[iPart].pt()/(pDR*pDR));
+    if(pPupId == 1) var += near_particles[iPart].pt();
+    if(pPupId == 2) var += (1./pDR)*(1./pDR);
+    if(pPupId == 3) var += (1./pDR)*(1./pDR);
+    if(pPupId == 4) var += near_particles[iPart].pt();
+    if(pPupId == 5) var += (near_particles[iPart].pt()/pDR)*(near_particles[iPart].pt()/pDR);
+  }
+
+  if(pPupId == 0 && var != 0) var = log(var);
+  if(pPupId == 3 && var != 0) var = log(var);
+  if(pPupId == 5 && var != 0) var = log(var);
+  return var;
+ 
+}
+
+
+float puppiCleanContainer::pt_within_R(const std::vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet & centre, const float & R){
+
+  fastjet::Selector sel = fastjet::SelectorCircle(R);
+  sel.set_reference(centre);
+  std::vector<fastjet::PseudoJet> near_particles = sel(particles);
+  double answer = 0.0;
+  for(size_t iPart = 0; iPart<near_particles.size(); iPart++){
+    answer += near_particles[iPart].pt();
+  }
+
+  return answer;
+}
+
+fastjet::PseudoJet puppiCleanContainer::flow_within_R(const vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, const float & R){
+
+  fastjet::Selector sel = fastjet::SelectorCircle(R);
+  sel.set_reference(centre);
+  std::vector<fastjet::PseudoJet> near_particles = sel(particles);
+  fastjet::PseudoJet flow;
+  for(unsigned int i=0; i<near_particles.size(); i++){
+    flow += near_particles[i];
+  }
+  return flow;
+
+}
+
+// compute median, mean value and RMS for puppi
+void puppiCleanContainer::computeMedRMS(const int & puppiAlgo) {
+
+  if(puppiAlgo > int(puppiAlgo_.size())  ) return;
+  if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size() == 0) return;
+
+  // sort in pVal increasing order
+  std::sort(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.begin(),puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.end(),puppiValSort());
+
+  // if apply correction
+  float lCorr = 1.;
+  if(puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_) lCorr *= 1.-fPVFrac_;
+
+  // count the position of the last particle with pval zero coming from PU
+  int lNum0 = 0;
+  for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size(); i0++) {
+    if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_[i0].fPval_ == 0) lNum0 = i0;
+  }
+
+  // take the median value on PU particles
+  int lNHalfway = lNum0 + int(float(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size()-lNum0)*0.50*lCorr);
+  puppiAlgo_.at(puppiAlgo).fMedian_ = puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(lNHalfway).fPval_;
+  float lMed = puppiAlgo_.at(puppiAlgo).fMedian_; //Just to make the readability easier
+
+  // take the RMS
+  int lNRMS = 0;
+  for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size(); i0++) {
+    puppiAlgo_.at(puppiAlgo).fMean_ += puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_;
+    if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_ == 0) continue;
+    if(!puppiAlgo_.at(puppiAlgo).fUseCharged_ && puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_ && puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_ > lMed) continue;
+    lNRMS++;
+    puppiAlgo_.at(puppiAlgo).fRMS_ += (puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_-lMed)*( puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_-lMed);
+  }
+
+  puppiAlgo_.at(puppiAlgo).fMean_ /= puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size();
+  if(lNRMS > 0) puppiAlgo_.at(puppiAlgo).fRMS_/=lNRMS;
+  if(puppiAlgo_.at(puppiAlgo).fRMS_ == 0) puppiAlgo_.at(puppiAlgo).fRMS_ = 1e-5;
+  puppiAlgo_.at(puppiAlgo).fRMS_  = sqrt(puppiAlgo_.at(puppiAlgo).fRMS_);
+  puppiAlgo_.at(puppiAlgo).fRMS_ *= puppiAlgo_.at(puppiAlgo).fRMSScaleFactor_;
+
+  if(!puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_) return;
+
+  //Adjust the p-value to correspond to the median
+  std::sort(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.begin(),puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.end(),puppiValSort());
+  int lNPV = 0; 
+  for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.size(); i0++){ 
+    if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_[i0].fPval_ <= lMed ) lNPV++;
+  }
+
+  // it helps in puppi the median value close to the mean one when a lot of pval 0 are present
+  float lAdjust = 1.5*float(lNPV)/float(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.size()+puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size());
+  if(lAdjust > 0) puppiAlgo_.at(puppiAlgo).fMedian_ -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*puppiAlgo_.at(puppiAlgo).fRMS_);
+
+}
+
+float puppiCleanContainer::getNeutralPtCut(const float & fNeutralMinE, const float & fNeutralPtSlope, const int & fNPV) {
+  return fNeutralMinE + fNPV * fNeutralPtSlope;
+}
+
+// take the type of algorithm : return a vector since more than one algo can be defined for the same eta region
+std::vector<int> puppiCleanContainer::getPuppiId(const float & pt, const float & eta, const std::vector<puppiAlgoBin> & puppiAlgos){
+  std::vector<int> PuppiId ;
+  for(size_t iPuppiAlgo = 0; iPuppiAlgo < puppiAlgos.size() ; iPuppiAlgo++){
+    if(fabs(eta) <= puppiAlgos[iPuppiAlgo].fEtaMin_) continue;
+    if(fabs(eta) > puppiAlgos[iPuppiAlgo].fEtaMax_) continue;
+    PuppiId.push_back(int(iPuppiAlgo));
+  }
+  return PuppiId;  
+}
+
+//check if a particle is good for an Algo definition
+bool puppiCleanContainer::isGoodPuppiId(const float & pt, const float & eta, const puppiAlgoBin & puppiAlgo){
+  if(fabs(eta) <= puppiAlgo.fEtaMin_) return false;
+  if(fabs(eta) > puppiAlgo.fEtaMax_) return false;
+  return true;
+
+}
+
+
+
+// ----------------------
+float puppiCleanContainer::compute(const float & chi2, const std::vector<puppiParticle> & particles, const std::vector<puppiAlgoBin> & puppiAlgos, const std::vector<int> & pPupId) {
+
+  if(particles.size() != pPupId.size() ) return 0; // default check
+
+  float lVal  = 0.;
+  float lPVal = 1.;
+  int   lNDOF = 0;
+
+  for( size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){
+
+    if(puppiAlgos.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.size() + puppiAlgos.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.size()  == 0) return 1; 
+
+    if(iAlgo > 0 ){
+      float pPVal = ROOT::Math::chisquared_cdf(lVal,lNDOF); // take a chi2 value since the should be multiplied (multiply weight and not summing chi2)
+      lPVal *= pPVal;
+      lNDOF  = 0;
+      lVal   = 0; 
+    }
+
+    if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == -1) continue;
+
+    float pVal = particles.at(iAlgo).fPval_ ;
+
+    if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 0 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
+    if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 3 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
+    if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 5 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
+    
+    lVal += (pVal-puppiAlgos.at(pPupId.at(iAlgo)).fMedian_)*(fabs(pVal-puppiAlgos.at(pPupId.at(iAlgo)).fMedian_))/puppiAlgos.at(pPupId.at(iAlgo)).fRMS_/puppiAlgos.at(pPupId.at(iAlgo)).fRMS_;
+    lNDOF++;
+    if(chi2 != 0) lNDOF++; 
+    if(chi2 != 0) lVal+=chi2; //Add external Chi2 to first element
+  }
+
+  lPVal *= ROOT::Math::chisquared_cdf(lVal,lNDOF);
+  return lPVal;
+
+}
+
+float puppiCleanContainer::getChi2FromdZ(float iDZ) {
+   //We need to obtain prob of PU + (1-Prob of LV)
+   // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm (its really more like 1mm)
+   //double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),0.2)*2.; //*2 is to do it double sided
+   //Take iDZ to be corrected by sigma already
+   double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),1.)*2.; //*2 is to do it double sided
+   double lProbPU = 1-lProbLV;
+   if(lProbPU <= 0) lProbPU = 1e-16; //Quick Trick to through out infs
+   if(lProbPU >= 0) lProbPU = 1-1e-16; //Ditto
+   double lChi2PU = TMath::ChisquareQuantile(lProbPU,1);
+   lChi2PU*=lChi2PU;
+   return lChi2PU;
+}
Index: external/PUPPI/puppiCleanContainer.hh
===================================================================
--- external/PUPPI/puppiCleanContainer.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
+++ external/PUPPI/puppiCleanContainer.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -0,0 +1,83 @@
+#ifndef PUPPICLEANCONTAINER_HH
+#define PUPPICLEANCONTAINER_HH
+
+
+#include "PUPPI/RecoObj.hh"
+#include "PUPPI/puppiParticle.hh"
+#include "PUPPI/puppiAlgoBin.hh"
+
+#include "fastjet/internal/base.hh"
+#include "fastjet/PseudoJet.hh"
+#include <algorithm>
+
+using namespace std;
+
+//......................
+
+class puppiCleanContainer{
+
+ public:
+
+  // basic constructor which takes input particles as RecoObj, the tracker eta extension and other two boolean info
+  puppiCleanContainer(std::vector<RecoObj> inParticles,    // incoming particles of the event
+                      std::vector<puppiAlgoBin> puppiAlgo, // vector with the definition of the puppi algorithm in different eta region (one for each eta) 
+                      float minPuppiWeight  = 0.01,        // min puppi weight cut
+                      bool  useExp = false                 // useDz vertex probability
+  ); 
+
+  ~puppiCleanContainer(); 
+
+  // ----- get methods
+
+  // get all the PF particles
+  std::vector<fastjet::PseudoJet> pfParticles()  { return  fPFParticles_; }    
+  // get all the PF charged from PV
+  std::vector<fastjet::PseudoJet> pvParticles()  { return  fChargedPV_; }        
+  // get all the PF charged from PU
+  std::vector<fastjet::PseudoJet> puParticles()  { return  fChargedNoPV_; }    
+  // get CHS particle collection
+  std::vector<fastjet::PseudoJet> pfchsParticles(){ return fPFchsParticles_; }    
+  // get puppi weight for all particles 
+  std::vector<float> getPuppiWeights() { return fPuppiWeights_; };
+
+  // process puppi
+  std::vector<fastjet::PseudoJet> puppiEvent();
+ 
+ protected:
+
+   void    getRMSAvg(const int &, std::vector<fastjet::PseudoJet> &, std::vector<fastjet::PseudoJet> &);        
+   float   goodVar  (const fastjet::PseudoJet &, const std::vector<fastjet::PseudoJet> &, const int &, const float &);    
+   void    computeMedRMS(const int &);  
+   float   compute(const float &, const std::vector<puppiParticle> &, const std::vector<puppiAlgoBin> &, const std::vector<int> &);
+
+   // some get functions
+   float getNeutralPtCut(const float&, const float&, const int&);
+   std::vector<int> getPuppiId(const float &, const float &, const std::vector<puppiAlgoBin> &);
+   bool  isGoodPuppiId(const float &, const float &, const puppiAlgoBin &);
+   float getChi2FromdZ(float);
+   // other functions
+   float  var_within_R(const int &, const vector<fastjet::PseudoJet> &, const fastjet::PseudoJet &, const float &);
+   float  pt_within_R(const std::vector<fastjet::PseudoJet> &, const fastjet::PseudoJet &, const float &);
+   fastjet::PseudoJet flow_within_R(const vector<fastjet::PseudoJet> &, const fastjet::PseudoJet &, const float &);
+   
+  
+ private:    
+ 
+  std::vector<RecoObj>            fRecoParticles_;
+  std::vector<fastjet::PseudoJet> fPFParticles_;
+  std::vector<fastjet::PseudoJet> fPFchsParticles_;    
+  std::vector<fastjet::PseudoJet> fChargedPV_;
+  std::vector<fastjet::PseudoJet> fChargedNoPV_;
+
+  std::vector<puppiAlgoBin> puppiAlgo_;
+  std::vector<float> fPuppiWeights_;
+
+  float  fMinPuppiWeight_;
+  float  fPVFrac_;
+
+  int    fNPV_;  
+  bool   fUseExp_ ;
+    
+};
+
+#endif
Index: external/PUPPI/puppiParticle.hh
===================================================================
--- external/PUPPI/puppiParticle.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
+++ external/PUPPI/puppiParticle.hh	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -0,0 +1,76 @@
+#ifndef PUPPIPARTICLE_HH
+#define PUPPIPARTICLE_HH
+
+#include <functional>
+#include <algorithm>
+
+
+//............ class to keep track of particles used in the puppi code storing pt,eta,metric value, id (LV,PU) and position in the orignal vector of particles 
+class puppiParticle {
+
+ public:
+
+  puppiParticle(){
+    fPt_         = 0.;
+    fEta_        = 0.;
+    fPval_       = 1.;
+    fParticleId_ = -1;
+    fPosition_   = -1; 
+  }; 
+
+  puppiParticle(const float & pt, const float & eta, const float & pval, const int & particleId, const int & position):
+    fPt_(pt),
+    fEta_(eta),
+    fPval_(pval),
+    fParticleId_(particleId),
+    fPosition_(position)
+  {};
+
+  puppiParticle(const puppiParticle & particle) {
+    fPt_         = particle.fPt_;
+    fEta_        = particle.fEta_;
+    fPval_       = particle.fPval_;
+    fParticleId_ = particle.fParticleId_;
+    fPosition_   = particle.fPosition_;
+
+  }
+
+  virtual ~puppiParticle(){};
+
+  // sort asaf of pt
+  bool operator < (const puppiParticle & particle2) const {
+    if(fPt_ < particle2.fPt_) return true;
+    else return false;
+  }; 
+
+
+  // oper operator =
+  bool operator == (const puppiParticle & particle2) const {
+    if( fPt_ == particle2.fPt_     and fEta_ == particle2.fEta_ and 
+        fPval_ == particle2.fPval_ and fParticleId_ == particle2.fParticleId_ and 
+        fPosition_ == particle2.fPosition_) return true;
+    else return false ;
+  };
+
+  float fPt_;         // pt of the candidate                                                                                                                                              
+  float fEta_;        // eta of the candidate                                                                                                                                             
+  float fPval_;       // metric value                                                                                                                                                     
+  int   fParticleId_; // particle id means user_index                                                                                                                                     
+  int   fPosition_;   // position in the original particle vector                                                                                                                        
+
+};
+
+class puppiValSort : public std::binary_function<int,int,bool> {
+ public:
+
+  puppiValSort(){};
+
+  ~puppiValSort(){};
+
+  bool operator() (const puppiParticle & x, const puppiParticle & y){
+    if(x.fPval_ < y.fPval_ ) return true;
+    return false;
+  }
+};
+
+#endif
Index: modules/RunPUPPI.cc
===================================================================
--- modules/RunPUPPI.cc	(revision 7b45ff5e3840b1a7e6b83f7376cdb74e900f182c)
+++ modules/RunPUPPI.cc	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -1,8 +1,8 @@
 #include "modules/RunPUPPI.h"
 
-#include "PUPPI/RecoObj2.hh"
-#include "PUPPI/AlgoObj.hh"
-//#include "PUPPI/puppiParticle.hh"
-//#include "PUPPI/puppiAlgoBin.hh"
+#include "PUPPI/puppiCleanContainer.hh"
+#include "PUPPI/RecoObj.hh"
+#include "PUPPI/puppiParticle.hh"
+#include "PUPPI/puppiAlgoBin.hh"
 
 #include "classes/DelphesClasses.h"
@@ -30,4 +30,5 @@
 
 void RunPUPPI::Init(){
+
   // input collection
   fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
@@ -37,101 +38,74 @@
   fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
   fPVItInputArray      = fPVInputArray->MakeIterator();
-  // puppi parameters                                 
-  fApplyNoLep     = GetBool("UseNoLep", true);    
+
+
+  // puppi parameters                                     
   fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
   fUseExp         = GetBool("UseExp", false);
-  // read eta min ranges
+
+  // 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
+
+  // 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
+
+  // 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
+
+  // 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
+
+  // 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
+
+  // 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("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
+
+  // 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());
-  // scheme for combining
-  param = GetParam("CombId");
-  fCombId.clear();
-  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
+
   // create output array
   fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
   fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
   fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
-  // Create algorithm list for puppi
-  std::vector<AlgoObj> 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() == 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++){
-    AlgoObj algoTmp ;
-    algoTmp.etaMin            = fEtaMinBin.at(iAlgo);
-    algoTmp.etaMax            = fEtaMaxBin.at(iAlgo);
-    algoTmp.ptMin             = fPtMinBin.at(iAlgo);
-    algoTmp.minNeutralPt      = fNeutralMinEBin.at(iAlgo);
-    algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
-    //Eta Extrapolation stuff is missing
-    //Loop through file requiring algos for same bins to be adjacent
-    while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) { 
-      AlgoSubObj algoSubTmp;
-      algoSubTmp.metricId          = fMetricId.at(iAlgo);
-      algoSubTmp.useCharged        = fUseCharged.at(iAlgo);
-      algoSubTmp.applyLowPUCorr    = fApplyLowPUCorr.at(iAlgo);
-      algoSubTmp.combId            = fCombId.at(iAlgo);
-      algoSubTmp.coneSize          = fConeSizeBin.at(iAlgo);
-      algoSubTmp.rmsPtMin          = fRMSPtMinBin.at(iAlgo);
-      algoSubTmp.rmsScaleFactor    = fRMSScaleFactorBin.at(iAlgo);
-      algoTmp.subAlgos.push_back(algoSubTmp);
-      iAlgo++;
-    }
-    iAlgo--;
-    //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;    
-    puppiAlgo.push_back(algoTmp);     
-  }
-  fPuppi  = new PuppiContainer(true,fUseExp,fMinPuppiWeight,puppiAlgo);
 }
 
@@ -150,10 +124,10 @@
   TLorentzVector momentum;
 
-  //DelphesFactory *factory = GetFactory();
+  DelphesFactory *factory = GetFactory();
 
   // loop over input objects
-  fItTrackInputArray   ->Reset();
-  fItNeutralInputArray ->Reset();
-  fPVItInputArray      ->Reset();
+  fItTrackInputArray->Reset();
+  fItNeutralInputArray->Reset();
+  fPVItInputArray->Reset();
 
   std::vector<Candidate *> InputParticles;
@@ -168,8 +142,10 @@
   std::vector<RecoObj> puppiInputVector;
   puppiInputVector.clear();
-  int lNBad  = 0; 
+
   // Loop on charge track candidate
   while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
+
       momentum = candidate->Momentum;
+      
       RecoObj curRecoObj;
       curRecoObj.pt  = momentum.Pt();
@@ -178,10 +154,7 @@
       curRecoObj.m   = momentum.M();  
       particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
-      //if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
-      //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
       if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
-	lNBad++;
 	curRecoObj.id    = 2;
-	curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
+	curRecoObj.vtxId = candidate->IsPU;
 	if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
 	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
@@ -210,5 +183,7 @@
   // Loop on neutral calo cells 
   while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
+
       momentum = candidate->Momentum;
+
       RecoObj curRecoObj;
       curRecoObj.pt  = momentum.Pt();
@@ -216,6 +191,6 @@
       curRecoObj.phi = momentum.Phi();
       curRecoObj.m   = momentum.M();
-      curRecoObj.charge = 0;
       particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
+
 
       if(candidate->Charge == 0){
@@ -235,8 +210,38 @@
       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
-  fPuppi->initialize(puppiInputVector);
-  fPuppi->puppiWeights();
-  std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
+  puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
+  std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
 
   // Loop on final particles
Index: modules/RunPUPPI.h
===================================================================
--- modules/RunPUPPI.h	(revision 7b45ff5e3840b1a7e6b83f7376cdb74e900f182c)
+++ modules/RunPUPPI.h	(revision 0e11b5ca37ecb51f9641a37e627795b31038660f)
@@ -3,5 +3,4 @@
 
 #include "classes/DelphesModule.h"
-#include "PUPPI/PuppiContainer.hh"
 #include <vector>
 
@@ -30,8 +29,7 @@
   const TObjArray *fNeutralInputArray; //!
   const TObjArray *fPVInputArray; //!                                                                                                                                                     
-  PuppiContainer* fPuppi;
+ 
   // puppi parameters
-  bool fApplyNoLep;
-  double fMinPuppiWeight;
+  float fMinPuppiWeight;
   bool fUseExp;
   
@@ -48,5 +46,4 @@
   std::vector<bool>  fApplyLowPUCorr;
   std::vector<int>   fMetricId;
-  std::vector<int>   fCombId;
 
   TObjArray *fOutputArray;
