| 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 = mdef->chisquare_error(nu) < sig_cut |
| 109 |
> |
double sigma = fabs(mdef->chisquare_error(nu)) < sig_cut |
| 110 |
|
? sig_cut : mdef->chisquare_error(nu); |
| 111 |
|
|
| 112 |
|
fsum += xstep * ystep * mdef->get_formula()->Eval(x, y) |
| 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 |
< |
double x = (static_cast<double>(i) + 0.5) * (energy.Xhi - energy.Xlo) |
| 133 |
> |
double x = (static_cast<double>(i) + 0.5) * (energy.Xhi - energy.Xlo) |
| 134 |
|
/ 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 |
< |
if (mdef->chisquare_error(nu) > sig_cut || !ignorezero) |
| 139 |
> |
if (fabs(mdef->chisquare_error(nu)) > sig_cut || !ignorezero) { |
| 140 |
|
chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0) |
| 141 |
|
/ pow(mdef->chisquare_error(nu), 2.0); |
| 142 |
< |
else |
| 142 |
> |
} |
| 143 |
> |
else { |
| 144 |
|
chisquare += pow(energy.bins.at(i).at(j) - nu, 2.0) |
| 145 |
|
/ pow(sig_cut, 2.0); |
| 146 |
+ |
} |
| 147 |
|
} |
| 148 |
|
} |
| 149 |
|
|
| 150 |
|
// add errors due to Gaussians outside histogram |
| 151 |
|
double eps = 0.1; // accuracy set for this function |
| 152 |
|
for (int i = 0; i < ngauss; i++) { |
| 153 |
< |
double *pval = xval+i*mdef->get_formula()->GetNpar(); |
| 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]; |
| 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 |
< |
bool bottom = par_y > energy.Ylo + cutoff_rad; |
| 163 |
> |
bool bottom = par_y < energy.Ylo + cutoff_rad; |
| 164 |
|
|
| 165 |
|
if (left) { |
| 166 |
|
double xlo = par_x - cutoff_rad; |
| 511 |
|
vector< vector<double> > sub_energy(energy.bins.size()); |
| 512 |
|
double max_sub, max_x = 0.0, max_y = 0.0; |
| 513 |
|
double pval_other[256]; |
| 514 |
< |
max_sub = -(numeric_limits<double>::infinity()); |
| 514 |
> |
max_sub = -(numeric_limits<double>::max()); |
| 515 |
|
for (int i = 0; i < energy.bins.size(); i++) { |
| 516 |
|
sub_energy[i].resize(energy.bins[i].size()); |
| 517 |
|
for (int j = 0; j < energy.bins[i].size(); j++) { |
| 534 |
|
double nu = fit_fcn(x, y, pval_other) * XbinSize * YbinSize; |
| 535 |
|
ngauss++; |
| 536 |
|
sub_energy[i][j] = energy.bins[i][j] - nu |
| 537 |
< |
- 4.0*mdef->chisquare_error(nu); |
| 537 |
> |
- 2.2*mdef->chisquare_error(nu); |
| 538 |
|
} |
| 539 |
|
catch (exception ex) { |
| 540 |
|
cerr << "Exception in par_init" << endl; |
| 555 |
|
if (mdef->get_indiv_max_E() == i - npar1) { |
| 556 |
|
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 ? 4.0*mdef->chisquare_error(nu) |
| 563 |
< |
: 0.0); |
| 564 |
< |
perr[i] = mdef->chisquare_error(pval[i])*0.1; |
| 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; |
| 566 |
> |
perr[i] = mdef->chisquare_error(pval[i])*0.5; |
| 567 |
|
plo[i] = 0.0; |
| 568 |
|
phi[i] = 1.0e6; |
| 569 |
|
} |
| 678 |
|
trouble t; |
| 679 |
|
t.occ = T_NULL; |
| 680 |
|
t.istat = 3; |
| 681 |
< |
// fix the N values and sigmas of the Gaussians |
| 681 |
> |
// fix sigmas of the Gaussians |
| 682 |
|
for (int i = 0; i < ngauss; i++) { |
| 683 |
|
double parno[2]; |
| 684 |
< |
parno[0] = i*npar_indiv + 1; |
| 685 |
< |
parno[1] = i*npar_indiv + 4; |
| 680 |
< |
gMinuit->mnexcm("FIX", parno, 2, erflg); |
| 684 |
> |
parno[0] = i*npar_indiv + 4; |
| 685 |
> |
gMinuit->mnexcm("FIX", parno, 1, erflg); |
| 686 |
|
} |
| 687 |
|
gMinuit->mnexcm("SIMPLEX", 0, 0, erflg); |
| 688 |
|
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); |
| 689 |
|
gMinuit->mnstat(chisquare, edm, errdef, nvpar, nparx, istat); |
| 690 |
|
if (istat == 3) { |
| 691 |
|
/* we're not concerned about MINOS errors right now |
| 750 |
|
} |
| 751 |
|
ndof -= npar; |
| 752 |
|
|
| 753 |
+ |
r.chisquare.push_back(chisquare); |
| 754 |
|
double P = prob(chisquare, ndof); |
| 770 |
– |
cout << "P = "<<P << endl; |
| 755 |
|
if (P > P_cutoff_val) { |
| 756 |
|
break; |
| 757 |
|
} |
