ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc

Comparing UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc (file contents):
Revision 1.4 by yjlee, Wed Nov 2 19:27:48 2011 UTC vs.
Revision 1.13 by yilmaz, Thu May 10 10:47:17 2012 UTC

#	Line 21 | Line 21
21		#include "DataFormats/HeavyIonEvent/interface/Centrality.h"
22
23		#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
24	–	#include "DataFormats/PatCandidates/interface/Jet.h"
24		#include "DataFormats/JetReco/interface/CaloJetCollection.h"
25		#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
26		#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
27	+	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
28		#include "DataFormats/JetReco/interface/GenJetCollection.h"
29		#include "DataFormats/VertexReco/interface/Vertex.h"
30		#include "DataFormats/Math/interface/deltaPhi.h"
#	Line 46 | Line 46 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
46
47
48		jetTag_ = iConfig.getParameter<InputTag>("jetTag");
49	<	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
49	>	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
50
51		isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
52	<
52	>
53		if(isMC_){
54		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
55		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 61 | Line 61 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
61		useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
62		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
63
64	+	doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
65	+	doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
66	+	saveBfragments_  = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
67	+
68	+	pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
69	+
70		if(!isMC_){
71		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
72		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
#	Line 78 | Line 84 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
84		}
85
86		cout<<" jet collection : "<<jetTag_<<endl;
87	<	if(isMC_)cout<<" genjet collection : "<<genjetTag_<<endl;
87	>	doSubEvent_ = 0;
88
89	+	if(isMC_){
90	+	cout<<" genjet collection : "<<genjetTag_<<endl;
91	+	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
92	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
93	+	}
94
95
96		}
#	Line 115 | Line 126 | HiInclusiveJetAnalyzer::beginJob() {
126		}
127		if (useCentrality_) {
128		t->Branch("hf",&jets_.hf,"hf/F");
118	–	t->Branch("nref",&jets_.nref,"nref/I");
129		t->Branch("bin",&jets_.bin,"bin/I");
130		}
131	+
132	+	t->Branch("nref",&jets_.nref,"nref/I");
133		t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
134		t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
135		t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
#	Line 125 | Line 137 | HiInclusiveJetAnalyzer::beginJob() {
137		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
138		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
139
140	+	// b-jet discriminators
141	+	if (doLifeTimeTagging_) {
142	+
143	+	t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
144	+	t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
145	+	t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
146	+	t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
147	+	t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
148	+	t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
149	+	t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
150	+	t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
151	+
152	+	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
153	+	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
154	+	t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
155	+	t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
156	+	t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
157	+	t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
158	+
159	+	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
160	+	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
161	+
162	+	if (doLifeTimeTaggingExtras_) {
163	+	t->Branch("nIP",&jets_.nIP,"nIP/I");
164	+	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
165	+	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
166	+	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
167	+	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
168	+	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
169	+	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
170	+	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
171	+	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
172	+	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
173	+	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
174	+	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
175	+
176	+	}
177	+
178	+	t->Branch("mue",     jets_.mue,     "mue[nref]/F");
179	+	t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
180	+	t->Branch("mueta",   jets_.mueta,   "mueta[nref]/F");
181	+	t->Branch("muphi",   jets_.muphi,   "muphi[nref]/F");
182	+	t->Branch("mudr",    jets_.mudr,    "mudr[nref]/F");
183	+	t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
184	+	t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
185	+	}
186	+
187		if(isMC_){
188	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
189	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
190	+
191		t->Branch("pthat",&jets_.pthat,"pthat/F");
192
193		// Only matched gen jets
#	Line 137 | Line 199 | HiInclusiveJetAnalyzer::beginJob() {
199		t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
200		// matched parton
201		t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
202	<	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
202	>	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
203	>	t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
204
205		// For all gen jets, matched or unmatched
206		t->Branch("ngen",&jets_.ngen,"ngen/I");
#	Line 148 | Line 211 | HiInclusiveJetAnalyzer::beginJob() {
211		t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
212		t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
213		t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
214	+
215	+	if(doSubEvent_){
216	+	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
217	+	}
218	+
219	+	if(saveBfragments_  ) {
220	+	t->Branch("bMult",&jets_.bMult,"bMult/I");
221	+	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
222	+	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
223	+	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
224	+	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
225	+	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
226	+	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
227	+	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
228	+	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
229	+	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
230	+	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
231	+	}
232	+
233		}
234	<
234	>	/*
235		if(!isMC_){
236		t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
237		t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
#	Line 161 | Line 243 | HiInclusiveJetAnalyzer::beginJob() {
243		t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
244
245		}
246	<
246	>	*/
247		TH1D::SetDefaultSumw2();
248
249
#	Line 201 | Line 283 | HiInclusiveJetAnalyzer::analyze(const Ev
283
284
285
204	–	//double jetPtMin = 35.;
205	–
286
287		// loop the events
288
#	Line 227 | Line 307 | HiInclusiveJetAnalyzer::analyze(const Ev
307		edm::Handle<reco::JetView> jets;
308		iEvent.getByLabel(jetTag_, jets);
309
310	+	edm::Handle<reco::PFCandidateCollection> pfCandidates;
311	+	if(doLifeTimeTagging_){
312	+	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
313	+	}
314		// FILL JRA TREE
315
316		jets_.b = b;
#	Line 240 | Line 324 | HiInclusiveJetAnalyzer::analyze(const Ev
324		for(unsigned int j = 0 ; j < jets->size(); ++j){
325		const reco::Jet& jet = (*jets)[j];
326
327	<	//cout<<" jet pt "<<jet.pt()<<endl;
244	<	//if(jet.pt() < jetPtMin) continue;
327	>
328		if (useJEC_ && usePat_){
329		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
330		}
331	+
332	+	if(doLifeTimeTagging_){
333	+
334	+	if(jetTag_.label()=="icPu5patJets"){
335	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
336	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
337	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
338	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
339	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
340	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
341	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
342	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
343	+	}
344	+	else if(jetTag_.label()=="akPu3PFpatJets"){
345	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
346	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
347	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
348	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
349	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
350	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
351	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
352	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
353	+	}
354	+	else{
355	+	cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
356	+	}
357	+
358	+	const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
359	+
360	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
361	+
362	+	if (tagInfoSV.nVertices()>0) {
363	+	jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
364	+	// this is the 3d flight distance, for 2-D use (0,true)
365	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
366	+	jets_.svtxdl[jets_.nref]    = m1D.value();
367	+	jets_.svtxdls[jets_.nref]   = m1D.significance();
368	+
369	+	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
370	+	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
371	+	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
372	+	jets_.svtxm[jets_.nref]    = svtx.p4().mass();
373	+	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
374	+	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
375	+	}
376	+
377	+	const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
378	+
379	+	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
380	+	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
381	+
382	+	if (doLifeTimeTaggingExtras_) {
383	+
384	+	TrackRefVector selTracks=tagInfoIP.selectedTracks();
385	+
386	+	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
387	+
388	+	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
389	+	{
390	+	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
391	+	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
392	+	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
393	+	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
394	+	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
395	+	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
396	+	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
397	+	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
398	+	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
399	+	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
400	+	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
401	+	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
402	+	}
403	+
404	+	jets_.nIP += jets_.nselIPtrk[jets_.nref];
405	+
406	+	}
407	+
408	+	const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
409	+	int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
410	+	if(pfMuonIndex >=0){
411	+	const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
412	+	jets_.mupt[jets_.nref]    =  muon.pt();
413	+	jets_.mueta[jets_.nref]   =  muon.eta();
414	+	jets_.muphi[jets_.nref]   =  muon.phi();
415	+	jets_.mue[jets_.nref]     =  muon.energy();
416	+	jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
417	+	jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
418	+	jets_.muchg[jets_.nref]   =  muon.charge();
419	+	}
420	+	}
421	+
422		jets_.jtpt[jets_.nref] = jet.pt();
423		jets_.jteta[jets_.nref] = jet.eta();
424		jets_.jtphi[jets_.nref] = jet.phi();
#	Line 252 | Line 426 | HiInclusiveJetAnalyzer::analyze(const Ev
426		jets_.jtpu[jets_.nref] = jet.pileup();
427
428		if(isMC_ && usePat_){
429	<	const reco::GenJet* genjet = (*patjets)[j].genJet();
429	>
430	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
431	>
432		if(genjet){
433		jets_.refpt[jets_.nref] = genjet->pt();
434		jets_.refeta[jets_.nref] = genjet->eta();
#	Line 269 | Line 445 | HiInclusiveJetAnalyzer::analyze(const Ev
445		jets_.refdrjt[jets_.nref] = -999.;
446		}
447
448	<	// matched partons
449	<	const reco::GenParticle * parton = (*patjets)[j].genParton();
450	<	if(parton){
451	<	jets_.refparton_pt[jets_.nref] = parton->pt();
452	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
453	<	} else {
454	<	jets_.refparton_pt[jets_.nref] = -999;
455	<	jets_.refparton_flavor[jets_.nref] = -999;
456	<	}
448	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
449	>
450	>	// matched partons
451	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
452	>
453	>	if((*patjets)[j].genParton()){
454	>	jets_.refparton_pt[jets_.nref] = parton.pt();
455	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
456	>
457	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
458	>
459	>	usedStringPts.clear();
460	>
461	>	// uncomment this if you want to know the ugly truth about parton matching -matt
462	>	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
463	>	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
464	>
465	>	jets_.bJetIndex[jets_.bMult] = jets_.nref;
466	>	jets_.bStatus[jets_.bMult] = parton.status();
467	>	jets_.bVx[jets_.bMult] = parton.vx();
468	>	jets_.bVy[jets_.bMult] = parton.vy();
469	>	jets_.bVz[jets_.bMult] = parton.vz();
470	>	jets_.bPt[jets_.bMult] = parton.pt();
471	>	jets_.bEta[jets_.bMult] = parton.eta();
472	>	jets_.bPhi[jets_.bMult] = parton.phi();
473	>	jets_.bPdg[jets_.bMult] = parton.pdgId();
474	>	jets_.bChg[jets_.bMult] = parton.charge();
475	>	jets_.bMult++;
476	>	saveDaughters(parton);
477	>	}
478	>
479	>
480	>	} else {
481	>	jets_.refparton_pt[jets_.nref] = -999;
482	>	jets_.refparton_flavor[jets_.nref] = -999;
483	>	}
484	>
485	>
486		}
487
488		jets_.nref++;
#	Line 288 | Line 493 | HiInclusiveJetAnalyzer::analyze(const Ev
493
494		if(isMC_){
495
496	+	edm::Handle<HepMCProduct> hepMCProduct;
497	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
498	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
499	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
500	+	jets_.beamId1 = beamParticles.first->pdg_id();
501	+	jets_.beamId2 = beamParticles.second->pdg_id();
502	+
503		edm::Handle<GenEventInfoProduct> hEventInfo;
504		iEvent.getByLabel(eventInfoTag_,hEventInfo);
505		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 305 | Line 517 | HiInclusiveJetAnalyzer::analyze(const Ev
517		float genjet_pt = genjet.pt();
518
519		// threshold to reduce size of output in minbias PbPb
520	<	if(genjet_pt>20.){
520	>	if(genjet_pt>genPtMin_){
521
522		jets_.genpt[jets_.ngen] = genjet_pt;
523		jets_.geneta[jets_.ngen] = genjet.eta();
524		jets_.genphi[jets_.ngen] = genjet.phi();
525		jets_.geny[jets_.ngen] = genjet.eta();
526
527	+	if(doSubEvent_){
528	+	const GenParticle* gencon = genjet.getGenConstituent(0);
529	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
530	+	}
531
532		// find matching patJet if there is one
533
534		jets_.gendrjt[jets_.ngen] = -1.0;
535		jets_.genmatchindex[jets_.ngen] = -1;
536
321	–
537		for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
538		const pat::Jet& jet = (*jets)[ijet];
539	<
540	<	if(jet.genJet()){
541	<	if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
542	<	fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
543	<	(fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
539	>	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
540	>	if(matchedGenJet){
541	>	// poor man's matching, someone fix please
542	>	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
543	>	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
544	>	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
545
546		jets_.genmatchindex[jets_.ngen] = (int)ijet;
547		jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
#	Line 342 | Line 558 | HiInclusiveJetAnalyzer::analyze(const Ev
558
559		}
560		t->Fill();
561	+	memset(&jets_,0,sizeof jets_);
562		}
563
564
#	Line 441 | Line 658 | inline bool HiInclusiveJetAnalyzer::getP
658		}
659
660
661	+	int
662	+	HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
663	+	{
664	+
665	+	int pfMuonIndex = -1;
666	+	float ptMax = 0.;
667	+
668	+
669	+	for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
670	+	const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
671	+
672	+	int id = pfCandidate.particleId();
673	+	if(abs(id) != 3) continue;
674	+
675	+	if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
676	+
677	+	double pt =  pfCandidate.pt();
678	+	if(pt>ptMax){
679	+	ptMax = pt;
680	+	pfMuonIndex = (int) icand;
681	+	}
682	+	}
683	+
684	+	return pfMuonIndex;
685	+
686	+	}
687	+
688	+
689	+	double
690	+	HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
691	+	{
692	+
693	+	float lj_x = jet.p4().px();
694	+	float lj_y = jet.p4().py();
695	+	float lj_z = jet.p4().pz();
696	+
697	+	// absolute values squared
698	+	float lj2  = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
699	+	float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
700	+
701	+	// projection vec(mu) to lepjet axis
702	+	float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
703	+
704	+	// absolute value squared and normalized
705	+	float pLrel2 = lepXlj*lepXlj/lj2;
706	+
707	+	// lep2 = pTrel2 + pLrel2
708	+	float pTrel2 = lep2-pLrel2;
709	+
710	+	return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
711	+	}
712	+
713	+	// Recursive function, but this version gets called only the first time
714	+	void
715	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
716	+
717	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
718	+	const reco::Candidate & daughter = *gen.daughter(i);
719	+	double daughterPt = daughter.pt();
720	+	if(daughterPt<1.) continue;
721	+	double daughterEta = daughter.eta();
722	+	if(fabs(daughterEta)>3.) continue;
723	+	int daughterPdgId = daughter.pdgId();
724	+	int daughterStatus = daughter.status();
725	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
726	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
727	+
728	+	// cheesy way of finding strings which were already used
729	+	if(daughter.pdgId()==92){
730	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
731	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
732	+	}
733	+	usedStringPts.push_back(daughter.pt());
734	+	}
735	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
736	+	jets_.bStatus[jets_.bMult] = daughterStatus;
737	+	jets_.bVx[jets_.bMult] = daughter.vx();
738	+	jets_.bVy[jets_.bMult] = daughter.vy();
739	+	jets_.bVz[jets_.bMult] = daughter.vz();
740	+	jets_.bPt[jets_.bMult] = daughterPt;
741	+	jets_.bEta[jets_.bMult] = daughterEta;
742	+	jets_.bPhi[jets_.bMult] = daughter.phi();
743	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
744	+	jets_.bChg[jets_.bMult] = daughter.charge();
745	+	jets_.bMult++;
746	+	saveDaughters(daughter);
747	+	}
748	+	}
749	+
750	+	// This version called for all subsequent calls
751	+	void
752	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
753	+
754	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
755	+	const reco::Candidate & daughter = *gen.daughter(i);
756	+	double daughterPt = daughter.pt();
757	+	if(daughterPt<1.) continue;
758	+	double daughterEta = daughter.eta();
759	+	if(fabs(daughterEta)>3.) continue;
760	+	int daughterPdgId = daughter.pdgId();
761	+	int daughterStatus = daughter.status();
762	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
763	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
764	+
765	+	// cheesy way of finding strings which were already used
766	+	if(daughter.pdgId()==92){
767	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
768	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
769	+	}
770	+	usedStringPts.push_back(daughter.pt());
771	+	}
772
773	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
774	+	jets_.bStatus[jets_.bMult] = daughterStatus;
775	+	jets_.bVx[jets_.bMult] = daughter.vx();
776	+	jets_.bVy[jets_.bMult] = daughter.vy();
777	+	jets_.bVz[jets_.bMult] = daughter.vz();
778	+	jets_.bPt[jets_.bMult] = daughterPt;
779	+	jets_.bEta[jets_.bMult] = daughterEta;
780	+	jets_.bPhi[jets_.bMult] = daughter.phi();
781	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
782	+	jets_.bChg[jets_.bMult] = daughter.charge();
783	+	jets_.bMult++;
784	+	saveDaughters(daughter);
785	+	}
786	+	}
787
788		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines