Hi, ROOT experts,
I am trying to switch from using TLinearFitter to using TMinuit2 in this piece of code:
static const char * formula[1+max_voltage_calibration_fit_order] = {"pol0","pol1","pol2","pol3","pol4","pol5","pol6","pol7","pol8","pol9"};
void recalculateFits(int order, double min, double max, double vref, uint32_t mask, int turnover_threshold)
{
fit_order = order < max_voltage_calibration_fit_order ? order : max_voltage_calibration_fit_order;
order = fit_order;
fit_min = min;
this->turnover_threshold = turnover_threshold;
fit_max = max;
fit_vref = vref;
if ( min > max)
{
fit_min = max;
fit_max = min;
}
TLinearFitter fit(1, formula[order]);
for (int ichan = 0; ichan < num_radiant_channels; ichan++)
{
fit_coeffs[ichan].clear();
if ( (mask & (1 << ichan)) == 0) continue;
//check if channel seems broken by looking for zeroes
fit_coeffs[ichan].resize((order+1) * num_lab4_samples, 0);
int nbroken = 0;
for (int i = 0; i < num_lab4_samples; I++)
{
fit.ClearPoints();
int jmin = 0;
int jmax = 0;
if (vref)
{
double last_v = 0;
for (unsigned j = 0; j < scanSize(); j++)
{
double v = ichan < num_radiant_channels / 2 ? vbias[0][j]: vbias[1][j];
if (v >= vref)
{
if (v == vref || j == 0)
{
adc_offset[ichan] = scan_result[j][ichan][I];
}
else //have too interpolate
{
double frac_low = (v - vref) / (v-last_v);
adc_offset[ichan] = scan_result[j][ichan][i] * (1-frac_low) + frac_low * scan_result[j-1][ichan][I];
}
break;
}
last_v = v;
}
}
double last_adc = -4096;
int nzero = 0;
for (unsigned j = 0; j < scanSize(); j++)
{
double v = ichan < num_radiant_channels / 2 ? vbias[0][j]: vbias[1][j];
if (v >= fit_min && v <= fit_max)
{
if (vref)
{
v-=vref;
}
double adc = scan_result[j][ichan][i] - adc_offset[ichan];
if (adc == 0) nzero++;
//don't include if it turns over
if (!jmin || adc > last_adc - turnover_threshold)
{
fit.AddPoint(&adc, v);
jmax = j;
}
last_adc = adc;
if (!jmin) jmin = j;
}
}
if (vref)
{
fit.FixParameter(0, 0);
}
if (nzero > 1)
{
nbroken++;
}
else
{
fit.Eval();
}
if (vref) fit.ReleaseParameter(0);
for (int iorder = 0; iorder <= fit_order; iorder++)
{
fit_coeffs[ichan][i * (order+1) + iorder] = fit.GetParameter(iorder);
}
fit_chisq[ichan][i] = 0;
fit_maxerr[ichan][i] = 0;
fit_ndof[ichan][i] = (jmax - jmin) + 1; // number of points including start and end
fit_ndof[ichan][i] -= (fit_order + 1); // subtract number of parameters of the fit function
if (vref) fit_ndof[ichan][i] += 1; // if parameter 0 is fixed, add back one degree of freedom
turnover_index[ichan][i] = jmax+1;
//calculate max deviation and chi square
for (int j = jmin ; j <= jmax; j++)
{
double adc = scan_result[j][ichan][i] - adc_offset[ichan];
double meas = ichan < num_radiant_channels / 2 ? vbias[0][j]: vbias[1][j];
meas -= vref;
double pred = evalPars(adc, fit_order, &fit_coeffs[ichan][i * (order+1)]);
double delta = fabs(meas-pred);
if (delta > fit_maxerr[ichan][i]) fit_maxerr[ichan][i] = delta;
fit_chisq[ichan][i] +=delta*delta / (0.010*0.010); // estimate 10 mV equivalent RMS?
}
}
if (nbroken) printf("WARNING: Channel %d seems to have %d broken samples?\n", ichan, nbroken);
}
}
Here are a couple of questions:
(1) Is there a similar way to call out a fit object for using TMinuit2 like the one using TLinearFitter in the code above?
TLinearFitter fit(1, formula)
(2) If TMinuit2 had been used, how could one clear the fitting points in a graph for preparing to do the fitting for the next graph when there were many different graphs that need to be fitted? In the TLinearFitter case there’s the fit.ClearPoints()
.
(3) The function AddPoint()
can be used with the TLinearFitter object but is this feature can also be used for other fitters? For example, if a TMinuit2 fit object was called, would there be something like this that can be used to add fitting points?
fit.AddPoint(&adc, v)
(4) I am confused with the ReleaseParameter()
, could you please tell me when to use it?
I know there are some tutorials about TMinuit2, but the ways it is used in those tutorials are more straightforward to me such as using ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2")
, and then a TF1 function can be used to fit a graph directly; however, I don’t know how to call out a fit object for TMinuit2 like we do for TLinearFitter ( TLinearFitter fit(1, formula)
), and I am really lost in this project, so I hope I can get some help from the experts.
Thank you!