Index: external/PUPPI/AlgoObj.hh
===================================================================
--- external/PUPPI/AlgoObj.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
+++ external/PUPPI/AlgoObj.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
@@ -0,0 +1,36 @@
+#ifndef ALGOOBJ_HH
+#define ALGOOBJ_HH
+#include <vector>
+
+class AlgoSubObj
+{
+public:
+  int    metricId;
+  bool   useCharged;
+  bool   applyLowPUCorr;
+  int    combId;
+  double coneSize;
+  double rmsPtMin;
+  double rmsScaleFactor;
+};
+class AlgoObj 
+{
+public:
+      AlgoObj():
+	etaMin(0), etaMax(0), 
+	ptMin(0),minNeutralPt(0),minNeutralPtSlope(0),
+	rmsEtaSF(0),medEtaSF(0),etaMaxExtrap(0)
+    {}
+    ~AlgoObj(){}
+  
+    float etaMin;
+    float etaMax;
+    float ptMin;
+    double minNeutralPt;
+    double minNeutralPtSlope;
+    double rmsEtaSF;
+    double medEtaSF;
+    double etaMaxExtrap;  
+    std::vector<AlgoSubObj> subAlgos; 
+};
+#endif
Index: external/PUPPI/PuppiAlgo.cc
===================================================================
--- external/PUPPI/PuppiAlgo.cc	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
+++ external/PUPPI/PuppiAlgo.cc	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
@@ -0,0 +1,190 @@
+//#include "FWCore/ParameterSet/interface/ParameterSet.h"
+//#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
+#include "PuppiAlgo.hh"
+#include "external/fastjet/internal/base.hh"
+#include "Math/QuantFuncMathCore.h"
+#include "Math/SpecFuncMathCore.h"
+#include "Math/ProbFunc.h"
+#include "TMath.h"
+
+PuppiAlgo::PuppiAlgo(AlgoObj &iAlgo) { 
+  fEtaMin             = iAlgo.etaMin;
+  fEtaMax             = iAlgo.etaMax;
+  fPtMin              = iAlgo.ptMin;
+  fNeutralPtMin       = iAlgo.minNeutralPt;
+  fNeutralPtSlope     = iAlgo.minNeutralPtSlope;
+  fRMSEtaSF           = iAlgo.rmsEtaSF;
+  fMedEtaSF           = iAlgo.medEtaSF;
+  fEtaMaxExtrap       = iAlgo.etaMaxExtrap;
+  std::vector<AlgoSubObj> lAlgos = iAlgo.subAlgos;
+  fNAlgos = lAlgos.size();
+  std::cout << "==> "  << fEtaMin << " - " << fEtaMax << " - " << fPtMin  << " - " << fNeutralPtMin  << " - " << fNeutralPtSlope  << " - " << fRMSEtaSF  << " - " << std::endl;
+  for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++)  { 
+    int    pAlgoId      = lAlgos[i0].metricId;
+    bool   pCharged     = lAlgos[i0].useCharged;
+    bool   pWeight0     = lAlgos[i0].applyLowPUCorr;
+    int    pComb        = lAlgos[i0].combId;
+    double pConeSize    = lAlgos[i0].coneSize;
+    double pRMSPtMin    = lAlgos[i0].rmsPtMin;
+    double pRMSSF       = lAlgos[i0].rmsScaleFactor;
+    std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
+    fAlgoId        .push_back(pAlgoId);
+    fCharged       .push_back(pCharged);
+    fAdjust        .push_back(pWeight0);
+    fCombId        .push_back(pComb);
+    fConeSize      .push_back(pConeSize);
+    fRMSPtMin      .push_back(pRMSPtMin);
+    fRMSScaleFactor.push_back(pRMSSF);
+    double pRMS  = 0; 
+    double pMed  = 0; 
+    double pMean = 0;
+    int    pNCount = 0; 
+    fRMS   .push_back(pRMS);
+    fMedian.push_back(pMed);
+    fMean  .push_back(pMean);
+    fNCount.push_back(pNCount);
+    /*
+      tmprms .clear();
+      tmpmed .clear();
+      for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++){
+      tmprms.push_back(pRMS);
+      tmpmed.push_back(pMed);
+      }
+      fRMS_perEta.push_back(tmprms);
+      fMedian_perEta.push_back(tmpmed);
+    */
+  }
+}
+PuppiAlgo::~PuppiAlgo() { 
+  fPups  .clear();
+  fPupsPV.clear();
+}
+void PuppiAlgo::reset() { 
+  fPups  .clear();
+  fPupsPV.clear();
+  for(unsigned int i0 = 0; i0 < fNAlgos; i0++) { 
+    fMedian[i0] =  0; 
+    fRMS   [i0] =  0;
+    fMean  [i0] =  0;
+    fNCount[i0] =  0;
+  }
+}
+void PuppiAlgo::add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo) { 
+  if(iParticle.pt() < fRMSPtMin[iAlgo]) return;
+  if(fCharged[iAlgo] && fabs(iParticle.user_index())  < 1) return;
+  if(fCharged[iAlgo] && (fabs(iParticle.user_index()) >=1 && fabs(iParticle.user_index()) <=2)) fPupsPV.push_back(iVal);
+  if(fCharged[iAlgo] && fabs(iParticle.user_index()) < 3) return;
+  fPups.push_back(iVal);
+  fNCount[iAlgo]++;
+  /* eta extrap (cmssw version)
+     if ((std::abs(iParticle.eta()) < fEtaMaxExtrap) && (std::abs(puppi_register) >= 3)){
+        fPups.push_back(iVal);
+        // fPupsPV.push_back(iVal);        
+        fNCount[iAlgo]++;
+	}
+    // for the low PU case, correction.  for checking that the PU-only median will be below the PV particles
+    if(std::abs(iParticle.eta()) < fEtaMaxExtrap && (std::abs(puppi_register) >=1 && std::abs(puppi_register) <=2)) fPupsPV.push_back(iVal);
+  */
+}
+void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) { 
+  if(iAlgo >= fNAlgos   ) return;
+  if(fNCount[iAlgo] == 0) return;
+  int lNBefore = 0; 
+  for(unsigned int i0 = 0; i0 < iAlgo; i0++) lNBefore += fNCount[i0];
+  std::sort(fPups.begin()+lNBefore,fPups.begin()+lNBefore+fNCount[iAlgo]);
+  int lNum0 = 0;
+  for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) { 
+    if(fPups[i0] == 0) lNum0 = i0-lNBefore; 
+  }
+  //lNum0 = 0; 
+  int lNHalfway = lNBefore + lNum0 + int( double( fNCount[iAlgo]-lNum0 )*0.50);
+  fMedian[iAlgo] = fPups[lNHalfway];
+  double lMed = fMedian[iAlgo];  //Just to make the readability easier
+
+  int lNRMS = 0; 
+  for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) {
+    fMean[iAlgo] += fPups[i0];
+    if(fPups[i0] == 0) continue;
+    if(fAdjust[iAlgo] && fPups[i0] > lMed) continue;
+    lNRMS++;
+    fRMS [iAlgo] += (fPups[i0]-lMed)*(fPups[i0]-lMed);
+  }
+  fMean[iAlgo]/=fNCount[iAlgo];
+  if(lNRMS > 0) fRMS [iAlgo]/=lNRMS;
+  if(fRMS[iAlgo] == 0) fRMS[iAlgo] = 1e-5;
+
+  fRMS [iAlgo] = sqrt(fRMS[iAlgo]);
+  fRMS [iAlgo] *= fRMSScaleFactor[iAlgo];
+  if(fAdjust[iAlgo]){ 
+    //Adjust the p-value to correspond to the median
+    std::sort(fPupsPV.begin(),fPupsPV.end());
+    int lNPV = 0; 
+    for(unsigned int i0 = 0; i0 < fPupsPV.size(); i0++) if(fPupsPV[i0] <= lMed ) lNPV++;
+    double lAdjust = double(lNPV)/double(lNPV+0.5*fNCount[iAlgo]);
+    if(lAdjust > 0) {
+      fMedian[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*fRMS[iAlgo]);
+      fRMS[iAlgo]    -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*fRMS[iAlgo]);
+    }        
+  }
+  /* eta extrapolated version
+  for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++){
+    fRMS_perEta[iAlgo][j0]    = fRMS[iAlgo]*fRMSEtaSF[j0];
+    fMedian_perEta[iAlgo][j0] = fMedian[iAlgo]*fMedEtaSF[j0];
+  }    
+  */
+}
+//This code is probably a bit confusing
+double PuppiAlgo::compute(std::vector<double> &iVals,double iChi2) { 
+  if(fAlgoId[0] == -1) return 1;
+  double lVal  = 0.;
+  double lPVal = 1.;
+  int    lNDOF = 0; 
+  for(unsigned int i0 = 0; i0 < fNAlgos; i0++) { 
+    if(fNCount[i0] == 0) return 1.;   //in the NoPU case return 1.
+    if(fCombId[i0] == 1 && i0 > 0) {  //Compute the previous p-value so that p-values can be multiplieed
+      double pPVal = ROOT::Math::chisquared_cdf(lVal,lNDOF);
+      lPVal *= pPVal;
+      lNDOF = 0; 
+      lVal  = 0; 
+    }
+    double pVal = iVals[i0];
+    //Special Check for any algo with log(0) 
+    if(fAlgoId[i0] == 0 && iVals[i0] == 0) pVal = fMedian[i0];
+    if(fAlgoId[i0] == 3 && iVals[i0] == 0) pVal = fMedian[i0];
+    if(fAlgoId[i0] == 5 && iVals[i0] == 0) pVal = fMedian[i0];
+    lVal += (pVal-fMedian[i0])*(fabs(pVal-fMedian[i0]))/fRMS[i0]/fRMS[i0];
+    lNDOF++;
+    if(i0 == 0 && iChi2 != 0) lNDOF++;      //Add external Chi2 to first element
+    if(i0 == 0 && iChi2 != 0) lVal+=iChi2;  //Add external Chi2 to first element
+  }
+  //Top it off with the last calc
+  lPVal *= ROOT::Math::chisquared_cdf(lVal,lNDOF);
+  return lPVal;
+}
+double PuppiAlgo::neutralPt(int iNPV) { 
+  return fNeutralPtMin + iNPV * fNeutralPtSlope;
+}
+int PuppiAlgo::numAlgos() { 
+  return fNAlgos;
+}
+double PuppiAlgo::ptMin() { 
+  return fPtMin;
+}
+double PuppiAlgo::etaMin() { 
+  return fEtaMin;
+}
+double PuppiAlgo::etaMax() { 
+  return fEtaMax;
+}
+int PuppiAlgo::algoId(const unsigned int &iAlgo) { 
+  assert(iAlgo < fNAlgos);
+  return fAlgoId[iAlgo];
+}
+bool PuppiAlgo::isCharged(const unsigned int &iAlgo) { 
+  assert(iAlgo < fNAlgos);
+  return fCharged[iAlgo];
+}
+double PuppiAlgo::coneSize(const unsigned int &iAlgo) { 
+  assert(iAlgo < fNAlgos);
+  return fConeSize[iAlgo];
+}
Index: external/PUPPI/PuppiAlgo.hh
===================================================================
--- external/PUPPI/PuppiAlgo.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
+++ external/PUPPI/PuppiAlgo.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
@@ -0,0 +1,49 @@
+#include "external/fastjet/internal/base.hh"
+#include "external/fastjet/PseudoJet.hh"
+#include "AlgoObj.hh"
+#include <vector>
+
+class PuppiAlgo{ 
+public:
+  PuppiAlgo(AlgoObj &iAlgoObj);
+  ~PuppiAlgo();
+  //Computing Mean and RMS
+  void   reset();
+  void   add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo);
+  void   computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac);
+  //Get the Weight
+  double compute(std::vector<double> &iVals,double iChi2);
+  //Helpers
+  double ptMin();
+  double etaMin();
+  double etaMax();
+  int    numAlgos ();
+  int    algoId   (const unsigned int &iAlgo);
+  bool   isCharged(const unsigned int &iAlgo);
+  double coneSize (const unsigned int &iAlgo);
+  double neutralPt(int iNPV);
+
+private:  
+  unsigned int   fNAlgos;
+  float  fEtaMax;
+  float  fEtaMin;
+  float  fPtMin ;
+  double fNeutralPtMin;
+  double fNeutralPtSlope;
+  double fRMSEtaSF;
+  double fMedEtaSF;
+  double fEtaMaxExtrap;
+  std::vector<float>  fPups;
+  std::vector<float>  fPupsPV;
+  std::vector<int>    fAlgoId;
+  std::vector<bool>   fCharged;
+  std::vector<bool>   fAdjust;
+  std::vector<int>    fCombId;
+  std::vector<double> fConeSize;
+  std::vector<double> fRMSPtMin;
+  std::vector<double> fRMSScaleFactor;
+  std::vector<double> fRMS;
+  std::vector<double> fMedian;
+  std::vector<double> fMean;
+  std::vector<int>    fNCount;
+};
Index: external/PUPPI/RecoObj2.hh
===================================================================
--- external/PUPPI/RecoObj2.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
+++ external/PUPPI/RecoObj2.hh	(revision 07b7064a30564da84ad3a9b5fe62e48301005d6d)
@@ -0,0 +1,34 @@
+#ifndef RECOOBJ_HH
+#define 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;
+    int           charge;
+};
+#endif
+
+// Only need 4-vector and ID
+// id = 0:   neutral
+// id = 1; charged from PV
+// id = 2; charged not from PV
