Create a function based on a simulation

[lamya86]

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

---

[Wile_E_Coyote]

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

[lamya86]

@Wile_E_Coyote

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.

[Wile_E_Coyote]

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.:

Compute uncertainty of integral of a component of a TF1

Hi, I am wondering whether there is functionality in standard ROOT (not RooFit) to compute the uncertainty of a component of a TF1. E.g. if fitting a function TF1 f1(“f1”,“gaus(0)+pol2(3)”) to a spectrum, I’d like to get the integral and its uncertainty only for the gaus(0) component. If there is no function available, how can I properly can do this myself? Another related question: If using a signal component, where the “amplitude” parameter already represents the integral (like it should b…

[lamya86]

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 

[Wile_E_Coyote]

Wikipedia → Covariance matrix

[lamya86]

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

[Wile_E_Coyote]

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;

[lamya86]

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

[Wile_E_Coyote]

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:

IntegralError give impossible result

Hello, I want to use IntegralError function for fit integral two TH1. For first TH1 IntegralError is ok. For second IntegralError return 0, but it is impossible. What to do? A cdoe: auto fitResult = MPi0Rho->Fit("tot", "LS", "", inf, sup); auto covMatrix = fitResult->GetCovarianceMatrix(); auto BW2_Phi = [](const double* x, const double* p) -> double { double E = x[0], E2 = E*E; double M = p[5], M2 = M*M; double G = p[6], G2 = G*G; double gamma …

[lamya86]

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);

[Wile_E_Coyote]

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(); }

[lamya86]

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

[Wile_E_Coyote]

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

[lamya86]

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);