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.9 by yilmaz, Sun Apr 22 19:12:44 2012 UTC vs.
Revision 1.19 by yilmaz, Fri May 11 19:07:29 2012 UTC

#	Line 24 | Line 24
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 42 | Line 43 | using namespace edm;
43		using namespace reco;
44
45		HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
45	–
46
47		jetTag_ = iConfig.getParameter<InputTag>("jetTag");
48	+	matchTag_ = iConfig.getUntrackedParameter<InputTag>("matchTag",jetTag_);
49	+
50		vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
51	+	trackTag_ = iConfig.getParameter<InputTag>("trackTag");
52	+	useQuality_ = iConfig.getUntrackedParameter<bool>("useQuality",0);
53	+	string trackQuality_ = iConfig.getUntrackedParameter<string>("trackQuality","highPurity");
54
55		isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
56	<
56	>	doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
57	>
58	>	rParam = iConfig.getParameter<double>("rParam");
59	>	hardPtMin_ = iConfig.getUntrackedParameter<double>("hardPtMin",4);
60	>
61		if(isMC_){
62		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
63		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 61 | Line 70 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
70		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
71
72		doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
73	+	doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
74	+	saveBfragments_  = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
75	+
76		pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
77
78	<	if(!isMC_){
78	>	if(doTrigger_){
79		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
80		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
81
#	Line 80 | Line 92 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
92		}
93
94		cout<<" jet collection : "<<jetTag_<<endl;
95	+	doSubEvent_ = 0;
96	+
97		if(isMC_){
98		cout<<" genjet collection : "<<genjetTag_<<endl;
99		genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
100	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
101		}
102
103
#	Line 130 | Line 145 | HiInclusiveJetAnalyzer::beginJob() {
145		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
146		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
147
148	+	// jet ID information, jet composition
149	+	t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
150	+
151	+	t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
152	+	t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
153	+	t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
154	+	t->Branch("trackHardSum", jets_.trackHardSum,"trackHardSum[nref]/F");
155	+	t->Branch("trackHardN", jets_.trackHardN,"trackHardN[nref]/I");
156	+
157	+	t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
158	+	t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
159	+	t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
160	+	t->Branch("chargedHardSum", jets_.chargedHardSum,"chargedHardSum[nref]/F");
161	+	t->Branch("chargedHardN", jets_.chargedHardN,"chargedHardN[nref]/I");
162	+
163	+	t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
164	+	t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
165	+	t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
166	+	t->Branch("photonHardSum", jets_.photonHardSum,"photonHardSum[nref]/F");
167	+	t->Branch("photonHardN", jets_.photonHardN,"photonHardN[nref]/I");
168	+
169	+	t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
170	+	t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
171	+	t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
172	+
173	+	t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
174	+	t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
175	+	t->Branch("eN", jets_.eN,"eN[nref]/I");
176	+
177	+	t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
178	+	t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
179	+	t->Branch("muN", jets_.muN,"muN[nref]/I");
180	+
181	+	t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
182	+	t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
183	+
184		// b-jet discriminators
185		if (doLifeTimeTagging_) {
186
#	Line 142 | Line 193 | HiInclusiveJetAnalyzer::beginJob() {
193		t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
194		t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
195
196	<	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/b");
197	<	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/b");
196	>	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
197	>	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
198		t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
199		t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
200		t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
201		t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
202
203	<	t->Branch("nIPtracks",jets_.nIPtracks,"nIPTracks[nref]/b");
204	<	t->Branch("nselIPtracks",jets_.nselIPtracks,"nselIPTracks[nref]/b");
203	>	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
204	>	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
205	>
206	>	if (doLifeTimeTaggingExtras_) {
207	>	t->Branch("nIP",&jets_.nIP,"nIP/I");
208	>	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
209	>	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
210	>	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
211	>	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
212	>	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
213	>	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
214	>	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
215	>	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
216	>	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
217	>	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
218	>	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
219	>
220	>	}
221
222		t->Branch("mue",     jets_.mue,     "mue[nref]/F");
223		t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
#	Line 160 | Line 227 | HiInclusiveJetAnalyzer::beginJob() {
227		t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
228		t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
229		}
230	+
231
232		if(isMC_){
233	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
234	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
235	+
236		t->Branch("pthat",&jets_.pthat,"pthat/F");
237
238		// Only matched gen jets
#	Line 185 | Line 256 | HiInclusiveJetAnalyzer::beginJob() {
256		t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
257		t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
258		t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
259	+
260	+	if(doSubEvent_){
261	+	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
262	+	}
263	+
264	+	if(saveBfragments_  ) {
265	+	t->Branch("bMult",&jets_.bMult,"bMult/I");
266	+	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
267	+	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
268	+	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
269	+	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
270	+	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
271	+	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
272	+	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
273	+	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
274	+	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
275	+	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
276	+	}
277	+
278		}
279		/*
280		if(!isMC_){
#	Line 238 | Line 328 | HiInclusiveJetAnalyzer::analyze(const Ev
328
329
330
241	–	//double jetPtMin = 35.;
242	–
331
332		// loop the events
333
#	Line 260 | Line 348 | HiInclusiveJetAnalyzer::analyze(const Ev
348
349		edm::Handle<pat::JetCollection> patjets;
350		if(usePat_)iEvent.getByLabel(jetTag_, patjets);
351	+
352	+	edm::Handle<pat::JetCollection> matchedjets;
353	+	iEvent.getByLabel(matchTag_, matchedjets);
354
355		edm::Handle<reco::JetView> jets;
356		iEvent.getByLabel(jetTag_, jets);
357
358		edm::Handle<reco::PFCandidateCollection> pfCandidates;
359	<	if(doLifeTimeTagging_){
360	<	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
361	<	}
359	>	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
360	>
361	>	edm::Handle<reco::TrackCollection> tracks;
362	>	iEvent.getByLabel(trackTag_,tracks);
363	>
364		// FILL JRA TREE
365
366		jets_.b = b;
367		jets_.nref = 0;
368
369	<	if(!isMC_){
369	>	if(doTrigger_){
370		fillL1Bits(iEvent);
371		fillHLTBits(iEvent);
372		}
#	Line 281 | Line 374 | HiInclusiveJetAnalyzer::analyze(const Ev
374		for(unsigned int j = 0 ; j < jets->size(); ++j){
375		const reco::Jet& jet = (*jets)[j];
376
377	<	//cout<<" jet pt "<<jet.pt()<<endl;
285	<	//if(jet.pt() < jetPtMin) continue;
377	>
378		if (useJEC_ && usePat_){
379		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
380		}
#	Line 310 | Line 402 | HiInclusiveJetAnalyzer::analyze(const Ev
402		jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
403		}
404		else{
405	<	cout<<" you fail at btagging "<<endl;
405	>	cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
406		}
407	<
407	>
408		const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
409	+
410	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
411	+
412		if (tagInfoSV.nVertices()>0) {
318	–	Measurement1D m1D = tagInfoSV.flightDistance(0);
319	–	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
413		jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
414	+	// this is the 3d flight distance, for 2-D use (0,true)
415	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
416		jets_.svtxdl[jets_.nref]    = m1D.value();
417		jets_.svtxdls[jets_.nref]   = m1D.significance();
418
419	<	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
420	<
419	>	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
420	>	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
421	>	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
422		jets_.svtxm[jets_.nref]    = svtx.p4().mass();
423	<	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
424	<
329	<	}
330	<	else {
331	<	jets_.nsvtx[jets_.nref]    =   0;
332	<	jets_.svtxntrk[jets_.nref] =   0;
333	<	jets_.svtxdl[jets_.nref]   = 0.0;
334	<	jets_.svtxdls[jets_.nref]  = 0.0;
335	<	jets_.svtxm[jets_.nref]    = 0.0;
336	<	jets_.svtxpt[jets_.nref]   = 0.0;
423	>	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
424	>	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
425		}
426	+
427		const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
428	+
429	+	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
430	+	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
431
432	<	jets_.nIPtracks[jets_.nref] = tagInfoIP.tracks().size();
341	<	jets_.nselIPtracks[jets_.nref] = tagInfoIP.selectedTracks().size();
432	>	if (doLifeTimeTaggingExtras_) {
433
434	<	// would be nice to save the info for the tracks, but requires a more sophisticated tree
435	<	/*
436	<	//cout << "Jet pt: " << tagInfoIP.jet()->pt() << endl;
437	<	cout << "Tot tracks: " << tagInfoIP.tracks().size() << endl;
438	<	TrackRefVector selTracks=tagInfoIP.selectedTracks();
439	<	int n=selTracks.size();
440	<	cout << "Sel tracks: " << n << endl;
441	<	// false      cout << " Pt  \t d len \t jet dist \t p3d \t p2d\t ip3d \t ip2d " << endl;
442	<	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
443	<	cout << pv << " vs " << tagInfoIP.primaryVertex()->position()   << endl;
444	<	for(int j=0;j< n;j++)
445	<	{
446	<	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[j];
447	<	cout << selTracks[j]->pt() << "\t";
448	<	cout << tagInfoIP.probabilities(0)[j] << "\t";
449	<	cout << tagInfoIP.probabilities(1)[j] << "\t";
450	<	cout << data.ip3d.value() << "\t";
451	<	cout << data.ip3d.significance() << "\t";
452	<	cout << data.distanceToJetAxis.value() << "\t";
453	<	cout << data.distanceToJetAxis.significance() << "\t";
454	<	cout << data.distanceToGhostTrack.value() << "\t";
364	<	cout << data.distanceToGhostTrack.significance() << "\t";
365	<	cout << data.closestToJetAxis << "\t";
366	<	cout << (data.closestToJetAxis - pv).mag() << "\t";
367	<	cout << data.closestToGhostTrack << "\t";
368	<	cout << (data.closestToGhostTrack - pv).mag() << "\t";
369	<	cout << data.ip2d.value() << "\t";
370	<	cout << data.ip2d.significance() <<  endl;
371	<	}
372	<	*/
434	>	TrackRefVector selTracks=tagInfoIP.selectedTracks();
435	>
436	>	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
437	>
438	>	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
439	>	{
440	>	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
441	>	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
442	>	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
443	>	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
444	>	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
445	>	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
446	>	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
447	>	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
448	>	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
449	>	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
450	>	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
451	>	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
452	>	}
453	>
454	>	jets_.nIP += jets_.nselIPtrk[jets_.nref];
455
456	+	}
457	+
458		const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
459		int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
460		if(pfMuonIndex >=0){
#	Line 382 | Line 466 | HiInclusiveJetAnalyzer::analyze(const Ev
466		jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
467		jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
468		jets_.muchg[jets_.nref]   =  muon.charge();
469	<	}
386	<	else{
469	>	}else{
470		jets_.mupt[jets_.nref]    =  0.0;
471		jets_.mueta[jets_.nref]   =  0.0;
472		jets_.muphi[jets_.nref]   =  0.0;
#	Line 392 | Line 475 | HiInclusiveJetAnalyzer::analyze(const Ev
475		jets_.muptrel[jets_.nref] =  0.0;
476		jets_.muchg[jets_.nref]   = 0;
477		}
478	+
479		}
480
481	+
482	+
483	+	// Jet ID variables
484	+
485	+	jets_.muMax[jets_.nref] = 0;
486	+	jets_.muSum[jets_.nref] = 0;
487	+	jets_.muN[jets_.nref] = 0;
488	+
489	+	jets_.eMax[jets_.nref] = 0;
490	+	jets_.eSum[jets_.nref] = 0;
491	+	jets_.eN[jets_.nref] = 0;
492	+
493	+	jets_.neutralMax[jets_.nref] = 0;
494	+	jets_.neutralSum[jets_.nref] = 0;
495	+	jets_.neutralN[jets_.nref] = 0;
496	+
497	+	jets_.photonMax[jets_.nref] = 0;
498	+	jets_.photonSum[jets_.nref] = 0;
499	+	jets_.photonN[jets_.nref] = 0;
500	+	jets_.photonHardSum[jets_.nref] = 0;
501	+	jets_.photonHardN[jets_.nref] = 0;
502	+
503	+	jets_.chargedMax[jets_.nref] = 0;
504	+	jets_.chargedSum[jets_.nref] = 0;
505	+	jets_.chargedN[jets_.nref] = 0;
506	+	jets_.chargedHardSum[jets_.nref] = 0;
507	+	jets_.chargedHardN[jets_.nref] = 0;
508	+
509	+	jets_.trackMax[jets_.nref] = 0;
510	+	jets_.trackSum[jets_.nref] = 0;
511	+	jets_.trackN[jets_.nref] = 0;
512	+	jets_.trackHardSum[jets_.nref] = 0;
513	+	jets_.trackHardN[jets_.nref] = 0;
514	+
515	+
516	+	for(unsigned int icand = 0; icand < tracks->size(); ++icand){
517	+	const reco::Track& track = (*tracks)[icand];
518	+	if(useQuality_ && !(track.quality(reco::TrackBase::qualityByName(trackQuality_)))) continue;
519	+
520	+	double dr = deltaR(jet,track);
521	+	if(dr < rParam){
522	+	double ptcand = track.pt();
523	+	jets_.trackSum[jets_.nref] += ptcand;
524	+	jets_.trackN[jets_.nref] += 1;
525	+
526	+	if(ptcand > hardPtMin_){
527	+	jets_.trackHardSum[jets_.nref] += ptcand;
528	+	jets_.trackHardN[jets_.nref] += 1;
529	+
530	+	}
531	+	if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
532	+
533	+	}
534	+	}
535	+
536	+	for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
537	+	const reco::PFCandidate& track = (*pfCandidates)[icand];
538	+	double dr = deltaR(jet,track);
539	+	if(dr < rParam){
540	+	double ptcand = track.pt();
541	+	int pfid = track.particleId();
542	+
543	+	switch(pfid){
544	+
545	+	case 1:
546	+	jets_.chargedSum[jets_.nref] += ptcand;
547	+	jets_.chargedN[jets_.nref] += 1;
548	+	if(ptcand > hardPtMin_){
549	+	jets_.chargedHardSum[jets_.nref] += ptcand;
550	+	jets_.chargedHardN[jets_.nref] += 1;
551	+	}
552	+	if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
553	+	break;
554	+
555	+	case 2:
556	+	jets_.eSum[jets_.nref] += ptcand;
557	+	jets_.eN[jets_.nref] += 1;
558	+	if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
559	+	break;
560	+
561	+	case 3:
562	+	jets_.muSum[jets_.nref] += ptcand;
563	+	jets_.muN[jets_.nref] += 1;
564	+	if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
565	+	break;
566	+
567	+	case 4:
568	+	jets_.photonSum[jets_.nref] += ptcand;
569	+	jets_.photonN[jets_.nref] += 1;
570	+	if(ptcand > hardPtMin_){
571	+	jets_.photonHardSum[jets_.nref] += ptcand;
572	+	jets_.photonHardN[jets_.nref] += 1;
573	+	}
574	+	if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
575	+	break;
576	+
577	+	case 5:
578	+	jets_.neutralSum[jets_.nref] += ptcand;
579	+	jets_.neutralN[jets_.nref] += 1;
580	+	if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
581	+	break;
582	+
583	+	}
584	+	}
585	+	}
586	+
587	+	double drMin = 100;
588	+	for(unsigned int j = 0 ; j < matchedjets->size(); ++j){
589	+	const reco::Jet& mjet = (*matchedjets)[j];
590	+
591	+	double dr = deltaR(jet,mjet);
592	+	if(dr < drMin){
593	+	jets_.matchedPt[jets_.nref] = mjet.pt();
594	+	jets_.matchedR[jets_.nref] = dr;
595	+
596	+	}
597	+	}
598	+
599	+
600	+
601	+
602	+	//     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
603	+
604	+	/////////////////////////////////////////////////////////////////
605	+	// Jet core pt^2 discriminant for fake jets
606	+	// Edited by Yue Shi Lai <ylai@mit.edu>
607	+
608	+	// Initial value is 0
609	+	jets_.discr_fr01[jets_.nref] = 0;
610	+	// Start with no directional adaption, i.e. the fake rejection
611	+	// axis is the jet axis
612	+	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
613	+	float azimuth_adapt = jets_.jtphi[jets_.nref];
614	+
615	+	// Unadapted discriminant with adaption search
616	+	for (size_t iteration = 0; iteration < 2; iteration++) {
617	+	float pseudorapidity_adapt_new = pseudorapidity_adapt;
618	+	float azimuth_adapt_new = azimuth_adapt;
619	+	float max_weighted_perp = 0;
620	+	float perp_square_sum = 0;
621	+
622	+	for (size_t index_pf_candidate = 0;
623	+	index_pf_candidate < pfCandidates->size();
624	+	index_pf_candidate++) {
625	+	const reco::PFCandidate &p =
626	+	(*pfCandidates)[index_pf_candidate];
627	+
628	+	switch (p.particleId()) {
629	+	case 1: // Charged hadron
630	+	case 3: // Muon
631	+	case 4: // Photon
632	+	{
633	+	const float dpseudorapidity =
634	+	p.eta() - pseudorapidity_adapt;
635	+	const float dazimuth =
636	+	reco::deltaPhi(p.phi(), azimuth_adapt);
637	+	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
638	+	// = 50
639	+	const float angular_weight =
640	+	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
641	+	dazimuth * dazimuth));
642	+	const float weighted_perp =
643	+	angular_weight * p.pt() * p.pt();
644	+	const float weighted_perp_square =
645	+	weighted_perp * p.pt();
646	+
647	+	perp_square_sum += weighted_perp_square;
648	+	if (weighted_perp >= max_weighted_perp) {
649	+	pseudorapidity_adapt_new = p.eta();
650	+	azimuth_adapt_new = p.phi();
651	+	max_weighted_perp = weighted_perp;
652	+	}
653	+	}
654	+	}
655	+	}
656	+	// Update the fake rejection value
657	+	jets_.discr_fr01[jets_.nref] = std::max(
658	+	jets_.discr_fr01[jets_.nref], perp_square_sum);
659	+	// Update the directional adaption
660	+	pseudorapidity_adapt = pseudorapidity_adapt_new;
661	+	azimuth_adapt = azimuth_adapt_new;
662	+	}
663	+
664	+
665		jets_.jtpt[jets_.nref] = jet.pt();
666		jets_.jteta[jets_.nref] = jet.eta();
667		jets_.jtphi[jets_.nref] = jet.phi();
#	Line 401 | Line 669 | HiInclusiveJetAnalyzer::analyze(const Ev
669		jets_.jtpu[jets_.nref] = jet.pileup();
670
671		if(isMC_ && usePat_){
672	<	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
673	<	const reco::GenJet* genjet = (*patjets)[j].genJet();
672	>
673	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
674
675		if(genjet){
676		jets_.refpt[jets_.nref] = genjet->pt();
#	Line 420 | Line 688 | HiInclusiveJetAnalyzer::analyze(const Ev
688		jets_.refdrjt[jets_.nref] = -999.;
689		}
690
691	<	// matched partons
692	<	const reco::GenParticle * parton = (*patjets)[j].genParton();
693	<	if(parton){
694	<	jets_.refparton_pt[jets_.nref] = parton->pt();
695	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
691	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
692	>
693	>	// matched partons
694	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
695	>
696	>	if((*patjets)[j].genParton()){
697	>	jets_.refparton_pt[jets_.nref] = parton.pt();
698	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
699	>
700	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
701	>
702	>	usedStringPts.clear();
703	>
704	>	// uncomment this if you want to know the ugly truth about parton matching -matt
705	>	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
706	>	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
707	>
708	>	jets_.bJetIndex[jets_.bMult] = jets_.nref;
709	>	jets_.bStatus[jets_.bMult] = parton.status();
710	>	jets_.bVx[jets_.bMult] = parton.vx();
711	>	jets_.bVy[jets_.bMult] = parton.vy();
712	>	jets_.bVz[jets_.bMult] = parton.vz();
713	>	jets_.bPt[jets_.bMult] = parton.pt();
714	>	jets_.bEta[jets_.bMult] = parton.eta();
715	>	jets_.bPhi[jets_.bMult] = parton.phi();
716	>	jets_.bPdg[jets_.bMult] = parton.pdgId();
717	>	jets_.bChg[jets_.bMult] = parton.charge();
718	>	jets_.bMult++;
719	>	saveDaughters(parton);
720	>	}
721	>
722	>
723		} else {
724	<	jets_.refparton_pt[jets_.nref] = -999;
725	<	jets_.refparton_flavor[jets_.nref] = -999;
724	>	jets_.refparton_pt[jets_.nref] = -999;
725	>	jets_.refparton_flavor[jets_.nref] = -999;
726		}
727	+
728	+
729		}
730
731		jets_.nref++;
#	Line 439 | Line 736 | HiInclusiveJetAnalyzer::analyze(const Ev
736
737		if(isMC_){
738
739	+	edm::Handle<HepMCProduct> hepMCProduct;
740	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
741	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
742	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
743	+	jets_.beamId1 = beamParticles.first->pdg_id();
744	+	jets_.beamId2 = beamParticles.second->pdg_id();
745	+
746		edm::Handle<GenEventInfoProduct> hEventInfo;
747		iEvent.getByLabel(eventInfoTag_,hEventInfo);
748		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 463 | Line 767 | HiInclusiveJetAnalyzer::analyze(const Ev
767		jets_.genphi[jets_.ngen] = genjet.phi();
768		jets_.geny[jets_.ngen] = genjet.eta();
769
770	+	if(doSubEvent_){
771	+	const GenParticle* gencon = genjet.getGenConstituent(0);
772	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
773	+	}
774
775		// find matching patJet if there is one
776
#	Line 493 | Line 801 | HiInclusiveJetAnalyzer::analyze(const Ev
801
802		}
803		t->Fill();
804	+	memset(&jets_,0,sizeof jets_);
805		}
806
807
#	Line 643 | Line 952 | HiInclusiveJetAnalyzer::getPtRel(const r
952
953		return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
954		}
955	+
956	+	// Recursive function, but this version gets called only the first time
957	+	void
958	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
959	+
960	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
961	+	const reco::Candidate & daughter = *gen.daughter(i);
962	+	double daughterPt = daughter.pt();
963	+	if(daughterPt<1.) continue;
964	+	double daughterEta = daughter.eta();
965	+	if(fabs(daughterEta)>3.) continue;
966	+	int daughterPdgId = daughter.pdgId();
967	+	int daughterStatus = daughter.status();
968	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
969	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
970	+
971	+	// cheesy way of finding strings which were already used
972	+	if(daughter.pdgId()==92){
973	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
974	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
975	+	}
976	+	usedStringPts.push_back(daughter.pt());
977	+	}
978	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
979	+	jets_.bStatus[jets_.bMult] = daughterStatus;
980	+	jets_.bVx[jets_.bMult] = daughter.vx();
981	+	jets_.bVy[jets_.bMult] = daughter.vy();
982	+	jets_.bVz[jets_.bMult] = daughter.vz();
983	+	jets_.bPt[jets_.bMult] = daughterPt;
984	+	jets_.bEta[jets_.bMult] = daughterEta;
985	+	jets_.bPhi[jets_.bMult] = daughter.phi();
986	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
987	+	jets_.bChg[jets_.bMult] = daughter.charge();
988	+	jets_.bMult++;
989	+	saveDaughters(daughter);
990	+	}
991	+	}
992	+
993	+	// This version called for all subsequent calls
994	+	void
995	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
996	+
997	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
998	+	const reco::Candidate & daughter = *gen.daughter(i);
999	+	double daughterPt = daughter.pt();
1000	+	if(daughterPt<1.) continue;
1001	+	double daughterEta = daughter.eta();
1002	+	if(fabs(daughterEta)>3.) continue;
1003	+	int daughterPdgId = daughter.pdgId();
1004	+	int daughterStatus = daughter.status();
1005	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
1006	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
1007	+
1008	+	// cheesy way of finding strings which were already used
1009	+	if(daughter.pdgId()==92){
1010	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
1011	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
1012	+	}
1013	+	usedStringPts.push_back(daughter.pt());
1014	+	}
1015	+
1016	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
1017	+	jets_.bStatus[jets_.bMult] = daughterStatus;
1018	+	jets_.bVx[jets_.bMult] = daughter.vx();
1019	+	jets_.bVy[jets_.bMult] = daughter.vy();
1020	+	jets_.bVz[jets_.bMult] = daughter.vz();
1021	+	jets_.bPt[jets_.bMult] = daughterPt;
1022	+	jets_.bEta[jets_.bMult] = daughterEta;
1023	+	jets_.bPhi[jets_.bMult] = daughter.phi();
1024	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
1025	+	jets_.bChg[jets_.bMult] = daughter.charge();
1026	+	jets_.bMult++;
1027	+	saveDaughters(daughter);
1028	+	}
1029	+	}
1030
1031		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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