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.4 by dnisson, Wed Sep 9 15:17:45 2009 UTC vs.
Revision 1.9 by dnisson, Tue Nov 10 17:03:15 2009 UTC

#	Line 83 | Line 83 | namespace jetfit {
83		return val;
84		}
85
86	<	// Integral of model formula / chisquare sigma
86	>	// Integral of (model formula)^2 / chisquare sigma
87		double formula_int(double xlo, double xhi, double ylo, double yhi,
88		double *pval, double XbinSize, double YbinSize,
89		double *xval) {
#	Line 106 | Line 106 | namespace jetfit {
106		double y = (static_cast<double>(j) + 0.5)*ystep + ylo;
107		double sig_cut = 1.0e-5;
108		double nu = XbinSize * YbinSize * fit_fcn(x, y, xval);
109	<	double sigma = fabs(mdef->chisquare_error(nu)) < sig_cut
110	<	? sig_cut : mdef->chisquare_error(nu);
109	>	double sigma = fabs(mdef->chisquare_error(0.0)) < sig_cut
110	>	? sig_cut : mdef->chisquare_error(0.0);
111
112	<	fsum += xstep * ystep * mdef->get_formula()->Eval(x, y)
112	>	fsum += pow(xstep * ystep * mdef->get_formula()->Eval(x, y), 2.0)
113		/ sigma / sigma;
114		}
115		}
#	Line 136 | Line 136 | namespace jetfit {
136		/ static_cast<double>(energy.bins.at(i).size()) + energy.Ylo;
137		double nu = fit_fcn(x, y, xval) * XbinSize * YbinSize;
138		double sig_cut = 1.0e-5;
139	<	if (fabs(mdef->chisquare_error(nu)) > sig_cut || !ignorezero) {
139	>	if (fabs(mdef->chisquare_error(energy.bins.at(i).at(j))) > sig_cut || !ignorezero) {
140		chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
141	<	/ pow(mdef->chisquare_error(nu), 2.0);
141	>	/ pow(mdef->chisquare_error(energy.bins.at(i).at(j)), 2.0);
142		}
143		else {
144		chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
#	Line 146 | Line 146 | namespace jetfit {
146		}
147		}
148		}
149	<
149	>
150		// add errors due to Gaussians outside histogram
151	<	double eps = 0.1; // accuracy set for this function
151	>	double eps = 0.01; // accuracy set for this function
152		for (int i = 0; i < ngauss; i++) {
153		double *pval = xval+i*(mdef->get_formula()->GetNpar());
154		double par_x = pval[mdef->get_indiv_max_x()];
155		double par_y = pval[mdef->get_indiv_max_y()];
156		double par_sig = pval[3];
157	<	double cutoff_rad = -log(eps * 2.0 * PI * par_sig * par_sig
158	<	/ pval[0]);
157	>	double cutoff_rad = sqrt(-log(eps * 2.0 * PI * par_sig * par_sig
158	>	/ pval[0]) * 2.0 * par_sig * par_sig);
159
160		bool left = par_x < energy.Xlo + cutoff_rad;
161		bool right = par_x > energy.Xhi - cutoff_rad;
#	Line 167 | Line 167 | namespace jetfit {
167		double xhi = energy.Xlo;
168		double ylo = energy.Ylo;
169		double yhi = energy.Yhi;
170	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
170	>	chisquare += 2.0*formula_int(xlo, xhi, ylo, yhi, pval,
171		XbinSize, YbinSize, xval);
172		}
173		if (right) {
#	Line 175 | Line 175 | namespace jetfit {
175		double xhi = par_x + cutoff_rad;
176		double ylo = energy.Ylo;
177		double yhi = energy.Yhi;
178	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
178	>	chisquare += 2.0*formula_int(xlo, xhi, ylo, yhi, pval,
179		XbinSize, YbinSize, xval);
180		}
181		if (top) {
#	Line 183 | Line 183 | namespace jetfit {
183		double xhi = right ? par_x + cutoff_rad : energy.Xhi;
184		double ylo = energy.Yhi;
185		double yhi = par_y + cutoff_rad;
186	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
186	>	chisquare += 2.0*formula_int(xlo, xhi, ylo, yhi, pval,
187		XbinSize, YbinSize, xval);
188		}
189		if (bottom) {
#	Line 191 | Line 191 | namespace jetfit {
191		double xhi = right ? par_x + cutoff_rad : energy.Xhi;
192		double ylo = par_y - cutoff_rad;
193		double yhi = energy.Ylo;
194	<	chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
194	>	chisquare += 2.0*formula_int(xlo, xhi, ylo, yhi, pval,
195		XbinSize, YbinSize, xval);
196		}
197		}
198	<
198	>
199		fcnval = chisquare;
200		}
201		catch (exception ex) {
#	Line 483 | Line 483 | namespace jetfit {
483		}
484		}
485
486	+	void fill_histo_with_vec(TH2 *hist, vector< vector<double> > &vec) {
487	+	for (int i = 0; i < vec.size(); i++) {
488	+	for (int j = 0; j < vec[i].size(); j++) {
489	+	hist->SetBinContent(i+1, j+1, vec[i][j]);
490	+	}
491	+	}
492	+	}
493	+
494		void par_init(TH2 *hist, TMinuit *gMinuit, vector<TString> &pars,
495		double *pval, double *perr, double *plo, double *phi,
496		int npar, results r) {
#	Line 491 | Line 499 | namespace jetfit {
499		if (ngauss <= mdef->get_n_special_par_sets()) {
500		init_spec_pars = true;
501		for (int k = 0; k < ngauss; k++) {
502	<	int ipar0 = k*mdef->get_formula()->GetNpar(); // index of par 0
502	>	int ipar0 = k*mdef->get_formula()->GetNpar(); // index of indiv par 0
503		for (int l = 0; l < mdef->get_formula()->GetNpar(); l++) {
504		int ipar = ipar0 + l;
505		if (ipar < npar)
#	Line 520 | Line 528 | namespace jetfit {
528		double y = static_cast<double>(j)*(energy.Yhi - energy.Ylo)
529		/ static_cast<double>(energy.bins[i].size()) + energy.Ylo;
530		if (ngauss > 1) {
531	<	// subtract 4*sigma plus integral of fit function
531	>	// subtract integral of fit function
532		try {
533		int npval_other = r.pval.at(r.pval.size()-1).size();
534		if (npval_other > 256) {
#	Line 533 | Line 541 | namespace jetfit {
541		ngauss--;
542		double nu = fit_fcn(x, y, pval_other) * XbinSize * YbinSize;
543		ngauss++;
544	<	sub_energy[i][j] = energy.bins[i][j] - nu
537	<	- 2.2*mdef->chisquare_error(nu);
544	>	sub_energy[i][j] = energy.bins[i][j] - nu;
545		}
546		catch (exception ex) {
547		cerr << "Exception in par_init" << endl;
#	Line 543 | Line 550 | namespace jetfit {
550		else {
551		sub_energy[i][j] = energy.bins[i][j];
552		}
546	–	if (sub_energy[i][j] > max_sub) {
547	–	max_sub = sub_energy[i][j];
548	–	max_x = x;
549	–	max_y = y;
550	–	}
553		}
554		}
555	+	TH2D *hist_sub = new TH2D("hist_sub", "Subtracted histo",
556	+	sub_energy.size(), energy.Xlo,
557	+	energy.Xhi,
558	+	sub_energy[0].size(), energy.Ylo,
559	+	energy.Yhi);
560	+	fill_histo_with_vec(hist_sub, sub_energy);
561	+	max_sub = 0.0;
562	+	max_x = 0.0;
563	+	max_y = 0.0;
564	+	for (int i = 1; i <= hist_sub->GetXaxis()->GetNbins(); i++) {
565	+	for (int j = 1; j <= hist_sub->GetYaxis()->GetNbins(); j++) {
566	+	double nu = energy.bins[i-1][j-1] - sub_energy[i-1][j-1];
567	+	if (hist_sub->GetBinContent(i, j)
568	+	- 3.0*pow(mdef->chisquare_error(nu), 2.0) > max_sub) {
569	+	max_sub = hist_sub->GetBinContent(i, j);
570	+	max_x = static_cast<double>(i-1)*XbinSize
571	+	+ hist_sub->GetXaxis()->GetXmin();
572	+	max_y = static_cast<double>(j-1)*YbinSize
573	+	+ hist_sub->GetYaxis()->GetXmin();
574	+	}
575	+	}
576	+	}
577
578		for (int i = npar1; i < npar; i++) {
579		if (mdef->get_indiv_max_E() == i - npar1) {
580	<	double nu = 0.0;
557	<	if (ngauss > 1) {
558	<	ngauss--;
559	<	nu = fit_fcn(max_x, max_y, pval_other) * XbinSize * YbinSize;
560	<	ngauss++;
561	<	}
562	<	pval[i] = (max_sub + (ngauss > 1 ? 2.2*mdef->chisquare_error(nu)
563	<	: 0.0)
564	<	- mdef->chisquare_error(nu))
565	<	* 2.0 * PI * 0.01 / XbinSize / YbinSize;
580	>	pval[i] = max_sub / XbinSize / YbinSize;
581		perr[i] = mdef->chisquare_error(pval[i])*0.5;
582		plo[i] = 0.0;
583		phi[i] = 1.0e6;
#	Line 671 | Line 686 | namespace jetfit {
686
687		// initialize parameters
688		par_init(hist, gMinuit, pars, pval, perr, plo, phi, npar, r);
689	<
689	>
690		// minimize
691		double chisquare, edm, errdef;
692		int nvpar, nparx;

Diff Legend

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