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.20 by dnisson, Wed Nov 11 01:46:56 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	+	/* 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>
25		#include <cmath>
26		#include <cstdlib>
27	+	#include <iostream>
28		#include <fstream>
29		#include <sstream>
30
#	Line 39 | Line 52 | private:
52
53		static double phi_cutoff_;
54
55	<	static double g2int(double xlo, double xhi, double ylo, double yhi,
56	<	double *pval) {
44	<	double sum1 = 0.0;
45	<	double sum2 = 0.0;
46	<	double xmid = 0.5 * (xlo + xhi);
47	<	double ymid = 0.5 * (ylo + yhi);
48	<	double xstep = (xhi - xlo) / 50.0;
49	<	double ystep = (yhi - ylo) / 50.0;
50	<	for (int i = 0; i < 50; i++) {
51	<	double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
52	<	sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
53	<	}
54	<	for (int i = 0; i < 50; i++) {
55	<	double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
56	<	sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
57	<	}
58	<	return sum1 * sum2;
59	<	}
60	<
61	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
62	<	double dist_sq = numeric_limits<double>::infinity();
63	<	unique_jets[hist].resize(ngauss);
64	<	int nbinsX = hist->GetXaxis()->GetNbins();
65	<	int nbinsY = hist->GetYaxis()->GetNbins();
66	<	double XbinSize = (hist->GetXaxis()->GetXmax()
67	<	- hist->GetXaxis()->GetXmin())
68	<	/ static_cast<double>(nbinsX);
69	<	double YbinSize = (hist->GetYaxis()->GetXmax()
70	<	- hist->GetYaxis()->GetXmin())
71	<	/ static_cast<double>(nbinsY);
72	<	for (int i = 0; i < ngauss; i++) {
73	<	double N, mu_x, mu_y, sig, err, lo, hi;
74	<	int iuint;
75	<	TString name;
76	<	gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
77	<	gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
78	<	gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
79	<	gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
80	<	for (int j = 0; j < i; j++) {
81	<	double N2, mu_x2, mu_y2, sig2;
82	<	gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
83	<	gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
84	<	gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
85	<	gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
86	<	double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
87	<	+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
88	<	if (_dist_sq < dist_sq)
89	<	dist_sq = _dist_sq;
90	<	}
91	<
92	<	jet j;
93	<	j.energy = N;
94	<	j.eta = mu_x; j.phi = mu_y;
95	<	unique_jets[hist][ngauss-1].push_back(j);
96	<	}
97	<	}
98	<
99	<	virtual void beginJob(const edm::EventSetup&);
100	<	virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
55	>	virtual void beginJob(const edm::EventSetup &es);
56	>	virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
57		virtual jetfit::model_def& make_model_def(const edm::Event&,
58		const edm::EventSetup&,
59		TH2 *);
60		virtual void analyze_results(jetfit::results r,
61		std::vector<jetfit::trouble> t, TH2 *);
62	+	vector<reco::Candidate *> get_particles(const edm::Event&);
63
64		fstream ofs;
65	+	edm::InputTag inputTagPFCandidates_;
66	+	edm::InputTag inputTagGenParticles_;
67	+	int info_type_;
68		double smear_;
69		int smear_coord_;
70	+	string jet_algo_;
71	+	string reclustered_jet_algo_;
72		};
73
74		map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
#	Line 115 | Line 77 | JetFinderAnalyzer::unique_jets;
77		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
78		: JetFitAnalyzer(pSet) // this is important!
79		{
80	+	info_type_ = pSet.getUntrackedParameter("info_type", 0);
81	+
82	+	if (info_type_ == 0) {
83	+	inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
84	+	}
85	+	if (info_type_ == 1) {
86	+	inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
87	+	}
88	+
89		smear_ = pSet.getUntrackedParameter("smear", 0.02);
90		smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
91		// 0 = eta-phi smear
92		// 1 = proper angle smear
93	<	set_user_minuit(jetfinder);
93	>	jet_algo_ = pSet.getParameter<string>("jet_algo");
94	>	set_user_minuit(0); // we don't need this feature--it's deprecated
95		}
96
97		TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
98	<	ostringstream oss;
99	<	oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
100	<	TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
101	<	600, -2.5, 2.5, 600, -PI, PI);
102	<
103	<	// fill unreduced histo
104	<	edm::Handle<edm::HepMCProduct> hRaw;
105	<	evt.getByLabel("source",
106	<	hRaw);
107	<	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);
98	>	const reco::Jet * highest_e_jet = 0;
99	>
100	>	if (info_type_ == 0) {
101	>	edm::Handle< vector<reco::GenJet> > evtJets;
102	>	evt.getByLabel(jet_algo_, evtJets);
103	>	for (unsigned i = 0; i < evtJets->size(); i++) {
104	>	if (highest_e_jet == 0
105	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
106	>	highest_e_jet = &((*evtJets)[i]);
107	>	}
108		}
109		}
110	<
111	<	for (int i = 0; i < particles.size(); i++) {
112	<	unred_histo->Fill(particles[i]->momentum().eta(),
113	<	particles[i]->momentum().phi(),
114	<	particles[i]->momentum().e());
115	<	}
116	<
151	<	// reduce histo
152	<	ostringstream oss2;
153	<	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
154	<	edm::Service<TFileService> fs;
155	<	// draw cone of radius 0.5 around highest energy bin, reduce
156	<	double maxE = 0.0;
157	<	int max_i = 29, max_j = 29;
158	<	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
159	<	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
160	<	double E = unred_histo->GetBinContent(i+1, j+1);
161	<	if (E > maxE) {
162	<	maxE = E;
163	<	max_i = i;
164	<	max_j = j;
110	>	else if (info_type_ == 1) {
111	>	edm::Handle< vector<reco::PFJet> > evtJets;
112	>	evt.getByLabel(jet_algo_, evtJets);
113	>	for (unsigned i = 0; i < evtJets->size(); i++) {
114	>	if (highest_e_jet == 0
115	>	|| (*evtJets)[i].energy() > highest_e_jet->energy()) {
116	>	highest_e_jet = &((*evtJets)[i]);
117		}
118		}
119		}
120	<
121	<	double rcone = 0.5;
122	<	double Xlo = unred_histo->GetXaxis()->GetXmin();
123	<	double Xhi = unred_histo->GetXaxis()->GetXmax();
172	<	double Ylo = unred_histo->GetYaxis()->GetXmin();
173	<	double Yhi = unred_histo->GetYaxis()->GetXmax();
174	<	double XbinSize = (Xhi - Xlo) /
175	<	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
176	<	double YbinSize = (Yhi - Ylo) /
177	<	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
178	<	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
179	<	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
180	<	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
181	<	60, max_x-rcone, max_x+rcone,
182	<	60, max_y-rcone, max_y+rcone);
183	<
184	<	// create an unsmeared reduced histo
185	<	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
186	<	(oss2.str()+"_unsmeared").c_str(),
187	<	60, max_x-rcone, max_x+rcone,
188	<	60, max_y-rcone, max_y+rcone);
189	<	for (int i = 0; i < particles.size(); i++) {
190	<	double N = particles[i]->momentum().e();
191	<	double x = particles[i]->momentum().eta();
192	<	double y = particles[i]->momentum().phi();
193	<	histo_unsmeared->Fill(x, y, N);
120	>
121	>	if (highest_e_jet == 0) {
122	>	cerr << "No fat jets found!" << endl;
123	>	return 0;
124		}
125
126	<	// create a smeared reduced histo
126	>	vector<const reco::Candidate *> particles =
127	>	highest_e_jet->getJetConstituentsQuick();
128	>
129	>	ostringstream oss;
130	>	oss << "eta_phi_energy"<<evt.id().event() << flush;
131	>	edm::Service<TFileService> fs;
132	>	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
133	>	30,
134	>	highest_e_jet->eta()-0.5,
135	>	highest_e_jet->eta()+0.5,
136	>	30,
137	>	highest_e_jet->phi()-0.5,
138	>	highest_e_jet->phi()+0.5);
139	>
140	>	if (smear_ > 0.0) {
141	>	// create a smeared histo
142		// create a temporary 2D vector for smeared energies
143	<	XbinSize = (histo->GetXaxis()->GetXmax()
143	>	double XbinSize = (histo->GetXaxis()->GetXmax()
144		- histo->GetXaxis()->GetXmin()) /
145		static_cast<double>(histo->GetXaxis()->GetNbins());
146	<	YbinSize = (histo->GetYaxis()->GetXmax()
146	>	double YbinSize = (histo->GetYaxis()->GetXmax()
147		- histo->GetYaxis()->GetXmin()) /
148		static_cast<double>(histo->GetYaxis()->GetNbins());
149	<	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
149	>	double Xlo = histo->GetXaxis()->GetXmin();
150	>	double Xhi = histo->GetXaxis()->GetXmax();
151	>	double Ylo = histo->GetYaxis()->GetXmin();
152	>	double Yhi = histo->GetYaxis()->GetXmax();
153	>	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
154		switch (smear_coord_) {
155		case 1:
156		for (int i = 0; i < particles.size(); i++) {
157	<	double N = particles[i]->momentum().e();
158	<	double x = particles[i]->momentum().eta();
159	<	double y = particles[i]->momentum().phi();
157	>	double N = particles[i]->energy();
158	>	double x = particles[i]->eta();
159	>	double y = particles[i]->phi();
160	>	if (y < Ylo) {
161	>	y += 2.0*PI;
162	>	}
163	>	if (y > Yhi) {
164	>	y -= 2.0*PI;
165	>	}
166		// loop over bins and add Gaussian in proper angle to smeared
167	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
168	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
167	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
168	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
169		double eta = static_cast<double>((signed int)i2) * XbinSize +
170	<	max_x - rcone - x;
171	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
172	<	max_y - rcone - y));
173	<	phi = sin(phi) > 0 ? phi : -phi;
219	<
170	>	Xlo - x;
171	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
172	>	Ylo - y;
173	>
174		// transform eta, phi to proper angle
175		double theta = 2.0*atan(exp(-eta));
176		double iota = asin(sin(theta)*sin(phi));
#	Line 230 | Line 184 | TH2D * JetFinderAnalyzer::make_histo(con
184		case 0:
185		default:
186		for (int i = 0; i < particles.size(); i++) {
187	<	double N = particles[i]->momentum().e();
188	<	double x = particles[i]->momentum().eta();
189	<	double y = particles[i]->momentum().phi();
187	>	double N = particles[i]->energy();
188	>	double x = particles[i]->eta();
189	>	double y = particles[i]->phi();
190	>	if (y < Ylo) {
191	>	y += 2.0*PI;
192	>	}
193	>	if (y > Yhi) {
194	>	y -= 2.0*PI;
195	>	}
196		// loop over bins and add Gaussian to smeared
197	<	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
198	<	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
197	>	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
198	>	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
199		double eta = static_cast<double>((signed int)i2) * XbinSize
200	<	+ max_x - rcone - x;
201	<	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
202	<	+ max_y - rcone - y));
203	<	phi = sin(phi) > 0 ? phi : -phi;
200	>	+ Xlo - x;
201	>	double phi = static_cast<double>((signed int)j2) * YbinSize +
202	>	Ylo - y;
203	>
204		smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
205		* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
206		}
#	Line 248 | Line 208 | TH2D * JetFinderAnalyzer::make_histo(con
208		}
209		}
210		// set histogram to match smear vector
211	<	for (int i = 1; i <= 60; i++) {
212	<	for (int j = 1; j <= 60; j++) {
211	>	for (int i = 1; i <= 30; i++) {
212	>	for (int j = 1; j <= 30; j++) {
213		histo->SetBinContent(i, j, smeared[i-1][j-1]);
214		}
215		}
216	<
217	<	return histo;
218	<	}
219	<
220	<	void seed_with_CA(const edm::Event &evt, TH2 *histo,
221	<	jetfit::model_def &_mdef) {
262	<	edm::Handle<edm::HepMCProduct> hMC;
263	<	evt.getByLabel("source", hMC);
264	<
265	<	// create a PseudoJet vector
266	<	vector<PseudoJet> particles;
267	<	const HepMC::GenEvent *hmcEvt = hMC->GetEvent();
268	<	for (HepMC::GenEvent::particle_const_iterator
269	<	pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
270	<	pit++) {
271	<	if ((*pit)->status() == 1) {
272	<	double x_max = (histo->GetXaxis()->GetXmax()
273	<	+ histo->GetXaxis()->GetXmin()) / 2.0;
274	<	double y_max = (histo->GetYaxis()->GetXmax()
275	<	+ histo->GetYaxis()->GetXmin()) / 2.0;
276	<	valarray<double> pmom(4);
277	<	pmom[0] = (*pit)->momentum().px();
278	<	pmom[1] = (*pit)->momentum().py();
279	<	pmom[2] = (*pit)->momentum().pz();
280	<	pmom[3] = (*pit)->momentum().e();
281	<	double eta = (*pit)->momentum().eta();
282	<	double phi = (*pit)->momentum().phi();
283	<	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
284	<	PseudoJet j(pmom);
285	<	particles.push_back(j);
286	<	}
216	>	}
217	>	else {
218	>	// don't smear--just fill with particles
219	>	for (int i = 0; i < particles.size(); i++) {
220	>	histo->Fill(particles[i]->eta(), particles[i]->phi(),
221	>	particles[i]->energy());
222		}
223		}
224
225	<	// choose a jet definition
226	<	double R = 0.2;
292	<	JetDefinition jet_def(cambridge_algorithm, R);
293	<
294	<	// run clustering and extract the jets
295	<	ClusterSequence cs(particles, jet_def);
296	<	vector<PseudoJet> jets = cs.inclusive_jets();
297	<
298	<	double XbinSize = (histo->GetXaxis()->GetXmax()
299	<	- histo->GetXaxis()->GetXmin()) /
300	<	static_cast<double>(histo->GetXaxis()->GetNbins());
301	<	double YbinSize = (histo->GetYaxis()->GetXmax()
302	<	- histo->GetYaxis()->GetXmin()) /
303	<	static_cast<double>(histo->GetYaxis()->GetNbins());
225	>	return histo;
226	>	}
227
228	<	// seed with C-A jets
228	>	template <class Jet>
229	>	void seed_with_jetcoll(vector<Jet> jets,
230	>	jetfit::model_def &_mdef, double phi_lo = -PI,
231	>	double phi_hi = PI) {
232	>	// seed with jet collection
233		int ijset = 0;
234		for (unsigned ij = 0; ij < jets.size(); ij++) {
235	<	double N = jets[ij].e();
235	>	double N = jets[ij].energy();
236	>	double eta = jets[ij].eta();
237	>	double phi = jets[ij].phi();
238		if (N > 10.0) {
239		_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
240		0.0, 1.0e6);
241	<	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
241	>	_mdef.set_special_par(ijset, 1, eta, 0.01,
242		0.0, 0.0);
243	<	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
244	<	: jets[ij].phi();
245	<	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
243	>	if (phi < phi_lo) {
244	>	phi += 2.0*PI;
245	>	}
246	>	if (phi > phi_hi) {
247	>	phi -= 2.0*PI;
248	>	}
249	>	_mdef.set_special_par(ijset, 2, phi, 0.01,
250		0.0, 0.0);
251		_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
252		ijset++;
#	Line 327 | Line 260 | jetfit::model_def& JetFinderAnalyzer::ma
260		class jf_model_def : public jetfit::model_def {
261		public:
262		virtual double chisquare_error(double E) {
263	<	return 0.97*E + 14.0;
331	<	// study from 08-27-09
263	>	return 0.298 - 0.9557*E; // study 11-09-09
264		}
265		};
266
#	Line 344 | Line 276 | jetfit::model_def& JetFinderAnalyzer::ma
276		_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
277		_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
278
279	<	seed_with_CA(evt, histo, *_mdef);
280	<
281	<	jetfit::set_model_def(_mdef);
279	>	// get jetcoll from event file and select highest e jet
280	>	if (info_type_ == 0) {
281	>	edm::Handle< vector<reco::GenJet> > jet_collection;
282	>	evt.getByLabel(jet_algo_, jet_collection);
283	>	reco::GenJet highest_e_jet;
284	>	bool found_jet = false;
285	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
286	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
287	>	{
288	>	highest_e_jet = (*jet_collection)[i];
289	>	found_jet = true;
290	>	}
291	>	}
292	>	vector<reco::GenJet> highest_e_jet_coll;
293	>	highest_e_jet_coll.push_back(highest_e_jet);
294	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
295	>	histo->GetYaxis()->GetXmin(),
296	>	histo->GetYaxis()->GetXmax());
297	>	}
298	>	if (info_type_ == 1) {
299	>	edm::Handle< vector<reco::PFJet> > jet_collection;
300	>	evt.getByLabel(jet_algo_, jet_collection);
301	>	reco::PFJet highest_e_jet;
302	>	bool found_jet = false;
303	>	for (unsigned i = 0; i < jet_collection->size(); i++) {
304	>	if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
305	>	{
306	>	highest_e_jet = (*jet_collection)[i];
307	>	found_jet = true;
308	>	}
309	>	}
310	>	vector<reco::PFJet> highest_e_jet_coll;
311	>	highest_e_jet_coll.push_back(highest_e_jet);
312	>	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
313	>	histo->GetYaxis()->GetXmin(),
314	>	histo->GetYaxis()->GetXmax());
315	>	}
316
317		// generate initial fit histogram
318		edm::Service<TFileService> fs;
#	Line 371 | Line 337 | jetfit::model_def& JetFinderAnalyzer::ma
337		double Ylo = histo->GetYaxis()->GetXmin();
338		double Yhi = histo->GetYaxis()->GetXmax();
339
340	<	for (int i = 0; i < 60; i++) {
341	<	for (int j = 0; j < 60; j++) {
340	>	for (int i = 0; i < 30; i++) {
341	>	for (int j = 0; j < 30; j++) {
342		double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
343		double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
344		double pval[256];
#	Line 389 | Line 355 | jetfit::model_def& JetFinderAnalyzer::ma
355		}
356		}
357		}
358	<	jetfit::set_ngauss(_mdef->get_n_special_par_sets());
358	>	_mdef->set_ngauss(_mdef->get_n_special_par_sets());
359		init_fit_histo->SetBinContent(i+1, j+1,
360	<	jetfit::fit_fcn(x, y, pval));
360	>	_mdef->fit_fcn(x, y, pval)
361	>	* XbinSize * YbinSize);
362		}
363		}
364
#	Line 444 | Line 411 | ostream& operator<<(ostream &out, jetfit
411		else
412		out << "Not calculated" << endl;
413		}
414	+	else {
415	+	out << "Error matrix accurate" << endl;
416	+	}
417
418		return out;
419		}
#	Line 452 | Line 422 | void JetFinderAnalyzer::analyze_results(
422		std::vector<jetfit::trouble> t,
423		TH2 *hist_orig) {
424		ofs << "Histogram "<<hist_orig->GetName() << endl;
425	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
425	>	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
426	>	i < unique_jets[hist_orig].size(); i++) {
427	>	int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
428		ofs << "For "<<i+1<<" gaussians: " << endl
429	<	<< t.at(i) << endl
430	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
429	>	<< t.at(i) << endl;
430	>	ofs << "chisquare="<<r.chisquare.at(ir)
431	>	<< endl;
432	>	ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
433		for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
434		jet _jet = unique_jets[hist_orig][i][j];
435		ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
#	Line 469 | Line 443 | void JetFinderAnalyzer::analyze_results(
443		edm::Service<TFileService> fs;
444		for (vector< vector<double> >::size_type i = 0;
445		i < r.pval.size(); i++) {
446	<	jetfit::set_ngauss(r.pval[i].size() / 4);
446	>	r.moddef->set_ngauss(r.pval[i].size() / 4);
447	>	jetfit::set_model_def(r.moddef);
448		TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
449		hist_orig->GetXaxis()->GetXmin(),
450		hist_orig->GetXaxis()->GetXmax(),
#	Line 491 | Line 466 | void JetFinderAnalyzer::analyze_results(
466		hist_orig->GetYaxis()->GetNbins(),
467		hist_orig->GetYaxis()->GetXmin(),
468		hist_orig->GetYaxis()->GetXmax());
469	<	TH1 *tf2_histo = tf2->CreateHistogram();
470	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
471	<	- fit_histo->GetXaxis()->GetXmin())
472	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
473	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
474	<	- fit_histo->GetYaxis()->GetXmin())
475	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
476	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
477	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
478	<	fit_histo->SetBinContent(ih+1, jh+1,
479	<	tf2_histo->GetBinContent(ih+1, jh+1)
480	<	* XbinSize * YbinSize);
469	>	double XbinSize = (hist_orig->GetXaxis()->GetXmax()
470	>	- hist_orig->GetXaxis()->GetXmin())
471	>	/ hist_orig->GetXaxis()->GetNbins();
472	>	double YbinSize = (hist_orig->GetYaxis()->GetXmax()
473	>	- hist_orig->GetYaxis()->GetXmin())
474	>	/ hist_orig->GetYaxis()->GetNbins();
475	>	for (int i = 1; i <= fit_histo->GetXaxis()->GetNbins(); i++) {
476	>	for (int j = 1; j <= fit_histo->GetYaxis()->GetNbins(); j++) {
477	>	double x = static_cast<double>(i - 1)
478	>	* (fit_histo->GetXaxis()->GetXmax()
479	>	- fit_histo->GetXaxis()->GetXmin())
480	>	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins())
481	>	+ fit_histo->GetXaxis()->GetXmin();;
482	>	double y = static_cast<double>(j - 1)
483	>	* (fit_histo->GetYaxis()->GetXmax()
484	>	- fit_histo->GetYaxis()->GetXmin())
485	>	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins())
486	>	+ fit_histo->GetYaxis()->GetXmin();
487	>	fit_histo->SetBinContent(i, j, tf2->Eval(x, y) * XbinSize * YbinSize);
488		}
489		}
490		}

Diff Legend

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