CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandBreitWigner.cc,v 1.6 2010/06/16 17:24: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 namespace CLHEP { 00030 00031 std::string RandBreitWigner::name() const {return "RandBreitWigner";} 00032 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;} 00033 00034 RandBreitWigner::~RandBreitWigner() { 00035 } 00036 00037 double RandBreitWigner::operator()() { 00038 return fire( defaultA, defaultB ); 00039 } 00040 00041 double RandBreitWigner::operator()( double a, double b ) { 00042 return fire( a, b ); 00043 } 00044 00045 double RandBreitWigner::operator()( double a, double b, double c ) { 00046 return fire( a, b, c ); 00047 } 00048 00049 double RandBreitWigner::shoot(double mean, double gamma) 00050 { 00051 double rval, displ; 00052 00053 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 00054 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 00055 00056 return mean + displ; 00057 } 00058 00059 double RandBreitWigner::shoot(double mean, double gamma, double cut) 00060 { 00061 double val, rval, displ; 00062 00063 if ( gamma == 0.0 ) return mean; 00064 val = std::atan(2.0*cut/gamma); 00065 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 00066 displ = 0.5*gamma*std::tan(rval*val); 00067 00068 return mean + displ; 00069 } 00070 00071 double RandBreitWigner::shootM2(double mean, double gamma ) 00072 { 00073 double val, rval, displ; 00074 00075 if ( gamma == 0.0 ) return mean; 00076 val = std::atan(-mean/gamma); 00077 rval = RandFlat::shoot(val, CLHEP::halfpi); 00078 displ = gamma*std::tan(rval); 00079 00080 return std::sqrt(mean*mean + mean*displ); 00081 } 00082 00083 double RandBreitWigner::shootM2(double mean, double gamma, double cut ) 00084 { 00085 double rval, displ; 00086 double lower, upper, tmp; 00087 00088 if ( gamma == 0.0 ) return mean; 00089 tmp = std::max(0.0,(mean-cut)); 00090 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00091 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00092 rval = RandFlat::shoot(lower, upper); 00093 displ = gamma*std::tan(rval); 00094 00095 return std::sqrt(std::max(0.0, mean*mean + mean*displ)); 00096 } 00097 00098 void RandBreitWigner::shootArray ( const int size, double* vect ) 00099 { 00100 for( double* v = vect; v != vect + size; ++v ) 00101 *v = shoot( 1.0, 0.2 ); 00102 } 00103 00104 void RandBreitWigner::shootArray ( const int size, double* vect, 00105 double a, double b ) 00106 { 00107 for( double* v = vect; v != vect + size; ++v ) 00108 *v = shoot( a, b ); 00109 } 00110 00111 void RandBreitWigner::shootArray ( const int size, double* vect, 00112 double a, double b, 00113 double c ) 00114 { 00115 for( double* v = vect; v != vect + size; ++v ) 00116 *v = shoot( a, b, c ); 00117 } 00118 00119 //---------------- 00120 00121 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 00122 double mean, double gamma) 00123 { 00124 double rval, displ; 00125 00126 rval = 2.0*anEngine->flat()-1.0; 00127 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 00128 00129 return mean + displ; 00130 } 00131 00132 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 00133 double mean, double gamma, double cut ) 00134 { 00135 double val, rval, displ; 00136 00137 if ( gamma == 0.0 ) return mean; 00138 val = std::atan(2.0*cut/gamma); 00139 rval = 2.0*anEngine->flat()-1.0; 00140 displ = 0.5*gamma*std::tan(rval*val); 00141 00142 return mean + displ; 00143 } 00144 00145 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 00146 double mean, double gamma ) 00147 { 00148 double val, rval, displ; 00149 00150 if ( gamma == 0.0 ) return mean; 00151 val = std::atan(-mean/gamma); 00152 rval = RandFlat::shoot(anEngine,val, CLHEP::halfpi); 00153 displ = gamma*std::tan(rval); 00154 00155 return std::sqrt(mean*mean + mean*displ); 00156 } 00157 00158 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 00159 double mean, double gamma, double cut ) 00160 { 00161 double rval, displ; 00162 double lower, upper, tmp; 00163 00164 if ( gamma == 0.0 ) return mean; 00165 tmp = std::max(0.0,(mean-cut)); 00166 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00167 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00168 rval = RandFlat::shoot(anEngine, lower, upper); 00169 displ = gamma*std::tan(rval); 00170 00171 return std::sqrt( std::max(0.0, mean*mean+mean*displ) ); 00172 } 00173 00174 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00175 const int size, double* vect ) 00176 { 00177 for( double* v = vect; v != vect + size; ++v ) 00178 *v = shoot( anEngine, 1.0, 0.2 ); 00179 } 00180 00181 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00182 const int size, double* vect, 00183 double a, double b ) 00184 { 00185 for( double* v = vect; v != vect + size; ++v ) 00186 *v = shoot( anEngine, a, b ); 00187 } 00188 00189 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 00190 const int size, double* vect, 00191 double a, double b, double c ) 00192 { 00193 for( double* v = vect; v != vect + size; ++v ) 00194 *v = shoot( anEngine, a, b, c ); 00195 } 00196 00197 //---------------- 00198 00199 double RandBreitWigner::fire() 00200 { 00201 return fire( defaultA, defaultB ); 00202 } 00203 00204 double RandBreitWigner::fire(double mean, double gamma) 00205 { 00206 double rval, displ; 00207 00208 rval = 2.0*localEngine->flat()-1.0; 00209 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 00210 00211 return mean + displ; 00212 } 00213 00214 double RandBreitWigner::fire(double mean, double gamma, double cut) 00215 { 00216 double val, rval, displ; 00217 00218 if ( gamma == 0.0 ) return mean; 00219 val = std::atan(2.0*cut/gamma); 00220 rval = 2.0*localEngine->flat()-1.0; 00221 displ = 0.5*gamma*std::tan(rval*val); 00222 00223 return mean + displ; 00224 } 00225 00226 double RandBreitWigner::fireM2() 00227 { 00228 return fireM2( defaultA, defaultB ); 00229 } 00230 00231 double RandBreitWigner::fireM2(double mean, double gamma ) 00232 { 00233 double val, rval, displ; 00234 00235 if ( gamma == 0.0 ) return mean; 00236 val = std::atan(-mean/gamma); 00237 rval = RandFlat::shoot(localEngine.get(),val, CLHEP::halfpi); 00238 displ = gamma*std::tan(rval); 00239 00240 return std::sqrt(mean*mean + mean*displ); 00241 } 00242 00243 double RandBreitWigner::fireM2(double mean, double gamma, double cut ) 00244 { 00245 double rval, displ; 00246 double lower, upper, tmp; 00247 00248 if ( gamma == 0.0 ) return mean; 00249 tmp = std::max(0.0,(mean-cut)); 00250 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 00251 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 00252 rval = RandFlat::shoot(localEngine.get(),lower, upper); 00253 displ = gamma*std::tan(rval); 00254 00255 return std::sqrt(std::max(0.0, mean*mean + mean*displ)); 00256 } 00257 00258 void RandBreitWigner::fireArray ( const int size, double* vect) 00259 { 00260 for( double* v = vect; v != vect + size; ++v ) 00261 *v = fire(defaultA, defaultB ); 00262 } 00263 00264 void RandBreitWigner::fireArray ( const int size, double* vect, 00265 double a, double b ) 00266 { 00267 for( double* v = vect; v != vect + size; ++v ) 00268 *v = fire( a, b ); 00269 } 00270 00271 void RandBreitWigner::fireArray ( const int size, double* vect, 00272 double a, double b, double c ) 00273 { 00274 for( double* v = vect; v != vect + size; ++v ) 00275 *v = fire( a, b, c ); 00276 } 00277 00278 00279 std::ostream & RandBreitWigner::put ( std::ostream & os ) const { 00280 int pr=os.precision(20); 00281 std::vector<unsigned long> t(2); 00282 os << " " << name() << "\n"; 00283 os << "Uvec" << "\n"; 00284 t = DoubConv::dto2longs(defaultA); 00285 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 00286 t = DoubConv::dto2longs(defaultB); 00287 os << defaultB << " " << t[0] << " " << t[1] << "\n"; 00288 os.precision(pr); 00289 return os; 00290 #ifdef REMOVED 00291 int pr=os.precision(20); 00292 os << " " << name() << "\n"; 00293 os << defaultA << " " << defaultB << "\n"; 00294 os.precision(pr); 00295 return os; 00296 #endif 00297 } 00298 00299 std::istream & RandBreitWigner::get ( std::istream & is ) { 00300 std::string inName; 00301 is >> inName; 00302 if (inName != name()) { 00303 is.clear(std::ios::badbit | is.rdstate()); 00304 std::cerr << "Mismatch when expecting to read state of a " 00305 << name() << " distribution\n" 00306 << "Name found was " << inName 00307 << "\nistream is left in the badbit state\n"; 00308 return is; 00309 } 00310 if (possibleKeywordInput(is, "Uvec", defaultA)) { 00311 std::vector<unsigned long> t(2); 00312 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 00313 is >> defaultB >> t[0] >> t[1]; defaultB = DoubConv::longs2double(t); 00314 return is; 00315 } 00316 // is >> defaultA encompassed by possibleKeywordInput 00317 is >> defaultB; 00318 return is; 00319 } 00320 00321 00322 } // namespace CLHEP 00323