Switching a fit object from TLinearFitter to TMinuit2

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!

I think @moneta can help you.

Hi,
I don’t understand why you want to move from TLinearFitter to TMinuit2. If you are fitting a linear function, it is better using the linear fitter. Otherwise, you should use either TGraph::Fit or using the ROOT::Math::Fitter class, which works for linear and non-linear fitter. The class TMinui2 is also deprecated and I recommend to not use it.
When using the Fitter class you either provide the model function (the one you provide as a formula in TLinearFitter) or the objective function (e.g. the chi-square function). If you want to do multiple fits adding fitting points, you can do adding points to the ROOT::Fit::BinData class.

If you need help in this conversion, please let me know

Cheers

Lorenzo

Hi, Lorenzo,

Thanks a lot for your response! I want to switch from using TLinearFitter to using another fitter because I am going to add other model function options such as the Chebyshev polynomials and some nonlinear user-defined functions, and it seems to me that the TLinearFitter is a little bit limited on that. Would you mind telling me how you would do to use another fitter that can work with other more complex model functions in my case or maybe giving me an example?

Thank you so much!

Hi,

Here is an example, similar to your code, on how to use the Fitter class.

Cheers

Lorenzo

exampleFitter.C (2.2 KB)

Hi, Lorenzo,

Thanks so much for the example, and it really helped me a lot to understand! I’ve changed the fitter from TLinearFitter to ROOT::Fit::Fitter, but it takes so much longer to run the program now, I know the TLinearFitter is supposed to be faster, but now it’s super slow to do the fit with the ROOT::Fit::Fitter, is it normal or maybe I am doing something wrong here?

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();
      TF1 f1("f1", formula[order]);
      ROOT::Fit::Fitter fitter;
      ROOT::Math::WrappedMultiTF1 modelFunc(f1, 1);
      fitter.SetFunction(modelFunc);
      ROOT::Fit::BinData data(scanSize(), f1.GetNdim(), ROOT::Fit::BinData::kNoError);

      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);
            data.Add(&adc, v);
            jmax = j;
          }
          last_adc = adc;

          if (!jmin) jmin = j;
        }
      }

      // Set minimizer
      fitter.Config().SetMinimizer("Minuit2");

      if (vref)
      {
        //fit.FixParameter(0, 0);
        fitter.Config().ParSettings(0).SetValue(0.);
        fitter.Config().ParSettings(0).Fix();
      }
      if (nzero > 1)
      {
        nbroken++;
      }
      else
      {
        //fit.Eval();
        bool ret = fitter.Fit(data);
      }


      //if (vref) fit.ReleaseParameter(0);
      if (vref) fitter.Config().ParSettings(0).Release();

      auto &result = fitter.Result();

      for (int iorder = 0; iorder <= fit_order; iorder++)
      {
        //fit_coeffs[ichan][i * (order+1) + iorder] = fit.GetParameter(iorder);
        fit_coeffs[ichan][i * (order+1) + iorder] = result.Parameter(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);
  }
}

Thank you!

Hi,
It is expected that Minuit2, which is an iterative minimisation, is slower than the linear fitter, but the difference should not be too large. To have an idea what is the size of num_radiant_channels, num_lab4_samples, (i.e. the number of fits you are doing) and scanSize (the size of the data ) ?

You can also try to using (with the ROOT::Fit::Fitter) the linear fitter by doing:

fitter.Config().SetMinimizer("Linear");

This should give a time similar to using the TLinearFitter.

Cheers

Lorenzo

The size of num_radiant_channels is 24, the size of num_lab4_samples is 4096, and the size of scanSize() is 155.

I tried fitter.Config().SetMinimizer("Linear"), and the fitting didn’t seem completely right, it gave me an error because I fixed one of the parameters to zero:
ROOT::Math::FitResult: FitConfiguration and Minimizer result are not consistent
Number of free parameters from FitConfig = 5
Number of free parameters from Minimizer = 6

Is there any way that I can fix this and make it work?

Thank you!

Another question I have is that I want to apply a triangle wave function to the fitting, because there seems to be some periodic behaviors in my data. I wrote a single program to deal with just one graph to see if I can do the fitting well, but the result is not ideal. Could you please also give me some suggestions on this? Also, if I want to use my C++ triangle wave function on different sub-ranges and then combine them to do the fit, because in different sections of the data, the periods are not the same for the triangle waves, is there a good way to do this? I am sorry that I am asking a lot of questions at once.

newFitter.C (3.7 KB)

sampleDataPoints.root (6.9 KB)

Thank you!

Hi,

Concerning your first question, this is a problem I have discovered last week after your post. It is now fixed in the master version. The problem is caused when fixing the parameters when using the Linear Fitter. I think if you are fixing the first parameter (for ipar=0), then you can ignore that error message.

Since you are doing many fits, it is possible to see a slowdown when using a non-linear fitter.

I will look now at your second question,

Cheers

Lorenzo

Hi,

The fit is not ideal probably because the data are different than the function. I see you have fixed the parameter of the triage function, maybe these parameter are not optimal.
If you want to change depending on the data, you can either perform a combined fit as in this example tutorial, ROOT: tutorials/fit/combinedFit.C File Reference, where you consider separate data sets, or you could also defining the function differently, where you apply different parameter depending on the data range.

Lorenzo

Hi, Lorenzo,

Thanks a lot for answering all my questions. This might be my last question in this post, why is the object ROOT::Fit::Fitter fitter cannot be declared outside of the loop? Because it seems like it has to be called every time for each fitting, and I know it wouldn’t work if it was outside of the loop, but I want to know the reason behind it.

Thank you!

Hi,
You can declare the Fitter object outside the loop. When you are setting a new function and new data it should automatically reset and use the new information. Creating the Fitter class is a rather light operation, so you should not gain very much by declaring it outside the loop.
Cheers
Lorenzo

Hi, Lorenzo,

I think I got confused by the automatic output information from the minimizer itself, if the fitter is declared outside of the loop, many parameters will show fixed and not show the errors, but if we look at the parameters and the errors by using result.Parameter(ipar) and result.ParError(ipar), the fitting works the same in both ways.

You’ve answered all my questions and I really appreciate it.

Thank you so much!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.