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.3 by dnisson, Thu Sep 3 18:02:35 2009 UTC vs.
Revision 1.6 by dnisson, Fri Sep 4 17:57:32 2009 UTC

#	Line 7 | Line 7
7		#include "FWCore/MessageLogger/interface/MessageLogger.h"
8		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
9
10	<	#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
11	<	#include "DataFormats/Candidate/interface/Particle.h"
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"
#	Line 29 | Line 29
29		using namespace std;
30		using namespace fastjet;
31
32	–	// Class to represent a "generic" particle, whether raw or reconstructed
33	–	class GenericParticle {
34	–	public:
35	–	GenericParticle(double __px, double __py, double __pz, double __e,
36	–	double __charge = 0.0)
37	–	: _px(__px), _py(__py), _pz(__pz), _e(__e), _charge(__charge) {
38	–	}
39	–	GenericParticle(const HepMC::GenParticle &genParticle)
40	–	: _charge(0.0) {
41	–	_px = genParticle.momentum().px();
42	–	_py = genParticle.momentum().py();
43	–	_pz = genParticle.momentum().pz();
44	–	_e = genParticle.momentum().e();
45	–	}
46	–	GenericParticle(const reco::PFCandidate &pfCandidate)
47	–	: _charge(0.0) {
48	–	_px = pfCandidate.px();
49	–	_py = pfCandidate.py();
50	–	_pz = pfCandidate.pz();
51	–	_e = pfCandidate.energy();
52	–	}
53	–	double px() {
54	–	return _px;
55	–	}
56	–	double py() {
57	–	return _py;
58	–	}
59	–	double pz() {
60	–	return _pz;
61	–	}
62	–	double e() {
63	–	return _e;
64	–	}
65	–	double charge() {
66	–	return _charge;
67	–	}
68	–	double eta() {
69	–	double theta = acos(_pz/sqrt(_px*_px + _py*_py + _pz*_pz));
70	–	return -log(tan(theta*0.5));
71	–	}
72	–	double phi() {
73	–	double phi0 = acos(_px/sqrt(_px*_px + _py*_py));
74	–	return _py > 0.0 ? phi0 : -phi0;
75	–	}
76	–	double m() {
77	–	return sqrt(_e*_e - _px*_px - _py*_py - _pz*_pz);
78	–	}
79	–	private:
80	–	double _px;
81	–	double _py;
82	–	double _pz;
83	–	double _e;
84	–	double _charge;
85	–	};
86	–
32		class JetFinderAnalyzer : public JetFitAnalyzer {
33		public:
34		struct jet {
#	Line 164 | Line 109 | private:
109		TH2 *);
110		virtual void analyze_results(jetfit::results r,
111		std::vector<jetfit::trouble> t, TH2 *);
112	<	vector<GenericParticle *> get_particles(const edm::Event&);
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_;
#	Line 183 | Line 131 | JetFinderAnalyzer::JetFinderAnalyzer(con
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		}
#	Line 211 | Line 162 | JetFinderAnalyzer::fetchCandidateCollect
162
163		}
164
165	<	vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
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<edm::HepMCProduct> hRaw;
184	>	edm::Handle<reco::GenParticleCollection> hRaw;
185		edm::Handle<reco::PFCandidateCollection> hPFlow;
186		if (info_type_ == 0) {
187	<	evt.getByLabel("source", hRaw);
187	>	fetchCandidateCollection(hRaw,
188	>	inputTagGenParticles_,
189	>	evt);
190		}
191		if (info_type_ == 1) {
192		fetchCandidateCollection(hPFlow,
#	Line 224 | Line 194 | vector<GenericParticle *> JetFinderAnaly
194		evt);
195		}
196
197	<	vector<GenericParticle *> particles;
197	>	vector<reco::Candidate *> particles;
198
199		switch (info_type_) {
200		case 0:
201	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
202	<	for (HepMC::GenEvent::particle_const_iterator
203	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
234	<	pit++) {
235	<	if ((*pit)->status() == 1) {
236	<	particles.push_back(new GenericParticle(**pit));
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		}
239	–
206		break;
207		case 1:
208		for (unsigned i = 0; i < hPFlow->size(); i++) {
209	<	particles.push_back(new GenericParticle((*hPFlow)[i]));
209	>	particles.push_back((reco::Candidate *)&((*hPFlow)[i]));
210		}
211		break;
212		default:
#	Line 256 | Line 222 | TH2D * JetFinderAnalyzer::make_histo(con
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	<	vector<GenericParticle *> particles = get_particles(evt);
225	>	vector<reco::Candidate *> particles = get_particles(evt);
226		for (int i = 0; i < particles.size(); i++) {
227		unred_histo->Fill(particles[i]->eta(),
228		particles[i]->phi(),
229	<	particles[i]->e());
229	>	particles[i]->energy());
230		}
231
232		// reduce histo
#	Line 302 | 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]->e();
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);
#	Line 320 | 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]->e();
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
#	Line 328 | Line 294 | TH2D * JetFinderAnalyzer::make_histo(con
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 338 | 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 345 | 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]->e();
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
#	Line 353 | Line 337 | TH2D * JetFinderAnalyzer::make_histo(con
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 372 | Line 372 | TH2D * JetFinderAnalyzer::make_histo(con
372		return histo;
373		}
374
375	<	void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
375	>	void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo,
376		jetfit::model_def &_mdef) {
377		// create a PseudoJet vector
378		vector<PseudoJet> particles;
#	Line 385 | Line 385 | void seed_with_CA(vector<GenericParticle
385		pmom[0] = gParticles[i]->px();
386		pmom[1] = gParticles[i]->py();
387		pmom[2] = gParticles[i]->pz();
388	<	pmom[3] = gParticles[i]->e();
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);
#	Line 414 | Line 417 | void seed_with_CA(vector<GenericParticle
417		for (unsigned ij = 0; ij < jets.size(); ij++) {
418		double N = jets[ij].e();
419		if (N > 50.0) {
420	+	cout << N << endl;
421		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
422		0.0, 1.0e6);
423		_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
#	Line 434 | Line 438 | jetfit::model_def& JetFinderAnalyzer::ma
438		class jf_model_def : public jetfit::model_def {
439		public:
440		virtual double chisquare_error(double E) {
441	<	return 0.97*E + 14.0;
442	<	// study from 08-27-09
441	>	return E > 2.7 ? 0.55*E - 1.0 : 0.5;
442	>	// study from 09-03-09
443		}
444		};
445

Diff Legend

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