CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandBreitWigner.cc,v 1.4.4.2 2005/04/15 16:32:53 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandBreitWigner --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 // This file is part of Geant4 (simulation toolkit for HEP). 00010 00011 // ======================================================================= 00012 // Gabriele Cosmo - Created: 5th September 1995 00013 // - Added methods to shoot arrays: 28th July 1997 00014 // J.Marraffino - Added default arguments as attributes and 00015 // operator() with arguments: 16th Feb 1998 00016 // M Fischler - put and get to/from streams 12/10/04 00017 // M Fischler - put/get to/from streams uses pairs of ulongs when 00018 // + storing doubles avoid problems with precision 00019 // 4/14/05 00020 // ======================================================================= 00021 00022 #include "CLHEP/Random/defs.h" 00023 #include "CLHEP/Random/RandBreitWigner.h" 00024 #include "CLHEP/Units/PhysicalConstants.h" 00025 #include "CLHEP/Random/DoubConv.hh" 00026 #include <algorithm> // for min() and max() 00027 #include <cmath> 00028 00029 using namespace std; 00030 00031 namespace CLHEP { 00032 00033 std::string RandBreitWigner::name() const {return "RandBreitWigner";} 00034 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;} 00035 00036 RandBreitWigner::~RandBreitWigner() { 00037 if ( deleteEngine ) delete localEngine; 00038 } 00039 00040 RandBreitWigner::RandBreitWigner(const RandBreitWigner& right) 00041 : defaultA(right.defaultA), defaultB(right.defaultB) 00042 {;} 00043 00044 double RandBreitWigner::operator()() { 00045 return fire( defaultA, defaultB ); 00046 } 00047 00048 double RandBreitWigner::operator()( double a, double b ) { 00049 return fire( a, b ); 00050 } 00051 00052 double RandBreitWigner::operator()( double a, double b, double c ) { 00053 return fire( a, b, c ); 00054 } 00055 00056 double RandBreitWigner::shoot(double mean, double gamma) 00057 { 00058 double rval, displ; 00059 00060 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 00061 displ = 0.5*gamma*tan(rval*CLHEP::halfpi); 00062 00063 return mean + displ; 00064 } 00065 00066 double RandBreitWigner::shoot(double mean, double gamma, double cut) 00067 { 00068 double val, rval, displ; 00069 00070 if ( gamma == 0.0 ) return mean; 00071 val = atan(2.0*cut/gamma); 00072 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 00073 displ = 0.5*gamma*tan(rval*val); 00074 00075 return mean + displ; 00076 } 00077 00078 double RandBreitWigner::shootM2(double mean, double gamma ) 00079 { 00080 double val, rval, displ; 00081 00082 if ( gamma == 0.0 ) return mean; 00083 val = atan(-mean/gamma); 00084 rval = RandFlat::shoot(val, CLHEP::halfpi); 00085 displ = gamma*tan(rval); 00086 00087 return sqrt(mean*mean + mean*displ); 00088 } 00089 00090 double RandBreitWigner::shootM2(double mean, double gamma, double cut ) 00091 { 00092 double rval, displ; 00093 double lower, upper, tmp; 00094 00095 if ( gamma == 0.0 ) return mean; 00096 tmp = max(0.0,(mean-cut)); 00097 lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00098 upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00099 rval = RandFlat::shoot(lower, upper); 00100 displ = gamma*tan(rval); 00101 00102 return sqrt(max(0.0, mean*mean + mean*displ)); 00103 } 00104 00105 void RandBreitWigner::shootArray ( const int size, double* vect ) 00106 { 00107 int i; 00108 00109 for (i=0; i<size; ++i) 00110 vect[i] = shoot( 1.0, 0.2 ); 00111 } 00112 00113 void RandBreitWigner::shootArray ( const int size, double* vect, 00114 double a, double b ) 00115 { 00116 int i; 00117 00118 for (i=0; i<size; ++i) 00119 vect[i] = shoot( a, b ); 00120 } 00121 00122 void RandBreitWigner::shootArray ( const int size, double* vect, 00123 double a, double b, 00124 double c ) 00125 { 00126 int i; 00127 00128 for (i=0; i<size; ++i) 00129 vect[i] = shoot( a, b, c ); 00130 } 00131 00132 //---------------- 00133 00134 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 00135 double mean, double gamma) 00136 { 00137 double rval, displ; 00138 00139 rval = 2.0*anEngine->flat()-1.0; 00140 displ = 0.5*gamma*tan(rval*CLHEP::halfpi); 00141 00142 return mean + displ; 00143 } 00144 00145 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 00146 double mean, double gamma, double cut ) 00147 { 00148 double val, rval, displ; 00149 00150 if ( gamma == 0.0 ) return mean; 00151 val = atan(2.0*cut/gamma); 00152 rval = 2.0*anEngine->flat()-1.0; 00153 displ = 0.5*gamma*tan(rval*val); 00154 00155 return mean + displ; 00156 } 00157 00158 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 00159 double mean, double gamma ) 00160 { 00161 double val, rval, displ; 00162 00163 if ( gamma == 0.0 ) return mean; 00164 val = atan(-mean/gamma); 00165 rval = RandFlat::shoot(anEngine,val, CLHEP::halfpi); 00166 displ = gamma*tan(rval); 00167 00168 return sqrt(mean*mean + mean*displ); 00169 } 00170 00171 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 00172 double mean, double gamma, double cut ) 00173 { 00174 double rval, displ; 00175 double lower, upper, tmp; 00176 00177 if ( gamma == 0.0 ) return mean; 00178 tmp = max(0.0,(mean-cut)); 00179 lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00180 upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00181 rval = RandFlat::shoot(anEngine, lower, upper); 00182 displ = gamma*tan(rval); 00183 00184 return sqrt( max(0.0, mean*mean+mean*displ) ); 00185 } 00186 00187 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00188 const int size, double* vect ) 00189 { 00190 int i; 00191 00192 for (i=0; i<size; ++i) 00193 vect[i] = shoot( anEngine, 1.0, 0.2 ); 00194 } 00195 00196 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00197 const int size, double* vect, 00198 double a, double b ) 00199 { 00200 int i; 00201 00202 for (i=0; i<size; ++i) 00203 vect[i] = shoot( anEngine, a, b ); 00204 } 00205 00206 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00207 const int size, double* vect, 00208 double a, double b, double c ) 00209 { 00210 int i; 00211 00212 for (i=0; i<size; ++i) 00213 vect[i] = shoot( anEngine, a, b, c ); 00214 } 00215 00216 //---------------- 00217 00218 double RandBreitWigner::fire() 00219 { 00220 return fire( defaultA, defaultB ); 00221 } 00222 00223 double RandBreitWigner::fire(double mean, double gamma) 00224 { 00225 double rval, displ; 00226 00227 rval = 2.0*localEngine->flat()-1.0; 00228 displ = 0.5*gamma*tan(rval*CLHEP::halfpi); 00229 00230 return mean + displ; 00231 } 00232 00233 double RandBreitWigner::fire(double mean, double gamma, double cut) 00234 { 00235 double val, rval, displ; 00236 00237 if ( gamma == 0.0 ) return mean; 00238 val = atan(2.0*cut/gamma); 00239 rval = 2.0*localEngine->flat()-1.0; 00240 displ = 0.5*gamma*tan(rval*val); 00241 00242 return mean + displ; 00243 } 00244 00245 double RandBreitWigner::fireM2() 00246 { 00247 return fireM2( defaultA, defaultB ); 00248 } 00249 00250 double RandBreitWigner::fireM2(double mean, double gamma ) 00251 { 00252 double val, rval, displ; 00253 00254 if ( gamma == 0.0 ) return mean; 00255 val = atan(-mean/gamma); 00256 rval = RandFlat::shoot(localEngine,val, CLHEP::halfpi); 00257 displ = gamma*tan(rval); 00258 00259 return sqrt(mean*mean + mean*displ); 00260 } 00261 00262 double RandBreitWigner::fireM2(double mean, double gamma, double cut ) 00263 { 00264 double rval, displ; 00265 double lower, upper, tmp; 00266 00267 if ( gamma == 0.0 ) return mean; 00268 tmp = max(0.0,(mean-cut)); 00269 lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00270 upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00271 rval = RandFlat::shoot(localEngine,lower, upper); 00272 displ = gamma*tan(rval); 00273 00274 return sqrt(max(0.0, mean*mean + mean*displ)); 00275 } 00276 00277 void RandBreitWigner::fireArray ( const int size, double* vect) 00278 { 00279 int i; 00280 00281 for (i=0; i<size; ++i) 00282 vect[i] = fire(defaultA, defaultB ); 00283 } 00284 00285 void RandBreitWigner::fireArray ( const int size, double* vect, 00286 double a, double b ) 00287 { 00288 int i; 00289 00290 for (i=0; i<size; ++i) 00291 vect[i] = fire( a, b ); 00292 } 00293 00294 void RandBreitWigner::fireArray ( const int size, double* vect, 00295 double a, double b, double c ) 00296 { 00297 int i; 00298 00299 for (i=0; i<size; ++i) 00300 vect[i] = fire( a, b, c ); 00301 } 00302 00303 00304 std::ostream & RandBreitWigner::put ( std::ostream & os ) const { 00305 int pr=os.precision(20); 00306 std::vector<unsigned long> t(2); 00307 os << " " << name() << "\n"; 00308 os << "Uvec" << "\n"; 00309 t = DoubConv::dto2longs(defaultA); 00310 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00311 t = DoubConv::dto2longs(defaultB); 00312 os << defaultB << " " << t[0] << " " << t[1] << "\n"; 00313 os.precision(pr); 00314 return os; 00315 #ifdef REMOVED 00316 int pr=os.precision(20); 00317 os << " " << name() << "\n"; 00318 os << defaultA << " " << defaultB << "\n"; 00319 os.precision(pr); 00320 return os; 00321 #endif 00322 } 00323 00324 std::istream & RandBreitWigner::get ( std::istream & is ) { 00325 std::string inName; 00326 is >> inName; 00327 if (inName != name()) { 00328 is.clear(std::ios::badbit | is.rdstate()); 00329 std::cerr << "Mismatch when expecting to read state of a " 00330 << name() << " distribution\n" 00331 << "Name found was " << inName 00332 << "\nistream is left in the badbit state\n"; 00333 return is; 00334 } 00335 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00336 std::vector<unsigned long> t(2); 00337 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00338 is >> defaultB >> t[0] >> t[1]; defaultB = DoubConv::longs2double(t); 00339 return is; 00340 } 00341 // is >> defaultA encompassed by possibleKeywordInput 00342 is >> defaultB; 00343 return is; 00344 } 00345 00346 00347 } // namespace CLHEP 00348