ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/JetFitAnalyzer/src/jetfit.cpp
Revision: 1.6
Committed: Wed Sep 23 17:52:05 2009 UTC (15 years, 7 months ago) by dnisson
Branch: MAIN
CVS Tags: V00-01-00
Changes since 1.5: +5 -6 lines
Log Message:
Fixed bug in subtraction code.

File Contents

# User Rev Content
1 dnisson 1.1 /* jetfit.cpp - Package to fit multi-Gaussian distributions to histograms
2     * rewrite after accidental deletion 07-24-09
3     * Author: David Nisson
4     */
5    
6     #include <cstdlib>
7     #include <cmath>
8     #include <ctime>
9     #include <cstdio>
10     #include <cstring>
11     #include <cctype>
12     #include <iostream>
13     #include <fstream>
14     #include <sstream>
15     #include <string>
16     #include <vector>
17     #include <set>
18     #include <map>
19     #include <limits>
20     #include <exception>
21    
22     #include <TMinuit.h>
23     #include <TH1.h>
24     #include <TH2.h>
25     #include <TF2.h>
26     #include <TRandom3.h>
27     #include <Rtypes.h>
28     #include <TFormula.h>
29     #include <TSystem.h>
30    
31     #include "UserCode/JetFitAnalyzer/interface/jetfit.h"
32    
33     #define PI 3.141592653
34    
35     #define MAX_GAUSS 3
36    
37     using namespace std;
38    
39     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    
51     model_def *mdef;
52    
53     model_def& curr_model_def() {
54     return *mdef;
55     }
56    
57     void set_model_def(model_def *_mdef) {
58     mdef = _mdef;
59     }
60    
61     int get_ngauss() {
62     return ngauss;
63     }
64    
65     void set_ngauss(int _ngauss) {
66     ngauss = _ngauss;
67     }
68    
69     // fit function
70     double fit_fcn(double x, double y, double *xval) {
71     int npar_indiv = mdef->get_formula()->GetNpar();
72     double val = 0.0;
73     for (int i = 0; i < ngauss; i++) {
74     mdef->get_formula()->SetParameters(xval + i*npar_indiv);
75     val += mdef->get_formula()->Eval(x, y);
76     }
77     return val;
78     }
79    
80     // TF2-compatible fit function
81     double fit_fcn_TF2(double *x, double *pval) {
82     double val = fit_fcn(x[0], x[1], pval);
83     return val;
84     }
85    
86     // Integral of model formula / chisquare sigma
87     double formula_int(double xlo, double xhi, double ylo, double yhi,
88     double *pval, double XbinSize, double YbinSize,
89     double *xval) {
90     double xstep = (xhi - xlo) / static_cast<double>(20);
91     double ystep = (yhi - ylo) / static_cast<double>(20);
92    
93     double fsum = 0.0;
94     double pval_old[256];
95     int npar = mdef->get_formula()->GetNpar();
96     if (npar > 256) {
97     cerr << "Parameter overload" << endl;
98     return 0.0;
99     }
100     mdef->get_formula()->GetParameters(pval_old);
101     mdef->get_formula()->SetParameters(pval);
102    
103     for (int i = 0; i < 20; i++) {
104     for (int j = 0; j < 20; j++) {
105     double x = (static_cast<double>(i) + 0.5)*xstep + xlo;
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 dnisson 1.2 double sigma = fabs(mdef->chisquare_error(nu)) < sig_cut
110 dnisson 1.1 ? sig_cut : mdef->chisquare_error(nu);
111    
112     fsum += xstep * ystep * mdef->get_formula()->Eval(x, y)
113     / sigma / sigma;
114     }
115     }
116    
117     mdef->get_formula()->SetParameters(pval_old);
118     return fsum;
119     }
120    
121     // MINUIT objective function
122     void fcn(int &npar, double *grad, double &fcnval, double *xval, int iflag)
123     {
124     double chisquare = 0.0;
125     double XbinSize = (energy.Xhi - energy.Xlo)
126     / static_cast<double>(energy.bins.size());
127     double YbinSize = (energy.Yhi - energy.Ylo)
128     / static_cast<double>(energy.bins.at(0).size());
129     try {
130     // add errors in data points in histogram
131     for (int i = 0; i < energy.bins.size(); i++) {
132     for (int j = 0; j < energy.bins.at(i).size(); j++) {
133 dnisson 1.4 double x = (static_cast<double>(i) + 0.5) * (energy.Xhi - energy.Xlo)
134 dnisson 1.1 / static_cast<double>(energy.bins.size()) + energy.Xlo;
135     double y = (static_cast<double>(j) + 0.5) * (energy.Yhi - energy.Ylo)
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 dnisson 1.3 if (fabs(mdef->chisquare_error(nu)) > sig_cut || !ignorezero) {
140 dnisson 1.1 chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
141     / pow(mdef->chisquare_error(nu), 2.0);
142 dnisson 1.3 }
143     else {
144 dnisson 1.1 chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0)
145     / pow(sig_cut, 2.0);
146 dnisson 1.3 }
147 dnisson 1.1 }
148     }
149    
150     // add errors due to Gaussians outside histogram
151 dnisson 1.6 double eps = 0.01; // accuracy set for this function
152 dnisson 1.1 for (int i = 0; i < ngauss; i++) {
153 dnisson 1.3 double *pval = xval+i*(mdef->get_formula()->GetNpar());
154 dnisson 1.1 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]);
159    
160     bool left = par_x < energy.Xlo + cutoff_rad;
161     bool right = par_x > energy.Xhi - cutoff_rad;
162     bool top = par_y > energy.Yhi - cutoff_rad;
163 dnisson 1.2 bool bottom = par_y < energy.Ylo + cutoff_rad;
164 dnisson 1.1
165     if (left) {
166     double xlo = par_x - cutoff_rad;
167     double xhi = energy.Xlo;
168     double ylo = energy.Ylo;
169     double yhi = energy.Yhi;
170     chisquare += formula_int(xlo, xhi, ylo, yhi, pval,
171     XbinSize, YbinSize, xval);
172     }
173     if (right) {
174     double xlo = energy.Xhi;
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,
179     XbinSize, YbinSize, xval);
180     }
181     if (top) {
182     double xlo = left ? par_x - cutoff_rad : energy.Xlo;
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,
187     XbinSize, YbinSize, xval);
188     }
189     if (bottom) {
190     double xlo = left ? par_x - cutoff_rad : energy.Xlo;
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,
195     XbinSize, YbinSize, xval);
196     }
197     }
198    
199     fcnval = chisquare;
200     }
201     catch (exception ex) {
202     cerr << "Exception in jetfit::fcn" << endl;
203     }
204     }
205    
206     bool get_ignorezero() {
207     return ignorezero;
208     }
209    
210     void set_ignorezero(bool _iz) {
211     ignorezero = _iz;
212     }
213    
214     TFormula * model_def::get_formula() {
215     return indiv_formula;
216     }
217    
218     void model_def::set_formula(TFormula *new_formula) {
219     indiv_formula = new_formula;
220     }
221    
222     void model_def::get_indiv_par(int i, string &name, double &start, double &step,
223     double &lo, double &hi) {
224     try {
225     name = indiv_par_names.at(i);
226     start = indiv_par_starts.at(i);
227     step = indiv_par_start_steps.at(i);
228     lo = indiv_par_lo.at(i);
229     hi = indiv_par_hi.at(i);
230     } catch (exception ex) {
231     cerr << "exception in jetfit::model_def::get_indiv_par" << endl;
232     }
233     }
234    
235     void model_def::set_indiv_par(int i, string name, double start, double step,
236     double lo, double hi) {
237     if (indiv_par_names.size() < i+1) {
238     indiv_par_names.resize(i+1);
239     indiv_par_starts.resize(i+1);
240     indiv_par_start_steps.resize(i+1);
241     indiv_par_lo.resize(i+1);
242     indiv_par_hi.resize(i+1);
243     }
244     try {
245     indiv_par_names.at(i) = name;
246     indiv_par_starts.at(i) = start;
247     indiv_par_start_steps.at(i) = step;
248     indiv_par_lo.at(i) = lo;
249     indiv_par_hi.at(i) = hi;
250     }
251     catch (exception ex) {
252     cerr << "exception in jetfit::model_def::set_indiv_par" << endl;
253     }
254     }
255    
256     int model_def::get_indiv_max_E() {
257     return indiv_max_E;
258     }
259    
260     void model_def::set_indiv_max_E(int _new_max_E) {
261     indiv_max_E = _new_max_E;
262     }
263    
264     int model_def::get_indiv_max_x() {
265     return indiv_max_x;
266     }
267    
268     void model_def::set_indiv_max_x(int new_max_x) {
269     indiv_max_x = new_max_x;
270     }
271    
272     int model_def::get_indiv_max_y() {
273     return indiv_max_y;
274     }
275    
276     void model_def::set_indiv_max_y(int new_max_y) {
277     indiv_max_y = new_max_y;
278     }
279    
280     unsigned model_def::get_n_special_par_sets() {
281     return special_par_starts.size();
282     }
283    
284     void model_def::get_special_par(int g, int i, double &pstart,
285     double &perr, double &plo, double &phi) {
286     try {
287     pstart = special_par_starts.at(g).at(i);
288     perr = special_par_start_steps.at(g).at(i);
289     plo = special_par_lo.at(g).at(i);
290     phi = special_par_hi.at(g).at(i);
291     } catch (exception ex) {
292     cerr << "Exception in model_def::get_special_par" << endl;
293     }
294     }
295    
296     void model_def::set_special_par(int g, int i, double pstart,
297     double perr, double plo, double phi) {
298     if (g+1 > special_par_starts.size()) {
299     special_par_starts.resize(g+1);
300     special_par_start_steps.resize(g+1);
301     special_par_lo.resize(g+1);
302     special_par_hi.resize(g+1);
303     }
304    
305     try {
306     if (i+1 > special_par_starts.at(g).size()) {
307     special_par_starts.at(g).resize(i+1);
308     special_par_start_steps.at(g).resize(i+1);
309     special_par_lo.at(g).resize(i+1);
310     special_par_hi.at(g).resize(i+1);
311     }
312    
313     special_par_starts.at(g).at(i) = pstart;
314     special_par_start_steps.at(g).at(i) = perr;
315     special_par_lo.at(g).at(i) = plo;
316     special_par_hi.at(g).at(i) = phi;
317     }
318     catch (exception ex) {
319     cerr << "exception in jetfit::model_def::set_special_par" << endl;
320     }
321     }
322    
323     // translation of CERN's ERFC function
324     double erfc(double x) {
325     bool lef = false;
326     double p2[5], q2[5];
327    
328     const double z1 = 1.0, hf = z1/2.0, c1 = 0.56418958;
329     const double vmax = 7.0;
330     const double switch_ = 4.0;
331    
332     double p10 = 3.6767877, q10 = 3.2584593, p11 = -9.7970465e-2;
333     p2[0] = 7.3738883;
334     q2[0] = 7.3739609;
335     p2[1] = 6.8650185;
336     q2[1] = 1.5184908e1;
337     p2[2] = 3.0317993;
338     q2[2] = 1.2795530e1;
339     p2[3] = 5.6316962e-1;
340     q2[3] = 5.3542168;
341     p2[4] = 4.3187787e-5;
342     q2[4] = 1.0;
343    
344     double p30 = -1.2436854e-1, q30 = 4.4091706e-1, p31 = -9.6821036e-2;
345    
346     double v = fabs(x);
347     double y, h, hc, ap, aq;
348     if (v < hf) {
349     y = v*v;
350     h = x*(p10 + p11*y)/(q10 + y);
351     hc = 1.0 - h;
352     }
353     else {
354     if (v < switch_) {
355     ap = p2[4];
356     aq = q2[4];
357     for (int i = 3; i >= 0; i--) {
358     ap = p2[i] + v*ap;
359     aq = q2[i] + v*aq;
360     }
361     hc = exp(-v*v)*ap/aq;
362     h = 1.0 - hc;
363     }
364     else if (v < vmax) {
365     y = 1/v/v;
366     hc = exp(-v*v)*(c1 + y*(p30 + p31*y)/(q30 + y))/v;
367     h = 1.0 - hc;
368     }
369     // for very big values we can save us any calculation
370     // and the FP-exceptions from exp
371     else {
372     h = 1.0;
373     hc = 0.0;
374     }
375     if (x < 0) {
376     h = -h;
377     hc = 2.0 - hc;
378     }
379     }
380     if (lef) {
381     return h;
382     }
383     else {
384     return hc;
385     }
386     }
387    
388     // translation of CERN's PROB function
389     double prob(double x, int n) {
390     const char *name = "PROB";
391     char errtxt[80];
392    
393     const double r1 = 1.0,
394     hf = r1/2.0, th = r1/3.0, f1 = 2.0*r1/9.0;
395     const double c1 = 1.128379167095513;
396     const int nmax = 300;
397     // maximum chi2 per df for df >= 2., if chi2/df > chipdf prob=0
398     const double chipdf = 100.0;
399     const double xmax = 174.673, xmax2 = 2.0*xmax;
400     const double xlim = 24.0;
401     const double eps = 1e-30;
402    
403     double y = x;
404     double u = hf*y;
405     double h, w, s, t, fi, e;
406     int m;
407     if (n < 0) {
408     h = 0.0;
409     sprintf(errtxt, "n = %d < 1", n);
410     cerr << "PROB: G100.1: "<<errtxt;
411     }
412     else if (y < 0.0) {
413     h = 0.0;
414     sprintf(errtxt, "x = %f < 0", n);
415     cerr << "PROB: G100.2: "<<errtxt;
416     }
417     else if (y == 0.0 || n/20 > y) {
418     h = 1.0;
419     }
420     else if (n == 1) {
421     w = sqrt(u);
422     if (w < xlim) {
423     h = erfc(w);
424     }
425     else {
426     h = 0.0;
427     }
428     }
429     else if (n > nmax) {
430     s = r1 / ((double)n);
431     t = f1 * s;
432     w = (pow(y*s, th) - (1.0 - t)) / sqrt(2.0*t);
433     if (w < -xlim) {
434     h = 1.0;
435     }
436     else if (w < xlim) {
437     h = hf * erfc(w);
438     }
439     else {
440     h = 0.0;
441     }
442     }
443     else {
444     m = n/2;
445     if (u < xmax2 && (y/n) <= chipdf) {
446     s = exp(-hf*u);
447     t = s;
448     e = s;
449     if (2*m == n) {
450     fi = 0.0;
451     for (int i = 1; i < m; i++) {
452     fi += 1.0;
453     t = u*t/fi;
454     s += t;
455     }
456     h = s*e;
457     }
458     else {
459     fi = 1.0;
460     for (int i = 1; i < m; i++) {
461     fi += 2.0;
462     t = t*y/fi;
463     s += t;
464     }
465     w = sqrt(u);
466     if (w < xlim) {
467     h = c1*w*s*e + erfc(w);
468     }
469     else {
470     h = 0.0;
471     }
472     }
473     }
474     else {
475     h = 0.0;
476     }
477     }
478     if (h > eps) {
479     return h;
480     }
481     else {
482     return 0.0;
483     }
484     }
485    
486 dnisson 1.5 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 dnisson 1.6 hist->SetBinContent(i+1, j+1, vec[i][j]);
490 dnisson 1.5 }
491     }
492     }
493    
494 dnisson 1.1 void par_init(TH2 *hist, TMinuit *gMinuit, vector<TString> &pars,
495     double *pval, double *perr, double *plo, double *phi,
496     int npar, results r) {
497     int npar1 = npar - mdef->get_formula()->GetNpar();
498     bool init_spec_pars = false;
499     if (ngauss <= mdef->get_n_special_par_sets()) {
500     init_spec_pars = true;
501     for (int k = 0; k < ngauss; k++) {
502 dnisson 1.6 int ipar0 = k*mdef->get_formula()->GetNpar(); // index of indiv par 0
503 dnisson 1.1 for (int l = 0; l < mdef->get_formula()->GetNpar(); l++) {
504     int ipar = ipar0 + l;
505     if (ipar < npar)
506     mdef->get_special_par(k, l, pval[ipar], perr[ipar], plo[ipar],
507     phi[ipar]);
508     else
509     cerr << "WARNING: Attempt to set parameter out of index range!"
510     << endl;
511     }
512     }
513     }
514     else {
515     double XbinSize = (energy.Xhi - energy.Xlo)
516     / static_cast<double>(energy.bins.size());
517     double YbinSize = (energy.Xhi - energy.Xlo)
518     / static_cast<double>(energy.bins.at(0).size());
519     vector< vector<double> > sub_energy(energy.bins.size());
520     double max_sub, max_x = 0.0, max_y = 0.0;
521     double pval_other[256];
522 dnisson 1.3 max_sub = -(numeric_limits<double>::max());
523 dnisson 1.1 for (int i = 0; i < energy.bins.size(); i++) {
524     sub_energy[i].resize(energy.bins[i].size());
525     for (int j = 0; j < energy.bins[i].size(); j++) {
526     double x = static_cast<double>(i)*(energy.Xhi - energy.Xlo)
527     / static_cast<double>(energy.bins.size()) + energy.Xlo;
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 dnisson 1.6 // subtract integral of fit function
532 dnisson 1.1 try {
533     int npval_other = r.pval.at(r.pval.size()-1).size();
534     if (npval_other > 256) {
535     cerr << "Parameter overload" << endl;
536     return;
537     }
538     for (int k = 0; k < npval_other; k++) {
539     pval_other[k] = r.pval[r.pval.size()-1][k];
540     }
541     ngauss--;
542     double nu = fit_fcn(x, y, pval_other) * XbinSize * YbinSize;
543     ngauss++;
544 dnisson 1.6 sub_energy[i][j] = energy.bins[i][j] - nu;
545 dnisson 1.1 }
546     catch (exception ex) {
547     cerr << "Exception in par_init" << endl;
548     }
549     }
550     else {
551     sub_energy[i][j] = energy.bins[i][j];
552     }
553     }
554     }
555 dnisson 1.5 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 = hist_sub->GetSumOfWeights();
562     max_x = hist_sub->GetMean(1);
563     max_y = hist_sub->GetMean(2);
564 dnisson 1.1
565     for (int i = npar1; i < npar; i++) {
566     if (mdef->get_indiv_max_E() == i - npar1) {
567     double nu = 0.0;
568     if (ngauss > 1) {
569 dnisson 1.3 ngauss--;
570 dnisson 1.1 nu = fit_fcn(max_x, max_y, pval_other) * XbinSize * YbinSize;
571 dnisson 1.3 ngauss++;
572 dnisson 1.1 }
573 dnisson 1.5 pval[i] = max_sub;
574 dnisson 1.4 perr[i] = mdef->chisquare_error(pval[i])*0.5;
575 dnisson 1.1 plo[i] = 0.0;
576     phi[i] = 1.0e6;
577     }
578     else if (mdef->get_indiv_max_x() == i - npar1) {
579     pval[i] = max_x;
580     perr[i] = (energy.Xhi - energy.Xlo) / static_cast<double>(hist->GetXaxis()->GetNbins());
581     plo[i] = 0.0;
582     phi[i] = 0.0;
583     }
584     else if (mdef->get_indiv_max_y() == i - npar1) {
585     pval[i] = max_y;
586     perr[i] = (energy.Yhi - energy.Ylo) / static_cast<double>(hist->GetYaxis()->GetNbins());
587     plo[i] = 0.0;
588     phi[i] = 0.0;
589     }
590     else {
591     string dummy;
592     mdef->get_indiv_par(i-npar1, dummy, pval[i], perr[i], plo[i], phi[i]);
593     }
594     }
595     }
596     int n_pars_to_set = init_spec_pars ? ngauss*mdef->get_formula()->GetNpar()
597     : npar - npar1;
598     int init_par_to_set = init_spec_pars ? 0 : npar1;
599     for (int i = 0; i < n_pars_to_set; i++) {
600     double dummy;
601     string s;
602     int ipar = i + init_par_to_set;
603     mdef->get_indiv_par(i % mdef->get_formula()->GetNpar(), s,
604     dummy, dummy, dummy, dummy);
605     ostringstream par_oss;
606     par_oss << s << ipar / (mdef->get_formula()->GetNpar()) << flush;
607     try {
608     pars.at(ipar) = TString(par_oss.str());
609     }
610     catch (exception ex) {
611     cerr << "exception 2 in par_init" << endl;
612     }
613     int error_flag;
614     try {
615     gMinuit->mnparm(ipar, pars.at(ipar),
616     pval[ipar], perr[ipar], plo[ipar], phi[ipar], error_flag);
617     } catch (exception ex) {
618     cerr << "exception 3 in par_init" << endl;
619     }
620     }
621     }
622    
623     results fit_histo(TH2 *hist, vector<trouble> &t_vec,
624     void (*cc_minuit)(TMinuit *, TH2 *, int),
625     int start_ngauss,
626     int rebinX, int rebinY,
627     double P_cutoff_val) {
628     TMinuit *gMinuit = new TMinuit();
629     int npar_indiv = mdef->get_formula()->GetNpar();
630     int istat, erflg;
631     results r;
632    
633     gMinuit->SetFCN(fcn);
634     gMinuit->mninit(5, 6, 7);
635    
636     // set error def and machine accuracy
637     double cs_err_def = 1.0;
638     double fcn_eps = 0.2;
639     gMinuit->mnexcm("SET ERR", &cs_err_def, 1, erflg);
640     gMinuit->mnexcm("SET EPS", &fcn_eps, 1, erflg);
641    
642     // rebin histogram
643     TH2 *hist_rebinned;
644     if (rebinX > 1 && rebinY > 1)
645     hist_rebinned = hist->Rebin2D(rebinX, rebinY);
646     else
647     hist_rebinned = hist;
648    
649     r.hist_rebinned = hist_rebinned;
650     // load histogram into energies vector
651     energy.bins.resize(hist_rebinned->GetXaxis()->GetNbins());
652     for (int i = 0; i < energy.bins.size(); i++) {
653     energy.bins[i].resize(hist_rebinned->GetYaxis()->GetNbins());
654     for (int j = 0; j < energy.bins[i].size(); j++) {
655     energy.bins[i][j] = hist_rebinned->GetBinContent(i+1, j+1);
656     }
657     }
658    
659     energy.Xlo = hist->GetXaxis()->GetXmin();
660     energy.Xhi = hist->GetXaxis()->GetXmax();
661     energy.Ylo = hist->GetYaxis()->GetXmin();
662     energy.Yhi = hist->GetYaxis()->GetXmax();
663    
664     if (start_ngauss < 1) {
665     start_ngauss = 1;
666     }
667    
668     for (ngauss = start_ngauss;
669     ngauss <= (MAX_GAUSS > start_ngauss ? MAX_GAUSS : start_ngauss);
670     ngauss++) {
671     t_vec.resize(t_vec.size() + 1);
672     int npar = npar_indiv*ngauss;
673     double pval[256], perr[256], plo[256], phi[256];
674     if (npar > 256) {
675     cerr << "Parameter overload" << endl;
676     return r;
677     }
678     vector<TString> pars(npar);
679    
680     // initialize parameters
681     par_init(hist, gMinuit, pars, pval, perr, plo, phi, npar, r);
682 dnisson 1.5
683 dnisson 1.1 // minimize
684     double chisquare, edm, errdef;
685     int nvpar, nparx;
686     trouble t;
687     t.occ = T_NULL;
688     t.istat = 3;
689 dnisson 1.4 // fix sigmas of the Gaussians
690 dnisson 1.1 for (int i = 0; i < ngauss; i++) {
691     double parno[2];
692 dnisson 1.4 parno[0] = i*npar_indiv + 4;
693     gMinuit->mnexcm("FIX", parno, 1, erflg);
694 dnisson 1.1 }
695     gMinuit->mnexcm("SIMPLEX", 0, 0, erflg);
696     gMinuit->mnexcm("MIGRAD", 0, 0, erflg);
697     gMinuit->mnstat(chisquare, edm, errdef, nvpar, nparx, istat);
698     if (istat == 3) {
699     /* we're not concerned about MINOS errors right now
700     gMinuit->mnexcm("MINOS", 0, 0, erflg);
701     gMinuit->mnstat(chisquare, edm, errdef, nvpar, nparx, istat);
702     if (istat != 3) {
703     t.occ = T_MINOS;
704     t.istat = istat;
705     }
706     */
707     }
708     else {
709     t.occ = T_MIGRAD;
710     t.istat = istat;
711     }
712    
713     try {
714     if (t_vec.size() < ngauss) {
715     t_vec.resize(ngauss);
716     }
717     t_vec.at(ngauss-1) = t;
718     }
719     catch (exception ex) {
720     cerr << "exception in fit_histo" << endl;
721     }
722    
723     // put parameters in map
724     for (int i = 0; i < npar && i < 256; i++) {
725     int iuint; // internal parameter number
726     gMinuit->mnpout(i, pars[i], pval[i], perr[i], plo[i], phi[i], iuint);
727     }
728     vector<double> pval_copy(npar);
729     vector<double> perr_copy(npar);
730     vector<double> plo_copy(npar);
731     vector<double> phi_copy(npar);
732     for (int i = 0; i < npar && i < 256; i++) {
733     pval_copy[i] = pval[i];
734     perr_copy[i] = perr[i];
735     plo_copy[i] = plo[i];
736     phi_copy[i] = phi[i];
737     }
738     r.pars.push_back(pars);
739     r.pval.push_back(pval_copy);
740     r.perr.push_back(perr_copy);
741     r.plo.push_back(plo_copy);
742     r.phi.push_back(phi_copy);
743    
744     // execute user minuit code
745     if (cc_minuit != 0)
746     (*cc_minuit)(gMinuit, hist_rebinned, ngauss);
747    
748     int ndof = 0;
749     if (!ignorezero)
750     ndof = energy.bins.size() * energy.bins[0].size();
751     else {
752     for (int i = 0; i < energy.bins.size(); i++) {
753     for (int j = 0; j < energy.bins[i].size(); j++) {
754     if (energy.bins[i][j] > 1.0e-30)
755     ndof++;
756     }
757     }
758     }
759     ndof -= npar;
760    
761 dnisson 1.2 r.chisquare.push_back(chisquare);
762 dnisson 1.1 double P = prob(chisquare, ndof);
763     if (P > P_cutoff_val) {
764     break;
765     }
766     }
767    
768     delete gMinuit;
769     return r;
770     }
771     }