cbmroot/NumericalIntegration_8h_source.html

#ifndef NUMERICAL_INTEGRATION_H

#define NUMERICAL_INTEGRATION_H

double Integrate2DLaguerre32Legendre32(double (*func)(double, double),

                                       double ay,

                                       double by) {

  // Integrate 2D function from 0 to infinity using Gauss-Laguerre integration

  // with 32 points and from ay to by using Gauss-Legendre integration with 32 points

  //


  Double_t xleg[32];

  Double_t wleg[32];

  Double_t x[32];

  Double_t w[32];


  xleg[0]  = -0.997263861849;

  xleg[1]  = -0.985611511545;

  xleg[2]  = -0.964762255588;

  xleg[3]  = -0.934906075938;

  xleg[4]  = -0.896321155766;

  xleg[5]  = -0.849367613733;

  xleg[6]  = -0.794483795968;

  xleg[7]  = -0.732182118740;

  xleg[8]  = -0.663044266930;

  xleg[9]  = -0.587715757241;

  xleg[10] = -0.506899908932;

  xleg[11] = -0.421351276131;

  xleg[12] = -0.331868602282;

  xleg[13] = -0.239287362252;

  xleg[14] = -0.144471961583;

  xleg[15] = -0.048307665688;

  xleg[16] = 0.048307665688;

  xleg[17] = 0.144471961583;

  xleg[18] = 0.239287362252;

  xleg[19] = 0.331868602282;

  xleg[20] = 0.421351276131;

  xleg[21] = 0.506899908932;

  xleg[22] = 0.587715757241;

  xleg[23] = 0.663044266930;

  xleg[24] = 0.732182118740;

  xleg[25] = 0.794483795968;

  xleg[26] = 0.849367613733;

  xleg[27] = 0.896321155766;

  xleg[28] = 0.934906075938;

  xleg[29] = 0.964762255588;

  xleg[30] = 0.985611511545;

  xleg[31] = 0.997263861849;


  wleg[0]  = 0.007018610009;

  wleg[1]  = 0.016274394716;

  wleg[2]  = 0.025392065309;

  wleg[3]  = 0.034273862913;

  wleg[4]  = 0.042835898022;

  wleg[5]  = 0.050998059262;

  wleg[6]  = 0.058684093479;

  wleg[7]  = 0.065822222776;

  wleg[8]  = 0.072345794109;

  wleg[9]  = 0.078193895787;

  wleg[10] = 0.083311924227;

  wleg[11] = 0.087652093004;

  wleg[12] = 0.091173878696;

  wleg[13] = 0.093844399081;

  wleg[14] = 0.095638720079;

  wleg[15] = 0.096540088515;

  wleg[16] = 0.096540088515;

  wleg[17] = 0.095638720079;

  wleg[18] = 0.093844399081;

  wleg[19] = 0.091173878696;

  wleg[20] = 0.087652093004;

  wleg[21] = 0.083311924227;

  wleg[22] = 0.078193895787;

  wleg[23] = 0.072345794109;

  wleg[24] = 0.065822222776;

  wleg[25] = 0.058684093479;

  wleg[26] = 0.050998059262;

  wleg[27] = 0.042835898022;

  wleg[28] = 0.034273862913;

  wleg[29] = 0.025392065309;

  wleg[30] = 0.016274394716;

  wleg[31] = 0.007018610009;


  Double_t xlag[32];

  Double_t wlag[32];


  xlag[0]  = 0.044489365833;

  xlag[1]  = 0.234526109520;

  xlag[2]  = 0.576884629302;

  xlag[3]  = 1.072448753818;

  xlag[4]  = 1.722408776445;

  xlag[5]  = 2.528336706426;

  xlag[6]  = 3.492213273022;

  xlag[7]  = 4.616456769750;

  xlag[8]  = 5.903958504174;

  xlag[9]  = 7.358126733186;

  xlag[10] = 8.982940924213;

  xlag[11] = 10.783018632540;

  xlag[12] = 12.763697986743;

  xlag[13] = 14.931139755523;

  xlag[14] = 17.292454336715;

  xlag[15] = 19.855860940336;

  xlag[16] = 22.630889013197;

  xlag[17] = 25.628636022459;

  xlag[18] = 28.862101816323;

  xlag[19] = 32.346629153965;

  xlag[20] = 36.100494805752;

  xlag[21] = 40.145719771539;

  xlag[22] = 44.509207995755;

  xlag[23] = 49.224394987309;

  xlag[24] = 54.333721333397;

  xlag[25] = 59.892509162134;

  xlag[26] = 65.975377287935;

  xlag[27] = 72.687628090663;

  xlag[28] = 80.187446977914;

  xlag[29] = 88.735340417892;

  xlag[30] = 98.829542868284;

  xlag[31] = 111.751398097938;


  wlag[0]  = 0.114187105768;

  wlag[1]  = 0.266065216898;

  wlag[2]  = 0.418793137325;

  wlag[3]  = 0.572532846500;

  wlag[4]  = 0.727648788381;

  wlag[5]  = 0.884536719340;

  wlag[6]  = 1.043618875892;

  wlag[7]  = 1.205349274152;

  wlag[8]  = 1.370221338522;

  wlag[9]  = 1.538777256469;

  wlag[10] = 1.711619352686;

  wlag[11] = 1.889424063449;

  wlag[12] = 2.072959340247;

  wlag[13] = 2.263106633997;

  wlag[14] = 2.460889072488;

  wlag[15] = 2.667508126397;

  wlag[16] = 2.884392092922;

  wlag[17] = 3.113261327040;

  wlag[18] = 3.356217692596;

  wlag[19] = 3.615869856484;

  wlag[20] = 3.895513044949;

  wlag[21] = 4.199394104712;

  wlag[22] = 4.533114978534;

  wlag[23] = 4.904270287611;

  wlag[24] = 5.323500972024;

  wlag[25] = 5.806333214234;

  wlag[26] = 6.376614674160;

  wlag[27] = 7.073526580707;

  wlag[28] = 7.967693509296;

  wlag[29] = 9.205040331278;

  wlag[30] = 11.163013090768;

  wlag[31] = 15.390180415261;


  Double_t sum = 0.;


  for (Int_t i = 0; i < 32; i++) {

    for (Int_t j = 0; j < 32; j++) {

      x[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;

      w[j] = (by - ay) / 2. * wleg[j];


      sum += wlag[i] * w[j] * func(xlag[i], x[j]);

    }

  }

  return sum;

}


void GetCoefs2DLaguerre32Legendre32(double ay,

                                    double by,

                                    std::vector<double>& xlag,

                                    std::vector<double>& wlag,

                                    std::vector<double>& xleg,

                                    std::vector<double>& wleg) {

  xlag.resize(32);

  wlag.resize(32);

  xleg.resize(32);

  wleg.resize(32);


  xleg[0]  = -0.997263861849;

  xleg[1]  = -0.985611511545;

  xleg[2]  = -0.964762255588;

  xleg[3]  = -0.934906075938;

  xleg[4]  = -0.896321155766;

  xleg[5]  = -0.849367613733;

  xleg[6]  = -0.794483795968;

  xleg[7]  = -0.732182118740;

  xleg[8]  = -0.663044266930;

  xleg[9]  = -0.587715757241;

  xleg[10] = -0.506899908932;

  xleg[11] = -0.421351276131;

  xleg[12] = -0.331868602282;

  xleg[13] = -0.239287362252;

  xleg[14] = -0.144471961583;

  xleg[15] = -0.048307665688;

  xleg[16] = 0.048307665688;

  xleg[17] = 0.144471961583;

  xleg[18] = 0.239287362252;

  xleg[19] = 0.331868602282;

  xleg[20] = 0.421351276131;

  xleg[21] = 0.506899908932;

  xleg[22] = 0.587715757241;

  xleg[23] = 0.663044266930;

  xleg[24] = 0.732182118740;

  xleg[25] = 0.794483795968;

  xleg[26] = 0.849367613733;

  xleg[27] = 0.896321155766;

  xleg[28] = 0.934906075938;

  xleg[29] = 0.964762255588;

  xleg[30] = 0.985611511545;

  xleg[31] = 0.997263861849;


  wleg[0]  = 0.007018610009;

  wleg[1]  = 0.016274394716;

  wleg[2]  = 0.025392065309;

  wleg[3]  = 0.034273862913;

  wleg[4]  = 0.042835898022;

  wleg[5]  = 0.050998059262;

  wleg[6]  = 0.058684093479;

  wleg[7]  = 0.065822222776;

  wleg[8]  = 0.072345794109;

  wleg[9]  = 0.078193895787;

  wleg[10] = 0.083311924227;

  wleg[11] = 0.087652093004;

  wleg[12] = 0.091173878696;

  wleg[13] = 0.093844399081;

  wleg[14] = 0.095638720079;

  wleg[15] = 0.096540088515;

  wleg[16] = 0.096540088515;

  wleg[17] = 0.095638720079;

  wleg[18] = 0.093844399081;

  wleg[19] = 0.091173878696;

  wleg[20] = 0.087652093004;

  wleg[21] = 0.083311924227;

  wleg[22] = 0.078193895787;

  wleg[23] = 0.072345794109;

  wleg[24] = 0.065822222776;

  wleg[25] = 0.058684093479;

  wleg[26] = 0.050998059262;

  wleg[27] = 0.042835898022;

  wleg[28] = 0.034273862913;

  wleg[29] = 0.025392065309;

  wleg[30] = 0.016274394716;

  wleg[31] = 0.007018610009;


  xlag[0]  = 0.044489365833;

  xlag[1]  = 0.234526109520;

  xlag[2]  = 0.576884629302;

  xlag[3]  = 1.072448753818;

  xlag[4]  = 1.722408776445;

  xlag[5]  = 2.528336706426;

  xlag[6]  = 3.492213273022;

  xlag[7]  = 4.616456769750;

  xlag[8]  = 5.903958504174;

  xlag[9]  = 7.358126733186;

  xlag[10] = 8.982940924213;

  xlag[11] = 10.783018632540;

  xlag[12] = 12.763697986743;

  xlag[13] = 14.931139755523;

  xlag[14] = 17.292454336715;

  xlag[15] = 19.855860940336;

  xlag[16] = 22.630889013197;

  xlag[17] = 25.628636022459;

  xlag[18] = 28.862101816323;

  xlag[19] = 32.346629153965;

  xlag[20] = 36.100494805752;

  xlag[21] = 40.145719771539;

  xlag[22] = 44.509207995755;

  xlag[23] = 49.224394987309;

  xlag[24] = 54.333721333397;

  xlag[25] = 59.892509162134;

  xlag[26] = 65.975377287935;

  xlag[27] = 72.687628090663;

  xlag[28] = 80.187446977914;

  xlag[29] = 88.735340417892;

  xlag[30] = 98.829542868284;

  xlag[31] = 111.751398097938;


  wlag[0]  = 0.114187105768;

  wlag[1]  = 0.266065216898;

  wlag[2]  = 0.418793137325;

  wlag[3]  = 0.572532846500;

  wlag[4]  = 0.727648788381;

  wlag[5]  = 0.884536719340;

  wlag[6]  = 1.043618875892;

  wlag[7]  = 1.205349274152;

  wlag[8]  = 1.370221338522;

  wlag[9]  = 1.538777256469;

  wlag[10] = 1.711619352686;

  wlag[11] = 1.889424063449;

  wlag[12] = 2.072959340247;

  wlag[13] = 2.263106633997;

  wlag[14] = 2.460889072488;

  wlag[15] = 2.667508126397;

  wlag[16] = 2.884392092922;

  wlag[17] = 3.113261327040;

  wlag[18] = 3.356217692596;

  wlag[19] = 3.615869856484;

  wlag[20] = 3.895513044949;

  wlag[21] = 4.199394104712;

  wlag[22] = 4.533114978534;

  wlag[23] = 4.904270287611;

  wlag[24] = 5.323500972024;

  wlag[25] = 5.806333214234;

  wlag[26] = 6.376614674160;

  wlag[27] = 7.073526580707;

  wlag[28] = 7.967693509296;

  wlag[29] = 9.205040331278;

  wlag[30] = 11.163013090768;

  wlag[31] = 15.390180415261;


  //Double_t sum = 0.;


  //for(Int_t i = 0 ; i < 32 ; i++){

  for (Int_t j = 0; j < 32; j++) {

    xleg[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;

    wleg[j] = (by - ay) / 2. * wleg[j];


    //sum += wlag[i]*w[j]*func(xlag[i],x[j]);

  }

  //}

}


void GetCoefs2DLegendre32Legendre32(double ay,

                                    double by,

                                    double a2y,

                                    double b2y,

                                    std::vector<double>& xlag,

                                    std::vector<double>& wlag,

                                    std::vector<double>& xleg,

                                    std::vector<double>& wleg) {

  xlag.resize(32);

  wlag.resize(32);

  xleg.resize(32);

  wleg.resize(32);


  xleg[0]  = -0.997263861849;

  xleg[1]  = -0.985611511545;

  xleg[2]  = -0.964762255588;

  xleg[3]  = -0.934906075938;

  xleg[4]  = -0.896321155766;

  xleg[5]  = -0.849367613733;

  xleg[6]  = -0.794483795968;

  xleg[7]  = -0.732182118740;

  xleg[8]  = -0.663044266930;

  xleg[9]  = -0.587715757241;

  xleg[10] = -0.506899908932;

  xleg[11] = -0.421351276131;

  xleg[12] = -0.331868602282;

  xleg[13] = -0.239287362252;

  xleg[14] = -0.144471961583;

  xleg[15] = -0.048307665688;

  xleg[16] = 0.048307665688;

  xleg[17] = 0.144471961583;

  xleg[18] = 0.239287362252;

  xleg[19] = 0.331868602282;

  xleg[20] = 0.421351276131;

  xleg[21] = 0.506899908932;

  xleg[22] = 0.587715757241;

  xleg[23] = 0.663044266930;

  xleg[24] = 0.732182118740;

  xleg[25] = 0.794483795968;

  xleg[26] = 0.849367613733;

  xleg[27] = 0.896321155766;

  xleg[28] = 0.934906075938;

  xleg[29] = 0.964762255588;

  xleg[30] = 0.985611511545;

  xleg[31] = 0.997263861849;


  wleg[0]  = 0.007018610009;

  wleg[1]  = 0.016274394716;

  wleg[2]  = 0.025392065309;

  wleg[3]  = 0.034273862913;

  wleg[4]  = 0.042835898022;

  wleg[5]  = 0.050998059262;

  wleg[6]  = 0.058684093479;

  wleg[7]  = 0.065822222776;

  wleg[8]  = 0.072345794109;

  wleg[9]  = 0.078193895787;

  wleg[10] = 0.083311924227;

  wleg[11] = 0.087652093004;

  wleg[12] = 0.091173878696;

  wleg[13] = 0.093844399081;

  wleg[14] = 0.095638720079;

  wleg[15] = 0.096540088515;

  wleg[16] = 0.096540088515;

  wleg[17] = 0.095638720079;

  wleg[18] = 0.093844399081;

  wleg[19] = 0.091173878696;

  wleg[20] = 0.087652093004;

  wleg[21] = 0.083311924227;

  wleg[22] = 0.078193895787;

  wleg[23] = 0.072345794109;

  wleg[24] = 0.065822222776;

  wleg[25] = 0.058684093479;

  wleg[26] = 0.050998059262;

  wleg[27] = 0.042835898022;

  wleg[28] = 0.034273862913;

  wleg[29] = 0.025392065309;

  wleg[30] = 0.016274394716;

  wleg[31] = 0.007018610009;


  xlag[0]  = 0.044489365833;

  xlag[1]  = 0.234526109520;

  xlag[2]  = 0.576884629302;

  xlag[3]  = 1.072448753818;

  xlag[4]  = 1.722408776445;

  xlag[5]  = 2.528336706426;

  xlag[6]  = 3.492213273022;

  xlag[7]  = 4.616456769750;

  xlag[8]  = 5.903958504174;

  xlag[9]  = 7.358126733186;

  xlag[10] = 8.982940924213;

  xlag[11] = 10.783018632540;

  xlag[12] = 12.763697986743;

  xlag[13] = 14.931139755523;

  xlag[14] = 17.292454336715;

  xlag[15] = 19.855860940336;

  xlag[16] = 22.630889013197;

  xlag[17] = 25.628636022459;

  xlag[18] = 28.862101816323;

  xlag[19] = 32.346629153965;

  xlag[20] = 36.100494805752;

  xlag[21] = 40.145719771539;

  xlag[22] = 44.509207995755;

  xlag[23] = 49.224394987309;

  xlag[24] = 54.333721333397;

  xlag[25] = 59.892509162134;

  xlag[26] = 65.975377287935;

  xlag[27] = 72.687628090663;

  xlag[28] = 80.187446977914;

  xlag[29] = 88.735340417892;

  xlag[30] = 98.829542868284;

  xlag[31] = 111.751398097938;


  wlag[0]  = 0.114187105768;

  wlag[1]  = 0.266065216898;

  wlag[2]  = 0.418793137325;

  wlag[3]  = 0.572532846500;

  wlag[4]  = 0.727648788381;

  wlag[5]  = 0.884536719340;

  wlag[6]  = 1.043618875892;

  wlag[7]  = 1.205349274152;

  wlag[8]  = 1.370221338522;

  wlag[9]  = 1.538777256469;

  wlag[10] = 1.711619352686;

  wlag[11] = 1.889424063449;

  wlag[12] = 2.072959340247;

  wlag[13] = 2.263106633997;

  wlag[14] = 2.460889072488;

  wlag[15] = 2.667508126397;

  wlag[16] = 2.884392092922;

  wlag[17] = 3.113261327040;

  wlag[18] = 3.356217692596;

  wlag[19] = 3.615869856484;

  wlag[20] = 3.895513044949;

  wlag[21] = 4.199394104712;

  wlag[22] = 4.533114978534;

  wlag[23] = 4.904270287611;

  wlag[24] = 5.323500972024;

  wlag[25] = 5.806333214234;

  wlag[26] = 6.376614674160;

  wlag[27] = 7.073526580707;

  wlag[28] = 7.967693509296;

  wlag[29] = 9.205040331278;

  wlag[30] = 11.163013090768;

  wlag[31] = 15.390180415261;


  //Double_t sum = 0.;


  //for(Int_t i = 0 ; i < 32 ; i++){

  for (Int_t j = 0; j < 32; j++) {

    xlag[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;

    wlag[j] = (by - ay) / 2. * wleg[j];

    xleg[j] = (b2y - a2y) / 2. * xleg[j] + (b2y + a2y) / 2.;

    wleg[j] = (b2y - a2y) / 2. * wleg[j];


    //sum += wlag[i]*w[j]*func(xlag[i],x[j]);

  }

  //}

}


void GetCoefsIntegrateLegendre32(double a,

                                 double b,

                                 std::vector<double>& x,

                                 std::vector<double>& w) {

  // Integrate function from a to b using Legendre-Gaussian integration

  // with 32 points.

  //


  x.resize(32);

  w.resize(32);


  x[0]  = -0.997263861849;

  x[1]  = -0.985611511545;

  x[2]  = -0.964762255588;

  x[3]  = -0.934906075938;

  x[4]  = -0.896321155766;

  x[5]  = -0.849367613733;

  x[6]  = -0.794483795968;

  x[7]  = -0.732182118740;

  x[8]  = -0.663044266930;

  x[9]  = -0.587715757241;

  x[10] = -0.506899908932;

  x[11] = -0.421351276131;

  x[12] = -0.331868602282;

  x[13] = -0.239287362252;

  x[14] = -0.144471961583;

  x[15] = -0.048307665688;

  x[16] = 0.048307665688;

  x[17] = 0.144471961583;

  x[18] = 0.239287362252;

  x[19] = 0.331868602282;

  x[20] = 0.421351276131;

  x[21] = 0.506899908932;

  x[22] = 0.587715757241;

  x[23] = 0.663044266930;

  x[24] = 0.732182118740;

  x[25] = 0.794483795968;

  x[26] = 0.849367613733;

  x[27] = 0.896321155766;

  x[28] = 0.934906075938;

  x[29] = 0.964762255588;

  x[30] = 0.985611511545;

  x[31] = 0.997263861849;


  w[0]  = 0.007018610009;

  w[1]  = 0.016274394716;

  w[2]  = 0.025392065309;

  w[3]  = 0.034273862913;

  w[4]  = 0.042835898022;

  w[5]  = 0.050998059262;

  w[6]  = 0.058684093479;

  w[7]  = 0.065822222776;

  w[8]  = 0.072345794109;

  w[9]  = 0.078193895787;

  w[10] = 0.083311924227;

  w[11] = 0.087652093004;

  w[12] = 0.091173878696;

  w[13] = 0.093844399081;

  w[14] = 0.095638720079;

  w[15] = 0.096540088515;

  w[16] = 0.096540088515;

  w[17] = 0.095638720079;

  w[18] = 0.093844399081;

  w[19] = 0.091173878696;

  w[20] = 0.087652093004;

  w[21] = 0.083311924227;

  w[22] = 0.078193895787;

  w[23] = 0.072345794109;

  w[24] = 0.065822222776;

  w[25] = 0.058684093479;

  w[26] = 0.050998059262;

  w[27] = 0.042835898022;

  w[28] = 0.034273862913;

  w[29] = 0.025392065309;

  w[30] = 0.016274394716;

  w[31] = 0.007018610009;


  for (int i = 0; i < 32; i++) {

    w[i] = (b - a) / 2. * w[i];

    x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;

  }


  /*Double_t sum = 0.;


  for(Int_t i = 0 ; i < 32 ; i++){

    w[i] = (b-a)/2.*w[i];

    x[i] = (b-a)/2.*x[i] + (b+a)/2.;

    sum += w[i]*func->Eval(x[i]);

  }


  return sum;*/

}


void GetCoefsIntegrateLegendre10(double a,

                                 double b,

                                 std::vector<double>& x,

                                 std::vector<double>& w) {

  // Integrate function from a to b using Legendre-Gaussian integration

  // with 32 points.

  //


  x.resize(10);

  w.resize(10);


  x[0] = -0.973906529;

  x[1] = -0.8650633667;

  x[2] = -0.6794095683;

  x[3] = -0.4333953941;

  x[4] = -0.148874339;

  x[5] = 0.148874339;

  x[6] = 0.433395394;

  x[7] = 0.6794095683;

  x[8] = 0.8650633667;

  x[9] = 0.973906529;


  w[0] = 0.0666713443;

  w[1] = 0.1494513491;

  w[2] = 0.219086363;

  w[3] = 0.2692667193;

  w[4] = 0.2955242247;

  w[5] = 0.295524225;

  w[6] = 0.2692667193;

  w[7] = 0.219086363;

  w[8] = 0.149451349;

  w[9] = 0.06667134431;


  for (int i = 0; i < 10; i++) {

    w[i] = (b - a) / 2. * w[i];

    x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;

  }


  /*Double_t sum = 0.;


  for(Int_t i = 0 ; i < 32 ; i++){

    w[i] = (b-a)/2.*w[i];

    x[i] = (b-a)/2.*x[i] + (b+a)/2.;

    sum += w[i]*func->Eval(x[i]);

  }


  return sum;*/

}


void GetCoefsIntegrateLegendre5(double a,

                                double b,

                                std::vector<double>& x,

                                std::vector<double>& w) {

  // Integrate function from a to b using Legendre-Gaussian integration

  // with 32 points.

  //


  x.resize(5);

  w.resize(5);


  x[0] = -0.906179846;

  x[1] = -0.5384693101;

  x[2] = 0;

  x[3] = 0.5384693101;

  x[4] = 0.9061798459;


  w[0] = 0.2369268851;

  w[1] = 0.4786286705;

  w[2] = 0.568888889;

  w[3] = 0.47862867;

  w[4] = 0.2369268851;


  for (int i = 0; i < 5; i++) {

    w[i] = (b - a) / 2. * w[i];

    x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;

  }


  /*Double_t sum = 0.;


  for(Int_t i = 0 ; i < 32 ; i++){

    w[i] = (b-a)/2.*w[i];

    x[i] = (b-a)/2.*x[i] + (b+a)/2.;

    sum += w[i]*func->Eval(x[i]);

  }


  return sum;*/

}


void GetCoefsIntegrateLegendre40(double a,

                                 double b,

                                 std::vector<double>& x,

                                 std::vector<double>& w) {

  // Integrate function from a to b using Legendre-Gaussian integration

  // with 40 points.

  //


  x.resize(40);

  w.resize(40);


  x[0]  = -0.998237709711;

  x[1]  = -0.990726238699;

  x[2]  = -0.977259949984;

  x[3]  = -0.957916819214;

  x[4]  = -0.932812808279;

  x[5]  = -0.902098806969;

  x[6]  = -0.865959503212;

  x[7]  = -0.824612230833;

  x[8]  = -0.778305651427;

  x[9]  = -0.727318255190;

  x[10] = -0.671956684614;

  x[11] = -0.612553889668;

  x[12] = -0.549467125095;

  x[13] = -0.483075801686;

  x[14] = -0.413779204372;

  x[15] = -0.341994090826;

  x[16] = -0.268152185007;

  x[17] = -0.192697580701;

  x[18] = -0.116084070675;

  x[19] = -0.038772417506;

  x[20] = 0.038772417506;

  x[21] = 0.116084070675;

  x[22] = 0.192697580701;

  x[23] = 0.268152185007;

  x[24] = 0.341994090826;

  x[25] = 0.413779204372;

  x[26] = 0.483075801686;

  x[27] = 0.549467125095;

  x[28] = 0.612553889668;

  x[29] = 0.671956684614;

  x[30] = 0.727318255190;

  x[31] = 0.778305651427;

  x[32] = 0.824612230833;

  x[33] = 0.865959503212;

  x[34] = 0.902098806969;

  x[35] = 0.932812808279;

  x[36] = 0.957916819214;

  x[37] = 0.977259949984;

  x[38] = 0.990726238699;

  x[39] = 0.998237709711;


  w[0]  = 0.004521277099;

  w[1]  = 0.010498284521;

  w[2]  = 0.016421058382;

  w[3]  = 0.022245849194;

  w[4]  = 0.027937006980;

  w[5]  = 0.033460195283;

  w[6]  = 0.038782167974;

  w[7]  = 0.043870908186;

  w[8]  = 0.048695807635;

  w[9]  = 0.053227846984;

  w[10] = 0.057439769099;

  w[11] = 0.061306242493;

  w[12] = 0.064804013457;

  w[13] = 0.067912045815;

  w[14] = 0.070611647391;

  w[15] = 0.072886582396;

  w[16] = 0.074723169058;

  w[17] = 0.076110361901;

  w[18] = 0.077039818164;

  w[19] = 0.077505947978;

  w[20] = 0.077505947978;

  w[21] = 0.077039818164;

  w[22] = 0.076110361901;

  w[23] = 0.074723169058;

  w[24] = 0.072886582396;

  w[25] = 0.070611647391;

  w[26] = 0.067912045815;

  w[27] = 0.064804013457;

  w[28] = 0.061306242493;

  w[29] = 0.057439769099;

  w[30] = 0.053227846984;

  w[31] = 0.048695807635;

  w[32] = 0.043870908186;

  w[33] = 0.038782167974;

  w[34] = 0.033460195283;

  w[35] = 0.027937006980;

  w[36] = 0.022245849194;

  w[37] = 0.016421058382;

  w[38] = 0.010498284521;

  w[39] = 0.004521277099;


  //Double_t sum = 0.;


  for (int i = 0; i < 40; i++) {

    w[i] = (b - a) / 2. * w[i];

    x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;

    //sum += w1*func->Eval(x1);

  }


  //return sum;

}


void GetCoefsIntegrateLaguerre32(std::vector<double>& x,

                                 std::vector<double>& w) {

  // Integrate function from 0 to infinity using Gauss-Laguerre integration

  // with 32 points

  //


  x.resize(32);

  w.resize(32);


  x[0]  = 0.044489365833;

  x[1]  = 0.234526109520;

  x[2]  = 0.576884629302;

  x[3]  = 1.072448753818;

  x[4]  = 1.722408776445;

  x[5]  = 2.528336706426;

  x[6]  = 3.492213273022;

  x[7]  = 4.616456769750;

  x[8]  = 5.903958504174;

  x[9]  = 7.358126733186;

  x[10] = 8.982940924213;

  x[11] = 10.783018632540;

  x[12] = 12.763697986743;

  x[13] = 14.931139755523;

  x[14] = 17.292454336715;

  x[15] = 19.855860940336;

  x[16] = 22.630889013197;

  x[17] = 25.628636022459;

  x[18] = 28.862101816323;

  x[19] = 32.346629153965;

  x[20] = 36.100494805752;

  x[21] = 40.145719771539;

  x[22] = 44.509207995755;

  x[23] = 49.224394987309;

  x[24] = 54.333721333397;

  x[25] = 59.892509162134;

  x[26] = 65.975377287935;

  x[27] = 72.687628090663;

  x[28] = 80.187446977914;

  x[29] = 88.735340417892;

  x[30] = 98.829542868284;

  x[31] = 111.751398097938;


  w[0]  = 0.114187105768;

  w[1]  = 0.266065216898;

  w[2]  = 0.418793137325;

  w[3]  = 0.572532846500;

  w[4]  = 0.727648788381;

  w[5]  = 0.884536719340;

  w[6]  = 1.043618875892;

  w[7]  = 1.205349274152;

  w[8]  = 1.370221338522;

  w[9]  = 1.538777256469;

  w[10] = 1.711619352686;

  w[11] = 1.889424063449;

  w[12] = 2.072959340247;

  w[13] = 2.263106633997;

  w[14] = 2.460889072488;

  w[15] = 2.667508126397;

  w[16] = 2.884392092922;

  w[17] = 3.113261327040;

  w[18] = 3.356217692596;

  w[19] = 3.615869856484;

  w[20] = 3.895513044949;

  w[21] = 4.199394104712;

  w[22] = 4.533114978534;

  w[23] = 4.904270287611;

  w[24] = 5.323500972024;

  w[25] = 5.806333214234;

  w[26] = 6.376614674160;

  w[27] = 7.073526580707;

  w[28] = 7.967693509296;

  w[29] = 9.205040331278;

  w[30] = 11.163013090768;

  w[31] = 15.390180415261;


  /*Double_t sum = 0.;


  for(Int_t i = 0 ; i < 32 ; i++){

    sum += w[i]*func->Eval(x[i]);

  }


  return sum;

  */

}


#endif