Create a function based on a simulation

Hello Experts,
I am using Spline to create a function based on a simulation data. Then I use this function to fit my real data.
The function I create has 2 parameters (amplitude for signal and amplitude for background).

Now the fitting function looks wired and it is not smooth (see attached). I am wondering if there is a way to improve it ? OR if there is a better method that I can use to make a function based on simulation then use it to fit my data.
Also I am interest on yield for both signal and background with their uncertainties. How I can get them correctly.

I have also attached my code and root files ( sim and data).

ndeFitsexp.C (1.0 KB)
mySplineFit.C (967 Bytes)
data.root (52.9 KB)
sim.root (140.3 KB)


Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided


Try: f1->SetNpx(3 * histname->GetNbinsX()); // e.g., 3 or 5

Thanks @Wile_E_Coyote for your reply, now it looks better. Can help me on a correct way to get a signal yield with uncertainty.

If you just want the integral of your total “f1” ("myfuncTest10") function, you can simply use the TF1::Integral and TF1::IntegralError methods.

If you want integrals of its “parts”, see, e.g.:

Thanks for quick reply @Wile_E_Coyote .
Once I replace

histname->Fit("myfuncTest1");

To:

TFitResultPtr r1 = histname->Fit(f1, "WEMRS");

I notice a covariance matrix is negative, is that correct?
Here is the output:

****************************************
Minimizer is Minuit / MigradImproved
Chi2                      =  1.58713e+06
NDf                       =          119
Edm                       =   1.0232e-16
NCalls                    =           91
peak                      =      2.02951   +/-   0.184018      -0.0015934   +0.0015934    (Minos) 
background                =      1.64114   +/-   0.0172051     -0.000148979 +0.000148979  (Minos) 
r1) = 4000

2x2 matrix is as follows

     |      0    |      1    |
-------------------------------
   0 |    0.03386  -0.0006664 
   1 | -0.0006664    0.000296 

Wikipedia → Covariance matrix

Thanks @Wile_E_Coyote, Once I calculate the integral of total function f1 as:

Double_t bin_width = histname->GetXaxis()->GetBinWidth(1);
   Double_t neutron_n_total = f1->Integral(0.45, 2.5)/bin_width;
   Double_t neutron_n_totalError = f1->IntegralError(0.45, 2.5)/bin_width;

I got the following error:

Error in <GSLError>: Error 18 in qags.c at 548 : cannot reach tolerance because of roundoff error
Warning in <TF1::IntegralOneDim>: Error found in integrating function myfuncTest10 in [0.450000,2.500000] using AdaptiveSingular. Result = 1082.182463 +/- 0.000129  - status = 18
Info in <TF1::IntegralOneDim>: 		Function Parameters = { peak =  2.029512 , background =  1.641142 }
For Yields = 86574.6 +/- 945.771

Try to play with the “epsrel” parameter (e.g, 1.e-4 ... 1.e-6 ... 1.e-8 ... 1e-10 ... 1.e-12):
Double_t neutron_n_total = f1->Integral(0.45, 2.5, 1.e-6) / bin_width;

Thanks a lot @Wile_E_Coyote for your help, now it works.
The last thing I want to check it with you is calculating integral error for only signal

TMatrixD c1 = r1->GetCovarianceMatrix();
   c1.Print();
   
   TMatrixD c1g = c1.GetSub(0, 0, 0, 1);
   c1g.Print();
   
   TF1 *sig_gaus = new TF1("sig_gaus",peak,0.45, 2.5, 2); 
   sig_gaus->SetParameter(0,par0);
   Double_t neutron_n_signal = sig_gaus->Integral(0.45, 2.5, 1.e-6)/bin_width;
   Double_t neutron_n_signalError = sig_gaus->IntegralError(0.45, 2.5, sig_gaus->GetParameters(), c1g.GetMatrixArray())/bin_width;

The output looks fine for me but I would like to make sure

2x2 matrix is as follows

     |      0    |      1    |
-------------------------------
   0 |    0.03386  -0.0006664 
   1 | -0.0006664    0.000296 


1x2 matrix is as follows

     |      0    |      1    |
-------------------------------
   0 |    0.03386  -0.0006664 

For Yields = 7261.59 +/- 658.415

Something like this will only work if “sig_gaus” is the same function that was used to get “r1” (full “c1” should be used).
Otherwise, you need to create a (square!) “c1g” covariance matrix with elements appropriate for the “sig_gaus” function (i.e., for all its parameters).

See also:

I really appreciate you help @Wile_E_Coyote but I am still confused of how I can do this. So, since I have only 2 parameters
I replaced this part

TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
    c1g.Print();

    if (use_partial_covariances) c1g -=
                                 c1.GetSub(0, 2, 3, 5) *
                                 c1.GetSub(3, 5, 3, 5).InvertFast() *
                                 c1.GetSub(3, 5, 0, 2);

To :

TMatrixD c1g = c1.GetSub(0, 1, 0, 1);
   c1g.Print();
   
   c1g -=c1.GetSub(0, 1, 1, 1) *
         c1.GetSub(1, 1, 1, 1).InvertFast() *
         c1.GetSub(1, 1, 0, 1);

If you want to “extract” the first parameter only (out of the two) …

TMatrixD c1g = c1.GetSub(0, 0, 0, 0); // c1 is the "full" 2x2 matrix
c1g.Print();

if (use_partial_covariances) { c1g -=
                                   c1.GetSub(0, 0, 1, 1) *
                                   c1.GetSub(1, 1, 1, 1).InvertFast() *
                                   c1.GetSub(1, 1, 0, 0);
							   c1g.Print(); }

So then for extract the second parameter only, it should look like this:

TMatrixD c1bc = c1.GetSub(1, 1, 1, 1);
   c1bc.Print();
   if (use_partial_covariances_bg) c1bc -=
                               c1.GetSub(1, 1, 0, 0) *
                               c1.GetSub(1, 1, 1, 1).InvertFast() *
                               c1.GetSub(0, 0, 1, 1);

Correct? Thanks a lot @Wile_E_Coyote for your helping :slight_smile:

The optional “partial covariance” correction is wrong.
Talk to some mathematician. :wink:

OK how about now @Wile_E_Coyote

if (use_partial_covariances_bg) c1bc -=
                               c1.GetSub(1, 1, 0, 1) *
                               c1.GetSub(0, 1, 0, 1).InvertFast() *
                               c1.GetSub(0, 1, 1, 1);

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