diff --git a/PDF/LeptonLeptonPDF.cc b/PDF/LeptonLeptonPDF.cc --- a/PDF/LeptonLeptonPDF.cc +++ b/PDF/LeptonLeptonPDF.cc @@ -1,114 +1,115 @@ // -*- C++ -*- // // LeptonLeptonPDF.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the LeptonLeptonPDF class. // #include "LeptonLeptonPDF.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/Maths.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Interface/ClassDocumentation.h" +#include using namespace ThePEG; IBPtr LeptonLeptonPDF::clone() const { return new_ptr(*this); } IBPtr LeptonLeptonPDF::fullclone() const { return new_ptr(*this); } bool LeptonLeptonPDF::canHandleParticle(tcPDPtr particle) const { return ( abs(particle->id()) < 20 || abs(particle->id()) > 10 ); } bool LeptonLeptonPDF::hasPoleIn1(tcPDPtr particle, tcPDPtr parton) const { return particle == parton; } cPDVector LeptonLeptonPDF::partons(tcPDPtr p) const { cPDVector ret; if ( canHandleParticle(p) ) ret.push_back(p); return ret; } double LeptonLeptonPDF:: xfl(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double l, Energy2) const { using namespace Constants; static const double gf = 3.0/4.0 - EulerGamma; static const double cut = 1.0e-300; static const double del = 1000.0*cut; if ( parton != particle || l <= cut ) return 0.0; const double x = exp(-l); const double eps = Math::exp1m(-l); const double lx = -l; const double l1x = log(eps); const double beta2 = SM().alphaEM()*(log(partonScale/sqr(particle->mass()))-1.0)/pi; double corr = 0.0; - double integral = exp(beta2*gf)*pow(cut,beta2)/Math::gamma(beta2+1.0) + double integral = exp(beta2*gf)*pow(cut,beta2)/std::tgamma(beta2+1.0) - beta2*cut*(4.0 - 3.0*beta2 + 4.0*beta2*log(cut))/4.0; if ( l < cut + del ) corr = integral/del; - return x*(exp(l1x*(beta2 - 1.0) + log(beta2) + gf*beta2)/Math::gamma(beta2+1.0) + return x*(exp(l1x*(beta2 - 1.0) + log(beta2) + gf*beta2)/std::tgamma(beta2+1.0) - 0.5*beta2*(1.0+x) + sqr(beta2)*((1.0+x)*(-4.0*l1x+3.0*lx) - 4.0*lx/eps - 5.0 - x)/8.0) + corr; } double LeptonLeptonPDF:: xfvl(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double l, Energy2 particleScale) const { return xfl(particle, parton, partonScale, l, particleScale); } double LeptonLeptonPDF:: xfvx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double x, double eps, Energy2 particleScale) const { using Math::log1m; return xfl(particle, parton, partonScale, (x < 0.5 || eps <= 0.0)? -log(x): -log1m(eps), particleScale); } double LeptonLeptonPDF:: flattenL(tcPDPtr particle, tcPDPtr, const PDFCuts & c, double z, double & jacobian) const { using namespace Constants; using Math::log1m; Energy2 scale = min(c.sMax(), c.scaleMax()); if ( c.scaleMin() > ZERO ) scale = min(scale, c.scaleMin()); double beta2 = SM().alphaEM()*(log(scale/sqr(particle->mass()))-1.0)/pi; // double zpow = pow(z, 1.0/beta2); // jacobian = zpow/(1.0/beta2+1); // return z*zpow; // double l = pow(z, 1.0/beta2); // jacobian = l/(z*beta2); // z = max(z, pow(std::numeric_limits::min()*1.5, beta2)); double eps = pow(z, 1.0/beta2); double l = -log1m(eps); jacobian *= eps/(beta2*z*(1.0 - eps)); return l; } NoPIOClassDescription LeptonLeptonPDF::initLeptonLeptonPDF; void LeptonLeptonPDF::Init() { static ClassDocumentation documentation ("Describes the distribution of leptons within leptons, ie. the energy " "loss of leptons due to radiation of soft photons."); } diff --git a/Utilities/Maths.cc b/Utilities/Maths.cc --- a/Utilities/Maths.cc +++ b/Utilities/Maths.cc @@ -1,100 +1,48 @@ // -*- C++ -*- // // Maths.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #include #include "Maths.h" - -double ThePEG::Math::atanh(double x) { -#ifndef ThePEG_HAS_ATANH - return 0.5*(ThePEG::Math::log1m(-x) - ThePEG::Math::log1m(x)); -#else - return std::atanh(x); -#endif -} +#include double ThePEG::Math::exp1m(double x) { using namespace std; #ifndef ThePEG_HAS_EXPM1 return 2.0*std::exp(x/2.0)*std::sinh(-x/2.0); #else return -expm1(x); #endif } double ThePEG::Math::log1m(double x) { using namespace std; #ifndef ThePEG_HAS_LOG1P #ifndef ThePEG_HAS_ATANH volatile double y = 1.0 - x; return log(y) - ((y - 1.0) + x)/y ; /* cancels errors with IEEE arithmetic */ #else return 2.0*atanh(x/(x-2.0)); #endif #else return log1p(-x); #endif } double ThePEG::Math::powi(double x, int i) { switch ( i ) { case 0: return 1.0; case -1: return 1/x; case -2: return 1/x/x; case -3: return 1/x/x/x; case 1: return x; case 2: return x*x; case 3: return x*x*x; default: return i > 0? powi(x, i - 1)*x: powi(x, i + 1)/x; } } - -double ThePEG::Math::gamma(double x) { - static const double b[8] = - {-0.577191652, 0.988205891, -0.897056937, 0.918206857, - -0.756704078, 0.482199394, -0.193527818, 0.035868343}; - const long nx = long(x); - const double dx = x - double(nx); - double res = 1.0; - double dxp = 1.0; - for ( int i = 0; i < 8; ++i ) { - dxp*=dx; - res += b[i]*dxp; - } - if ( x < 0.0 ) { - res /= x; - } else { - for ( long ix = 1; ix < nx; ++ix ) res *= x - double(ix); - } - return res; -} - -double ThePEG::Math::lngamma(double x) { - -#ifndef ThePEG_HAS_LGAMMA - - static const double b[6]={76.18009173, -86.50532033, 24.01409822, - -1.231739516, 0.120858003e-2, -0.536382e-5}; - x -= 1.0; - double t = x + 5.5; - t -= (x + 0.5)*log(t); - double sum = 1.0; - for (int j = 0; j < 6; j++) { - x += 1.0; - sum += b[j]/x; - } - return -t + log(2.50662827465*sum); - -#else - - return lgamma(x); - -#endif - -} - diff --git a/Utilities/Maths.h b/Utilities/Maths.h --- a/Utilities/Maths.h +++ b/Utilities/Maths.h @@ -1,349 +1,340 @@ // -*- C++ -*- // // Maths.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_Math_H #define ThePEG_Math_H #include namespace ThePEG { /** The Math namespace includes the declaration of some useful * mathematical functions. */ namespace Math { /** * MathType is an empty non-polymorphic base class for all * mathematical function types. */ struct MathType {}; -/** The gamma function */ -double gamma(double); - -/** The log of the gamma function */ -double lngamma(double); - -/** Return \f${\rm atanh}(x)\f$ */ -double atanh(double); - /** Return \f$1-e^x\f$, with highest possible precision for * \f$x\rightarrow 0\f$. */ double exp1m(double x); /** Return \f$1\log(1-x)\f$, with highest possible precision for * \f$x\rightarrow 0\f$. */ double log1m(double); /** Return x rased to the integer power p, using recursion. */ double powi(double x, int p); /** Return the integral of \f$x^p dx\f$ between xl and xu. */ inline double pIntegrate(double p, double xl, double xu) { return p == -1.0? log(xu/xl): (pow(xu, p + 1.0) - pow(xl, p + 1.0))/(p + 1.0); } /** Return the integral of \f$x^p dx\f$ between xl and xu. */ inline double pIntegrate(int p, double xl, double xu) { return p == -1? log(xu/xl): (powi(xu, p + 1) - powi(xl, p + 1))/double(p + 1); } /** Return the integral of \f$x^{e-1} dx\f$ between xl and xl+dx with * highest possible precision for \f$dx\rightarrow 0\f$ and/or * \f$e\rightarrow 0\f$. */ inline double pXIntegrate(double e, double xl, double dx) { return e == 0.0? log1m(-dx/xl): -pow(xl, e)*exp1m(e*log1m(-dx/xl))/e; } /** Generate an x between xl and xu distributed as \f$x^p\f$. */ inline double pGenerate(double p, double xl, double xu, double rnd) { return p == -1.0? xl*pow(xu/xl, rnd): pow((1.0 - rnd)*pow(xl, p + 1.0) + rnd*pow(xu, p + 1.0), 1.0/(1.0 + p)); } /** Generate an x between xl and xu distributed as \f$x^p\f$. */ inline double pGenerate(int p, double xl, double xu, double rnd) { return p == -1? xl*pow(xu/xl, rnd): pow((1.0 - rnd)*powi(xl, p + 1) + rnd*powi(xu, p + 1), 1.0/double(1 + p)); } /** Generate an x between xl and xl + dx distributed as \f$x^{e-1}\f$ * with highest possible precision for\f$dx\rightarrow 0\f$ and/or * * \f$e\rightarrow 0\f$. * @param e the parameter defining the power in \f$x^{e-1}\f$. * @param xl the lower bound of the generation interval. * @param dx the interval. * @param rnd a flat random number in the interval ]0,1[. */ inline double pXGenerate(double e, double xl, double dx, double rnd) { return e == 0.0? -xl*exp1m(rnd*log1m(-dx/xl)): -exp1m(log1m(rnd*exp1m(e*log1m(-dx/xl)))/e)*xl; } /** Returns (x - y)/(|x| + |y|). */ template inline double relativeError(FloatType x, FloatType y) { return ( x == y ? 0.0 : double((x - y)/(abs(x) + abs(y))) ); } /** Return x if |x|<|y|, else return y. */ template inline T absmin(const T & x, const T & y) { return abs(x) < abs(y)? x: y; } /** Return x if |x|>|y|, else return y. */ template inline T absmax(const T & x, const T & y) { return abs(x) > abs(y)? x: y; } /** Transfer the sign of the second argument to the first. * @return \f$|x|\f$ if \f$y>0\f$ otherwise return \f$-|x|\f$. */ template inline T sign(T x, U y) { return y > U()? abs(x): -abs(x); } /** Templated class for calculating integer powers. */ //@{ /** * Struct for powers */ template struct Power: public MathType {}; /** * Struct for powers */ template struct Power { /** Member for the power*/ static double pow(double x) { return x*Power::pow(x); } }; /** * Struct for powers */ template struct Power { /** Member for the power*/ static double pow(double x) { return Power::pow(x)/x; } }; /** * Struct for powers */ template <> struct Power<0,true> { /** Member for the power*/ static double pow(double) { return 1.0; } }; /** * Struct for powers */ template <> struct Power<0,false> { /** Member for the power*/ static double pow(double) { return 1.0; } }; //@} /** Templated function to calculate integer powers known at * compile-time. */ template inline double Pow(double x) { return Power::pow(x); } /** This namespace introduces some useful function classes with known * primitive and inverse primitive functions. Useful to sample * corresponding distributions.*/ namespace Functions { /** Class corresponding to functions of the form \f$x^N\f$ with integer N. */ template struct PowX: public MathType { /** The primitive function. */ static double primitive(double x) { return Pow(x)/double(N+1); } /** Integrate function in a given interval. */ static double integrate(double x0, double x1) { return primitive(x1) - primitive(x0); } /** Sample a distribution in a given interval given a flat random * number R in the interval ]0,1[. */ static double generate(double x0, double x1, double R) { return pow(primitive(x0) + R*integrate(x0, x1), 1.0/double(N+1)); } }; /** @cond TRAITSPECIALIZATIONS */ /** * Template for generating according to a specific power */ template <> inline double PowX<1>::generate(double x0, double x1, double R) { return std::sqrt(x0*x0 + R*(x1*x1 - x0*x0)); } /** * Template for generating according to a specific power */ template <> inline double PowX<0>::generate(double x0, double x1, double R) { return x0 + R*(x1 - x0); } /** * Template for generating according to a specific power */ template<> inline double PowX<-1>::primitive(double x) { return log(x); } /** * Template for generating according to a specific power */ template <> inline double PowX<-1>::integrate(double x0, double x1) { return log(x1/x0); } /** * Template for generating according to a specific power */ template <> inline double PowX<-1>::generate(double x0, double x1, double R) { return x0*pow(x1/x0, R); } /** * Template for generating according to a specific power */ template <> inline double PowX<-2>::generate(double x0, double x1, double R) { return x0*x1/(x1 - R*(x1 - x0)); } /** * Template for generating according to a specific power */ template <> inline double PowX<-3>::generate(double x0, double x1, double R) { return x0*x1/std::sqrt(x1*x1 - R*(x1*x1 - x0*x0)); } /** @endcond */ /** Class corresponding to functions of the form \f$(1-x)^N\f$ * with integer N. */ template struct Pow1mX: public MathType { /** The primitive function. */ static double primitive(double x) { return -PowX::primitive(1.0 - x); } /** Integrate function in a given interval. */ static double integrate(double x0, double x1) { return PowX::integrate(1.0 - x1, 1.0 - x0); } /** Sample a distribution in a given interval given a flat random * number R in the interval ]0,1[. */ static double generate(double x0, double x1, double R) { return 1.0 - PowX::generate(1.0 - x1, 1.0 - x0, R); } }; /** Class corresponding to functions of the form \f$1/(x(1-x))\f$ */ struct InvX1mX: public MathType { /** The primitive function. */ static double primitive(double x) { return log(x/(1.0 - x)); } /** Integrate function in a given interval. */ static double integrate(double x0, double x1) { return log(x1*(1.0 - x0)/(x0*(1.0 - x1))); } /** Sample a distribution in a given interval given a flat random * number R in the interval ]0,1[. */ static double generate(double x0, double x1, double R) { double r = pow(x1*(1.0 - x0)/(x0*(1.0 - x1)), R)*x0/(1.0 - x0); return r/(1.0 + r); } }; /** Class corresponding to functions of the form \f$e^x\f$ */ struct ExpX: public MathType { /** The primitive function. */ static double primitive(double x) { return exp(x); } /** Integrate function in a given interval. */ static double integrate(double x0, double x1) { return exp(x1) - exp(x0); } /** Sample a distribution in a given interval given a flat random * number R in the interval ]0,1[. */ static double generate(double x0, double x1, double R) { return log(exp(x0) + R*(exp(x1) - exp(x0))); } }; /** Class corresponding to functions of the form \f$x^{N/D}\f$ * with integer N and D. */ template struct FracPowX: public MathType { /** The primitive function. */ static double primitive(double x) { double r = double(N)/double(D) + 1.0; return pow(x, r)/r; } /** Integrate function in a given interval. */ static double integrate(double x0, double x1) { return primitive(x1) - primitive(x0); } /** Sample a distribution in a given interval given a flat random * number R in the interval ]0,1[. */ static double generate(double x0, double x1, double R) { double r = double(N)/double(D) + 1.0; return pow(primitive(x0) + R*integrate(x0, x1), 1.0/r); } }; } } } #endif /* ThePEG_Math_H */