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.2 by yilmaz, Tue Sep 20 20:06:06 2011 UTC vs.
Revision 1.16 by ylai, Fri May 11 14:39:32 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	>	doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
53	>
54		if(isMC_){
55		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
56		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 57 | Line 58 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
58		verbose_ = iConfig.getUntrackedParameter<bool>("verbose",false);
59
60		useCentrality_ = iConfig.getUntrackedParameter<bool>("useCentrality",false);
61	<	useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",true);
61	>	useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",false);
62		useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
63		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
64
65	<	if(!isMC_){
65	>	doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
66	>	doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
67	>	saveBfragments_  = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
68	>
69	>	pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
70	>
71	>	if(doTrigger_){
72		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
73		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
74
#	Line 78 | Line 85 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
85		}
86
87		cout<<" jet collection : "<<jetTag_<<endl;
88	<	if(isMC_)cout<<" genjet collection : "<<genjetTag_<<endl;
88	>	doSubEvent_ = 0;
89
90	+	if(isMC_){
91	+	cout<<" genjet collection : "<<genjetTag_<<endl;
92	+	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
93	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
94	+	}
95
96
97		}
#	Line 104 | Line 116 | HiInclusiveJetAnalyzer::beginJob() {
116		t = fs1->make<TTree>("t",jetTagTitle.c_str());
117
118		//  TTree* t= new TTree("t","Jet Response Analyzer");
119	<	t->Branch("run",&jets_.run,"run/I");
119	>	//t->Branch("run",&jets_.run,"run/I");
120		t->Branch("evt",&jets_.evt,"evt/I");
121	<	t->Branch("lumi",&jets_.lumi,"lumi/I");
121	>	//t->Branch("lumi",&jets_.lumi,"lumi/I");
122		t->Branch("b",&jets_.b,"b/F");
123	<	t->Branch("vx",&jets_.vx,"vx/F");
124	<	t->Branch("vy",&jets_.vy,"vy/F");
125	<	t->Branch("vz",&jets_.vz,"vz/F");
126	<	t->Branch("hf",&jets_.hf,"hf/F");
123	>	if (useVtx_) {
124	>	t->Branch("vx",&jets_.vx,"vx/F");
125	>	t->Branch("vy",&jets_.vy,"vy/F");
126	>	t->Branch("vz",&jets_.vz,"vz/F");
127	>	}
128	>	if (useCentrality_) {
129	>	t->Branch("hf",&jets_.hf,"hf/F");
130	>	t->Branch("bin",&jets_.bin,"bin/I");
131	>	}
132	>
133		t->Branch("nref",&jets_.nref,"nref/I");
116	–	t->Branch("bin",&jets_.bin,"bin/I");
134		t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
135		t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
136		t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
#	Line 121 | Line 138 | HiInclusiveJetAnalyzer::beginJob() {
138		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
139		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
140
141	+	// b-jet discriminators
142	+	if (doLifeTimeTagging_) {
143	+
144	+	t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
145	+	t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
146	+	t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
147	+	t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
148	+	t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
149	+	t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
150	+	t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
151	+	t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
152	+
153	+	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
154	+	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
155	+	t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
156	+	t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
157	+	t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
158	+	t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
159	+
160	+	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
161	+	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
162	+
163	+	if (doLifeTimeTaggingExtras_) {
164	+	t->Branch("nIP",&jets_.nIP,"nIP/I");
165	+	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
166	+	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
167	+	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
168	+	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
169	+	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
170	+	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
171	+	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
172	+	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
173	+	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
174	+	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
175	+	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
176	+
177	+	}
178	+
179	+	t->Branch("mue",     jets_.mue,     "mue[nref]/F");
180	+	t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
181	+	t->Branch("mueta",   jets_.mueta,   "mueta[nref]/F");
182	+	t->Branch("muphi",   jets_.muphi,   "muphi[nref]/F");
183	+	t->Branch("mudr",    jets_.mudr,    "mudr[nref]/F");
184	+	t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
185	+	t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
186	+	}
187	+
188	+	t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
189	+
190		if(isMC_){
191	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
192	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
193	+
194		t->Branch("pthat",&jets_.pthat,"pthat/F");
195
196		// Only matched gen jets
#	Line 133 | Line 202 | HiInclusiveJetAnalyzer::beginJob() {
202		t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
203		// matched parton
204		t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
205	<	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
205	>	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
206	>	t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
207
208		// For all gen jets, matched or unmatched
209		t->Branch("ngen",&jets_.ngen,"ngen/I");
#	Line 144 | Line 214 | HiInclusiveJetAnalyzer::beginJob() {
214		t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
215		t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
216		t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
217	+
218	+	if(doSubEvent_){
219	+	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
220	+	}
221	+
222	+	if(saveBfragments_  ) {
223	+	t->Branch("bMult",&jets_.bMult,"bMult/I");
224	+	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
225	+	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
226	+	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
227	+	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
228	+	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
229	+	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
230	+	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
231	+	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
232	+	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
233	+	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
234	+	}
235	+
236		}
237	<
237	>	/*
238		if(!isMC_){
239		t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
240		t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
#	Line 157 | Line 246 | HiInclusiveJetAnalyzer::beginJob() {
246		t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
247
248		}
249	<
249	>	*/
250		TH1D::SetDefaultSumw2();
251
252
#	Line 178 | Line 267 | HiInclusiveJetAnalyzer::analyze(const Ev
267
268		LogDebug("HiInclusiveJetAnalyzer")<<"START event: "<<event<<" in run "<<run<<endl;
269
181	–
270		int bin = -1;
271		double hf = 0.;
272		double b = 999.;
273
274		if(useCentrality_){
187	–	//if(!isMC_){
275		if(!centrality_) centrality_ = new CentralityProvider(iSetup);
276		centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
277		//double c = centrality_->centralityValue();
#	Line 194 | Line 281 | HiInclusiveJetAnalyzer::analyze(const Ev
281
282		bin = centrality_->getBin();
283		b = centrality_->bMean();
197	–	//}
198	–	/*
199	–	else{
200	–
201	–	TFile * centFile = new TFile("/net/hidsk0001/d00/scratch/mnguyen/CMSSW_3_9_1_patch1/src/macros/Hydjet_CentTable.root");
202	–
203	–	edm::Handle<reco::Centrality> cent;
204	–	iEvent.getByLabel(edm::InputTag("hiCentrality"),cent);
205	–	//cout<<" grabbed centrality "<<endl;
206	–	CentralityBins::RunMap cmap = getCentralityFromFile(centFile, "makeCentralityTableTFile", "HFhitsHydjet_2760GeV", 1, 1);
207	–
208	–	// Still getting cent from hits.  When tower tables become available, we need to switch
209	–	hf = cent->EtHFhitSum();
210	–	//cout<<" hf "<<hf<<endl;
211	–	bin = cmap[run]->getBin(hf);
212	–	b = cmap[run]->bMeanOfBin(bin);
213	–	}
214	–	*/
284		}
285
286
287
288
220	–	//double jetPtMin = 35.;
221	–
289
290		// loop the events
291
#	Line 243 | Line 310 | HiInclusiveJetAnalyzer::analyze(const Ev
310		edm::Handle<reco::JetView> jets;
311		iEvent.getByLabel(jetTag_, jets);
312
313	+	edm::Handle<reco::PFCandidateCollection> pfCandidates;
314	+	if(doLifeTimeTagging_){
315	+	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
316	+	}
317		// FILL JRA TREE
318
319		jets_.b = b;
320		jets_.nref = 0;
321
322	<	if(!isMC_){
322	>	if(doTrigger_){
323		fillL1Bits(iEvent);
324		fillHLTBits(iEvent);
325		}
#	Line 256 | Line 327 | HiInclusiveJetAnalyzer::analyze(const Ev
327		for(unsigned int j = 0 ; j < jets->size(); ++j){
328		const reco::Jet& jet = (*jets)[j];
329
330	<	//cout<<" jet pt "<<jet.pt()<<endl;
260	<	//if(jet.pt() < jetPtMin) continue;
330	>
331		if (useJEC_ && usePat_){
332		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
333		}
334	+
335	+	if(doLifeTimeTagging_){
336	+
337	+	if(jetTag_.label()=="icPu5patJets"){
338	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
339	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
340	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
341	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
342	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
343	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
344	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
345	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
346	+	}
347	+	else if(jetTag_.label()=="akPu3PFpatJets"){
348	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
349	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
350	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
351	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
352	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
353	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
354	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
355	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
356	+	}
357	+	else{
358	+	cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
359	+	}
360	+
361	+	const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
362	+
363	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
364	+
365	+	if (tagInfoSV.nVertices()>0) {
366	+	jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
367	+	// this is the 3d flight distance, for 2-D use (0,true)
368	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
369	+	jets_.svtxdl[jets_.nref]    = m1D.value();
370	+	jets_.svtxdls[jets_.nref]   = m1D.significance();
371	+
372	+	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
373	+	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
374	+	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
375	+	jets_.svtxm[jets_.nref]    = svtx.p4().mass();
376	+	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
377	+	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
378	+	}
379	+
380	+	const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
381	+
382	+	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
383	+	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
384	+
385	+	if (doLifeTimeTaggingExtras_) {
386	+
387	+	TrackRefVector selTracks=tagInfoIP.selectedTracks();
388	+
389	+	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
390	+
391	+	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
392	+	{
393	+	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
394	+	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
395	+	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
396	+	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
397	+	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
398	+	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
399	+	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
400	+	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
401	+	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
402	+	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
403	+	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
404	+	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
405	+	}
406	+
407	+	jets_.nIP += jets_.nselIPtrk[jets_.nref];
408	+
409	+	}
410	+
411	+	const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
412	+	int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
413	+	if(pfMuonIndex >=0){
414	+	const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
415	+	jets_.mupt[jets_.nref]    =  muon.pt();
416	+	jets_.mueta[jets_.nref]   =  muon.eta();
417	+	jets_.muphi[jets_.nref]   =  muon.phi();
418	+	jets_.mue[jets_.nref]     =  muon.energy();
419	+	jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
420	+	jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
421	+	jets_.muchg[jets_.nref]   =  muon.charge();
422	+	}
423	+	<<<<<<< HiInclusiveJetAnalyzer.cc
424	+	else{
425	+	jets_.mupt[jets_.nref]    =  0.0;
426	+	jets_.mueta[jets_.nref]   =  0.0;
427	+	jets_.muphi[jets_.nref]   =  0.0;
428	+	jets_.mue[jets_.nref]     =  0.0;
429	+	jets_.mudr[jets_.nref]    =  9.9;
430	+	jets_.muptrel[jets_.nref] =  0.0;
431	+	jets_.muchg[jets_.nref]   = 0;
432	+	}
433	+
434	+	=======
435	+	>>>>>>> 1.15
436	+	}
437	+
438	+	/////////////////////////////////////////////////////////////////
439	+	// Jet core pt^2 discriminant for fake jets
440	+	// Edited by Yue Shi Lai <ylai@mit.edu>
441	+
442	+	// Initial value is 0
443	+	jets_.discr_fr01[jets_.nref] = 0;
444	+	// Start with no directional adaption, i.e. the fake rejection
445	+	// axis is the jet axis
446	+	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
447	+	float azimuth_adapt = jets_.jtphi[jets_.nref];
448	+
449	+	// Unadapted discriminant with adaption search
450	+	for (size_t iteration = 0; iteration < 2; iteration++) {
451	+	float pseudorapidity_adapt_new = pseudorapidity_adapt;
452	+	float azimuth_adapt_new = azimuth_adapt;
453	+	float max_weighted_perp = 0;
454	+	float perp_square_sum = 0;
455	+
456	+	for (size_t index_pf_candidate = 0;
457	+	index_pf_candidate < pfCandidates->size();
458	+	index_pf_candidate++) {
459	+	const reco::PFCandidate &p =
460	+	(*pfCandidates)[index_pf_candidate];
461	+
462	+	switch (p.particleId()) {
463	+	case 1: // Charged hadron
464	+	case 3: // Muon
465	+	case 4: // Photon
466	+	{
467	+	const float dpseudorapidity =
468	+	p.eta() - pseudorapidity_adapt;
469	+	const float dazimuth =
470	+	reco::deltaPhi(p.phi(), azimuth_adapt);
471	+	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
472	+	// = 50
473	+	const float angular_weight =
474	+	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
475	+	dazimuth * dazimuth));
476	+	const float weighted_perp =
477	+	angular_weight * p.pt() * p.pt();
478	+	const float weighted_perp_square =
479	+	weighted_perp * p.pt();
480	+
481	+	perp_square_sum += weighted_perp_square;
482	+	if (weighted_perp >= max_weighted_perp) {
483	+	pseudorapidity_adapt_new = p.eta();
484	+	azimuth_adapt_new = p.phi();
485	+	max_weighted_perp = weighted_perp;
486	+	}
487	+	}
488	+	}
489	+	}
490	+	// Update the fake rejection value
491	+	jets_.discr_fr01[jets_.nref] = std::max(
492	+	jets_.discr_fr01[jets_.nref], perp_square_sum);
493	+	// Update the directional adaption
494	+	pseudorapidity_adapt = pseudorapidity_adapt_new;
495	+	azimuth_adapt = azimuth_adapt_new;
496	+	}
497	+
498	+
499		jets_.jtpt[jets_.nref] = jet.pt();
500		jets_.jteta[jets_.nref] = jet.eta();
501		jets_.jtphi[jets_.nref] = jet.phi();
#	Line 268 | Line 503 | HiInclusiveJetAnalyzer::analyze(const Ev
503		jets_.jtpu[jets_.nref] = jet.pileup();
504
505		if(isMC_ && usePat_){
506	<	const reco::GenJet* genjet = (*patjets)[j].genJet();
506	>
507	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
508	>
509		if(genjet){
510		jets_.refpt[jets_.nref] = genjet->pt();
511		jets_.refeta[jets_.nref] = genjet->eta();
#	Line 285 | Line 522 | HiInclusiveJetAnalyzer::analyze(const Ev
522		jets_.refdrjt[jets_.nref] = -999.;
523		}
524
525	<	// matched partons
526	<	const reco::GenParticle * parton = (*patjets)[j].genParton();
527	<	if(parton){
528	<	jets_.refparton_pt[jets_.nref] = parton->pt();
529	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
530	<	} else {
531	<	jets_.refparton_pt[jets_.nref] = -999;
532	<	jets_.refparton_flavor[jets_.nref] = -999;
533	<	}
525	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
526	>
527	>	// matched partons
528	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
529	>
530	>	if((*patjets)[j].genParton()){
531	>	jets_.refparton_pt[jets_.nref] = parton.pt();
532	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
533	>
534	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
535	>
536	>	usedStringPts.clear();
537	>
538	>	// uncomment this if you want to know the ugly truth about parton matching -matt
539	>	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
540	>	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
541	>
542	>	jets_.bJetIndex[jets_.bMult] = jets_.nref;
543	>	jets_.bStatus[jets_.bMult] = parton.status();
544	>	jets_.bVx[jets_.bMult] = parton.vx();
545	>	jets_.bVy[jets_.bMult] = parton.vy();
546	>	jets_.bVz[jets_.bMult] = parton.vz();
547	>	jets_.bPt[jets_.bMult] = parton.pt();
548	>	jets_.bEta[jets_.bMult] = parton.eta();
549	>	jets_.bPhi[jets_.bMult] = parton.phi();
550	>	jets_.bPdg[jets_.bMult] = parton.pdgId();
551	>	jets_.bChg[jets_.bMult] = parton.charge();
552	>	jets_.bMult++;
553	>	saveDaughters(parton);
554	>	}
555	>
556	>
557	>	} else {
558	>	jets_.refparton_pt[jets_.nref] = -999;
559	>	jets_.refparton_flavor[jets_.nref] = -999;
560	>	}
561	>
562	>
563		}
564
565		jets_.nref++;
#	Line 304 | Line 570 | HiInclusiveJetAnalyzer::analyze(const Ev
570
571		if(isMC_){
572
573	+	edm::Handle<HepMCProduct> hepMCProduct;
574	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
575	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
576	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
577	+	jets_.beamId1 = beamParticles.first->pdg_id();
578	+	jets_.beamId2 = beamParticles.second->pdg_id();
579	+
580		edm::Handle<GenEventInfoProduct> hEventInfo;
581		iEvent.getByLabel(eventInfoTag_,hEventInfo);
582		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 321 | Line 594 | HiInclusiveJetAnalyzer::analyze(const Ev
594		float genjet_pt = genjet.pt();
595
596		// threshold to reduce size of output in minbias PbPb
597	<	if(genjet_pt>20.){
597	>	if(genjet_pt>genPtMin_){
598
599		jets_.genpt[jets_.ngen] = genjet_pt;
600		jets_.geneta[jets_.ngen] = genjet.eta();
601		jets_.genphi[jets_.ngen] = genjet.phi();
602		jets_.geny[jets_.ngen] = genjet.eta();
603
604	+	if(doSubEvent_){
605	+	const GenParticle* gencon = genjet.getGenConstituent(0);
606	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
607	+	}
608
609		// find matching patJet if there is one
610
611		jets_.gendrjt[jets_.ngen] = -1.0;
612		jets_.genmatchindex[jets_.ngen] = -1;
613
337	–
614		for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
615		const pat::Jet& jet = (*jets)[ijet];
616	<
617	<	if(jet.genJet()){
618	<	if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
619	<	fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
620	<	(fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
616	>	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
617	>	if(matchedGenJet){
618	>	// poor man's matching, someone fix please
619	>	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
620	>	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
621	>	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
622
623		jets_.genmatchindex[jets_.ngen] = (int)ijet;
624		jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
#	Line 358 | Line 635 | HiInclusiveJetAnalyzer::analyze(const Ev
635
636		}
637		t->Fill();
638	+	memset(&jets_,0,sizeof jets_);
639		}
640
641
#	Line 457 | Line 735 | inline bool HiInclusiveJetAnalyzer::getP
735		}
736
737
738	+	int
739	+	HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
740	+	{
741	+
742	+	int pfMuonIndex = -1;
743	+	float ptMax = 0.;
744	+
745	+
746	+	for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
747	+	const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
748	+
749	+	int id = pfCandidate.particleId();
750	+	if(abs(id) != 3) continue;
751
752	+	if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
753	+
754	+	double pt =  pfCandidate.pt();
755	+	if(pt>ptMax){
756	+	ptMax = pt;
757	+	pfMuonIndex = (int) icand;
758	+	}
759	+	}
760	+
761	+	return pfMuonIndex;
762	+
763	+	}
764	+
765	+
766	+	double
767	+	HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
768	+	{
769	+
770	+	float lj_x = jet.p4().px();
771	+	float lj_y = jet.p4().py();
772	+	float lj_z = jet.p4().pz();
773	+
774	+	// absolute values squared
775	+	float lj2  = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
776	+	float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
777	+
778	+	// projection vec(mu) to lepjet axis
779	+	float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
780	+
781	+	// absolute value squared and normalized
782	+	float pLrel2 = lepXlj*lepXlj/lj2;
783	+
784	+	// lep2 = pTrel2 + pLrel2
785	+	float pTrel2 = lep2-pLrel2;
786	+
787	+	return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
788	+	}
789	+
790	+	// Recursive function, but this version gets called only the first time
791	+	void
792	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
793	+
794	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
795	+	const reco::Candidate & daughter = *gen.daughter(i);
796	+	double daughterPt = daughter.pt();
797	+	if(daughterPt<1.) continue;
798	+	double daughterEta = daughter.eta();
799	+	if(fabs(daughterEta)>3.) continue;
800	+	int daughterPdgId = daughter.pdgId();
801	+	int daughterStatus = daughter.status();
802	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
803	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
804	+
805	+	// cheesy way of finding strings which were already used
806	+	if(daughter.pdgId()==92){
807	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
808	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
809	+	}
810	+	usedStringPts.push_back(daughter.pt());
811	+	}
812	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
813	+	jets_.bStatus[jets_.bMult] = daughterStatus;
814	+	jets_.bVx[jets_.bMult] = daughter.vx();
815	+	jets_.bVy[jets_.bMult] = daughter.vy();
816	+	jets_.bVz[jets_.bMult] = daughter.vz();
817	+	jets_.bPt[jets_.bMult] = daughterPt;
818	+	jets_.bEta[jets_.bMult] = daughterEta;
819	+	jets_.bPhi[jets_.bMult] = daughter.phi();
820	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
821	+	jets_.bChg[jets_.bMult] = daughter.charge();
822	+	jets_.bMult++;
823	+	saveDaughters(daughter);
824	+	}
825	+	}
826	+
827	+	// This version called for all subsequent calls
828	+	void
829	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
830	+
831	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
832	+	const reco::Candidate & daughter = *gen.daughter(i);
833	+	double daughterPt = daughter.pt();
834	+	if(daughterPt<1.) continue;
835	+	double daughterEta = daughter.eta();
836	+	if(fabs(daughterEta)>3.) continue;
837	+	int daughterPdgId = daughter.pdgId();
838	+	int daughterStatus = daughter.status();
839	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
840	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
841	+
842	+	// cheesy way of finding strings which were already used
843	+	if(daughter.pdgId()==92){
844	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
845	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
846	+	}
847	+	usedStringPts.push_back(daughter.pt());
848	+	}
849	+
850	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
851	+	jets_.bStatus[jets_.bMult] = daughterStatus;
852	+	jets_.bVx[jets_.bMult] = daughter.vx();
853	+	jets_.bVy[jets_.bMult] = daughter.vy();
854	+	jets_.bVz[jets_.bMult] = daughter.vz();
855	+	jets_.bPt[jets_.bMult] = daughterPt;
856	+	jets_.bEta[jets_.bMult] = daughterEta;
857	+	jets_.bPhi[jets_.bMult] = daughter.phi();
858	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
859	+	jets_.bChg[jets_.bMult] = daughter.charge();
860	+	jets_.bMult++;
861	+	saveDaughters(daughter);
862	+	}
863	+	}
864
865		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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