ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Morgan/src/JetAnalyzer.cc

Comparing UserCode/Morgan/src/JetAnalyzer.cc (file contents):
Revision 1.2 by lethuill, Mon Dec 1 15:58:05 2008 UTC vs.
Revision 1.7 by lethuill, Wed Mar 11 12:44:56 2009 UTC

#	Line 20 | Line 20 | JetAnalyzer::JetAnalyzer(const edm::Para
20		{
21		dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
22		jetProducer_ = producersNames.getParameter<edm::InputTag>("jetProducer");
23	<	useMC_ = myConfig.getParameter<bool>("doMC");
23	>	useMC_ = myConfig.getUntrackedParameter<bool>("doJetMC");
24		}
25
26		JetAnalyzer::~JetAnalyzer()
27		{
28		}
29
30	<	bool Rsortrule (std::pair <double,double> p1, std::pair <double,double> p2 ) {
31	<	return p1.second<p2.second;
30	>	bool Rsortrule (std::pair <double,double> p1, std::pair <double,double> p2 )
31	>	{
32	>	return p1.second<p2.second;
33		}
34
35		void JetAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootJets)
36		{
37
38	<	edm::Handle< vector< pat::Jet > > recoJets;
39	<	iEvent.getByLabel(jetProducer_, recoJets);
40	<	if(verbosity_>1) std::cout <<"  Number of jets = "<< recoJets->size()<< "  Label: "<<jetProducer_.label()<<"  Instance: "<<jetProducer_.instance()<<std::endl;
41	<	for (unsigned int j=0; j<recoJets->size(); j++)
42	<	{
43	<	TRootJet localJet(TLorentzVector( (*recoJets)[j].px(),(*recoJets)[j].py(),(*recoJets)[j].pz(),(*recoJets)[j].energy() ) ,
44	<	TVector3( (*recoJets)[j].vx() ,(*recoJets)[j].vy() ,(*recoJets)[j].vz() ),
45	<	abs((*recoJets)[j].pdgId()),
46	<	(*recoJets)[j].charge()
47	<	//      (*recoJets)[j].genJet()
48	<	);
49	<	localJet.setBtag_trackCountingHighEff((*recoJets)[j].bDiscriminator("trackCountingHighEffBJetTags"));
50	<	localJet.setBtag_trackCountingHighPur((*recoJets)[j].bDiscriminator("trackCountingHighPurBJetTags"));
51	<	localJet.setBtag_jetProbability((*recoJets)[j].bDiscriminator("jetProbabilityBJetTags"));
52	<
53	<	reco::TrackRefVector tracks =  (*recoJets)[j].associatedTracks() ;
54	<	Float_t pTrackerTot = 0;
55	<	Float_t ptTrackerTot = 0;
56	<	//  Float_t eTrackerTot = 0;
57	<	Float_t pTrackerCone01 = 0;
57	<	Float_t pTrackerCone02 = 0;
58	<	Float_t pTrackerCone03 = 0;
59	<	Float_t pTrackerCone04 = 0;
60	<	Float_t pTrackerCone05 = 0;
61	<	vector < pair < Float_t , Float_t > > tracksVPair ;
62	<	pair < Float_t , Float_t > tracksPair;
63	<	for(unsigned int iTracks = 0; iTracks< (*recoJets)[j].associatedTracks().size() ; iTracks ++)
38	>	unsigned int nJets=0;
39	>
40	>	// reco::CaloJet or reco::PFJet ?
41	>	std::string jetType = "BASIC";
42	>	if( jetProducer_.label()=="kt4CaloJets"
43	>	|| jetProducer_.label()=="kt6CaloJets"
44	>	|| jetProducer_.label()=="iterativeCone5CaloJets"
45	>	|| jetProducer_.label()=="sisCone5CaloJets"
46	>	|| jetProducer_.label()=="sisCone7CaloJets"
47	>	) jetType="CALO";
48	>
49	>	if( jetProducer_.label()=="kt4PFJets"
50	>	|| jetProducer_.label()=="kt6PFJets"
51	>	|| jetProducer_.label()=="iterativeCone5PFJets"
52	>	|| jetProducer_.label()=="sisCone5PFJets"
53	>	|| jetProducer_.label()=="sisCone7PFJets"
54	>	) jetType="PF";
55	>
56	>	edm::Handle < std::vector <reco::CaloJet> > recoCaloJets;
57	>	if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="CALO" )
58		{
59	<	//      Float_t buzz = tracks[iTracks]->p();
60	<	pTrackerTot += tracks[iTracks]->p();
67	<	ptTrackerTot += tracks[iTracks]->pt();
68	<	tracksPair.first = tracks[iTracks]->p();
69	<	tracksPair.second = localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) ;
70	<	tracksVPair.push_back(tracksPair);
71	<	// eTrackerTot += tracks[iTracks]->energy();
72	<	if(localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) < 0.5)
73	<	{
74	<	pTrackerCone05 += tracks[iTracks]->p();
75	<	if(localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) < 0.4)
76	<	{
77	<	pTrackerCone04 += tracks[iTracks]->p();
78	<	if(localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) < 0.3)
79	<	{
80	<	pTrackerCone03 += tracks[iTracks]->p();
81	<	if(localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) < 0.2)
82	<	{
83	<	pTrackerCone02 += tracks[iTracks]->p();
84	<	if(localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0)) < 0.1)
85	<	{
86	<	pTrackerCone01 += tracks[iTracks]->p();
87	<	}
88	<	}
89	<	}
90	<	}
91	<	}
59	>	iEvent.getByLabel(jetProducer_, recoCaloJets);
60	>	nJets = recoCaloJets->size();
61		}
62	<	sort(tracksVPair.begin(),tracksVPair.end(), Rsortrule);
63	<	Float_t Rmin = 0;
64	<	Float_t pDummy = 0;
65	<	for(std::vector<std::pair< Float_t,Float_t > >::iterator i = tracksVPair.begin(); i !=tracksVPair.end(); i++) {
66	<	pDummy+=(*i).first;
67	<	if (pDummy>0.75*(pTrackerTot)) {
99	<	Rmin = (*i).second;
100	<	break;
62	>
63	>	edm::Handle < std::vector <reco::PFJet> > recoPFJets;
64	>	if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="PF" )
65	>	{
66	>	iEvent.getByLabel(jetProducer_, recoPFJets);
67	>	nJets = recoPFJets->size();
68		}
102	–	}
103	–	if (Rmin<1e-5) {Rmin=0.;}
69
70	<	localJet.setBroadness(Rmin);
71	<	localJet.setBroadnessDeltaR01(pTrackerCone01/pTrackerTot);
72	<	localJet.setBroadnessDeltaR02(pTrackerCone02/pTrackerTot);
73	<	localJet.setBroadnessDeltaR03(pTrackerCone03/pTrackerTot);
74	<	localJet.setBroadnessDeltaR04(pTrackerCone04/pTrackerTot);
75	<	localJet.setBroadnessDeltaR05(pTrackerCone05/pTrackerTot);
76	<	localJet.setChargedPt(ptTrackerTot);
77	<	localJet.setPartonFlavour((*recoJets)[j].partonFlavour());
78	<	localJet.setCorrB((*recoJets)[j].correctionFactor(pat::Jet::bCorrection));
79	<	localJet.setCorrC((*recoJets)[j].correctionFactor(pat::Jet::cCorrection));
80	<	localJet.setCorrUDS((*recoJets)[j].correctionFactor(pat::Jet::udsCorrection));
81	<	localJet.setCorrGlu((*recoJets)[j].correctionFactor(pat::Jet::gCorrection));
82	<	localJet.setN90((*recoJets)[j].n90());
83	<
84	<	// we still need to implemente the matching or the referencing of the MC part of the Jet see the boolean "useMC_"
85	<
86	<	if((*recoJets)[j].isCaloJet()) {
87	<	// do anything related to a CaloJet
88	<	localJet.setJetType(1);
89	<	localJet.setEcalEnergy((*recoJets)[j].emEnergyFraction());
90	<	localJet.setHcalEnergy((*recoJets)[j].energyFractionHadronic());
91	<	localJet.setChargedEnergy(pTrackerTot);
92	<	localJet.setChargedMultiplicity((*recoJets)[j].associatedTracks().size()) ;
93	<	localJet.setNeutralMultiplicity((*recoJets)[j].getCaloConstituents().size()) ;
94	<	}
95	<	if((*recoJets)[j].isPFJet()) {
96	<	// do anything related to a PFJet
97	<	localJet.setJetType(2);
98	<	localJet.setEcalEnergy((*recoJets)[j].chargedEmEnergyFraction() + (*recoJets)[j].neutralEmEnergyFraction());
99	<	localJet.setHcalEnergy((*recoJets)[j].chargedHadronEnergyFraction() + (*recoJets)[j].neutralHadronEnergyFraction());
100	<	localJet.setChargedEnergy((*recoJets)[j].chargedEmEnergy() + (*recoJets)[j].chargedHadronEnergy() + (*recoJets)[j].chargedMuEnergy() );
101	<	localJet.setChargedMultiplicity((int)(*recoJets)[j].chargedMultiplicity()) ;
102	<	localJet.setNeutralMultiplicity((int)(*recoJets)[j].neutralMultiplicity()) ;
103	<	}
104	<	// if it is neither CaloJet nor PFJet we are bad
105	<	new( (*rootJets)[j] ) TRootJet(localJet);
106	<	if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localJet << endl;
107	<	}
70	>	edm::Handle < std::vector <pat::Jet> > patJets;
71	>	if( dataType_=="PAT" || dataType_=="PATAOD" )
72	>	{
73	>	iEvent.getByLabel(jetProducer_, patJets);
74	>	nJets = patJets->size();
75	>	}
76	>
77	>	if(verbosity_>1) std::cout << "   Number of jets = " << nJets << "   Label: " << jetProducer_.label() << "   Instance: " << jetProducer_.instance() << std::endl;
78	>
79	>	for (unsigned int j=0; j<nJets; j++)
80	>	{
81	>	const reco::Jet* jet = 0;
82	>	if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="CALO" ) jet = (const reco::Jet*) ( & ((*recoCaloJets)[j]) );
83	>	if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="PF" ) jet = (const reco::Jet*) ( & ((*recoPFJets)[j]) );
84	>	if( dataType_=="PAT" || dataType_=="PATAOD" )
85	>	{
86	>	jet = (const reco::Jet*) ( & ((*patJets)[j]) );
87	>	if( (*patJets)[j].isCaloJet() ) jetType="CALO";
88	>	if( (*patJets)[j].isPFJet() ) jetType="PF";
89	>	}
90	>
91	>	TRootJet localJet(
92	>	jet->px()
93	>	,jet->py()
94	>	,jet->pz()
95	>	,jet->energy()
96	>	,jet->vx()
97	>	,jet->vy()
98	>	,jet->vz()
99	>	,jet->pdgId()
100	>	,jet->charge()
101	>	);
102	>
103	>	localJet.setNConstituents(jet->nConstituents());
104	>	localJet.setN90(jet->nCarrying(0.9));
105	>	localJet.setN60(jet->nCarrying(0.6));
106	>	localJet.setJetArea(jet->jetArea());
107	>	localJet.setPileupEnergy(jet->pileup());
108	>	localJet.setMaxDistance(jet->maxDistance());
109	>	float totalEnergy = jet->etInAnnulus(0.,999.);
110	>	if(totalEnergy!=0)
111	>	{
112	>	localJet.setDR01EnergyFraction( jet->etInAnnulus(0.,0.1) / totalEnergy );
113	>	localJet.setDR02EnergyFraction( jet->etInAnnulus(0.,0.2) / totalEnergy );
114	>	localJet.setDR03EnergyFraction( jet->etInAnnulus(0.,0.3) / totalEnergy );
115	>	localJet.setDR04EnergyFraction( jet->etInAnnulus(0.,0.4) / totalEnergy );
116	>	localJet.setDR05EnergyFraction( jet->etInAnnulus(0.,0.5) / totalEnergy );
117	>	}
118	>
119	>
120	>	if( dataType_=="RECO" || dataType_=="AOD" )
121	>	{
122	>	// Some specific methods to reco::Jet
123	>	// Do association to genParticle ?
124	>
125	>	if( jetType=="CALO" )
126	>	{
127	>	// Variables from reco::CaloJet
128	>	const reco::CaloJet *caloJet = dynamic_cast<const reco::CaloJet*>(&*jet);
129	>	localJet.setJetType(1);
130	>	localJet.setEcalEnergyFraction(caloJet->emEnergyFraction());
131	>	localJet.setHcalEnergyFraction(caloJet->energyFractionHadronic());
132	>	//std::vector<CaloTowerPtr> getCaloConstituents () const;
133	>	}
134	>
135	>	if( jetType=="PF" )
136	>	{
137	>	// Variables from reco::PFJet
138	>	const reco::PFJet *pfJet = dynamic_cast<const reco::PFJet*>(&*jet);
139	>	localJet.setJetType(2);
140	>	localJet.setEcalEnergyFraction(pfJet->chargedEmEnergyFraction() + pfJet->neutralEmEnergyFraction());
141	>	localJet.setHcalEnergyFraction(pfJet->chargedHadronEnergyFraction() + pfJet->neutralHadronEnergyFraction());
142	>	if(pfJet->energy() != 0) localJet.setChargedEnergyFraction( (pfJet->chargedEmEnergy() + pfJet->chargedHadronEnergy() + pfJet->chargedMuEnergy() ) / pfJet->energy());
143	>	localJet.setChargedMultiplicity((int)pfJet->chargedMultiplicity()) ;
144	>	//std::vector <const reco::PFCandidate*> getPFConstituents () const;
145	>	}
146	>
147	>	}
148	>
149	>
150	>	if( dataType_=="PATAOD" || dataType_=="PAT" )
151	>	{
152	>	// Some specific methods to pat::Jet
153	>	const pat::Jet *patJet = dynamic_cast<const pat::Jet*>(&*jet);
154
155	+	// Variables from pat::Jet (Basic)
156	+	localJet.setBtag_trackCountingHighEff(patJet->bDiscriminator("trackCountingHighEffBJetTags"));
157	+	localJet.setBtag_trackCountingHighPur(patJet->bDiscriminator("trackCountingHighPurBJetTags"));
158	+	localJet.setBtag_jetProbability(patJet->bDiscriminator("jetProbabilityBJetTags"));
159	+
160	+	// Jet energy correction changed in 22x
161	+	// see DataFormats/PatCandidates/interface/JetCorrFactors.h
162	+	// j'ai pas tout compris, mais je choisis le niveau L7 pour tout le monde
163	+	// A verifier. S.P.
164	+	localJet.setBCorrection(patJet->jetCorrFactors().correction(pat::JetCorrFactors::L7b));
165	+	localJet.setCCorrection(patJet->jetCorrFactors().correction(pat::JetCorrFactors::L7c));
166	+	localJet.setUDSCorrection(patJet->jetCorrFactors().correction(pat::JetCorrFactors::L7uds));
167	+	localJet.setGCorrection(patJet->jetCorrFactors().correction(pat::JetCorrFactors::L7g));
168	+
169	+	// Use  associated tracks to calculate charged broadness of the jet
170	+	// FIXME - Check generalTracks collection is present
171	+	reco::TrackRefVector tracks =  patJet->associatedTracks() ;
172	+	Float_t deltaR = 0.;
173	+	Float_t pTrackerTot = 0.;
174	+	Float_t pTrackerCone01 = 0.;
175	+	Float_t pTrackerCone02 = 0.;
176	+	Float_t pTrackerCone03 = 0.;
177	+	Float_t pTrackerCone04 = 0.;
178	+	Float_t pTrackerCone05 = 0.;
179	+	// TODO - Use std::map....
180	+	vector < pair < Float_t , Float_t > > tracksVPair ;
181	+	pair < Float_t , Float_t > tracksPair;
182	+
183	+	for(unsigned int iTracks = 0; iTracks< patJet->associatedTracks().size() ; iTracks++)
184	+	{
185	+	deltaR = localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0));
186	+	pTrackerTot += tracks[iTracks]->p();
187	+	tracksPair.first = tracks[iTracks]->p();
188	+	tracksPair.second = deltaR;
189	+	tracksVPair.push_back(tracksPair);
190	+	if(deltaR < 0.5) pTrackerCone05 += tracks[iTracks]->p();
191	+	if(deltaR < 0.4) pTrackerCone04 += tracks[iTracks]->p();
192	+	if(deltaR < 0.3) pTrackerCone03 += tracks[iTracks]->p();
193	+	if(deltaR < 0.2) pTrackerCone02 += tracks[iTracks]->p();
194	+	if(deltaR < 0.1) pTrackerCone01 += tracks[iTracks]->p();
195	+	}
196	+	sort(tracksVPair.begin(), tracksVPair.end(), Rsortrule);
197	+	Float_t Rmin = 0;
198	+	Float_t pDummy = 0;
199	+	for(std::vector<std::pair< Float_t,Float_t > >::iterator i = tracksVPair.begin(); i!=tracksVPair.end(); i++)
200	+	{
201	+	pDummy+=(*i).first;
202	+	if (pDummy>0.75*(pTrackerTot))
203	+	{
204	+	Rmin = (*i).second;
205	+	break;
206	+	}
207	+	}
208	+	if (Rmin<1e-5) {Rmin=0.;}
209	+	localJet.setChargedBroadness(Rmin);
210	+	if(pTrackerTot !=0)
211	+	{
212	+	localJet.setChargedBroadnessDR01(pTrackerCone01/pTrackerTot);
213	+	localJet.setChargedBroadnessDR02(pTrackerCone02/pTrackerTot);
214	+	localJet.setChargedBroadnessDR03(pTrackerCone03/pTrackerTot);
215	+	localJet.setChargedBroadnessDR04(pTrackerCone04/pTrackerTot);
216	+	localJet.setChargedBroadnessDR05(pTrackerCone05/pTrackerTot);
217	+	}
218	+
219	+
220	+	if ( patJet->isCaloJet() ) {
221	+	// Variables from pat::Jet (CaloSpecific)
222	+	localJet.setJetType(1);
223	+	localJet.setEcalEnergyFraction(patJet->emEnergyFraction());
224	+	localJet.setHcalEnergyFraction(patJet->energyFractionHadronic());
225	+	localJet.setChargedMultiplicity(patJet->associatedTracks().size()) ;
226	+	//std::vector<CaloTowerPtr> getCaloConstituents () const;
227	+	}
228	+
229	+	if(patJet->isPFJet()) {
230	+	// Variables from pat::Jet (PFSpecific)
231	+	localJet.setJetType(2);
232	+	localJet.setEcalEnergyFraction(patJet->chargedEmEnergyFraction() + patJet->neutralEmEnergyFraction());
233	+	localJet.setHcalEnergyFraction(patJet->chargedHadronEnergyFraction() + patJet->neutralHadronEnergyFraction());
234	+	if(patJet->energy() != 0) localJet.setChargedEnergyFraction( (patJet->chargedEmEnergy() + patJet->chargedHadronEnergy() + patJet->chargedMuEnergy() ) / patJet->energy());
235	+	localJet.setChargedMultiplicity((int)patJet->chargedMultiplicity()) ;
236	+	//std::vector <const reco::PFCandidate*> getPFConstituents () const;
237	+	}
238	+
239	+	// Matched genParticle
240	+	if (useMC_)
241	+	{
242	+	// MC truth associator index
243	+	if ((patJet->genParticleRef()).isNonnull()) {
244	+	localJet.setGenParticleIndex((patJet->genParticleRef()).index());
245	+	} else {
246	+	localJet.setGenParticleIndex(-1);
247	+	}
248	+
249	+	// Matched generated jet
250	+	// FIXME - Keep genJet index or UID instead
251	+	// genJet stocked in MCParticles branch
252	+	/*
253	+	if ( patJet->genJet() != NULL )
254	+	{
255	+	localJet.setMomentumMCJet(TLorentzVector(patJet->genJet()->px(),patJet->genJet()->py(),patJet->genJet()->pz(),patJet->genJet()->energy()));
256	+	localJet.setVertexMCJet(TVector3(patJet->genJet()->vx(),patJet->genJet()->vy(),patJet->genJet()->vz() ));
257	+	localJet.setPdgIdMCJet(patJet->genJet()->pdgId());
258	+	}
259	+	*/
260	+
261	+	// check if jet comes from a top
262	+	bool IsTopJet =  false;
263	+	if(patJet->genParton())
264	+	{
265	+	const reco::Candidate* gen = patJet->genParton();
266	+	while(gen->mother())
267	+	{
268	+	if(abs((gen->mother())->pdgId()) == 6)
269	+	{
270	+	IsTopJet =  true;
271	+	break;
272	+	}
273	+	else
274	+	{
275	+	gen = (gen->mother() );
276	+	}
277	+	}
278	+	}
279	+	localJet.setIsTopJet(IsTopJet);
280	+	}
281	+	}
282	+
283	+	new( (*rootJets)[j] ) TRootJet(localJet);
284	+	if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localJet << endl;
285	+	}
286		}

Diff Legend

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