Skip to content

Commit

Permalink
rebase to head -- Merge branch 'CMSSW_8_0_X' into from-CMSSW_8_0_X_20…
Browse files Browse the repository at this point in the history 
…15-10-22-2300_fixEraTypo
  • Loading branch information
@hengne
hengne committed Oct 23, 2015
2 parents ef84225 + 8b1004f commit 3f23db9
Show file tree
Hide file tree
Showing 6 changed files with 148 additions and 87 deletions.
44 changes: 28 additions & 16 deletions CommonTools/PileupAlgos/interface/PuppiAlgo.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,38 +12,49 @@ class PuppiAlgo{
~PuppiAlgo();
//Computing Mean and RMS
void reset();
void fixAlgoEtaBin(int i_eta );
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> const &iVals,double iChi2) const;
const std::vector<float> & alphas(){ return fPups; }
//Helpers
inline double ptMin() const { return fPtMin; }
inline double etaMin() const { return fEtaMin; }
inline double etaMax() const { return fEtaMax; }
inline int etaBins() const {return fEtaMin.size(); }
inline double etaMin(int i) const { return fEtaMin[i]; }
inline double etaMax(int i) const { return fEtaMax[i]; }
inline double ptMin() const { return cur_PtMin; }

inline int numAlgos () const { return fNAlgos;}
inline int algoId ( unsigned int iAlgo) const { return fAlgoId.at(iAlgo); }
inline bool isCharged ( unsigned int iAlgo) const { return fCharged.at(iAlgo); }
inline double coneSize ( unsigned int iAlgo) const { return fConeSize.at(iAlgo); }
inline double neutralPt (int iNPV) const { return fNeutralPtMin + iNPV * fNeutralPtSlope; }
inline double neutralPt (int iNPV) const { return cur_NeutralPtMin + iNPV * cur_NeutralPtSlope; }

inline double rms( unsigned int i ) const {return fRMS[i];}
inline double median( unsigned int i ) const {return fMedian[i];}
inline double rms() const {return cur_RMS;}
inline double median() const {return cur_Med;}

private:
unsigned int fNAlgos;
float fEtaMax;
float fEtaMin;
float fPtMin ;
double fNeutralPtMin;
double fNeutralPtSlope;

double fRMSEtaSF;
double fMedEtaSF;
std::vector<double> fEtaMax;
std::vector<double> fEtaMin;
std::vector<double> fPtMin ;
std::vector<double> fNeutralPtMin;
std::vector<double> fNeutralPtSlope;

std::vector<double> fRMSEtaSF;
std::vector<double> fMedEtaSF;
double fEtaMaxExtrap;

std::vector<double> fRMS;
std::vector<double> fMedian;
double cur_PtMin;
double cur_NeutralPtMin;
double cur_NeutralPtSlope;
double cur_RMS;
double cur_Med;

std::vector<double> fRMS; // this is the raw RMS per algo
std::vector<double> fMedian; // this is the raw Median per algo
std::vector< std::vector<double> > fRMS_perEta; // this is the final RMS used after eta corrections
std::vector< std::vector<double> > fMedian_perEta; // this is the final Med used after eta corrections

std::vector<float> fPups;
std::vector<float> fPupsPV;
Expand All @@ -56,6 +67,7 @@ class PuppiAlgo{
std::vector<double> fRMSScaleFactor;
std::vector<double> fMean;
std::vector<int> fNCount;

};


Expand Down
58 changes: 28 additions & 30 deletions CommonTools/PileupAlgos/python/Puppi_cff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,41 +43,39 @@
vtxZCut = cms.double(24),
algos = cms.VPSet(
cms.PSet(
etaMin = cms.double(0.),
etaMax = cms.double(2.5),
ptMin = cms.double(0.),
MinNeutralPt = cms.double(0.1),
MinNeutralPtSlope = cms.double(0.015),
RMSEtaSF = cms.double(1.0),
MedEtaSF = cms.double(1.0),
etaMin = cms.vdouble(0.),
etaMax = cms.vdouble(2.5),
ptMin = cms.vdouble(0.),
MinNeutralPt = cms.vdouble(0.1),
MinNeutralPtSlope = cms.vdouble(0.015),
RMSEtaSF = cms.vdouble(1.0),
MedEtaSF = cms.vdouble(1.0),
EtaMaxExtrap = cms.double(2.0),
puppiAlgos = puppiCentral
),
cms.PSet(
etaMin = cms.double(2.5),
etaMax = cms.double(3.0),
ptMin = cms.double(0.0),
MinNeutralPt = cms.double(1.7),
MinNeutralPtSlope = cms.double(0.07),
# RMSEtaSF = cms.double(1.545),
# MedEtaSF = cms.double(0.845),
RMSEtaSF = cms.double(1.30),
MedEtaSF = cms.double(1.05),
EtaMaxExtrap = cms.double(2.0),
etaMin = cms.vdouble( 2.5, 3.0),
etaMax = cms.vdouble( 3.0, 10.0),
ptMin = cms.vdouble( 0.0, 0.0),
MinNeutralPt = cms.vdouble( 1.7, 2.0),
MinNeutralPtSlope = cms.vdouble(0.07, 0.07),
RMSEtaSF = cms.vdouble(1.30, 1.10),
MedEtaSF = cms.vdouble(1.05, 0.90),
EtaMaxExtrap = cms.double( 2.0),
puppiAlgos = puppiForward
),
cms.PSet(
etaMin = cms.double(3.0),
etaMax = cms.double(10.0),
ptMin = cms.double(0.0),
MinNeutralPt = cms.double(2.0),
MinNeutralPtSlope = cms.double(0.07),
# RMSEtaSF = cms.double(1.18),
# MedEtaSF = cms.double(0.4397),
RMSEtaSF = cms.double(1.10),
MedEtaSF = cms.double(0.90),
EtaMaxExtrap = cms.double(2.0),
puppiAlgos = puppiForward
)
# cms.PSet(
# etaMin = cms.double(3.0),
# etaMax = cms.double(10.0),
# ptMin = cms.double(0.0),
# MinNeutralPt = cms.double(2.0),
# MinNeutralPtSlope = cms.double(0.07),
# # RMSEtaSF = cms.double(1.18),
# # MedEtaSF = cms.double(0.4397),
# RMSEtaSF = cms.double(1.10),
# MedEtaSF = cms.double(0.90),
# EtaMaxExtrap = cms.double(2.0),
# puppiAlgos = puppiForward
# )
)
)
84 changes: 61 additions & 23 deletions CommonTools/PileupAlgos/src/PuppiAlgo.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,21 @@


PuppiAlgo::PuppiAlgo(edm::ParameterSet &iConfig) {
fEtaMin = iConfig.getParameter<double>("etaMin");
fEtaMax = iConfig.getParameter<double>("etaMax");
fPtMin = iConfig.getParameter<double>("ptMin");
fNeutralPtMin = iConfig.getParameter<double>("MinNeutralPt"); // Weighted Neutral Pt Cut
fNeutralPtSlope = iConfig.getParameter<double>("MinNeutralPtSlope"); // Slope vs #pv
fRMSEtaSF = iConfig.getParameter<double>("RMSEtaSF");
fMedEtaSF = iConfig.getParameter<double>("MedEtaSF");
fEtaMin = iConfig.getParameter<std::vector< double > >("etaMin");
fEtaMax = iConfig.getParameter<std::vector< double > >("etaMax");
fPtMin = iConfig.getParameter<std::vector< double > >("ptMin");
fNeutralPtMin = iConfig.getParameter<std::vector< double > >("MinNeutralPt"); // Weighted Neutral Pt Cut
fNeutralPtSlope = iConfig.getParameter<std::vector< double > >("MinNeutralPtSlope"); // Slope vs #pv
fRMSEtaSF = iConfig.getParameter<std::vector< double > >("RMSEtaSF");
fMedEtaSF = iConfig.getParameter<std::vector< double > >("MedEtaSF");
fEtaMaxExtrap = iConfig.getParameter<double>("EtaMaxExtrap");

std::vector<edm::ParameterSet> lAlgos = iConfig.getParameter<std::vector<edm::ParameterSet> >("puppiAlgos");
fNAlgos = lAlgos.size();
//Uber Configurable Puppi
std::vector<double> tmprms;
std::vector<double> tmpmed;

for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
int pAlgoId = lAlgos[i0].getParameter<int > ("algoId");
bool pCharged = lAlgos[i0].getParameter<bool> ("useCharged");
Expand All @@ -44,7 +47,22 @@ PuppiAlgo::PuppiAlgo(edm::ParameterSet &iConfig) {
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);
}

cur_PtMin = -99.;
cur_NeutralPtMin = -99.;
cur_NeutralPtSlope = -99.;
cur_RMS = -99.;
cur_Med = -99.;
}
PuppiAlgo::~PuppiAlgo() {
fPups .clear();
Expand All @@ -60,6 +78,15 @@ void PuppiAlgo::reset() {
fNCount[i0] = 0;
}
}

void PuppiAlgo::fixAlgoEtaBin(int i_eta) {
cur_PtMin = fPtMin[i_eta];
cur_NeutralPtMin = fNeutralPtMin[i_eta];
cur_NeutralPtSlope = fNeutralPtSlope[i_eta];
cur_RMS = fRMS_perEta[0][i_eta]; // 0 is number of algos within this eta bin
cur_Med = fMedian_perEta[0][i_eta]; // 0 is number of algos within this eta bin
}

void PuppiAlgo::add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo) {
if(iParticle.pt() < fRMSPtMin[iAlgo]) return;
// Change from SRR : Previously used fastjet::PseudoJet::user_index to decide the particle type.
Expand Down Expand Up @@ -126,8 +153,8 @@ void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {
for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) {
fMean[iAlgo] += fPups[i0];
if(fPups[i0] == 0) continue;
if(!fCharged[iAlgo] && fAdjust[iAlgo] && fPups[i0] > lMed) continue;
//if(fAdjust[iAlgo] && fPups[i0] > lMed) continue;
// if(!fCharged[iAlgo] && fAdjust[iAlgo] && fPups[i0] > lMed) continue;
if(fAdjust[iAlgo] && fPups[i0] > lMed) continue;
lNRMS++;
fRMS [iAlgo] += (fPups[i0]-lMed)*(fPups[i0]-lMed);
}
Expand All @@ -140,22 +167,32 @@ void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {
// some ways to do corrections to fRMS and fMedian
fRMS [iAlgo] *= fRMSScaleFactor[iAlgo];

fRMS[iAlgo] *= fRMSEtaSF;
fMedian[iAlgo] *= fMedEtaSF;

if(!fAdjust[iAlgo]) return;
//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(fPupsPV.size()+fNCount[iAlgo]);
if(lAdjust > 0) fMedian[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*fRMS[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]);
}
}

// fRMS_perEta[iAlgo] *= cur_RMSEtaSF;
// fMedian_perEta[iAlgo] *= cur_MedEtaSF;

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> const &iVals,double iChi2) const {

if(fAlgoId[0] == -1) return 1;
double lVal = 0.;
double lPVal = 1.;
Expand All @@ -170,15 +207,16 @@ double PuppiAlgo::compute(std::vector<double> const &iVals,double iChi2) const {
}
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];
if(fAlgoId[i0] == 0 && iVals[i0] == 0) pVal = cur_Med;
if(fAlgoId[i0] == 3 && iVals[i0] == 0) pVal = cur_Med;
if(fAlgoId[i0] == 5 && iVals[i0] == 0) pVal = cur_Med;
lVal += (pVal-cur_Med)*(fabs(pVal-cur_Med))/cur_RMS/cur_RMS;
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;

}
31 changes: 18 additions & 13 deletions CommonTools/PileupAlgos/src/PuppiContainer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,11 +47,11 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
// float nom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt)*(cosh(fRecoParticle.eta))*(cosh(fRecoParticle.eta))) + (fRecoParticle.pt)*sinh(fRecoParticle.eta);//hacked
// float denom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt));//hacked
// float rapidity = log(nom/denom);//hacked
if (edm::isFinite(fRecoParticle.rapidity)){
curPseudoJet.reset_PtYPhiM(fRecoParticle.pt,fRecoParticle.rapidity,fRecoParticle.phi,fRecoParticle.m);//hacked
} else {
curPseudoJet.reset_PtYPhiM(0, 99., 0, 0);//skipping may have been a better choice
}
if (edm::isFinite(fRecoParticle.rapidity)){
curPseudoJet.reset_PtYPhiM(fRecoParticle.pt,fRecoParticle.rapidity,fRecoParticle.phi,fRecoParticle.m);//hacked
} else {
curPseudoJet.reset_PtYPhiM(0, 99., 0, 0);//skipping may have been a better choice
}
//curPseudoJet.reset_PtYPhiM(fRecoParticle.pt,fRecoParticle.eta,fRecoParticle.phi,fRecoParticle.m);
int puppi_register = 0;
if(fRecoParticle.id == 0 or fRecoParticle.charge == 0) puppi_register = 0; // zero is neutral hadron
Expand Down Expand Up @@ -168,11 +168,16 @@ void PuppiContainer::getRawAlphas(int iOpt,std::vector<fastjet::PseudoJet> const
int PuppiContainer::getPuppiId( float iPt, float iEta) {
int lId = -1;
for(int i0 = 0; i0 < fNAlgos; i0++) {
if(std::abs(iEta) < fPuppiAlgo[i0].etaMin()) continue;
if(std::abs(iEta) > fPuppiAlgo[i0].etaMax()) continue;
if(iPt < fPuppiAlgo[i0].ptMin()) continue;
lId = i0;
break;
int nEtaBinsPerAlgo = fPuppiAlgo[i0].etaBins();
for (int i1 = 0; i1 < nEtaBinsPerAlgo; i1++){
if ( (std::abs(iEta) > fPuppiAlgo[i0].etaMin(i1)) && (std::abs(iEta) < fPuppiAlgo[i0].etaMax(i1)) ){
fPuppiAlgo[i0].fixAlgoEtaBin( i1 );
if(iPt > fPuppiAlgo[i0].ptMin()){
lId = i0;
break;
}
}
}
}
//if(lId == -1) std::cerr << "Error : Full fiducial range is not defined " << std::endl;
return lId;
Expand Down Expand Up @@ -246,16 +251,16 @@ std::vector<double> const & PuppiContainer::puppiWeights() {
//std::cout << "fRecoParticles[i0].pt = " << fRecoParticles[i0].pt << ", fRecoParticles[i0].charge = " << fRecoParticles[i0].charge << ", fRecoParticles[i0].id = " << fRecoParticles[i0].id << ", weight = " << pWeight << std::endl;

fWeights .push_back(pWeight);
fAlphaMed.push_back(fPuppiAlgo[pPupId].median(0));
fAlphaRMS.push_back(fPuppiAlgo[pPupId].rms(0));
fAlphaMed.push_back(fPuppiAlgo[pPupId].median());
fAlphaRMS.push_back(fPuppiAlgo[pPupId].rms());
//Now get rid of the thrown out weights for the particle collection

// leave these lines in, in case want to move eventually to having no 1-to-1 correspondence between puppi and pf cands
// if( std::abs(pWeight) < std::numeric_limits<double>::denorm_min() ) continue; // this line seems not to work like it's supposed to...
// if(std::abs(pWeight) <= 0. ) continue;

//Produce
PseudoJet curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e());
PseudoJet curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e() );
curjet.set_user_index(i0);
fPupParticles.push_back(curjet);
}
Expand Down
2 changes: 2 additions & 0 deletions RecoJets/Configuration/python/RecoJets_EventContent_cff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@
'keep *_fixedGridRhoFastjetAll_*_*',
'keep *_fixedGridRhoFastjetAllTmp_*_*',
'keep *_fixedGridRhoFastjetAllCalo_*_*',
'keep *_fixedGridRhoFastjetCentral_*_*',
'keep *_fixedGridRhoFastjetCentralCalo_*_*',
'keep *_fixedGridRhoFastjetCentralChargedPileUp_*_*',
'keep *_fixedGridRhoFastjetCentralNeutral_*_*',
Expand Down Expand Up @@ -117,6 +118,7 @@
'keep *_fixedGridRhoAll_*_*',
'keep *_fixedGridRhoFastjetAll_*_*',
'keep *_fixedGridRhoFastjetAllTmp_*_*',
'keep *_fixedGridRhoFastjetCentral_*_*',
'keep *_fixedGridRhoFastjetAllCalo_*_*',
'keep *_fixedGridRhoFastjetCentralCalo_*_*',
'keep *_fixedGridRhoFastjetCentralChargedPileUp_*_*',
Expand Down
Loading

0 comments on commit 3f23db9

Please sign in to comment.