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.6 by yjlee, Wed Nov 2 19:41:13 2011 UTC vs.
Revision 1.28 by yilmaz, Fri Jan 25 17:08:00 2013 UTC

#	Line 19 | Line 19
19
20		#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
21		#include "DataFormats/HeavyIonEvent/interface/Centrality.h"
22	+	#include "Geometry/Records/interface/CaloGeometryRecord.h"
23
24		#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
24	–	#include "DataFormats/PatCandidates/interface/Jet.h"
25		#include "DataFormats/JetReco/interface/CaloJetCollection.h"
26		#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
27		#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
28	+	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
29		#include "DataFormats/JetReco/interface/GenJetCollection.h"
30		#include "DataFormats/VertexReco/interface/Vertex.h"
31		#include "DataFormats/Math/interface/deltaPhi.h"
#	Line 38 | Line 39
39		#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
40		#include "L1Trigger/GlobalTrigger/interface/L1GlobalTrigger.h"
41
42	+	#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
43	+	#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
44	+
45		using namespace std;
46		using namespace edm;
47		using namespace reco;
48
49	<	HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
50	<
49	>	HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) :
50	>	geo(0)
51	>	{
52
53		jetTag_ = iConfig.getParameter<InputTag>("jetTag");
54	<	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
54	>	matchTag_ = iConfig.getUntrackedParameter<InputTag>("matchTag",jetTag_);
55	>
56	>	vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
57	>	trackTag_ = iConfig.getParameter<InputTag>("trackTag");
58	>	useQuality_ = iConfig.getUntrackedParameter<bool>("useQuality",1);
59	>	trackQuality_ = iConfig.getUntrackedParameter<string>("trackQuality","highPurity");
60
61		isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
62	+	fillGenJets_ = iConfig.getUntrackedParameter<bool>("fillGenJets",false);
63	+
64	+	doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
65	+
66	+	rParam = iConfig.getParameter<double>("rParam");
67	+	hardPtMin_ = iConfig.getUntrackedParameter<double>("hardPtMin",4);
68	+	jetPtMin_ = iConfig.getUntrackedParameter<double>("jetPtMin",0);
69
70		if(isMC_){
71		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
#	Line 61 | Line 78 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
78		useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
79		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
80
81	<	if(!isMC_){
81	>	doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
82	>	doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
83	>	saveBfragments_  = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
84	>	skipCorrections_  = iConfig.getUntrackedParameter<bool>("skipCorrections",false);
85	>
86	>	pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
87	>
88	>	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
89	>	EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
90	>	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
91	>	HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
92	>
93	>	genParticleSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("genParticles",edm::InputTag("hiGenParticles"));
94	>
95	>	if(doTrigger_){
96		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
97		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
98
#	Line 78 | Line 109 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
109		}
110
111		cout<<" jet collection : "<<jetTag_<<endl;
112	<	if(isMC_)cout<<" genjet collection : "<<genjetTag_<<endl;
112	>	doSubEvent_ = 0;
113
114	+	if(isMC_){
115	+	cout<<" genjet collection : "<<genjetTag_<<endl;
116	+	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
117	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
118	+	}
119
120
121		}
#	Line 115 | Line 151 | HiInclusiveJetAnalyzer::beginJob() {
151		}
152		if (useCentrality_) {
153		t->Branch("hf",&jets_.hf,"hf/F");
118	–	t->Branch("nref",&jets_.nref,"nref/I");
154		t->Branch("bin",&jets_.bin,"bin/I");
155		}
156	+
157	+	t->Branch("nref",&jets_.nref,"nref/I");
158		t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
159	<	t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
159	>	if(!skipCorrections_) t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
160		t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
161		t->Branch("jty",jets_.jty,"jty[nref]/F");
162		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
163		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
164	+	t->Branch("jtm",jets_.jtm,"jtm[nref]/F");
165	+
166	+	// jet ID information, jet composition
167	+	t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
168	+
169	+	t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
170	+	t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
171	+	t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
172	+	t->Branch("trackHardSum", jets_.trackHardSum,"trackHardSum[nref]/F");
173	+	t->Branch("trackHardN", jets_.trackHardN,"trackHardN[nref]/I");
174	+
175	+	t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
176	+	t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
177	+	t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
178	+	t->Branch("chargedHardSum", jets_.chargedHardSum,"chargedHardSum[nref]/F");
179	+	t->Branch("chargedHardN", jets_.chargedHardN,"chargedHardN[nref]/I");
180	+
181	+	t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
182	+	t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
183	+	t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
184	+	t->Branch("photonHardSum", jets_.photonHardSum,"photonHardSum[nref]/F");
185	+	t->Branch("photonHardN", jets_.photonHardN,"photonHardN[nref]/I");
186	+
187	+	t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
188	+	t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
189	+	t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
190	+
191	+	t->Branch("hcalSum", jets_.hcalSum,"hcalSum[nref]/F");
192	+	t->Branch("ecalSum", jets_.ecalSum,"ecalSum[nref]/F");
193	+
194	+	t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
195	+	t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
196	+	t->Branch("eN", jets_.eN,"eN[nref]/I");
197	+
198	+	t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
199	+	t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
200	+	t->Branch("muN", jets_.muN,"muN[nref]/I");
201	+
202	+	if(!skipCorrections_) t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
203	+	t->Branch("matchedRawPt", jets_.matchedRawPt,"matchedRawPt[nref]/F");
204	+	t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
205	+
206	+	// b-jet discriminators
207	+	if (doLifeTimeTagging_) {
208	+
209	+	t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
210	+	t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
211	+	t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
212	+	t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
213	+	t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
214	+	t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
215	+	t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
216	+	t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
217	+
218	+	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
219	+	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
220	+	t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
221	+	t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
222	+	t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
223	+	t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
224	+
225	+	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
226	+	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
227
228	+	if (doLifeTimeTaggingExtras_) {
229	+	t->Branch("nIP",&jets_.nIP,"nIP/I");
230	+	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
231	+	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
232	+	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
233	+	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
234	+	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
235	+	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
236	+	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
237	+	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
238	+	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
239	+	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
240	+	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
241	+
242	+	}
243	+
244	+	t->Branch("mue",     jets_.mue,     "mue[nref]/F");
245	+	t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
246	+	t->Branch("mueta",   jets_.mueta,   "mueta[nref]/F");
247	+	t->Branch("muphi",   jets_.muphi,   "muphi[nref]/F");
248	+	t->Branch("mudr",    jets_.mudr,    "mudr[nref]/F");
249	+	t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
250	+	t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
251	+	}
252	+
253	+
254		if(isMC_){
255	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
256	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
257	+
258		t->Branch("pthat",&jets_.pthat,"pthat/F");
259
260		// Only matched gen jets
#	Line 137 | Line 266 | HiInclusiveJetAnalyzer::beginJob() {
266		t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
267		// matched parton
268		t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
269	<	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
269	>	t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
270	>	t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
271	>
272	>	t->Branch("genChargedSum", jets_.genChargedSum,"genChargedSum[nref]/F");
273	>	t->Branch("genHardSum", jets_.genHardSum,"genHardSum[nref]/F");
274	>	t->Branch("signalChargedSum", jets_.signalChargedSum,"signalChargedSum[nref]/F");
275	>	t->Branch("signalHardSum", jets_.signalHardSum,"signalHardSum[nref]/F");
276	>
277	>	if(doSubEvent_){
278	>	t->Branch("subid",jets_.subid,"subid[nref]/I");
279	>	}
280	>
281	>	if(fillGenJets_){
282	>	// For all gen jets, matched or unmatched
283	>	t->Branch("ngen",&jets_.ngen,"ngen/I");
284	>	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
285	>	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
286	>	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
287	>	t->Branch("geny",jets_.geny,"geny[ngen]/F");
288	>	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
289	>	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
290	>	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
291	>
292	>	if(doSubEvent_){
293	>	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
294	>	}
295	>	}
296	>
297	>	if(saveBfragments_  ) {
298	>	t->Branch("bMult",&jets_.bMult,"bMult/I");
299	>	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
300	>	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
301	>	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
302	>	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
303	>	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
304	>	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
305	>	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
306	>	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
307	>	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
308	>	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
309	>	}
310
142	–	// For all gen jets, matched or unmatched
143	–	t->Branch("ngen",&jets_.ngen,"ngen/I");
144	–	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
145	–	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
146	–	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
147	–	t->Branch("geny",jets_.geny,"geny[ngen]/F");
148	–	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
149	–	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
150	–	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
311		}
312	<
312	>	/*
313		if(!isMC_){
314		t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
315		t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
#	Line 161 | Line 321 | HiInclusiveJetAnalyzer::beginJob() {
321		t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
322
323		}
324	<
324	>	*/
325		TH1D::SetDefaultSumw2();
326
327
#	Line 186 | Line 346 | HiInclusiveJetAnalyzer::analyze(const Ev
346		double hf = 0.;
347		double b = 999.;
348
349	+	if(!geo){
350	+	edm::ESHandle<CaloGeometry> pGeo;
351	+	iSetup.get<CaloGeometryRecord>().get(pGeo);
352	+	geo = pGeo.product();
353	+	}
354		if(useCentrality_){
355		if(!centrality_) centrality_ = new CentralityProvider(iSetup);
356		centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
#	Line 198 | Line 363 | HiInclusiveJetAnalyzer::analyze(const Ev
363		b = centrality_->bMean();
364		}
365
201	–
202	–
203	–
204	–	//double jetPtMin = 35.;
205	–
206	–
366		// loop the events
367
368		jets_.bin = bin;
369		jets_.hf = hf;
370
371	<
372	<	if (useVtx_) {
373	<	edm::Handle<vector<reco::Vertex> >vertex;
374	<	iEvent.getByLabel(vtxTag_, vertex);
375	<
376	<	if(vertex->size()>0) {
377	<	jets_.vx=vertex->begin()->x();
378	<	jets_.vy=vertex->begin()->y();
379	<	jets_.vz=vertex->begin()->z();
380	<	}
381	<	}
371	>	reco::Vertex::Point vtx(0,0,0);
372	>	if (useVtx_) {
373	>	edm::Handle<vector<reco::Vertex> >vertex;
374	>	iEvent.getByLabel(vtxTag_, vertex);
375	>
376	>	if(vertex->size()>0) {
377	>	jets_.vx=vertex->begin()->x();
378	>	jets_.vy=vertex->begin()->y();
379	>	jets_.vz=vertex->begin()->z();
380	>	vtx = vertex->begin()->position();
381	>	}
382	>	}
383
384		edm::Handle<pat::JetCollection> patjets;
385		if(usePat_)iEvent.getByLabel(jetTag_, patjets);
386	+
387	+	edm::Handle<pat::JetCollection> matchedjets;
388	+	iEvent.getByLabel(matchTag_, matchedjets);
389
390		edm::Handle<reco::JetView> jets;
391		iEvent.getByLabel(jetTag_, jets);
392
393	<	// FILL JRA TREE
393	>	edm::Handle<reco::PFCandidateCollection> pfCandidates;
394	>	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
395	>
396	>	edm::Handle<reco::TrackCollection> tracks;
397	>	iEvent.getByLabel(trackTag_,tracks);
398	>
399	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
400	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
401	>
402	>	edm::Handle<HFRecHitCollection> hfHits;
403	>	edm::Handle<HBHERecHitCollection> hbheHits;
404	>
405	>	iEvent.getByLabel(HcalRecHitHBHESrc_,hbheHits);
406	>	iEvent.getByLabel(HcalRecHitHFSrc_,hfHits);
407	>	iEvent.getByLabel(EBSrc_,ebHits);
408	>	iEvent.getByLabel(EESrc_,eeHits);
409	>
410	>	edm::Handle<reco::GenParticleCollection> genparts;
411	>	iEvent.getByLabel(genParticleSrc_,genparts);
412
413	+
414	+	// FILL JRA TREE
415		jets_.b = b;
416		jets_.nref = 0;
417	+
418
419	<	if(!isMC_){
419	>	if(doTrigger_){
420		fillL1Bits(iEvent);
421		fillHLTBits(iEvent);
422		}
423
424	<	for(unsigned int j = 0 ; j < jets->size(); ++j){
425	<	const reco::Jet& jet = (*jets)[j];
242	<
243	<	//cout<<" jet pt "<<jet.pt()<<endl;
244	<	//if(jet.pt() < jetPtMin) continue;
424	>	for(unsigned int j = 0; j < jets->size(); ++j){
425	>	const reco::Jet& jet = (*jets)[j];
426		if (useJEC_ && usePat_){
427		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
428		}
429	+
430	+	if(doLifeTimeTagging_){
431	+
432	+	if(jetTag_.label()=="icPu5patJets"){
433	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
434	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
435	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
436	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
437	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
438	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
439	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
440	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
441	+	}
442	+	else if(jetTag_.label()=="akPu3PFpatJets"){
443	+	jets_.discr_csvMva[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
444	+	jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
445	+	jets_.discr_muByIp3[jets_.nref]   = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
446	+	jets_.discr_muByPt[jets_.nref]    = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
447	+	jets_.discr_prob[jets_.nref]      = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
448	+	jets_.discr_probb[jets_.nref]     = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
449	+	jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
450	+	jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
451	+	}
452	+	else{
453	+	cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
454	+	}
455	+
456	+	const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
457	+
458	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
459	+
460	+	if (tagInfoSV.nVertices()>0) {
461	+	jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
462	+	// this is the 3d flight distance, for 2-D use (0,true)
463	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
464	+	jets_.svtxdl[jets_.nref]    = m1D.value();
465	+	jets_.svtxdls[jets_.nref]   = m1D.significance();
466	+
467	+	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
468	+	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
469	+	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
470	+	jets_.svtxm[jets_.nref]    = svtx.p4().mass();
471	+	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
472	+	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
473	+	}
474	+
475	+	const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
476	+
477	+	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
478	+	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
479	+
480	+	if (doLifeTimeTaggingExtras_) {
481	+
482	+	TrackRefVector selTracks=tagInfoIP.selectedTracks();
483	+
484	+	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
485	+
486	+	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
487	+	{
488	+	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
489	+	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
490	+	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
491	+	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
492	+	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
493	+	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
494	+	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
495	+	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
496	+	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
497	+	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
498	+	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
499	+	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
500	+	}
501	+
502	+	jets_.nIP += jets_.nselIPtrk[jets_.nref];
503	+
504	+	}
505	+
506	+	const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
507	+	int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
508	+	if(pfMuonIndex >=0){
509	+	const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
510	+	jets_.mupt[jets_.nref]    =  muon.pt();
511	+	jets_.mueta[jets_.nref]   =  muon.eta();
512	+	jets_.muphi[jets_.nref]   =  muon.phi();
513	+	jets_.mue[jets_.nref]     =  muon.energy();
514	+	jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
515	+	jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
516	+	jets_.muchg[jets_.nref]   =  muon.charge();
517	+	}else{
518	+	jets_.mupt[jets_.nref]    =  0.0;
519	+	jets_.mueta[jets_.nref]   =  0.0;
520	+	jets_.muphi[jets_.nref]   =  0.0;
521	+	jets_.mue[jets_.nref]     =  0.0;
522	+	jets_.mudr[jets_.nref]    =  9.9;
523	+	jets_.muptrel[jets_.nref] =  0.0;
524	+	jets_.muchg[jets_.nref]   = 0;
525	+	}
526	+
527	+	}
528	+
529	+	// Jet ID variables
530	+
531	+	jets_.muMax[jets_.nref] = 0;
532	+	jets_.muSum[jets_.nref] = 0;
533	+	jets_.muN[jets_.nref] = 0;
534	+
535	+	jets_.eMax[jets_.nref] = 0;
536	+	jets_.eSum[jets_.nref] = 0;
537	+	jets_.eN[jets_.nref] = 0;
538	+
539	+	jets_.neutralMax[jets_.nref] = 0;
540	+	jets_.neutralSum[jets_.nref] = 0;
541	+	jets_.neutralN[jets_.nref] = 0;
542	+
543	+	jets_.photonMax[jets_.nref] = 0;
544	+	jets_.photonSum[jets_.nref] = 0;
545	+	jets_.photonN[jets_.nref] = 0;
546	+	jets_.photonHardSum[jets_.nref] = 0;
547	+	jets_.photonHardN[jets_.nref] = 0;
548	+
549	+	jets_.chargedMax[jets_.nref] = 0;
550	+	jets_.chargedSum[jets_.nref] = 0;
551	+	jets_.chargedN[jets_.nref] = 0;
552	+	jets_.chargedHardSum[jets_.nref] = 0;
553	+	jets_.chargedHardN[jets_.nref] = 0;
554	+
555	+	jets_.trackMax[jets_.nref] = 0;
556	+	jets_.trackSum[jets_.nref] = 0;
557	+	jets_.trackN[jets_.nref] = 0;
558	+	jets_.trackHardSum[jets_.nref] = 0;
559	+	jets_.trackHardN[jets_.nref] = 0;
560	+
561	+	jets_.hcalSum[jets_.nref] = 0;
562	+	jets_.ecalSum[jets_.nref] = 0;
563	+
564	+	jets_.genChargedSum[jets_.nref] = 0;
565	+	jets_.genHardSum[jets_.nref] = 0;
566	+
567	+	jets_.signalChargedSum[jets_.nref] = 0;
568	+	jets_.signalHardSum[jets_.nref] = 0;
569	+
570	+	jets_.subid[jets_.nref] = -1;
571	+
572	+	for(unsigned int icand = 0; icand < tracks->size(); ++icand){
573	+	const reco::Track& track = (*tracks)[icand];
574	+	if(useQuality_ ){
575	+	bool goodtrack = track.quality(reco::TrackBase::qualityByName(trackQuality_));
576	+	if(!goodtrack) continue;
577	+	}
578	+
579	+	double dr = deltaR(jet,track);
580	+	if(dr < rParam){
581	+	double ptcand = track.pt();
582	+	jets_.trackSum[jets_.nref] += ptcand;
583	+	jets_.trackN[jets_.nref] += 1;
584	+
585	+	if(ptcand > hardPtMin_){
586	+	jets_.trackHardSum[jets_.nref] += ptcand;
587	+	jets_.trackHardN[jets_.nref] += 1;
588	+
589	+	}
590	+	if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
591	+
592	+	}
593	+	}
594	+
595	+	for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
596	+	const reco::PFCandidate& track = (*pfCandidates)[icand];
597	+	double dr = deltaR(jet,track);
598	+	if(dr < rParam){
599	+	double ptcand = track.pt();
600	+	int pfid = track.particleId();
601	+
602	+	switch(pfid){
603	+
604	+	case 1:
605	+	jets_.chargedSum[jets_.nref] += ptcand;
606	+	jets_.chargedN[jets_.nref] += 1;
607	+	if(ptcand > hardPtMin_){
608	+	jets_.chargedHardSum[jets_.nref] += ptcand;
609	+	jets_.chargedHardN[jets_.nref] += 1;
610	+	}
611	+	if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
612	+	break;
613	+
614	+	case 2:
615	+	jets_.eSum[jets_.nref] += ptcand;
616	+	jets_.eN[jets_.nref] += 1;
617	+	if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
618	+	break;
619	+
620	+	case 3:
621	+	jets_.muSum[jets_.nref] += ptcand;
622	+	jets_.muN[jets_.nref] += 1;
623	+	if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
624	+	break;
625	+
626	+	case 4:
627	+	jets_.photonSum[jets_.nref] += ptcand;
628	+	jets_.photonN[jets_.nref] += 1;
629	+	if(ptcand > hardPtMin_){
630	+	jets_.photonHardSum[jets_.nref] += ptcand;
631	+	jets_.photonHardN[jets_.nref] += 1;
632	+	}
633	+	if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
634	+	break;
635	+
636	+	case 5:
637	+	jets_.neutralSum[jets_.nref] += ptcand;
638	+	jets_.neutralN[jets_.nref] += 1;
639	+	if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
640	+	break;
641	+
642	+	default:
643	+	break;
644	+
645	+	}
646	+	}
647	+	}
648	+
649	+	// Calorimeter fractions
650	+
651	+	for(unsigned int i = 0; i < hbheHits->size(); ++i){
652	+	const HBHERecHit & hit= (*hbheHits)[i];
653	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
654	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
655	+	if(dr < rParam){
656	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
657	+	}
658	+	}
659	+
660	+	for(unsigned int i = 0; i < hfHits->size(); ++i){
661	+	const HFRecHit & hit= (*hfHits)[i];
662	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
663	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
664	+	if(dr < rParam){
665	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
666	+	}
667	+	}
668	+
669	+
670	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
671	+	const EcalRecHit & hit= (*ebHits)[i];
672	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
673	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
674	+	if(dr < rParam){
675	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
676	+	}
677	+	}
678	+
679	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
680	+	const EcalRecHit & hit= (*eeHits)[i];
681	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
682	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
683	+	if(dr < rParam){
684	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
685	+	}
686	+	}
687	+
688	+	// Alternative reconstruction matching (PF for calo, calo for PF)
689	+
690	+	double drMin = 100;
691	+	for(unsigned int imatch = 0 ; imatch < matchedjets->size(); ++imatch){
692	+	const pat::Jet& mjet = (*matchedjets)[imatch];
693	+
694	+	double dr = deltaR(jet,mjet);
695	+	if(dr < drMin){
696	+	jets_.matchedPt[jets_.nref] = mjet.pt();
697	+	jets_.matchedRawPt[jets_.nref] = mjet.correctedJet("Uncorrected").pt();
698	+	jets_.matchedR[jets_.nref] = dr;
699	+	drMin = dr;
700	+	}
701	+	}
702	+
703	+	//     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
704	+
705	+	/////////////////////////////////////////////////////////////////
706	+	// Jet core pt^2 discriminant for fake jets
707	+	// Edited by Yue Shi Lai <ylai@mit.edu>
708	+
709	+	// Initial value is 0
710	+	jets_.discr_fr01[jets_.nref] = 0;
711	+	// Start with no directional adaption, i.e. the fake rejection
712	+	// axis is the jet axis
713	+	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
714	+	float azimuth_adapt = jets_.jtphi[jets_.nref];
715	+
716	+	// Unadapted discriminant with adaption search
717	+	for (size_t iteration = 0; iteration < 2; iteration++) {
718	+	float pseudorapidity_adapt_new = pseudorapidity_adapt;
719	+	float azimuth_adapt_new = azimuth_adapt;
720	+	float max_weighted_perp = 0;
721	+	float perp_square_sum = 0;
722	+
723	+	for (size_t index_pf_candidate = 0;
724	+	index_pf_candidate < pfCandidates->size();
725	+	index_pf_candidate++) {
726	+	const reco::PFCandidate &p =
727	+	(*pfCandidates)[index_pf_candidate];
728	+
729	+	switch (p.particleId()) {
730	+	//case 1:     // Charged hadron
731	+	//case 3:     // Muon
732	+	case 4: // Photon
733	+	{
734	+	const float dpseudorapidity =
735	+	p.eta() - pseudorapidity_adapt;
736	+	const float dazimuth =
737	+	reco::deltaPhi(p.phi(), azimuth_adapt);
738	+	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
739	+	// = 50
740	+	const float angular_weight =
741	+	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
742	+	dazimuth * dazimuth));
743	+	const float weighted_perp =
744	+	angular_weight * p.pt() * p.pt();
745	+	const float weighted_perp_square =
746	+	weighted_perp * p.pt();
747	+
748	+	perp_square_sum += weighted_perp_square;
749	+	if (weighted_perp >= max_weighted_perp) {
750	+	pseudorapidity_adapt_new = p.eta();
751	+	azimuth_adapt_new = p.phi();
752	+	max_weighted_perp = weighted_perp;
753	+	}
754	+	}
755	+	default:
756	+	break;
757	+	}
758	+	}
759	+	// Update the fake rejection value
760	+	jets_.discr_fr01[jets_.nref] = std::max(
761	+	jets_.discr_fr01[jets_.nref], perp_square_sum);
762	+	// Update the directional adaption
763	+	pseudorapidity_adapt = pseudorapidity_adapt_new;
764	+	azimuth_adapt = azimuth_adapt_new;
765	+	}
766	+
767		jets_.jtpt[jets_.nref] = jet.pt();
768		jets_.jteta[jets_.nref] = jet.eta();
769		jets_.jtphi[jets_.nref] = jet.phi();
770		jets_.jty[jets_.nref] = jet.eta();
771		jets_.jtpu[jets_.nref] = jet.pileup();
772	+	jets_.jtm[jets_.nref] = jet.mass();
773
774		if(isMC_ && usePat_){
775	<	const reco::GenJet* genjet = (*patjets)[j].genJet();
775	>
776	>	for(UInt_t i = 0; i < genparts->size(); ++i){
777	>	const reco::GenParticle& p = (*genparts)[i];
778	>	if (p.status()!=1) continue;
779	>	if (p.charge()==0) continue;
780	>	double dr = deltaR(jet,p);
781	>	if(dr < rParam){
782	>	double ppt = p.pt();
783	>	jets_.genChargedSum[jets_.nref] += ppt;
784	>	if(ppt > hardPtMin_) jets_.genHardSum[jets_.nref] += ppt;
785	>	if(p.collisionId() == 0){
786	>	jets_.signalChargedSum[jets_.nref] += ppt;
787	>	if(ppt > hardPtMin_) jets_.signalHardSum[jets_.nref] += ppt;
788	>	}
789	>
790	>	}
791	>	}
792	>
793	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
794	>
795		if(genjet){
796		jets_.refpt[jets_.nref] = genjet->pt();
797		jets_.refeta[jets_.nref] = genjet->eta();
#	Line 260 | Line 799 | HiInclusiveJetAnalyzer::analyze(const Ev
799		jets_.refy[jets_.nref] = genjet->eta();
800		jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
801		jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
802	+
803	+	if(doSubEvent_){
804	+	const GenParticle* gencon = genjet->getGenConstituent(0);
805	+	jets_.subid[jets_.nref] = gencon->collisionId();
806	+	}
807	+
808		}else{
809		jets_.refpt[jets_.nref] = -999.;
810		jets_.refeta[jets_.nref] = -999.;
#	Line 269 | Line 814 | HiInclusiveJetAnalyzer::analyze(const Ev
814		jets_.refdrjt[jets_.nref] = -999.;
815		}
816
817	<	// matched partons
818	<	const reco::GenParticle * parton = (*patjets)[j].genParton();
819	<	if(parton){
820	<	jets_.refparton_pt[jets_.nref] = parton->pt();
821	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
822	<	} else {
823	<	jets_.refparton_pt[jets_.nref] = -999;
824	<	jets_.refparton_flavor[jets_.nref] = -999;
825	<	}
817	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
818	>
819	>	// matched partons
820	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
821	>
822	>	if((*patjets)[j].genParton()){
823	>	jets_.refparton_pt[jets_.nref] = parton.pt();
824	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
825	>
826	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
827	>
828	>	usedStringPts.clear();
829	>
830	>	// uncomment this if you want to know the ugly truth about parton matching -matt
831	>	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
832	>	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
833	>
834	>	jets_.bJetIndex[jets_.bMult] = jets_.nref;
835	>	jets_.bStatus[jets_.bMult] = parton.status();
836	>	jets_.bVx[jets_.bMult] = parton.vx();
837	>	jets_.bVy[jets_.bMult] = parton.vy();
838	>	jets_.bVz[jets_.bMult] = parton.vz();
839	>	jets_.bPt[jets_.bMult] = parton.pt();
840	>	jets_.bEta[jets_.bMult] = parton.eta();
841	>	jets_.bPhi[jets_.bMult] = parton.phi();
842	>	jets_.bPdg[jets_.bMult] = parton.pdgId();
843	>	jets_.bChg[jets_.bMult] = parton.charge();
844	>	jets_.bMult++;
845	>	saveDaughters(parton);
846	>	}
847	>
848	>
849	>	} else {
850	>	jets_.refparton_pt[jets_.nref] = -999;
851	>	jets_.refparton_flavor[jets_.nref] = -999;
852	>	}
853	>
854	>
855		}
856
857		jets_.nref++;
#	Line 288 | Line 862 | HiInclusiveJetAnalyzer::analyze(const Ev
862
863		if(isMC_){
864
865	+	edm::Handle<HepMCProduct> hepMCProduct;
866	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
867	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
868	+
869	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
870	+	if(beamParticles.first != 0)jets_.beamId1 = beamParticles.first->pdg_id();
871	+	if(beamParticles.second != 0)jets_.beamId2 = beamParticles.second->pdg_id();
872	+
873		edm::Handle<GenEventInfoProduct> hEventInfo;
874		iEvent.getByLabel(eventInfoTag_,hEventInfo);
875		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 305 | Line 887 | HiInclusiveJetAnalyzer::analyze(const Ev
887		float genjet_pt = genjet.pt();
888
889		// threshold to reduce size of output in minbias PbPb
890	<	if(genjet_pt>20.){
890	>	if(genjet_pt>genPtMin_){
891
892		jets_.genpt[jets_.ngen] = genjet_pt;
893		jets_.geneta[jets_.ngen] = genjet.eta();
894		jets_.genphi[jets_.ngen] = genjet.phi();
895		jets_.geny[jets_.ngen] = genjet.eta();
896
897	+	if(doSubEvent_){
898	+	const GenParticle* gencon = genjet.getGenConstituent(0);
899	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
900	+	}
901
902		// find matching patJet if there is one
903
904		jets_.gendrjt[jets_.ngen] = -1.0;
905		jets_.genmatchindex[jets_.ngen] = -1;
906
321	–
907		for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
908		const pat::Jet& jet = (*jets)[ijet];
909	<
910	<	if(jet.genJet()){
911	<	if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
912	<	fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
913	<	(fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
909	>	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
910	>	if(matchedGenJet){
911	>	// poor man's matching, someone fix please
912	>	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
913	>	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
914	>	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
915
916		jets_.genmatchindex[jets_.ngen] = (int)ijet;
917		jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
#	Line 341 | Line 927 | HiInclusiveJetAnalyzer::analyze(const Ev
927		}
928
929		}
930	+
931		t->Fill();
932	+	memset(&jets_,0,sizeof jets_);
933	+
934		}
935
936
#	Line 441 | Line 1030 | inline bool HiInclusiveJetAnalyzer::getP
1030		}
1031
1032
1033	+	int
1034	+	HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
1035	+	{
1036	+
1037	+	int pfMuonIndex = -1;
1038	+	float ptMax = 0.;
1039	+
1040	+
1041	+	for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
1042	+	const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
1043	+
1044	+	int id = pfCandidate.particleId();
1045	+	if(abs(id) != 3) continue;
1046	+
1047	+	if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
1048	+
1049	+	double pt =  pfCandidate.pt();
1050	+	if(pt>ptMax){
1051	+	ptMax = pt;
1052	+	pfMuonIndex = (int) icand;
1053	+	}
1054	+	}
1055	+
1056	+	return pfMuonIndex;
1057	+
1058	+	}
1059	+
1060	+
1061	+	double
1062	+	HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
1063	+	{
1064	+
1065	+	float lj_x = jet.p4().px();
1066	+	float lj_y = jet.p4().py();
1067	+	float lj_z = jet.p4().pz();
1068	+
1069	+	// absolute values squared
1070	+	float lj2  = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
1071	+	float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
1072	+
1073	+	// projection vec(mu) to lepjet axis
1074	+	float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
1075	+
1076	+	// absolute value squared and normalized
1077	+	float pLrel2 = lepXlj*lepXlj/lj2;
1078	+
1079	+	// lep2 = pTrel2 + pLrel2
1080	+	float pTrel2 = lep2-pLrel2;
1081	+
1082	+	return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
1083	+	}
1084	+
1085	+	// Recursive function, but this version gets called only the first time
1086	+	void
1087	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
1088	+
1089	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
1090	+	const reco::Candidate & daughter = *gen.daughter(i);
1091	+	double daughterPt = daughter.pt();
1092	+	if(daughterPt<1.) continue;
1093	+	double daughterEta = daughter.eta();
1094	+	if(fabs(daughterEta)>3.) continue;
1095	+	int daughterPdgId = daughter.pdgId();
1096	+	int daughterStatus = daughter.status();
1097	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
1098	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
1099	+
1100	+	// cheesy way of finding strings which were already used
1101	+	if(daughter.pdgId()==92){
1102	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
1103	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
1104	+	}
1105	+	usedStringPts.push_back(daughter.pt());
1106	+	}
1107	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
1108	+	jets_.bStatus[jets_.bMult] = daughterStatus;
1109	+	jets_.bVx[jets_.bMult] = daughter.vx();
1110	+	jets_.bVy[jets_.bMult] = daughter.vy();
1111	+	jets_.bVz[jets_.bMult] = daughter.vz();
1112	+	jets_.bPt[jets_.bMult] = daughterPt;
1113	+	jets_.bEta[jets_.bMult] = daughterEta;
1114	+	jets_.bPhi[jets_.bMult] = daughter.phi();
1115	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
1116	+	jets_.bChg[jets_.bMult] = daughter.charge();
1117	+	jets_.bMult++;
1118	+	saveDaughters(daughter);
1119	+	}
1120	+	}
1121	+
1122	+	// This version called for all subsequent calls
1123	+	void
1124	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
1125	+
1126	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
1127	+	const reco::Candidate & daughter = *gen.daughter(i);
1128	+	double daughterPt = daughter.pt();
1129	+	if(daughterPt<1.) continue;
1130	+	double daughterEta = daughter.eta();
1131	+	if(fabs(daughterEta)>3.) continue;
1132	+	int daughterPdgId = daughter.pdgId();
1133	+	int daughterStatus = daughter.status();
1134	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
1135	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
1136	+
1137	+	// cheesy way of finding strings which were already used
1138	+	if(daughter.pdgId()==92){
1139	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
1140	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
1141	+	}
1142	+	usedStringPts.push_back(daughter.pt());
1143	+	}
1144	+
1145	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
1146	+	jets_.bStatus[jets_.bMult] = daughterStatus;
1147	+	jets_.bVx[jets_.bMult] = daughter.vx();
1148	+	jets_.bVy[jets_.bMult] = daughter.vy();
1149	+	jets_.bVz[jets_.bMult] = daughter.vz();
1150	+	jets_.bPt[jets_.bMult] = daughterPt;
1151	+	jets_.bEta[jets_.bMult] = daughterEta;
1152	+	jets_.bPhi[jets_.bMult] = daughter.phi();
1153	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
1154	+	jets_.bChg[jets_.bMult] = daughter.charge();
1155	+	jets_.bMult++;
1156	+	saveDaughters(daughter);
1157	+	}
1158	+	}
1159	+
1160	+	double HiInclusiveJetAnalyzer::getEt(math::XYZPoint pos, double energy){
1161	+	double et = energy*sin(pos.theta());
1162	+	return et;
1163	+	}
1164	+
1165	+	math::XYZPoint HiInclusiveJetAnalyzer::getPosition(const DetId &id, reco::Vertex::Point vtx){
1166	+	const GlobalPoint& pos=geo->getPosition(id);
1167	+	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
1168	+	return posV;
1169	+	}
1170	+
1171
1172
1173		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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