RandomFuns.C

Go to the documentation of this file.

#include "jbltools/sfh/RandomFuns.h"
#include "jbltools/sfh/IntFun.h"
#include "jbltools/sfh/FloatFun.h"

#include "TH1.h"

#include <cmath>
#include <cassert>
 

static const char *ident="@(#)$Id: RandomFuns.C,v 1.5 2005/11/29 20:03:31 blist Exp $";

RandomFun::RandomFun (const IntFun& seed_, int index_, const char *name_ ) 
   : FloatFun (name_ ? name_ : str("Random(")+(seed_.getName() ? seed_.getName() : "?") + ", ",
     str(index_)+")"), 
     seed(seed_),
     index(index_),
     cleanup (false)
   {}

RandomFun::RandomFun (const IntFun& seed_, int index_, const std::string &name_) 
    : FloatFun (name_), 
      seed(seed_),
      index(index_),
     cleanup (false)
      {}

// Large primes taken from http://primes.utm.edu/lists/small/millions/primes1.zip:
// List of 1st million primes

RandomFun::RandomFun (const IntFun& seed1_, const IntFun& seed2_, int index_, const char *name_ ) 
   : FloatFun (name_ ? name_ : str("Random(")+(seed1_.getName() ? seed1_.getName() : "?") + ", " +
     (seed2_.getName() ? seed2_.getName() : "?") + ", ",str(index_)+")"), 
     seed(1327*seed1_+1997*seed2_),
     index(index_),
     cleanup (true)
   {}

RandomFun::RandomFun (const IntFun& seed1_, const IntFun& seed2_, int index_, const std::string &name_) 
    : FloatFun (name_), 
      seed(1327*seed1_+1997*seed2_),
      index(index_),
      cleanup (true)
      {}
      
const long long RandomFun::a = 16807LL;
const long long RandomFun::c = 0LL;
const long long RandomFun::m = 2147483647LL;

Float_FF RandomFun::operator() () const {
  long long i = 3498493*seed(); 
  i ^= 123456789LL;
  i = std::abs(i);
  i = iterate (i);
  i = ((i & 65535)<<16) | ((i >> 16) & 65535);
  i = iterate (i);
  for (int j = 0; j < 3*index; j++) i = iterate (i);
  return toDouble (i);
//   long long i = 3498493*seed()+7369133*index; 
//   i ^= 123456789LL;
//   i = std::abs(i);
//   i = iterate (i);
//   i = ((i & 65535)<<16) | ((i >> 16) & 65535);
//   i = iterate (i);
//   i = ((i & 65535)<<16) | ((i >> 16) & 65535);
//   i = iterate (i);
//   i = iterate (i);
//   return toDouble (i);
}

RandomFun::~RandomFun() {
  if (cleanup) const_cast<IntFun&>(seed).destroy();
}

const FillIterator *RandomFun::getIterator() const {
  return seed.getIterator();
}

GaussRandomFun::GaussRandomFun (const FloatFun& rnd_, double mean_, double sigma_, const char *name_ ) 
   : FloatFun (name_ ? 
               str(name_) : 
               str("GaussRandom(")+(rnd_.getName() ? rnd_.getName() : "?") + ", " + str(mean_) + ", " + str(sigma_)+")"), 
     rnd(rnd_),
     mean(mean_),
     sigma (sigma_)
   {}

GaussRandomFun::GaussRandomFun (const FloatFun& rnd_, double mean_, double sigma_, const std::string &name_) 
    : FloatFun (name_), 
      rnd(rnd_),
      mean(mean_),
      sigma (sigma_)
      {}

const FillIterator *GaussRandomFun::getIterator() const {
  return rnd.getIterator();
}

      
Float_FF GaussRandomFun::operator() () const {
  double x1 = rnd();
  double x2 = RandomFun::generateNext(x1);
  double sq = sqrt(-2*log(x2));
  return sigma*sin(2*3.141592653589793238*x1)*sq + mean;
}


Poisson1Fun::Poisson1Fun (const FloatFun& rnd_, double mean_, const char *name_) 
  : IntFun (name_ ? name_ : str("Poisson(")+rnd_.getName()+")"),
    rnd(rnd_),
    mean (mean_)
{}

int Poisson1Fun::operator() () const {
  double r = rnd();
  double e = exp (-mean);
  double x = 0;
  double s = exp (-mean);
  for (int result = 0; result < 126; s *= mean/(++result)) {
    x += s;
    if (r < x) return result;
  }
  return 126;
}

const FillIterator *Poisson1Fun::getIterator() const {
  return rnd.getIterator();
}

Poisson2Fun::Poisson2Fun (const FloatFun& rnd_, double mean_, const char *name_) 
  : IntFun (name_ ? name_ : str("Poisson(")+rnd_.getName()+")"),
    gaussrnd(*new GaussRandomFun (rnd_, mean_, std::sqrt(std::abs(mean_))))
{}

Poisson2Fun::Poisson2Fun (const FloatFun& rnd_, double mean_, const std::string& name_) 
  : IntFun (name_),
    gaussrnd(*new GaussRandomFun (rnd_, mean_, std::sqrt(std::abs(mean_))))
{}

int Poisson2Fun::operator() () const {
  double g = gaussrnd();
  return (g < 0) ? 0 : (int)(g+0.5);
}

const FillIterator *Poisson2Fun::getIterator() const {
  return gaussrnd.getIterator();
}

Poisson2Fun::~Poisson2Fun() {
  gaussrnd.destroy();
}

HistoRandomFun::HistoRandomFun (const FloatFun& rnd_, const TH1& histo, const char *name_ ) 
  : FloatFun (name_ ? str(name_) : str("HistoRandomFun(") + rnd_.getName() + ")"),
    rnd (rnd_),
    nbins (0),
    edges (0),
    values (0)
{
  init (histo);
}
HistoRandomFun::HistoRandomFun (const FloatFun& rnd_, const TH1& histo, const std::string& name_ ) 
  : FloatFun (name_),
    rnd (rnd_),
    nbins (0),
    edges (0),
    values (0)
{
  init (histo);
}

Float_FF HistoRandomFun::operator() () const {
  if (nbins <= 0) return 0;
  assert (edges);
  assert (values);
  double x1 = rnd();
  double x2 = RandomFun::generateNext(x1);
  int ibin = 0;
  for (ibin=0; (ibin < nbins-1) && values[ibin] <= x1; ibin++);
  assert (ibin >=0);
  assert (ibin < nbins);
  return edges[ibin] + x2*(edges[ibin+1]-edges[ibin]);
}

const FillIterator *HistoRandomFun::getIterator() const {
  return rnd.getIterator();
}

void HistoRandomFun::init (const TH1& histo) {
  nbins = histo.GetNbinsX()+2;
  edges = new double[nbins+1];
  values = new double[nbins];
  values[0] = histo.GetBinContent(0); // underflow bin
  for (int ibin = 1; ibin < nbins; ibin++) {
    edges[ibin] = histo.GetBinLowEdge(ibin);
    values[ibin] = values[ibin-1] + histo.GetBinContent(ibin);
  }
  edges[0] = nbins>2 ? 2*edges[1]-edges[2] : edges[1]-1; // make underflow bin as wide as first bin
  edges[nbins] = 2*edges[nbins-1]-edges[nbins-2];        // make overflow bin as wide as last bin
  double f = values[nbins-1] ? 1./values[nbins-1] : 1;
  for (int ibin = 0; ibin < nbins; ibin++) values[ibin] *= f;
}

HistoRandomFun::~HistoRandomFun() {
  delete[] edges;
  delete[] values;
}


HistoRandomIntFun::HistoRandomIntFun (const FloatFun& rnd_, const TH1& histo, const char *name_ ) 
  : IntFun (name_ ? str(name_) : str("HistoRandomIntFun(") + rnd_.getName() + ")"),
    rnd (rnd_),
    nbins (0),
    edges (0),
    values (0)
{
  init (histo);
}
HistoRandomIntFun::HistoRandomIntFun (const FloatFun& rnd_, const TH1& histo, const std::string& name_ ) 
  : IntFun (name_),
    rnd (rnd_),
    nbins (0),
    edges (0),
    values (0)
{
  init (histo);
}

int HistoRandomIntFun::operator() () const {
  if (nbins <= 0) return 0;
  assert (edges);
  assert (values);
  double x1 = rnd();
  double x2 = RandomFun::generateNext(x1);
  int ibin = 0;
  for (ibin=0; (ibin < nbins+1) && values[ibin] <= x1; ibin++);
  assert (ibin >=0);
  assert (ibin < nbins);
  return (int)(std::floor(edges[ibin] + x2*(edges[ibin+1]-edges[ibin]) + 0.5));
}

const FillIterator *HistoRandomIntFun::getIterator() const {
  return rnd.getIterator();
}

void HistoRandomIntFun::init (const TH1& histo) {
  nbins = histo.GetNbinsX()+2;
  edges = new double[nbins+1];
  values = new double[nbins];
  values[0] = histo.GetBinContent(0); // underflow bin
  for (int ibin = 1; ibin < nbins; ibin++) {
    edges[ibin] = histo.GetBinLowEdge(ibin);
    values[ibin] = values[ibin-1] + histo.GetBinContent(ibin);
  }
  edges[0] = nbins>2 ? 2*edges[1]-edges[2] : edges[1]-1;
  edges[nbins] = 2*edges[nbins-1]-edges[nbins-2];
  double f = values[nbins-1] ? 1./values[nbins-1] : 1;
  for (int ibin = 0; ibin < nbins; ibin++) values[ibin] *= f;
  // Round edges
  for (int ibin = 0; ibin <= nbins; ibin++) {
    if (edges[ibin] == std::floor(edges[ibin])) {
      edges[ibin] = edges[ibin]-0.5;
    }
    else { 
      edges[ibin] = std::floor(edges[ibin]) + 0.5;
    }
  }
}

HistoRandomIntFun::~HistoRandomIntFun() {
  delete[] edges;
  delete[] values;
}

FloatFun& gauss (const FloatFunPoR& rnd, double mean, double sigma) {
  if (!rnd.pff) return *new ConstFun(0);
  return *new GaussRandomFun (*(rnd.pff), mean, sigma);
}


IntFun& poisson (const FloatFunPoR& rnd,       
                 double mean                   
                ) {
  if (!rnd.pff) return *new ConstIntFun(0);
  return (mean < 81) ? 
          static_cast<IntFun&>(*new Poisson1Fun (*(rnd.pff), mean)) : 
          static_cast<IntFun&>(*new Poisson2Fun (*(rnd.pff), mean));
}

---