Analysis of the example 'fitConvolution.C' in Fit Tutorial


I’m not sure that I am writing to the right place, but I hope you
correct me.

The example ‘fitConvolution.C’ in the fitting tutorial is very vague
and, moreover, it seems to me contains an error.
The result of the fit in the example is very different from the
expected values: (0.073 +/- 0.002) instead of -0.3, (-2.26 +/- 0.05)
instead of 0 and 1.13 +/- 0.06 instead of 3. At the same time, the
goodness of the fit is very poor: chi^2/ndf = 298/96.

First, the example uses the convolution of the functions specified
through TFormula: the ctr TF1Convolution is given through the
functions “expo” and “gaus”.
In my opinion, this is an outdated and extremely unfortunate
interface that contributes to the appearance of errors. It is
completely unclear how the normalization constants of each of these
functions are specified. But okay, let’s assume that it gives
convolution of exp(a+bx) (+) exp(-0.5((x-m)/s)^2). Here (+) denotes
the convolution given by the integral from minus infinity to plus
infinity. It is easy to show that by changing the variables in the
integral, the expression is transformed to the form:
exp((a+bm) + bx) (+) exp(-0.5((x/s)^2).
This shows that the parameters a and m are fully correlated.

The second part of the problem is that the Monte Carlo simulation
uses gRandom->Exp() that generates values only for positive values.
So instead of exp(a + b*x) we have to use a function that is zero
for negative x.

Below is a modification of the example:

void fitConvolution() {
   gStyle -> SetOptStat(1111111); //name, entrs, mean, rms, under, over, int
   gStyle -> SetOptFit(112); // print fixed parameters

   gRandom -> SetSeed(65539); // for repeatable results

   // histograms to fit
   TH1F* h_Exp = new TH1F("h_Exp", "Exp", 100,-1.,10.);
   TH1F* h_Gauss = new TH1F("h_Gauss", "Gauss", 100,-9.,9.);
   TH1F* h_ExpGauss = new TH1F("h_ExpGauss",
         "Exponential convoluted by gaussian", 100,0.,5.);
   for (int i=0;i<1e6;i++) {
      // 1/0.3 gives a alpha of -0.3 in the exp
      Double_t x = gRandom->Exp(1./0.3);

      // probability density function of the addition of two
      // variables is the convolution of two density functions
//       double y = gRandom->Gaus(0.,3.);
      double y = gRandom->Gaus(0.,2.);

   // clone histogram for second fit
   TH1F* h_ExpGaussC = (TH1F*)h_ExpGauss->Clone("h_ExpGaussC");
//    h_ExpGaussC -> SetTitle("example with lambda functions");

   // 1) original example from root site:

   bool useFFT = false;
//    bool useFFT = true;

   TF1Convolution *f_conv = new TF1Convolution("expo","gaus",-1,6,useFFT);
   TF1   *f = new TF1("f",*f_conv, 0., 5., f_conv->GetNpar());

   // 2) the same example modified using lambda functions
   auto Lexp = [](const double* x,const double* p) -> double {
//       return exp(p[0] + x[0]*p[1]);
      return (x[0]>=0) ? exp(p[0] + x[0]*p[1]) : 0;
   TF1* tex = new TF1("tex",Lexp,-1.,6.,2);

   auto Lgs = [](const double* x,const double* p) -> double {
      return TMath::Gaus(x[0],p[0],p[1],false);
   TF1* tgs = new TF1("tgs",Lgs,-15.,15.,2);

   TF1Convolution* conv_eg = new TF1Convolution(tex,tgs,useFFT);
   int Npar = conv_eg -> GetNpar();
//    cout << " Npar= " << Npar << endl;
   conv_eg -> SetRange(-1.,6.);
   conv_eg -> SetNofPointsFFT(10000);

   TF1* tconv_eg = new TF1("tconv_eg", *conv_eg, 0., 5., Npar);

   // Commented out lines crashes CLING:
   //  -> cling_runtime_internal_throwIfInvalidPointer
   // Man got to sit and wonder ‘why, why, why?’
//    auto Lconv = [&conv_eg](const double* x,const double* p) -> double {
//       return conv_eg -> operator()(x,p);
//    };
//    TF1* tconv_eg = new TF1("tconv_eg", Lconv, 0., 5., Npar);

   tconv_eg -> SetParameters( 1.,-0.3, 0., 2.);
   tconv_eg -> FixParameter( 2, 0. );

   // fit and draw
   auto c = new TCanvas("c","...",0,0,1200,800);
   c -> Divide(2,1);
   auto pad = c -> cd(1);
   pad -> Divide(1,2);
   pad -> cd(1);
   h_Exp -> SetLineWidth(2);
   h_Exp -> Fit("expo");
   h_Exp -> Draw("EP");
   pad -> cd(2);
   h_Gauss -> SetLineWidth(2);
   h_Gauss -> Fit("gaus");
   h_Gauss -> Draw("EP");

   c -> cd(2);
//    h_ExpGauss -> Draw("EP");
//    h_ExpGauss -> Fit("f","E");
   h_ExpGaussC -> Fit("tconv_eg","E");
   h_ExpGaussC -> Draw("EP");

Better agreement with expected parameter values, but there is some
shift. Apparently, this is due to a very wide Gaussian compared to
the fitting range, which leads to significant smoothing and leakage
of events outside the fitting range.
Therefore, for this example, I would suggest reducing the sigma of
Gaussian to two.

I hope it will be use useful.

Please read tips for efficient and successful posting and posting code

_ROOT Version: v6.24.00-gcc10.2
Platform: x86_64 GNU/Linux
Compiler: building CLING

Thanks. I’m sure @moneta will comment on this

One more comment on the reasons why the fit gives biased values for
the wide Gaussian case. My last point in previous letter was wrong.
The reason lies in the wrong choice of the range of convolution. For
the case of such a wide Gaussian (sigma=3) the offset by one is
unreasonably small. The right border of this range should be at least
two sigmas from the fitting area. The left border can be chosen very
close to zero due to the fact that the first function vanishes at
negative values.

Therefore the following modification of the above example gives
a perfect fit for sigma=3.

 conv_eg -> SetRange(-0.1,11.);

The values obtained in fit:

 FCN=102.795 FROM MINOS     STATUS=SUCCESSFUL     45 CALLS         323 TOTAL
                     EDM=5.38663e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           7.59554e+00   7.24513e-03   1.57851e-05  -2.15537e-01
   2  p1          -2.97369e-01   5.11931e-03  -4.92313e-05  -5.10843e-01
   3  p2           0.00000e+00     fixed    
   4  p3           2.97040e+00   5.49297e-02   5.49297e-02  -8.30937e-03