CbmRoot
HRGModel/NumericalIntegration.h
Go to the documentation of this file.
1 #ifndef NUMERICAL_INTEGRATION_H
2 #define NUMERICAL_INTEGRATION_H
3 double Integrate2DLaguerre32Legendre32(double (*func)(double, double),
4  double ay,
5  double by) {
6  // Integrate 2D function from 0 to infinity using Gauss-Laguerre integration
7  // with 32 points and from ay to by using Gauss-Legendre integration with 32 points
8  //
9 
10  double xleg[32];
11  double wleg[32];
12  double x[32];
13  double w[32];
14 
15  xleg[0] = -0.997263861849;
16  xleg[1] = -0.985611511545;
17  xleg[2] = -0.964762255588;
18  xleg[3] = -0.934906075938;
19  xleg[4] = -0.896321155766;
20  xleg[5] = -0.849367613733;
21  xleg[6] = -0.794483795968;
22  xleg[7] = -0.732182118740;
23  xleg[8] = -0.663044266930;
24  xleg[9] = -0.587715757241;
25  xleg[10] = -0.506899908932;
26  xleg[11] = -0.421351276131;
27  xleg[12] = -0.331868602282;
28  xleg[13] = -0.239287362252;
29  xleg[14] = -0.144471961583;
30  xleg[15] = -0.048307665688;
31  xleg[16] = 0.048307665688;
32  xleg[17] = 0.144471961583;
33  xleg[18] = 0.239287362252;
34  xleg[19] = 0.331868602282;
35  xleg[20] = 0.421351276131;
36  xleg[21] = 0.506899908932;
37  xleg[22] = 0.587715757241;
38  xleg[23] = 0.663044266930;
39  xleg[24] = 0.732182118740;
40  xleg[25] = 0.794483795968;
41  xleg[26] = 0.849367613733;
42  xleg[27] = 0.896321155766;
43  xleg[28] = 0.934906075938;
44  xleg[29] = 0.964762255588;
45  xleg[30] = 0.985611511545;
46  xleg[31] = 0.997263861849;
47 
48  wleg[0] = 0.007018610009;
49  wleg[1] = 0.016274394716;
50  wleg[2] = 0.025392065309;
51  wleg[3] = 0.034273862913;
52  wleg[4] = 0.042835898022;
53  wleg[5] = 0.050998059262;
54  wleg[6] = 0.058684093479;
55  wleg[7] = 0.065822222776;
56  wleg[8] = 0.072345794109;
57  wleg[9] = 0.078193895787;
58  wleg[10] = 0.083311924227;
59  wleg[11] = 0.087652093004;
60  wleg[12] = 0.091173878696;
61  wleg[13] = 0.093844399081;
62  wleg[14] = 0.095638720079;
63  wleg[15] = 0.096540088515;
64  wleg[16] = 0.096540088515;
65  wleg[17] = 0.095638720079;
66  wleg[18] = 0.093844399081;
67  wleg[19] = 0.091173878696;
68  wleg[20] = 0.087652093004;
69  wleg[21] = 0.083311924227;
70  wleg[22] = 0.078193895787;
71  wleg[23] = 0.072345794109;
72  wleg[24] = 0.065822222776;
73  wleg[25] = 0.058684093479;
74  wleg[26] = 0.050998059262;
75  wleg[27] = 0.042835898022;
76  wleg[28] = 0.034273862913;
77  wleg[29] = 0.025392065309;
78  wleg[30] = 0.016274394716;
79  wleg[31] = 0.007018610009;
80 
81  double xlag[32];
82  double wlag[32];
83 
84  xlag[0] = 0.044489365833;
85  xlag[1] = 0.234526109520;
86  xlag[2] = 0.576884629302;
87  xlag[3] = 1.072448753818;
88  xlag[4] = 1.722408776445;
89  xlag[5] = 2.528336706426;
90  xlag[6] = 3.492213273022;
91  xlag[7] = 4.616456769750;
92  xlag[8] = 5.903958504174;
93  xlag[9] = 7.358126733186;
94  xlag[10] = 8.982940924213;
95  xlag[11] = 10.783018632540;
96  xlag[12] = 12.763697986743;
97  xlag[13] = 14.931139755523;
98  xlag[14] = 17.292454336715;
99  xlag[15] = 19.855860940336;
100  xlag[16] = 22.630889013197;
101  xlag[17] = 25.628636022459;
102  xlag[18] = 28.862101816323;
103  xlag[19] = 32.346629153965;
104  xlag[20] = 36.100494805752;
105  xlag[21] = 40.145719771539;
106  xlag[22] = 44.509207995755;
107  xlag[23] = 49.224394987309;
108  xlag[24] = 54.333721333397;
109  xlag[25] = 59.892509162134;
110  xlag[26] = 65.975377287935;
111  xlag[27] = 72.687628090663;
112  xlag[28] = 80.187446977914;
113  xlag[29] = 88.735340417892;
114  xlag[30] = 98.829542868284;
115  xlag[31] = 111.751398097938;
116 
117  wlag[0] = 0.114187105768;
118  wlag[1] = 0.266065216898;
119  wlag[2] = 0.418793137325;
120  wlag[3] = 0.572532846500;
121  wlag[4] = 0.727648788381;
122  wlag[5] = 0.884536719340;
123  wlag[6] = 1.043618875892;
124  wlag[7] = 1.205349274152;
125  wlag[8] = 1.370221338522;
126  wlag[9] = 1.538777256469;
127  wlag[10] = 1.711619352686;
128  wlag[11] = 1.889424063449;
129  wlag[12] = 2.072959340247;
130  wlag[13] = 2.263106633997;
131  wlag[14] = 2.460889072488;
132  wlag[15] = 2.667508126397;
133  wlag[16] = 2.884392092922;
134  wlag[17] = 3.113261327040;
135  wlag[18] = 3.356217692596;
136  wlag[19] = 3.615869856484;
137  wlag[20] = 3.895513044949;
138  wlag[21] = 4.199394104712;
139  wlag[22] = 4.533114978534;
140  wlag[23] = 4.904270287611;
141  wlag[24] = 5.323500972024;
142  wlag[25] = 5.806333214234;
143  wlag[26] = 6.376614674160;
144  wlag[27] = 7.073526580707;
145  wlag[28] = 7.967693509296;
146  wlag[29] = 9.205040331278;
147  wlag[30] = 11.163013090768;
148  wlag[31] = 15.390180415261;
149 
150  double sum = 0.;
151 
152  for (int i = 0; i < 32; i++) {
153  for (int j = 0; j < 32; j++) {
154  x[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;
155  w[j] = (by - ay) / 2. * wleg[j];
156 
157  sum += wlag[i] * w[j] * func(xlag[i], x[j]);
158  }
159  }
160  return sum;
161 }
162 
163 void GetCoefs2DLaguerre32Legendre32(double ay,
164  double by,
165  std::vector<double>& xlag,
166  std::vector<double>& wlag,
167  std::vector<double>& xleg,
168  std::vector<double>& wleg) {
169  xlag.resize(32);
170  wlag.resize(32);
171  xleg.resize(32);
172  wleg.resize(32);
173 
174  xleg[0] = -0.997263861849;
175  xleg[1] = -0.985611511545;
176  xleg[2] = -0.964762255588;
177  xleg[3] = -0.934906075938;
178  xleg[4] = -0.896321155766;
179  xleg[5] = -0.849367613733;
180  xleg[6] = -0.794483795968;
181  xleg[7] = -0.732182118740;
182  xleg[8] = -0.663044266930;
183  xleg[9] = -0.587715757241;
184  xleg[10] = -0.506899908932;
185  xleg[11] = -0.421351276131;
186  xleg[12] = -0.331868602282;
187  xleg[13] = -0.239287362252;
188  xleg[14] = -0.144471961583;
189  xleg[15] = -0.048307665688;
190  xleg[16] = 0.048307665688;
191  xleg[17] = 0.144471961583;
192  xleg[18] = 0.239287362252;
193  xleg[19] = 0.331868602282;
194  xleg[20] = 0.421351276131;
195  xleg[21] = 0.506899908932;
196  xleg[22] = 0.587715757241;
197  xleg[23] = 0.663044266930;
198  xleg[24] = 0.732182118740;
199  xleg[25] = 0.794483795968;
200  xleg[26] = 0.849367613733;
201  xleg[27] = 0.896321155766;
202  xleg[28] = 0.934906075938;
203  xleg[29] = 0.964762255588;
204  xleg[30] = 0.985611511545;
205  xleg[31] = 0.997263861849;
206 
207  wleg[0] = 0.007018610009;
208  wleg[1] = 0.016274394716;
209  wleg[2] = 0.025392065309;
210  wleg[3] = 0.034273862913;
211  wleg[4] = 0.042835898022;
212  wleg[5] = 0.050998059262;
213  wleg[6] = 0.058684093479;
214  wleg[7] = 0.065822222776;
215  wleg[8] = 0.072345794109;
216  wleg[9] = 0.078193895787;
217  wleg[10] = 0.083311924227;
218  wleg[11] = 0.087652093004;
219  wleg[12] = 0.091173878696;
220  wleg[13] = 0.093844399081;
221  wleg[14] = 0.095638720079;
222  wleg[15] = 0.096540088515;
223  wleg[16] = 0.096540088515;
224  wleg[17] = 0.095638720079;
225  wleg[18] = 0.093844399081;
226  wleg[19] = 0.091173878696;
227  wleg[20] = 0.087652093004;
228  wleg[21] = 0.083311924227;
229  wleg[22] = 0.078193895787;
230  wleg[23] = 0.072345794109;
231  wleg[24] = 0.065822222776;
232  wleg[25] = 0.058684093479;
233  wleg[26] = 0.050998059262;
234  wleg[27] = 0.042835898022;
235  wleg[28] = 0.034273862913;
236  wleg[29] = 0.025392065309;
237  wleg[30] = 0.016274394716;
238  wleg[31] = 0.007018610009;
239 
240  xlag[0] = 0.044489365833;
241  xlag[1] = 0.234526109520;
242  xlag[2] = 0.576884629302;
243  xlag[3] = 1.072448753818;
244  xlag[4] = 1.722408776445;
245  xlag[5] = 2.528336706426;
246  xlag[6] = 3.492213273022;
247  xlag[7] = 4.616456769750;
248  xlag[8] = 5.903958504174;
249  xlag[9] = 7.358126733186;
250  xlag[10] = 8.982940924213;
251  xlag[11] = 10.783018632540;
252  xlag[12] = 12.763697986743;
253  xlag[13] = 14.931139755523;
254  xlag[14] = 17.292454336715;
255  xlag[15] = 19.855860940336;
256  xlag[16] = 22.630889013197;
257  xlag[17] = 25.628636022459;
258  xlag[18] = 28.862101816323;
259  xlag[19] = 32.346629153965;
260  xlag[20] = 36.100494805752;
261  xlag[21] = 40.145719771539;
262  xlag[22] = 44.509207995755;
263  xlag[23] = 49.224394987309;
264  xlag[24] = 54.333721333397;
265  xlag[25] = 59.892509162134;
266  xlag[26] = 65.975377287935;
267  xlag[27] = 72.687628090663;
268  xlag[28] = 80.187446977914;
269  xlag[29] = 88.735340417892;
270  xlag[30] = 98.829542868284;
271  xlag[31] = 111.751398097938;
272 
273  wlag[0] = 0.114187105768;
274  wlag[1] = 0.266065216898;
275  wlag[2] = 0.418793137325;
276  wlag[3] = 0.572532846500;
277  wlag[4] = 0.727648788381;
278  wlag[5] = 0.884536719340;
279  wlag[6] = 1.043618875892;
280  wlag[7] = 1.205349274152;
281  wlag[8] = 1.370221338522;
282  wlag[9] = 1.538777256469;
283  wlag[10] = 1.711619352686;
284  wlag[11] = 1.889424063449;
285  wlag[12] = 2.072959340247;
286  wlag[13] = 2.263106633997;
287  wlag[14] = 2.460889072488;
288  wlag[15] = 2.667508126397;
289  wlag[16] = 2.884392092922;
290  wlag[17] = 3.113261327040;
291  wlag[18] = 3.356217692596;
292  wlag[19] = 3.615869856484;
293  wlag[20] = 3.895513044949;
294  wlag[21] = 4.199394104712;
295  wlag[22] = 4.533114978534;
296  wlag[23] = 4.904270287611;
297  wlag[24] = 5.323500972024;
298  wlag[25] = 5.806333214234;
299  wlag[26] = 6.376614674160;
300  wlag[27] = 7.073526580707;
301  wlag[28] = 7.967693509296;
302  wlag[29] = 9.205040331278;
303  wlag[30] = 11.163013090768;
304  wlag[31] = 15.390180415261;
305 
306  //double sum = 0.;
307 
308  //for(int i = 0 ; i < 32 ; i++){
309  for (int j = 0; j < 32; j++) {
310  xleg[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;
311  wleg[j] = (by - ay) / 2. * wleg[j];
312 
313  //sum += wlag[i]*w[j]*func(xlag[i],x[j]);
314  }
315  //}
316 }
317 
318 void GetCoefs2DLegendre32Legendre32(double ay,
319  double by,
320  double a2y,
321  double b2y,
322  std::vector<double>& xlag,
323  std::vector<double>& wlag,
324  std::vector<double>& xleg,
325  std::vector<double>& wleg) {
326  xlag.resize(32);
327  wlag.resize(32);
328  xleg.resize(32);
329  wleg.resize(32);
330 
331  xleg[0] = -0.997263861849;
332  xleg[1] = -0.985611511545;
333  xleg[2] = -0.964762255588;
334  xleg[3] = -0.934906075938;
335  xleg[4] = -0.896321155766;
336  xleg[5] = -0.849367613733;
337  xleg[6] = -0.794483795968;
338  xleg[7] = -0.732182118740;
339  xleg[8] = -0.663044266930;
340  xleg[9] = -0.587715757241;
341  xleg[10] = -0.506899908932;
342  xleg[11] = -0.421351276131;
343  xleg[12] = -0.331868602282;
344  xleg[13] = -0.239287362252;
345  xleg[14] = -0.144471961583;
346  xleg[15] = -0.048307665688;
347  xleg[16] = 0.048307665688;
348  xleg[17] = 0.144471961583;
349  xleg[18] = 0.239287362252;
350  xleg[19] = 0.331868602282;
351  xleg[20] = 0.421351276131;
352  xleg[21] = 0.506899908932;
353  xleg[22] = 0.587715757241;
354  xleg[23] = 0.663044266930;
355  xleg[24] = 0.732182118740;
356  xleg[25] = 0.794483795968;
357  xleg[26] = 0.849367613733;
358  xleg[27] = 0.896321155766;
359  xleg[28] = 0.934906075938;
360  xleg[29] = 0.964762255588;
361  xleg[30] = 0.985611511545;
362  xleg[31] = 0.997263861849;
363 
364  wleg[0] = 0.007018610009;
365  wleg[1] = 0.016274394716;
366  wleg[2] = 0.025392065309;
367  wleg[3] = 0.034273862913;
368  wleg[4] = 0.042835898022;
369  wleg[5] = 0.050998059262;
370  wleg[6] = 0.058684093479;
371  wleg[7] = 0.065822222776;
372  wleg[8] = 0.072345794109;
373  wleg[9] = 0.078193895787;
374  wleg[10] = 0.083311924227;
375  wleg[11] = 0.087652093004;
376  wleg[12] = 0.091173878696;
377  wleg[13] = 0.093844399081;
378  wleg[14] = 0.095638720079;
379  wleg[15] = 0.096540088515;
380  wleg[16] = 0.096540088515;
381  wleg[17] = 0.095638720079;
382  wleg[18] = 0.093844399081;
383  wleg[19] = 0.091173878696;
384  wleg[20] = 0.087652093004;
385  wleg[21] = 0.083311924227;
386  wleg[22] = 0.078193895787;
387  wleg[23] = 0.072345794109;
388  wleg[24] = 0.065822222776;
389  wleg[25] = 0.058684093479;
390  wleg[26] = 0.050998059262;
391  wleg[27] = 0.042835898022;
392  wleg[28] = 0.034273862913;
393  wleg[29] = 0.025392065309;
394  wleg[30] = 0.016274394716;
395  wleg[31] = 0.007018610009;
396 
397  xlag[0] = 0.044489365833;
398  xlag[1] = 0.234526109520;
399  xlag[2] = 0.576884629302;
400  xlag[3] = 1.072448753818;
401  xlag[4] = 1.722408776445;
402  xlag[5] = 2.528336706426;
403  xlag[6] = 3.492213273022;
404  xlag[7] = 4.616456769750;
405  xlag[8] = 5.903958504174;
406  xlag[9] = 7.358126733186;
407  xlag[10] = 8.982940924213;
408  xlag[11] = 10.783018632540;
409  xlag[12] = 12.763697986743;
410  xlag[13] = 14.931139755523;
411  xlag[14] = 17.292454336715;
412  xlag[15] = 19.855860940336;
413  xlag[16] = 22.630889013197;
414  xlag[17] = 25.628636022459;
415  xlag[18] = 28.862101816323;
416  xlag[19] = 32.346629153965;
417  xlag[20] = 36.100494805752;
418  xlag[21] = 40.145719771539;
419  xlag[22] = 44.509207995755;
420  xlag[23] = 49.224394987309;
421  xlag[24] = 54.333721333397;
422  xlag[25] = 59.892509162134;
423  xlag[26] = 65.975377287935;
424  xlag[27] = 72.687628090663;
425  xlag[28] = 80.187446977914;
426  xlag[29] = 88.735340417892;
427  xlag[30] = 98.829542868284;
428  xlag[31] = 111.751398097938;
429 
430  wlag[0] = 0.114187105768;
431  wlag[1] = 0.266065216898;
432  wlag[2] = 0.418793137325;
433  wlag[3] = 0.572532846500;
434  wlag[4] = 0.727648788381;
435  wlag[5] = 0.884536719340;
436  wlag[6] = 1.043618875892;
437  wlag[7] = 1.205349274152;
438  wlag[8] = 1.370221338522;
439  wlag[9] = 1.538777256469;
440  wlag[10] = 1.711619352686;
441  wlag[11] = 1.889424063449;
442  wlag[12] = 2.072959340247;
443  wlag[13] = 2.263106633997;
444  wlag[14] = 2.460889072488;
445  wlag[15] = 2.667508126397;
446  wlag[16] = 2.884392092922;
447  wlag[17] = 3.113261327040;
448  wlag[18] = 3.356217692596;
449  wlag[19] = 3.615869856484;
450  wlag[20] = 3.895513044949;
451  wlag[21] = 4.199394104712;
452  wlag[22] = 4.533114978534;
453  wlag[23] = 4.904270287611;
454  wlag[24] = 5.323500972024;
455  wlag[25] = 5.806333214234;
456  wlag[26] = 6.376614674160;
457  wlag[27] = 7.073526580707;
458  wlag[28] = 7.967693509296;
459  wlag[29] = 9.205040331278;
460  wlag[30] = 11.163013090768;
461  wlag[31] = 15.390180415261;
462 
463  //double sum = 0.;
464 
465  //for(int i = 0 ; i < 32 ; i++){
466  for (int j = 0; j < 32; j++) {
467  xlag[j] = (by - ay) / 2. * xleg[j] + (by + ay) / 2.;
468  wlag[j] = (by - ay) / 2. * wleg[j];
469  xleg[j] = (b2y - a2y) / 2. * xleg[j] + (b2y + a2y) / 2.;
470  wleg[j] = (b2y - a2y) / 2. * wleg[j];
471 
472  //sum += wlag[i]*w[j]*func(xlag[i],x[j]);
473  }
474  //}
475 }
476 
477 void GetCoefsIntegrateLegendre32(double a,
478  double b,
479  std::vector<double>& x,
480  std::vector<double>& w) {
481  // Integrate function from a to b using Legendre-Gaussian integration
482  // with 32 points.
483  //
484 
485  x.resize(32);
486  w.resize(32);
487 
488  x[0] = -0.997263861849;
489  x[1] = -0.985611511545;
490  x[2] = -0.964762255588;
491  x[3] = -0.934906075938;
492  x[4] = -0.896321155766;
493  x[5] = -0.849367613733;
494  x[6] = -0.794483795968;
495  x[7] = -0.732182118740;
496  x[8] = -0.663044266930;
497  x[9] = -0.587715757241;
498  x[10] = -0.506899908932;
499  x[11] = -0.421351276131;
500  x[12] = -0.331868602282;
501  x[13] = -0.239287362252;
502  x[14] = -0.144471961583;
503  x[15] = -0.048307665688;
504  x[16] = 0.048307665688;
505  x[17] = 0.144471961583;
506  x[18] = 0.239287362252;
507  x[19] = 0.331868602282;
508  x[20] = 0.421351276131;
509  x[21] = 0.506899908932;
510  x[22] = 0.587715757241;
511  x[23] = 0.663044266930;
512  x[24] = 0.732182118740;
513  x[25] = 0.794483795968;
514  x[26] = 0.849367613733;
515  x[27] = 0.896321155766;
516  x[28] = 0.934906075938;
517  x[29] = 0.964762255588;
518  x[30] = 0.985611511545;
519  x[31] = 0.997263861849;
520 
521  w[0] = 0.007018610009;
522  w[1] = 0.016274394716;
523  w[2] = 0.025392065309;
524  w[3] = 0.034273862913;
525  w[4] = 0.042835898022;
526  w[5] = 0.050998059262;
527  w[6] = 0.058684093479;
528  w[7] = 0.065822222776;
529  w[8] = 0.072345794109;
530  w[9] = 0.078193895787;
531  w[10] = 0.083311924227;
532  w[11] = 0.087652093004;
533  w[12] = 0.091173878696;
534  w[13] = 0.093844399081;
535  w[14] = 0.095638720079;
536  w[15] = 0.096540088515;
537  w[16] = 0.096540088515;
538  w[17] = 0.095638720079;
539  w[18] = 0.093844399081;
540  w[19] = 0.091173878696;
541  w[20] = 0.087652093004;
542  w[21] = 0.083311924227;
543  w[22] = 0.078193895787;
544  w[23] = 0.072345794109;
545  w[24] = 0.065822222776;
546  w[25] = 0.058684093479;
547  w[26] = 0.050998059262;
548  w[27] = 0.042835898022;
549  w[28] = 0.034273862913;
550  w[29] = 0.025392065309;
551  w[30] = 0.016274394716;
552  w[31] = 0.007018610009;
553 
554  for (int i = 0; i < 32; i++) {
555  w[i] = (b - a) / 2. * w[i];
556  x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;
557  }
558 
559  /*double sum = 0.;
560 
561  for(int i = 0 ; i < 32 ; i++){
562  w[i] = (b-a)/2.*w[i];
563  x[i] = (b-a)/2.*x[i] + (b+a)/2.;
564  sum += w[i]*func->Eval(x[i]);
565  }
566 
567  return sum;*/
568 }
569 
570 void GetCoefsIntegrateLegendre10(double a,
571  double b,
572  std::vector<double>& x,
573  std::vector<double>& w) {
574  // Integrate function from a to b using Legendre-Gaussian integration
575  // with 32 points.
576  //
577 
578  x.resize(10);
579  w.resize(10);
580 
581  x[0] = -0.973906529;
582  x[1] = -0.8650633667;
583  x[2] = -0.6794095683;
584  x[3] = -0.4333953941;
585  x[4] = -0.148874339;
586  x[5] = 0.148874339;
587  x[6] = 0.433395394;
588  x[7] = 0.6794095683;
589  x[8] = 0.8650633667;
590  x[9] = 0.973906529;
591 
592  w[0] = 0.0666713443;
593  w[1] = 0.1494513491;
594  w[2] = 0.219086363;
595  w[3] = 0.2692667193;
596  w[4] = 0.2955242247;
597  w[5] = 0.295524225;
598  w[6] = 0.2692667193;
599  w[7] = 0.219086363;
600  w[8] = 0.149451349;
601  w[9] = 0.06667134431;
602 
603  for (int i = 0; i < 10; i++) {
604  w[i] = (b - a) / 2. * w[i];
605  x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;
606  }
607 
608  /*double sum = 0.;
609 
610  for(int i = 0 ; i < 32 ; i++){
611  w[i] = (b-a)/2.*w[i];
612  x[i] = (b-a)/2.*x[i] + (b+a)/2.;
613  sum += w[i]*func->Eval(x[i]);
614  }
615 
616  return sum;*/
617 }
618 
619 void GetCoefsIntegrateLegendre5(double a,
620  double b,
621  std::vector<double>& x,
622  std::vector<double>& w) {
623  // Integrate function from a to b using Legendre-Gaussian integration
624  // with 32 points.
625  //
626 
627  x.resize(5);
628  w.resize(5);
629 
630  x[0] = -0.906179846;
631  x[1] = -0.5384693101;
632  x[2] = 0;
633  x[3] = 0.5384693101;
634  x[4] = 0.9061798459;
635 
636  w[0] = 0.2369268851;
637  w[1] = 0.4786286705;
638  w[2] = 0.568888889;
639  w[3] = 0.47862867;
640  w[4] = 0.2369268851;
641 
642  for (int i = 0; i < 5; i++) {
643  w[i] = (b - a) / 2. * w[i];
644  x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;
645  }
646 
647  /*double sum = 0.;
648 
649  for(int i = 0 ; i < 32 ; i++){
650  w[i] = (b-a)/2.*w[i];
651  x[i] = (b-a)/2.*x[i] + (b+a)/2.;
652  sum += w[i]*func->Eval(x[i]);
653  }
654 
655  return sum;*/
656 }
657 
658 void GetCoefsIntegrateLegendre40(double a,
659  double b,
660  std::vector<double>& x,
661  std::vector<double>& w) {
662  // Integrate function from a to b using Legendre-Gaussian integration
663  // with 40 points.
664  //
665 
666  x.resize(40);
667  w.resize(40);
668 
669  x[0] = -0.998237709711;
670  x[1] = -0.990726238699;
671  x[2] = -0.977259949984;
672  x[3] = -0.957916819214;
673  x[4] = -0.932812808279;
674  x[5] = -0.902098806969;
675  x[6] = -0.865959503212;
676  x[7] = -0.824612230833;
677  x[8] = -0.778305651427;
678  x[9] = -0.727318255190;
679  x[10] = -0.671956684614;
680  x[11] = -0.612553889668;
681  x[12] = -0.549467125095;
682  x[13] = -0.483075801686;
683  x[14] = -0.413779204372;
684  x[15] = -0.341994090826;
685  x[16] = -0.268152185007;
686  x[17] = -0.192697580701;
687  x[18] = -0.116084070675;
688  x[19] = -0.038772417506;
689  x[20] = 0.038772417506;
690  x[21] = 0.116084070675;
691  x[22] = 0.192697580701;
692  x[23] = 0.268152185007;
693  x[24] = 0.341994090826;
694  x[25] = 0.413779204372;
695  x[26] = 0.483075801686;
696  x[27] = 0.549467125095;
697  x[28] = 0.612553889668;
698  x[29] = 0.671956684614;
699  x[30] = 0.727318255190;
700  x[31] = 0.778305651427;
701  x[32] = 0.824612230833;
702  x[33] = 0.865959503212;
703  x[34] = 0.902098806969;
704  x[35] = 0.932812808279;
705  x[36] = 0.957916819214;
706  x[37] = 0.977259949984;
707  x[38] = 0.990726238699;
708  x[39] = 0.998237709711;
709 
710  w[0] = 0.004521277099;
711  w[1] = 0.010498284521;
712  w[2] = 0.016421058382;
713  w[3] = 0.022245849194;
714  w[4] = 0.027937006980;
715  w[5] = 0.033460195283;
716  w[6] = 0.038782167974;
717  w[7] = 0.043870908186;
718  w[8] = 0.048695807635;
719  w[9] = 0.053227846984;
720  w[10] = 0.057439769099;
721  w[11] = 0.061306242493;
722  w[12] = 0.064804013457;
723  w[13] = 0.067912045815;
724  w[14] = 0.070611647391;
725  w[15] = 0.072886582396;
726  w[16] = 0.074723169058;
727  w[17] = 0.076110361901;
728  w[18] = 0.077039818164;
729  w[19] = 0.077505947978;
730  w[20] = 0.077505947978;
731  w[21] = 0.077039818164;
732  w[22] = 0.076110361901;
733  w[23] = 0.074723169058;
734  w[24] = 0.072886582396;
735  w[25] = 0.070611647391;
736  w[26] = 0.067912045815;
737  w[27] = 0.064804013457;
738  w[28] = 0.061306242493;
739  w[29] = 0.057439769099;
740  w[30] = 0.053227846984;
741  w[31] = 0.048695807635;
742  w[32] = 0.043870908186;
743  w[33] = 0.038782167974;
744  w[34] = 0.033460195283;
745  w[35] = 0.027937006980;
746  w[36] = 0.022245849194;
747  w[37] = 0.016421058382;
748  w[38] = 0.010498284521;
749  w[39] = 0.004521277099;
750 
751  //double sum = 0.;
752 
753  for (int i = 0; i < 40; i++) {
754  w[i] = (b - a) / 2. * w[i];
755  x[i] = (b - a) / 2. * x[i] + (b + a) / 2.;
756  //sum += w1*func->Eval(x1);
757  }
758 
759  //return sum;
760 }
761 
762 void GetCoefsIntegrateLaguerre32(std::vector<double>& x,
763  std::vector<double>& w) {
764  // Integrate function from 0 to infinity using Gauss-Laguerre integration
765  // with 32 points
766  //
767 
768  x.resize(32);
769  w.resize(32);
770 
771  x[0] = 0.044489365833;
772  x[1] = 0.234526109520;
773  x[2] = 0.576884629302;
774  x[3] = 1.072448753818;
775  x[4] = 1.722408776445;
776  x[5] = 2.528336706426;
777  x[6] = 3.492213273022;
778  x[7] = 4.616456769750;
779  x[8] = 5.903958504174;
780  x[9] = 7.358126733186;
781  x[10] = 8.982940924213;
782  x[11] = 10.783018632540;
783  x[12] = 12.763697986743;
784  x[13] = 14.931139755523;
785  x[14] = 17.292454336715;
786  x[15] = 19.855860940336;
787  x[16] = 22.630889013197;
788  x[17] = 25.628636022459;
789  x[18] = 28.862101816323;
790  x[19] = 32.346629153965;
791  x[20] = 36.100494805752;
792  x[21] = 40.145719771539;
793  x[22] = 44.509207995755;
794  x[23] = 49.224394987309;
795  x[24] = 54.333721333397;
796  x[25] = 59.892509162134;
797  x[26] = 65.975377287935;
798  x[27] = 72.687628090663;
799  x[28] = 80.187446977914;
800  x[29] = 88.735340417892;
801  x[30] = 98.829542868284;
802  x[31] = 111.751398097938;
803 
804  w[0] = 0.114187105768;
805  w[1] = 0.266065216898;
806  w[2] = 0.418793137325;
807  w[3] = 0.572532846500;
808  w[4] = 0.727648788381;
809  w[5] = 0.884536719340;
810  w[6] = 1.043618875892;
811  w[7] = 1.205349274152;
812  w[8] = 1.370221338522;
813  w[9] = 1.538777256469;
814  w[10] = 1.711619352686;
815  w[11] = 1.889424063449;
816  w[12] = 2.072959340247;
817  w[13] = 2.263106633997;
818  w[14] = 2.460889072488;
819  w[15] = 2.667508126397;
820  w[16] = 2.884392092922;
821  w[17] = 3.113261327040;
822  w[18] = 3.356217692596;
823  w[19] = 3.615869856484;
824  w[20] = 3.895513044949;
825  w[21] = 4.199394104712;
826  w[22] = 4.533114978534;
827  w[23] = 4.904270287611;
828  w[24] = 5.323500972024;
829  w[25] = 5.806333214234;
830  w[26] = 6.376614674160;
831  w[27] = 7.073526580707;
832  w[28] = 7.967693509296;
833  w[29] = 9.205040331278;
834  w[30] = 11.163013090768;
835  w[31] = 15.390180415261;
836 
837  /*double sum = 0.;
838 
839  for(int i = 0 ; i < 32 ; i++){
840  sum += w[i]*func->Eval(x[i]);
841  }
842 
843  return sum;
844  */
845 }
846 
847 #endif
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
Integrate2DLaguerre32Legendre32
double Integrate2DLaguerre32Legendre32(double(*func)(double, double), double ay, double by)
Definition: HRGModel/NumericalIntegration.h:3
GetCoefsIntegrateLaguerre32
void GetCoefsIntegrateLaguerre32(std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:762
GetCoefsIntegrateLegendre5
void GetCoefsIntegrateLegendre5(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:619
GetCoefsIntegrateLegendre32
void GetCoefsIntegrateLegendre32(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:477
GetCoefsIntegrateLegendre10
void GetCoefsIntegrateLegendre10(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:570
GetCoefsIntegrateLegendre40
void GetCoefsIntegrateLegendre40(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:658
GetCoefs2DLegendre32Legendre32
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)
Definition: HRGModel/NumericalIntegration.h:318
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
GetCoefs2DLaguerre32Legendre32
void GetCoefs2DLaguerre32Legendre32(double ay, double by, std::vector< double > &xlag, std::vector< double > &wlag, std::vector< double > &xleg, std::vector< double > &wleg)
Definition: HRGModel/NumericalIntegration.h:163