The deviation for the area under the Gaussian function and the difference between the 2 areas for 2 overlapped Gaussians

Hi everybody,
Could you remind me how we can cover entire %100 of the Gaussian area by using Integral(…,…) for the function we defined earlier? I guess, we can’t write -infinity to + infinity as an option.
For instance:

TF1*f1= new TF1("f1","gaus", -10,10);
f1->SetParameter(0,10);
f1->SetParameter(1,0);
f1->SetParameter(2,2);
f1->Draw();
f1->Integral(-6, 6);  

All we know, for example, we can cover 99.73% if we go range from -3sigma to +3sigma as I chose above from -6 to 6. Thus, I can retrieve the full area by dividing the result of this integral to 0.9973.

Also, the formulation says the area of a Gaussian is height x sigma x sqrt(2*pi) . When I compare the ROOT’s Integral result, they don’t match %100. For instance integral gives 50.128 after percentage correction, but formula gives 50.132.

Another, more important question for me, is there any example to look at to get the blue area I badly marked in the picture I attached. It shows the simple case when 2 Gaussian functions are overlapped.
Mathematically, we can just get the integral of the (f1-f2) in the specified range on X axis. That will give immediately, the blue area. How should I implement that to my code, so ROOT understands? Probably, I will need errors, too. However, let’s get the 1st bit, first.

Appreciate to all in advance.

HI @spyhunter01,
when I compute the Integral in the range (-3sigma,3sigma) with the line below

I obtain 49.997218, that corrected is for the percentage 50.132576.

In any case to integrate gaussian I suggest you to integrate between -5sigma and 5 sigma, the result you obtain is already correct.

Regarding the second point of your question you can define a function as the difference of two gaussians

 TF1*f1= new TF1("f1","gaus(0)-gaus(3)", -15,15);

Then I set the parameters of both the gaussians with SetParameters.

f1->SetParameters(10,0,2,10,2,2)

here I have two gaussian with the same sigma=2 and mean 0 and 2.
The using Integral on such function you obtain the value of the area you are searching for.

Stefano

First of all, thanks Stefano.

Is that correct that your line of " TF1*f1= new TF1(“f1”,“gaus(0)-gaus(3)”, -15,15); " actually represents the red area in my shape in ROOT coding. If so, that’s good news because how simple is that.

Why I got negative value on my screen for the simple example I did? I attached.

Code lines:

TF1 * f1= new TF1(“f1”,“gaus”, 10 , 45) ;
f1->SetParameters(10,20,2) ;
TF1 * f2= new TF1(“f2”,“gaus”, 10 , 45) ;
f2->SetParameters(10,30,3) ;
TF1 * f3= new TF1(“f3”,“gaus(0)-gaus(3)”, 10 , 45) ;
f3->SetParameters(10,20,2,10,30,3)
f1->Draw();
f2->Draw("same);
f3->Draw("same);
f3->Integral(10,45); // the answer is (Double_t)(-2.50662755600367500e+001)
if I say,
f3->Integral(10,24); //(Double_t)4.72812449233806120e+001
I calculated the intersection as 24. The graphs also shows 24. By coincidence, maybe this f3 ends around 24.
Is that because the max point that I gave to integral to calculate?

I think the range from -10 to 12 will be 5 sigma left from the 1st peak and 5 sigma right from the 2nd peak as you said instead of -15 to 15. However, I guess it won’t hurt anything at this stage, right? The main point will be to put the correct range when we write the integral part for this subtraction as f1->Integral(-10,12). Eventually, they go to infinity anyway, and that’s their common limits on both edges.

Draw f3 separately and see what it really looks like.

As dastudillo said the difference of the 2 gaussians is the one showed in the figure below


So a certain point you start to integrate negative value.

If instead you want to integrate just up to the intersection point you should do something like this

TF1*f2= new TF1("f2","gaus(0)>gaus(3)?gaus(0)-gaus(3):0", -10,10);
f2->SetParameters(10,0,2,10,2,2)
f2->Draw()

This is the difference of the two functions if the y of first gaussian is larger of the y of the second one, and is 0 otherwise.

This will give you the right answer.

Stefano

PS

Maybe try to use f2->SetNpx(200) just to have more points in the curve, the curve will look smoother

1 Like

I see. Of course, I haven’t set the Y values to show negative parts. Thus, I missed the point. Also, I haven’t set the integral limits narrower as you said. Therefore, I missed the possibility to have negative value after subtraction in terms of which integral has the most counts underneath it.

I guess I was tired. Thanks, man.

Now, I guess I start to work on the integral errors. Before I start, will that IntegralError(…) method work in this subtraction case as it was working only for 1 Gaussian case?

Cheers.

IntegralError should work fine