Commit de30ff96 authored by Enrico Bothmann's avatar Enrico Bothmann Committed by Frank Krauss

Fix (limited) Gaussian random number generation

... for the primordial kT determination. The Marsaglia algorithm was
wrongly implemented (and missed an opportunity for optimisation by not
storing the second Gaussian random number for later use), and the
LimitedWeight function did not rescale its hit-or-miss comparison
function to a maximum of 1, which was extremely inefficient for low
kTmax, where basically every kT candidate was unnecessarily rejected by
LimitedWeight.

Other changes
===

The Gaussian random number algorithm is now implemented in Random.H and
supersedes the previous one, which was less efficient. At that
opportunity, I have also removed superfluous `inline` keyword uses in
front of member function definitions.

Observations in performance measurements
===

- When running 10k events for Examples/CI/LO_Z, this will increase event
  generation speed by ~50% when using -O3. This measurement will vary
  strongly, because events with very low kTmax seem to occur only every
  few thousand events.

- There is however a significant increase in Retried events from
  Hadron_Decays (before 9, after 38).

Further remarks
===

- Note that there is a chance that this fix renders previous tuning partly
  invalid.

- Note that the Primordial kT generation using a Gaussian with a fixed
  width and an a-posteriori hit-or-miss Monte-Carlo that can have a
  width that is way smaller might still be very inefficient. The Gaussian
  generates a lot of large numbers before being accepted by the
  LimitedWeight. The Gaussian takes a kTmax argument but does not use it.
  This raises the question if perhaps this is not intended behaviour in
  the first place.
parent d9efa01c
......@@ -54,9 +54,8 @@ double Matter_Overlap::SelectB() const {
break;
}
do {
ran->Gaussian(b,b2);
b = dabs(b)*radius;
if (b>m_bmax) b = dabs(b2)*radius;
auto b = dabs(ran->GetGaussian())*radius;
if (b>m_bmax) b = dabs(ran->GetGaussian())*radius;
} while (b>m_bmax);
return b;
}
......
......@@ -418,16 +418,6 @@ ptrdiff_t ATOOLS::Random::operator() (ptrdiff_t max) {
return Min(static_cast<ptrdiff_t>(Get() * max),max-1);
}
void ATOOLS::Random::Gaussian(double & x,double & y)
{
double phi(2.*M_PI*Get()), random(Get());
while (random==0.) random = Get();
double r(sqrt(-2.*log(random)));
x = r*std::cos(phi);
y = r*std::sin(phi);
}
External_RNG::~External_RNG()
{
}
......@@ -69,7 +69,7 @@ namespace ATOOLS {
void SetSeed(long nid); // seed for Rnd2()
void SetSeed(unsigned int i1,unsigned int i2,
unsigned int i3,unsigned int i4);
inline long int GetSeed() { return m_id; }
long int GetSeed() { return m_id; }
int WriteOutStatus(const char *outfile);
int WriteOutStatus(std::ostream &os,const size_t &idx);
......@@ -80,17 +80,35 @@ namespace ATOOLS {
void RestoreStatus();
void FastForward(const size_t &n);
void ResetToLastIncrementedSeed();
inline void EraseLastIncrementedSeed()
void EraseLastIncrementedSeed()
{ m_lastincrementedseed.str(std::string()); }
inline void SetSeedStorageIncrement(size_t inc) { m_increment=inc; }
void SetSeedStorageIncrement(size_t inc) { m_increment=inc; }
// return uniformly distributed random number in [0,1] using active Generator
double Get();
ptrdiff_t operator() (ptrdiff_t max);
// produce two Gaussian distributed random number using active Generator
void Gaussian(double &,double &);
// produce Gaussian distributed random number using active Generator
// according to the Marsaglia method
double GetGaussian() {
static auto hasSpare = false;
static double spare;
if (hasSpare) {
hasSpare = false;
return spare;
}
double ran1, ran2, R;
do {
ran1 = 2.*Get()-1.;
ran2 = 2.*Get()-1.;
R = ran1*ran1+ran2*ran2;
} while (R>1. || R==0.);
R = sqrt(-2.*log(R)/R);
hasSpare = true;
spare = ran2 * R;
return ran1 * R;
}
// produce Poissonian distributed random number using active Generator
inline double Poissonian(const double & lambda) {
double Poissonian(const double & lambda) {
if(lambda>500.) {
double u = Get();
double v = Get();
......@@ -102,10 +120,10 @@ namespace ATOOLS {
return N;
}
inline double Theta() { return acos(2.*Get()-1.); }
double GetNZ();
double Theta() { return acos(2.*Get()-1.); }
double GetNZ();
inline External_RNG* GetExternalRng() { return p_external; }
External_RNG* GetExternalRng() { return p_external; }
};// end of class Random
......@@ -156,7 +174,7 @@ namespace ATOOLS {
*/
/*!
\fn inline double Random::Get()
\fn double Random::Get()
\brief is the main routine, returns a single random number in [0,1]
The number is determined either by using Ran2() or Ran4(), depending on
......@@ -164,12 +182,12 @@ namespace ATOOLS {
*/
/*!
\fn inline double Random::GetNZ()
\fn double Random::GetNZ()
\brief retrun a not zero random number
*/
/*!
\fn inline long Random::GetSeed()
\fn long Random::GetSeed()
\brief returns a the seed
No corresponding method for Ran4() exists so far.
......@@ -186,7 +204,7 @@ namespace ATOOLS {
*/
/*!
\fn inline double Random::Theta()
\fn double Random::Theta()
\brief returns an angle \f$\phi\f$ for a uniform \f$cos(\phi)\f$ distribution
*/
......
......@@ -187,10 +187,9 @@ void IsotropicSpectator::GeneratePoint(ATOOLS::Vec4D * p,
double lambda_qcd=0.2, pspat, critE2(sqr(m_residual_mass));
do {
do {
double gauss1, gauss2;
ran->Gaussian(gauss1,gauss2);
pspat = lambda_qcd+lambda_qcd/m_decayer_mass*gauss1;
if(pspat<1e-6) pspat = lambda_qcd+lambda_qcd/m_decayer_mass*gauss2;
pspat = lambda_qcd+lambda_qcd/m_decayer_mass*ran->GetGaussian();
if (pspat < 1e-6)
pspat = lambda_qcd+lambda_qcd/m_decayer_mass*ran->GetGaussian();
} while(pspat<1e-6);
} while (sqr(p[0][0]-sqrt(sqr(m_spectator_mass)+sqr(pspat))) - sqr(pspat) < critE2);
......
......@@ -119,22 +119,13 @@ Vec4D Primordial_KPerp::KT(const Vec4D & mom) {
}
double Primordial_KPerp::KT_Gauss(const double & ktmax) const {
// Generate normalised Gaussian random numbers according to the
// Marsaglia method
double ran1, ran2, R;
do {
ran1 = 2.*ran->Get()-1.;
ran2 = 2.*ran->Get()-1.;
R = ran1*ran1+ran2*ran2;
} while (R>1. || R==0.);
R = sqrt(-2.*log(R)/R);
// shift the Gaussian random numbers
return m_mean[m_beam] + R * m_sigma[m_beam];
return m_mean[m_beam] + dabs(ran->GetGaussian()) * m_sigma[m_beam];
}
double Primordial_KPerp::KT_Gauss_Limited(const double & ktmax) const {
// Generate normalised Gaussian random numbers according to the
// Marsaglia method, with an additional polynomial limitation
// Generate normalised Gaussian random numbers
// with an additional polynomial limitation
double ran1, ran2, R, kt;
do {
kt = KT_Gauss(ktmax);
......@@ -165,7 +156,10 @@ double Primordial_KPerp::DipoleWeight(const double & kt) const {
}
double Primordial_KPerp::LimitedWeight(const double & kt) const {
return Max(0.,pow(m_ktmax[m_beam],m_eta[m_beam])-pow(kt,m_eta[m_beam]));
if (kt < m_ktmax[m_beam])
return 1.0 - pow(kt/m_ktmax[m_beam], m_eta[m_beam]);
else
return 0.0;
}
void Primordial_KPerp::InitAnalysis() {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment