Simple Integral

Roofit and RooStats
1

Hi there,

This might result being a very dumb mistake from myself, but I do not understand what is happening with a simple piece of code. I am trying to get the integral of a Gaussian distribution (this is not my final goal, but it is a test to ensure my final answer will be correct) using the following structure:

    double pi = TMath::Pi();
    double Lon = 1.0;
    double sig = 0.1;
    double Lmin = -15.0+Lon;
    double Lmax = +15.0+Lon;
    
    TF1 *Prob = new TF1("Prob",Form("((1.0/(%f*sqrt(2*%f)))*exp(-0.5*pow((x-%f)/(%f),2)))",sig,pi,Lon,sig),Lmin,Lmax);
    
    double average = 0.0;
    average = Prob->Integral(Lmin, Lmax);

Then, when I print the value of average, I get crazy values (I am expecting the result to be one).
Am I making a mistake? Or something I don’t understand…?

And just to remark: when I test the defined function by printing some values (using the Eval function), it actually gives a Gaussian.

Thank you for your help and apologise for this “trivial” issue.

Mario AAO.

2

I guess it’s a deficiency in ROOT.
Your function is always 0 except for a very sharp peak around x=1. This fools the “integrator” (it fails to “detect” this peak).

3

Hi,

I think I disagree. You see, if I print out some values of the function by using

    for (int i = 0 ; i < 10 ; i++) {
        double fac = i/5.0;
        double testVal = Prob->Eval(fac*Lon);
        std::cout << "L = " << fac*Lon << " The test value is = " << testVal << std::endl;
    }

the result is

L = 0 The test value is = 7.6946e-22
L = 0.2 The test value is = 5.05227e-14
L = 0.4 The test value is = 6.07588e-08
L = 0.6 The test value is = 0.0013383
L = 0.8 The test value is = 0.53991
L = 1 The test value is = 3.98942
L = 1.2 The test value is = 0.53991
L = 1.4 The test value is = 0.0013383
L = 1.6 The test value is = 6.07588e-08
L = 1.8 The test value is = 5.05227e-14

I could also plot it, but I think these values make clear that the function is working as expected, i.e., as a Gaussian with a peak in L=1 and a width of 0.1.

What do you think?

Thanks,

Mario AAO.

4

Hello again,

So, I solved the problem, although I must say I don’t know the reason for this.

What happens is that

does not accept Lmin and Lmax the way I passed them to it, even though I have defined them as doubles, as you can see in my initial message.

I finally changed to

and it works (everything else is just the same as in my initial message).

Does anyone know why this happen?

Cheers,

Mario AAO.

5

You “solved” nothing. Try simply:
double Lon = 2; // 2, 3, 4, …

6

[quote=“Pepe Le Pew”]You “solved” nothing. Try simply:
double Lon = 2; // 2, 3, 4, …[/quote]

You are right!

But it only makes things more unclear, don’t you think?
In order to get it working properly, I have to use either

    double mx = +15.0;
    double mn = -15.0;
    double testInt = Prob->Integral(mn, mx);

or

    double testInt = Prob->Integral(-15.0, 15.0);

As long as the value of Lon is in the integration interval, it will give the correct answer. So, I did not solve the problem and my question remains: what is the problem with my initial procedure (and the second one, for that matters)?

Is there any kind of “limitations” on the Integral() function? May be I cannot pass any operator as a parameter to this function (and why does it work for Lon=1.0)?

Thanks,

Mario AAO.