Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions JetAnalyzers/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
<use name="DataFormats/HepMCCandidate"/>
<use name="DataFormats/ParticleFlowCandidate"/>
<use name="DataFormats/Common"/>
<use name="DataFormats/PatCandidates"/>
<use name="SimDataFormats/GeneratorProducts"/>
<use name="SimDataFormats/JetMatching"/>
<use name="PhysicsTools/Utilities"/>
Expand Down
1 change: 1 addition & 0 deletions JetAnalyzers/bin/jet_correction_analyzer_x.cc
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
#include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
#if __has_include("xrootd/XrdCl/XrdClFileSystem.hh")
#include "xrootd/XProtocol/XProtocol.hh"
#include "xrootd/XrdCl/XrdClFileSystem.hh"
#define has_xrdcl 1
#else
Expand Down
4 changes: 2 additions & 2 deletions JetAnalyzers/interface/JetResponseAnalyzer.hh
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/PatCandidates/interface/PackedGenParticle.h"

#include "JetMETCorrections/Objects/interface/JetCorrector.h"

Expand Down Expand Up @@ -101,10 +102,9 @@ private:
edm::EDGetTokenT<reco::VertexCollection> srcVtx_;
edm::EDGetTokenT<GenEventInfoProduct> srcGenInfo_;
edm::EDGetTokenT<vector<PileupSummaryInfo> > srcPileupInfo_;
//edm::EDGetTokenT<vector<reco::PFCandidate> > srcPFCandidates_;
edm::EDGetTokenT<PFCandidateView> srcPFCandidates_;
edm::EDGetTokenT<std::vector<edm::FwdPtr<reco::PFCandidate> > > srcPFCandidatesAsFwdPtr_;
edm::EDGetTokenT<vector<reco::GenParticle> > srcGenParticles_;
edm::EDGetTokenT<vector<pat::PackedGenParticle> > srcGenParticles_;

std::string jecLabel_;

Expand Down
2 changes: 1 addition & 1 deletion JetAnalyzers/python/Defaults_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
JetResponseParameters = cms.PSet(
# record flavor information, consider both RefPt and JetPt
doComposition = cms.bool(True),
doFlavor = cms.bool(True),
doFlavor = cms.bool(False),
doRefPt = cms.bool(True),
doJetPt = cms.bool(True),
saveCandidates = cms.bool(False),
Expand Down
2 changes: 1 addition & 1 deletion JetAnalyzers/python/JetReconstruction_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@
#!
#! PF JETS CHS
#!
ak1PFCHSJets = ak1PFJets.clone( src = 'pfNoPileUpJME' )
ak1PFCHSJets = ak1PFJets.clone( src = 'pfCHS' )
ak2PFCHSJets = ak1PFCHSJets.clone( rParam=0.2 )
ak3PFCHSJets = ak1PFCHSJets.clone( rParam=0.3 )
ak4PFCHSJets = ak1PFCHSJets.clone( rParam=0.4 )
Expand Down
45 changes: 25 additions & 20 deletions JetAnalyzers/python/addAlgorithm.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
################################################################################

partons = cms.EDProducer('PartonSelector',
src = cms.InputTag('genParticles'),
src = cms.InputTag('packedGenParticles'),
withLeptons = cms.bool(False),
skipFirstN = cms.uint32(0)
)
Expand All @@ -20,6 +20,9 @@
from JetMETAnalysis.JetAnalyzers.JetCorrection_cff import *
from RecoTauTag.TauTagTools.tauDecayModes_cfi import *
from CommonTools.PileupAlgos.Puppi_cff import *
from JetMETAnalysis.JetAnalyzers.customizePuppiTune_cff import *

genParticlesForJetsNoNu.src = cms.InputTag("packedGenParticles")

stdClusteringAlgorithms = ['ak'] #Options: {ak,kt}
stdJetTypes = ['calo','pf','pfchs','puppi'] #Options: {'calo','pf','pfchs','puppi'}
Expand Down Expand Up @@ -382,7 +385,12 @@ def addAlgorithm(process, alg_size_type_corr, Defaults, reco, doProducer):
sequence = cms.Sequence(recJets * sequence)
if type == 'PUPPI':
process.load('CommonTools.PileupAlgos.Puppi_cff')
process.load('JetMETAnalysis.JetAnalyzers.customizePuppiTune_cff')
#puppi.candName = cms.InputTag("particleFlow")
#80x change
puppiCentral[0].applyLowPUCorr = cms.bool(False)
puppiForward[0].applyLowPUCorr = cms.bool(False)
puppi.vertexName = "offlineSlimmedPrimaryVertices"
sequence = cms.Sequence(puppi * sequence)
if type == 'Track':
process.load('JetMETAnalysis.JetAnalyzers.TrackJetReconstruction_cff')
Expand Down Expand Up @@ -505,10 +513,10 @@ def addAlgorithm(process, alg_size_type_corr, Defaults, reco, doProducer):
srcRhos = cms.InputTag(''),
srcRho = cms.InputTag(''),
srcRhoHLT = cms.InputTag(''),
srcVtx = cms.InputTag('offlinePrimaryVertices'),
srcVtx = cms.InputTag('offlineSlimmedPrimaryVertices'),
srcJetToUncorJetMap = cms.InputTag(jetToUncorJet.label(), 'rec2gen'),
srcPFCandidates = cms.InputTag(''),
srcGenParticles = cms.InputTag('genParticles')
srcGenParticles = cms.InputTag('packedGenParticles')
)
if doProducer:
jraAnalyzer = 'JetResponseAnalyzerProducer'
Expand All @@ -519,10 +527,10 @@ def addAlgorithm(process, alg_size_type_corr, Defaults, reco, doProducer):
jecLabel = cms.string(''),
srcRho = cms.InputTag(''),
srcRhoHLT = cms.InputTag(''),
srcVtx = cms.InputTag('offlinePrimaryVertices'),
srcVtx = cms.InputTag('offlineSlimmedPrimaryVertices'),
srcJetToUncorJetMap = cms.InputTag(jetToUncorJet.label(), 'rec2gen'),
srcPFCandidates = cms.InputTag(''),
srcGenParticles = cms.InputTag('genParticles')
srcGenParticles = cms.InputTag('packedGenParticles')
)

if type == 'CaloHLT':
Expand All @@ -531,35 +539,38 @@ def addAlgorithm(process, alg_size_type_corr, Defaults, reco, doProducer):
elif type == 'Calo':
jra.srcRho = cms.InputTag("fixedGridRhoFastjetAllCalo")
elif type == 'PFchs':
process.kt6PFchsJetsRhos = kt6PFJets.clone(src = 'pfNoPileUpJME',
process.pfCHS = cms.EDFilter("CandPtrSelector", src = cms.InputTag("packedPFCandidates"), cut = cms.string("fromPV"))
process.kt6PFchsJetsRhos = kt6PFJets.clone(src = 'pfCHS',
doFastJetNonUniform = cms.bool(True),
puCenters = cms.vdouble(-5,-4,-3,-2,-1,0,1,2,3,4,5),
puWidth = cms.double(.8),
nExclude = cms.uint32(2))
sequence = cms.Sequence(process.kt6PFchsJetsRhos * sequence)
sequence = cms.Sequence(process.pfCHS * process.kt6PFchsJetsRhos * sequence)
jra.srcRhos = cms.InputTag("kt6PFchsJetsRhos", "rhos")
jra.srcRho = cms.InputTag("fixedGridRhoFastjetAll")
jra.srcPFCandidates = cms.InputTag('pfNoPileUpJME')
jra.srcPFCandidates = cms.InputTag('pfCHS')
elif type == 'PFHLT':
jra.srcRho = ak4PFL1Fastjet.srcRho #added 02/15/2012
jra.srcRho = ak4PFL1Fastjet.srcRho
jra.srcRhoHLT = ak5PFHLTL1Fastjet.srcRho
elif type == 'PFchsHLT':
jra.srcRho = ak4PFchsL1Fastjet.srcRho #added 02/15/2012
jra.srcRho = ak4PFchsL1Fastjet.srcRho
jra.srcRhoHLT = ak5PFchsHLTL1Fastjet.srcRho
elif type == 'PF':
process.particleFlow = cms.EDFilter("CandPtrSelector", src = cms.InputTag("packedPFCandidates"), cut = cms.string(""))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are you recreating the particleFlow collection without any cuts? Why not just use the packedPFCandidates directly as inputs on line 583? Same comment for lines 585-589.

process.kt6PFJetsRhos = kt6PFJets.clone(doFastJetNonUniform = cms.bool(True),
puCenters = cms.vdouble(-5,-4,-3,-2,-1,0,1,2,3,4,5),
puWidth = cms.double(.8),
nExclude = cms.uint32(2))
sequence = cms.Sequence(process.kt6PFJetsRhos * sequence)
sequence = cms.Sequence(process.particleFlow * process.kt6PFJetsRhos * sequence)
jra.srcRhos = cms.InputTag("kt6PFJetsRhos", "rhos")
jra.srcRho = cms.InputTag("fixedGridRhoFastjetAll")
jra.srcPFCandidates = cms.InputTag('particleFlow')
jra.srcPFCandidates = cms.InputTag('packedPFCandidates')
elif type == 'PUPPI':
process.particleFlow = cms.EDFilter("CandPtrSelector", src = cms.InputTag("packedPFCandidates"), cut = cms.string(""))
process.kt6PFJetsRhos = kt6PFJets.clone(doFastJetNonUniform = cms.bool(True),
puCenters = cms.vdouble(-5,-4,-3,-2,-1,0,1,2,3,4,5),
puWidth = cms.double(.8), nExclude = cms.uint32(2))
sequence = cms.Sequence(process.kt6PFJetsRhos * sequence)
sequence = cms.Sequence(process.particleFlow * process.kt6PFJetsRhos * sequence)
jra.srcRhos = cms.InputTag("kt6PFJetsRhos", "rhos")
jra.srcRho = cms.InputTag("fixedGridRhoFastjetAll")
jra.srcPFCandidates = cms.InputTag('puppi')
Expand All @@ -574,14 +585,8 @@ def addAlgorithm(process, alg_size_type_corr, Defaults, reco, doProducer):

setattr(process,alg_size_type_corr,jra)
sequence = cms.Sequence(sequence * jra)

## add chs to path is needed
if type == 'PFchs':
sequence = cms.Sequence(process.pfNoPileUpJMESequence * sequence)

## create the path and put in the sequence
sequence = cms.Sequence(sequence)
setattr(process, alg_size_type_corr + 'Sequence', sequence)
path = cms.Path( sequence )
setattr(process, alg_size_type_corr + 'Path', path)
setattr(process, alg_size_type_corr + 'Path', path)
print alg_size_type_corr
74 changes: 58 additions & 16 deletions JetAnalyzers/src/JetResponseAnalyzer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,11 @@ JetResponseAnalyzer::JetResponseAnalyzer(const edm::ParameterSet& iConfig)
, srcRhoHLT_ (consumes<double>(iConfig.getParameter<edm::InputTag> ("srcRhoHLT")))
, srcVtx_ (consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag> ("srcVtx")))
, srcGenInfo_ (consumes<GenEventInfoProduct>(edm::InputTag("generator")) )
, srcPileupInfo_ (consumes<vector<PileupSummaryInfo> >(edm::InputTag("addPileupInfo")) )
, srcPileupInfo_ (consumes<vector<PileupSummaryInfo> >(edm::InputTag("slimmedAddPileupInfo")) )
//, srcPFCandidates_ (consumes<vector<reco::PFCandidate> >(iConfig.getParameter<edm::InputTag>("srcPFCandidates")))
, srcPFCandidates_ (consumes<PFCandidateView>(iConfig.getParameter<edm::InputTag>("srcPFCandidates")))
, srcPFCandidatesAsFwdPtr_(consumes<std::vector<edm::FwdPtr<reco::PFCandidate> > >(iConfig.getParameter<edm::InputTag>("srcPFCandidates")))
, srcGenParticles_ (consumes<vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("srcGenParticles")))
, srcGenParticles_ (consumes<vector<pat::PackedGenParticle> >(iConfig.getParameter<edm::InputTag>("srcGenParticles")))
, jecLabel_ (iConfig.getParameter<std::string> ("jecLabel"))
, doComposition_ (iConfig.getParameter<bool> ("doComposition"))
, doFlavor_ (iConfig.getParameter<bool> ("doFlavor"))
Expand Down Expand Up @@ -115,6 +115,46 @@ void JetResponseAnalyzer::beginJob()
}


void getMult( vector<reco::CandidatePtr> const & particles, int* nMult, int* chMult ) {

vector<reco::CandidatePtr>::const_iterator itParticle;
for (itParticle=particles.begin();itParticle!=particles.end();++itParticle){
const reco::Candidate* pfCand = itParticle->get();

switch (std::abs(pfCand->pdgId())) {

case 211: //PFCandidate::h: // charged hadron
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't use unnamed constants. Instead, you should define enums as in

enum Flavor{all=0, h, e, mu, gamma, h0, h_HF, egamma_HF, chs, numFlavors, X}; //X=undefined
and use those as your cases.

You also only really have two cases:

  1. You add one to the nMult variable
  2. You add one to the chMult variable
    You should combine these cases using and OR of the conditions inside if, else if, else statements. This will reduce the number of comparisons. You can use the std::any_of function to reduce the verbosity of the conditional (https://www.cplusplus.com/reference/algorithm/any_of/).

(*chMult)++;
break;

case 130: //PFCandidate::h0 : // neutral hadron
(*nMult)++;
break;

case 22: //PFCandidate::gamma: // photon
(*nMult)++;
break;

case 11: // PFCandidate::e: // electron
(*chMult)++;
break;

case 13: //PFCandidate::mu: // muon
(*chMult)++;
break;

case 1: // PFCandidate::h_HF : // hadron in HF
(*nMult)++;
break;

case 2: //PFCandidate::egamma_HF : // electromagnetic in HF
(*nMult)++;
break;
}
}
}

//_______________
//______________________________________________________________________________
void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
const edm::EventSetup& iSetup)
Expand All @@ -134,7 +174,7 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
edm::Handle<reco::VertexCollection> vtx;
edm::Handle<PFCandidateView> pfCandidates;
edm::Handle<std::vector<edm::FwdPtr<reco::PFCandidate> > > pfCandidatesAsFwdPtr;
edm::Handle<vector<reco::GenParticle> > genParticles;
edm::Handle<vector<pat::PackedGenParticle> > genParticles;

// Jet CORRECTOR
jetCorrector_ = (jecLabel_.empty()) ? 0 : JetCorrector::getJetCorrector(jecLabel_,iSetup);
Expand Down Expand Up @@ -189,8 +229,8 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
JRAEvt_->refpvz = -1000.0;
iEvent.getByToken(srcGenParticles_, genParticles);
for (size_t i = 0; i < genParticles->size(); ++i) {
const reco::GenParticle & genIt = (*genParticles)[i];
if ( genIt.isHardProcess() ) {
const pat::PackedGenParticle & genIt = (*genParticles)[i];
if ( genIt.fromHardProcessFinalState() ) {
JRAEvt_->refpvz = genIt.vz();
break;
}
Expand Down Expand Up @@ -264,11 +304,10 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
else JRAEvt_->refdrjt->pop_back();
continue;
}

JRAEvt_->refpdgid->push_back(0);
JRAEvt_->refpdgid_algorithmicDef->push_back(0);
JRAEvt_->refpdgid_physicsDef->push_back(0);
if (getFlavorFromMap_) {
JRAEvt_->refpdgid_algorithmicDef->push_back(0);
JRAEvt_->refpdgid_physicsDef->push_back(0);
reco::JetMatchedPartonsCollection::const_iterator itPartonMatch;
itPartonMatch=refToPartonMap->begin();
for (;itPartonMatch!=refToPartonMap->end();++itPartonMatch) {
Expand Down Expand Up @@ -315,10 +354,6 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
}
}
}
else {
JRAEvt_->refpdgid_algorithmicDef->at(JRAEvt_->nref)=0;
JRAEvt_->refpdgid_physicsDef->at(JRAEvt_->nref)=0;
}
JRAEvt_->refpdgid->at(JRAEvt_->nref)=ref->pdgId();

// Beta/Beta Star Calculation
Expand Down Expand Up @@ -439,12 +474,21 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
JRAEvt_->jtmuf ->push_back(pfJetRef->muonEnergyFraction() *JRAEvt_->jtjec->at(JRAEvt_->nref));
JRAEvt_->jthfhf->push_back(pfJetRef->HFHadronEnergyFraction() *JRAEvt_->jtjec->at(JRAEvt_->nref));
JRAEvt_->jthfef->push_back(pfJetRef->HFEMEnergyFraction() *JRAEvt_->jtjec->at(JRAEvt_->nref));
int chMult=0, nMult=0;
getMult( ref.castTo<reco::GenJetRef>()->getJetConstituents(), &nMult, &chMult );
JRAEvt_->refnMult ->push_back( nMult );
JRAEvt_->refchMult->push_back( chMult );

//this method exists for pfjets (neutralMultiplicity()), but not for genjets
//original i thought since genjet didn't have it i should make this method
chMult=0; nMult=0;
getMult( jet.castTo<reco::PFJetRef>()->getJetConstituents(), &nMult, &chMult );
JRAEvt_->jtnMult ->push_back( nMult );
JRAEvt_->jtchMult->push_back( chMult );
}
}

JRAEvt_->nref++;
}

// PFCANDIDATE INFORMATION
//Dual handle idea from https://github.com/aperloff/cmssw/blob/CMSSW_7_6_X/RecoJets/JetProducers/plugins/VirtualJetProducer.cc
//Random-Cone algo from https://github.com/cihar29/OffsetAnalysis/blob/master/run_offset.py
Expand Down Expand Up @@ -483,8 +527,6 @@ void JetResponseAnalyzer::analyze(const edm::Event& iEvent,
}
}
}


tree_->Fill();

return;
Expand Down
8 changes: 8 additions & 0 deletions JetUtilities/interface/JRAEvent.h
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,8 @@ public :
vector<Int_t>* refpdgid_physicsDef;
vector<Float_t>* refe;
vector<Float_t>* refpt;
vector<Int_t>* refnMult;
vector<Int_t>* refchMult;
vector<Float_t>* refeta;
vector<Float_t>* refphi;
vector<Float_t>* refy;
Expand All @@ -88,6 +90,8 @@ public :
vector<Float_t>* refarea;
vector<Float_t>* jte;
vector<Float_t>* jtpt;
vector<Int_t>* jtnMult;
vector<Int_t>* jtchMult;
vector<Float_t>* jteta;
vector<Float_t>* jtphi;
vector<Float_t>* jty;
Expand Down Expand Up @@ -140,6 +144,8 @@ public :
TBranch *b_refpdgid_physicsDef; //!
TBranch *b_refe; //!
TBranch *b_refpt; //!
TBranch *b_refnMult; //!
TBranch *b_refchMult; //!
TBranch *b_refeta; //!
TBranch *b_refphi; //!
TBranch *b_refy; //!
Expand All @@ -148,6 +154,8 @@ public :
TBranch *b_refarea; //!
TBranch *b_jte; //!
TBranch *b_jtpt; //!
TBranch *b_jtnMult; //!
TBranch *b_jtchMult; //!
TBranch *b_jteta; //!
TBranch *b_jtphi; //!
TBranch *b_jty; //!
Expand Down
17 changes: 15 additions & 2 deletions JetUtilities/src/HistogramUtilities.cc
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ namespace HistUtil {
b[i] = hist->GetBinLowEdge(i+1);
}
Double_t median = Median(nbins,x,y,b,0,debug);
// Set the error equal to the difference between the median and the bin boundary
/* // Set the error equal to the difference between the median and the bin boundary
// Choose the difference which is biggest
// We cannot determine the exact median from binned values
// The best we can do is a linear approximation of where inside the bin the median resides
Expand All @@ -126,14 +126,27 @@ namespace HistUtil {
Double_t width = hist->GetBinWidth(hist->FindBin(median));
Double_t high_edge = hist->GetBinLowEdge(hist->FindBin(median)) + width;
Double_t error = std::max(std::abs(median-low_edge),std::abs(high_edge-median));
*/

// The error on the median is non-trivial
// We use the "standard" relationship http://davidmlane.com/hyperstat/A106993.html
// for the error. But this only works if your distribution is ~normal.
// Fine binning is better to more precisely calculate the median & RMS
// EffectiveEntries is used to account for weights
Double_t error = 1.253 * hist->GetRMS() / TMath::Sqrt(hist->GetEffectiveEntries());

if (debug) {
cout << "\tlow_dege: " << low_edge << endl << "\thigh_edge: " << high_edge << endl
/* cout << "\tlow_dege: " << low_edge << endl << "\thigh_edge: " << high_edge << endl
<< "\twidth: " << width << endl << "\terror: " << error << endl;
cout << "\tdouble check error: " << std::max(std::abs(median-low_edge),std::abs(high_edge-median)) << endl
<< "\tabs(median-low_edge): " << std::abs(median-low_edge) << endl
<< "\tabs(high_edge-median): " << std::abs(high_edge-median) << endl
<< "\tmedian-low_edge: " << median-low_edge << endl
<< "\thigh_edge-median: " << high_edge-median << endl;
*/
cout << "\terror: " << error << endl;
cout << "\tRMS: " << hist->GetRMS() << endl
<< "\tentries: " << hist->GetEffectiveEntries() << endl;
}
delete [] x;
delete [] y;
Expand Down
Loading