ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc

Comparing UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc (file contents):
Revision 1.9 by yilmaz, Sun Apr 22 19:12:44 2012 UTC vs.
Revision 1.35 by yjlee, Sat Jan 26 17:06:06 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"
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 37 | 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	+	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	<
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");
72		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 61 | Line 79 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
79		usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
80
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	<	if(!isMC_){
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 79 | Line 108 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
108		}
109		}
110
111	<	cout<<" jet collection : "<<jetTag_<<endl;
111	>	//  cout<<" jet collection : "<<jetTag_<<endl;
112	>	doSubEvent_ = 0;
113	>
114		if(isMC_){
115	<	cout<<" genjet collection : "<<genjetTag_<<endl;
115	>	//     cout<<" genjet collection : "<<genjetTag_<<endl;
116		genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
117	+	doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
118		}
119
120
#	Line 124 | Line 156 | HiInclusiveJetAnalyzer::beginJob() {
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	+
203	+	t->Branch("fHPD",jets_.fHPD,"fHPD[nref]/F");
204	+	t->Branch("fRBX",jets_.fRBX,"fRBX[nref]/F");
205	+	t->Branch("n90",jets_.n90,"n90[nref]/I");
206	+	t->Branch("fSubDet1",jets_.fSubDet1,"fSubDet1[nref]/F");
207	+	t->Branch("fSubDet2",jets_.fSubDet2,"fSubDet2[nref]/F");
208	+	t->Branch("fSubDet3",jets_.fSubDet3,"fSubDet3[nref]/F");
209	+	t->Branch("fSubDet4",jets_.fSubDet4,"fSubDet4[nref]/F");
210	+	t->Branch("restrictedEMF",jets_.restrictedEMF,"restrictedEMF[nref]/F");
211	+	t->Branch("nHCAL",jets_.nHCAL,"nHCAL[nref]/I");
212	+	t->Branch("nECAL",jets_.nECAL,"nECAL[nref]/I");
213	+	t->Branch("apprHPD",jets_.apprHPD,"apprHPD[nref]/F");
214	+	t->Branch("apprRBX",jets_.apprRBX,"apprRBX[nref]/F");
215	+
216	+	//  t->Branch("hitsInN90",jets_.n90,"hitsInN90[nref]");
217	+	t->Branch("n2RPC",jets_.n2RPC,"n2RPC[nref]/I");
218	+	t->Branch("n3RPC",jets_.n3RPC,"n3RPC[nref]/I");
219	+	t->Branch("nRPC",jets_.nRPC,"nRPC[nref]/I");
220	+
221	+	t->Branch("fEB",jets_.fEB,"fEB[nref]/F");
222	+	t->Branch("fEE",jets_.fEE,"fEE[nref]/F");
223	+	t->Branch("fHB",jets_.fHB,"fHB[nref]/F");
224	+	t->Branch("fHE",jets_.fHE,"fHE[nref]/F");
225	+	t->Branch("fHO",jets_.fHO,"fHO[nref]/F");
226	+	t->Branch("fLong",jets_.fLong,"fLong[nref]/F");
227	+	t->Branch("fShort",jets_.fShort,"fShort[nref]/F");
228	+	t->Branch("fLS",jets_.fLS,"fLS[nref]/F");
229	+	t->Branch("fHFOOT",jets_.fHFOOT,"fHFOOT[nref]/F");
230	+
231	+
232	+	// Jet ID
233	+
234	+	if(!skipCorrections_) t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
235	+	t->Branch("matchedRawPt", jets_.matchedRawPt,"matchedRawPt[nref]/F");
236	+	t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
237
238		// b-jet discriminators
239		if (doLifeTimeTagging_) {
#	Line 142 | Line 247 | HiInclusiveJetAnalyzer::beginJob() {
247		t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
248		t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
249
250	<	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/b");
251	<	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/b");
250	>	t->Branch("nsvtx",    jets_.nsvtx,    "nsvtx[nref]/I");
251	>	t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
252		t->Branch("svtxdl",   jets_.svtxdl,   "svtxdl[nref]/F");
253		t->Branch("svtxdls",  jets_.svtxdls,  "svtxdls[nref]/F");
254		t->Branch("svtxm",    jets_.svtxm,    "svtxm[nref]/F");
255		t->Branch("svtxpt",   jets_.svtxpt,   "svtxpt[nref]/F");
256
257	<	t->Branch("nIPtracks",jets_.nIPtracks,"nIPTracks[nref]/b");
258	<	t->Branch("nselIPtracks",jets_.nselIPtracks,"nselIPTracks[nref]/b");
257	>	t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
258	>	t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
259	>
260	>	if (doLifeTimeTaggingExtras_) {
261	>	t->Branch("nIP",&jets_.nIP,"nIP/I");
262	>	t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
263	>	t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
264	>	t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
265	>	t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
266	>	t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
267	>	t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
268	>	t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
269	>	t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
270	>	t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
271	>	t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
272	>	t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
273	>
274	>	}
275
276		t->Branch("mue",     jets_.mue,     "mue[nref]/F");
277		t->Branch("mupt",    jets_.mupt,    "mupt[nref]/F");
#	Line 160 | Line 281 | HiInclusiveJetAnalyzer::beginJob() {
281		t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
282		t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
283		}
284	+
285
286		if(isMC_){
287	+	t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
288	+	t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
289	+
290		t->Branch("pthat",&jets_.pthat,"pthat/F");
291
292		// Only matched gen jets
#	Line 176 | Line 301 | HiInclusiveJetAnalyzer::beginJob() {
301		t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
302		t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
303
304	<	// For all gen jets, matched or unmatched
305	<	t->Branch("ngen",&jets_.ngen,"ngen/I");
306	<	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
307	<	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
308	<	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
309	<	t->Branch("geny",jets_.geny,"geny[ngen]/F");
310	<	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
311	<	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
312	<	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
304	>	t->Branch("genChargedSum", jets_.genChargedSum,"genChargedSum[nref]/F");
305	>	t->Branch("genHardSum", jets_.genHardSum,"genHardSum[nref]/F");
306	>	t->Branch("signalChargedSum", jets_.signalChargedSum,"signalChargedSum[nref]/F");
307	>	t->Branch("signalHardSum", jets_.signalHardSum,"signalHardSum[nref]/F");
308	>
309	>	if(doSubEvent_){
310	>	t->Branch("subid",jets_.subid,"subid[nref]/I");
311	>	}
312	>
313	>	if(fillGenJets_){
314	>	// For all gen jets, matched or unmatched
315	>	t->Branch("ngen",&jets_.ngen,"ngen/I");
316	>	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
317	>	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
318	>	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
319	>	t->Branch("geny",jets_.geny,"geny[ngen]/F");
320	>	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
321	>	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
322	>	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
323	>
324	>	if(doSubEvent_){
325	>	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
326	>	}
327	>	}
328	>
329	>	if(saveBfragments_  ) {
330	>	t->Branch("bMult",&jets_.bMult,"bMult/I");
331	>	t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
332	>	t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
333	>	t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
334	>	t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
335	>	t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
336	>	t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
337	>	t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
338	>	t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
339	>	t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
340	>	t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
341	>	}
342	>
343		}
344		/*
345		if(!isMC_){
#	Line 223 | Line 378 | HiInclusiveJetAnalyzer::analyze(const Ev
378		double hf = 0.;
379		double b = 999.;
380
381	+	if(!geo){
382	+	edm::ESHandle<CaloGeometry> pGeo;
383	+	iSetup.get<CaloGeometryRecord>().get(pGeo);
384	+	geo = pGeo.product();
385	+	}
386		if(useCentrality_){
387		if(!centrality_) centrality_ = new CentralityProvider(iSetup);
388		centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
#	Line 235 | Line 395 | HiInclusiveJetAnalyzer::analyze(const Ev
395		b = centrality_->bMean();
396		}
397
238	–
239	–
240	–
241	–	//double jetPtMin = 35.;
242	–
243	–
398		// loop the events
399
400		jets_.bin = bin;
401		jets_.hf = hf;
402
403	<
404	<	if (useVtx_) {
405	<	edm::Handle<vector<reco::Vertex> >vertex;
406	<	iEvent.getByLabel(vtxTag_, vertex);
407	<
408	<	if(vertex->size()>0) {
409	<	jets_.vx=vertex->begin()->x();
410	<	jets_.vy=vertex->begin()->y();
411	<	jets_.vz=vertex->begin()->z();
412	<	}
413	<	}
403	>	reco::Vertex::Point vtx(0,0,0);
404	>	if (useVtx_) {
405	>	edm::Handle<vector<reco::Vertex> >vertex;
406	>	iEvent.getByLabel(vtxTag_, vertex);
407	>
408	>	if(vertex->size()>0) {
409	>	jets_.vx=vertex->begin()->x();
410	>	jets_.vy=vertex->begin()->y();
411	>	jets_.vz=vertex->begin()->z();
412	>	vtx = vertex->begin()->position();
413	>	}
414	>	}
415
416		edm::Handle<pat::JetCollection> patjets;
417		if(usePat_)iEvent.getByLabel(jetTag_, patjets);
418	+
419	+	edm::Handle<pat::JetCollection> matchedjets;
420	+	iEvent.getByLabel(matchTag_, matchedjets);
421
422		edm::Handle<reco::JetView> jets;
423		iEvent.getByLabel(jetTag_, jets);
424
425		edm::Handle<reco::PFCandidateCollection> pfCandidates;
426	<	if(doLifeTimeTagging_){
427	<	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
428	<	}
429	<	// FILL JRA TREE
426	>	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
427	>
428	>	edm::Handle<reco::TrackCollection> tracks;
429	>	iEvent.getByLabel(trackTag_,tracks);
430	>
431	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
432	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
433	>
434	>	edm::Handle<HFRecHitCollection> hfHits;
435	>	edm::Handle<HBHERecHitCollection> hbheHits;
436	>
437	>	iEvent.getByLabel(HcalRecHitHBHESrc_,hbheHits);
438	>	iEvent.getByLabel(HcalRecHitHFSrc_,hfHits);
439	>	iEvent.getByLabel(EBSrc_,ebHits);
440	>	iEvent.getByLabel(EESrc_,eeHits);
441	>
442	>	edm::Handle<reco::GenParticleCollection> genparts;
443	>	iEvent.getByLabel(genParticleSrc_,genparts);
444	>
445
446	+	// FILL JRA TREE
447		jets_.b = b;
448		jets_.nref = 0;
449	+
450
451	<	if(!isMC_){
451	>	if(doTrigger_){
452		fillL1Bits(iEvent);
453		fillHLTBits(iEvent);
454		}
455
456	<	for(unsigned int j = 0 ; j < jets->size(); ++j){
457	<	const reco::Jet& jet = (*jets)[j];
458	<
459	<	//cout<<" jet pt "<<jet.pt()<<endl;
285	<	//if(jet.pt() < jetPtMin) continue;
456	>	for(unsigned int j = 0; j < jets->size(); ++j){
457	>	const reco::Jet& jet = (*jets)[j];
458	>
459	>	if(jet.pt() < jetPtMin_) continue;
460		if (useJEC_ && usePat_){
461		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
462		}
#	Line 310 | Line 484 | HiInclusiveJetAnalyzer::analyze(const Ev
484		jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
485		}
486		else{
487	<	cout<<" you fail at btagging "<<endl;
487	>	//      cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
488		}
489	<
489	>
490		const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
491	+
492	+	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
493	+
494		if (tagInfoSV.nVertices()>0) {
318	–	Measurement1D m1D = tagInfoSV.flightDistance(0);
319	–	jets_.nsvtx[jets_.nref]     = tagInfoSV.nVertices();
495		jets_.svtxntrk[jets_.nref]  = tagInfoSV.nVertexTracks(0);
496	+	// this is the 3d flight distance, for 2-D use (0,true)
497	+	Measurement1D m1D = tagInfoSV.flightDistance(0);
498		jets_.svtxdl[jets_.nref]    = m1D.value();
499		jets_.svtxdls[jets_.nref]   = m1D.significance();
500
501	<	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
502	<
501	>	const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
502	>	//cout<<" SV:  vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
503	>	//cout<<" PV:  vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
504		jets_.svtxm[jets_.nref]    = svtx.p4().mass();
505	<	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
506	<
329	<	}
330	<	else {
331	<	jets_.nsvtx[jets_.nref]    =   0;
332	<	jets_.svtxntrk[jets_.nref] =   0;
333	<	jets_.svtxdl[jets_.nref]   = 0.0;
334	<	jets_.svtxdls[jets_.nref]  = 0.0;
335	<	jets_.svtxm[jets_.nref]    = 0.0;
336	<	jets_.svtxpt[jets_.nref]   = 0.0;
505	>	jets_.svtxpt[jets_.nref]   = svtx.p4().pt();
506	>	//cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
507		}
508	+
509		const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
510	+
511	+	jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
512	+	jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
513
514	<	jets_.nIPtracks[jets_.nref] = tagInfoIP.tracks().size();
341	<	jets_.nselIPtracks[jets_.nref] = tagInfoIP.selectedTracks().size();
514	>	if (doLifeTimeTaggingExtras_) {
515
516	<	// would be nice to save the info for the tracks, but requires a more sophisticated tree
517	<	/*
518	<	//cout << "Jet pt: " << tagInfoIP.jet()->pt() << endl;
519	<	cout << "Tot tracks: " << tagInfoIP.tracks().size() << endl;
520	<	TrackRefVector selTracks=tagInfoIP.selectedTracks();
521	<	int n=selTracks.size();
522	<	cout << "Sel tracks: " << n << endl;
523	<	// false      cout << " Pt  \t d len \t jet dist \t p3d \t p2d\t ip3d \t ip2d " << endl;
524	<	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
525	<	cout << pv << " vs " << tagInfoIP.primaryVertex()->position()   << endl;
526	<	for(int j=0;j< n;j++)
527	<	{
528	<	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[j];
529	<	cout << selTracks[j]->pt() << "\t";
530	<	cout << tagInfoIP.probabilities(0)[j] << "\t";
531	<	cout << tagInfoIP.probabilities(1)[j] << "\t";
532	<	cout << data.ip3d.value() << "\t";
533	<	cout << data.ip3d.significance() << "\t";
534	<	cout << data.distanceToJetAxis.value() << "\t";
535	<	cout << data.distanceToJetAxis.significance() << "\t";
536	<	cout << data.distanceToGhostTrack.value() << "\t";
364	<	cout << data.distanceToGhostTrack.significance() << "\t";
365	<	cout << data.closestToJetAxis << "\t";
366	<	cout << (data.closestToJetAxis - pv).mag() << "\t";
367	<	cout << data.closestToGhostTrack << "\t";
368	<	cout << (data.closestToGhostTrack - pv).mag() << "\t";
369	<	cout << data.ip2d.value() << "\t";
370	<	cout << data.ip2d.significance() <<  endl;
371	<	}
372	<	*/
516	>	TrackRefVector selTracks=tagInfoIP.selectedTracks();
517	>
518	>	GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
519	>
520	>	for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
521	>	{
522	>	jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
523	>	TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
524	>	jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
525	>	jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
526	>	jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
527	>	jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
528	>	jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
529	>	jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
530	>	jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
531	>	jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
532	>	jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
533	>	jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag();  //decay length
534	>	}
535	>
536	>	jets_.nIP += jets_.nselIPtrk[jets_.nref];
537
538	+	}
539	+
540		const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
541		int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
542		if(pfMuonIndex >=0){
#	Line 382 | Line 548 | HiInclusiveJetAnalyzer::analyze(const Ev
548		jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
549		jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
550		jets_.muchg[jets_.nref]   =  muon.charge();
551	<	}
386	<	else{
551	>	}else{
552		jets_.mupt[jets_.nref]    =  0.0;
553		jets_.mueta[jets_.nref]   =  0.0;
554		jets_.muphi[jets_.nref]   =  0.0;
#	Line 392 | Line 557 | HiInclusiveJetAnalyzer::analyze(const Ev
557		jets_.muptrel[jets_.nref] =  0.0;
558		jets_.muchg[jets_.nref]   = 0;
559		}
560	+
561		}
562
563	+	// Jet ID variables
564	+
565	+	jets_.muMax[jets_.nref] = 0;
566	+	jets_.muSum[jets_.nref] = 0;
567	+	jets_.muN[jets_.nref] = 0;
568	+
569	+	jets_.eMax[jets_.nref] = 0;
570	+	jets_.eSum[jets_.nref] = 0;
571	+	jets_.eN[jets_.nref] = 0;
572	+
573	+	jets_.neutralMax[jets_.nref] = 0;
574	+	jets_.neutralSum[jets_.nref] = 0;
575	+	jets_.neutralN[jets_.nref] = 0;
576	+
577	+	jets_.photonMax[jets_.nref] = 0;
578	+	jets_.photonSum[jets_.nref] = 0;
579	+	jets_.photonN[jets_.nref] = 0;
580	+	jets_.photonHardSum[jets_.nref] = 0;
581	+	jets_.photonHardN[jets_.nref] = 0;
582	+
583	+	jets_.chargedMax[jets_.nref] = 0;
584	+	jets_.chargedSum[jets_.nref] = 0;
585	+	jets_.chargedN[jets_.nref] = 0;
586	+	jets_.chargedHardSum[jets_.nref] = 0;
587	+	jets_.chargedHardN[jets_.nref] = 0;
588	+
589	+	jets_.trackMax[jets_.nref] = 0;
590	+	jets_.trackSum[jets_.nref] = 0;
591	+	jets_.trackN[jets_.nref] = 0;
592	+	jets_.trackHardSum[jets_.nref] = 0;
593	+	jets_.trackHardN[jets_.nref] = 0;
594	+
595	+	jets_.hcalSum[jets_.nref] = 0;
596	+	jets_.ecalSum[jets_.nref] = 0;
597	+
598	+	jets_.genChargedSum[jets_.nref] = 0;
599	+	jets_.genHardSum[jets_.nref] = 0;
600	+
601	+	jets_.signalChargedSum[jets_.nref] = 0;
602	+	jets_.signalHardSum[jets_.nref] = 0;
603	+
604	+	jets_.subid[jets_.nref] = -1;
605	+
606	+	for(unsigned int icand = 0; icand < tracks->size(); ++icand){
607	+	const reco::Track& track = (*tracks)[icand];
608	+	if(useQuality_ ){
609	+	bool goodtrack = track.quality(reco::TrackBase::qualityByName(trackQuality_));
610	+	if(!goodtrack) continue;
611	+	}
612	+
613	+	double dr = deltaR(jet,track);
614	+	if(dr < rParam){
615	+	double ptcand = track.pt();
616	+	jets_.trackSum[jets_.nref] += ptcand;
617	+	jets_.trackN[jets_.nref] += 1;
618	+
619	+	if(ptcand > hardPtMin_){
620	+	jets_.trackHardSum[jets_.nref] += ptcand;
621	+	jets_.trackHardN[jets_.nref] += 1;
622	+
623	+	}
624	+	if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
625	+
626	+	}
627	+	}
628	+
629	+	for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
630	+	const reco::PFCandidate& track = (*pfCandidates)[icand];
631	+	double dr = deltaR(jet,track);
632	+	if(dr < rParam){
633	+	double ptcand = track.pt();
634	+	int pfid = track.particleId();
635	+
636	+	switch(pfid){
637	+
638	+	case 1:
639	+	jets_.chargedSum[jets_.nref] += ptcand;
640	+	jets_.chargedN[jets_.nref] += 1;
641	+	if(ptcand > hardPtMin_){
642	+	jets_.chargedHardSum[jets_.nref] += ptcand;
643	+	jets_.chargedHardN[jets_.nref] += 1;
644	+	}
645	+	if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
646	+	break;
647	+
648	+	case 2:
649	+	jets_.eSum[jets_.nref] += ptcand;
650	+	jets_.eN[jets_.nref] += 1;
651	+	if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
652	+	break;
653	+
654	+	case 3:
655	+	jets_.muSum[jets_.nref] += ptcand;
656	+	jets_.muN[jets_.nref] += 1;
657	+	if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
658	+	break;
659	+
660	+	case 4:
661	+	jets_.photonSum[jets_.nref] += ptcand;
662	+	jets_.photonN[jets_.nref] += 1;
663	+	if(ptcand > hardPtMin_){
664	+	jets_.photonHardSum[jets_.nref] += ptcand;
665	+	jets_.photonHardN[jets_.nref] += 1;
666	+	}
667	+	if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
668	+	break;
669	+
670	+	case 5:
671	+	jets_.neutralSum[jets_.nref] += ptcand;
672	+	jets_.neutralN[jets_.nref] += 1;
673	+	if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
674	+	break;
675	+
676	+	default:
677	+	break;
678	+
679	+	}
680	+	}
681	+	}
682	+
683	+	// Calorimeter fractions
684	+
685	+	for(unsigned int i = 0; i < hbheHits->size(); ++i){
686	+	const HBHERecHit & hit= (*hbheHits)[i];
687	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
688	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
689	+	if(dr < rParam){
690	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
691	+	}
692	+	}
693	+
694	+	for(unsigned int i = 0; i < hfHits->size(); ++i){
695	+	const HFRecHit & hit= (*hfHits)[i];
696	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
697	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
698	+	if(dr < rParam){
699	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
700	+	}
701	+	}
702	+
703	+
704	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
705	+	const EcalRecHit & hit= (*ebHits)[i];
706	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
707	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
708	+	if(dr < rParam){
709	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
710	+	}
711	+	}
712	+
713	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
714	+	const EcalRecHit & hit= (*eeHits)[i];
715	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
716	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
717	+	if(dr < rParam){
718	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
719	+	}
720	+	}
721	+
722	+	// Jet ID for CaloJets
723	+
724	+
725	+
726	+
727	+
728	+
729	+
730	+	// Alternative reconstruction matching (PF for calo, calo for PF)
731	+
732	+	double drMin = 100;
733	+	for(unsigned int imatch = 0 ; imatch < matchedjets->size(); ++imatch){
734	+	const pat::Jet& mjet = (*matchedjets)[imatch];
735	+
736	+	double dr = deltaR(jet,mjet);
737	+	if(dr < drMin){
738	+	jets_.matchedPt[jets_.nref] = mjet.pt();
739	+	jets_.matchedRawPt[jets_.nref] = mjet.correctedJet("Uncorrected").pt();
740	+	jets_.matchedR[jets_.nref] = dr;
741	+	drMin = dr;
742	+	}
743	+	}
744	+
745	+	//     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
746	+
747	+	/////////////////////////////////////////////////////////////////
748	+	// Jet core pt^2 discriminant for fake jets
749	+	// Edited by Yue Shi Lai <ylai@mit.edu>
750	+
751	+	// Initial value is 0
752	+	jets_.discr_fr01[jets_.nref] = 0;
753	+	// Start with no directional adaption, i.e. the fake rejection
754	+	// axis is the jet axis
755	+	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
756	+	float azimuth_adapt = jets_.jtphi[jets_.nref];
757	+
758	+	// Unadapted discriminant with adaption search
759	+	for (size_t iteration = 0; iteration < 2; iteration++) {
760	+	float pseudorapidity_adapt_new = pseudorapidity_adapt;
761	+	float azimuth_adapt_new = azimuth_adapt;
762	+	float max_weighted_perp = 0;
763	+	float perp_square_sum = 0;
764	+
765	+	for (size_t index_pf_candidate = 0;
766	+	index_pf_candidate < pfCandidates->size();
767	+	index_pf_candidate++) {
768	+	const reco::PFCandidate &p =
769	+	(*pfCandidates)[index_pf_candidate];
770	+
771	+	switch (p.particleId()) {
772	+	//case 1:     // Charged hadron
773	+	//case 3:     // Muon
774	+	case 4: // Photon
775	+	{
776	+	const float dpseudorapidity =
777	+	p.eta() - pseudorapidity_adapt;
778	+	const float dazimuth =
779	+	reco::deltaPhi(p.phi(), azimuth_adapt);
780	+	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
781	+	// = 50
782	+	const float angular_weight =
783	+	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
784	+	dazimuth * dazimuth));
785	+	const float weighted_perp =
786	+	angular_weight * p.pt() * p.pt();
787	+	const float weighted_perp_square =
788	+	weighted_perp * p.pt();
789	+
790	+	perp_square_sum += weighted_perp_square;
791	+	if (weighted_perp >= max_weighted_perp) {
792	+	pseudorapidity_adapt_new = p.eta();
793	+	azimuth_adapt_new = p.phi();
794	+	max_weighted_perp = weighted_perp;
795	+	}
796	+	}
797	+	default:
798	+	break;
799	+	}
800	+	}
801	+	// Update the fake rejection value
802	+	jets_.discr_fr01[jets_.nref] = std::max(
803	+	jets_.discr_fr01[jets_.nref], perp_square_sum);
804	+	// Update the directional adaption
805	+	pseudorapidity_adapt = pseudorapidity_adapt_new;
806	+	azimuth_adapt = azimuth_adapt_new;
807	+	}
808	+
809		jets_.jtpt[jets_.nref] = jet.pt();
810		jets_.jteta[jets_.nref] = jet.eta();
811		jets_.jtphi[jets_.nref] = jet.phi();
812		jets_.jty[jets_.nref] = jet.eta();
813		jets_.jtpu[jets_.nref] = jet.pileup();
814	+	jets_.jtm[jets_.nref] = jet.mass();
815
816	+	if(usePat_){
817	+
818	+	jets_.fHPD[jets_.nref] = (*patjets)[j].jetID().fHPD;
819	+	jets_.fRBX[jets_.nref] = (*patjets)[j].jetID().fRBX;
820	+	jets_.n90[jets_.nref] = (*patjets)[j].n90();
821	+
822	+	jets_.fSubDet1[jets_.nref] = (*patjets)[j].jetID().fSubDetector1;
823	+	jets_.fSubDet2[jets_.nref] = (*patjets)[j].jetID().fSubDetector2;
824	+	jets_.fSubDet3[jets_.nref] = (*patjets)[j].jetID().fSubDetector3;
825	+	jets_.fSubDet4[jets_.nref] = (*patjets)[j].jetID().fSubDetector4;
826	+	jets_.restrictedEMF[jets_.nref] = (*patjets)[j].jetID().restrictedEMF;
827	+	jets_.nHCAL[jets_.nref] = (*patjets)[j].jetID().nHCALTowers;
828	+	jets_.nECAL[jets_.nref] = (*patjets)[j].jetID().nECALTowers;
829	+	jets_.apprHPD[jets_.nref] = (*patjets)[j].jetID().approximatefHPD;
830	+	jets_.apprRBX[jets_.nref] = (*patjets)[j].jetID().approximatefRBX;
831	+
832	+	//       jets_.n90[jets_.nref] = (*patjets)[j].jetID().hitsInN90;
833	+	jets_.n2RPC[jets_.nref] = (*patjets)[j].jetID().numberOfHits2RPC;
834	+	jets_.n3RPC[jets_.nref] = (*patjets)[j].jetID().numberOfHits3RPC;
835	+	jets_.nRPC[jets_.nref] = (*patjets)[j].jetID().numberOfHitsRPC;
836	+
837	+	jets_.fEB[jets_.nref] = (*patjets)[j].jetID().fEB;
838	+	jets_.fEE[jets_.nref] = (*patjets)[j].jetID().fEE;
839	+	jets_.fHB[jets_.nref] = (*patjets)[j].jetID().fHB;
840	+	jets_.fHE[jets_.nref] = (*patjets)[j].jetID().fHE;
841	+	jets_.fHO[jets_.nref] = (*patjets)[j].jetID().fHO;
842	+	jets_.fLong[jets_.nref] = (*patjets)[j].jetID().fLong;
843	+	jets_.fShort[jets_.nref] = (*patjets)[j].jetID().fShort;
844	+	jets_.fLS[jets_.nref] = (*patjets)[j].jetID().fLS;
845	+	jets_.fHFOOT[jets_.nref] = (*patjets)[j].jetID().fHFOOT;
846	+
847	+
848	+	}
849	+
850	+	if(isMC_){
851	+
852	+	for(UInt_t i = 0; i < genparts->size(); ++i){
853	+	const reco::GenParticle& p = (*genparts)[i];
854	+	if (p.status()!=1) continue;
855	+	if (p.charge()==0) continue;
856	+	double dr = deltaR(jet,p);
857	+	if(dr < rParam){
858	+	double ppt = p.pt();
859	+	jets_.genChargedSum[jets_.nref] += ppt;
860	+	if(ppt > hardPtMin_) jets_.genHardSum[jets_.nref] += ppt;
861	+	if(p.collisionId() == 0){
862	+	jets_.signalChargedSum[jets_.nref] += ppt;
863	+	if(ppt > hardPtMin_) jets_.signalHardSum[jets_.nref] += ppt;
864	+	}
865	+
866	+	}
867	+	}
868	+
869	+	}
870	+
871		if(isMC_ && usePat_){
872	<	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
873	<	const reco::GenJet* genjet = (*patjets)[j].genJet();
872	>
873	>
874	>	const reco::GenJet * genjet = (*patjets)[j].genJet();
875
876		if(genjet){
877		jets_.refpt[jets_.nref] = genjet->pt();
#	Line 411 | Line 880 | HiInclusiveJetAnalyzer::analyze(const Ev
880		jets_.refy[jets_.nref] = genjet->eta();
881		jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
882		jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
883	+
884	+	if(doSubEvent_){
885	+	const GenParticle* gencon = genjet->getGenConstituent(0);
886	+	jets_.subid[jets_.nref] = gencon->collisionId();
887	+	}
888	+
889		}else{
890		jets_.refpt[jets_.nref] = -999.;
891		jets_.refeta[jets_.nref] = -999.;
#	Line 420 | Line 895 | HiInclusiveJetAnalyzer::analyze(const Ev
895		jets_.refdrjt[jets_.nref] = -999.;
896		}
897
898	<	// matched partons
899	<	const reco::GenParticle * parton = (*patjets)[j].genParton();
900	<	if(parton){
901	<	jets_.refparton_pt[jets_.nref] = parton->pt();
902	<	jets_.refparton_flavor[jets_.nref] = parton->pdgId();
898	>	jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
899	>
900	>	// matched partons
901	>	const reco::GenParticle & parton = *(*patjets)[j].genParton();
902	>
903	>	if((*patjets)[j].genParton()){
904	>	jets_.refparton_pt[jets_.nref] = parton.pt();
905	>	jets_.refparton_flavor[jets_.nref] = parton.pdgId();
906	>
907	>	if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
908	>
909	>	usedStringPts.clear();
910	>
911	>	// uncomment this if you want to know the ugly truth about parton matching -matt
912	>	//if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
913	>	// cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
914	>
915	>	jets_.bJetIndex[jets_.bMult] = jets_.nref;
916	>	jets_.bStatus[jets_.bMult] = parton.status();
917	>	jets_.bVx[jets_.bMult] = parton.vx();
918	>	jets_.bVy[jets_.bMult] = parton.vy();
919	>	jets_.bVz[jets_.bMult] = parton.vz();
920	>	jets_.bPt[jets_.bMult] = parton.pt();
921	>	jets_.bEta[jets_.bMult] = parton.eta();
922	>	jets_.bPhi[jets_.bMult] = parton.phi();
923	>	jets_.bPdg[jets_.bMult] = parton.pdgId();
924	>	jets_.bChg[jets_.bMult] = parton.charge();
925	>	jets_.bMult++;
926	>	saveDaughters(parton);
927	>	}
928	>
929	>
930		} else {
931	<	jets_.refparton_pt[jets_.nref] = -999;
932	<	jets_.refparton_flavor[jets_.nref] = -999;
931	>	jets_.refparton_pt[jets_.nref] = -999;
932	>	jets_.refparton_flavor[jets_.nref] = -999;
933		}
934	+
935	+
936		}
937
938		jets_.nref++;
#	Line 439 | Line 943 | HiInclusiveJetAnalyzer::analyze(const Ev
943
944		if(isMC_){
945
946	+	edm::Handle<HepMCProduct> hepMCProduct;
947	+	iEvent.getByLabel(eventInfoTag_,hepMCProduct);
948	+	const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
949	+
950	+	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
951	+	if(beamParticles.first != 0)jets_.beamId1 = beamParticles.first->pdg_id();
952	+	if(beamParticles.second != 0)jets_.beamId2 = beamParticles.second->pdg_id();
953	+
954		edm::Handle<GenEventInfoProduct> hEventInfo;
955		iEvent.getByLabel(eventInfoTag_,hEventInfo);
956		//jets_.pthat = hEventInfo->binningValues()[0];
#	Line 458 | Line 970 | HiInclusiveJetAnalyzer::analyze(const Ev
970		// threshold to reduce size of output in minbias PbPb
971		if(genjet_pt>genPtMin_){
972
973	+
974		jets_.genpt[jets_.ngen] = genjet_pt;
975		jets_.geneta[jets_.ngen] = genjet.eta();
976		jets_.genphi[jets_.ngen] = genjet.phi();
977		jets_.geny[jets_.ngen] = genjet.eta();
978
979	+	if(doSubEvent_){
980	+	const GenParticle* gencon = genjet.getGenConstituent(0);
981	+	jets_.gensubid[jets_.ngen] = gencon->collisionId();
982	+	}
983
984		// find matching patJet if there is one
985
986		jets_.gendrjt[jets_.ngen] = -1.0;
987		jets_.genmatchindex[jets_.ngen] = -1;
988
989	<	for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
473	<	const pat::Jet& jet = (*jets)[ijet];
474	<	const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
475	<	if(matchedGenJet){
989	>	for(int ijet = 0 ; ijet < jets_.nref; ++ijet){
990		// poor man's matching, someone fix please
991	<	if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
992	<	fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
993	<	(fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
994	<
995	<	jets_.genmatchindex[jets_.ngen] = (int)ijet;
996	<	jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
483	<	jets_.gendrjt[jets_.ngen] = reco::deltaR(jet,genjet);
991	>	if(fabs(genjet.pt()-jets_.refpt[ijet])<0.00001 &&
992	>	fabs(genjet.eta()-jets_.refeta[ijet])<0.00001){
993	>
994	>	jets_.genmatchindex[jets_.ngen] = (int)ijet;
995	>	jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jets_.refphi[ijet],genjet.phi());
996	>	jets_.gendrjt[jets_.ngen] = sqrt(pow(jets_.gendphijt[jets_.ngen],2)+pow(fabs(genjet.eta()-jets_.refeta[ijet]),2));
997
998		break;
999	<	}
999	>	}
1000		}
1001		}
1002	+
1003		jets_.ngen++;
490	–	}
491	–
1004		}
1005	<
1005	>
1006		}
1007	+
1008	+
1009	+
1010	+
1011	+
1012		t->Fill();
1013	+	memset(&jets_,0,sizeof jets_);
1014	+
1015		}
1016
1017
#	Line 643 | Line 1162 | HiInclusiveJetAnalyzer::getPtRel(const r
1162
1163		return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
1164		}
1165	+
1166	+	// Recursive function, but this version gets called only the first time
1167	+	void
1168	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
1169	+
1170	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
1171	+	const reco::Candidate & daughter = *gen.daughter(i);
1172	+	double daughterPt = daughter.pt();
1173	+	if(daughterPt<1.) continue;
1174	+	double daughterEta = daughter.eta();
1175	+	if(fabs(daughterEta)>3.) continue;
1176	+	int daughterPdgId = daughter.pdgId();
1177	+	int daughterStatus = daughter.status();
1178	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
1179	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
1180	+
1181	+	// cheesy way of finding strings which were already used
1182	+	if(daughter.pdgId()==92){
1183	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
1184	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
1185	+	}
1186	+	usedStringPts.push_back(daughter.pt());
1187	+	}
1188	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
1189	+	jets_.bStatus[jets_.bMult] = daughterStatus;
1190	+	jets_.bVx[jets_.bMult] = daughter.vx();
1191	+	jets_.bVy[jets_.bMult] = daughter.vy();
1192	+	jets_.bVz[jets_.bMult] = daughter.vz();
1193	+	jets_.bPt[jets_.bMult] = daughterPt;
1194	+	jets_.bEta[jets_.bMult] = daughterEta;
1195	+	jets_.bPhi[jets_.bMult] = daughter.phi();
1196	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
1197	+	jets_.bChg[jets_.bMult] = daughter.charge();
1198	+	jets_.bMult++;
1199	+	saveDaughters(daughter);
1200	+	}
1201	+	}
1202	+
1203	+	// This version called for all subsequent calls
1204	+	void
1205	+	HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
1206	+
1207	+	for(unsigned i=0;i<gen.numberOfDaughters();i++){
1208	+	const reco::Candidate & daughter = *gen.daughter(i);
1209	+	double daughterPt = daughter.pt();
1210	+	if(daughterPt<1.) continue;
1211	+	double daughterEta = daughter.eta();
1212	+	if(fabs(daughterEta)>3.) continue;
1213	+	int daughterPdgId = daughter.pdgId();
1214	+	int daughterStatus = daughter.status();
1215	+	// Special case when b->b+string, both b and string contain all daughters, so only take the string
1216	+	if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
1217	+
1218	+	// cheesy way of finding strings which were already used
1219	+	if(daughter.pdgId()==92){
1220	+	for(unsigned ist=0;ist<usedStringPts.size();ist++){
1221	+	if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
1222	+	}
1223	+	usedStringPts.push_back(daughter.pt());
1224	+	}
1225	+
1226	+	jets_.bJetIndex[jets_.bMult] = jets_.nref;
1227	+	jets_.bStatus[jets_.bMult] = daughterStatus;
1228	+	jets_.bVx[jets_.bMult] = daughter.vx();
1229	+	jets_.bVy[jets_.bMult] = daughter.vy();
1230	+	jets_.bVz[jets_.bMult] = daughter.vz();
1231	+	jets_.bPt[jets_.bMult] = daughterPt;
1232	+	jets_.bEta[jets_.bMult] = daughterEta;
1233	+	jets_.bPhi[jets_.bMult] = daughter.phi();
1234	+	jets_.bPdg[jets_.bMult] = daughterPdgId;
1235	+	jets_.bChg[jets_.bMult] = daughter.charge();
1236	+	jets_.bMult++;
1237	+	saveDaughters(daughter);
1238	+	}
1239	+	}
1240	+
1241	+	double HiInclusiveJetAnalyzer::getEt(math::XYZPoint pos, double energy){
1242	+	double et = energy*sin(pos.theta());
1243	+	return et;
1244	+	}
1245	+
1246	+	math::XYZPoint HiInclusiveJetAnalyzer::getPosition(const DetId &id, reco::Vertex::Point vtx){
1247	+	const GlobalPoint& pos=geo->getPosition(id);
1248	+	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
1249	+	return posV;
1250	+	}
1251	+
1252	+
1253
1254		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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