ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/JetFitAnalyzer/src/jetfit.cpp

Comparing UserCode/JetFitAnalyzer/src/jetfit.cpp (file contents):
Revision 1.1 by dnisson, Wed Sep 2 21:44:14 2009 UTC vs.
Revision 1.13 by dnisson, Fri Nov 13 03:57:48 2009 UTC

#	Line 1 | Line 1
1		/* jetfit.cpp - Package to fit multi-Gaussian distributions to histograms
2	–	* rewrite after accidental deletion 07-24-09
2		* Author: David Nisson
3		*/
4
#	Line 27 | Line 26
26		#include <Rtypes.h>
27		#include <TFormula.h>
28		#include <TSystem.h>
29	+	#include <TMath.h>
30
31		#include "UserCode/JetFitAnalyzer/interface/jetfit.h"
32
#	Line 40 | Line 40 | namespace jetfit {
40		// configurable options
41		bool ignorezero = false; // ignore zero bins when fitting
42
43	–	int ngauss; // number of Gaussians in current model
44	–
45	–	struct histo {
46	–	vector< vector<double> > bins;
47	–	double Xlo, Xhi, Ylo, Yhi;
48	–	} energy;
49	–
50	–
43		model_def *mdef;
44
45		model_def& curr_model_def() {
#	Line 58 | Line 50 | namespace jetfit {
50		mdef = _mdef;
51		}
52
53	<	int get_ngauss() {
53	>	int model_def::get_ngauss() {
54		return ngauss;
55		}
56
57	<	void set_ngauss(int _ngauss) {
57	>	void model_def::set_ngauss(int _ngauss) {
58		ngauss = _ngauss;
59		}
60
61		// fit function
62	<	double fit_fcn(double x, double y, double *xval) {
63	<	int npar_indiv = mdef->get_formula()->GetNpar();
62	>	double model_def::fit_fcn(double x, double y, double *xval) {
63	>	int npar_indiv = get_formula()->GetNpar();
64		double val = 0.0;
65		for (int i = 0; i < ngauss; i++) {
66	<	mdef->get_formula()->SetParameters(xval + i*npar_indiv);
66	>	get_formula()->SetParameters(xval + i*npar_indiv);
67		val += mdef->get_formula()->Eval(x, y);
68		}
69		return val;
#	Line 79 | Line 71 | namespace jetfit {
71
72		// TF2-compatible fit function
73		double fit_fcn_TF2(double *x, double *pval) {
74	<	double val = fit_fcn(x[0], x[1], pval);
74	>	double val = mdef->fit_fcn(x[0], x[1], pval);
75		return val;
76		}
77
78	<	// Integral of model formula / chisquare sigma
79	<	double formula_int(double xlo, double xhi, double ylo, double yhi,
78	>	// Integral of (model formula)^2 / chisquare sigma
79	>	double model_def::formula_int(double xlo, double xhi, double ylo, double yhi,
80		double *pval, double XbinSize, double YbinSize,
81		double *xval) {
82		double xstep = (xhi - xlo) / static_cast<double>(20);
#	Line 92 | Line 84 | namespace jetfit {
84
85		double fsum = 0.0;
86		double pval_old[256];
87	<	int npar = mdef->get_formula()->GetNpar();
87	>	int npar = get_formula()->GetNpar();
88		if (npar > 256) {
89		cerr << "Parameter overload" << endl;
90		return 0.0;
91		}
92	<	mdef->get_formula()->GetParameters(pval_old);
93	<	mdef->get_formula()->SetParameters(pval);
92	>	get_formula()->GetParameters(pval_old);
93	>	get_formula()->SetParameters(pval);
94
95		for (int i = 0; i < 20; i++) {
96		for (int j = 0; j < 20; j++) {
#	Line 106 | Line 98 | namespace jetfit {
98		double y = (static_cast<double>(j) + 0.5)*ystep + ylo;
99		double sig_cut = 1.0e-5;
100		double nu = XbinSize * YbinSize * fit_fcn(x, y, xval);
101	<	double sigma = mdef->chisquare_error(nu) < sig_cut
102	<	? sig_cut : mdef->chisquare_error(nu);
101	>	double sigma = fabs(chisquare_error(0.0)) < sig_cut
102	>	? sig_cut : chisquare_error(0.0);
103
104	<	fsum += xstep * ystep * mdef->get_formula()->Eval(x, y)
104	>	fsum += pow(xstep * ystep * get_formula()->Eval(x, y), 2.0)
105		/ sigma / sigma;
106		}
107		}
108
109	<	mdef->get_formula()->SetParameters(pval_old);
109	>	get_formula()->SetParameters(pval_old);
110		return fsum;
111		}
112
#	Line 130 | Line 122 | namespace jetfit {
122		// add errors in data points in histogram
123		for (int i = 0; i < energy.bins.size(); i++) {
124		for (int j = 0; j < energy.bins.at(i).size(); j++) {
125	<	double x = (static_cast<double>(i) + 0.5) * (energy.Xhi - energy.Xlo)
125	>	double x = (static_cast<double>(i) + 0.5) * (energy.Xhi - energy.Xlo)
126		/ static_cast<double>(energy.bins.size()) + energy.Xlo;
127		double y = (static_cast<double>(j) + 0.5) * (energy.Yhi - energy.Ylo)
128		/ static_cast<double>(energy.bins.at(i).size()) + energy.Ylo;
129	<	double nu = fit_fcn(x, y, xval) * XbinSize * YbinSize;
129	>	double nu = mdef->fit_fcn(x, y, xval) * XbinSize * YbinSize;
130		double sig_cut = 1.0e-5;
131	<	if (mdef->chisquare_error(nu) > sig_cut || !ignorezero)
131	>	if (fabs(mdef->chisquare_error(energy.bins.at(i).at(j))) > sig_cut || !ignorezero) {
132		chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
133	<	/ pow(mdef->chisquare_error(nu), 2.0);
134	<	else
133	>	/ pow(mdef->chisquare_error(energy.bins.at(i).at(j)), 2.0);
134	>	}
135	>	else {
136		chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
137		/ pow(sig_cut, 2.0);
138	+	}
139		}
140		}
141	<
141	>
142		// add errors due to Gaussians outside histogram
143	<	double eps = 0.1; // accuracy set for this function
144	<	for (int i = 0; i < ngauss; i++) {
145	<	double *pval = xval+i*mdef->get_formula()->GetNpar();
143	>	double eps = 0.01; // accuracy set for this function
144	>	for (int i = 0; i < mdef->get_ngauss(); i++) {
145	>	double *pval = xval+i*(mdef->get_formula()->GetNpar());
146		double par_x = pval[mdef->get_indiv_max_x()];
147		double par_y = pval[mdef->get_indiv_max_y()];
148		double par_sig = pval[3];
149	<	double cutoff_rad = -log(eps * 2.0 * PI * par_sig * par_sig
150	<	/ pval[0]);
149	>	double cutoff_rad = sqrt(-log(eps * 2.0 * PI * par_sig * par_sig
150	>	/ pval[0]) * 2.0 * par_sig * par_sig);
151
152		bool left = par_x < energy.Xlo + cutoff_rad;
153		bool right = par_x > energy.Xhi - cutoff_rad;
154		bool top = par_y > energy.Yhi - cutoff_rad;
155	<	bool bottom = par_y > energy.Ylo + cutoff_rad;
155	>	bool bottom = par_y < energy.Ylo + cutoff_rad;
156
157		if (left) {
158		double xlo = par_x - cutoff_rad;
159		double xhi = energy.Xlo;
160		double ylo = energy.Ylo;
161		double yhi = energy.Yhi;
162	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
162	>	chisquare += 2.0*mdef->formula_int(xlo, xhi, ylo, yhi, pval,
163		XbinSize, YbinSize, xval);
164		}
165		if (right) {
#	Line 173 | Line 167 | namespace jetfit {
167		double xhi = par_x + cutoff_rad;
168		double ylo = energy.Ylo;
169		double yhi = energy.Yhi;
170	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
170	>	chisquare += 2.0*mdef->formula_int(xlo, xhi, ylo, yhi, pval,
171		XbinSize, YbinSize, xval);
172		}
173		if (top) {
#	Line 181 | Line 175 | namespace jetfit {
175		double xhi = right ? par_x + cutoff_rad : energy.Xhi;
176		double ylo = energy.Yhi;
177		double yhi = par_y + cutoff_rad;
178	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
178	>	chisquare += 2.0*mdef->formula_int(xlo, xhi, ylo, yhi, pval,
179		XbinSize, YbinSize, xval);
180		}
181		if (bottom) {
#	Line 189 | Line 183 | namespace jetfit {
183		double xhi = right ? par_x + cutoff_rad : energy.Xhi;
184		double ylo = par_y - cutoff_rad;
185		double yhi = energy.Ylo;
186	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
186	>	chisquare += 2.0*mdef->formula_int(xlo, xhi, ylo, yhi, pval,
187		XbinSize, YbinSize, xval);
188		}
189		}
190	<
190	>
191		fcnval = chisquare;
192		}
193		catch (exception ex) {
#	Line 318 | Line 312 | namespace jetfit {
312		}
313		}
314
315	<	// translation of CERN's ERFC function
316	<	double erfc(double x) {
317	<	bool lef = false;
318	<	double p2[5], q2[5];
325	<
326	<	const double z1 = 1.0, hf = z1/2.0, c1 = 0.56418958;
327	<	const double vmax = 7.0;
328	<	const double switch_ = 4.0;
329	<
330	<	double p10 = 3.6767877, q10 = 3.2584593, p11 = -9.7970465e-2;
331	<	p2[0] = 7.3738883;
332	<	q2[0] = 7.3739609;
333	<	p2[1] = 6.8650185;
334	<	q2[1] = 1.5184908e1;
335	<	p2[2] = 3.0317993;
336	<	q2[2] = 1.2795530e1;
337	<	p2[3] = 5.6316962e-1;
338	<	q2[3] = 5.3542168;
339	<	p2[4] = 4.3187787e-5;
340	<	q2[4] = 1.0;
341	<
342	<	double p30 = -1.2436854e-1, q30 = 4.4091706e-1, p31 = -9.6821036e-2;
343	<
344	<	double v = fabs(x);
345	<	double y, h, hc, ap, aq;
346	<	if (v < hf) {
347	<	y = v*v;
348	<	h = x*(p10 + p11*y)/(q10 + y);
349	<	hc = 1.0 - h;
350	<	}
351	<	else {
352	<	if (v < switch_) {
353	<	ap = p2[4];
354	<	aq = q2[4];
355	<	for (int i = 3; i >= 0; i--) {
356	<	ap = p2[i] + v*ap;
357	<	aq = q2[i] + v*aq;
358	<	}
359	<	hc = exp(-v*v)*ap/aq;
360	<	h = 1.0 - hc;
361	<	}
362	<	else if (v < vmax) {
363	<	y = 1/v/v;
364	<	hc = exp(-v*v)*(c1 + y*(p30 + p31*y)/(q30 + y))/v;
365	<	h = 1.0 - hc;
366	<	}
367	<	// for very big values we can save us any calculation
368	<	// and the FP-exceptions from exp
369	<	else {
370	<	h = 1.0;
371	<	hc = 0.0;
315	>	void fill_histo_with_vec(TH2 *hist, vector< vector<double> > &vec) {
316	>	for (int i = 0; i < vec.size(); i++) {
317	>	for (int j = 0; j < vec[i].size(); j++) {
318	>	hist->SetBinContent(i+1, j+1, vec[i][j]);
319		}
373	–	if (x < 0) {
374	–	h = -h;
375	–	hc = 2.0 - hc;
376	–	}
377	–	}
378	–	if (lef) {
379	–	return h;
380	–	}
381	–	else {
382	–	return hc;
320		}
321		}
322
323	<	// translation of CERN's PROB function
387	<	double prob(double x, int n) {
388	<	const char *name = "PROB";
389	<	char errtxt[80];
390	<
391	<	const double r1 = 1.0,
392	<	hf = r1/2.0, th = r1/3.0, f1 = 2.0*r1/9.0;
393	<	const double c1 = 1.128379167095513;
394	<	const int nmax = 300;
395	<	// maximum chi2 per df for df >= 2., if chi2/df > chipdf prob=0
396	<	const double chipdf = 100.0;
397	<	const double xmax = 174.673, xmax2 = 2.0*xmax;
398	<	const double xlim = 24.0;
399	<	const double eps = 1e-30;
400	<
401	<	double y = x;
402	<	double u = hf*y;
403	<	double h, w, s, t, fi, e;
404	<	int m;
405	<	if (n < 0) {
406	<	h = 0.0;
407	<	sprintf(errtxt, "n = %d < 1", n);
408	<	cerr << "PROB: G100.1: "<<errtxt;
409	<	}
410	<	else if (y < 0.0) {
411	<	h = 0.0;
412	<	sprintf(errtxt, "x = %f < 0", n);
413	<	cerr << "PROB: G100.2: "<<errtxt;
414	<	}
415	<	else if (y == 0.0 || n/20 > y) {
416	<	h = 1.0;
417	<	}
418	<	else if (n == 1) {
419	<	w = sqrt(u);
420	<	if (w < xlim) {
421	<	h = erfc(w);
422	<	}
423	<	else {
424	<	h = 0.0;
425	<	}
426	<	}
427	<	else if (n > nmax) {
428	<	s = r1 / ((double)n);
429	<	t = f1 * s;
430	<	w = (pow(y*s, th) - (1.0 - t)) / sqrt(2.0*t);
431	<	if (w < -xlim) {
432	<	h = 1.0;
433	<	}
434	<	else if (w < xlim) {
435	<	h = hf * erfc(w);
436	<	}
437	<	else {
438	<	h = 0.0;
439	<	}
440	<	}
441	<	else {
442	<	m = n/2;
443	<	if (u < xmax2 && (y/n) <= chipdf) {
444	<	s = exp(-hf*u);
445	<	t = s;
446	<	e = s;
447	<	if (2*m == n) {
448	<	fi = 0.0;
449	<	for (int i = 1; i < m; i++) {
450	<	fi += 1.0;
451	<	t = u*t/fi;
452	<	s += t;
453	<	}
454	<	h = s*e;
455	<	}
456	<	else {
457	<	fi = 1.0;
458	<	for (int i = 1; i < m; i++) {
459	<	fi += 2.0;
460	<	t = t*y/fi;
461	<	s += t;
462	<	}
463	<	w = sqrt(u);
464	<	if (w < xlim) {
465	<	h = c1*w*s*e + erfc(w);
466	<	}
467	<	else {
468	<	h = 0.0;
469	<	}
470	<	}
471	<	}
472	<	else {
473	<	h = 0.0;
474	<	}
475	<	}
476	<	if (h > eps) {
477	<	return h;
478	<	}
479	<	else {
480	<	return 0.0;
481	<	}
482	<	}
483	<
484	<	void par_init(TH2 *hist, TMinuit *gMinuit, vector<TString> &pars,
323	>	void model_def::par_init(TH2 *hist, TMinuit *gMinuit, vector<TString> &pars,
324		double *pval, double *perr, double *plo, double *phi,
325	<	int npar, results r) {
326	<	int npar1 = npar - mdef->get_formula()->GetNpar();
325	>	int npar, FitResults r) {
326	>	int npar1 = npar - get_formula()->GetNpar();
327		bool init_spec_pars = false;
328	<	if (ngauss <= mdef->get_n_special_par_sets()) {
328	>	if (ngauss <= get_n_special_par_sets()) {
329		init_spec_pars = true;
330		for (int k = 0; k < ngauss; k++) {
331	<	int ipar0 = k*mdef->get_formula()->GetNpar(); // index of par 0
332	<	for (int l = 0; l < mdef->get_formula()->GetNpar(); l++) {
331	>	int ipar0 = k*get_formula()->GetNpar(); // index of indiv par 0
332	>	for (int l = 0; l < get_formula()->GetNpar(); l++) {
333		int ipar = ipar0 + l;
334		if (ipar < npar)
335	<	mdef->get_special_par(k, l, pval[ipar], perr[ipar], plo[ipar],
335	>	get_special_par(k, l, pval[ipar], perr[ipar], plo[ipar],
336		phi[ipar]);
337		else
338		cerr << "WARNING: Attempt to set parameter out of index range!"
#	Line 509 | Line 348 | namespace jetfit {
348		vector< vector<double> > sub_energy(energy.bins.size());
349		double max_sub, max_x = 0.0, max_y = 0.0;
350		double pval_other[256];
351	<	max_sub = -(numeric_limits<double>::infinity());
351	>	max_sub = -(numeric_limits<double>::max());
352		for (int i = 0; i < energy.bins.size(); i++) {
353		sub_energy[i].resize(energy.bins[i].size());
354		for (int j = 0; j < energy.bins[i].size(); j++) {
#	Line 518 | Line 357 | namespace jetfit {
357		double y = static_cast<double>(j)*(energy.Yhi - energy.Ylo)
358		/ static_cast<double>(energy.bins[i].size()) + energy.Ylo;
359		if (ngauss > 1) {
360	<	// subtract 4*sigma plus integral of fit function
360	>	// subtract integral of fit function
361		try {
362		int npval_other = r.pval.at(r.pval.size()-1).size();
363		if (npval_other > 256) {
#	Line 531 | Line 370 | namespace jetfit {
370		ngauss--;
371		double nu = fit_fcn(x, y, pval_other) * XbinSize * YbinSize;
372		ngauss++;
373	<	sub_energy[i][j] = energy.bins[i][j] - nu
535	<	- 4.0*mdef->chisquare_error(nu);
373	>	sub_energy[i][j] = energy.bins[i][j] - nu;
374		}
375		catch (exception ex) {
376		cerr << "Exception in par_init" << endl;
#	Line 541 | Line 379 | namespace jetfit {
379		else {
380		sub_energy[i][j] = energy.bins[i][j];
381		}
544	–	if (sub_energy[i][j] > max_sub) {
545	–	max_sub = sub_energy[i][j];
546	–	max_x = x;
547	–	max_y = y;
548	–	}
382		}
383		}
384	+	TH2D *hist_sub = new TH2D("hist_sub", "Subtracted histo",
385	+	sub_energy.size(), energy.Xlo,
386	+	energy.Xhi,
387	+	sub_energy[0].size(), energy.Ylo,
388	+	energy.Yhi);
389	+	fill_histo_with_vec(hist_sub, sub_energy);
390	+	max_sub = 0.0;
391	+	max_x = 0.0;
392	+	max_y = 0.0;
393	+	for (int i = 1; i <= hist_sub->GetXaxis()->GetNbins(); i++) {
394	+	for (int j = 1; j <= hist_sub->GetYaxis()->GetNbins(); j++) {
395	+	double nu = energy.bins[i-1][j-1] - sub_energy[i-1][j-1];
396	+	if (hist_sub->GetBinContent(i, j)
397	+	- 3.0*pow(chisquare_error(nu), 2.0) > max_sub) {
398	+	max_sub = hist_sub->GetBinContent(i, j);
399	+	max_x = static_cast<double>(i-1)*XbinSize
400	+	+ hist_sub->GetXaxis()->GetXmin();
401	+	max_y = static_cast<double>(j-1)*YbinSize
402	+	+ hist_sub->GetYaxis()->GetXmin();
403	+	}
404	+	}
405	+	}
406
407		for (int i = npar1; i < npar; i++) {
408	<	if (mdef->get_indiv_max_E() == i - npar1) {
409	<	double nu = 0.0;
410	<	if (ngauss > 1) {
556	<	nu = fit_fcn(max_x, max_y, pval_other) * XbinSize * YbinSize;
557	<	}
558	<	pval[i] = max_sub + (ngauss > 1 ? 4.0*mdef->chisquare_error(nu)
559	<	: 0.0);
560	<	perr[i] = mdef->chisquare_error(pval[i])*0.1;
408	>	if (get_indiv_max_E() == i - npar1) {
409	>	pval[i] = max_sub / XbinSize / YbinSize;
410	>	perr[i] = mdef->chisquare_error(pval[i])*0.5;
411		plo[i] = 0.0;
412		phi[i] = 1.0e6;
413		}
414	<	else if (mdef->get_indiv_max_x() == i - npar1) {
414	>	else if (get_indiv_max_x() == i - npar1) {
415		pval[i] = max_x;
416		perr[i] = (energy.Xhi - energy.Xlo) / static_cast<double>(hist->GetXaxis()->GetNbins());
417		plo[i] = 0.0;
418		phi[i] = 0.0;
419		}
420	<	else if (mdef->get_indiv_max_y() == i - npar1) {
420	>	else if (get_indiv_max_y() == i - npar1) {
421		pval[i] = max_y;
422		perr[i] = (energy.Yhi - energy.Ylo) / static_cast<double>(hist->GetYaxis()->GetNbins());
423		plo[i] = 0.0;
#	Line 575 | Line 425 | namespace jetfit {
425		}
426		else {
427		string dummy;
428	<	mdef->get_indiv_par(i-npar1, dummy, pval[i], perr[i], plo[i], phi[i]);
428	>	get_indiv_par(i-npar1, dummy, pval[i], perr[i], plo[i], phi[i]);
429		}
430		}
431		}
432	<	int n_pars_to_set = init_spec_pars ? ngauss*mdef->get_formula()->GetNpar()
432	>	int n_pars_to_set = init_spec_pars ? ngauss*get_formula()->GetNpar()
433		: npar - npar1;
434		int init_par_to_set = init_spec_pars ? 0 : npar1;
435		for (int i = 0; i < n_pars_to_set; i++) {
436		double dummy;
437		string s;
438		int ipar = i + init_par_to_set;
439	<	mdef->get_indiv_par(i % mdef->get_formula()->GetNpar(), s,
439	>	get_indiv_par(i % get_formula()->GetNpar(), s,
440		dummy, dummy, dummy, dummy);
441		ostringstream par_oss;
442	<	par_oss << s << ipar / (mdef->get_formula()->GetNpar()) << flush;
442	>	par_oss << s << ipar / (get_formula()->GetNpar()) << flush;
443		try {
444		pars.at(ipar) = TString(par_oss.str());
445		}
#	Line 606 | Line 456 | namespace jetfit {
456		}
457		}
458
459	<	results fit_histo(TH2 *hist, vector<trouble> &t_vec,
460	<	void (*cc_minuit)(TMinuit *, TH2 *, int),
461	<	int start_ngauss,
462	<	int rebinX, int rebinY,
463	<	double P_cutoff_val) {
459	>	FitResults fit_histo(TH2 *hist, vector<trouble> &t_vec,
460	>	void (*cc_minuit)(TMinuit *, TH2 *, int),
461	>	int start_ngauss,
462	>	int rebinX, int rebinY,
463	>	double P_cutoff_val) {
464		TMinuit *gMinuit = new TMinuit();
465		int npar_indiv = mdef->get_formula()->GetNpar();
466		int istat, erflg;
467	<	results r;
467	>	FitResults r;
468
469		gMinuit->SetFCN(fcn);
470		gMinuit->mninit(5, 6, 7);
#	Line 651 | Line 501 | namespace jetfit {
501		start_ngauss = 1;
502		}
503
504	<	for (ngauss = start_ngauss;
505	<	ngauss <= (MAX_GAUSS > start_ngauss ? MAX_GAUSS : start_ngauss);
506	<	ngauss++) {
504	>	for (mdef->set_ngauss(start_ngauss);
505	>	mdef->get_ngauss() <=
506	>	(MAX_GAUSS > start_ngauss ? MAX_GAUSS : start_ngauss);
507	>	mdef->set_ngauss(mdef->get_ngauss()+1)) {
508	>
509	>	int ngauss = mdef->get_ngauss();
510		t_vec.resize(t_vec.size() + 1);
511	<	int npar = npar_indiv*ngauss;
511	>	int npar = npar_indiv*mdef->get_ngauss();
512		double pval[256], perr[256], plo[256], phi[256];
513		if (npar > 256) {
514		cerr << "Parameter overload" << endl;
#	Line 664 | Line 517 | namespace jetfit {
517		vector<TString> pars(npar);
518
519		// initialize parameters
520	<	par_init(hist, gMinuit, pars, pval, perr, plo, phi, npar, r);
521	<
520	>	mdef->par_init(hist, gMinuit, pars, pval, perr, plo, phi, npar, r);
521	>
522		// minimize
523		double chisquare, edm, errdef;
524		int nvpar, nparx;
525		trouble t;
526		t.occ = T_NULL;
527		t.istat = 3;
528	<	// fix the N values and sigmas of the Gaussians
528	>	// fix sigmas of the Gaussians
529		for (int i = 0; i < ngauss; i++) {
530		double parno[2];
531	<	parno[0] = i*npar_indiv + 1;
532	<	parno[1] = i*npar_indiv + 4;
680	<	gMinuit->mnexcm("FIX", parno, 2, erflg);
531	>	parno[0] = i*npar_indiv + 4;
532	>	gMinuit->mnexcm("FIX", parno, 1, erflg);
533		}
534		gMinuit->mnexcm("SIMPLEX", 0, 0, erflg);
535		gMinuit->mnexcm("MIGRAD", 0, 0, erflg);
684	–	// release N values, fix mu_x and mu_y
685	–	for (int i = 0; i < ngauss; i++) {
686	–	double parno[4];
687	–	parno[0] = i*npar_indiv + 1;
688	–	parno[1] = i*npar_indiv + 2;
689	–	parno[2] = i*npar_indiv + 3;
690	–	parno[3] = i*npar_indiv + 4;
691	–	gMinuit->mnexcm("RELEASE", parno, 1, erflg);
692	–	gMinuit->mnexcm("FIX", parno+1, 2, erflg);
693	–	}
694	–	gMinuit->mnexcm("MIGRAD", 0, 0, erflg);
695	–	// release mu_x and mu_y
696	–	for (int i = 0; i < ngauss; i++) {
697	–	double parno[4];
698	–	parno[0] = i*npar_indiv + 1;
699	–	parno[1] = i*npar_indiv + 2;
700	–	parno[2] = i*npar_indiv + 3;
701	–	parno[3] = i*npar_indiv + 4;
702	–	gMinuit->mnexcm("RELEASE", parno+1, 2, erflg);
703	–	}
704	–	gMinuit->mnexcm("MIGRAD", 0, 0, erflg);
536		gMinuit->mnstat(chisquare, edm, errdef, nvpar, nparx, istat);
537		if (istat == 3) {
538		/* we're not concerned about MINOS errors right now
#	Line 766 | Line 597 | namespace jetfit {
597		}
598		ndof -= npar;
599
600	<	double P = prob(chisquare, ndof);
601	<	cout << "P = "<<P << endl;
600	>	r.chisquare.push_back(chisquare);
601	>	double P = TMath::Prob(chisquare, ndof);
602		if (P > P_cutoff_val) {
603		break;
604		}

Diff Legend

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