GooFit  v2.1.3
landau.h
Go to the documentation of this file.
1 // Taken from LCG ROOT MathLib
2 // License info:
3 // Authors: Andras Zsenei & Lorenzo Moneta 06/2005
4 
5 /**********************************************************************
6  * *
7  * Copyright (c) 2005 , LCG ROOT MathLib Team *
8  * *
9  * *
10  **********************************************************************/
11 
12 // Licenced under LGPL 2.1, see https://github.com/SiLab-Bonn/pyLandau/blob/master/LICENSE
13 
14 // Langaus authors:
15 // Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
16 // Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and
17 // Markus Friedl (Markus.Friedl@cern.ch)
18 //
19 // Adaption for Python by David-Leon Pohl, pohl@physik.uni-bonn.de
20 // Adaption for pybind11 by Henry Schreiner, henry.schreiner@cern.ch
21 
22 #pragma once
23 
24 #include <cmath>
25 #include <limits>
26 
27 const double invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
28 
29 inline double landauPDF(const double &x, const double &x0, const double &xi) // same algorithm is used in GSL
30 {
31  static double p1[5] = {0.4259894875, -0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
32  static double q1[5] = {1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
33 
34  static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
35  static double q2[5] = {1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
36 
37  static double p3[5] = {0.1788544503, 0.09359161662, 0.006325387654, 0.00006611667319, -0.000002031049101};
38  static double q3[5] = {1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
39 
40  static double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
41  static double q4[5] = {1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511};
42 
43  static double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
44  static double q5[5] = {1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357};
45 
46  static double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
47  static double q6[5] = {1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939};
48 
49  static double a1[3] = {0.04166666667, -0.01996527778, 0.02709538966};
50 
51  static double a2[2] = {-1.845568670, -4.284640743};
52 
53  if(xi <= 0)
54  return 0;
55  double v = (x - x0) / xi;
56  double u, ue, us, denlan;
57  if(v < -5.5) {
58  u = std::exp(v + 1.0);
59  if(u < 1e-10)
60  return 0.0;
61  ue = std::exp(-1 / u);
62  us = std::sqrt(u);
63  denlan = 0.3989422803 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
64  } else if(v < -1) {
65  u = std::exp(-v - 1);
66  denlan = std::exp(-u) * std::sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v)
67  / (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
68  } else if(v < 1) {
69  denlan = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v)
70  / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
71  } else if(v < 5) {
72  denlan = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v)
73  / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
74  } else if(v < 12) {
75  u = 1 / v;
76  denlan = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u)
77  / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
78  } else if(v < 50) {
79  u = 1 / v;
80  denlan = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u)
81  / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
82  } else if(v < 300) {
83  u = 1 / v;
84  denlan = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u)
85  / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
86  } else {
87  u = 1 / (v - v * std::log(v) / (v + 1));
88  denlan = u * u * (1 + (a2[0] + a2[1] * u) * u);
89  }
90  return denlan / xi;
91 }
92 
93 inline double gaussPDF(const double &x, const double &mu, const double &sigma) {
94  return invsq2pi / sigma * exp(-pow((x - mu), 2) / (2 * pow(sigma, 2)));
95 }
96 
97 inline double landauGaussPDF(const double &x, const double &mu, const double &eta, const double &sigma) {
98  // Numeric constants
99  double mpshift = 0.; // -0.22278298; // Landau maximum location shift in original code is wrong, since the shift
100  // does not depend on mu only
101 
102  // Control constants
103  unsigned int np = 100; // number of convolution steps
104  const unsigned int sc = 8; // convolution extends to +-sc Gaussian sigmas
105 
106  // Convolution steps have to be increased if sigma > eta * 5 to get stable solution that does not oscillate,
107  // addresses #1
108  if(sigma > 3 * eta)
109  np *= int(sigma / eta / 3.);
110  if(np > 100000) // Do not use too many convolution steps to save time
111  np = 100000;
112 
113  // Variables
114  double xx;
115  double mpc;
116  double fland;
117  double sum = 0.0;
118  double xlow, xupp;
119  double step;
120 
121  // MP shift correction
122  mpc = mu - mpshift; // * eta;
123 
124  // Range of convolution integral
125  xlow = x - sc * sigma;
126  xupp = x + sc * sigma;
127 
128  step = (xupp - xlow) / (double)np;
129 
130  // Discrete linear convolution of Landau and Gaussian
131  for(unsigned int i = 1; i <= np / 2; ++i) {
132  xx = xlow + double(i - 0.5) * step;
133  fland = landauPDF(xx, mpc, eta) / eta;
134  sum += fland * gaussPDF(x, xx, sigma);
135 
136  xx = xupp - double(i - 0.5) * step;
137  fland = landauPDF(xx, mpc, eta) / eta;
138  sum += fland * gaussPDF(x, xx, sigma);
139  }
140 
141  return (step * sum);
142 }
143 
144 double landau_quantile(double z, double xi) {
145  // LANDAU quantile : algorithm from CERNLIB G110 ranlan
146  // with scale parameter xi
147  // Converted by Rene Brun from CERNLIB routine ranlan(G110),
148  // Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
149  // Moved once again by Henry Schreiner to GooFit
150 
151  static double f[982]
152  = {0, 0, 0, 0, 0, -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
153  -2.076740, -2.050397, -2.025605, -2.002150, -1.979866, -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
154  -1.864385, -1.847451, -1.831030, -1.815083, -1.799574, -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
155  -1.714187, -1.701029, -1.688130, -1.675477, -1.663057, -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
156  -1.592811, -1.581729, -1.570806, -1.560034, -1.549407, -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
157  -1.488386, -1.478628, -1.468976, -1.459428, -1.449979, -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
158  -1.395194, -1.386356, -1.377594, -1.368906, -1.360291, -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
159  -1.310006, -1.301843, -1.293737, -1.285688, -1.277693, -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
160  -1.230800, -1.223153, -1.215550, -1.207990, -1.200474, -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
161  -1.156220, -1.148977, -1.141770, -1.134598, -1.127459, -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
162  -1.085306, -1.078388, -1.071498, -1.064636, -1.057802, -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
163  -1.017350, -1.010695, -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247, -.964749, -.958271,
164  -.951813, -.945375, -.938957, -.932558, -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
165  -.888272, -.882014, -.875773, -.869547, -.863337, -.857142, -.850963, -.844798, -.838648, -.832512,
166  -.826390, -.820282, -.814187, -.808106, -.802038, -.795982, -.789940, -.783909, -.777891, -.771884,
167  -.765889, -.759906, -.753934, -.747973, -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
168  -.706541, -.700661, -.694791, -.688931, -.683079, -.677236, -.671402, -.665576, -.659759, -.653950,
169  -.648149, -.642356, -.636570, -.630793, -.625022, -.619259, -.613503, -.607754, -.602012, -.596276,
170  -.590548, -.584825, -.579109, -.573399, -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
171  -.533592, -.527926, -.522266, -.516611, -.510961, -.505315, -.499674, -.494037, -.488405, -.482777,
172  -.477153, -.471533, -.465917, -.460305, -.454697, -.449092, -.443491, -.437893, -.432299, -.426707,
173  -.421119, -.415534, -.409951, -.404372, -.398795, -.393221, -.387649, -.382080, -.376513, -.370949,
174  -.365387, -.359826, -.354268, -.348712, -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
175  -.309863, -.304318, -.298775, -.293233, -.287692, -.282152, -.276613, -.271074, -.265536, -.259999,
176  -.254462, -.248926, -.243389, -.237854, -.232318, -.226783, -.221247, -.215712, -.210176, -.204641,
177  -.199105, -.193568, -.188032, -.182495, -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
178  -.143717, -.138173, -.132629, -.127083, -.121537, -.115989, -.110439, -.104889, -.099336, -.093782,
179  -.088227, -.082670, -.077111, -.071550, -.065987, -.060423, -.054856, -.049288, -.043717, -.038144,
180  -.032569, -.026991, -.021411, -.015828, -.010243, -.004656, .000934, .006527, .012123, .017722,
181  .023323, .028928, .034535, .040146, .045759, .051376, .056997, .062620, .068247, .073877,
182  .079511, .085149, .090790, .096435, .102083, .107736, .113392, .119052, .124716, .130385,
183  .136057, .141734, .147414, .153100, .158789, .164483, .170181, .175884, .181592, .187304,
184  .193021, .198743, .204469, .210201, .215937, .221678, .227425, .233177, .238933, .244696,
185  .250463, .256236, .262014, .267798, .273587, .279382, .285183, .290989, .296801, .302619,
186  .308443, .314273, .320109, .325951, .331799, .337654, .343515, .349382, .355255, .361135,
187  .367022, .372915, .378815, .384721, .390634, .396554, .402481, .408415, .414356, .420304,
188  .426260, .432222, .438192, .444169, .450153, .456145, .462144, .468151, .474166, .480188,
189  .486218, .492256, .498302, .504356, .510418, .516488, .522566, .528653, .534747, .540850,
190  .546962, .553082, .559210, .565347, .571493, .577648, .583811, .589983, .596164, .602355,
191  .608554, .614762, .620980, .627207, .633444, .639689, .645945, .652210, .658484, .664768,
192  .671062, .677366, .683680, .690004, .696338, .702682, .709036, .715400, .721775, .728160,
193  .734556, .740963, .747379, .753807, .760246, .766695, .773155, .779627, .786109, .792603,
194  .799107, .805624, .812151, .818690, .825241, .831803, .838377, .844962, .851560, .858170,
195  .864791, .871425, .878071, .884729, .891399, .898082, .904778, .911486, .918206, .924940,
196  .931686, .938446, .945218, .952003, .958802, .965614, .972439, .979278, .986130, .992996,
197  .999875, 1.006769, 1.013676, 1.020597, 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
198  1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778, 1.111889, 1.119016, 1.126159, 1.133316,
199  1.140490, 1.147679, 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149, 1.198451, 1.205770,
200  1.213105, 1.220457, 1.227826, 1.235211, 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
201  1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203, 1.332819, 1.340454, 1.348108, 1.355780,
202  1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216, 1.410024, 1.417851, 1.425698, 1.433565,
203  1.441453, 1.449360, 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240, 1.505293, 1.513368,
204  1.521465, 1.529583, 1.537723, 1.545885, 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
205  1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593, 1.654053, 1.662538, 1.671047, 1.679581,
206  1.688139, 1.696721, 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746, 1.757506, 1.766293,
207  1.775106, 1.783945, 1.792810, 1.801703, 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
208  1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572, 1.919830, 1.929117, 1.938434, 1.947781,
209  1.957158, 1.966566, 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659, 2.033285, 2.042943,
210  2.052633, 2.062355, 2.072110, 2.081899, 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
211  2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018, 2.212257, 2.222533, 2.232845, 2.243195,
212  2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084, 2.316701, 2.327356, 2.338051, 2.348786,
213  2.359562, 2.370377, 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138, 2.447246, 2.458397,
214  2.469591, 2.480828, 2.492110, 2.503436, 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
215  2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936, 2.654871, 2.666855, 2.678890, 2.690975,
216  2.703110, 2.715297, 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514, 2.802070, 2.814681,
217  2.827347, 2.840069, 2.852846, 2.865680, 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
218  2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263, 3.037875, 3.051551, 3.065290, 3.079095,
219  3.092965, 3.106900, 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928, 3.206340, 3.220824,
220  3.235378, 3.250005, 3.264704, 3.279477, 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
221  3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711, 3.478500, 3.494372, 3.510328, 3.526370,
222  3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450, 3.641111, 3.657863, 3.674708, 3.691646,
223  3.708680, 3.725809, 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645, 3.848470, 3.866400,
224  3.884434, 3.902574, 3.920821, 3.939176, 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
225  4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292, 4.188160, 4.208154, 4.228275, 4.248524,
226  4.268903, 4.289413, 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310, 4.436781, 4.458395,
227  4.480154, 4.502060, 4.524114, 4.546319, 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
228  4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179, 4.849511, 4.874020, 4.898710, 4.923582,
229  4.948639, 4.973883, 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411, 5.156034, 5.182864,
230  5.209903, 5.237156, 5.264625, 5.292312, 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
231  5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734, 5.673552, 5.704634, 5.735986, 5.767610,
232  5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316, 5.996967, 6.030925, 6.065194, 6.099780,
233  6.134687, 6.169921, 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465, 6.426130, 6.464166,
234  6.502578, 6.541371, 6.580553, 6.620130, 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
235  6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398, 7.173764, 7.219632, 7.266011, 7.312910,
236  7.360339, 7.408308, 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035, 7.760077, 7.812747,
237  7.866057, 7.920019, 7.974647, 8.029953, 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
238  8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955, 8.817481, 8.883916, 8.951282, 9.019600,
239  9.088889, 9.159174, 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118, 9.681054, 9.760191,
240  9.840558, 9.922186, 10.005107, 10.089353, 10.174959, 10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
241  10.719188, 10.815362, 10.913214, 11.012789, 11.114137, 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
242  11.762390, 11.877664, 11.995170, 12.114979, 12.237161, 12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
243  13.024498, 13.165570, 13.309711, 13.457026, 13.607625, 13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
244  14.587072, 14.764233, 14.945778, 15.131877, 15.322712, 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
245  16.578520, 16.808433, 17.044929, 17.288305, 17.538873, 17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
246  19.213574, 19.525133, 19.847249, 20.180480, 20.525429, 20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
247  22.880933, 23.329017, 23.795634, 24.281981, 24.789364, 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
248  28.365274, 29.068370, 29.808638, 30.589157, 31.413354, 32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
249  37.527131, 38.796172, 40.157721, 41.622399, 43.202525, 44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
250  56.123356, 59.103894};
251 
252  if(xi <= 0)
253  return 0;
254  if(z <= 0)
255  return -std::numeric_limits<double>::infinity();
256  if(z >= 1)
257  return std::numeric_limits<double>::infinity();
258 
259  double ranlan, u, v;
260  u = 1000 * z;
261  int i = int(u);
262  u -= i;
263  if(i >= 70 && i < 800) {
264  ranlan = f[i - 1] + u * (f[i] - f[i - 1]);
265  } else if(i >= 7 && i <= 980) {
266  ranlan = f[i - 1] + u * (f[i] - f[i - 1] - 0.25 * (1 - u) * (f[i + 1] - f[i] - f[i - 1] + f[i - 2]));
267  } else if(i < 7) {
268  v = std::log(z);
269  u = 1 / v;
270  ranlan = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * u) * u) / (1 + (3.41760202E1 + 4.01244582 * u) * u))
271  * (-std::log(-0.91893853 - v) - 1);
272  } else {
273  u = 1 - z;
274  v = u * u;
275  if(z <= 0.999) {
276  ranlan
277  = (1.00060006 + 2.63991156E2 * u + 4.37320068E3 * v) / ((1 + 2.57368075E2 * u + 3.41448018E3 * v) * u);
278  } else {
279  ranlan
280  = (1.00001538 + 6.07514119E3 * u + 7.34266409E5 * v) / ((1 + 6.06511919E3 * u + 6.94021044E5 * v) * u);
281  }
282  }
283  return xi * ranlan;
284 }
285 
286 double landau_quantile_c(double z, double xi) { return landau_quantile(1. - z, xi); }
landauGaussPDF
double landauGaussPDF(const double &x, const double &mu, const double &eta, const double &sigma)
Definition: landau.h:97
invsq2pi
const double invsq2pi
Definition: landau.h:27
sigma
Observable * sigma
Definition: pipipi0DPFit.cpp:75
landau_quantile
double landau_quantile(double z, double xi)
Definition: landau.h:144
gaussPDF
double gaussPDF(const double &x, const double &mu, const double &sigma)
Definition: landau.h:93
landauPDF
double landauPDF(const double &x, const double &x0, const double &xi)
Definition: landau.h:29
landau_quantile_c
double landau_quantile_c(double z, double xi)
Definition: landau.h:286
Generated on Sat Apr 21 2018 20:13:17 for GooFit by   doxygen 1.8.13