RdGeoCeLDFFitter/LikelihoodFunction.cc
Go to the documentation of this file.
1 #include "LikelihoodFunction.h"
2 
3 #include <iostream>
4 #include <cmath>
5 #include <vector>
6 #include <limits>
7 
8 #include <Minuit2/FCNBase.h>
9 
10 #include <utl/CoordinateSystemPtr.h>
11 #include <utl/Vector.h>
12 #include <utl/BasicVector.h>
13 #include <utl/Point.h>
14 #include <utl/Math.h>
15 #include <utl/PhysicalFunctions.h>
16 #include <utl/RadioGeometryUtilities.h>
17 #include <utl/Probability.h>
18 #include <utl/AugerUnits.h>
19 
20 #include "TF1.h"
21 #include "TF2.h"
22 #include "TMath.h"
23 #include <math.h>
24 
25 
26 namespace RdGeoCeLDFFitter {
27 
28  // returns the value of a BSpline function at position x for internal knots t, coefficients c and order k
29  double BSpline(const double x, std::vector<double> t, const std::vector<double> c, const int k)
30  {
31  for (int i = 0; i < k; ++i) { // add outer knots
32  t.insert(t.begin(), t[0]);
33  t.push_back(t[t.size() - 1]);
34  }
35  const int m = t.size();
36  double B[m - 1][k + 1];
37 
38 
39  for (int i = 0; i < (m - 1); ++i) {
40  if ((x >= t[i]) && (x < t[i + 1])) {
41  B[i][0] = 1.;
42  } else {
43  B[i][0] = 0.;
44  }
45  }
46 
47  for (int j = 1; j < (k + 1); ++j) {
48  for (int i = 0; i < (m - j - 1); ++i) {
49  double first_term = 0;
50  double second_term = 0;
51  if (t[i + j] - t[i] == 0.0) {
52  first_term = 0.0;
53  } else {
54  first_term = ((x - t[i]) / (t[i + j] - t[i])) * B[i][j - 1];
55  }
56  if (t[i + j + 1] - t[i + 1] == 0.0) {
57  second_term = 0.0;
58  } else {
59  second_term = ((t[i + j + 1] - x) / (t[i + j + 1] - t[i + 1])) * B[i + 1][j - 1];
60  }
61  B[i][j] = first_term + second_term;
62  }
63  }
64  double y = 0;
65  for (int i = 0; i < (m - k - 1); ++i) {
66  y += c[i] * B[i][k];
67  }
68 
69  return y;
70  }
71 
72  double getRcutGeo(double x)
73  {
74  static const double t[] = { 466.81311596, 500. , 750. , 1000. ,
75  1250. , 1500. , 1750. , 2000. ,
76  2250. , 2500. , 2750. , 3000. ,
77  3250. , 3500. , 3750. , 3885.68693916 };
78  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
79  static const double c[] = { 230.76422954, 292.23358156, 153.16581136, 195.99815817,
80  289.37971316, 499.73283119, 788.21596404, 1088.74340871,
81  1288.02204404, 1398.97369355, 1498.09425301, 1621.72311365,
82  1749.19342988, 1849.69730923, 2005.52546444, 2110.19363873,
83  2208.87505369, 2160.6098839};
84  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
85  return BSpline(x, knots, coeffs, 3);
86  }
87 
88  double getBGeo(double x)
89  {
90  static const double t[] = {466.81311596, 500. , 750. , 1000. ,
91  1250. , 1500. , 1750. , 2000. ,
92  2250. , 2500. , 2750. , 3000. ,
93  3250. , 3500. , 3750. , 3986.00008755 };
94  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
95  static const double c[] = {155.97606476, 182.2254283 , 182.86021036, 178.86262811,
96  207.31477709, 146.36606353, 147.3021935 , 119.22748464,
97  107.44908171, 149.15858921, 148.00772949, 125.11434689,
98  133.93007258, 125.64835967, 118.2237385 , 140.58832422,
99  83.65478038, 170.77145449};
100  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
101  return BSpline(x, knots, coeffs, 3);
102  }
103 
104  double getGeoR(double x)
105  {
106  static const double t[] = {-200. , 130.11929063, 164.97522784, 200.62942539,
107  224.205947 , 305.22361845, 416.77956938, 795.48741779,
108  1530.28993967, 1866.76266188, 1965.43098131, 2109.17204951,
109  2279.55497793, 2419.10847863, 2563.25426508, 2654.53207932,
110  2751.82229916, 2874.81031606, 2988.51488067, 3055.12705453,
111  3085.36752724, 3152.97188299, 3289.48859202, 3396.7478351 ,
112  3500.33523494, 3537.07041433, 3602.07484135, 3673.28026584,
113  3769.15406 , 3773.79577228, 3785.86008755, 3820.94789224,
114  3946.57706295, 4027.80175655, 4101.69012082, 4485.40709032,
115  4711.07709032, 6755.69846645 };
116  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
117  static const double c[] = { -12.52965737, -31.59360635, -25.24425533, -67.83585992,
118  -79.71324199, -93.60518414, -94.27396049, 70.49357834,
119  147.75729402, 234.2965183 , 295.39764637, 326.80549038,
120  347.32762753, 369.38657812, 394.63663292, 410.28574442,
121  423.9399281 , 443.66582118, 459.37314999, 473.35999024,
122  483.05109551, 490.09444485, 500.16219506, 520.04475797,
123  528.52336809, 546.69229895, 554.58563345, 559.9322956 ,
124  577.37871356, 589.77842791, 594.22847378, 596.00265837,
125  599.98206278, 607.75915513, 619.6370265 , 634.92896525,
126  713.91896463, 762.78472157, 995.47281507, 1051.81132903 };
127  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
128  return BSpline(x, knots, coeffs, 3);
129  }
130 
131  double getGeoSigma(double x)
132  {
133  static const double t[] = { -101.45678401, 0. , 250. , 500. ,
134  750. , 1000. , 1250. , 1500. ,
135  1750. , 2000. , 2250. , 2500. ,
136  2750. , 3000. , 3250. , 3500. ,
137  4750. , 5000. , 5500. , 6755.69846645 };
138  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
139  static const double c[] = { 13.49120468, 39.52387312, -5.50210824, 83.25407145,
140  92.48704203, 87.44313674, 121.88717453, 146.24929624,
141  178.96255878, 205.85276253, 235.74854479, 262.40281119,
142  287.51330043, 314.39461196, 339.95264153, 362.62280653,
143  427.29130168, 468.13556834, 558.41447367, 690.04158918,
144  610.01095483, 701.48353634 };
145  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
146  return BSpline(x, knots, coeffs, 3);
147  }
148 
149  double getCeSigma(double x)
150  {
151  static const double t[] = { -101.45678401, 1521.56069445, 2405.56898675, 3105.11956325,
152  3370.44648263, 3628.78839557, 3715.55057002, 3785.86008755,
153  3918.7967475 , 4070.43012082, 4093.84957495, 4101.69012082,
154  4160.78677462, 4677.60910462, 4711.07709032, 4864.08796619,
155  5895.47312434, 6065.52312434};
156  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
157  static const double c[] = { -58.42692309, 122.54642736, 334.89738829, 687.39838556,
158  956.19509815, 1256.96223705, 1449.1703743 , 1613.29719523,
159  1733.37464689, 1855.41536503, 1867.27344587, 2041.3608888 ,
160  2141.90965386, 2969.56208685, 2354.8500548 , 3525.00539393,
161  9805.42066704, 34273.37238895, 17841.04620772, 18300.542827 };
162  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
163  return BSpline(x, knots, coeffs, 3);
164  }
165 
166  double getGeoEcorr(double x)
167  {
168  static const double t[] = { -200.0, -187.878787879, -181.818181818, -175.757575758,
169  -169.696969697, -163.636363636, -157.575757576, -151.515151515, -145.454545455,
170  -139.393939394, -133.333333333, -127.272727273, -121.212121212, -115.151515152,
171  -109.090909091, -103.03030303, -96.9696969697, -90.9090909091, -84.8484848485, -78.7878787879,
172  -72.7272727273, -66.6666666667, -60.6060606061, -54.5454545455, -48.4848484848,
173  -42.4242424242, -36.3636363636, -30.303030303, -24.2424242424, -18.1818181818, -12.1212121212,
174  -6.06060606061, 0.0, 6.06060606061, 12.1212121212, 18.1818181818, 24.2424242424, 30.303030303,
175  36.3636363636, 42.4242424242, 48.4848484848, 54.5454545455, 60.6060606061, 66.6666666667,
176  72.7272727273, 78.7878787879, 84.8484848485, 90.9090909091, 96.9696969697, 103.03030303,
177  109.090909091, 115.151515152, 121.212121212, 127.272727273, 133.333333333, 139.393939394,
178  145.454545455, 151.515151515, 157.575757576, 163.636363636, 169.696969697, 175.757575758,
179  181.818181818, 187.878787879, 193.939393939, 200.0, 206.060606061, 212.121212121,
180  218.181818182, 224.242424242, 230.303030303, 236.363636364, 242.424242424, 248.484848485,
181  254.545454545, 260.606060606, 266.666666667, 272.727272727, 278.787878788, 284.848484848,
182  290.909090909, 296.96969697, 303.03030303, 309.090909091, 315.151515152, 321.212121212,
183  327.272727273, 333.333333333, 339.393939394, 345.454545455, 351.515151515, 357.575757576,
184  363.636363636, 369.696969697, 375.757575758, 381.818181818, 387.878787879, 393.939393939,
185  400.0, 405.0, 405.555555556, 406.111111111, 406.666666667, 407.222222222, 407.777777778,
186  408.333333333, 408.888888889, 409.444444444, 410.0, 410.555555556, 411.111111111,
187  411.666666667, 412.222222222, 412.777777778, 413.333333333, 413.888888889, 414.444444444,
188  415.0, 415.555555556, 416.111111111, 416.666666667, 417.222222222, 417.777777778,
189  418.333333333, 418.888888889, 419.444444444, 420.0, 420.555555556, 421.111111111,
190  421.666666667, 422.222222222, 422.777777778, 423.333333333, 423.888888889, 424.444444444,
191  425.0, 425.555555556, 426.111111111, 426.666666667, 427.222222222, 427.777777778,
192  428.333333333, 428.888888889, 429.444444444, 430.0, 430.555555556, 431.111111111,
193  431.666666667, 432.222222222, 432.777777778, 433.333333333, 433.888888889, 434.444444444,
194  435.0, 435.555555556, 436.111111111, 436.666666667, 437.222222222, 437.777777778,
195  438.333333333, 438.888888889, 439.444444444, 440.0, 440.555555556, 441.111111111,
196  441.666666667, 442.222222222, 442.777777778, 443.333333333, 443.888888889, 444.444444444,
197  445.0, 445.555555556, 446.111111111, 446.666666667, 447.222222222, 447.777777778,
198  448.333333333, 448.888888889, 449.444444444, 450.0, 450.555555556, 451.111111111,
199  451.666666667, 452.222222222, 452.777777778, 453.333333333, 453.888888889, 454.444444444,
200  455.0, 455.555555556, 456.111111111, 456.666666667, 457.222222222, 457.777777778,
201  458.333333333, 458.888888889, 459.444444444, 460.0, 470.0, 505.656565657, 541.313131313,
202  576.96969697, 612.626262626, 648.282828283, 683.939393939, 719.595959596, 755.252525253,
203  790.909090909, 826.565656566, 862.222222222, 897.878787879, 933.535353535, 969.191919192,
204  1004.84848485, 1040.50505051, 1076.16161616, 1111.81818182, 1147.47474747, 1183.13131313,
205  1218.78787879, 1254.44444444, 1290.1010101, 1325.75757576, 1361.41414141, 1397.07070707,
206  1432.72727273, 1468.38383838, 1504.04040404, 1539.6969697, 1575.35353535, 1611.01010101,
207  1646.66666667, 1682.32323232, 1717.97979798, 1753.63636364, 1789.29292929, 1824.94949495,
208  1860.60606061, 1896.26262626, 1931.91919192, 1967.57575758, 2003.23232323, 2038.88888889,
209  2074.54545455, 2110.2020202, 2145.85858586, 2181.51515152, 2217.17171717, 2252.82828283,
210  2288.48484848, 2324.14141414, 2359.7979798, 2395.45454545, 2431.11111111, 2466.76767677,
211  2502.42424242, 2538.08080808, 2573.73737374, 2609.39393939, 2645.05050505, 2680.70707071,
212  2716.36363636, 2752.02020202, 2787.67676768, 2823.33333333, 2858.98989899, 2894.64646465,
213  2930.3030303, 2965.95959596, 3001.61616162, 3037.27272727, 3072.92929293, 3108.58585859,
214  3144.24242424, 3179.8989899, 3215.55555556, 3251.21212121, 3286.86868687, 3322.52525253,
215  3358.18181818, 3393.83838384, 3429.49494949, 3465.15151515, 3500.80808081, 3536.46464646,
216  3572.12121212, 3607.77777778, 3643.43434343, 3679.09090909, 3714.74747475, 3750.4040404,
217  3786.06060606, 3821.71717172, 3857.37373737, 3893.03030303, 3928.68686869, 4000.0 };
218  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
219  static const double c[] = { 4.09579894127, 3.95814473798, 3.78869310835, 3.62182870224,
220  3.52333396852, 3.44930116982, 3.39853812511, 3.37004718055, 3.36613087593, 3.27958211762,
221  3.24711391337, 3.23104096523, 3.24157759279, 3.2771187285, 3.34014359932, 3.43287033812,
222  3.55830046632, 3.71972378416, 3.92045692834, 4.16326499461, 4.44943837368, 4.77774665142,
223  5.14226438824, 5.53326618721, 5.93685895075, 6.33354150807, 6.69973125636, 7.01032448944,
224  7.24243064894, 7.37924796732, 7.41295259226, 7.34568883638, 7.18845091121, 6.95839122567,
225  6.67551251138, 6.35966206781, 6.02839152064, 5.69582618326, 5.37238795788, 5.06508757689,
226  4.77810706649, 4.51346343496, 4.27162771941, 4.05204051047, 3.85350942738, 3.67449753973,
227  3.51332239939, 3.36828603327, 3.23775500519, 3.12020612565, 3.01424679054, 2.9186279676,
228  2.83222708497, 2.7540610952, 2.68326017336, 2.61906331421, 2.56080510984, 2.50790542096,
229  2.45985927161, 2.41622790323, 2.37663080539, 2.3407387498, 2.30826774855, 2.27897386308,
230  2.25264877734, 2.22911606907, 2.208228059, 2.18986335259, 2.17392435522, 2.16033719872,
231  2.14904478176, 2.14002747103, 2.13322242738, 2.12882306666, 2.12564698083, 2.12201854712,
232  2.11659479604, 2.10814411591, 2.09561156502, 2.07811922302, 2.05497933105, 2.02575567259,
233  1.99011688437, 1.94731353305, 1.89712016167, 1.83975298662, 1.77587590106, 1.70663059988,
234  1.63331965101, 1.55831076986, 1.4880592269, 1.42559586641, 1.37403563446, 1.32979694416,
235  1.29794465966, 1.27637988648, 1.25716851693, 1.24380954831, 1.23431580706, 1.22877429087,
236  1.22691295001, 1.22633190091, 1.22625708453, 1.22597486687, 1.2265865964, 1.22397704742,
237  1.2200475633, 1.21672310841, 1.2134646402, 1.21039658889, 1.20746811043, 1.20467765397,
238  1.20201203592, 1.19946290309, 1.19702183843, 1.19468165525, 1.19243566184, 1.19027794408,
239  1.18820322767, 1.18620654321, 1.18428339196, 1.18242971297, 1.18064179239, 1.17891600774,
240  1.17724959741, 1.17563891971, 1.17408216555, 1.17257592298, 1.17111842564, 1.16970713015,
241  1.16834017981, 1.16701555272, 1.16573146854, 1.16448627925, 1.16327822823, 1.1621062931,
242  1.16096818419, 1.1598633316, 1.15879036712, 1.15774777453, 1.15673468874, 1.15574989515,
243  1.15479244601, 1.15386128518, 1.15295567141, 1.15207416768, 1.15121713069, 1.15038193538,
244  1.14956997758, 1.14877874166, 1.14800836415, 1.14725788309, 1.14652695962, 1.14581371366,
245  1.14511558509, 1.14443000596, 1.14375494903, 1.1430883285, 1.14242833366, 1.14177328466,
246  1.14112159039, 1.1404718697, 1.13982281972, 1.1391731397, 1.13852176773, 1.13786831308,
247  1.13721074488, 1.13654863867, 1.13588125486, 1.13520766702, 1.13452736293, 1.13383972348,
248  1.133144317, 1.13244047239, 1.13172853577, 1.13100742358, 1.13027735118, 1.12953799486,
249  1.12878931292, 1.12803116572, 1.12726364402, 1.12648636764, 1.12570018097, 1.12490435674,
250  1.12409943828, 1.12328543091, 1.12246249816, 1.12163095977, 1.12079070775, 1.11994283996,
251  1.11908465193, 1.1182241988, 1.11734341384, 1.11649933489, 1.1154871854, 1.11507125651,
252  1.11454695471, 1.11403505851, 1.11350381317, 1.10988194726, 1.09258636011, 1.02657968908,
253  0.99999799562, 1.02372864707, 1.03501526097, 1.03933604676, 1.03789067022, 1.03508270779,
254  1.03206701319, 1.0303774479, 1.02945044936, 1.02883845486, 1.0281511747, 1.02710878311,
255  1.02574021888, 1.02415999223, 1.0225106602, 1.02091588421, 1.01943647656, 1.01807294318,
256  1.01683382549, 1.01572353563, 1.01473937946, 1.0138658263, 1.01309016187, 1.01240209451,
257  1.01178984348, 1.01124491706, 1.01075883617, 1.01032437874, 1.00993513598, 1.0095853978,
258  1.00927021739, 1.00898533932, 1.0087266415, 1.00849109498, 1.00827536522, 1.00807707971,
259  1.00789378823, 1.00772346284, 1.00756427991, 1.00741458382, 1.00727294338, 1.00713808878,
260  1.00700892182, 1.00688432631, 1.0067634952, 1.0066455815, 1.00652990595, 1.00641579296,
261  1.00630274815, 1.00619024639, 1.0060778877, 1.00596529875, 1.00585215927, 1.00573822944,
262  1.00562321621, 1.00550723156, 1.00538947765, 1.00527052728, 1.00514998281, 1.00502784383,
263  1.00490410573, 1.0047787678, 1.00465185826, 1.00452342566, 1.00439354928, 1.0042623099,
264  1.00412981708, 1.00399619245, 1.00386158205, 1.00372610402, 1.00358997035, 1.00345327396,
265  1.00331635713, 1.00317918347, 1.00304216319, 1.00290537135, 1.0027691235, 1.00263357544,
266  1.00249898172, 1.00236561414, 1.00223363118, 1.00210337885, 1.00197504385, 1.00184891781,
267  1.00172525483, 1.00160431226, 1.00148632168, 1.00137172692, 1.00126042581, 1.00115303192,
268  1.00104961784, 1.00095053713, 1.00085586797, 1.00076599552, 1.00068100448, 1.00057451822,
269  1.00050253578, 1.00045819906 };
270  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
271  return BSpline(x, knots, coeffs, 3);
272  }
273 
274  double getCeEcorr(double x)
275  {
276  static const double t[] = { -200.0, -187.878787879, -181.818181818, -175.757575758,
277  -169.696969697, -163.636363636, -157.575757576, -151.515151515, -145.454545455,
278  -139.393939394, -133.333333333, -127.272727273, -121.212121212, -115.151515152,
279  -109.090909091, -103.03030303, -96.9696969697, -90.9090909091, -84.8484848485, -78.7878787879,
280  -72.7272727273, -66.6666666667, -60.6060606061, -54.5454545455, -48.4848484848,
281  -42.4242424242, -36.3636363636, -30.303030303, -24.2424242424, -18.1818181818, -12.1212121212,
282  -6.06060606061, 0.0, 6.06060606061, 12.1212121212, 18.1818181818, 24.2424242424, 30.303030303,
283  36.3636363636, 42.4242424242, 48.4848484848, 54.5454545455, 60.6060606061, 66.6666666667,
284  72.7272727273, 78.7878787879, 84.8484848485, 90.9090909091, 96.9696969697, 103.03030303,
285  109.090909091, 115.151515152, 121.212121212, 127.272727273, 133.333333333, 139.393939394,
286  145.454545455, 151.515151515, 157.575757576, 163.636363636, 169.696969697, 175.757575758,
287  181.818181818, 187.878787879, 193.939393939, 200.0, 206.060606061, 212.121212121,
288  218.181818182, 224.242424242, 230.303030303, 236.363636364, 242.424242424, 248.484848485,
289  254.545454545, 260.606060606, 266.666666667, 272.727272727, 278.787878788, 284.848484848,
290  290.909090909, 296.96969697, 303.03030303, 309.090909091, 315.151515152, 321.212121212,
291  327.272727273, 333.333333333, 339.393939394, 345.454545455, 351.515151515, 357.575757576,
292  363.636363636, 369.696969697, 375.757575758, 381.818181818, 387.878787879, 393.939393939,
293  400.0, 405.0, 405.555555556, 406.111111111, 406.666666667, 407.222222222, 407.777777778,
294  408.333333333, 408.888888889, 409.444444444, 410.0, 410.555555556, 411.111111111,
295  411.666666667, 412.222222222, 412.777777778, 413.333333333, 413.888888889, 414.444444444,
296  415.0, 415.555555556, 416.111111111, 416.666666667, 417.222222222, 417.777777778,
297  418.333333333, 418.888888889, 419.444444444, 420.0, 420.555555556, 421.111111111,
298  421.666666667, 422.222222222, 422.777777778, 423.333333333, 423.888888889, 424.444444444,
299  425.0, 425.555555556, 426.111111111, 426.666666667, 427.222222222, 427.777777778,
300  428.333333333, 428.888888889, 429.444444444, 430.0, 430.555555556, 431.111111111,
301  431.666666667, 432.222222222, 432.777777778, 433.333333333, 433.888888889, 434.444444444,
302  435.0, 435.555555556, 436.111111111, 436.666666667, 437.222222222, 437.777777778,
303  438.333333333, 438.888888889, 439.444444444, 440.0, 440.555555556, 441.111111111,
304  441.666666667, 442.222222222, 442.777777778, 443.333333333, 443.888888889, 444.444444444,
305  445.0, 445.555555556, 446.111111111, 446.666666667, 447.222222222, 447.777777778,
306  448.333333333, 448.888888889, 449.444444444, 450.0, 450.555555556, 451.111111111,
307  451.666666667, 452.222222222, 452.777777778, 453.333333333, 453.888888889, 454.444444444,
308  455.0, 455.555555556, 456.111111111, 456.666666667, 457.222222222, 457.777777778,
309  458.333333333, 458.888888889, 459.444444444, 460.0, 470.0, 505.656565657, 541.313131313,
310  576.96969697, 612.626262626, 648.282828283, 683.939393939, 719.595959596, 755.252525253,
311  790.909090909, 826.565656566, 862.222222222, 897.878787879, 933.535353535, 969.191919192,
312  1004.84848485, 1040.50505051, 1076.16161616, 1111.81818182, 1147.47474747, 1183.13131313,
313  1218.78787879, 1254.44444444, 1290.1010101, 1325.75757576, 1361.41414141, 1397.07070707,
314  1432.72727273, 1468.38383838, 1504.04040404, 1539.6969697, 1575.35353535, 1611.01010101,
315  1646.66666667, 1682.32323232, 1717.97979798, 1753.63636364, 1789.29292929, 1824.94949495,
316  1860.60606061, 1896.26262626, 1931.91919192, 1967.57575758, 2003.23232323, 2038.88888889,
317  2074.54545455, 2110.2020202, 2145.85858586, 2181.51515152, 2217.17171717, 2252.82828283,
318  2288.48484848, 2324.14141414, 2359.7979798, 2395.45454545, 2431.11111111, 2466.76767677,
319  2502.42424242, 2538.08080808, 2573.73737374, 2609.39393939, 2645.05050505, 2680.70707071,
320  2716.36363636, 2752.02020202, 2787.67676768, 2823.33333333, 2858.98989899, 2894.64646465,
321  2930.3030303, 2965.95959596, 3001.61616162, 3037.27272727, 3072.92929293, 3108.58585859,
322  3144.24242424, 3179.8989899, 3215.55555556, 3251.21212121, 3286.86868687, 3322.52525253,
323  3358.18181818, 3393.83838384, 3429.49494949, 3465.15151515, 3500.80808081, 3536.46464646,
324  3572.12121212, 3607.77777778, 3643.43434343, 3679.09090909, 3714.74747475, 3750.4040404,
325  3786.06060606, 3821.71717172, 3857.37373737, 3893.03030303, 3928.68686869, 4000.0 };
326  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
327  static const double c[] = { 1.17248741329, 1.17159509349, 1.17021760552, 1.16832452168,
328  1.16688005928, 1.16541357182, 1.16392686764, 1.16242122524, 1.16089758654, 1.15935672261,
329  1.15779935581, 1.15622624747, 1.15463825518, 1.15303636618, 1.15142171256, 1.14979557385,
330  1.14815937098, 1.14651465512, 1.14486309358, 1.14320645407, 1.1415465882, 1.13988541437,
331  1.13822490018, 1.13656704419, 1.13491385693, 1.13326734114, 1.13162947115, 1.13000217158,
332  1.12838729559, 1.12678660297, 1.12520173849, 1.12363421107, 1.12208537438, 1.12055640932,
333  1.11904830921, 1.11756186808, 1.11609767266, 1.11465609837, 1.11323730964, 1.1118412645,
334  1.1104677234, 1.10911626191, 1.10778628691, 1.1064770555, 1.10518769613, 1.10391723101,
335  1.1026645991, 1.10142867894, 1.10020831053, 1.09900231589, 1.09780951757, 1.09662875495,
336  1.0954588981, 1.09429885906, 1.09314760057, 1.09200414244, 1.09086756564, 1.08973701441,
337  1.08861169661, 1.08749088282, 1.08637390331, 1.08526014765, 1.08414905069, 1.0830401261,
338  1.08193282054, 1.08082685231, 1.07972173557, 1.07861713607, 1.07751271358, 1.07640815898,
339  1.07530317934, 1.07419749872, 1.07309084976, 1.07198298853, 1.07087362543, 1.06976229978,
340  1.06864862248, 1.06753221826, 1.06641276676, 1.06528994426, 1.06416345938, 1.06303303303,
341  1.06189840029, 1.06075930714, 1.05961550919, 1.05846677044, 1.05731286226, 1.05615356342,
342  1.05498865754, 1.05381794095, 1.05264118756, 1.05145815123, 1.05026860659, 1.04907233963,
343  1.04786914975, 1.04665884898, 1.04544126198, 1.04421622597, 1.0429835908, 1.04181557392,
344  1.04101853101, 1.04059799177, 1.04048323233, 1.0403684059, 1.04025351238, 1.04013855171,
345  1.0400235238, 1.03990842858, 1.03979326596, 1.03967803588, 1.03956273825, 1.039447373,
346  1.03933194006, 1.03921643936, 1.03910087082, 1.03898523437, 1.03886952995, 1.03875375748,
347  1.0386379169, 1.03852200814, 1.03840603114, 1.03828998582, 1.03817387213, 1.03805769,
348  1.03794143936, 1.03782512017, 1.03770873233, 1.03759227594, 1.03747575028, 1.03735915781,
349  1.03724248906, 1.03712577909, 1.03700889681, 1.03689233119, 1.03677425713, 1.03666148574,
350  1.03652859734, 1.03647045926, 1.03613302153, 1.03683761642, 1.03365296457, 1.0449828508,
351  1.00214341491, 1.0136517575, 1.01059730539, 1.0114449521, 1.0112469551, 1.01132902545,
352  1.01133595742, 1.01136292444, 1.01138440336, 1.01140732715, 1.01142965643, 1.01145207716,
353  1.01147441225, 1.01149659755, 1.01151874666, 1.01154079308, 1.01156276351, 1.01158464331,
354  1.01160643123, 1.01162813428, 1.01164974664, 1.01167126877, 1.01169270014, 1.01171402917,
355  1.01173528592, 1.01175643526, 1.01177749143, 1.01179846011, 1.01181931681, 1.01184008778,
356  1.01186076095, 1.01188133378, 1.01190180736, 1.01192218157, 1.0119424565, 1.01196262632,
357  1.01198270008, 1.01200265969, 1.01202255057, 1.01204228519, 1.01206194636, 1.01208149263,
358  1.01210093126, 1.0121202808, 1.01213949234, 1.01215861588, 1.01217762174, 1.01219652453,
359  1.01221530582, 1.01223398695, 1.01225254768, 1.01227103915, 1.01228918678, 1.01230813354,
360  1.01232355356, 1.01235156281, 1.01235955046, 1.01249157308, 1.01292215263, 1.01354761021,
361  1.01401718484, 1.01429901759, 1.01446365434, 1.01458161403, 1.01469928711, 1.01485456397,
362  1.01503743576, 1.01522416515, 1.01539298556, 1.01551990397, 1.01557991501, 1.01554677101,
363  1.01539371619, 1.0150934031, 1.01461879578, 1.01394341921, 1.01304181242, 1.01188993492,
364  1.01046545744, 1.00874800102, 1.00671930301, 1.00436325956, 1.00166603765, 0.998615907577,
365  0.995203316753, 0.991420600303, 0.987262248801, 0.982724085905, 0.977804245822,
366  0.972500274262, 0.966820615778, 0.960750737617, 0.954310731513, 0.947496919269,
367  0.940317271328, 0.932778212118, 0.924887881518, 0.91665585042, 0.908092200656, 0.89920892333,
368  0.890017570325, 0.88053403005, 0.870770779421, 0.860744499424, 0.850471451621, 0.83996873897,
369  0.829254221802, 0.818346175119, 0.80726329626, 0.796024525776, 0.784648938313, 0.773155620944,
370  0.76156354159, 0.749891454836, 0.738157764473, 0.726380487174, 0.714577088157, 0.702764450704,
371  0.690958886861, 0.679175799766, 0.667430103811, 0.655735774222, 0.644105951423,
372  0.632553086439, 0.621088672858, 0.609723434709, 0.598467310765, 0.587329287617,
373  0.576317706003, 0.565440106057, 0.554703156508, 0.544112965094, 0.533674800998,
374  0.523393388487, 0.513272786428, 0.503316405571, 0.493527154979, 0.483907383486,
375  0.474458950635, 0.465183253611, 0.45608126232, 0.447153504999, 0.438400292467, 0.429821073933,
376  0.421416172907, 0.413184128601, 0.40512399896, 0.397235075951, 0.389515330909, 0.381963353722,
377  0.374577278565, 0.36735513376, 0.360294813063, 0.353394100519, 0.346650623364, 0.337865804291,
378  0.331505676562, 0.327363727255 };
379  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
380  return BSpline(x, knots, coeffs, 3);
381  }
382 
383  double getbCe(const double k, const double dxmax)
384  {
385  if (k < 1e-5) {
386  return 146.92691815 - 0.25112664 * dxmax;
387  } else {
388  return 55.55667917 + 0.32392104 * dxmax;
389  }
390  }
391 
392  double getRcutCe(const double k, const double dxmax)
393  {
394  if (k < 1e-5) {
395  return 0.;
396  } else {
397  double b = getbCe(k, dxmax);
398  const double p0 = 2.90571462e+01;
399  const double p1 = 1.97413284e-01;
400  const double p2 = 1.80588511e-03;
401  return 0.5 * (-p1 + pow((4 * b - 4 * p0) * p2 + p1 * p1, 0.5)) / p2;
402  }
403  }
404 
405 
406  double get_k_ce(const double dxmax)
407  {
408  const double a = 5.80505613e+02;
409  const double b = -1.76588481e+00;
410  const double e = 3.12029983e+00;
411  const double f = 3.73038601e-03;
412  double g = dxmax - a;
413  double res = b + (e - b) / (1 + exp(-f * g));
414  if (res < 0) {
415  res = 0;
416  }
417 
418  static const double t[] = {-101.45678401, 3986.00008755};
419  std::vector<double> knots(t, t + sizeof(t) / sizeof(t[0]));
420  static const double c[] = {-3.11958152, 7.23549875, 1.6947458 , 2.52750041 };
421  std::vector<double> coeffs(c, c + sizeof(c) / sizeof(c[0]));
422  return BSpline(dxmax, knots, coeffs, 3);
423  }
424  double geta(double rho)
425  {
426  const double average_density = 648.18353008270035;
427  double a = -0.23604683 + 0.43426141 * exp(1.11141046 * 1e-3 * (rho - average_density));
428  return a;
429  }
430 
431  double getp(double r, double rcut, const double p2)
432  {
433  rcut = std::max(1., fabs(rcut));
434  double b = 1e-3 * p2;
435  double p_geo = 2. * pow(rcut, b);
436  if (r <= rcut) {
437  return 2.;
438  } else {
439  return p_geo * pow(r, (-1. * b));
440  }
441  }
442 
443  double LDF_vB_parts(const double r, const double sigma, const double R, const double p)
444  {
445  return exp(-1. * pow(fabs(r - R) / (pow(2, 0.5) * sigma), p));
446  }
447 
448  double LDFvB(const double r, const double sigma, const double R, const double E, const double p)
449  {
450  if (R < 0) {
451  double norm = fabs(
452  sigma * M_PI * pow(2, 0.5)
453  * ((erfc(-1. * R * pow(2, 0.5) / (2. * sigma))) * pow(M_PI, 0.5) * R
454  + pow(2, 0.5) * sigma * exp(-R * R / 2. / pow(sigma, 2))));
455  return E / norm * LDF_vB_parts(r, sigma, R, p);
456  } else {
457  double norm = 1. / pow(sigma, 2) * 0.5 * exp(0.5 * R * R / pow(sigma, 2)) * sigma
458  / ((erf(0.5 * R * pow(2, 0.5) / sigma) * pow(M_PI, 0.5) * pow(2, 0.5)
459  * exp(0.5 * R * R / pow(sigma, 2)) * R + 2 * sigma) * M_PI);
460  return E * norm * (LDF_vB_parts(r, sigma, R, p) + LDF_vB_parts(r, sigma, -R, p));
461  }
462  }
463 
464  double my_gamma2(const double xx, const double E, const double sigma, const double k,
465  const double rcut, const double b)
466  {
467  const double p = getp(fabs(xx), rcut, b);
468  if (k < 0) {
469  return NAN;
470  }
471  const double norm = (k + 1.) / (pow(2., k) * pow(2. * k + 2, (-0.5 * k))) / (pow(sigma, k + 2.))
472  / (2 * M_PI) / tgamma(0.5 * k + 1);
473  return norm * E * pow(fabs(xx), k) * exp(-(pow(fabs(xx), p) / (p / (k + 1.) * pow(sigma, p))));
474  }
475 
476  double LDFCeDxmax(const double r, const double E, const double dxmax)
477  {
478  const double k = get_k_ce(dxmax);
479  const double sigma = getCeSigma(dxmax);
480  const double Ecorr = getCeEcorr(dxmax);
481  const double rcut = getRcutCe(k, dxmax);
482  const double b = getbCe(k, dxmax);
483  return my_gamma2(r, E, sigma, k, rcut, b) / Ecorr;
484  }
485 
486  double LDFGeoDxmax(double r, double Erad, double dxmax)
487  {
488  double rcut = getRcutGeo(dxmax);
489  double b = getBGeo(dxmax);
490  double p = getp(r, rcut, b);
491  double R = getGeoR(dxmax);
492  double sigma = getGeoSigma(dxmax);
493  double Ecorr = getGeoEcorr(dxmax);
494  const double res = LDFvB(r, sigma, R, Erad, p) / Ecorr;
495  return res;
496  }
497 
498  std::pair<double, double> LDFGeoCe(const double x, const double y, const double Erad,
499  const double dxmax, const double coreX, const double coreY,
500  const double zenith, const double sinalpha)
501  {
502  // calculate geomagnetic and charge-excess radiation energies from full radiation energy, Dxmax and zenith angles.
503  const double Xatm2 = getAtmosphere(1564.);
504  const double h2 = getVerticalHeight(Xatm2 - dxmax * cos(zenith));
505  const double rho = getDensity(h2);
506 
507  const double a = geta(rho);
508  //double sinalpha = get_sine_angle_to_lorentz_force(zenith, azimuth,
509  // magnetic_field_vector);
510  //double sinalpha = 1;
511  //std::cout << " a = " << a << " sinalpha = " << sinalpha << std::endl;
512  const double Egeo = Erad / (1 + pow(a / sinalpha, 2));
513  const double Ece = Erad - Egeo;
514 
515  const double x2 = x - coreX;
516  const double y2 = y - coreY;
517 
518  const double r = pow(x2 * x2 + y2 * y2, 0.5);
519 
520  // get energy fluence of geomagnetic and charge-excess component
521  const double fce = LDFCeDxmax(r, Ece, dxmax);
522  const double fgeo = LDFGeoDxmax(r, Egeo, dxmax);
523 
524  // combine the two emission processes depending on the position relative to the shower axis
525  const double az = atan2(y2, x2);
526  const double fvB = pow(pow(fgeo, 0.5) + pow(fce, 0.5) * cos(az), 2);
527  const double fvvB = fce * pow(sin(az), 2);
528  //const double f = fvB + fvvB;
529 
530  //std::cout << "geoceLDF = " << x << ", " << y << ", " << r << ", " << Erad << ", " << Egeo << ", "
531  // << Ece << ", " << dxmax << ", " << coreX << ", " << coreY << ", " << zenith << ", "
532  // << sinalpha << ", " << f << ", " << fvB << ", " << fvvB << ", " << fgeo << ", " << fce
533  // << std::endl;
534 
535  return std::pair<double, double>(fvB, fvvB);
536  }
537 
538  const double a[] = { -133.13151125e4, -13.973209265e4, 0.8378263431e4, 3.111742176, 0.01128292e4 };
539  const double b[] = { 1176.9833473e4, 1244.234531e4, 1464.0120855e4, 622.11207419e4, 1e4 };
540  const double c[] = { 954151.404e-2, 692708.89816e-2, 615439.43936e-2, 747969.08133e-2, 1.e7 };
541  const double layers[] = { 9.5e3, 15.5e3, 36.5e3, 100e3 };
542  const double h_max = 112829.2; // height above sea level where the mass overburden vanishes
543 
544  double getAtmosphere(double h)
545  {
546  double result = 0;
547  if (h < layers[0]) {
548  result = a[0] + b[0] * exp(-1 * h / c[0]);
549  } else if (h < layers[1]) {
550  result = a[1] + b[1] * exp(-1 * h / c[1]);
551  } else if (h < layers[2]) {
552  result = a[2] + b[2] * exp(-1 * h / c[2]);
553  } else if (h < layers[3]) {
554  result = a[3] + b[3] * exp(-1 * h / c[3]);
555  } else if (h < h_max) {
556  result = a[4] - b[4] * h / c[4];
557  } else {
558  result = 0;
559  }
560  return result * 1e-4;
561  }
562 
563  double getVerticalHeight(double at)
564  {
565  at = at * 1e4;
566  double h = 0.;
567  int i = 0;
568  if (at > getAtmosphere(layers[0])) {
569  i = 0;
570  } else if (at > getAtmosphere(layers[1])) {
571  i = 1;
572  } else if (at > getAtmosphere(layers[2])) {
573  i = 2;
574  } else if (at > getAtmosphere(layers[3])) {
575  i = 3;
576  } else {
577  i = 4;
578  }
579  if (i == 4) {
580  h = -1. * c[i] * (at - a[i]) / b[i];
581  } else {
582  h = -1. * c[i] * log((at - a[i]) / b[i]);
583  }
584  return h;
585  }
586 
587  double getDensity(double h)
588  {
589  double result = 0;
590  if (h < layers[0]) {
591  result = b[0] * exp(-1 * h / c[0]) / c[0];
592  } else if (h < layers[1]) {
593  result = b[1] * exp(-1 * h / c[1]) / c[1];
594  } else if (h < layers[2]) {
595  result = b[2] * exp(-1 * h / c[2]) / c[2];
596  } else if (h < layers[3]) {
597  result = b[3] * exp(-1 * h / c[3]) / c[3];
598  } else if (h < h_max) {
599  result = b[4] / c[4];
600  } else {
601  result = 0;
602  }
603  return result;
604  }
605 
606 
607  LDFLikelihoodFunction::LDFLikelihoodFunction(const FitConfig fitconfig,
608  const EventFitData eventData,
609  const std::vector<StationFitData>& stationData) :
610  fFitConfig(fitconfig),
611  fEventData(eventData),
612  fStationData(stationData),
613  fTheErrorDef(1)
614  {
615  }
616 
617  /*
618  A different LDF model can be implemented easily by defining a new Likelihood of Chi2 function.
619  This contribution to the combined neg. log. Likelihood needs to be added in the station loop to
620  the variable "negLogLikelihood". Just follow the example for LDFModel==1.
621  The parameters of the LDF function needs to be added to the minimizer in RdLDFFitter.cc
622  Parameters are:
623  0, 1, 2 = (coreX, coreY, coreZ) in local coordinate system
624  3 = a (charge excess strength)
625  4, 5, 6 = Scintillator LDF parameters
626  7, 8,... = optional LDF parameters
627  */
628 
629  double LDFLikelihoodFunction::operator()(const std::vector<double>& pars) const
630  {
631  const double Erad = pars[0] * utl::MeV; // the fit parameter is the energy in MeV to be in a similar order of magnitude as the other fit parameters. The LDF function however requires the energy to be in eV (standard auger units)
632  const double dxmax = pars[1];
633  const double coreX = pars[2];
634  const double coreY = pars[3];
635 
636  double negLogLikelihood = 0;
637 
638  // std::cout << "station loop " << std::endl;
639 
640  // loop over all stations
641  for (std::vector<StationFitData>::const_iterator sIt = fStationData.begin(), end = fStationData
642  .end(); sIt != end; ++sIt) {
643 
644  negLogLikelihood += GetChi2LDFModel1(sIt->x, sIt->y, Erad, dxmax, coreX, coreY, sIt->f_vB,
645  sIt->ferror_vB, sIt->f_vvB, sIt->ferror_vvB,
646  fEventData.zenith, fEventData.sinalpha);
647 
648  } // loop over stations
649  // std::cout << std::endl;
650 
651  return std::isnan(negLogLikelihood) ? std::numeric_limits<double>::max() : negLogLikelihood;
652  }
653 
654  /*double Up() const { return 0.5; }*/
655 
656  double LDFLikelihoodFunction::GetChi2LDFModel1(const double x, const double y, const double Erad,
657  const double dxmax, const double coreX,
658  const double coreY, const double f_vB,
659  const double ferror_vB, const double f_vvB,
660  const double ferror_vvB, const double zenith,
661  const double sinalpha) const
662  {
663  std::pair<double, double> tmp = LDFGeoCe(x, y, Erad, dxmax, coreX, coreY, zenith, sinalpha);
664  //std::cout << x << ", " << y << " = " << tmp.first << ", " << tmp.second << std::endl;
665  const double chi1 = pow((tmp.first - f_vB) / ferror_vB, 2);
666  const double chi2 = pow((tmp.second - f_vvB) / ferror_vvB, 2);
667  double c = fabs(f_vvB / f_vB);
668  if (c > 1) {
669  c = 1;
670  }
671 
672  const double res = chi1 + c * chi2;
673 
674  //std::cout << "(" << tmp.first << " - " << f_vB << " / " << ferror_vB << ") + (" << tmp.second
675  // << " - " << f_vvB << " / " << ferror_vvB << ") = " << res << ", " << std::endl;
676  return res;
677  }
678 
679  double LDFLikelihoodFunction::GetChi2FullEnergyFluence(const std::vector<double>& pars) const
680  {
681  const double Erad = pars[0] * utl::MeV; // the fit parameter is the energy in MeV to be in a similar order of magnitude as the other fit parameters. The LDF function however requires the energy to be in eV (standard auger units)
682  const double dxmax = pars[1];
683  const double coreX = pars[2];
684  const double coreY = pars[3];
685 
686  double negLogLikelihood = 0;
687 
688  // loop over all stations
689  for (std::vector<StationFitData>::const_iterator sIt = fStationData.begin(), end = fStationData
690  .end(); sIt != end; ++sIt) {
691  std::pair<double, double> tmp = LDFGeoCe(sIt->x, sIt->y, Erad, dxmax, coreX, coreY,
692  fEventData.zenith, fEventData.sinalpha);
693  const double ffit = tmp.first + tmp.second;
694  const double f = sIt->f_vB + sIt->f_vvB;
695  const double f_error = sIt->ferror_vB + sIt->ferror_vvB;
696  negLogLikelihood += pow((ffit - f) / f_error, 2);
697 
698  } // loop over stations
699 
700  return std::isnan(negLogLikelihood) ? std::numeric_limits<double>::max() : negLogLikelihood;
701  }
702 
703 }
RdGeoCeLDFFitter::getRcutCe
double getRcutCe(const double k, const double dxmax)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:392
RdGeoCeLDFFitter::LDFLikelihoodFunction::LDFLikelihoodFunction
LDFLikelihoodFunction(const FitConfig ldfconfig, const EventFitData eventData, const std::vector< StationFitData > &stationData)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:607
RdGeoCeLDFFitter::LDFLikelihoodFunction::GetChi2FullEnergyFluence
double GetChi2FullEnergyFluence(const std::vector< double > &pars) const
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:679
RdGeoCeLDFFitter::LDFLikelihoodFunction::GetChi2LDFModel1
double GetChi2LDFModel1(const double x, const double y, const double Erad, const double dxmax, const double coreX, const double coreY, const double f_vB, const double ferror_vB, const double f_vvB, const double ferror_vvB, const double zenith, const double sinalpha) const
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:656
RdGeoCeLDFFitter::LDFCeDxmax
double LDFCeDxmax(const double r, const double E, const double dxmax)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:476
RdGeoCeLDFFitter::get_k_ce
double get_k_ce(const double dxmax)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:406
std::pow
double pow(const double x, const unsigned int i)
Definition: MuonArrivalTime.h:17
utl::MeV
constexpr double MeV
Definition: AugerUnits.h:184
max
#define max(a, b)
Definition: DLECorrectionGG/DLECorrection.cc:30
result
const Data result[]
Definition: testReadStream.cc:56
utl::g
constexpr double g
Definition: AugerUnits.h:200
RdGeoCeLDFFitter::getbCe
double getbCe(const double k, const double dxmax)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:383
RdGeoCeLDFFitter::LDFGeoDxmax
double LDFGeoDxmax(double r, double Erad, double dxmax)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:486
RdGeoCeLDFFitter::LDFLikelihoodFunction::operator()
double operator()(const std::vector< double > &pars) const override
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:629
RdBeamFormer::norm
double norm(utl::Vector x)
Definition: RdBeamFormer/RdWaveModel.cc:39
RdGeoCeLDFFitter::LDF_vB_parts
double LDF_vB_parts(const double r, const double sigma, const double R, const double p)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:443
RdGeoCeLDFFitter::LDFvB
double LDFvB(const double r, const double sigma, const double R, const double E, const double p)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:448
RdGeoCeLDFFitter::BSpline
double BSpline(const double x, std::vector< double > t, const std::vector< double > c, const int k)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:29
RdGeoCeLDFFitter::LDFGeoCe
std::pair< double, double > LDFGeoCe(const double x, const double y, const double Erad, const double dxmax, const double coreX, const double coreY, const double zenith, const double sinalpha)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:498
RdGeoCeLDFFitter::getp
double getp(double r, double rcut, const double p2)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:431
utl::m
constexpr double m
Definition: AugerUnits.h:121
RdGeoCeLDFFitter::my_gamma2
double my_gamma2(const double xx, const double E, const double sigma, const double k, const double rcut, const double b)
Definition: RdGeoCeLDFFitter/LikelihoodFunction.cc:464
, generated on Tue Sep 26 2023.