CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandGeneral.cc,v 1.4.4.2 2005/04/15 16:32: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 <assert.h> 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 : deleteEngine(false), 00069 nBins(theProbSize), 00070 InterpolationType(IntType) 00071 { 00072 localEngine = HepRandom::getTheEngine(); 00073 prepareTable(aProbFunc); 00074 } 00075 00076 RandGeneral::RandGeneral(HepRandomEngine& anEngine, 00077 const double* aProbFunc, 00078 int theProbSize, 00079 int IntType ) 00080 : localEngine(&anEngine), 00081 deleteEngine(false), 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 : localEngine(anEngine), 00093 deleteEngine(true), 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 if ( deleteEngine ) delete localEngine; 00175 } 00176 00177 00179 // mapRandom(rand) 00181 00182 double RandGeneral::mapRandom(double rand) const { 00183 // 00184 // Private method to take the random (however it is created) and map it 00185 // according to the distribution. 00186 // 00187 00188 int nbelow = 0; // largest k such that I[k] is known to be <= rand 00189 int nabove = nBins; // largest k such that I[k] is known to be > rand 00190 int middle; 00191 00192 while (nabove > nbelow+1) { 00193 middle = (nabove + nbelow+1)>>1; 00194 if (rand >= theIntegralPdf[middle]) { 00195 nbelow = middle; 00196 } else { 00197 nabove = middle; 00198 } 00199 } // after this loop, nabove is always nbelow+1 and they straddle rad: 00200 assert ( nabove == nbelow+1 ); 00201 assert ( theIntegralPdf[nbelow] <= rand ); 00202 assert ( theIntegralPdf[nabove] >= rand ); 00203 // If a defective engine produces rand=1, that will 00204 // still give sensible results so we relax the > rand assertion 00205 00206 if ( InterpolationType == 1 ) { 00207 00208 return nbelow * oneOverNbins; 00209 00210 } else { 00211 00212 double binMeasure = theIntegralPdf[nabove] - theIntegralPdf[nbelow]; 00213 // binMeasure is always aProbFunc[nbelow], 00214 // but we don't have aProbFunc any more so we subtract. 00215 00216 if ( binMeasure == 0 ) { 00217 // rand lies right in a bin of measure 0. Simply return the center 00218 // of the range of that bin. (Any value between k/N and (k+1)/N is 00219 // equally good, in this rare case.) 00220 return (nbelow + .5) * oneOverNbins; 00221 } 00222 00223 double binFraction = (rand - theIntegralPdf[nbelow]) / binMeasure; 00224 00225 return (nbelow + binFraction) * oneOverNbins; 00226 } 00227 00228 } // mapRandom(rand) 00229 00230 00231 00232 void RandGeneral::shootArray( HepRandomEngine* anEngine, 00233 const int size, double* vect ) 00234 { 00235 register int i; 00236 00237 for (i=0; i<size; ++i) { 00238 vect[i] = shoot(anEngine); 00239 } 00240 } 00241 00242 void RandGeneral::fireArray( const int size, double* vect ) 00243 { 00244 register int i; 00245 00246 for (i=0; i<size; ++i) { 00247 vect[i] = fire(); 00248 } 00249 } 00250 00251 std::ostream & RandGeneral::put ( std::ostream & os ) const { 00252 int pr=os.precision(20); 00253 std::vector<unsigned long> t(2); 00254 os << " " << name() << "\n"; 00255 os << "Uvec" << "\n"; 00256 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n"; 00257 t = DoubConv::dto2longs(oneOverNbins); 00258 os << t[0] << " " << t[1] << "\n"; 00259 assert (static_cast<int>(theIntegralPdf.size())==nBins+1); 00260 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) { 00261 t = DoubConv::dto2longs(theIntegralPdf[i]); 00262 os << theIntegralPdf[i] << " " << t[0] << " " << t[1] << "\n"; 00263 } 00264 os.precision(pr); 00265 return os; 00266 #ifdef REMOVED 00267 int pr=os.precision(20); 00268 os << " " << name() << "\n"; 00269 os << nBins << " " << oneOverNbins << " " << InterpolationType << "\n"; 00270 assert (static_cast<int>(theIntegralPdf.size())==nBins+1); 00271 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) 00272 os << theIntegralPdf[i] << "\n"; 00273 os.precision(pr); 00274 return os; 00275 #endif 00276 } 00277 00278 std::istream & RandGeneral::get ( std::istream & is ) { 00279 std::string inName; 00280 is >> inName; 00281 if (inName != name()) { 00282 is.clear(std::ios::badbit | is.rdstate()); 00283 std::cerr << "Mismatch when expecting to read state of a " 00284 << name() << " distribution\n" 00285 << "Name found was " << inName 00286 << "\nistream is left in the badbit state\n"; 00287 return is; 00288 } 00289 if (possibleKeywordInput(is, "Uvec", nBins)) { 00290 std::vector<unsigned long> t(2); 00291 is >> nBins >> oneOverNbins >> InterpolationType; 00292 is >> t[0] >> t[1]; oneOverNbins = DoubConv::longs2double(t); 00293 theIntegralPdf.resize(nBins+1); 00294 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) { 00295 is >> theIntegralPdf[i] >> t[0] >> t[1]; 00296 theIntegralPdf[i] = DoubConv::longs2double(t); 00297 } 00298 return is; 00299 } 00300 // is >> nBins encompassed by possibleKeywordInput 00301 is >> oneOverNbins >> InterpolationType; 00302 theIntegralPdf.resize(nBins+1); 00303 for (unsigned int i=0; i<theIntegralPdf.size(); ++i) is >> theIntegralPdf[i]; 00304 return is; 00305 } 00306 00307 } // namespace CLHEP