GooFit  v2.1.3
Functions | Variables
landau.h File Reference
#include <cmath>
#include <limits>

Go to the source code of this file.

Functions

double landauPDF (const double &x, const double &x0, const double &xi)
 
double gaussPDF (const double &x, const double &mu, const double &sigma)
 
double landauGaussPDF (const double &x, const double &mu, const double &eta, const double &sigma)
 
double landau_quantile (double z, double xi)
 
double landau_quantile_c (double z, double xi)
 

Variables

const double invsq2pi = 0.3989422804014
 

Function Documentation

◆ gaussPDF()

double gaussPDF ( const double &  x,
const double &  mu,
const double &  sigma 
)
inline

Definition at line 93 of file landau.h.

References invsq2pi.

Referenced by landauGaussPDF().

93  {
94  return invsq2pi / sigma * exp(-pow((x - mu), 2) / (2 * pow(sigma, 2)));
95 }
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 at line 144 of file landau.h.

Referenced by landau_quantile_c().

144  {
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 }

◆ landau_quantile_c()

double landau_quantile_c ( double  z,
double  xi 
)

Definition at line 286 of file landau.h.

References landau_quantile().

286 { return landau_quantile(1. - z, xi); }
landau_quantile
double landau_quantile(double z, double xi)
Definition: landau.h:144

◆ landauGaussPDF()

double landauGaussPDF ( const double &  x,
const double &  mu,
const double &  eta,
const double &  sigma 
)
inline

Definition at line 97 of file landau.h.

References gaussPDF(), landauPDF(), and sigma.

97  {
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 }
sigma
Observable * sigma
Definition: pipipi0DPFit.cpp:75
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

◆ landauPDF()

double landauPDF ( const double &  x,
const double &  x0,
const double &  xi 
)
inline

Definition at line 29 of file landau.h.

Referenced by landauGaussPDF().

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 }

Variable Documentation

◆ invsq2pi

const double invsq2pi = 0.3989422804014

Definition at line 27 of file landau.h.

Referenced by gaussPDF().

Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13