Sampler.cxx
Go to the documentation of this file.
1 // GENIE includes
2 #include "GENIE/Framework/Messenger/Messenger.h"
3 #include "GENIE/Framework/Numerical/RandomGen.h"
4 
5 // NOvASoft includes
7 
8 using namespace genie::supernova;
9 
10 
11 //.......................................................................
13 
14 
15 //.......................................................................
17 
18 
19 //.......................................................................
21 {
22  fHist = (TH1*)hist->Clone();
23  fHist->SetDirectory(0);
24 
25  LOG("Sampler", pDEBUG) << "init from hist \"" << fHist->GetName()
26  << "\", int=" << fHist->Integral();
27 }
28 
29 
30 //.......................................................................
32 {
33  LOG("Sampler", pDEBUG) << "get random from hist=" << fHist
34  << "(name=" << fHist->GetName()
35  << "), int=" << fHist->Integral();
36 
37  return fHist->GetRandom();
38 }
39 
40 
41 
42 
43 //.......................................................................
44 OrderedSampler::OrderedSampler(): fProbScale(1) {}
45 
46 
47 //.......................................................................
49 
50 
51 //.......................................................................
53 {
54  LOG("Sampler", pDEBUG) << "init";
55 
56  fBin = 0;
57  int Nbins = hist->GetNbinsX() + 1;
58 
59  fBinEdges.resize(Nbins);
60  fBinValues.resize(Nbins);
61 
62  for (int i = 0; i < Nbins; ++i) {
63  fBinEdges[i] = hist->GetBinLowEdge(i);
64  fBinValues[i] = hist->GetBinContent(i);
65  if(integral)
66  fBinValues[i] /= hist->GetBinWidth(i);
67  }
68 
69  // Start from bin#1, skipping underflow
70  this->StartNextBin();
71  LOG("Sampler", pNOTICE) << "Ordered sample with nbins=" << Nbins
72  << " [" << fBinEdges.front() << "," << fBinEdges.back() << "]";
73 }
74 
75 
76 //.......................................................................
78 {
79  return fBin >= fBinEdges.size() - 1;
80 }
81 
82 
83 //.......................................................................
85 {
86  LOG("Sampler", pDEBUG) << "from pos=" << fPos;
87 
88  if (this->End()) return 0;
89  do {
90  if (fPos >= this->BinEnd()) this->StartNextBin();
91  if (End()) return 0;
92 
93  // Generate random step, using the average distance in current bin
94  this->Step();
95  LOG("Sampler", pDEBUG) << "now pos=" << fPos
96  << " (bin =[" << this->BinStart()
97  << "," << this->BinEnd() << "])";
98  } while (fPos >= this->BinEnd());
99 
100  LOG("Sampler", pDEBUG) << "return pos=" << fPos;
101  return fPos;
102 }
103 
104 
105 // ......................................................................
107 {
108  return fBin;
109 }
110 
111 
112 // ......................................................................
114 {
115  fBin = 0;
116  this->StartNextBin();
117 }
118 
119 
120 // ......................................................................
122 {
123  fProbScale = scale;
124 }
125 
126 
127 // ......................................................................
128 const size_t OrderedSampler::Nbins()
129 {
130  return fBinValues.size();
131 }
132 
133 
134 // ......................................................................
136 {
137  return fBinEdges.data();
138 }
139 
140 
141 // ......................................................................
142 const double OrderedSampler::GetBinValue(size_t n) {
143  return fBinValues[n] * fProbScale;
144 }
145 
146 
147 // ......................................................................
149 {
150  return fBinEdges[fBin];
151 }
152 
153 
154 // ......................................................................
156 {
157  return fBinEdges[fBin + 1];
158 }
159 
160 
161 //.......................................................................
163 {
164  double binValue = this->GetBinValue(fBin);
165 
166  // If the bin value is too small, move on to the next bin.
167  if (binValue < 1e-15) {
168  fPos = this->BinEnd();
169  return;
170  }
171 
173 
174  float tau = 1. / binValue;
175  double step = rnd->RndFlux().Exp(tau);
176 
177  LOG("Sampler", pDEBUG) << "step by " << step << ": " << fPos << "->" << fPos + step;
178 
179  fPos += step;
180 
181  return;
182 }
183 
184 
185 //.......................................................................
187  {
188  // Go to next bin and update all the variables
189  ++fBin;
190  LOG("Sampler", pDEBUG) << "Starting bin#" << fBin << "/" << fBinValues.size()
191  << "(" << this->End() << ")";
192 
193  // Start the steps from beginning of this bin
194  fPos = this->BinStart();
195 
196  if (this->End()) return;
197 
198  LOG("Sampler", pDEBUG) << "Starting bin#" << fBin
199  << " [" << this->BinStart() << "," << this->BinEnd() << "]"
200  << " Freq=" << fBinValues[fBin];
201 }
void SetProbScale(double scale)
Definition: Sampler.cxx:121
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const double * GetBinEdges()
Definition: Sampler.cxx:135
const double GetBinValue(size_t n)
Definition: Sampler.cxx:142
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
Double_t scale
Definition: plot.C:25
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
std::vector< double > fBinValues
Definition: Sampler.h:49
std::vector< double > fBinEdges
Definition: Sampler.h:48
#define pNOTICE
Definition: Messenger.h:62
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
Float_t e
Definition: plot.C:35
void Init(TH1 *hist, bool integral=false)
Definition: Sampler.cxx:52
#define pDEBUG
Definition: Messenger.h:64