00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include "CLHEP/Random/RandGeneral.h"
00051 #include "CLHEP/Random/DoubConv.h"
00052 #include <cassert>
00053
00054 namespace CLHEP {
00055
00056 std::string RandGeneral::name() const {return "RandGeneral";}
00057 HepRandomEngine & RandGeneral::engine() {return *localEngine;}
00058
00059
00061
00063
00064 RandGeneral::RandGeneral( const double* aProbFunc,
00065 int theProbSize,
00066 int IntType )
00067 : HepRandom(),
00068 localEngine(HepRandom::getTheEngine(), do_nothing_deleter()),
00069 nBins(theProbSize),
00070 InterpolationType(IntType)
00071 {
00072 prepareTable(aProbFunc);
00073 }
00074
00075 RandGeneral::RandGeneral(HepRandomEngine& anEngine,
00076 const double* aProbFunc,
00077 int theProbSize,
00078 int IntType )
00079 : HepRandom(),
00080 localEngine(&anEngine, do_nothing_deleter()),
00081 nBins(theProbSize),
00082 InterpolationType(IntType)
00083 {
00084 prepareTable(aProbFunc);
00085 }
00086
00087 RandGeneral::RandGeneral(HepRandomEngine* anEngine,
00088 const double* aProbFunc,
00089 int theProbSize,
00090 int IntType )
00091 : HepRandom(),
00092 localEngine(anEngine),
00093 nBins(theProbSize),
00094 InterpolationType(IntType)
00095 {
00096 prepareTable(aProbFunc);
00097 }
00098
00099 void RandGeneral::prepareTable(const double* aProbFunc) {
00100
00101
00102
00103 if (nBins < 1) {
00104 std::cerr <<
00105 "RandGeneral constructed with no bins - will use flat distribution\n";
00106 useFlatDistribution();
00107 return;
00108 }
00109
00110 theIntegralPdf.resize(nBins+1);
00111 theIntegralPdf[0] = 0;
00112 register int ptn;
00113 register double weight;
00114
00115 for ( ptn = 0; ptn<nBins; ++ptn ) {
00116 weight = aProbFunc[ptn];
00117 if ( weight < 0 ) {
00118
00119
00120 std::cerr <<
00121 "RandGeneral constructed with negative-weight bin " << ptn <<
00122 " = " << weight << " \n -- will substitute 0 weight \n";
00123 weight = 0;
00124 }
00125
00126 theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight;
00127 }
00128
00129 if ( theIntegralPdf[nBins] <= 0 ) {
00130 std::cerr <<
00131 "RandGeneral constructed nothing in bins - will use flat distribution\n";
00132 useFlatDistribution();
00133 return;
00134 }
00135
00136 for ( ptn = 0; ptn < nBins+1; ++ptn ) {
00137 theIntegralPdf[ptn] /= theIntegralPdf[nBins];
00138
00139 }
00140
00141
00142 oneOverNbins = 1.0 / nBins;
00143
00144
00145
00146 if ( (InterpolationType != 0) && (InterpolationType != 1) ) {
00147 std::cerr <<
00148 "RandGeneral does not recognize IntType " << InterpolationType
00149 << "\n Will use type 0 (continuous linear interpolation \n";
00150 InterpolationType = 0;
00151 }
00152
00153 }
00154
00155 void RandGeneral::useFlatDistribution() {
00156
00157
00158
00159 nBins = 1;
00160 theIntegralPdf.resize(2);
00161 theIntegralPdf[0] = 0;
00162 theIntegralPdf[1] = 1;
00163 oneOverNbins = 1.0;
00164 return;
00165
00166 }
00167
00169
00171
00172 RandGeneral::~RandGeneral() {
00173 }
00174
00175
00177
00179
00180 double RandGeneral::mapRandom(double rand) const {
00181
00182
00183
00184
00185
00186 int nbelow = 0;
00187 int nabove = nBins;
00188 int middle;
00189
00190 while (nabove > nbelow+1) {
00191 middle = (nabove + nbelow+1)>>1;
00192 if (rand >= theIntegralPdf[middle]) {
00193 nbelow = middle;
00194 } else {
00195 nabove = middle;
00196 }
00197 }
00198 assert ( nabove == nbelow+1 );
00199 assert ( theIntegralPdf[nbelow] <= rand );
00200 assert ( theIntegralPdf[nabove] >= rand );
00201
00202
00203
00204 if ( InterpolationType == 1 ) {
00205
00206 return nbelow * oneOverNbins;
00207
00208 } else {
00209
00210 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow];
00211
00212
00213
00214 if ( binMeasure == 0 ) {
00215
00216
00217
00218 return (nbelow + .5) * oneOverNbins;
00219 }
00220
00221 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure;
00222
00223 return (nbelow + binFraction) * oneOverNbins;
00224 }
00225
00226 }
00227
00228 void RandGeneral::shootArray( HepRandomEngine* anEngine,
00229 const int size, double* vect )
00230 {
00231 register int i;
00232
00233 for (i=0; i<size; ++i) {
00234 vect[i] = shoot(anEngine);
00235 }
00236 }
00237
00238 void RandGeneral::fireArray( const int size, double* vect )
00239 {
00240 register int i;
00241
00242 for (i=0; i<size; ++i) {
00243 vect[i] = fire();
00244 }
00245 }
00246
00247 std::ostream & RandGeneral::put ( std::ostream & os ) const {
00248 int pr=os.precision(20);
00249 std::vector<unsigned long> t(2);
00250 os << " " << name() << "\n";
00251 os << "Uvec" << "\n";
00252 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n";
00253 t = DoubConv::dto2longs(oneOverNbins);
00254 os << t[0] << " " << t[1] << "\n";
00255 assert (static_cast<int>(theIntegralPdf.size())==nBins+1);
00256 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
00257 t = DoubConv::dto2longs(theIntegralPdf[i]);
00258 os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n";
00259 }
00260 os.precision(pr);
00261 return os;
00262 }
00263
00264 std::istream & RandGeneral::get ( std::istream & is ) {
00265 std::string inName;
00266 is >> inName;
00267 if (inName != name()) {
00268 is.clear(std::ios::badbit | is.rdstate());
00269 std::cerr << "Mismatch when expecting to read state of a "
00270 << name() << " distribution\n"
00271 << "Name found was " << inName
00272 << "\nistream is left in the badbit state\n";
00273 return is;
00274 }
00275 if (possibleKeywordInput(is, "Uvec", nBins)) {
00276 std::vector<unsigned long> t(2);
00277 is >> nBins >> oneOverNbins >> InterpolationType;
00278 is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t);
00279 theIntegralPdf.resize(nBins+1);
00280 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) {
00281 is >> theIntegralPdf[i] >> t[0] >> t[1];
00282 theIntegralPdf[i] = DoubConv::longs2double(t);
00283 }
00284 return is;
00285 }
00286
00287 is >> oneOverNbins >> InterpolationType;
00288 theIntegralPdf.resize(nBins+1);
00289 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i];
00290 return is;
00291 }
00292
00293 }