Re: Unexpected result from the integration over 2D gauss function

In your code you are creating a 1-D function TF1, but looking at your function integrand_PP, you are using it as if it was a 2-D function. You should create a TF2 instead with the correct range and parameters setting.
Your code produces an error message indicating that your call to TF1 is wrong.

Rene Brun

On
Sun, 22 Jul 2007, Andrey Kormilitzin wrote:

> Dear Rooter,
>
> I need to perform a very simple 2D integration over 2D gaussian function (x,y):
>
> \int_0^inf \int_0^inf dx dy 4*pi^2*x*y*exp(- (x^2+y^2)/S^2 )/pi^2/S^2
>
> The correct answer is 1, but I obtain another result (25.6846). Actually, there
> is no reason for that. I checked the 1D integration over the 1D gaussian
> function and it works well. I want to overcome the change of variables, since
> my real integrand looks more complicated as a function of 'x' and 'y'.
>
> The code is very simple:
>
> //**************************************************************************
> #define S_p 1.0
> #define PI TMath::Pi()
> #define eps 1e-8
>
> double proton_wave_function(double x, double y)
> {// normalized 2D gaussian
>
> return TMath::Exp( - (x*x + y*y)/S_p/S_p )/PI/PI/S_p/S_p;
> }
>
> double integrand_PP(double* var, double* par)
> {
> double x = var[0];
> double y = var[1];
>
> double t1 = ((1-x)/x);
> double J_t1 = 2*PI*(1/x/x)*t1;
>
> double t2 = ((1-y)/y);
> double J_t2 = 2*PI*(1/y/y)*t2;
>
> return J_t1*J_t2*proton_wave_function(t1, t2);
> }
>
> double sigma_tot_PP()
> {
> double limit_d[2] = {eps,eps} ,limit_u[2] = {1-eps,1-eps}, ifail=0.0,
> relerr=0.0, nfnevl=0.0;
> int min_p = 10, as1, as2;
>
> TF1* soft = new TF1("soft_int",integrand_PP,1);
>
> return (soft->IntegralMultiple(2, limit_d, limit_u, min_p, 5e4, 1e-20, relerr,
> as1, as2));
>
> }
>
> //**************************************************************************
>
> ** I made a change of variables (x,y)->(t1,t2) in order to change the limits of
> integration from 0->inf to 0->1. J_t1, J_t2 correspond to jacobian of the
> transformation.
>
>
> Is that problem can be solved without the change of
> variables x^2 + y^2 -> r^2 ?
>
>
> Sincerely,
>
> Andrey K.
>
> *************************************************************
> * Andrey Kormilitzin * *
> * * *
> * Department of * *
> * Particle Physics * *
> * School of Physics * Tel: ++ 972- 3 - 640 7954 (o) *
> * Tel Aviv University * *
> * Ramat Aviv,Tel Aviv * E-mail: andreyk1_at_post.tau.ac.il *
> * 69978, ISRAEL * *
> * * *
> *************************************************************
>
>
> ----------------------------------------------------------------
> This message was sent using IMP, the Internet Messaging Program.
>
>
Received on Sat Jul 28 2007 - 08:15:14 CEST