Skip to content

Commit

Permalink
Merge pull request cms-sw#67 from mverwe/addJetIDVarsMerge
Browse files Browse the repository at this point in the history 
Add jet id variables for PF jets
  • Loading branch information
@mverwe
mverwe committed Mar 13, 2016
2 parents 28f9801 + b4dbd5b commit fbfcdab
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 60 deletions.
12 changes: 12 additions & 0 deletions HeavyIonsAnalysis/JetAnalysis/interface/HiInclusiveJetAnalyzer.h
View file
Expand Up @@ -207,6 +207,18 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer {
float jtm[MAXJETS];
float jtarea[MAXJETS];

float jtPfCHF[MAXJETS];
float jtPfNHF[MAXJETS];
float jtPfCEF[MAXJETS];
float jtPfNEF[MAXJETS];
float jtPfMUF[MAXJETS];

int jtPfCHM[MAXJETS];
int jtPfNHM[MAXJETS];
int jtPfCEM[MAXJETS];
int jtPfNEM[MAXJETS];
int jtPfMUM[MAXJETS];

float jttau1[MAXJETS];
float jttau2[MAXJETS];
float jttau3[MAXJETS];
Expand Down
71 changes: 30 additions & 41 deletions HeavyIonsAnalysis/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
View file
Expand Up @@ -214,6 +214,18 @@ HiInclusiveJetAnalyzer::beginJob() {
t->Branch("jtm",jets_.jtm,"jtm[nref]/F");
t->Branch("jtarea",jets_.jtarea,"jtarea[nref]/F");

t->Branch("jtPfCHF",jets_.jtPfCHF,"jtPfCHF[nref]/F");
t->Branch("jtPfNHF",jets_.jtPfNHF,"jtPfNHF[nref]/F");
t->Branch("jtPfCEF",jets_.jtPfCEF,"jtPfCEF[nref]/F");
t->Branch("jtPfNEF",jets_.jtPfNEF,"jtPfNEF[nref]/F");
t->Branch("jtPfMUF",jets_.jtPfMUF,"jtPfMUF[nref]/F");

t->Branch("jtPfCHM",jets_.jtPfCHM,"jtPfCHM[nref]/I");
t->Branch("jtPfNHM",jets_.jtPfNHM,"jtPfNHM[nref]/I");
t->Branch("jtPfCEM",jets_.jtPfCEM,"jtPfCEM[nref]/I");
t->Branch("jtPfNEM",jets_.jtPfNEM,"jtPfNEM[nref]/I");
t->Branch("jtPfMUM",jets_.jtPfMUM,"jtPfMUM[nref]/I");

t->Branch("jttau1",jets_.jttau1,"jttau1[nref]/F");
t->Branch("jttau2",jets_.jttau2,"jttau2[nref]/F");
t->Branch("jttau3",jets_.jttau3,"jttau3[nref]/F");
Expand Down Expand Up @@ -872,44 +884,7 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
}

// Calorimeter fractions

// for(unsigned int i = 0; i < hbheHits->size(); ++i){
// const HBHERecHit & hit= (*hbheHits)[i];
// math::XYZPoint pos = getPosition(hit.id(),vtx);
// double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
// if(dr < rParam){
// jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
// }
// }

// for(unsigned int i = 0; i < hfHits->size(); ++i){
// const HFRecHit & hit= (*hfHits)[i];
// math::XYZPoint pos = getPosition(hit.id(),vtx);
// double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
// if(dr < rParam){
// jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
// }
// }


// for(unsigned int i = 0; i < ebHits->size(); ++i){
// const EcalRecHit & hit= (*ebHits)[i];
// math::XYZPoint pos = getPosition(hit.id(),vtx);
// double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
// if(dr < rParam){
// jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
// }
// }

// for(unsigned int i = 0; i < eeHits->size(); ++i){
// const EcalRecHit & hit= (*eeHits)[i];
// math::XYZPoint pos = getPosition(hit.id(),vtx);
// double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
// if(dr < rParam){
// jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
// }
// }

// Jet ID for CaloJets
if(doTower){
// changing it to take things from towers
for(unsigned int i = 0; i < towers->size(); ++i){
Expand All @@ -925,7 +900,8 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
}

}
// Jet ID for CaloJets




if(doMatch_){
Expand Down Expand Up @@ -1030,7 +1006,7 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
jets_.jttau3[jets_.nref] = -999.;

if(doSubJets_) analyzeSubjets(jet);

if(usePat_){
if( (*patjets)[j].hasUserFloat(jetName_+"Njettiness:tau1") )
jets_.jttau1[jets_.nref] = (*patjets)[j].userFloat(jetName_+"Njettiness:tau1");
Expand All @@ -1039,6 +1015,20 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
if( (*patjets)[j].hasUserFloat(jetName_+"Njettiness:tau3") )
jets_.jttau3[jets_.nref] = (*patjets)[j].userFloat(jetName_+"Njettiness:tau3");

if( (*patjets)[j].isPFJet()) {
jets_.jtPfCHF[jets_.nref] = (*patjets)[j].chargedHadronEnergyFraction();
jets_.jtPfNHF[jets_.nref] = (*patjets)[j].neutralHadronEnergyFraction();
jets_.jtPfCEF[jets_.nref] = (*patjets)[j].chargedEmEnergyFraction();
jets_.jtPfNEF[jets_.nref] = (*patjets)[j].neutralEmEnergyFraction();
jets_.jtPfMUF[jets_.nref] = (*patjets)[j].muonEnergyFraction();

jets_.jtPfCHM[jets_.nref] = (*patjets)[j].chargedHadronMultiplicity();
jets_.jtPfNHM[jets_.nref] = (*patjets)[j].neutralHadronMultiplicity();
jets_.jtPfCEM[jets_.nref] = (*patjets)[j].electronMultiplicity();
jets_.jtPfNEM[jets_.nref] = (*patjets)[j].photonMultiplicity();
jets_.jtPfMUM[jets_.nref] = (*patjets)[j].muonMultiplicity();
}

if(doStandardJetID_){
jets_.fHPD[jets_.nref] = (*patjets)[j].jetID().fHPD;
jets_.fRBX[jets_.nref] = (*patjets)[j].jetID().fRBX;
Expand Down Expand Up @@ -1530,5 +1520,4 @@ void HiInclusiveJetAnalyzer::analyzeSubjets(const reco::Jet jet) {

}


DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);
30 changes: 11 additions & 19 deletions RecoJets/JetProducers/plugins/VirtualJetProducer.cc
View file
Expand Up @@ -878,7 +878,7 @@ void VirtualJetProducer::writeJetsWithConstituents( edm::Event & iEvent, edm::E

// get a list of output jets MV: make this compatible with template
std::auto_ptr<reco::PFJetCollection> jetCollection( new reco::PFJetCollection() );
// this is the mapping of constituent to jet
// this is the mapping of jet to constituents
std::vector< std::vector<int> > indices;
// this is the list of jet 4-momenta
std::vector<math::XYZTLorentzVector> p4_Jets;
Expand Down Expand Up @@ -914,40 +914,32 @@ void VirtualJetProducer::writeJetsWithConstituents( edm::Event & iEvent, edm::E
constituents = it->pieces();
else if ( it->has_constituents() )
fastjet::SelectorIsPureGhost().sift(it->constituents(), ghosts, constituents); //filter out ghosts

//loop over constituents of jet (can be subjets or normal constituents)
std::vector<fastjet::PseudoJet>::const_iterator itConstBegin = constituents.begin(),
itConst = itConstBegin, itConstEnd = constituents.end();
for (; itConst != itConstEnd; ++itConst ) {
fastjet::PseudoJet const & constit = *itConst;

if ( verbosity_ >= 1 ) {
std::cout << "jet #" << jetIndex << " constituent #" << (itConst - itConstBegin) << ": Pt = " << constit.pt() << ", eta = " << constit.eta() << ", phi = " << constit.phi() << ", mass = " << constit.m() << ", uid: " << constit.user_index() << ")" << std::endl;
std::cout << "jet #" << jetIndex << " constituent #" << (itConst - itConstBegin) << ": Pt = " << constit.pt() << ", eta = " << constit.eta() << ", phi = " << constit.phi() << ", mass = " << constit.m() << ", uid: " << constit.user_index() << ", pos: " << constituentsSub.size() << ")" << std::endl;
}

indices[jetIndex].push_back( constituentsSub.size() );
constituentsSub.push_back(constit);
}
}

//Loop over original constituents and find constituent subtracted version
//Loop over constituents and store in the event
static const reco::PFCandidate dummySinceTranslateIsNotStatic;
for (unsigned int i=0;i<inputs_.size();i++) {

auto subMatched = find_if( constituentsSub.begin(), constituentsSub.end(), [&i]( fastjet::PseudoJet const & p ){ return p.user_index() == (int)i; } );
auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(inputs_[i]->pdgId());
reco::PFCandidate pCand( reco::PFCandidate(inputs_[i]->charge(), inputs_[i]->p4(), id) );

for (std::vector<fastjet::PseudoJet>::const_iterator itsub = constituentsSub.begin() ; itsub != constituentsSub.end(); ++itsub ) {
fastjet::PseudoJet const & constit = *itsub;
auto orig = inputs_[constit.user_index()];
auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(orig->pdgId());
reco::PFCandidate pCand( reco::PFCandidate(orig->charge(), orig->p4(), id) );
math::XYZTLorentzVector pVec;
if ( subMatched != constituentsSub.end() )
pVec.SetPxPyPzE(subMatched->px(),subMatched->py(),subMatched->pz(),subMatched->e());
else
pVec.SetPxPyPzE( 0, 0, 0, 0);
pVec.SetPxPyPzE(constit.px(),constit.py(),constit.pz(),constit.e());
pCand.setP4(pVec);
pCand.setSourceCandidatePtr( inputs_[i]->sourceCandidatePtr(0) );
pCand.setSourceCandidatePtr( orig->sourceCandidatePtr(0) );
constituentCollection->push_back(pCand);
}

// put constituents into event record
constituentHandleAfterPut = iEvent.put( constituentCollection, jetCollInstanceName_ );

Expand Down

0 comments on commit fbfcdab

Please sign in to comment.