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.1 by yilmaz, Tue Sep 20 19:45:06 2011 UTC vs.
Revision 1.18 by yilmaz, Fri May 11 18:43:10 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 43 | Line 43 | using namespace edm;
43		using namespace reco;
44
45		HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
46	–
46
47		jetTag_ = iConfig.getParameter<InputTag>("jetTag");
48	<	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
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	+	doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
57
58	+	rParam = iConfig.getParameter<double>("rParam");
59	+
60		if(isMC_){
61		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
62		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 57 | Line 64 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
64		verbose_ = iConfig.getUntrackedParameter<bool>("verbose",false);
65
66		useCentrality_ = iConfig.getUntrackedParameter<bool>("useCentrality",false);
67	<	useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",true);
67	>	useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",false);
68		useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
69	+	usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
70
71	<	if(!isMC_){
71	>	doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
72	>	doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
73	>	saveBfragments_  = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
74	>
75	>	pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
76	>
77	>	if(doTrigger_){
78		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
79		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
80
#	Line 77 | Line 91 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
91		}
92
93		cout<<" jet collection : "<<jetTag_<<endl;
94	<	if(isMC_)cout<<" genjet collection : "<<genjetTag_<<endl;
94	>	doSubEvent_ = 0;
95
96	+	if(isMC_){
97	+	cout<<" genjet collection : "<<genjetTag_<<endl;
98	+	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
99	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
100	+	}
101
102
103		}
#	Line 102 | Line 121 | HiInclusiveJetAnalyzer::beginJob() {
121		string jetTagTitle = jetTag_.label()+" Jet Analysis Tree";
122		t = fs1->make<TTree>("t",jetTagTitle.c_str());
123
105	–
124		//  TTree* t= new TTree("t","Jet Response Analyzer");
125	<	t->Branch("run",&jets_.run,"run/I");
125	>	//t->Branch("run",&jets_.run,"run/I");
126		t->Branch("evt",&jets_.evt,"evt/I");
127	<	t->Branch("lumi",&jets_.lumi,"lumi/I");
127	>	//t->Branch("lumi",&jets_.lumi,"lumi/I");
128		t->Branch("b",&jets_.b,"b/F");
129	<	t->Branch("vx",&jets_.vx,"vx/F");
130	<	t->Branch("vy",&jets_.vy,"vy/F");
131	<	t->Branch("vz",&jets_.vz,"vz/F");
132	<	t->Branch("hf",&jets_.hf,"hf/F");
129	>	if (useVtx_) {
130	>	t->Branch("vx",&jets_.vx,"vx/F");
131	>	t->Branch("vy",&jets_.vy,"vy/F");
132	>	t->Branch("vz",&jets_.vz,"vz/F");
133	>	}
134	>	if (useCentrality_) {
135	>	t->Branch("hf",&jets_.hf,"hf/F");
136	>	t->Branch("bin",&jets_.bin,"bin/I");
137	>	}
138	>
139		t->Branch("nref",&jets_.nref,"nref/I");
116	–	t->Branch("bin",&jets_.bin,"bin/I");
140		t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
141		t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
142		t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
143		t->Branch("jty",jets_.jty,"jty[nref]/F");
144		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
145	+	t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
146	+
147	+	// jet ID information, jet composition
148	+	t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
149
150	+	t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
151	+	t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
152	+	t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
153	+
154	+	t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
155	+	t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
156	+	t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
157	+
158	+	t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
159	+	t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
160	+	t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
161	+
162	+	t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
163	+	t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
164	+	t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
165	+
166	+	t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
167	+	t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
168	+	t->Branch("eN", jets_.eN,"eN[nref]/I");
169	+
170	+	t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
171	+	t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
172	+	t->Branch("muN", jets_.muN,"muN[nref]/I");
173	+
174	+	t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
175	+	t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
176	+
177	+	// b-jet discriminators
178	+	if (doLifeTimeTagging_) {
179	+
180	+	t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
181	+	t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
182	+	t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
183	+	t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
184	+	t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
185	+	t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
186	+	t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
187	+	t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
188	+
189	+	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
190	+	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
191	+	t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
192	+	t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
193	+	t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
194	+	t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
195	+
196	+	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
197	+	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
198	+
199	+	if (doLifeTimeTaggingExtras_) {
200	+	t->Branch("nIP",&jets_.nIP,"nIP/I");
201	+	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
202	+	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
203	+	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
204	+	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
205	+	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
206	+	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
207	+	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
208	+	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
209	+	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
210	+	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
211	+	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
212	+
213	+	}
214	+
215	+	t->Branch("mue",     jets_.mue,     "mue[nref]/F");
216	+	t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
217	+	t->Branch("mueta",   jets_.mueta,   "mueta[nref]/F");
218	+	t->Branch("muphi",   jets_.muphi,   "muphi[nref]/F");
219	+	t->Branch("mudr",    jets_.mudr,    "mudr[nref]/F");
220	+	t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
221	+	t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
222	+	}
223	+
224	+
225		if(isMC_){
226	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
227	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
228	+
229		t->Branch("pthat",&jets_.pthat,"pthat/F");
230
231		// Only matched gen jets
#	Line 132 | Line 237 | HiInclusiveJetAnalyzer::beginJob() {
237		t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
238		// matched parton
239		t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
240	<	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
240	>	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
241	>	t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
242
243		// For all gen jets, matched or unmatched
244		t->Branch("ngen",&jets_.ngen,"ngen/I");
#	Line 143 | Line 249 | HiInclusiveJetAnalyzer::beginJob() {
249		t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
250		t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
251		t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
252	+
253	+	if(doSubEvent_){
254	+	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
255	+	}
256	+
257	+	if(saveBfragments_  ) {
258	+	t->Branch("bMult",&jets_.bMult,"bMult/I");
259	+	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
260	+	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
261	+	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
262	+	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
263	+	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
264	+	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
265	+	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
266	+	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
267	+	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
268	+	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
269	+	}
270	+
271		}
272	<
272	>	/*
273		if(!isMC_){
274		t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
275		t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
#	Line 156 | Line 281 | HiInclusiveJetAnalyzer::beginJob() {
281		t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
282
283		}
284	<
284	>	*/
285		TH1D::SetDefaultSumw2();
286
287
#	Line 177 | Line 302 | HiInclusiveJetAnalyzer::analyze(const Ev
302
303		LogDebug("HiInclusiveJetAnalyzer")<<"START event: "<<event<<" in run "<<run<<endl;
304
180	–
305		int bin = -1;
306		double hf = 0.;
307		double b = 999.;
308
309		if(useCentrality_){
186	–	//if(!isMC_){
310		if(!centrality_) centrality_ = new CentralityProvider(iSetup);
311		centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
312		//double c = centrality_->centralityValue();
#	Line 193 | Line 316 | HiInclusiveJetAnalyzer::analyze(const Ev
316
317		bin = centrality_->getBin();
318		b = centrality_->bMean();
196	–	//}
197	–	/*
198	–	else{
199	–
200	–	TFile * centFile = new TFile("/net/hidsk0001/d00/scratch/mnguyen/CMSSW_3_9_1_patch1/src/macros/Hydjet_CentTable.root");
201	–
202	–	edm::Handle<reco::Centrality> cent;
203	–	iEvent.getByLabel(edm::InputTag("hiCentrality"),cent);
204	–	//cout<<" grabbed centrality "<<endl;
205	–	CentralityBins::RunMap cmap = getCentralityFromFile(centFile, "makeCentralityTableTFile", "HFhitsHydjet_2760GeV", 1, 1);
206	–
207	–	// Still getting cent from hits.  When tower tables become available, we need to switch
208	–	hf = cent->EtHFhitSum();
209	–	//cout<<" hf "<<hf<<endl;
210	–	bin = cmap[run]->getBin(hf);
211	–	b = cmap[run]->bMeanOfBin(bin);
212	–	}
213	–	*/
319		}
320
321
322
323
219	–	//double jetPtMin = 35.;
220	–
324
325		// loop the events
326
#	Line 236 | Line 339 | HiInclusiveJetAnalyzer::analyze(const Ev
339		}
340		}
341
342	+	edm::Handle<pat::JetCollection> patjets;
343	+	if(usePat_)iEvent.getByLabel(jetTag_, patjets);
344
345	<
346	<	edm::Handle<pat::JetCollection> jets;
345	>	edm::Handle<pat::JetCollection> matchedjets;
346	>	iEvent.getByLabel(matchTag_, matchedjets);
347	>
348	>	edm::Handle<reco::JetView> jets;
349		iEvent.getByLabel(jetTag_, jets);
350
351	<
351	>	edm::Handle<reco::PFCandidateCollection> pfCandidates;
352	>	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
353	>
354	>	edm::Handle<reco::TrackCollection> tracks;
355	>	iEvent.getByLabel(trackTag_,tracks);
356	>
357		// FILL JRA TREE
358
359		jets_.b = b;
360		jets_.nref = 0;
361
362	<	if(!isMC_){
362	>	if(doTrigger_){
363		fillL1Bits(iEvent);
364		fillHLTBits(iEvent);
365		}
366
367		for(unsigned int j = 0 ; j < jets->size(); ++j){
368	<	const pat::Jet& jet = (*jets)[j];
368	>	const reco::Jet& jet = (*jets)[j];
369	>
370	>
371	>	if (useJEC_ && usePat_){
372	>	jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
373	>	}
374
375	<	//cout<<" jet pt "<<jet.pt()<<endl;
376	<	//if(jet.pt() < jetPtMin) continue;
377	<	if (useJEC_) jets_.rawpt[jets_.nref]=jet.correctedJet("Uncorrected").pt();
375	>	if(doLifeTimeTagging_){
376	>
377	>	if(jetTag_.label()=="icPu5patJets"){
378	>	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
379	>	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
380	>	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
381	>	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
382	>	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
383	>	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
384	>	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
385	>	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
386	>	}
387	>	else if(jetTag_.label()=="akPu3PFpatJets"){
388	>	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
389	>	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
390	>	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
391	>	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
392	>	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
393	>	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
394	>	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
395	>	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
396	>	}
397	>	else{
398	>	cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
399	>	}
400	>
401	>	const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
402	>
403	>	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
404	>
405	>	if (tagInfoSV.nVertices()>0) {
406	>	jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
407	>	// this is the 3d flight distance, for 2-D use (0,true)
408	>	Measurement1D m1D = tagInfoSV.flightDistance(0);
409	>	jets_.svtxdl[jets_.nref]    = m1D.value();
410	>	jets_.svtxdls[jets_.nref]   = m1D.significance();
411	>
412	>	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
413	>	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
414	>	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
415	>	jets_.svtxm[jets_.nref]    = svtx.p4().mass();
416	>	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
417	>	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
418	>	}
419	>
420	>	const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
421	>
422	>	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
423	>	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
424	>
425	>	if (doLifeTimeTaggingExtras_) {
426	>
427	>	TrackRefVector selTracks=tagInfoIP.selectedTracks();
428	>
429	>	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
430	>
431	>	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
432	>	{
433	>	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
434	>	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
435	>	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
436	>	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
437	>	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
438	>	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
439	>	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
440	>	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
441	>	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
442	>	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
443	>	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
444	>	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
445	>	}
446	>
447	>	jets_.nIP += jets_.nselIPtrk[jets_.nref];
448	>
449	>	}
450	>
451	>	const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
452	>	int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
453	>	if(pfMuonIndex >=0){
454	>	const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
455	>	jets_.mupt[jets_.nref]    =  muon.pt();
456	>	jets_.mueta[jets_.nref]   =  muon.eta();
457	>	jets_.muphi[jets_.nref]   =  muon.phi();
458	>	jets_.mue[jets_.nref]     =  muon.energy();
459	>	jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
460	>	jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
461	>	jets_.muchg[jets_.nref]   =  muon.charge();
462	>	}else{
463	>	jets_.mupt[jets_.nref]    =  0.0;
464	>	jets_.mueta[jets_.nref]   =  0.0;
465	>	jets_.muphi[jets_.nref]   =  0.0;
466	>	jets_.mue[jets_.nref]     =  0.0;
467	>	jets_.mudr[jets_.nref]    =  9.9;
468	>	jets_.muptrel[jets_.nref] =  0.0;
469	>	jets_.muchg[jets_.nref]   = 0;
470	>	}
471	>
472	>	}
473	>
474	>
475	>
476	>	// Jet ID variables
477	>
478	>	jets_.muMax[jets_.nref] = 0;
479	>	jets_.muSum[jets_.nref] = 0;
480	>	jets_.muN[jets_.nref] = 0;
481	>
482	>	jets_.eMax[jets_.nref] = 0;
483	>	jets_.eSum[jets_.nref] = 0;
484	>	jets_.eN[jets_.nref] = 0;
485	>
486	>	jets_.neutralMax[jets_.nref] = 0;
487	>	jets_.neutralSum[jets_.nref] = 0;
488	>	jets_.neutralN[jets_.nref] = 0;
489	>
490	>	jets_.photonMax[jets_.nref] = 0;
491	>	jets_.photonSum[jets_.nref] = 0;
492	>	jets_.photonN[jets_.nref] = 0;
493	>
494	>	jets_.chargedMax[jets_.nref] = 0;
495	>	jets_.chargedSum[jets_.nref] = 0;
496	>	jets_.chargedN[jets_.nref] = 0;
497	>
498	>	jets_.trackMax[jets_.nref] = 0;
499	>	jets_.trackSum[jets_.nref] = 0;
500	>	jets_.trackN[jets_.nref] = 0;
501	>
502	>
503	>	for(unsigned int icand = 0; icand < tracks->size(); ++icand){
504	>	const reco::Track& track = (*tracks)[icand];
505	>	if(useQuality_ && !(track.quality(reco::TrackBase::qualityByName(trackQuality_)))) continue;
506	>
507	>	double dr = deltaR(jet,track);
508	>	if(dr < rParam){
509	>	double ptcand = track.pt();
510	>	jets_.trackSum[jets_.nref] += ptcand;
511	>	jets_.trackN[jets_.nref] += 1;
512	>	if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
513	>
514	>	}
515	>	}
516	>
517	>	for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
518	>	const reco::PFCandidate& track = (*pfCandidates)[icand];
519	>	double dr = deltaR(jet,track);
520	>	if(dr < rParam){
521	>	double ptcand = track.pt();
522	>	int pfid = track.particleId();
523	>
524	>	switch(pfid){
525	>
526	>	case 1:
527	>	jets_.chargedSum[jets_.nref] += ptcand;
528	>	jets_.chargedN[jets_.nref] += 1;
529	>	if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
530	>	break;
531	>
532	>	case 2:
533	>	jets_.eSum[jets_.nref] += ptcand;
534	>	jets_.eN[jets_.nref] += 1;
535	>	if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
536	>	break;
537	>
538	>	case 3:
539	>	jets_.muSum[jets_.nref] += ptcand;
540	>	jets_.muN[jets_.nref] += 1;
541	>	if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
542	>	break;
543	>
544	>	case 4:
545	>	jets_.photonSum[jets_.nref] += ptcand;
546	>	jets_.photonN[jets_.nref] += 1;
547	>	if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
548	>	break;
549	>
550	>	case 5:
551	>	jets_.neutralSum[jets_.nref] += ptcand;
552	>	jets_.neutralN[jets_.nref] += 1;
553	>	if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
554	>	break;
555	>
556	>	}
557	>	}
558	>	}
559	>
560	>	double drMin = 100;
561	>	for(unsigned int j = 0 ; j < matchedjets->size(); ++j){
562	>	const reco::Jet& mjet = (*matchedjets)[j];
563	>
564	>	double dr = deltaR(jet,mjet);
565	>	if(dr < drMin){
566	>	jets_.matchedPt[jets_.nref] = mjet.pt();
567	>	jets_.matchedR[jets_.nref] = dr;
568	>
569	>	}
570	>	}
571	>
572	>
573	>
574	>
575	>	//     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
576	>
577	>	/////////////////////////////////////////////////////////////////
578	>	// Jet core pt^2 discriminant for fake jets
579	>	// Edited by Yue Shi Lai <ylai@mit.edu>
580	>
581	>	// Initial value is 0
582	>	jets_.discr_fr01[jets_.nref] = 0;
583	>	// Start with no directional adaption, i.e. the fake rejection
584	>	// axis is the jet axis
585	>	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
586	>	float azimuth_adapt = jets_.jtphi[jets_.nref];
587	>
588	>	// Unadapted discriminant with adaption search
589	>	for (size_t iteration = 0; iteration < 2; iteration++) {
590	>	float pseudorapidity_adapt_new = pseudorapidity_adapt;
591	>	float azimuth_adapt_new = azimuth_adapt;
592	>	float max_weighted_perp = 0;
593	>	float perp_square_sum = 0;
594	>
595	>	for (size_t index_pf_candidate = 0;
596	>	index_pf_candidate < pfCandidates->size();
597	>	index_pf_candidate++) {
598	>	const reco::PFCandidate &p =
599	>	(*pfCandidates)[index_pf_candidate];
600	>
601	>	switch (p.particleId()) {
602	>	case 1: // Charged hadron
603	>	case 3: // Muon
604	>	case 4: // Photon
605	>	{
606	>	const float dpseudorapidity =
607	>	p.eta() - pseudorapidity_adapt;
608	>	const float dazimuth =
609	>	reco::deltaPhi(p.phi(), azimuth_adapt);
610	>	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
611	>	// = 50
612	>	const float angular_weight =
613	>	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
614	>	dazimuth * dazimuth));
615	>	const float weighted_perp =
616	>	angular_weight * p.pt() * p.pt();
617	>	const float weighted_perp_square =
618	>	weighted_perp * p.pt();
619	>
620	>	perp_square_sum += weighted_perp_square;
621	>	if (weighted_perp >= max_weighted_perp) {
622	>	pseudorapidity_adapt_new = p.eta();
623	>	azimuth_adapt_new = p.phi();
624	>	max_weighted_perp = weighted_perp;
625	>	}
626	>	}
627	>	}
628	>	}
629	>	// Update the fake rejection value
630	>	jets_.discr_fr01[jets_.nref] = std::max(
631	>	jets_.discr_fr01[jets_.nref], perp_square_sum);
632	>	// Update the directional adaption
633	>	pseudorapidity_adapt = pseudorapidity_adapt_new;
634	>	azimuth_adapt = azimuth_adapt_new;
635	>	}
636	>
637	>
638		jets_.jtpt[jets_.nref] = jet.pt();
639		jets_.jteta[jets_.nref] = jet.eta();
640		jets_.jtphi[jets_.nref] = jet.phi();
641		jets_.jty[jets_.nref] = jet.eta();
642	<
642	>	jets_.jtpu[jets_.nref] = jet.pileup();
643
644	<	if(jet.genJet()){
268	<	jets_.refpt[jets_.nref] = jet.genJet()->pt();
269	<	jets_.refeta[jets_.nref] = jet.genJet()->eta();
270	<	jets_.refphi[jets_.nref] = jet.genJet()->phi();
271	<	jets_.refy[jets_.nref] = jet.genJet()->eta();
272	<	jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(),jet.genJet()->phi());
273	<	jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),jet.genJet()->eta(),jet.genJet()->phi());
274	<	}
275	<	else{
276	<	jets_.refpt[jets_.nref] = -999.;
277	<	jets_.refeta[jets_.nref] = -999.;
278	<	jets_.refphi[jets_.nref] = -999.;
279	<	jets_.refy[jets_.nref] = -999.;
280	<	jets_.refdphijt[jets_.nref] = -999.;
281	<	jets_.refdrjt[jets_.nref] = -999.;
282	<	}
644	>	if(isMC_ && usePat_){
645
646	<	// matched partons
647	<	if (jet.genParton()) {
648	<	jets_.refparton_pt[jets_.nref] = jet.genParton()->pt();
649	<	jets_.refparton_flavor[jets_.nref] = jet.genParton()->pdgId();
650	<	} else {
651	<	jets_.refparton_pt[jets_.nref] = -999;
652	<	jets_.refparton_flavor[jets_.nref] = -999;
653	<	}
646	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
647	>
648	>	if(genjet){
649	>	jets_.refpt[jets_.nref] = genjet->pt();
650	>	jets_.refeta[jets_.nref] = genjet->eta();
651	>	jets_.refphi[jets_.nref] = genjet->phi();
652	>	jets_.refy[jets_.nref] = genjet->eta();
653	>	jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
654	>	jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
655	>	}else{
656	>	jets_.refpt[jets_.nref] = -999.;
657	>	jets_.refeta[jets_.nref] = -999.;
658	>	jets_.refphi[jets_.nref] = -999.;
659	>	jets_.refy[jets_.nref] = -999.;
660	>	jets_.refdphijt[jets_.nref] = -999.;
661	>	jets_.refdrjt[jets_.nref] = -999.;
662	>	}
663	>
664	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
665	>
666	>	// matched partons
667	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
668	>
669	>	if((*patjets)[j].genParton()){
670	>	jets_.refparton_pt[jets_.nref] = parton.pt();
671	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
672	>
673	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
674
675	+	usedStringPts.clear();
676	+
677	+	// uncomment this if you want to know the ugly truth about parton matching -matt
678	+	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
679	+	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
680	+
681	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
682	+	jets_.bStatus[jets_.bMult] = parton.status();
683	+	jets_.bVx[jets_.bMult] = parton.vx();
684	+	jets_.bVy[jets_.bMult] = parton.vy();
685	+	jets_.bVz[jets_.bMult] = parton.vz();
686	+	jets_.bPt[jets_.bMult] = parton.pt();
687	+	jets_.bEta[jets_.bMult] = parton.eta();
688	+	jets_.bPhi[jets_.bMult] = parton.phi();
689	+	jets_.bPdg[jets_.bMult] = parton.pdgId();
690	+	jets_.bChg[jets_.bMult] = parton.charge();
691	+	jets_.bMult++;
692	+	saveDaughters(parton);
693	+	}
694	+
695	+
696	+	} else {
697	+	jets_.refparton_pt[jets_.nref] = -999;
698	+	jets_.refparton_flavor[jets_.nref] = -999;
699	+	}
700	+
701	+
702	+	}
703
704		jets_.nref++;
705
#	Line 299 | Line 709 | HiInclusiveJetAnalyzer::analyze(const Ev
709
710		if(isMC_){
711
712	+	edm::Handle<HepMCProduct> hepMCProduct;
713	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
714	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
715	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
716	+	jets_.beamId1 = beamParticles.first->pdg_id();
717	+	jets_.beamId2 = beamParticles.second->pdg_id();
718	+
719		edm::Handle<GenEventInfoProduct> hEventInfo;
720		iEvent.getByLabel(eventInfoTag_,hEventInfo);
721		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 316 | Line 733 | HiInclusiveJetAnalyzer::analyze(const Ev
733		float genjet_pt = genjet.pt();
734
735		// threshold to reduce size of output in minbias PbPb
736	<	if(genjet_pt>20.){
736	>	if(genjet_pt>genPtMin_){
737
738		jets_.genpt[jets_.ngen] = genjet_pt;
739		jets_.geneta[jets_.ngen] = genjet.eta();
740		jets_.genphi[jets_.ngen] = genjet.phi();
741		jets_.geny[jets_.ngen] = genjet.eta();
742
743	+	if(doSubEvent_){
744	+	const GenParticle* gencon = genjet.getGenConstituent(0);
745	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
746	+	}
747
748		// find matching patJet if there is one
749
750		jets_.gendrjt[jets_.ngen] = -1.0;
751		jets_.genmatchindex[jets_.ngen] = -1;
752
332	–
753		for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
754		const pat::Jet& jet = (*jets)[ijet];
755	<
756	<	if(jet.genJet()){
757	<	if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
758	<	fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
759	<	(fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
755	>	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
756	>	if(matchedGenJet){
757	>	// poor man's matching, someone fix please
758	>	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
759	>	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
760	>	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
761
762		jets_.genmatchindex[jets_.ngen] = (int)ijet;
763		jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
#	Line 353 | Line 774 | HiInclusiveJetAnalyzer::analyze(const Ev
774
775		}
776		t->Fill();
777	+	memset(&jets_,0,sizeof jets_);
778		}
779
780
#	Line 452 | Line 874 | inline bool HiInclusiveJetAnalyzer::getP
874		}
875
876
877	+	int
878	+	HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
879	+	{
880	+
881	+	int pfMuonIndex = -1;
882	+	float ptMax = 0.;
883	+
884
885	+	for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
886	+	const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
887	+
888	+	int id = pfCandidate.particleId();
889	+	if(abs(id) != 3) continue;
890	+
891	+	if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
892	+
893	+	double pt =  pfCandidate.pt();
894	+	if(pt>ptMax){
895	+	ptMax = pt;
896	+	pfMuonIndex = (int) icand;
897	+	}
898	+	}
899	+
900	+	return pfMuonIndex;
901	+
902	+	}
903	+
904	+
905	+	double
906	+	HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
907	+	{
908	+
909	+	float lj_x = jet.p4().px();
910	+	float lj_y = jet.p4().py();
911	+	float lj_z = jet.p4().pz();
912	+
913	+	// absolute values squared
914	+	float lj2  = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
915	+	float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
916	+
917	+	// projection vec(mu) to lepjet axis
918	+	float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
919	+
920	+	// absolute value squared and normalized
921	+	float pLrel2 = lepXlj*lepXlj/lj2;
922	+
923	+	// lep2 = pTrel2 + pLrel2
924	+	float pTrel2 = lep2-pLrel2;
925	+
926	+	return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
927	+	}
928	+
929	+	// Recursive function, but this version gets called only the first time
930	+	void
931	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
932	+
933	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
934	+	const reco::Candidate & daughter = *gen.daughter(i);
935	+	double daughterPt = daughter.pt();
936	+	if(daughterPt<1.) continue;
937	+	double daughterEta = daughter.eta();
938	+	if(fabs(daughterEta)>3.) continue;
939	+	int daughterPdgId = daughter.pdgId();
940	+	int daughterStatus = daughter.status();
941	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
942	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
943	+
944	+	// cheesy way of finding strings which were already used
945	+	if(daughter.pdgId()==92){
946	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
947	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
948	+	}
949	+	usedStringPts.push_back(daughter.pt());
950	+	}
951	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
952	+	jets_.bStatus[jets_.bMult] = daughterStatus;
953	+	jets_.bVx[jets_.bMult] = daughter.vx();
954	+	jets_.bVy[jets_.bMult] = daughter.vy();
955	+	jets_.bVz[jets_.bMult] = daughter.vz();
956	+	jets_.bPt[jets_.bMult] = daughterPt;
957	+	jets_.bEta[jets_.bMult] = daughterEta;
958	+	jets_.bPhi[jets_.bMult] = daughter.phi();
959	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
960	+	jets_.bChg[jets_.bMult] = daughter.charge();
961	+	jets_.bMult++;
962	+	saveDaughters(daughter);
963	+	}
964	+	}
965	+
966	+	// This version called for all subsequent calls
967	+	void
968	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
969	+
970	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
971	+	const reco::Candidate & daughter = *gen.daughter(i);
972	+	double daughterPt = daughter.pt();
973	+	if(daughterPt<1.) continue;
974	+	double daughterEta = daughter.eta();
975	+	if(fabs(daughterEta)>3.) continue;
976	+	int daughterPdgId = daughter.pdgId();
977	+	int daughterStatus = daughter.status();
978	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
979	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
980	+
981	+	// cheesy way of finding strings which were already used
982	+	if(daughter.pdgId()==92){
983	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
984	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
985	+	}
986	+	usedStringPts.push_back(daughter.pt());
987	+	}
988	+
989	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
990	+	jets_.bStatus[jets_.bMult] = daughterStatus;
991	+	jets_.bVx[jets_.bMult] = daughter.vx();
992	+	jets_.bVy[jets_.bMult] = daughter.vy();
993	+	jets_.bVz[jets_.bMult] = daughter.vz();
994	+	jets_.bPt[jets_.bMult] = daughterPt;
995	+	jets_.bEta[jets_.bMult] = daughterEta;
996	+	jets_.bPhi[jets_.bMult] = daughter.phi();
997	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
998	+	jets_.bChg[jets_.bMult] = daughter.charge();
999	+	jets_.bMult++;
1000	+	saveDaughters(daughter);
1001	+	}
1002	+	}
1003
1004		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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