Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DRandom2.h
Go to the documentation of this file.
1 // $Id$
2 //
3 
4 // Random number generator used in mcsmear. All random numbers
5 // should come from the global "gDRandom" object declared here.
6 //
7 // Because we want to record the seeds used for every event,
8 // we use the TRandom2 class. This one has only 3 seed values
9 // (as opposed to 24 for TRandom1 and 624 for TRandom3) but
10 // is fast with a large period (10^26).
11 //
12 // Because the seeds in TRandom2 are declared protected in
13 // the class with no methods to access/set them, we derive
14 // a new class, DRandom2 from TRandom2. This allows us access
15 // to the numbers for easy recording/retrieving.
16 
17 #ifndef _DRANDOM2_H_
18 #define _DRANDOM2_H_
19 
20 #include <TRandom2.h>
21 #include <iostream>
22 using std::cerr;
23 using std::endl;
24 
25 class DRandom2:public TRandom2{
26  public:
27 
28  DRandom2(UInt_t seed=1):TRandom2(seed){}
29 
30  void GetSeeds(UInt_t &seed, UInt_t &seed1, UInt_t &seed2){
31  seed = this->fSeed;
32  seed1 = this->fSeed1;
33  seed2 = this->fSeed2;
34  }
35 
36  void SetSeeds(UInt_t &seed, UInt_t &seed1, UInt_t &seed2){
37 
38  // See the comments in TRandom2::SetSeed(int)
39  if( (seed<2) | (seed1<8) | (seed2<16) ){
40  cerr << endl;
41  cerr << "*********************************************************" << endl;
42  cerr << "WARNING: Random seeds passed to DRandom2::SetSeeds have" << endl;
43  cerr << "forbidden values: " << endl;
44  cerr << " seed = " << seed << " (must be at least 2)" <<endl;
45  cerr << " seed1 = " << seed1 << " (must be at least 8)" <<endl;
46  cerr << " seed1 = " << seed2 << " (must be at least 16)" <<endl;
47  cerr << "See comments in source for TRandom2::SetSeed(int)" <<endl;
48  cerr << "The seeds will all be adjusted to be in range." <<endl;
49  cerr << "*********************************************************" << endl;
50  cerr << endl;
51  seed += 2;
52  seed1 += 8;
53  seed2 += 15;
54  }
55 
56 
57  this->fSeed = seed;
58  this->fSeed1 = seed1;
59  this->fSeed2 = seed2;
60  }
61 
62  // legacy mcsmear interface
63  inline double SampleGaussian(double sigma) {
64  return Gaus(0.0, sigma);
65  }
66 
67  inline double SamplePoisson(double lambda) {
68  return Poisson(lambda);
69  }
70 
71  inline double SampleRange(double x1, double x2) {
72  double s, f;
73  double xlo, xhi;
74 
75  if(x1<x2){
76  xlo = x1;
77  xhi = x2;
78  }else{
79  xlo = x2;
80  xhi = x1;
81  }
82 
83  s = Rndm();
84  f = xlo + s*(xhi-xlo);
85 
86  return f;
87  }
88 
89  inline bool DecideToAcceptHit(double prob) {
90  // This function is used for sculpting simulation efficiencies to match those
91  // in the data. With the data/sim matching efficiency as an input parameter,
92  // this function decides if we should keep the hit or not
93 
94  // Tolerance for seeing if a number is near zero
95  // Could use std::numeric_limits::epsilon(), but that's probably too restrictive
96  const double maxAbsDiff = 1.e-8;
97 
98  // If the hit efficiency is zero, then always reject it
99  // For floating point numbers, using the absolute difference to see if a
100  // number is consistent with zero is preferred.
101  // Reference: https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
102  if( fabs(prob - 0.0) < maxAbsDiff )
103  return false;
104 
105  // If the effiency is less than 0, then we always reject it
106  // This really shouldn't happen, though
107  if( prob < 0.0 )
108  return false;
109 
110  // If the efficiency is greater than 1, then always accept it
111  // (though why would it be larger?)
112  if( prob > 1.0 )
113  return true;
114 
115  // If the efficiency is equal to one, then always accept it
116  if(AlmostEqual(prob, 1.0))
117  return true;
118 
119  // Otherwise, our efficiency should be some number in (0,1)
120  // Throw a random number in that range, and reject if the random
121  // number is larger than our efficiency
122  if(Uniform() > prob)
123  return false;
124  return true;
125  }
126 
127  private:
128  bool AlmostEqual(double a, double b) {
129  // Comparing floating point numbers is tough!
130  // This should work for numbers not near zero.
131  // Reference: https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
132  const double maxRelDiff = 1.e-8;
133 
134  double diff = fabs(a-b);
135  a = fabs(a);
136  b = fabs(b);
137  double largest = (b > a) ? b : a;
138 
139  if(diff <= largest*maxRelDiff)
140  return true;
141  return false;
142  }
143 
144 };
145 
146 #endif // _DRANDOM2_H_
147 
148 extern DRandom2 gDRandom;
149 
150 
bool AlmostEqual(double a, double b)
Definition: DRandom2.h:128
DRandom2 gDRandom
DRandom2(UInt_t seed=1)
Definition: DRandom2.h:28
TF1 * f
Definition: FitGains.C:21
void SetSeeds(UInt_t &seed, UInt_t &seed1, UInt_t &seed2)
Definition: DRandom2.h:36
double SampleGaussian(double sigma)
Definition: DRandom2.h:63
Double_t sigma[NCHANNELS]
Definition: st_tw_resols.C:37
double SampleRange(double x1, double x2)
Definition: DRandom2.h:71
bool DecideToAcceptHit(double prob)
Definition: DRandom2.h:89
void GetSeeds(UInt_t &seed, UInt_t &seed1, UInt_t &seed2)
Definition: DRandom2.h:30
int seed
double SamplePoisson(double lambda)
Definition: DRandom2.h:67