Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

update PUPPI based on DP note #10597

Merged
merged 7 commits into from Aug 20, 2015
Merged

update PUPPI based on DP note #10597

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4466938
update PUPPI based on DP note
nhanvtran Aug 5, 2015
7aa2807
reindenting
nhanvtran Aug 17, 2015
53a9d9f
a few line-by-liners and two mis-initializations
nhanvtran Aug 18, 2015
bd7e8f8
speed up PATPackedCandidateProducer.cc
nhanvtran Aug 19, 2015
ca7126c
more speeding up
nhanvtran Aug 19, 2015
35e96ca
one more small optimization
nhanvtran Aug 19, 2015
06dcff6
diagnostics turned off to speed it up
nhanvtran Aug 19, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 3 additions & 0 deletions CommonTools/PileupAlgos/interface/PuppiContainer.h
View file
Open in desktop
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 Down Expand Up @@ -68,6 +70,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
Open in desktop
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
Open in desktop
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
Open in desktop
@@ -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
Open in desktop
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
Open in desktop
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