CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGeneral.cc,v 1.6 2010/06/16 17:24:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandGeneral --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // S.Magni & G.Pieri - Created: 5th September 1995 00012 // G.Cosmo - Added constructor using default engine from the 00013 // static generator. Simplified shoot() and 00014 // shootArray() (not needed in principle!): 20th Aug 1998 00015 // M.G.Pia & G.Cosmo - Fixed bug in computation of theIntegralPdf in 00016 // two constructors: 5th Jan 1999 00017 // S.Magni & G.Pieri - Added linear interpolation: 24th Mar 1999 00018 // M. Fischler - General cleanup: 14th May 1999 00019 // + Eliminated constructor code replication by factoring 00020 // common code into prepareTable. 00021 // + Eliminated fire/shoot code replication by factoring 00022 // out common code into mapRandom. 00023 // + A couple of methods are moved inline to avoid a 00024 // speed cost for factoring out mapRandom: fire() 00025 // and shoot(anEngine). 00026 // + Inserted checks for negative weight and zero total 00027 // weight in the bins. 00028 // + Modified the binary search loop to avoid incorrect 00029 // behavior when rand is example on a boundary. 00030 // + Moved the check of InterpolationType up into 00031 // the constructor. A type other than 0 or 1 00032 // will give the interpolated distribution (instead of 00033 // a distribution that always returns 0). 00034 // + Modified the computation of the returned value 00035 // to use algeraic simplification to improve speed. 00036 // Eliminated two of the three divisionns, made 00037 // use of the fact that nabove-nbelow is always 1, etc. 00038 // + Inserted a check for rand hitting the boundary of a 00039 // zero-width bin, to avoid dividing 0/0. 00040 // M. Fischler - Minor correction in assert 31 July 2001 00041 // + changed from assert (above = below+1) to == 00042 // M Fischler - put and get to/from streams 12/15/04 00043 // + Modifications to use a vector as theIntegraPdf 00044 // M Fischler - put/get to/from streams uses pairs of ulongs when 00045 // + storing doubles avoid problems with precision 00046 // 4/14/05 00047 // 00048 // ======================================================================= 00049 00050 #include "CLHEP/Random/defs.h" 00051 #include "CLHEP/Random/RandGeneral.h" 00052 #include "CLHEP/Random/DoubConv.hh" 00053 #include <cassert> 00054 00055 namespace CLHEP { 00056 00057 std::string RandGeneral::name() const {return "RandGeneral";} 00058 HepRandomEngine & RandGeneral::engine() {return *localEngine;} 00059 00060 00062 // Constructors 00064 00065 RandGeneral::RandGeneral( const double* aProbFunc, 00066 int theProbSize, 00067 int IntType ) 00068 : HepRandom(), 00069 localEngine(HepRandom::getTheEngine(), do_nothing_deleter()), 00070 nBins(theProbSize), 00071 InterpolationType(IntType) 00072 { 00073 prepareTable(aProbFunc); 00074 } 00075 00076 RandGeneral::RandGeneral(HepRandomEngine& anEngine, 00077 const double* aProbFunc, 00078 int theProbSize, 00079 int IntType ) 00080 : HepRandom(), 00081 localEngine(&anEngine, do_nothing_deleter()), 00082 nBins(theProbSize), 00083 InterpolationType(IntType) 00084 { 00085 prepareTable(aProbFunc); 00086 } 00087 00088 RandGeneral::RandGeneral(HepRandomEngine* anEngine, 00089 const double* aProbFunc, 00090 int theProbSize, 00091 int IntType ) 00092 : HepRandom(), 00093 localEngine(anEngine), 00094 nBins(theProbSize), 00095 InterpolationType(IntType) 00096 { 00097 prepareTable(aProbFunc); 00098 } 00099 00100 void RandGeneral::prepareTable(const double* aProbFunc) { 00101 // 00102 // Private method called only by constructors. Prepares theIntegralPdf. 00103 // 00104 if (nBins < 1) { 00105 std::cerr << 00106 "RandGeneral constructed with no bins - will use flat distribution\n"; 00107 useFlatDistribution(); 00108 return; 00109 } 00110 00111 theIntegralPdf.resize(nBins+1); 00112 theIntegralPdf[0] = 0; 00113 register int ptn; 00114 register double weight; 00115 00116 for ( ptn = 0; ptn<nBins; ++ptn ) { 00117 weight = aProbFunc[ptn]; 00118 if ( weight < 0 ) { 00119 // We can't stomach negative bin contents, they invalidate the 00120 // search algorithm when the distribution is fired. 00121 std::cerr << 00122 "RandGeneral constructed with negative-weight bin " << ptn << 00123 " = " << weight << " \n -- will substitute 0 weight \n"; 00124 weight = 0; 00125 } 00126 // std::cout << ptn << " " << weight << " " << theIntegralPdf[ptn] << "\n"; 00127 theIntegralPdf[ptn+1] = theIntegralPdf[ptn] + weight; 00128 } 00129 00130 if ( theIntegralPdf[nBins] <= 0 ) { 00131 std::cerr << 00132 "RandGeneral constructed nothing in bins - will use flat distribution\n"; 00133 useFlatDistribution(); 00134 return; 00135 } 00136 00137 for ( ptn = 0; ptn < nBins+1; ++ptn ) { 00138 theIntegralPdf[ptn] /= theIntegralPdf[nBins]; 00139 // std::cout << ptn << " " << theIntegralPdf[ptn] << "\n"; 00140 } 00141 00142 // And another useful variable is ... 00143 oneOverNbins = 1.0 / nBins; 00144 00145 // One last chore: 00146 00147 if ( (InterpolationType != 0) && (InterpolationType != 1) ) { 00148 std::cerr << 00149 "RandGeneral does not recognize IntType " << InterpolationType 00150 << "\n Will use type 0 (continuous linear interpolation \n"; 00151 InterpolationType = 0; 00152 } 00153 00154 } // prepareTable() 00155 00156 void RandGeneral::useFlatDistribution() { 00157 // 00158 // Private method called only by prepareTables in case of user error. 00159 // 00160 nBins = 1; 00161 theIntegralPdf.resize(2); 00162 theIntegralPdf[0] = 0; 00163 theIntegralPdf[1] = 1; 00164 oneOverNbins = 1.0; 00165 return; 00166 00167 } // UseFlatDistribution() 00168 00170 // Destructor 00172 00173 RandGeneral::~RandGeneral() { 00174 } 00175 00176 00178 // mapRandom(rand) 00180 00181 double RandGeneral::mapRandom(double rand) const { 00182 // 00183 // Private method to take the random (however it is created) and map it 00184 // according to the distribution. 00185 // 00186 00187 int nbelow = 0; // largest k such that I[k] is known to be <= rand 00188 int nabove = nBins; // largest k such that I[k] is known to be > rand 00189 int middle; 00190 00191 while (nabove > nbelow+1) { 00192 middle = (nabove + nbelow+1)>>1; 00193 if (rand >= theIntegralPdf[middle]) { 00194 nbelow = middle; 00195 } else { 00196 nabove = middle; 00197 } 00198 } // after this loop, nabove is always nbelow+1 and they straddle rad: 00199 assert ( nabove == nbelow+1 ); 00200 assert ( theIntegralPdf[nbelow] <= rand ); 00201 assert ( theIntegralPdf[nabove] >= rand ); 00202 // If a defective engine produces rand=1, that will 00203 // still give sensible results so we relax the > rand assertion 00204 00205 if ( InterpolationType == 1 ) { 00206 00207 return nbelow * oneOverNbins; 00208 00209 } else { 00210 00211 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow]; 00212 // binMeasure is always aProbFunc[nbelow], 00213 // but we don't have aProbFunc any more so we subtract. 00214 00215 if ( binMeasure == 0 ) { 00216 // rand lies right in a bin of measure 0. Simply return the center 00217 // of the range of that bin. (Any value between k/N and (k+1)/N is 00218 // equally good, in this rare case.) 00219 return (nbelow + .5) * oneOverNbins; 00220 } 00221 00222 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure; 00223 00224 return (nbelow + binFraction) * oneOverNbins; 00225 } 00226 00227 } // mapRandom(rand) 00228 00229 00230 00231 void RandGeneral::shootArray( HepRandomEngine* anEngine, 00232 const int size, double* vect ) 00233 { 00234 register int i; 00235 00236 for (i=0; i<size; ++i) { 00237 vect[i] = shoot(anEngine); 00238 } 00239 } 00240 00241 void RandGeneral::fireArray( const int size, double* vect ) 00242 { 00243 register int i; 00244 00245 for (i=0; i<size; ++i) { 00246 vect[i] = fire(); 00247 } 00248 } 00249 00250 std::ostream & RandGeneral::put ( std::ostream & os ) const { 00251 int pr=os.precision(20); 00252 std::vector<unsigned long> t(2); 00253 os << " " << name() << "\n"; 00254 os << "Uvec" << "\n"; 00255 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n"; 00256 t = DoubConv::dto2longs(oneOverNbins); 00257 os << t[0] << " " << t[1] << "\n"; 00258 assert (static_cast<int>(theIntegralPdf.size())==nBins+1); 00259 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) { 00260 t = DoubConv::dto2longs(theIntegralPdf[i]); 00261 os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n"; 00262 } 00263 os.precision(pr); 00264 return os; 00265 #ifdef REMOVED 00266 int pr=os.precision(20); 00267 os << " " << name() << "\n"; 00268 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n"; 00269 assert (static_cast<int>(theIntegralPdf.size())==nBins+1); 00270 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) 00271 os << theIntegralPdf[i] << "\n"; 00272 os.precision(pr); 00273 return os; 00274 #endif 00275 } 00276 00277 std::istream & RandGeneral::get ( std::istream & is ) { 00278 std::string inName; 00279 is >> inName; 00280 if (inName != name()) { 00281 is.clear(std::ios::badbit | is.rdstate()); 00282 std::cerr << "Mismatch when expecting to read state of a " 00283 << name() << " distribution\n" 00284 << "Name found was " << inName 00285 << "\nistream is left in the badbit state\n"; 00286 return is; 00287 } 00288 if (possibleKeywordInput(is, "Uvec", nBins)) { 00289 std::vector<unsigned long> t(2); 00290 is >> nBins >> oneOverNbins >> InterpolationType; 00291 is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t); 00292 theIntegralPdf.resize(nBins+1); 00293 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) { 00294 is >> theIntegralPdf[i] >> t[0] >> t[1]; 00295 theIntegralPdf[i] = DoubConv::longs2double(t); 00296 } 00297 return is; 00298 } 00299 // is >> nBins encompassed by possibleKeywordInput 00300 is >> oneOverNbins >> InterpolationType; 00301 theIntegralPdf.resize(nBins+1); 00302 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i]; 00303 return is; 00304 } 00305 00306 } // namespace CLHEP