50#include "CLHEP/Random/defs.h"
51#include "CLHEP/Random/RandGeneral.h"
52#include "CLHEP/Random/DoubConv.hh"
71 InterpolationType(IntType)
73 prepareTable(aProbFunc);
77 const double* aProbFunc,
83 InterpolationType(IntType)
85 prepareTable(aProbFunc);
89 const double* aProbFunc,
93 localEngine(anEngine),
95 InterpolationType(IntType)
97 prepareTable(aProbFunc);
100void RandGeneral::prepareTable(
const double* aProbFunc) {
106 "RandGeneral constructed with no bins - will use flat distribution\n";
107 useFlatDistribution();
111 theIntegralPdf.resize(nBins+1);
112 theIntegralPdf[0] = 0;
114 register double weight;
116 for ( ptn = 0; ptn<nBins; ++ptn ) {
117 weight = aProbFunc[ptn];
122 "RandGeneral constructed with negative-weight bin " << ptn <<
123 " = " << weight <<
" \n -- will substitute 0 weight \n";
127 theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
130 if ( theIntegralPdf[nBins] <= 0 ) {
132 "RandGeneral constructed nothing in bins - will use flat distribution\n";
133 useFlatDistribution();
137 for ( ptn = 0; ptn < nBins+1; ++ptn ) {
138 theIntegralPdf[ptn] /= theIntegralPdf[nBins];
143 oneOverNbins = 1.0 / nBins;
147 if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
149 "RandGeneral does not recognize IntType " << InterpolationType
150 <<
"\n Will use type 0 (continuous linear interpolation \n";
151 InterpolationType = 0;
156void RandGeneral::useFlatDistribution() {
161 theIntegralPdf.resize(2);
162 theIntegralPdf[0] = 0;
163 theIntegralPdf[1] = 1;
181double RandGeneral::mapRandom(
double rand)
const {
191 while (nabove > nbelow+1) {
192 middle = (nabove + nbelow+1)>>1;
193 if (rand >= theIntegralPdf[middle]) {
199 assert ( nabove == nbelow+1 );
200 assert ( theIntegralPdf[nbelow] <= rand );
201 assert ( theIntegralPdf[nabove] >= rand );
205 if ( InterpolationType == 1 ) {
207 return nbelow * oneOverNbins;
211 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
215 if ( binMeasure == 0 ) {
219 return (nbelow + .5) * oneOverNbins;
222 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
224 return (nbelow + binFraction) * oneOverNbins;
232 const int size,
double* vect )
236 for (i=0; i<size; ++i) {
237 vect[i] =
shoot(anEngine);
245 for (i=0; i<size; ++i) {
251 int pr=os.precision(20);
252 std::vector<unsigned long> t(2);
253 os <<
" " <<
name() <<
"\n";
254 os <<
"Uvec" <<
"\n";
255 os << nBins <<
" " << oneOverNbins <<
" " << InterpolationType <<
"\n";
257 os << t[0] <<
" " << t[1] <<
"\n";
258 assert (
static_cast<int>(theIntegralPdf.size())==nBins+1);
259 for (
unsigned int i=0; i<theIntegralPdf.size(); ++i) {
261 os << theIntegralPdf[i] <<
" " << t[0] <<
" " << t[1] <<
"\n";
266 int pr=os.precision(20);
267 os <<
" " <<
name() <<
"\n";
268 os << nBins <<
" " << oneOverNbins <<
" " << InterpolationType <<
"\n";
269 assert (
static_cast<int>(theIntegralPdf.size())==nBins+1);
270 for (
unsigned int i=0; i<theIntegralPdf.size(); ++i)
271 os << theIntegralPdf[i] <<
"\n";
280 if (inName !=
name()) {
281 is.clear(std::ios::badbit | is.rdstate());
282 std::cerr <<
"Mismatch when expecting to read state of a "
283 <<
name() <<
" distribution\n"
284 <<
"Name found was " << inName
285 <<
"\nistream is left in the badbit state\n";
289 std::vector<unsigned long> t(2);
290 is >> nBins >> oneOverNbins >> InterpolationType;
292 theIntegralPdf.resize(nBins+1);
293 for (
unsigned int i=0; i<theIntegralPdf.size(); ++i) {
294 is >> theIntegralPdf[i] >> t[0] >> t[1];
300 is >> oneOverNbins >> InterpolationType;
301 theIntegralPdf.resize(nBins+1);
302 for (
unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i];
static double longs2double(const std::vector< unsigned long > &v)
static std::vector< unsigned long > dto2longs(double d)
RandGeneral(const double *aProbFunc, int theProbSize, int IntType=0)
std::istream & get(std::istream &is)
std::ostream & put(std::ostream &os) const
void fireArray(const int size, double *vect)
HepRandomEngine & engine()
void shootArray(const int size, double *vect)
bool possibleKeywordInput(IS &is, const std::string &key, T &t)