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.5 by dnisson, Fri Sep 4 15:30:38 2009 UTC vs.
Revision 1.27 by dnisson, Fri Jan 22 23:53:10 2010 UTC

#	Line 1 | Line 1
1	<	/* A simple jet-finding analyzer */
1	>	/* A JetFitAnalyzer that makes histograms with smearing */
2
3		#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
5	–	#include "fastjet/ClusterSequence.hh"
5		#include "FWCore/ServiceRegistry/interface/Service.h"
6		#include "FWCore/MessageLogger/interface/MessageLogger.h"
7		#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
8
9	<	#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>
19	<	#include <cmath>
20	<	#include <cstdlib>
21	<	#include <fstream>
9	>	#include <iostream>
10		#include <sstream>
11
24	–	#include "TFormula.h"
12		#include "TF2.h"
13	+	#include "TNtuple.h"
14
15		#define PI 3.141593
16
17		using namespace std;
30	–	using namespace fastjet;
18
19		class JetFinderAnalyzer : public JetFitAnalyzer {
20		public:
34	–	struct jet {
35	–	double energy;
36	–	double eta;
37	–	double phi;
38	–	};
39	–
21		explicit JetFinderAnalyzer( const edm::ParameterSet&);
22		~JetFinderAnalyzer() {}
23
24		private:
25	<	static map<TH2 *, vector< vector<jet> > > unique_jets;
26	<
27	<	static double phi_cutoff_;
47	<
48	<	static double g2int(double xlo, double xhi, double ylo, double yhi,
49	<	double *pval) {
50	<	double sum1 = 0.0;
51	<	double sum2 = 0.0;
52	<	double xmid = 0.5 * (xlo + xhi);
53	<	double ymid = 0.5 * (ylo + yhi);
54	<	double xstep = (xhi - xlo) / 50.0;
55	<	double ystep = (yhi - ylo) / 50.0;
56	<	for (int i = 0; i < 50; i++) {
57	<	double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
58	<	sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
59	<	}
60	<	for (int i = 0; i < 50; i++) {
61	<	double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
62	<	sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
63	<	}
64	<	return sum1 * sum2;
65	<	}
66	<
67	<	static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
68	<	double dist_sq = numeric_limits<double>::infinity();
69	<	unique_jets[hist].resize(ngauss);
70	<	int nbinsX = hist->GetXaxis()->GetNbins();
71	<	int nbinsY = hist->GetYaxis()->GetNbins();
72	<	double XbinSize = (hist->GetXaxis()->GetXmax()
73	<	- hist->GetXaxis()->GetXmin())
74	<	/ static_cast<double>(nbinsX);
75	<	double YbinSize = (hist->GetYaxis()->GetXmax()
76	<	- hist->GetYaxis()->GetXmin())
77	<	/ static_cast<double>(nbinsY);
78	<	for (int i = 0; i < ngauss; i++) {
79	<	double N, mu_x, mu_y, sig, err, lo, hi;
80	<	int iuint;
81	<	TString name;
82	<	gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
83	<	gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
84	<	gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
85	<	gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
86	<	for (int j = 0; j < i; j++) {
87	<	double N2, mu_x2, mu_y2, sig2;
88	<	gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
89	<	gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
90	<	gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
91	<	gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
92	<	double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
93	<	+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
94	<	if (_dist_sq < dist_sq)
95	<	dist_sq = _dist_sq;
96	<	}
97	<
98	<	jet j;
99	<	j.energy = N;
100	<	j.eta = mu_x; j.phi = mu_y;
101	<	unique_jets[hist][ngauss-1].push_back(j);
102	<	}
103	<	}
104	<
105	<	virtual void beginJob(const edm::EventSetup&);
106	<	virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
107	<	virtual jetfit::model_def& make_model_def(const edm::Event&,
108	<	const edm::EventSetup&,
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_;
25	>	virtual void beginJob(const edm::EventSetup &es);
26	>	virtual void analyze_results(HistoFitter::FitResults, std::vector<HistoFitter::Trouble>,
27	>	TH2 *);
28		};
29
126	–	map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
127	–	JetFinderAnalyzer::unique_jets;
128	–
30		JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
31		: JetFitAnalyzer(pSet) // this is important!
32		{
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
144	–	// 1 = proper angle smear
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	–	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]->energy());
230	–	}
231	–
232	–	// reduce histo
233	–	ostringstream oss2;
234	–	oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
235	–	edm::Service<TFileService> fs;
236	–	// draw cone of radius 0.5 around highest energy bin, reduce
237	–	double maxE = 0.0;
238	–	int max_i = 29, max_j = 29;
239	–	for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
240	–	for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
241	–	double E = unred_histo->GetBinContent(i+1, j+1);
242	–	if (E > maxE) {
243	–	maxE = E;
244	–	max_i = i;
245	–	max_j = j;
246	–	}
247	–	}
248	–	}
249	–
250	–	double rcone = 0.5;
251	–	double Xlo = unred_histo->GetXaxis()->GetXmin();
252	–	double Xhi = unred_histo->GetXaxis()->GetXmax();
253	–	double Ylo = unred_histo->GetYaxis()->GetXmin();
254	–	double Yhi = unred_histo->GetYaxis()->GetXmax();
255	–	double XbinSize = (Xhi - Xlo) /
256	–	static_cast<double>(unred_histo->GetXaxis()->GetNbins());
257	–	double YbinSize = (Yhi - Ylo) /
258	–	static_cast<double>(unred_histo->GetYaxis()->GetNbins());
259	–	double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
260	–	double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
261	–	TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
262	–	60, max_x-rcone, max_x+rcone,
263	–	60, max_y-rcone, max_y+rcone);
264	–
265	–	// create an unsmeared reduced histo
266	–	TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
267	–	(oss2.str()+"_unsmeared").c_str(),
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]->energy();
272	–	double x = particles[i]->eta();
273	–	double y = particles[i]->phi();
274	–	histo_unsmeared->Fill(x, y, N);
275	–	}
276	–
277	–	// create a smeared reduced histo
278	–	// create a temporary 2D vector for smeared energies
279	–	XbinSize = (histo->GetXaxis()->GetXmax()
280	–	- histo->GetXaxis()->GetXmin()) /
281	–	static_cast<double>(histo->GetXaxis()->GetNbins());
282	–	YbinSize = (histo->GetYaxis()->GetXmax()
283	–	- histo->GetYaxis()->GetXmin()) /
284	–	static_cast<double>(histo->GetYaxis()->GetNbins());
285	–	vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
286	–	switch (smear_coord_) {
287	–	case 1:
288	–	for (int i = 0; i < particles.size(); i++) {
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;
300	–
301	–	// transform eta, phi to proper angle
302	–	double theta = 2.0*atan(exp(-eta));
303	–	double iota = asin(sin(theta)*sin(phi));
304	–
305	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
306	–	* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
307	–	}
308	–	}
309	–	}
310	–	break;
311	–	case 0:
312	–	default:
313	–	for (int i = 0; i < particles.size(); i++) {
314	–	double N = particles[i]->energy();
315	–	double x = particles[i]->eta();
316	–	double y = particles[i]->phi();
317	–	// loop over bins and add Gaussian to smeared
318	–	for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
319	–	for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
320	–	double eta = static_cast<double>((signed int)i2) * XbinSize
321	–	+ max_x - rcone - x;
322	–	double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
323	–	+ max_y - rcone - y));
324	–	phi = sin(phi) > 0 ? phi : -phi;
325	–	smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
326	–	* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
327	–	}
328	–	}
329	–	}
330	–	}
331	–	// set histogram to match smear vector
332	–	for (int i = 1; i <= 60; i++) {
333	–	for (int j = 1; j <= 60; j++) {
334	–	histo->SetBinContent(i, j, smeared[i-1][j-1]);
335	–	}
336	–	}
337	–
338	–	return histo;
339	–	}
340	–
341	–	void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo,
342	–	jetfit::model_def &_mdef) {
343	–	// create a PseudoJet vector
344	–	vector<PseudoJet> particles;
345	–	for (unsigned i = 0; i < gParticles.size(); i++) {
346	–	double x_max = (histo->GetXaxis()->GetXmax()
347	–	+ histo->GetXaxis()->GetXmin()) / 2.0;
348	–	double y_max = (histo->GetYaxis()->GetXmax()
349	–	+ histo->GetYaxis()->GetXmin()) / 2.0;
350	–	valarray<double> pmom(4);
351	–	pmom[0] = gParticles[i]->px();
352	–	pmom[1] = gParticles[i]->py();
353	–	pmom[2] = gParticles[i]->pz();
354	–	pmom[3] = gParticles[i]->energy();
355	–	double eta = gParticles[i]->eta();
356	–	double phi = gParticles[i]->phi();
357	–	if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
358	–	PseudoJet j(pmom);
359	–	particles.push_back(j);
360	–	}
361	–	}
362	–
363	–	// choose a jet definition
364	–	double R = 0.2;
365	–	JetDefinition jet_def(cambridge_algorithm, R);
366	–
367	–	// run clustering and extract the jets
368	–	ClusterSequence cs(particles, jet_def);
369	–	vector<PseudoJet> jets = cs.inclusive_jets();
370	–
371	–	double XbinSize = (histo->GetXaxis()->GetXmax()
372	–	- histo->GetXaxis()->GetXmin()) /
373	–	static_cast<double>(histo->GetXaxis()->GetNbins());
374	–	double YbinSize = (histo->GetYaxis()->GetXmax()
375	–	- histo->GetYaxis()->GetXmin()) /
376	–	static_cast<double>(histo->GetYaxis()->GetNbins());
377	–
378	–	// seed with C-A jets
379	–	int ijset = 0;
380	–	for (unsigned ij = 0; ij < jets.size(); ij++) {
381	–	double N = jets[ij].e();
382	–	if (N > 50.0) {
383	–	cout << N << endl;
384	–	_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
385	–	0.0, 1.0e6);
386	–	_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
387	–	0.0, 0.0);
388	–	double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
389	–	: jets[ij].phi();
390	–	_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
391	–	0.0, 0.0);
392	–	_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
393	–	ijset++;
394	–	}
395	–	}
396	–	}
397	–
398	–	jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
399	–	const edm::EventSetup&,
400	–	TH2 *histo) {
401	–	class jf_model_def : public jetfit::model_def {
402	–	public:
403	–	virtual double chisquare_error(double E) {
404	–	return E > 2.7 ? 0.55*E - 1.0 : 0.5;
405	–	// study from 09-03-09
406	–	}
407	–	};
408	–
409	–	jf_model_def *_mdef = new jf_model_def();
410	–	TFormula *formula = new TFormula("gaus2d",
411	–	"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
412	–	_mdef->set_formula(formula);
413	–	_mdef->set_indiv_max_E(0);
414	–	_mdef->set_indiv_max_x(1);
415	–	_mdef->set_indiv_max_y(2);
416	–	_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
417	–	_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
418	–	_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
419	–	_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
420	–
421	–	seed_with_CA(get_particles(evt), histo, *_mdef);
422	–
423	–	jetfit::set_model_def(_mdef);
424	–
425	–	// generate initial fit histogram
426	–	edm::Service<TFileService> fs;
427	–	TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
428	–	.c_str(),
429	–	("Initial fit for "
430	–	+string(histo->GetName())).c_str(),
431	–	histo->GetXaxis()->GetNbins(),
432	–	histo->GetXaxis()->GetXmin(),
433	–	histo->GetXaxis()->GetXmax(),
434	–	histo->GetXaxis()->GetNbins(),
435	–	histo->GetXaxis()->GetXmin(),
436	–	histo->GetXaxis()->GetXmax());
437	–	double XbinSize = (histo->GetXaxis()->GetXmax()
438	–	- histo->GetXaxis()->GetXmin()) /
439	–	static_cast<double>(histo->GetXaxis()->GetNbins());
440	–	double YbinSize = (histo->GetYaxis()->GetXmax()
441	–	- histo->GetYaxis()->GetXmin()) /
442	–	static_cast<double>(histo->GetYaxis()->GetNbins());
443	–	double Xlo = histo->GetXaxis()->GetXmin();
444	–	double Xhi = histo->GetXaxis()->GetXmax();
445	–	double Ylo = histo->GetYaxis()->GetXmin();
446	–	double Yhi = histo->GetYaxis()->GetXmax();
447	–
448	–	for (int i = 0; i < 60; i++) {
449	–	for (int j = 0; j < 60; j++) {
450	–	double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
451	–	double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
452	–	double pval[256];
453	–	if (_mdef->get_n_special_par_sets() > 64) {
454	–	cerr << "Parameter overload" << endl;
455	–	return *_mdef;
456	–	}
457	–	else {
458	–	for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
459	–	for (int ii = 0; ii < 4; ii++) {
460	–	double spval, sperr, splo, sphi;
461	–	_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
462	–	pval[4*is + ii] = spval;
463	–	}
464	–	}
465	–	}
466	–	jetfit::set_ngauss(_mdef->get_n_special_par_sets());
467	–	init_fit_histo->SetBinContent(i+1, j+1,
468	–	jetfit::fit_fcn(x, y, pval));
469	–	}
470	–	}
471	–
472	–	return *_mdef;
33		}
34
35		void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
476	–	ofs.open("jetfindlog.txt", ios::out);
477	–	if (ofs.fail()) {
478	–	cerr << "Opening jetfindlog.txt FAILED" << endl;
479	–	}
480	–	ofs << "Jetfinder log" << endl
481	–	<< "=============" << endl << endl;
482	–	}
483	–
484	–	ostream& operator<<(ostream &out, jetfit::trouble t) {
485	–	string action, error_string;
486	–
487	–	if (t.istat != 3) {
488	–	switch(t.occ) {
489	–	case jetfit::T_NULL:
490	–	action = "Program"; break;
491	–	case jetfit::T_SIMPLEX:
492	–	action = "SIMPLEX"; break;
493	–	case jetfit::T_MIGRAD:
494	–	action = "MIGRAD"; break;
495	–	case jetfit::T_MINOS:
496	–	action = "MINOS"; break;
497	–	default:
498	–	action = "Program"; break;
499	–	}
36
37	<	switch (t.istat) {
502	<	case 0:
503	<	error_string = "Unable to calculate error matrix"; break;
504	<	case 1:
505	<	error_string = "Error matrix a diagonal approximation"; break;
506	<	case 2:
507	<	error_string = "Error matrix not positive definite"; break;
508	<	case 3:
509	<	error_string = "Converged successfully"; break;
510	<	default:
511	<	ostringstream oss;
512	<	oss<<"Unknown status code "<<t.istat << flush;
513	<	error_string = oss.str(); break;
514	<	}
37	>	}
38
39	<	if (t.occ != jetfit::T_NULL)
40	<	out << action<<" trouble: "<<error_string;
41	<	else
42	<	out << "Not calculated" << endl;
39	>	double evalFitFunction(HistoFitter::FitResults r, double x, double y) {
40	>	unsigned nFits = r.pars.size();
41	>	unsigned nGauss = r.pars[nFits-1].size() / 4;
42	>	double fitVal = 0.0;
43	>	for (unsigned i = 0; i < nGauss; i++) {
44	>	double N = r.pval[nFits-1][4*i];
45	>	double mu_x = r.pval[nFits-1][4*i + 1];
46	>	double mu_y = r.pval[nFits-1][4*i + 2];
47	>	double sig = r.pval[nFits-1][4*i + 3];
48	>
49	>	double rel_x = x - mu_x; double rel_y = y - mu_y;
50	>	fitVal += (N / 2.0 / M_PI / sig / sig)
51	>	* exp(-(rel_x * rel_x + rel_y * rel_y)/2.0/sig/sig);
52		}
53	<
522	<	return out;
53	>	return fitVal;
54		}
55
56	<	void JetFinderAnalyzer::analyze_results(jetfit::results r,
57	<	std::vector<jetfit::trouble> t,
56	>	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
57	>	std::vector<HistoFitter::Trouble> t,
58		TH2 *hist_orig) {
59	<	ofs << "Histogram "<<hist_orig->GetName() << endl;
529	<	for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
530	<	ofs << "For "<<i+1<<" gaussians: " << endl
531	<	<< t.at(i) << endl
532	<	<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
533	<	for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
534	<	jet _jet = unique_jets[hist_orig][i][j];
535	<	ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
536	<	<< ", phi = "<<_jet.phi << endl;
537	<	}
538	<	ofs << endl;
539	<	}
540	<	ofs << endl;
541	<
542	<	// save fit function histograms to root file
59	>	// perform analysis of fit results
60		edm::Service<TFileService> fs;
61	<	for (vector< vector<double> >::size_type i = 0;
62	<	i < r.pval.size(); i++) {
63	<	jetfit::set_ngauss(r.pval[i].size() / 4);
64	<	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
65	<	hist_orig->GetXaxis()->GetXmin(),
66	<	hist_orig->GetXaxis()->GetXmax(),
67	<	hist_orig->GetYaxis()->GetXmin(),
68	<	hist_orig->GetYaxis()->GetXmax(),
69	<	r.pval[i].size());
70	<	for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
71	<	tf2->SetParameter(j, r.pval[i][j]);
72	<	}
73	<	ostringstream fit_histo_oss;
74	<	fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
75	<	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
76	<	tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
77	<	TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
78	<	fit_histo_oss.str().c_str(),
79	<	hist_orig->GetXaxis()->GetNbins(),
80	<	hist_orig->GetXaxis()->GetXmin(),
81	<	hist_orig->GetXaxis()->GetXmax(),
82	<	hist_orig->GetYaxis()->GetNbins(),
566	<	hist_orig->GetYaxis()->GetXmin(),
567	<	hist_orig->GetYaxis()->GetXmax());
568	<	TH1 *tf2_histo = tf2->CreateHistogram();
569	<	double XbinSize = (fit_histo->GetXaxis()->GetXmax()
570	<	- fit_histo->GetXaxis()->GetXmin())
571	<	/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
572	<	double YbinSize = (fit_histo->GetYaxis()->GetXmax()
573	<	- fit_histo->GetYaxis()->GetXmin())
574	<	/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
575	<	for (int ih = 0; ih < tf2->GetNpx(); ih++) {
576	<	for (int jh = 0; jh < tf2->GetNpy(); jh++) {
577	<	fit_histo->SetBinContent(ih+1, jh+1,
578	<	tf2_histo->GetBinContent(ih+1, jh+1)
579	<	* XbinSize * YbinSize);
580	<	}
61	>	TH2D *fitHisto = fs->make<TH2D>((std::string(hist_orig->GetName())+"_fit").c_str(),
62	>	("Fitted distribution to "
63	>	+std::string(hist_orig->GetName())).c_str(),
64	>	hist_orig->GetNbinsX(),
65	>	hist_orig->GetXaxis()->GetXmin(),
66	>	hist_orig->GetXaxis()->GetXmax(),
67	>	hist_orig->GetNbinsY(),
68	>	hist_orig->GetYaxis()->GetXmin(),
69	>	hist_orig->GetYaxis()->GetXmax());
70	>
71	>	double Xlo = fitHisto->GetXaxis()->GetXmin();
72	>	double Xhi = fitHisto->GetXaxis()->GetXmax();
73	>	double Ylo = fitHisto->GetYaxis()->GetXmin();
74	>	double Yhi = fitHisto->GetYaxis()->GetXmax();
75	>	double XbinSize = (Xhi - Xlo) / static_cast<double>(fitHisto->GetNbinsX());
76	>	double YbinSize = (Yhi - Ylo) / static_cast<double>(fitHisto->GetNbinsY());
77	>
78	>	for (int i = 1; i <= fitHisto->GetNbinsX(); i++) {
79	>	for (int j = 1; j <= fitHisto->GetNbinsY(); j++) {
80	>	double x = (static_cast<double>(i) - 0.5) * XbinSize + Xlo;
81	>	double y = (static_cast<double>(j) - 0.5) * YbinSize + Ylo;
82	>	fitHisto->SetBinContent(i, j, evalFitFunction(r, x, y) * XbinSize * YbinSize);
83		}
84		}
85
86	<	// save results to file
87	<	ostringstream res_tree_oss, rt_title_oss;
88	<	res_tree_oss << hist_orig->GetName()<<"_results" << flush;
89	<	rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
86	>	// save fit results to an ntuple
87	>	TNtuple *rNtuple = fs->make<TNtuple>((std::string(hist_orig->GetName())+"_results").c_str(),
88	>	("Fit results for "+std::string(hist_orig->GetName())).c_str(),
89	>	"N:mu_x:mu_y:sigma");
90	>	unsigned nFits = r.pval.size();
91	>	unsigned nGauss = r.pval[nFits-1].size() / 4;
92	>	for (unsigned i = 0; i < nGauss; i++) {
93	>	rNtuple->Fill(r.pval[nFits-1][4*i], r.pval[nFits-1][4*i+1], r.pval[nFits-1][4*i+2],
94	>	r.pval[nFits-1][4*i+3]);
95	>	}
96	>
97	>	// save chisquares to ntuple
98	>	for (unsigned i = 0; i < r.chisquare.size(); i++) {
99	>	ostringstream csNtupleName, csNtupleTitle;
100	>	csNtupleName << hist_orig->GetName() << "_chi2_" << i << flush;
101	>	csNtupleTitle << "Chisquare "<<i<<" for histo "<<hist_orig->GetName()
102	>	<< flush;
103	>	TNtuple *csNtuple = fs->make<TNtuple>(csNtupleName.str().c_str(),
104	>	csNtupleTitle.str().c_str(),
105	>	"chisq");
106	>	csNtuple->Fill(r.chisquare[i]);
107	>	}
108		}
109
110		DEFINE_FWK_MODULE(JetFinderAnalyzer);

Diff Legend

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