Line data    Source code

       1             : // PythiaStdlib.cc is a part of the PYTHIA event generator.
       2             : // Copyright (C) 2015 Torbjorn Sjostrand.
       3             : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
       4             : // Please respect the MCnet Guidelines, see GUIDELINES for details.
       5             : 
       6             : // Function definitions (not found in the header) for the Gamma function.
       7             : 
       8             : #include "Pythia8/PythiaStdlib.h"
       9             : 
      10             : namespace Pythia8 {
      11             : 
      12             : //==========================================================================
      13             : 
      14             : // The Gamma function for real arguments, using the Lanczos approximation.
      15             : // Code based on http://en.wikipedia.org/wiki/Lanczos_approximation
      16             : 
      17             : double GammaCoef[9] = {
      18             :      0.99999999999980993,     676.5203681218851,   -1259.1392167224028,
      19             :       771.32342877765313,   -176.61502916214059,    12.507343278686905,
      20             :     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
      21             : 
      22             : double GammaReal(double x) {
      23             : 
      24             :   // Reflection formula (recursive!) for x < 0.5.
      25           0 :   if (x < 0.5) return M_PI / (sin(M_PI * x) * GammaReal(1 - x));
      26             : 
      27             :   // Iterate through terms.
      28           0 :   double z = x - 1.;
      29           0 :   double gamma = GammaCoef[0];
      30           0 :   for (int i = 1; i < 9; ++i) gamma += GammaCoef[i] / (z + i);
      31             : 
      32             :   // Answer.
      33           0 :   double t = z + 7.5;
      34           0 :   gamma *= sqrt(2. * M_PI) * pow(t, z + 0.5) * exp(-t);
      35             :   return gamma;
      36             : 
      37           0 : }
      38             : 
      39             : //==========================================================================
      40             : 
      41             : } // end namespace Pythia8