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.7 by yjlee, Wed Nov 2 19:44:09 2011 UTC vs.
Revision 1.8 by mnguyen, Tue Mar 6 15:03: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"
#	Line 46 | Line 45 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
45
46
47		jetTag_ = iConfig.getParameter<InputTag>("jetTag");
48	<	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
48	>	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
49
50		isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
51	<
51	>
52		if(isMC_){
53		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
54		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 61 | Line 60 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
60		useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
61		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
62
63	+	doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
64	+	pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
65	+
66		if(!isMC_){
67		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
68		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
#	Line 126 | Line 128 | HiInclusiveJetAnalyzer::beginJob() {
128		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
129		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
130
131	+	// b-jet discriminators
132	+	if (doLifeTimeTagging_) {
133	+
134	+	t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
135	+	t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
136	+	t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
137	+	t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
138	+	t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
139	+	t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
140	+	t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
141	+	t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
142	+
143	+	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/b");
144	+	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/b");
145	+	t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
146	+	t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
147	+	t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
148	+	t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
149	+
150	+	t->Branch("nIPtracks",jets_.nIPtracks,"nIPTracks[nref]/b");
151	+	t->Branch("nselIPtracks",jets_.nselIPtracks,"nselIPTracks[nref]/b");
152	+
153	+	t->Branch("mue",     jets_.mue,     "mue[nref]/F");
154	+	t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
155	+	t->Branch("mueta",   jets_.mueta,   "mueta[nref]/F");
156	+	t->Branch("muphi",   jets_.muphi,   "muphi[nref]/F");
157	+	t->Branch("mudr",    jets_.mudr,    "mudr[nref]/F");
158	+	t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
159	+	t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
160	+	}
161	+
162		if(isMC_){
163		t->Branch("pthat",&jets_.pthat,"pthat/F");
164
#	Line 138 | Line 171 | HiInclusiveJetAnalyzer::beginJob() {
171		t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
172		// matched parton
173		t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
174	<	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
174	>	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
175	>	t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
176
177		// For all gen jets, matched or unmatched
178		t->Branch("ngen",&jets_.ngen,"ngen/I");
#	Line 228 | Line 262 | HiInclusiveJetAnalyzer::analyze(const Ev
262		edm::Handle<reco::JetView> jets;
263		iEvent.getByLabel(jetTag_, jets);
264
265	+	edm::Handle<reco::PFCandidateCollection> pfCandidates;
266	+	if(doLifeTimeTagging_){
267	+	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
268	+	}
269		// FILL JRA TREE
270
271		jets_.b = b;
#	Line 246 | Line 284 | HiInclusiveJetAnalyzer::analyze(const Ev
284		if (useJEC_ && usePat_){
285		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
286		}
287	+
288	+	if(doLifeTimeTagging_){
289	+
290	+	if(jetTag_.label()=="icPu5patJets"){
291	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
292	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
293	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
294	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
295	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
296	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
297	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
298	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
299	+	}
300	+	else if(jetTag_.label()=="akPu3PFpatJets"){
301	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
302	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
303	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
304	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
305	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
306	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
307	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
308	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
309	+	}
310	+	else{
311	+	cout<<" you fail at btagging "<<endl;
312	+	}
313	+
314	+	const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
315	+	if (tagInfoSV.nVertices()>0) {
316	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
317	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
318	+	jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
319	+	jets_.svtxdl[jets_.nref]    = m1D.value();
320	+	jets_.svtxdls[jets_.nref]   = m1D.significance();
321	+
322	+	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
323	+
324	+	jets_.svtxm[jets_.nref]    = svtx.p4().mass();
325	+	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
326	+
327	+	}
328	+	else {
329	+	jets_.nsvtx[jets_.nref]    =   0;
330	+	jets_.svtxntrk[jets_.nref] =   0;
331	+	jets_.svtxdl[jets_.nref]   = 0.0;
332	+	jets_.svtxdls[jets_.nref]  = 0.0;
333	+	jets_.svtxm[jets_.nref]    = 0.0;
334	+	jets_.svtxpt[jets_.nref]   = 0.0;
335	+	}
336	+	const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
337	+
338	+	jets_.nIPtracks[jets_.nref] = tagInfoIP.tracks().size();
339	+	jets_.nselIPtracks[jets_.nref] = tagInfoIP.selectedTracks().size();
340	+
341	+	// would be nice to save the info for the tracks, but requires a more sophisticated tree
342	+	/*
343	+	//cout << "Jet pt: " << tagInfoIP.jet()->pt() << endl;
344	+	cout << "Tot tracks: " << tagInfoIP.tracks().size() << endl;
345	+	TrackRefVector selTracks=tagInfoIP.selectedTracks();
346	+	int n=selTracks.size();
347	+	cout << "Sel tracks: " << n << endl;
348	+	// false      cout << " Pt  \t d len \t jet dist \t p3d \t p2d\t ip3d \t ip2d " << endl;
349	+	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
350	+	cout << pv << " vs " << tagInfoIP.primaryVertex()->position()   << endl;
351	+	for(int j=0;j< n;j++)
352	+	{
353	+	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[j];
354	+	cout << selTracks[j]->pt() << "\t";
355	+	cout << tagInfoIP.probabilities(0)[j] << "\t";
356	+	cout << tagInfoIP.probabilities(1)[j] << "\t";
357	+	cout << data.ip3d.value() << "\t";
358	+	cout << data.ip3d.significance() << "\t";
359	+	cout << data.distanceToJetAxis.value() << "\t";
360	+	cout << data.distanceToJetAxis.significance() << "\t";
361	+	cout << data.distanceToGhostTrack.value() << "\t";
362	+	cout << data.distanceToGhostTrack.significance() << "\t";
363	+	cout << data.closestToJetAxis << "\t";
364	+	cout << (data.closestToJetAxis - pv).mag() << "\t";
365	+	cout << data.closestToGhostTrack << "\t";
366	+	cout << (data.closestToGhostTrack - pv).mag() << "\t";
367	+	cout << data.ip2d.value() << "\t";
368	+	cout << data.ip2d.significance() <<  endl;
369	+	}
370	+	*/
371	+
372	+	const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
373	+	int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
374	+	if(pfMuonIndex >=0){
375	+	const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
376	+	jets_.mupt[jets_.nref]    =  muon.pt();
377	+	jets_.mueta[jets_.nref]   =  muon.eta();
378	+	jets_.muphi[jets_.nref]   =  muon.phi();
379	+	jets_.mue[jets_.nref]     =  muon.energy();
380	+	jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
381	+	jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
382	+	jets_.muchg[jets_.nref]   =  muon.charge();
383	+	}
384	+	else{
385	+	jets_.mupt[jets_.nref]    =  0.0;
386	+	jets_.mueta[jets_.nref]   =  0.0;
387	+	jets_.muphi[jets_.nref]   =  0.0;
388	+	jets_.mue[jets_.nref]     =  0.0;
389	+	jets_.mudr[jets_.nref]    =  9.9;
390	+	jets_.muptrel[jets_.nref] =  0.0;
391	+	jets_.muchg[jets_.nref]   = 0;
392	+	}
393	+	}
394	+
395		jets_.jtpt[jets_.nref] = jet.pt();
396		jets_.jteta[jets_.nref] = jet.eta();
397		jets_.jtphi[jets_.nref] = jet.phi();
#	Line 253 | Line 399 | HiInclusiveJetAnalyzer::analyze(const Ev
399		jets_.jtpu[jets_.nref] = jet.pileup();
400
401		if(isMC_ && usePat_){
402	+	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
403		const reco::GenJet* genjet = (*patjets)[j].genJet();
404	+
405		if(genjet){
406		jets_.refpt[jets_.nref] = genjet->pt();
407		jets_.refeta[jets_.nref] = genjet->eta();
#	Line 273 | Line 421 | HiInclusiveJetAnalyzer::analyze(const Ev
421		// matched partons
422		const reco::GenParticle * parton = (*patjets)[j].genParton();
423		if(parton){
424	<	jets_.refparton_pt[jets_.nref] = parton->pt();
425	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
426	<	} else {
424	>	jets_.refparton_pt[jets_.nref] = parton->pt();
425	>	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
426	>	} else {
427		jets_.refparton_pt[jets_.nref] = -999;
428		jets_.refparton_flavor[jets_.nref] = -999;
429	<	}
429	>	}
430		}
431
432		jets_.nref++;
#	Line 319 | Line 467 | HiInclusiveJetAnalyzer::analyze(const Ev
467		jets_.gendrjt[jets_.ngen] = -1.0;
468		jets_.genmatchindex[jets_.ngen] = -1;
469
322	–
470		for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
471		const pat::Jet& jet = (*jets)[ijet];
472	<
473	<	if(jet.genJet()){
474	<	if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
475	<	fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
476	<	(fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
472	>	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
473	>	if(matchedGenJet){
474	>	// poor man's matching, someone fix please
475	>	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
476	>	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
477	>	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
478
479		jets_.genmatchindex[jets_.ngen] = (int)ijet;
480		jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
#	Line 442 | Line 590 | inline bool HiInclusiveJetAnalyzer::getP
590		}
591
592
593	+	int
594	+	HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
595	+	{
596	+
597	+	int pfMuonIndex = -1;
598	+	float ptMax = 0.;
599	+
600	+
601	+	for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
602	+	const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
603	+
604	+	int id = pfCandidate.particleId();
605	+	if(abs(id) != 3) continue;
606
607	+	if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
608	+
609	+	double pt =  pfCandidate.pt();
610	+	if(pt>ptMax){
611	+	ptMax = pt;
612	+	pfMuonIndex = (int) icand;
613	+	}
614	+	}
615	+
616	+	return pfMuonIndex;
617	+
618	+	}
619	+
620	+
621	+	double
622	+	HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
623	+	{
624	+
625	+	float lj_x = jet.p4().px();
626	+	float lj_y = jet.p4().py();
627	+	float lj_z = jet.p4().pz();
628	+
629	+	// absolute values squared
630	+	float lj2  = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
631	+	float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
632	+
633	+	// projection vec(mu) to lepjet axis
634	+	float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
635	+
636	+	// absolute value squared and normalized
637	+	float pLrel2 = lepXlj*lepXlj/lj2;
638	+
639	+	// lep2 = pTrel2 + pLrel2
640	+	float pTrel2 = lep2-pLrel2;
641	+
642	+	return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
643	+	}
644
645		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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