Hello,
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);
h_Exp->Fill(x);
// 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.);
h_Gauss->Fill(y);
h_ExpGauss->Fill(x+y);
}
// 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:
// https://root.cern/doc/master/fitConvolution_8C.html
bool useFFT = false;
// bool useFFT = true;
TF1Convolution *f_conv = new TF1Convolution("expo","gaus",-1,6,useFFT);
f_conv->SetRange(-1.,6.);
f_conv->SetNofPointsFFT(10000);
TF1 *f = new TF1("f",*f_conv, 0., 5., f_conv->GetNpar());
f->SetParameters(1.,-0.3,0.,3.);
// 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.
Yury.
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