ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/JetFitAnalyzer/test/JetFinderAnalyzer.cc

Comparing UserCode/JetFitAnalyzer/test/JetFinderAnalyzer.cc (file contents):
Revision 1.1 by dnisson, Wed Sep 2 21:48:13 2009 UTC vs.
Revision 1.9 by dnisson, Fri Sep 4 20:44:29 2009 UTC

#	Line 2 | Line 2
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
5	–	#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
5		#include "fastjet/ClusterSequence.hh"
6		#include "FWCore/ServiceRegistry/interface/Service.h"
7	+	#include "FWCore/MessageLogger/interface/MessageLogger.h"
8		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
9
10	+	#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
11	+	#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
12	+	#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
13	+	#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
14	+	#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
15	+
16		#include <map>
17		#include <vector>
18		#include <limits>
#	Line 103 | Line 109 | private:
109		TH2 *);
110		virtual void analyze_results(jetfit::results r,
111		std::vector<jetfit::trouble> t, TH2 *);
112	+	vector<reco::Candidate *> get_particles(const edm::Event&);
113	+	void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
114	+	const edm::InputTag&, const edm::Event&) const;
115	+	void fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>&,
116	+	const edm::InputTag&, const edm::Event&) const;
117
118		fstream ofs;
119	+	edm::InputTag inputTagPFCandidates_;
120	+	edm::InputTag inputTagGenParticles_;
121	+	int info_type_;
122		double smear_;
123		int smear_coord_;
124		};
#	Line 115 | Line 129 | JetFinderAnalyzer::unique_jets;
129		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
130		: JetFitAnalyzer(pSet) // this is important!
131		{
132	+	info_type_ = pSet.getUntrackedParameter("info_type", 0);
133	+
134	+	if (info_type_ == 0) {
135	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
136	+	}
137	+	if (info_type_ == 1) {
138	+	inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
139	+	}
140	+
141		smear_ = pSet.getUntrackedParameter("smear", 0.02);
142		smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
143		// 0 = eta-phi smear
#	Line 122 | Line 145 | JetFinderAnalyzer::JetFinderAnalyzer(con
145		set_user_minuit(jetfinder);
146		}
147
148	+	void
149	+	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c,
150	+	const edm::InputTag& tag,
151	+	const edm::Event& iEvent) const {
152	+
153	+	bool found = iEvent.getByLabel(tag, c);
154	+
155	+	if(!found ) {
156	+	ostringstream  err;
157	+	err<<" cannot get PFCandidates: "
158	+	<<tag<<endl;
159	+	edm::LogError("PFCandidates")<<err.str();
160	+	throw cms::Exception( "MissingProduct", err.str());
161	+	}
162	+
163	+	}
164	+
165	+	void
166	+	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>& c,
167	+	const edm::InputTag& tag,
168	+	const edm::Event& iEvent) const {
169	+
170	+	bool found = iEvent.getByLabel(tag, c);
171	+
172	+	if(!found ) {
173	+	ostringstream  err;
174	+	err<<" cannot get GenParticles: "
175	+	<<tag<<endl;
176	+	edm::LogError("GenParticles")<<err.str();
177	+	throw cms::Exception( "MissingProduct", err.str());
178	+	}
179	+
180	+	}
181	+
182	+	vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
183	+	// fill unreduced histo
184	+	edm::Handle<reco::GenParticleCollection> hRaw;
185	+	edm::Handle<reco::PFCandidateCollection> hPFlow;
186	+	if (info_type_ == 0) {
187	+	fetchCandidateCollection(hRaw,
188	+	inputTagGenParticles_,
189	+	evt);
190	+	}
191	+	if (info_type_ == 1) {
192	+	fetchCandidateCollection(hPFlow,
193	+	inputTagPFCandidates_,
194	+	evt);
195	+	}
196	+
197	+	vector<reco::Candidate *> particles;
198	+
199	+	switch (info_type_) {
200	+	case 0:
201	+	for (unsigned i = 0; i < hRaw->size(); i++) {
202	+	if ((*hRaw)[i].status() == 1) {
203	+	particles.push_back((reco::Candidate *)&((*hRaw)[i]));
204	+	}
205	+	}
206	+	break;
207	+	case 1:
208	+	for (unsigned i = 0; i < hPFlow->size(); i++) {
209	+	particles.push_back((reco::Candidate *)&((*hPFlow)[i]));
210	+	}
211	+	break;
212	+	default:
213	+	cerr << "Unknown event type" << endl; // TODO use MessageLogger
214	+	}
215	+
216	+	return particles;
217	+	}
218	+
219		TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
220		ostringstream oss;
221		oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
222		TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
223		600, -2.5, 2.5, 600, -PI, PI);
224
225	<	// fill unreduced histo
132	<	edm::Handle<edm::HepMCProduct> hRaw;
133	<	evt.getByLabel("source",
134	<	hRaw);
135	<	vector<HepMC::GenParticle *> particles;
136	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
137	<	for (HepMC::GenEvent::particle_const_iterator
138	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
139	<	pit++) {
140	<	if ((*pit)->status() == 1) {
141	<	particles.push_back(*pit);
142	<	}
143	<	}
144	<
225	>	vector<reco::Candidate *> particles = get_particles(evt);
226		for (int i = 0; i < particles.size(); i++) {
227	<	unred_histo->Fill(particles[i]->momentum().eta(),
228	<	particles[i]->momentum().phi(),
229	<	particles[i]->momentum().e());
227	>	unred_histo->Fill(particles[i]->eta(),
228	>	particles[i]->phi(),
229	>	particles[i]->energy());
230		}
231
232		// reduce histo
#	Line 187 | Line 268 | TH2D * JetFinderAnalyzer::make_histo(con
268		60, max_x-rcone, max_x+rcone,
269		60, max_y-rcone, max_y+rcone);
270		for (int i = 0; i < particles.size(); i++) {
271	<	double N = particles[i]->momentum().e();
272	<	double x = particles[i]->momentum().eta();
273	<	double y = particles[i]->momentum().phi();
271	>	double N = particles[i]->energy();
272	>	double x = particles[i]->eta();
273	>	double y = particles[i]->phi();
274		histo_unsmeared->Fill(x, y, N);
275		}
276
#	Line 205 | Line 286 | TH2D * JetFinderAnalyzer::make_histo(con
286		switch (smear_coord_) {
287		case 1:
288		for (int i = 0; i < particles.size(); i++) {
289	<	double N = particles[i]->momentum().e();
290	<	double x = particles[i]->momentum().eta();
291	<	double y = particles[i]->momentum().phi();
289	>	double N = particles[i]->energy();
290	>	double x = particles[i]->eta();
291	>	double y = particles[i]->phi();
292		// loop over bins and add Gaussian in proper angle to smeared
293		for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
294		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
295		double eta = static_cast<double>((signed int)i2) * XbinSize +
296		max_x - rcone - x;
297	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
298	<	max_y - rcone - y));
299	<	phi = sin(phi) > 0 ? phi : -phi;
297	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
298	>	max_y - rcone - y;
299	>	double phi = acos(cos(phi0));
300	>	phi = sin(phi0) > 0 ? phi : -phi;
301
302		// transform eta, phi to proper angle
303		double theta = 2.0*atan(exp(-eta));
#	Line 223 | Line 305 | TH2D * JetFinderAnalyzer::make_histo(con
305
306		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
307		* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
308	+
309	+	// correct for out-of-range phi
310	+	double transform = 0.0;
311	+	if (histo->GetYaxis()->GetXmin() < -PI) {
312	+	transform = -2.0*PI;
313	+	phi += transform;
314	+	iota = asin(sin(theta)*sin(phi));
315	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
316	+	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
317	+	}
318	+	if (histo->GetYaxis()->GetXmax() > PI) {
319	+	transform = 2.0*PI;
320	+	phi += transform;
321	+	iota = asin(sin(theta)*sin(phi));
322	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
323	+	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
324	+	}
325		}
326		}
327		}
#	Line 230 | Line 329 | TH2D * JetFinderAnalyzer::make_histo(con
329		case 0:
330		default:
331		for (int i = 0; i < particles.size(); i++) {
332	<	double N = particles[i]->momentum().e();
333	<	double x = particles[i]->momentum().eta();
334	<	double y = particles[i]->momentum().phi();
332	>	double N = particles[i]->energy();
333	>	double x = particles[i]->eta();
334	>	double y = particles[i]->phi();
335		// loop over bins and add Gaussian to smeared
336		for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
337		for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
338		double eta = static_cast<double>((signed int)i2) * XbinSize
339		+ max_x - rcone - x;
340	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
341	<	+ max_y - rcone - y));
342	<	phi = sin(phi) > 0 ? phi : -phi;
340	>	double phi0 = static_cast<double>((signed int)j2) * YbinSize +
341	>	max_y - rcone - y;
342	>	double phi = acos(cos(phi0));
343	>	phi = sin(phi0) > 0 ? phi : -phi;
344		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
345		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
346	+
347	+	// correct for out-of-range phi
348	+	double transform = 0.0;
349	+	if (histo->GetYaxis()->GetXmin() < -PI) {
350	+	transform = -2.0*PI;
351	+	phi += transform;
352	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
353	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
354	+	}
355	+	if (histo->GetYaxis()->GetXmax() > PI) {
356	+	transform = 2.0*PI;
357	+	phi += transform;
358	+	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
359	+	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
360	+	}
361		}
362		}
363		}
#	Line 257 | Line 372 | TH2D * JetFinderAnalyzer::make_histo(con
372		return histo;
373		}
374
375	<	void seed_with_CA(const edm::Event &evt, TH2 *histo,
375	>	void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo,
376		jetfit::model_def &_mdef) {
262	–	edm::Handle<edm::HepMCProduct> hMC;
263	–	evt.getByLabel("source", hMC);
264	–
377		// create a PseudoJet vector
378		vector<PseudoJet> particles;
379	<	const HepMC::GenEvent *hmcEvt = hMC->GetEvent();
380	<	for (HepMC::GenEvent::particle_const_iterator
381	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
382	<	pit++) {
383	<	if ((*pit)->status() == 1) {
384	<	double x_max = (histo->GetXaxis()->GetXmax()
385	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
386	<	double y_max = (histo->GetYaxis()->GetXmax()
387	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
388	<	valarray<double> pmom(4);
389	<	pmom[0] = (*pit)->momentum().px();
390	<	pmom[1] = (*pit)->momentum().py();
391	<	pmom[2] = (*pit)->momentum().pz();
392	<	pmom[3] = (*pit)->momentum().e();
393	<	double eta = (*pit)->momentum().eta();
394	<	double phi = (*pit)->momentum().phi();
395	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
396	<	PseudoJet j(pmom);
285	<	particles.push_back(j);
286	<	}
379	>	for (unsigned i = 0; i < gParticles.size(); i++) {
380	>	double x_max = (histo->GetXaxis()->GetXmax()
381	>	+ histo->GetXaxis()->GetXmin()) / 2.0;
382	>	double y_max = (histo->GetYaxis()->GetXmax()
383	>	+ histo->GetYaxis()->GetXmin()) / 2.0;
384	>	valarray<double> pmom(4);
385	>	pmom[0] = gParticles[i]->px();
386	>	pmom[1] = gParticles[i]->py();
387	>	pmom[2] = gParticles[i]->pz();
388	>	pmom[3] = gParticles[i]->energy();
389	>	double eta = gParticles[i]->eta();
390	>	double phi = gParticles[i]->phi();
391	>	// phi range corrections
392	>	if (phi - y_max > 2.0*PI) phi -= 2.0*PI;
393	>	if (phi - y_max < -2.0*PI) phi += 2.0*PI;
394	>	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
395	>	PseudoJet j(pmom);
396	>	particles.push_back(j);
397		}
398		}
399
#	Line 306 | Line 416 | void seed_with_CA(const edm::Event &evt,
416		int ijset = 0;
417		for (unsigned ij = 0; ij < jets.size(); ij++) {
418		double N = jets[ij].e();
419	<	if (N > 10.0) {
419	>	if (N > 50.0) {
420		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
421		0.0, 1.0e6);
422		_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
#	Line 327 | Line 437 | jetfit::model_def& JetFinderAnalyzer::ma
437		class jf_model_def : public jetfit::model_def {
438		public:
439		virtual double chisquare_error(double E) {
440	<	return 0.97*E + 14.0;
441	<	// study from 08-27-09
440	>	return exp(-(4.2 + 0.11*E));
441	>	// study from 09-04-09
442		}
443		};
444
#	Line 344 | Line 454 | jetfit::model_def& JetFinderAnalyzer::ma
454		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
455		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
456
457	<	seed_with_CA(evt, histo, *_mdef);
457	>	seed_with_CA(get_particles(evt), histo, *_mdef);
458
459		jetfit::set_model_def(_mdef);
460
#	Line 444 | Line 554 | ostream& operator<<(ostream &out, jetfit
554		else
555		out << "Not calculated" << endl;
556		}
557	+	else {
558	+	out << "Error matrix accurate" << endl;
559	+	}
560
561		return out;
562		}
#	Line 452 | Line 565 | void JetFinderAnalyzer::analyze_results(
565		std::vector<jetfit::trouble> t,
566		TH2 *hist_orig) {
567		ofs << "Histogram "<<hist_orig->GetName() << endl;
568	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
568	>	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
569	>	i < unique_jets[hist_orig].size(); i++) {
570	>	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
571		ofs << "For "<<i+1<<" gaussians: " << endl
572	<	<< t.at(i) << endl
573	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
572	>	<< t.at(i) << endl;
573	>	ofs << "chisquare="<<r.chisquare.at(ir)
574	>	<< endl;
575	>	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
576		for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
577		jet _jet = unique_jets[hist_orig][i][j];
578		ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta

Diff Legend

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