Page MenuHomeHEPForge
Files F7878538

No OneTemporary
Actions

View Options

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 <cmath>
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<double>::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> LeptonLeptonPDF::initLeptonLeptonPDF;
void LeptonLeptonPDF::Init() {
static ClassDocumentation<LeptonLeptonPDF> 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 <config.h>
#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 <cmath>
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 <cmath>
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 <typename FloatType>
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 <typename T>
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 <typename T>
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 <typename T, typename U>
inline T sign(T x, U y) {
return y > U()? abs(x): -abs(x);
}
/** Templated class for calculating integer powers. */
//@{
/**
* Struct for powers
*/
template <int N, bool Inv>
struct Power: public MathType {};
/**
* Struct for powers
*/
template <int N>
struct Power<N,false> {
/** Member for the power*/
static double pow(double x) { return x*Power<N-1,false>::pow(x); }
};
/**
* Struct for powers
*/
template <int N>
struct Power<N,true> {
/** Member for the power*/
static double pow(double x) { return Power<N+1,true>::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 <int N>
inline double Pow(double x) { return Power<N, (N < 0)>::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 <int N>
struct PowX: public MathType {
/** The primitive function. */
static double primitive(double x) {
return Pow<N+1>(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 <int N>
struct Pow1mX: public MathType {
/** The primitive function. */
static double primitive(double x) {
return -PowX<N>::primitive(1.0 - x);
}
/** Integrate function in a given interval. */
static double integrate(double x0, double x1) {
return PowX<N>::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<N>::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 <int N, int D>
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 */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:13 PM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805542
Default Alt Text
(15 KB)

Event Timeline