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