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.13 by yilmaz, Thu May 10 10:47:17 2012 UTC vs.
Revision 1.26 by yilmaz, Mon Jul 9 17:28:40 2012 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"
#	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	+	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	>
69		if(isMC_){
70		genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
71		eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
#	Line 67 | Line 83 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
83
84		pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
85
86	<	if(!isMC_){
86	>	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
87	>	EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
88	>	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
89	>	HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
90	>
91	>	genParticleSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("genParticles",edm::InputTag("hiGenParticles"));
92	>
93	>	if(doTrigger_){
94		L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
95		hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
96
#	Line 136 | Line 159 | HiInclusiveJetAnalyzer::beginJob() {
159		t->Branch("jty",jets_.jty,"jty[nref]/F");
160		t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
161		t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
162	+	t->Branch("jtm",jets_.jtm,"jtm[nref]/F");
163	+
164	+	// jet ID information, jet composition
165	+	t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
166	+
167	+	t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
168	+	t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
169	+	t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
170	+	t->Branch("trackHardSum", jets_.trackHardSum,"trackHardSum[nref]/F");
171	+	t->Branch("trackHardN", jets_.trackHardN,"trackHardN[nref]/I");
172	+
173	+	t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
174	+	t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
175	+	t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
176	+	t->Branch("chargedHardSum", jets_.chargedHardSum,"chargedHardSum[nref]/F");
177	+	t->Branch("chargedHardN", jets_.chargedHardN,"chargedHardN[nref]/I");
178	+
179	+	t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
180	+	t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
181	+	t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
182	+	t->Branch("photonHardSum", jets_.photonHardSum,"photonHardSum[nref]/F");
183	+	t->Branch("photonHardN", jets_.photonHardN,"photonHardN[nref]/I");
184	+
185	+	t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
186	+	t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
187	+	t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
188	+
189	+	t->Branch("hcalSum", jets_.hcalSum,"hcalSum[nref]/F");
190	+	t->Branch("ecalSum", jets_.ecalSum,"ecalSum[nref]/F");
191	+
192	+	t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
193	+	t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
194	+	t->Branch("eN", jets_.eN,"eN[nref]/I");
195	+
196	+	t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
197	+	t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
198	+	t->Branch("muN", jets_.muN,"muN[nref]/I");
199	+
200	+	t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
201	+	t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
202
203		// b-jet discriminators
204		if (doLifeTimeTagging_) {
#	Line 183 | Line 246 | HiInclusiveJetAnalyzer::beginJob() {
246		t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
247		t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
248		}
249	+
250
251		if(isMC_){
252		t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
#	Line 202 | Line 266 | HiInclusiveJetAnalyzer::beginJob() {
266		t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
267		t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
268
269	<	// For all gen jets, matched or unmatched
270	<	t->Branch("ngen",&jets_.ngen,"ngen/I");
271	<	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
272	<	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
209	<	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
210	<	t->Branch("geny",jets_.geny,"geny[ngen]/F");
211	<	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
212	<	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
213	<	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
269	>	t->Branch("genChargedSum", jets_.genChargedSum,"genChargedSum[nref]/F");
270	>	t->Branch("genHardSum", jets_.genHardSum,"genHardSum[nref]/F");
271	>	t->Branch("signalChargedSum", jets_.signalChargedSum,"signalChargedSum[nref]/F");
272	>	t->Branch("signalHardSum", jets_.signalHardSum,"signalHardSum[nref]/F");
273
274		if(doSubEvent_){
275	<	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
275	>	t->Branch("subid",jets_.subid,"subid[nref]/I");
276		}
277
278	+	if(fillGenJets_){
279	+	// For all gen jets, matched or unmatched
280	+	t->Branch("ngen",&jets_.ngen,"ngen/I");
281	+	t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
282	+	t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
283	+	t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
284	+	t->Branch("geny",jets_.geny,"geny[ngen]/F");
285	+	t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
286	+	t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
287	+	t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
288	+
289	+	if(doSubEvent_){
290	+	t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
291	+	}
292	+	}
293	+
294		if(saveBfragments_  ) {
295		t->Branch("bMult",&jets_.bMult,"bMult/I");
296		t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
#	Line 268 | Line 343 | HiInclusiveJetAnalyzer::analyze(const Ev
343		double hf = 0.;
344		double b = 999.;
345
346	+	if(!geo){
347	+	edm::ESHandle<CaloGeometry> pGeo;
348	+	iSetup.get<CaloGeometryRecord>().get(pGeo);
349	+	geo = pGeo.product();
350	+	}
351		if(useCentrality_){
352		if(!centrality_) centrality_ = new CentralityProvider(iSetup);
353		centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
#	Line 280 | Line 360 | HiInclusiveJetAnalyzer::analyze(const Ev
360		b = centrality_->bMean();
361		}
362
283	–
284	–
285	–
286	–
363		// loop the events
364
365		jets_.bin = bin;
366		jets_.hf = hf;
367
368	<
369	<	if (useVtx_) {
370	<	edm::Handle<vector<reco::Vertex> >vertex;
371	<	iEvent.getByLabel(vtxTag_, vertex);
372	<
373	<	if(vertex->size()>0) {
374	<	jets_.vx=vertex->begin()->x();
375	<	jets_.vy=vertex->begin()->y();
376	<	jets_.vz=vertex->begin()->z();
377	<	}
378	<	}
368	>	reco::Vertex::Point vtx(0,0,0);
369	>	if (useVtx_) {
370	>	edm::Handle<vector<reco::Vertex> >vertex;
371	>	iEvent.getByLabel(vtxTag_, vertex);
372	>
373	>	if(vertex->size()>0) {
374	>	jets_.vx=vertex->begin()->x();
375	>	jets_.vy=vertex->begin()->y();
376	>	jets_.vz=vertex->begin()->z();
377	>	vtx = vertex->begin()->position();
378	>	}
379	>	}
380
381		edm::Handle<pat::JetCollection> patjets;
382		if(usePat_)iEvent.getByLabel(jetTag_, patjets);
383	+
384	+	edm::Handle<pat::JetCollection> matchedjets;
385	+	iEvent.getByLabel(matchTag_, matchedjets);
386
387		edm::Handle<reco::JetView> jets;
388		iEvent.getByLabel(jetTag_, jets);
389
390		edm::Handle<reco::PFCandidateCollection> pfCandidates;
391	<	if(doLifeTimeTagging_){
392	<	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
393	<	}
394	<	// FILL JRA TREE
391	>	iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
392	>
393	>	edm::Handle<reco::TrackCollection> tracks;
394	>	iEvent.getByLabel(trackTag_,tracks);
395	>
396	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
397	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
398	>
399	>	edm::Handle<HFRecHitCollection> hfHits;
400	>	edm::Handle<HBHERecHitCollection> hbheHits;
401
402	+	iEvent.getByLabel(HcalRecHitHBHESrc_,hbheHits);
403	+	iEvent.getByLabel(HcalRecHitHFSrc_,hfHits);
404	+	iEvent.getByLabel(EBSrc_,ebHits);
405	+	iEvent.getByLabel(EESrc_,eeHits);
406	+
407	+	edm::Handle<reco::GenParticleCollection> genparts;
408	+	iEvent.getByLabel(genParticleSrc_,genparts);
409	+
410	+
411	+	// FILL JRA TREE
412		jets_.b = b;
413		jets_.nref = 0;
414	+
415
416	<	if(!isMC_){
416	>	if(doTrigger_){
417		fillL1Bits(iEvent);
418		fillHLTBits(iEvent);
419		}
420
421	<	for(unsigned int j = 0 ; j < jets->size(); ++j){
422	<	const reco::Jet& jet = (*jets)[j];
326	<
327	<
421	>	for(unsigned int j = 0; j < jets->size(); ++j){
422	>	const reco::Jet& jet = (*jets)[j];
423		if (useJEC_ && usePat_){
424		jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
425		}
#	Line 416 | Line 511 | HiInclusiveJetAnalyzer::analyze(const Ev
511		jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
512		jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
513		jets_.muchg[jets_.nref]   =  muon.charge();
514	+	}else{
515	+	jets_.mupt[jets_.nref]    =  0.0;
516	+	jets_.mueta[jets_.nref]   =  0.0;
517	+	jets_.muphi[jets_.nref]   =  0.0;
518	+	jets_.mue[jets_.nref]     =  0.0;
519	+	jets_.mudr[jets_.nref]    =  9.9;
520	+	jets_.muptrel[jets_.nref] =  0.0;
521	+	jets_.muchg[jets_.nref]   = 0;
522		}
523	+
524		}
525
526	+	// Jet ID variables
527	+
528	+	jets_.muMax[jets_.nref] = 0;
529	+	jets_.muSum[jets_.nref] = 0;
530	+	jets_.muN[jets_.nref] = 0;
531	+
532	+	jets_.eMax[jets_.nref] = 0;
533	+	jets_.eSum[jets_.nref] = 0;
534	+	jets_.eN[jets_.nref] = 0;
535	+
536	+	jets_.neutralMax[jets_.nref] = 0;
537	+	jets_.neutralSum[jets_.nref] = 0;
538	+	jets_.neutralN[jets_.nref] = 0;
539	+
540	+	jets_.photonMax[jets_.nref] = 0;
541	+	jets_.photonSum[jets_.nref] = 0;
542	+	jets_.photonN[jets_.nref] = 0;
543	+	jets_.photonHardSum[jets_.nref] = 0;
544	+	jets_.photonHardN[jets_.nref] = 0;
545	+
546	+	jets_.chargedMax[jets_.nref] = 0;
547	+	jets_.chargedSum[jets_.nref] = 0;
548	+	jets_.chargedN[jets_.nref] = 0;
549	+	jets_.chargedHardSum[jets_.nref] = 0;
550	+	jets_.chargedHardN[jets_.nref] = 0;
551	+
552	+	jets_.trackMax[jets_.nref] = 0;
553	+	jets_.trackSum[jets_.nref] = 0;
554	+	jets_.trackN[jets_.nref] = 0;
555	+	jets_.trackHardSum[jets_.nref] = 0;
556	+	jets_.trackHardN[jets_.nref] = 0;
557	+
558	+	jets_.hcalSum[jets_.nref] = 0;
559	+	jets_.ecalSum[jets_.nref] = 0;
560	+
561	+	jets_.genChargedSum[jets_.nref] = 0;
562	+	jets_.genHardSum[jets_.nref] = 0;
563	+
564	+	jets_.signalChargedSum[jets_.nref] = 0;
565	+	jets_.signalHardSum[jets_.nref] = 0;
566	+
567	+	jets_.subid[jets_.nref] = -1;
568	+
569	+	for(unsigned int icand = 0; icand < tracks->size(); ++icand){
570	+	const reco::Track& track = (*tracks)[icand];
571	+	if(useQuality_ ){
572	+	bool goodtrack = track.quality(reco::TrackBase::qualityByName(trackQuality_));
573	+	if(!goodtrack) continue;
574	+	}
575	+
576	+	double dr = deltaR(jet,track);
577	+	if(dr < rParam){
578	+	double ptcand = track.pt();
579	+	jets_.trackSum[jets_.nref] += ptcand;
580	+	jets_.trackN[jets_.nref] += 1;
581	+
582	+	if(ptcand > hardPtMin_){
583	+	jets_.trackHardSum[jets_.nref] += ptcand;
584	+	jets_.trackHardN[jets_.nref] += 1;
585	+
586	+	}
587	+	if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
588	+
589	+	}
590	+	}
591	+
592	+	for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
593	+	const reco::PFCandidate& track = (*pfCandidates)[icand];
594	+	double dr = deltaR(jet,track);
595	+	if(dr < rParam){
596	+	double ptcand = track.pt();
597	+	int pfid = track.particleId();
598	+
599	+	switch(pfid){
600	+
601	+	case 1:
602	+	jets_.chargedSum[jets_.nref] += ptcand;
603	+	jets_.chargedN[jets_.nref] += 1;
604	+	if(ptcand > hardPtMin_){
605	+	jets_.chargedHardSum[jets_.nref] += ptcand;
606	+	jets_.chargedHardN[jets_.nref] += 1;
607	+	}
608	+	if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
609	+	break;
610	+
611	+	case 2:
612	+	jets_.eSum[jets_.nref] += ptcand;
613	+	jets_.eN[jets_.nref] += 1;
614	+	if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
615	+	break;
616	+
617	+	case 3:
618	+	jets_.muSum[jets_.nref] += ptcand;
619	+	jets_.muN[jets_.nref] += 1;
620	+	if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
621	+	break;
622	+
623	+	case 4:
624	+	jets_.photonSum[jets_.nref] += ptcand;
625	+	jets_.photonN[jets_.nref] += 1;
626	+	if(ptcand > hardPtMin_){
627	+	jets_.photonHardSum[jets_.nref] += ptcand;
628	+	jets_.photonHardN[jets_.nref] += 1;
629	+	}
630	+	if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
631	+	break;
632	+
633	+	case 5:
634	+	jets_.neutralSum[jets_.nref] += ptcand;
635	+	jets_.neutralN[jets_.nref] += 1;
636	+	if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
637	+	break;
638	+
639	+	default:
640	+	break;
641	+
642	+	}
643	+	}
644	+	}
645	+
646	+	// Calorimeter fractions
647	+
648	+	for(unsigned int i = 0; i < hbheHits->size(); ++i){
649	+	const HBHERecHit & hit= (*hbheHits)[i];
650	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
651	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
652	+	if(dr < rParam){
653	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
654	+	}
655	+	}
656	+
657	+	for(unsigned int i = 0; i < hfHits->size(); ++i){
658	+	const HFRecHit & hit= (*hfHits)[i];
659	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
660	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
661	+	if(dr < rParam){
662	+	jets_.hcalSum[jets_.nref] += getEt(pos,hit.energy());
663	+	}
664	+	}
665	+
666	+
667	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
668	+	const EcalRecHit & hit= (*ebHits)[i];
669	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
670	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
671	+	if(dr < rParam){
672	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
673	+	}
674	+	}
675	+
676	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
677	+	const EcalRecHit & hit= (*eeHits)[i];
678	+	math::XYZPoint pos = getPosition(hit.id(),vtx);
679	+	double dr = deltaR(jet.eta(),jet.phi(),pos.eta(),pos.phi());
680	+	if(dr < rParam){
681	+	jets_.ecalSum[jets_.nref] += getEt(pos,hit.energy());
682	+	}
683	+	}
684	+
685	+	// Alternative reconstruction matching (PF for calo, calo for PF)
686	+
687	+	double drMin = 100;
688	+	for(unsigned int imatch = 0 ; imatch < matchedjets->size(); ++imatch){
689	+	const reco::Jet& mjet = (*matchedjets)[imatch];
690	+
691	+	double dr = deltaR(jet,mjet);
692	+	if(dr < drMin){
693	+	jets_.matchedPt[jets_.nref] = mjet.pt();
694	+	jets_.matchedR[jets_.nref] = dr;
695	+	drMin = dr;
696	+	}
697	+	}
698	+
699	+	//     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
700	+
701	+	/////////////////////////////////////////////////////////////////
702	+	// Jet core pt^2 discriminant for fake jets
703	+	// Edited by Yue Shi Lai <ylai@mit.edu>
704	+
705	+	// Initial value is 0
706	+	jets_.discr_fr01[jets_.nref] = 0;
707	+	// Start with no directional adaption, i.e. the fake rejection
708	+	// axis is the jet axis
709	+	float pseudorapidity_adapt = jets_.jteta[jets_.nref];
710	+	float azimuth_adapt = jets_.jtphi[jets_.nref];
711	+
712	+	// Unadapted discriminant with adaption search
713	+	for (size_t iteration = 0; iteration < 2; iteration++) {
714	+	float pseudorapidity_adapt_new = pseudorapidity_adapt;
715	+	float azimuth_adapt_new = azimuth_adapt;
716	+	float max_weighted_perp = 0;
717	+	float perp_square_sum = 0;
718	+
719	+	for (size_t index_pf_candidate = 0;
720	+	index_pf_candidate < pfCandidates->size();
721	+	index_pf_candidate++) {
722	+	const reco::PFCandidate &p =
723	+	(*pfCandidates)[index_pf_candidate];
724	+
725	+	switch (p.particleId()) {
726	+	//case 1:     // Charged hadron
727	+	//case 3:     // Muon
728	+	case 4: // Photon
729	+	{
730	+	const float dpseudorapidity =
731	+	p.eta() - pseudorapidity_adapt;
732	+	const float dazimuth =
733	+	reco::deltaPhi(p.phi(), azimuth_adapt);
734	+	// The Gaussian scale factor is 0.5 / (0.1 * 0.1)
735	+	// = 50
736	+	const float angular_weight =
737	+	exp(-50.0F * (dpseudorapidity * dpseudorapidity +
738	+	dazimuth * dazimuth));
739	+	const float weighted_perp =
740	+	angular_weight * p.pt() * p.pt();
741	+	const float weighted_perp_square =
742	+	weighted_perp * p.pt();
743	+
744	+	perp_square_sum += weighted_perp_square;
745	+	if (weighted_perp >= max_weighted_perp) {
746	+	pseudorapidity_adapt_new = p.eta();
747	+	azimuth_adapt_new = p.phi();
748	+	max_weighted_perp = weighted_perp;
749	+	}
750	+	}
751	+	default:
752	+	break;
753	+	}
754	+	}
755	+	// Update the fake rejection value
756	+	jets_.discr_fr01[jets_.nref] = std::max(
757	+	jets_.discr_fr01[jets_.nref], perp_square_sum);
758	+	// Update the directional adaption
759	+	pseudorapidity_adapt = pseudorapidity_adapt_new;
760	+	azimuth_adapt = azimuth_adapt_new;
761	+	}
762	+
763		jets_.jtpt[jets_.nref] = jet.pt();
764		jets_.jteta[jets_.nref] = jet.eta();
765		jets_.jtphi[jets_.nref] = jet.phi();
766		jets_.jty[jets_.nref] = jet.eta();
767		jets_.jtpu[jets_.nref] = jet.pileup();
768	+	jets_.jtm[jets_.nref] = jet.mass();
769
770		if(isMC_ && usePat_){
771
772	+	for(UInt_t i = 0; i < genparts->size(); ++i){
773	+	const reco::GenParticle& p = (*genparts)[i];
774	+	if (p.status()!=1) continue;
775	+	if (p.charge()==0) continue;
776	+	double dr = deltaR(jet,p);
777	+	if(dr < rParam){
778	+	double ppt = p.pt();
779	+	jets_.genChargedSum[jets_.nref] += ppt;
780	+	if(ppt > hardPtMin_) jets_.genHardSum[jets_.nref] += ppt;
781	+	if(p.collisionId() == 0){
782	+	jets_.signalChargedSum[jets_.nref] += ppt;
783	+	if(ppt > hardPtMin_) jets_.signalHardSum[jets_.nref] += ppt;
784	+	}
785	+
786	+	}
787	+	}
788	+
789		const reco::GenJet * genjet = (*patjets)[j].genJet();
790
791		if(genjet){
#	Line 436 | Line 795 | HiInclusiveJetAnalyzer::analyze(const Ev
795		jets_.refy[jets_.nref] = genjet->eta();
796		jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
797		jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
798	+
799	+	if(doSubEvent_){
800	+	const GenParticle* gencon = genjet->getGenConstituent(0);
801	+	jets_.subid[jets_.nref] = gencon->collisionId();
802	+	}
803	+
804		}else{
805		jets_.refpt[jets_.nref] = -999.;
806		jets_.refeta[jets_.nref] = -999.;
#	Line 496 | Line 861 | HiInclusiveJetAnalyzer::analyze(const Ev
861		edm::Handle<HepMCProduct> hepMCProduct;
862		iEvent.getByLabel(eventInfoTag_,hepMCProduct);
863		const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
864	<	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
865	<	jets_.beamId1 = beamParticles.first->pdg_id();
866	<	jets_.beamId2 = beamParticles.second->pdg_id();
864	>
865	>	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
866	>	if(beamParticles.first != 0)jets_.beamId1 = beamParticles.first->pdg_id();
867	>	if(beamParticles.second != 0)jets_.beamId2 = beamParticles.second->pdg_id();
868
869		edm::Handle<GenEventInfoProduct> hEventInfo;
870		iEvent.getByLabel(eventInfoTag_,hEventInfo);
#	Line 557 | Line 923 | HiInclusiveJetAnalyzer::analyze(const Ev
923		}
924
925		}
926	+
927		t->Fill();
928		memset(&jets_,0,sizeof jets_);
929	+
930		}
931
932
#	Line 784 | Line 1152 | HiInclusiveJetAnalyzer::saveDaughters(co
1152		saveDaughters(daughter);
1153		}
1154		}
1155	+
1156	+	double HiInclusiveJetAnalyzer::getEt(math::XYZPoint pos, double energy){
1157	+	double et = energy*sin(pos.theta());
1158	+	return et;
1159	+	}
1160	+
1161	+	math::XYZPoint HiInclusiveJetAnalyzer::getPosition(const DetId &id, reco::Vertex::Point vtx){
1162	+	const GlobalPoint& pos=geo->getPosition(id);
1163	+	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
1164	+	return posV;
1165	+	}
1166	+
1167	+
1168
1169		DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);

Diff Legend

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