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.17 by dnisson, Wed Oct 14 02:45:58 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"
15
16	+	/* Jet reco stuff */
17	+	#include "DataFormats/JetReco/interface/PFJetCollection.h"
18	+	#include "DataFormats/JetReco/interface/GenJetCollection.h"
19	+	#include "DataFormats/JetReco/interface/PFJet.h"
20	+	#include "DataFormats/JetReco/interface/GenJet.h"
21	+
22		#include <map>
23		#include <vector>
24		#include <limits>
#	Line 29 | Line 35
35		using namespace std;
36		using namespace fastjet;
37
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	–
38		class JetFinderAnalyzer : public JetFitAnalyzer {
39		public:
40		struct jet {
#	Line 119 | Line 70 | private:
70		return sum1 * sum2;
71		}
72
73	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
74	<	double dist_sq = numeric_limits<double>::infinity();
124	<	unique_jets[hist].resize(ngauss);
125	<	int nbinsX = hist->GetXaxis()->GetNbins();
126	<	int nbinsY = hist->GetYaxis()->GetNbins();
127	<	double XbinSize = (hist->GetXaxis()->GetXmax()
128	<	- hist->GetXaxis()->GetXmin())
129	<	/ static_cast<double>(nbinsX);
130	<	double YbinSize = (hist->GetYaxis()->GetXmax()
131	<	- hist->GetYaxis()->GetXmin())
132	<	/ static_cast<double>(nbinsY);
133	<	for (int i = 0; i < ngauss; i++) {
134	<	double N, mu_x, mu_y, sig, err, lo, hi;
135	<	int iuint;
136	<	TString name;
137	<	gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
138	<	gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
139	<	gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
140	<	gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
141	<	for (int j = 0; j < i; j++) {
142	<	double N2, mu_x2, mu_y2, sig2;
143	<	gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
144	<	gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
145	<	gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
146	<	gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
147	<	double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
148	<	+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
149	<	if (_dist_sq < dist_sq)
150	<	dist_sq = _dist_sq;
151	<	}
152	<
153	<	jet j;
154	<	j.energy = N;
155	<	j.eta = mu_x; j.phi = mu_y;
156	<	unique_jets[hist][ngauss-1].push_back(j);
157	<	}
158	<	}
159	<
160	<	virtual void beginJob(const edm::EventSetup&);
161	<	virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
73	>	virtual void beginJob(const edm::EventSetup &es);
74	>	virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
75		virtual jetfit::model_def& make_model_def(const edm::Event&,
76		const edm::EventSetup&,
77		TH2 *);
78		virtual void analyze_results(jetfit::results r,
79		std::vector<jetfit::trouble> t, TH2 *);
80	<	vector<GenericParticle *> get_particles(const edm::Event&);
168	<	void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
169	<	const edm::InputTag&, const edm::Event&) const;
80	>	vector<reco::Candidate *> get_particles(const edm::Event&);
81
82		fstream ofs;
83		edm::InputTag inputTagPFCandidates_;
84	+	edm::InputTag inputTagGenParticles_;
85		int info_type_;
86		double smear_;
87		int smear_coord_;
88	+	string jet_algo_;
89	+	string reclustered_jet_algo_;
90		};
91
92		map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
#	Line 183 | Line 97 | JetFinderAnalyzer::JetFinderAnalyzer(con
97		{
98		info_type_ = pSet.getUntrackedParameter("info_type", 0);
99
100	+	if (info_type_ == 0) {
101	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
102	+	}
103		if (info_type_ == 1) {
104		inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
105		}
#	Line 191 | Line 108 | JetFinderAnalyzer::JetFinderAnalyzer(con
108		smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
109		// 0 = eta-phi smear
110		// 1 = proper angle smear
111	<	set_user_minuit(jetfinder);
111	>	jet_algo_ = pSet.getParameter<string>("jet_algo");
112	>	set_user_minuit(0); // we don't need this feature--it's deprecated
113		}
114
115	<	void
116	<	JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c,
199	<	const edm::InputTag& tag,
200	<	const edm::Event& iEvent) const {
201	<
202	<	bool found = iEvent.getByLabel(tag, c);
203	<
204	<	if(!found ) {
205	<	ostringstream  err;
206	<	err<<" cannot get PFCandidates: "
207	<	<<tag<<endl;
208	<	edm::LogError("PFCandidates")<<err.str();
209	<	throw cms::Exception( "MissingProduct", err.str());
210	<	}
211	<
212	<	}
115	>	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
116	>	const reco::Jet * highest_e_jet = 0;
117
214	–	vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
215	–	// fill unreduced histo
216	–	edm::Handle<edm::HepMCProduct> hRaw;
217	–	edm::Handle<reco::PFCandidateCollection> hPFlow;
118		if (info_type_ == 0) {
119	<	evt.getByLabel("source", hRaw);
120	<	}
121	<	if (info_type_ == 1) {
122	<	fetchCandidateCollection(hPFlow,
123	<	inputTagPFCandidates_,
124	<	evt);
119	>	edm::Handle< vector<reco::GenJet> > evtJets;
120	>	evt.getByLabel(jet_algo_, evtJets);
121	>	for (unsigned i = 0; i < evtJets->size(); i++) {
122	>	if (highest_e_jet == 0
123	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
124	>	highest_e_jet = &((*evtJets)[i]);
125	>	}
126	>	}
127		}
128	<
129	<	vector<GenericParticle *> particles;
130	<
131	<	switch (info_type_) {
132	<	case 0:
133	<	const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
134	<	for (HepMC::GenEvent::particle_const_iterator
233	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
234	<	pit++) {
235	<	if ((*pit)->status() == 1) {
236	<	particles.push_back(new GenericParticle(**pit));
128	>	else if (info_type_ == 1) {
129	>	edm::Handle< vector<reco::PFJet> > evtJets;
130	>	evt.getByLabel(jet_algo_, evtJets);
131	>	for (unsigned i = 0; i < evtJets->size(); i++) {
132	>	if (highest_e_jet == 0
133	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
134	>	highest_e_jet = &((*evtJets)[i]);
135		}
136		}
137	+	}
138
139	<	break;
140	<	case 1:
141	<	for (unsigned i = 0; i < hPFlow->size(); i++) {
243	<	particles.push_back(new GenericParticle((*hPFlow)[i]));
244	<	}
245	<	break;
246	<	default:
247	<	cerr << "Unknown event type" << endl; // TODO use MessageLogger
139	>	if (highest_e_jet == 0) {
140	>	cerr << "No fat jets found!" << endl;
141	>	return 0;
142		}
143
144	<	return particles;
145	<	}
146	<
253	<	TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
144	>	vector<const reco::Candidate *> particles =
145	>	highest_e_jet->getJetConstituentsQuick();
146	>
147		ostringstream oss;
148	<	oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
256	<	TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
257	<	600, -2.5, 2.5, 600, -PI, PI);
258	<
259	<	vector<GenericParticle *> particles = get_particles(evt);
260	<	for (int i = 0; i < particles.size(); i++) {
261	<	unred_histo->Fill(particles[i]->eta(),
262	<	particles[i]->phi(),
263	<	particles[i]->e());
264	<	}
265	<
266	<	// reduce histo
267	<	ostringstream oss2;
268	<	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
148	>	oss << "eta_phi_energy"<<evt.id().event() << flush;
149		edm::Service<TFileService> fs;
150	<	// draw cone of radius 0.5 around highest energy bin, reduce
151	<	double maxE = 0.0;
152	<	int max_i = 29, max_j = 29;
153	<	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
154	<	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
155	<	double E = unred_histo->GetBinContent(i+1, j+1);
156	<	if (E > maxE) {
277	<	maxE = E;
278	<	max_i = i;
279	<	max_j = j;
280	<	}
281	<	}
282	<	}
283	<
284	<	double rcone = 0.5;
285	<	double Xlo = unred_histo->GetXaxis()->GetXmin();
286	<	double Xhi = unred_histo->GetXaxis()->GetXmax();
287	<	double Ylo = unred_histo->GetYaxis()->GetXmin();
288	<	double Yhi = unred_histo->GetYaxis()->GetXmax();
289	<	double XbinSize = (Xhi - Xlo) /
290	<	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
291	<	double YbinSize = (Yhi - Ylo) /
292	<	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
293	<	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
294	<	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
295	<	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
296	<	60, max_x-rcone, max_x+rcone,
297	<	60, max_y-rcone, max_y+rcone);
298	<
299	<	// create an unsmeared reduced histo
300	<	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
301	<	(oss2.str()+"_unsmeared").c_str(),
302	<	60, max_x-rcone, max_x+rcone,
303	<	60, max_y-rcone, max_y+rcone);
304	<	for (int i = 0; i < particles.size(); i++) {
305	<	double N = particles[i]->e();
306	<	double x = particles[i]->eta();
307	<	double y = particles[i]->phi();
308	<	histo_unsmeared->Fill(x, y, N);
309	<	}
150	>	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
151	>	30,
152	>	highest_e_jet->eta()-0.5,
153	>	highest_e_jet->eta()+0.5,
154	>	30,
155	>	highest_e_jet->phi()-0.5,
156	>	highest_e_jet->phi()+0.5);
157
158	<	// create a smeared reduced histo
158	>	if (smear_ > 0.0) {
159	>	// create a smeared histo
160		// create a temporary 2D vector for smeared energies
161	<	XbinSize = (histo->GetXaxis()->GetXmax()
161	>	double XbinSize = (histo->GetXaxis()->GetXmax()
162		- histo->GetXaxis()->GetXmin()) /
163		static_cast<double>(histo->GetXaxis()->GetNbins());
164	<	YbinSize = (histo->GetYaxis()->GetXmax()
164	>	double YbinSize = (histo->GetYaxis()->GetXmax()
165		- histo->GetYaxis()->GetXmin()) /
166		static_cast<double>(histo->GetYaxis()->GetNbins());
167	<	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
167	>	double Xlo = histo->GetXaxis()->GetXmin();
168	>	double Xhi = histo->GetXaxis()->GetXmax();
169	>	double Ylo = histo->GetYaxis()->GetXmin();
170	>	double Yhi = histo->GetYaxis()->GetXmax();
171	>	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
172		switch (smear_coord_) {
173		case 1:
174		for (int i = 0; i < particles.size(); i++) {
175	<	double N = particles[i]->e();
175	>	double N = particles[i]->energy();
176		double x = particles[i]->eta();
177		double y = particles[i]->phi();
178	+	if (y < Ylo) {
179	+	y += 2.0*PI;
180	+	}
181	+	if (y > Yhi) {
182	+	y -= 2.0*PI;
183	+	}
184		// loop over bins and add Gaussian in proper angle to smeared
185	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
186	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
185	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
186	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
187		double eta = static_cast<double>((signed int)i2) * XbinSize +
188	<	max_x - rcone - x;
189	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
190	<	max_y - rcone - y));
191	<	phi = sin(phi) > 0 ? phi : -phi;
334	<
188	>	Xlo - x;
189	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
190	>	Ylo - y;
191	>
192		// transform eta, phi to proper angle
193		double theta = 2.0*atan(exp(-eta));
194		double iota = asin(sin(theta)*sin(phi));
#	Line 345 | Line 202 | TH2D * JetFinderAnalyzer::make_histo(con
202		case 0:
203		default:
204		for (int i = 0; i < particles.size(); i++) {
205	<	double N = particles[i]->e();
205	>	double N = particles[i]->energy();
206		double x = particles[i]->eta();
207		double y = particles[i]->phi();
208	+	if (y < Ylo) {
209	+	y += 2.0*PI;
210	+	}
211	+	if (y > Yhi) {
212	+	y -= 2.0*PI;
213	+	}
214		// loop over bins and add Gaussian to smeared
215	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
216	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
215	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
216	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
217		double eta = static_cast<double>((signed int)i2) * XbinSize
218	<	+ max_x - rcone - x;
219	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
220	<	+ max_y - rcone - y));
221	<	phi = sin(phi) > 0 ? phi : -phi;
218	>	+ Xlo - x;
219	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
220	>	Ylo - y;
221	>
222		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
223		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
224		}
#	Line 363 | Line 226 | TH2D * JetFinderAnalyzer::make_histo(con
226		}
227		}
228		// set histogram to match smear vector
229	<	for (int i = 1; i <= 60; i++) {
230	<	for (int j = 1; j <= 60; j++) {
229	>	for (int i = 1; i <= 30; i++) {
230	>	for (int j = 1; j <= 30; j++) {
231		histo->SetBinContent(i, j, smeared[i-1][j-1]);
232		}
233		}
234	<
235	<	return histo;
236	<	}
237	<
238	<	void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
239	<	jetfit::model_def &_mdef) {
377	<	// create a PseudoJet vector
378	<	vector<PseudoJet> particles;
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]->e();
389	<	double eta = gParticles[i]->eta();
390	<	double phi = gParticles[i]->phi();
391	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
392	<	PseudoJet j(pmom);
393	<	particles.push_back(j);
234	>	}
235	>	else {
236	>	// don't smear--just fill with particles
237	>	for (int i = 0; i < particles.size(); i++) {
238	>	histo->Fill(particles[i]->eta(), particles[i]->phi(),
239	>	particles[i]->energy());
240		}
241		}
242
243	<	// choose a jet definition
244	<	double R = 0.2;
399	<	JetDefinition jet_def(cambridge_algorithm, R);
400	<
401	<	// run clustering and extract the jets
402	<	ClusterSequence cs(particles, jet_def);
403	<	vector<PseudoJet> jets = cs.inclusive_jets();
404	<
405	<	double XbinSize = (histo->GetXaxis()->GetXmax()
406	<	- histo->GetXaxis()->GetXmin()) /
407	<	static_cast<double>(histo->GetXaxis()->GetNbins());
408	<	double YbinSize = (histo->GetYaxis()->GetXmax()
409	<	- histo->GetYaxis()->GetXmin()) /
410	<	static_cast<double>(histo->GetYaxis()->GetNbins());
243	>	return histo;
244	>	}
245
246	<	// seed with C-A jets
246	>	template <class Jet>
247	>	void seed_with_jetcoll(vector<Jet> jets,
248	>	jetfit::model_def &_mdef, double phi_lo = -PI,
249	>	double phi_hi = PI) {
250	>	// seed with jet collection
251		int ijset = 0;
252		for (unsigned ij = 0; ij < jets.size(); ij++) {
253	<	double N = jets[ij].e();
254	<	if (N > 50.0) {
253	>	double N = jets[ij].energy();
254	>	double eta = jets[ij].eta();
255	>	double phi = jets[ij].phi();
256	>	if (N > 10.0) {
257		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
258		0.0, 1.0e6);
259	<	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
259	>	_mdef.set_special_par(ijset, 1, eta, 0.01,
260		0.0, 0.0);
261	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
262	<	: jets[ij].phi();
263	<	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
261	>	if (phi < phi_lo) {
262	>	phi += 2.0*PI;
263	>	}
264	>	if (phi > phi_hi) {
265	>	phi -= 2.0*PI;
266	>	}
267	>	_mdef.set_special_par(ijset, 2, phi, 0.01,
268		0.0, 0.0);
269		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
270		ijset++;
#	Line 434 | Line 278 | jetfit::model_def& JetFinderAnalyzer::ma
278		class jf_model_def : public jetfit::model_def {
279		public:
280		virtual double chisquare_error(double E) {
281	<	return 0.97*E + 14.0;
282	<	// study from 08-27-09
281	>	return exp(-(4.2 + 0.11*E));
282	>	// study from 09-04-09
283		}
284		};
285
#	Line 451 | Line 295 | jetfit::model_def& JetFinderAnalyzer::ma
295		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
296		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
297
298	<	seed_with_CA(get_particles(evt), histo, *_mdef);
298	>	// get jetcoll from event file and select highest e jet
299	>	if (info_type_ == 0) {
300	>	edm::Handle< vector<reco::GenJet> > jet_collection;
301	>	evt.getByLabel(jet_algo_, jet_collection);
302	>	reco::GenJet highest_e_jet;
303	>	bool found_jet = false;
304	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
305	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
306	>	{
307	>	highest_e_jet = (*jet_collection)[i];
308	>	found_jet = true;
309	>	}
310	>	}
311	>	vector<reco::GenJet> highest_e_jet_coll;
312	>	highest_e_jet_coll.push_back(highest_e_jet);
313	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
314	>	histo->GetYaxis()->GetXmin(),
315	>	histo->GetYaxis()->GetXmax());
316	>	}
317	>	if (info_type_ == 1) {
318	>	edm::Handle< vector<reco::PFJet> > jet_collection;
319	>	evt.getByLabel(jet_algo_, jet_collection);
320	>	reco::PFJet highest_e_jet;
321	>	bool found_jet = false;
322	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
323	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
324	>	{
325	>	highest_e_jet = (*jet_collection)[i];
326	>	found_jet = true;
327	>	}
328	>	}
329	>	vector<reco::PFJet> highest_e_jet_coll;
330	>	highest_e_jet_coll.push_back(highest_e_jet);
331	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
332	>	histo->GetYaxis()->GetXmin(),
333	>	histo->GetYaxis()->GetXmax());
334	>	}
335
336		jetfit::set_model_def(_mdef);
337
#	Line 498 | Line 378 | jetfit::model_def& JetFinderAnalyzer::ma
378		}
379		jetfit::set_ngauss(_mdef->get_n_special_par_sets());
380		init_fit_histo->SetBinContent(i+1, j+1,
381	<	jetfit::fit_fcn(x, y, pval));
381	>	jetfit::fit_fcn(x, y, pval)
382	>	* XbinSize * YbinSize);
383		}
384		}
385
#	Line 551 | Line 432 | ostream& operator<<(ostream &out, jetfit
432		else
433		out << "Not calculated" << endl;
434		}
435	+	else {
436	+	out << "Error matrix accurate" << endl;
437	+	}
438
439		return out;
440		}
#	Line 559 | Line 443 | void JetFinderAnalyzer::analyze_results(
443		std::vector<jetfit::trouble> t,
444		TH2 *hist_orig) {
445		ofs << "Histogram "<<hist_orig->GetName() << endl;
446	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
446	>	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
447	>	i < unique_jets[hist_orig].size(); i++) {
448	>	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
449		ofs << "For "<<i+1<<" gaussians: " << endl
450	<	<< t.at(i) << endl
451	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
450	>	<< t.at(i) << endl;
451	>	ofs << "chisquare="<<r.chisquare.at(ir)
452	>	<< endl;
453	>	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
454		for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
455		jet _jet = unique_jets[hist_orig][i][j];
456		ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
#	Line 577 | Line 465 | void JetFinderAnalyzer::analyze_results(
465		for (vector< vector<double> >::size_type i = 0;
466		i < r.pval.size(); i++) {
467		jetfit::set_ngauss(r.pval[i].size() / 4);
580	–	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
581	–	hist_orig->GetXaxis()->GetXmin(),
582	–	hist_orig->GetXaxis()->GetXmax(),
583	–	hist_orig->GetYaxis()->GetXmin(),
584	–	hist_orig->GetYaxis()->GetXmax(),
585	–	r.pval[i].size());
586	–	for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
587	–	tf2->SetParameter(j, r.pval[i][j]);
588	–	}
468		ostringstream fit_histo_oss;
469		fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
591	–	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
592	–	tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
470		TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
471		fit_histo_oss.str().c_str(),
472		hist_orig->GetXaxis()->GetNbins(),
#	Line 598 | Line 475 | void JetFinderAnalyzer::analyze_results(
475		hist_orig->GetYaxis()->GetNbins(),
476		hist_orig->GetYaxis()->GetXmin(),
477		hist_orig->GetYaxis()->GetXmax());
478	<	TH1 *tf2_histo = tf2->CreateHistogram();
479	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
480	<	- fit_histo->GetXaxis()->GetXmin())
481	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
482	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
483	<	- fit_histo->GetYaxis()->GetXmin())
484	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
485	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
486	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
487	<	fit_histo->SetBinContent(ih+1, jh+1,
488	<	tf2_histo->GetBinContent(ih+1, jh+1)
489	<	* XbinSize * YbinSize);
478	>	double XbinSize = (hist_orig->GetXaxis()->GetXmax()
479	>	- hist_orig->GetXaxis()->GetXmin())
480	>	/ static_cast<double>(hist_orig->GetXaxis()->GetNbins());
481	>	double YbinSize = (hist_orig->GetYaxis()->GetXmax()
482	>	- hist_orig->GetYaxis()->GetXmin())
483	>	/ static_cast<double>(hist_orig->GetYaxis()->GetNbins());
484	>
485	>	double Xlo = hist_orig->GetXaxis()->GetXmin();
486	>	double Xhi = hist_orig->GetXaxis()->GetXmax();
487	>	double Ylo = hist_orig->GetYaxis()->GetXmin();
488	>	double Yhi = hist_orig->GetYaxis()->GetXmax();
489	>	for (int bi = 0; bi < 30; bi++) {
490	>	for (int bj = 0; bj < 30; bj++) {
491	>	double x = (static_cast<double>(bi) + 0.5)*XbinSize + Xlo;
492	>	double y = (static_cast<double>(bj) + 0.5)*YbinSize + Ylo;
493	>	double pval[256];
494	>	for (int ii = 0; ii < r.pval[i].size(); ii++) {
495	>	pval[ii] = r.pval[i][ii];
496	>	}
497	>	fit_histo->SetBinContent(bi+1, bj+1,
498	>	jetfit::fit_fcn(x, y, pval)
499	>	* XbinSize * YbinSize);
500		}
501		}
502		}

Diff Legend

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