Skip to content

Commit

Permalink
Merge CMSSW_7_5_X into CMSSW_7_6_X.
Browse files Browse the repository at this point in the history
  • Loading branch information
@cmsbuild
cmsbuild committed Aug 20, 2015
2 parents 3249e5c + 2c91b74 commit 9d4aa93
Show file tree
Hide file tree
Showing 24 changed files with 554 additions and 103 deletions.
4 changes: 4 additions & 0 deletions CommonTools/PileupAlgos/interface/PuppiContainer.h
View file
Expand Up @@ -37,6 +37,8 @@ class PuppiContainer{
PuppiContainer(const edm::ParameterSet &iConfig);
~PuppiContainer();
void initialize(const std::vector<RecoObj> &iRecoObjects);
void setNPV(int iNPV){ fNPV = iNPV; }

std::vector<fastjet::PseudoJet> const & pfParticles() const { return fPFParticles; }
std::vector<fastjet::PseudoJet> const & pvParticles() const { return fChargedPV; }
std::vector<double> const & puppiWeights();
Expand All @@ -57,6 +59,7 @@ class PuppiContainer{
int getPuppiId ( float iPt, float iEta);
double var_within_R (int iId, const std::vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, double R);

bool fPuppiDiagnostics;
std::vector<RecoObj> fRecoParticles;
std::vector<fastjet::PseudoJet> fPFParticles;
std::vector<fastjet::PseudoJet> fChargedPV;
Expand All @@ -68,6 +71,7 @@ class PuppiContainer{
std::vector<double> fAlphaRMS;

bool fApplyCHS;
bool fInvert;
bool fUseExp;
double fNeutralMinPt;
double fNeutralSlope;
Expand Down
134 changes: 104 additions & 30 deletions CommonTools/PileupAlgos/plugins/PuppiProducer.cc
View file
Expand Up @@ -30,8 +30,13 @@
// ------------------------------------------------------------------------------------------
PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) {
fPuppiDiagnostics = iConfig.getParameter<bool>("puppiDiagnostics");
fPuppiForLeptons = iConfig.getParameter<bool>("puppiForLeptons");
fUseDZ = iConfig.getParameter<bool>("UseDeltaZCut");
fDZCut = iConfig.getParameter<double>("DeltaZCut");
fUseExistingWeights = iConfig.getParameter<bool>("useExistingWeights");
fUseWeightsNoLep = iConfig.getParameter<bool>("useWeightsNoLep");
fVtxNdofCut = iConfig.getParameter<int>("vtxNdofCut");
fVtxZCut = iConfig.getParameter<double>("vtxZCut");
fPuppiContainer = std::unique_ptr<PuppiContainer> ( new PuppiContainer(iConfig) );

tokenPFCandidates_
Expand Down Expand Up @@ -70,9 +75,17 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
iEvent.getByToken(tokenVertices_,hVertexProduct);
const reco::VertexCollection *pvCol = hVertexProduct.product();

int npv = 0;
const reco::VertexCollection::const_iterator vtxEnd = pvCol->end();
for (reco::VertexCollection::const_iterator vtxIter = pvCol->begin(); vtxEnd != vtxIter; ++vtxIter) {
if (!vtxIter->isFake() && vtxIter->ndof()>=fVtxNdofCut && fabs(vtxIter->z())<=fVtxZCut)
npv++;
}

//Fill the reco objects
fRecoObjCollection.clear();
for(CandidateView::const_iterator itPF = pfCol->begin(); itPF!=pfCol->end(); itPF++) {
// std::cout << "itPF->pdgId() = " << itPF->pdgId() << std::endl;
RecoObj pReco;
pReco.pt = itPF->pt();
pReco.eta = itPF->eta();
Expand All @@ -87,7 +100,10 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
bool lFirst = true;
const pat::PackedCandidate *lPack = dynamic_cast<const pat::PackedCandidate*>(&(*itPF));
if(lPack == 0 ) {

const reco::PFCandidate *pPF = dynamic_cast<const reco::PFCandidate*>(&(*itPF));
double curdz = 9999;
int closestVtxForUnassociateds = -9999;
for(reco::VertexCollection::const_iterator iV = pvCol->begin(); iV!=pvCol->end(); ++iV) {
if(lFirst) {
if ( pPF->trackRef().isNonnull() ) pDZ = pPF->trackRef() ->dz(iV->position());
Expand All @@ -100,43 +116,99 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
if(iV->trackWeight(pPF->trackRef())>0) {
closestVtx = &(*iV);
break;
}
}
// in case it's unassocciated, keep more info
double tmpdz = 99999;
if ( pPF->trackRef().isNonnull() ) tmpdz = pPF->trackRef() ->dz(iV->position());
else if ( pPF->gsfTrackRef().isNonnull() ) tmpdz = pPF->gsfTrackRef()->dz(iV->position());
if (fabs(tmpdz) < curdz){
curdz = fabs(tmpdz);
closestVtxForUnassociateds = pVtxId;
}
pVtxId++;

}
} else if(lPack->vertexRef().isNonnull() ) {
pDZ = lPack->dz();
pD0 = lPack->dxy();
closestVtx = &(*(lPack->vertexRef()));
pVtxId = (lPack->fromPV() != (pat::PackedCandidate::PVUsedInFit));
if( (lPack->fromPV() == pat::PackedCandidate::PVLoose) || (lPack->fromPV() == pat::PackedCandidate::PVTight) ){
pVtxId = 0;
}
}
int tmpFromPV = 0;
// mocking the miniAOD definitions
if (closestVtx != 0 && fabs(pReco.charge) > 0 && pVtxId > 0) tmpFromPV = 0;
if (closestVtx != 0 && fabs(pReco.charge) > 0 && pVtxId == 0) tmpFromPV = 3;
if (closestVtx == 0 && fabs(pReco.charge) > 0 && closestVtxForUnassociateds == 0) tmpFromPV = 2;
if (closestVtx == 0 && fabs(pReco.charge) > 0 && closestVtxForUnassociateds != 0) tmpFromPV = 1;
pReco.dZ = pDZ;
pReco.d0 = pD0;

if(closestVtx == 0) pReco.vtxId = -1;
if(closestVtx != 0) pReco.vtxId = pVtxId;
//if(closestVtx != 0) pReco.vtxChi2 = closestVtx->trackWeight(itPF->trackRef());
//Set the id for Puppi Algo: 0 is neutral pfCandidate, id = 1 for particles coming from PV and id = 2 for charged particles from non-leading vertex
pReco.id = 0;

if(closestVtx != 0 && pVtxId == 0 && fabs(pReco.charge) > 0) pReco.id = 1;
if(closestVtx != 0 && pVtxId > 0 && fabs(pReco.charge) > 0) pReco.id = 2;
//Add a dZ cut if wanted (this helps)
if(fUseDZ && pDZ > -9999 && closestVtx == 0 && (fabs(pDZ) < fDZCut) && fabs(pReco.charge) > 0) pReco.id = 1;
if(fUseDZ && pDZ > -9999 && closestVtx == 0 && (fabs(pDZ) > fDZCut) && fabs(pReco.charge) > 0) pReco.id = 2;

//std::cout << "pVtxId = " << pVtxId << ", and charge = " << itPF->charge() << ", and closestVtx = " << closestVtx << ", and id = " << pReco.id << std::endl;
//std::cout << "charge = " << itPF->charge() << ", pDZ = " << pDZ << ", pVtxId = " << pVtxId << ", closestVtx = " << closestVtx << ", fromPV() = " << lPack->fromPV() << ", pReco.id = " << pReco.id << std::endl;
pReco.id = 0;
if (fabs(pReco.charge) == 0){ pReco.id = 0; }
else{
if (tmpFromPV == 0){ pReco.id = 2; } // 0 is associated to PU vertex
if (tmpFromPV == 3){ pReco.id = 1; }
if (tmpFromPV == 1 || tmpFromPV == 2){
pReco.id = 0;
if (!fPuppiForLeptons && fUseDZ && (fabs(pDZ) < fDZCut)) pReco.id = 1;
if (!fPuppiForLeptons && fUseDZ && (fabs(pDZ) > fDZCut)) pReco.id = 2;
if (fPuppiForLeptons && tmpFromPV == 1) pReco.id = 2;
if (fPuppiForLeptons && tmpFromPV == 2) pReco.id = 1;
}
}
}
else if(lPack->vertexRef().isNonnull() ) {
pDZ = lPack->dz();
pD0 = lPack->dxy();
closestVtx = &(*(lPack->vertexRef()));
pReco.dZ = pDZ;
pReco.d0 = pD0;

pReco.id = 0;
if (fabs(pReco.charge) == 0){ pReco.id = 0; }
if (fabs(pReco.charge) > 0){
if (lPack->fromPV() == 0){ pReco.id = 2; } // 0 is associated to PU vertex
if (lPack->fromPV() == (pat::PackedCandidate::PVUsedInFit)){ pReco.id = 1; }
if (lPack->fromPV() == (pat::PackedCandidate::PVTight) || lPack->fromPV() == (pat::PackedCandidate::PVLoose)){
pReco.id = 0;
if (!fPuppiForLeptons && fUseDZ && (fabs(pDZ) < fDZCut)) pReco.id = 1;
if (!fPuppiForLeptons && fUseDZ && (fabs(pDZ) > fDZCut)) pReco.id = 2;
if (fPuppiForLeptons && lPack->fromPV() == (pat::PackedCandidate::PVLoose)) pReco.id = 2;
if (fPuppiForLeptons && lPack->fromPV() == (pat::PackedCandidate::PVTight)) pReco.id = 1;
}
}
}

fRecoObjCollection.push_back(pReco);
fRecoObjCollection.push_back(pReco);

}

fPuppiContainer->initialize(fRecoObjCollection);
fPuppiContainer->setNPV( npv );

std::vector<double> lWeights;
std::vector<fastjet::PseudoJet> lCandidates;
if (!fUseExistingWeights){
//Compute the weights and get the particles
lWeights = fPuppiContainer->puppiWeights();
lCandidates = fPuppiContainer->puppiParticles();
}
else{
//Use the existing weights
int lPackCtr = 0;
for(CandidateView::const_iterator itPF = pfCol->begin(); itPF!=pfCol->end(); itPF++) {
const pat::PackedCandidate *lPack = dynamic_cast<const pat::PackedCandidate*>(&(*itPF));
float curpupweight = -1.;
if(lPack == 0 ) {
// throw error
throw edm::Exception(edm::errors::LogicError,"PuppiProducer: cannot get weights since inputs are not packedPFCandidates");
}
else{
// if (fUseWeightsNoLep){ curpupweight = itPF->puppiWeightNoLep(); }
// else{ curpupweight = itPF->puppiWeight(); }
curpupweight = lPack->puppiWeight();
}
lWeights.push_back(curpupweight);
fastjet::PseudoJet curjet( curpupweight*lPack->px(), curpupweight*lPack->py(), curpupweight*lPack->pz(), curpupweight*lPack->energy());
curjet.set_user_index(lPackCtr);
lCandidates.push_back(curjet);
lPackCtr++;
}
}

//Compute the weights
const std::vector<double> lWeights = fPuppiContainer->puppiWeights();
//Fill it into the event
std::auto_ptr<edm::ValueMap<float> > lPupOut(new edm::ValueMap<float>());
edm::ValueMap<float>::Filler lPupFiller(*lPupOut);
Expand All @@ -152,7 +224,6 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
// Since the size of the ValueMap must be equal to the input collection, we need
// to search the "puppi" particles to find a match for each input. If none is found,
// the input is set to have a four-vector of 0,0,0,0
const std::vector<fastjet::PseudoJet> lCandidates = fPuppiContainer->puppiParticles();
fPuppiCandidates.reset( new PFOutputCollection );
std::auto_ptr<edm::ValueMap<LorentzVector> > p4PupOut(new edm::ValueMap<LorentzVector>());
LorentzVectorCollection puppiP4s;
Expand All @@ -174,11 +245,14 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
auto puppiMatched = find_if( lCandidates.begin(), lCandidates.end(), [&val]( fastjet::PseudoJet const & i ){ return i.user_index() == val; } );
if ( puppiMatched != lCandidates.end() ) {
pVec.SetPxPyPzE(puppiMatched->px(),puppiMatched->py(),puppiMatched->pz(),puppiMatched->E());
// fPuppiCandidates->push_back(pCand);
} else {
pVec.SetPxPyPzE( 0, 0, 0, 0);
}
pCand.setP4(pVec);
puppiP4s.push_back( pVec );

pCand.setSourceCandidatePtr( i0->sourceCandidatePtr(0) );
fPuppiCandidates->push_back(pCand);
}

Expand All @@ -203,7 +277,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {


//////////////////////////////////////////////
if (fPuppiDiagnostics){
if (fPuppiDiagnostics && !fUseExistingWeights){

// all the different alphas per particle
// THE alpha per particle
Expand Down
5 changes: 5 additions & 0 deletions CommonTools/PileupAlgos/plugins/PuppiProducer.h
View file
Expand Up @@ -47,8 +47,13 @@ class PuppiProducer : public edm::stream::EDProducer<> {
std::string fPFName;
std::string fPVName;
bool fPuppiDiagnostics;
bool fPuppiForLeptons;
bool fUseDZ;
float fDZCut;
bool fUseExistingWeights;
bool fUseWeightsNoLep;
int fVtxNdofCut;
double fVtxZCut;
std::unique_ptr<PuppiContainer> fPuppiContainer;
std::vector<RecoObj> fRecoObjCollection;
std::auto_ptr< PFOutputCollection > fPuppiCandidates;
Expand Down
7 changes: 7 additions & 0 deletions CommonTools/PileupAlgos/python/PUPuppi_cff.py
View file
@@ -0,0 +1,7 @@
import FWCore.ParameterSet.Config as cms

from CommonTools.PileupAlgos.Puppi_cff import *

pupuppi = puppi.clone()
pupuppi.invertPuppi = True

45 changes: 20 additions & 25 deletions CommonTools/PileupAlgos/python/Puppi_cff.py
View file
Expand Up @@ -6,7 +6,7 @@
useCharged = cms.bool(True),
applyLowPUCorr = cms.bool(True),
combOpt = cms.int32(0),
cone = cms.double(0.3),
cone = cms.double(0.4),
rmsPtMin = cms.double(0.1),
rmsScaleFactor = cms.double(1.0)
)
Expand All @@ -18,69 +18,64 @@
useCharged = cms.bool(False),
applyLowPUCorr = cms.bool(True),
combOpt = cms.int32(0),
cone = cms.double(0.3),
rmsPtMin = cms.double(0.1),
cone = cms.double(0.4),
rmsPtMin = cms.double(0.5),
rmsScaleFactor = cms.double(1.0)
)
)

puppi = cms.EDProducer("PuppiProducer",#cms.PSet(#"PuppiProducer",
puppiDiagnostics = cms.bool(False),
UseDeltaZCut = cms.bool(False),
puppiForLeptons = cms.bool(False),
UseDeltaZCut = cms.bool(True),
DeltaZCut = cms.double(0.3),
candName = cms.InputTag('particleFlow'),
vertexName = cms.InputTag('offlinePrimaryVertices'),
#candName = cms.string('packedPFCandidates'),
#vertexName = cms.string('offlineSlimmedPrimaryVertices'),
applyCHS = cms.bool (True),
invertPuppi = cms.bool (False),
useExp = cms.bool (False),
MinPuppiWeight = cms.double(0.01),
useExistingWeights = cms.bool(False),
useWeightsNoLep = cms.bool(False),
vtxNdofCut = cms.int32(4),
vtxZCut = cms.double(24),
algos = cms.VPSet(
cms.PSet(
etaMin = cms.double(0.),
etaMax = cms.double( 2.0),
etaMax = cms.double(2.5),
ptMin = cms.double(0.),
MinNeutralPt = cms.double(0.1),
MinNeutralPtSlope = cms.double(0.005),
MinNeutralPtSlope = cms.double(0.015),
RMSEtaSF = cms.double(1.0),
MedEtaSF = cms.double(1.0),
EtaMaxExtrap = cms.double(2.0),
puppiAlgos = puppiCentral
),
cms.PSet(
etaMin = cms.double(2.0),
etaMax = cms.double( 2.5),
ptMin = cms.double(0.),
MinNeutralPt = cms.double(0.1),
MinNeutralPtSlope = cms.double(0.005),
RMSEtaSF = cms.double(0.95),
MedEtaSF = cms.double(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.0),
MinNeutralPtSlope = cms.double(0.005),
MinNeutralPt = cms.double(1.7),
MinNeutralPtSlope = cms.double(0.07),
# RMSEtaSF = cms.double(1.545),
# MedEtaSF = cms.double(0.845),
RMSEtaSF = cms.double(1.484),
MedEtaSF = cms.double(0.845),
RMSEtaSF = cms.double(1.30),
MedEtaSF = cms.double(1.05),
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(1.0),
MinNeutralPtSlope = cms.double(0.005),
MinNeutralPt = cms.double(2.0),
MinNeutralPtSlope = cms.double(0.07),
# RMSEtaSF = cms.double(1.18),
# MedEtaSF = cms.double(0.4397),
RMSEtaSF = cms.double(1.25),
MedEtaSF = cms.double(0.4397),
RMSEtaSF = cms.double(1.10),
MedEtaSF = cms.double(0.90),
EtaMaxExtrap = cms.double(2.0),
puppiAlgos = puppiForward
)
Expand Down
4 changes: 2 additions & 2 deletions CommonTools/PileupAlgos/src/PuppiAlgo.cc
View file
Expand Up @@ -94,6 +94,7 @@ void PuppiAlgo::add(const fastjet::PseudoJet &iParticle,const double &iVal,const
}
// 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);

}

////////////////////////////////////////////////////////////////////////////////
Expand All @@ -109,7 +110,7 @@ void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {
int lNBefore = 0;
for(unsigned int i0 = 0; i0 < iAlgo; i0++) lNBefore += fNCount[i0];
std::sort(fPups.begin()+lNBefore,fPups.begin()+lNBefore+fNCount[iAlgo]);

// in case you have alphas == 0
int lNum0 = 0;
for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) {
Expand Down Expand Up @@ -141,7 +142,6 @@ void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {

fRMS[iAlgo] *= fRMSEtaSF;
fMedian[iAlgo] *= fMedEtaSF;
//if(!fCharged[iAlgo]) std::cout << " Process : " << iAlgo << " Median : " << fMedian[iAlgo] << " +/- " << fRMS[iAlgo] << " -- Begin : " << lNBefore << " -- Total : " << fNCount[iAlgo] << " -- 50% " << lNHalfway << " Fraction less than @ Median : " << std::endl;

if(!fAdjust[iAlgo]) return;
//Adjust the p-value to correspond to the median
Expand Down

0 comments on commit 9d4aa93

Please sign in to comment.