Geant4 Cross Reference |
>> 1 // $Id:$ 1 // -*- C++ -*- 2 // -*- C++ -*- 2 // 3 // 3 // ------------------------------------------- 4 // ----------------------------------------------------------------------- 4 // HEP Random 5 // HEP Random 5 // --- RandBreitWigner - 6 // --- RandBreitWigner --- 6 // class implementation f 7 // class implementation file 7 // ------------------------------------------- 8 // ----------------------------------------------------------------------- 8 // This file is part of Geant4 (simulation too 9 // This file is part of Geant4 (simulation toolkit for HEP). 9 10 10 // =========================================== 11 // ======================================================================= 11 // Gabriele Cosmo - Created: 5th September 199 12 // Gabriele Cosmo - Created: 5th September 1995 12 // - Added methods to shoot arr 13 // - Added methods to shoot arrays: 28th July 1997 13 // J.Marraffino - Added default arguments as 14 // J.Marraffino - Added default arguments as attributes and 14 // operator() with arguments: 15 // operator() with arguments: 16th Feb 1998 15 // M Fischler - put and get to/from stream 16 // M Fischler - put and get to/from streams 12/10/04 16 // M Fischler - put/get to/from streams 17 // M Fischler - put/get to/from streams uses pairs of ulongs when 17 // + storing doubles avoid problems with 18 // + storing doubles avoid problems with precision 18 // 4/14/05 19 // 4/14/05 19 // =========================================== 20 // ======================================================================= 20 21 21 #include "CLHEP/Random/RandBreitWigner.h" 22 #include "CLHEP/Random/RandBreitWigner.h" 22 #include "CLHEP/Units/PhysicalConstants.h" 23 #include "CLHEP/Units/PhysicalConstants.h" 23 #include "CLHEP/Random/DoubConv.h" 24 #include "CLHEP/Random/DoubConv.h" 24 #include <algorithm> // for min() and max() 25 #include <algorithm> // for min() and max() 25 #include <cmath> 26 #include <cmath> 26 #include <iostream> << 27 #include <string> << 28 #include <vector> << 29 27 30 namespace CLHEP { 28 namespace CLHEP { 31 29 32 std::string RandBreitWigner::name() const {ret 30 std::string RandBreitWigner::name() const {return "RandBreitWigner";} 33 HepRandomEngine & RandBreitWigner::engine() {r 31 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;} 34 32 35 RandBreitWigner::~RandBreitWigner() { 33 RandBreitWigner::~RandBreitWigner() { 36 } 34 } 37 35 38 double RandBreitWigner::operator()() { 36 double RandBreitWigner::operator()() { 39 return fire( defaultA, defaultB ); 37 return fire( defaultA, defaultB ); 40 } 38 } 41 39 42 double RandBreitWigner::operator()( double a, 40 double RandBreitWigner::operator()( double a, double b ) { 43 return fire( a, b ); 41 return fire( a, b ); 44 } 42 } 45 43 46 double RandBreitWigner::operator()( double a, 44 double RandBreitWigner::operator()( double a, double b, double c ) { 47 return fire( a, b, c ); 45 return fire( a, b, c ); 48 } 46 } 49 47 50 double RandBreitWigner::shoot(double mean, dou 48 double RandBreitWigner::shoot(double mean, double gamma) 51 { 49 { 52 double rval, displ; 50 double rval, displ; 53 51 54 rval = 2.0*HepRandom::getTheEngine()->flat( 52 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 55 displ = 0.5*gamma*std::tan(rval*CLHEP::half 53 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 56 54 57 return mean + displ; 55 return mean + displ; 58 } 56 } 59 57 60 double RandBreitWigner::shoot(double mean, dou 58 double RandBreitWigner::shoot(double mean, double gamma, double cut) 61 { 59 { 62 double val, rval, displ; 60 double val, rval, displ; 63 61 64 if ( gamma == 0.0 ) return mean; 62 if ( gamma == 0.0 ) return mean; 65 val = std::atan(2.0*cut/gamma); 63 val = std::atan(2.0*cut/gamma); 66 rval = 2.0*HepRandom::getTheEngine()->flat( 64 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0; 67 displ = 0.5*gamma*std::tan(rval*val); 65 displ = 0.5*gamma*std::tan(rval*val); 68 66 69 return mean + displ; 67 return mean + displ; 70 } 68 } 71 69 72 double RandBreitWigner::shootM2(double mean, d 70 double RandBreitWigner::shootM2(double mean, double gamma ) 73 { 71 { 74 double val, rval, displ; 72 double val, rval, displ; 75 73 76 if ( gamma == 0.0 ) return mean; 74 if ( gamma == 0.0 ) return mean; 77 val = std::atan(-mean/gamma); 75 val = std::atan(-mean/gamma); 78 rval = RandFlat::shoot(val, CLHEP::halfpi); 76 rval = RandFlat::shoot(val, CLHEP::halfpi); 79 displ = gamma*std::tan(rval); 77 displ = gamma*std::tan(rval); 80 78 81 return std::sqrt(mean*mean + mean*displ); 79 return std::sqrt(mean*mean + mean*displ); 82 } 80 } 83 81 84 double RandBreitWigner::shootM2(double mean, d 82 double RandBreitWigner::shootM2(double mean, double gamma, double cut ) 85 { 83 { 86 double rval, displ; 84 double rval, displ; 87 double lower, upper, tmp; 85 double lower, upper, tmp; 88 86 89 if ( gamma == 0.0 ) return mean; 87 if ( gamma == 0.0 ) return mean; 90 tmp = std::max(0.0,(mean-cut)); 88 tmp = std::max(0.0,(mean-cut)); 91 lower = std::atan( (tmp*tmp-mean*mean)/(mea 89 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 92 upper = std::atan( ((mean+cut)*(mean+cut)-m 90 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 93 rval = RandFlat::shoot(lower, upper); 91 rval = RandFlat::shoot(lower, upper); 94 displ = gamma*std::tan(rval); 92 displ = gamma*std::tan(rval); 95 93 96 return std::sqrt(std::max(0.0, mean*mean + 94 return std::sqrt(std::max(0.0, mean*mean + mean*displ)); 97 } 95 } 98 96 99 void RandBreitWigner::shootArray ( const int s 97 void RandBreitWigner::shootArray ( const int size, double* vect ) 100 { 98 { 101 for( double* v = vect; v != vect + size; ++v 99 for( double* v = vect; v != vect + size; ++v ) 102 *v = shoot( 1.0, 0.2 ); 100 *v = shoot( 1.0, 0.2 ); 103 } 101 } 104 102 105 void RandBreitWigner::shootArray ( const int s 103 void RandBreitWigner::shootArray ( const int size, double* vect, 106 double a, d 104 double a, double b ) 107 { 105 { 108 for( double* v = vect; v != vect + size; ++v 106 for( double* v = vect; v != vect + size; ++v ) 109 *v = shoot( a, b ); 107 *v = shoot( a, b ); 110 } 108 } 111 109 112 void RandBreitWigner::shootArray ( const int s 110 void RandBreitWigner::shootArray ( const int size, double* vect, 113 double a, d 111 double a, double b, 114 double c ) 112 double c ) 115 { 113 { 116 for( double* v = vect; v != vect + size; ++v 114 for( double* v = vect; v != vect + size; ++v ) 117 *v = shoot( a, b, c ); 115 *v = shoot( a, b, c ); 118 } 116 } 119 117 120 //---------------- 118 //---------------- 121 119 122 double RandBreitWigner::shoot(HepRandomEngine* 120 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 123 double mean, 121 double mean, double gamma) 124 { 122 { 125 double rval, displ; 123 double rval, displ; 126 124 127 rval = 2.0*anEngine->flat()-1.0; 125 rval = 2.0*anEngine->flat()-1.0; 128 displ = 0.5*gamma*std::tan(rval*CLHEP::half 126 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 129 127 130 return mean + displ; 128 return mean + displ; 131 } 129 } 132 130 133 double RandBreitWigner::shoot(HepRandomEngine* 131 double RandBreitWigner::shoot(HepRandomEngine* anEngine, 134 double mean, 132 double mean, double gamma, double cut ) 135 { 133 { 136 double val, rval, displ; 134 double val, rval, displ; 137 135 138 if ( gamma == 0.0 ) return mean; 136 if ( gamma == 0.0 ) return mean; 139 val = std::atan(2.0*cut/gamma); 137 val = std::atan(2.0*cut/gamma); 140 rval = 2.0*anEngine->flat()-1.0; 138 rval = 2.0*anEngine->flat()-1.0; 141 displ = 0.5*gamma*std::tan(rval*val); 139 displ = 0.5*gamma*std::tan(rval*val); 142 140 143 return mean + displ; 141 return mean + displ; 144 } 142 } 145 143 146 double RandBreitWigner::shootM2(HepRandomEngin 144 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 147 double mean 145 double mean, double gamma ) 148 { 146 { 149 double val, rval, displ; 147 double val, rval, displ; 150 148 151 if ( gamma == 0.0 ) return mean; 149 if ( gamma == 0.0 ) return mean; 152 val = std::atan(-mean/gamma); 150 val = std::atan(-mean/gamma); 153 rval = RandFlat::shoot(anEngine,val, CLHEP: 151 rval = RandFlat::shoot(anEngine,val, CLHEP::halfpi); 154 displ = gamma*std::tan(rval); 152 displ = gamma*std::tan(rval); 155 153 156 return std::sqrt(mean*mean + mean*displ); 154 return std::sqrt(mean*mean + mean*displ); 157 } 155 } 158 156 159 double RandBreitWigner::shootM2(HepRandomEngin 157 double RandBreitWigner::shootM2(HepRandomEngine* anEngine, 160 double mean 158 double mean, double gamma, double cut ) 161 { 159 { 162 double rval, displ; 160 double rval, displ; 163 double lower, upper, tmp; 161 double lower, upper, tmp; 164 162 165 if ( gamma == 0.0 ) return mean; 163 if ( gamma == 0.0 ) return mean; 166 tmp = std::max(0.0,(mean-cut)); 164 tmp = std::max(0.0,(mean-cut)); 167 lower = std::atan( (tmp*tmp-mean*mean)/(mea 165 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 168 upper = std::atan( ((mean+cut)*(mean+cut)-m 166 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 169 rval = RandFlat::shoot(anEngine, lower, upp 167 rval = RandFlat::shoot(anEngine, lower, upper); 170 displ = gamma*std::tan(rval); 168 displ = gamma*std::tan(rval); 171 169 172 return std::sqrt( std::max(0.0, mean*mean+m 170 return std::sqrt( std::max(0.0, mean*mean+mean*displ) ); 173 } 171 } 174 172 175 void RandBreitWigner::shootArray ( HepRandomEn 173 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 176 const int s 174 const int size, double* vect ) 177 { 175 { 178 for( double* v = vect; v != vect + size; ++v 176 for( double* v = vect; v != vect + size; ++v ) 179 *v = shoot( anEngine, 1.0, 0.2 ); 177 *v = shoot( anEngine, 1.0, 0.2 ); 180 } 178 } 181 179 182 void RandBreitWigner::shootArray ( HepRandomEn 180 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 183 const int s 181 const int size, double* vect, 184 double a, d 182 double a, double b ) 185 { 183 { 186 for( double* v = vect; v != vect + size; ++v 184 for( double* v = vect; v != vect + size; ++v ) 187 *v = shoot( anEngine, a, b ); 185 *v = shoot( anEngine, a, b ); 188 } 186 } 189 187 190 void RandBreitWigner::shootArray ( HepRandomEn 188 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine, 191 const int s 189 const int size, double* vect, 192 double a, d 190 double a, double b, double c ) 193 { 191 { 194 for( double* v = vect; v != vect + size; ++v 192 for( double* v = vect; v != vect + size; ++v ) 195 *v = shoot( anEngine, a, b, c ); 193 *v = shoot( anEngine, a, b, c ); 196 } 194 } 197 195 198 //---------------- 196 //---------------- 199 197 200 double RandBreitWigner::fire() 198 double RandBreitWigner::fire() 201 { 199 { 202 return fire( defaultA, defaultB ); 200 return fire( defaultA, defaultB ); 203 } 201 } 204 202 205 double RandBreitWigner::fire(double mean, doub 203 double RandBreitWigner::fire(double mean, double gamma) 206 { 204 { 207 double rval, displ; 205 double rval, displ; 208 206 209 rval = 2.0*localEngine->flat()-1.0; 207 rval = 2.0*localEngine->flat()-1.0; 210 displ = 0.5*gamma*std::tan(rval*CLHEP::half 208 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi); 211 209 212 return mean + displ; 210 return mean + displ; 213 } 211 } 214 212 215 double RandBreitWigner::fire(double mean, doub 213 double RandBreitWigner::fire(double mean, double gamma, double cut) 216 { 214 { 217 double val, rval, displ; 215 double val, rval, displ; 218 216 219 if ( gamma == 0.0 ) return mean; 217 if ( gamma == 0.0 ) return mean; 220 val = std::atan(2.0*cut/gamma); 218 val = std::atan(2.0*cut/gamma); 221 rval = 2.0*localEngine->flat()-1.0; 219 rval = 2.0*localEngine->flat()-1.0; 222 displ = 0.5*gamma*std::tan(rval*val); 220 displ = 0.5*gamma*std::tan(rval*val); 223 221 224 return mean + displ; 222 return mean + displ; 225 } 223 } 226 224 227 double RandBreitWigner::fireM2() 225 double RandBreitWigner::fireM2() 228 { 226 { 229 return fireM2( defaultA, defaultB ); 227 return fireM2( defaultA, defaultB ); 230 } 228 } 231 229 232 double RandBreitWigner::fireM2(double mean, do 230 double RandBreitWigner::fireM2(double mean, double gamma ) 233 { 231 { 234 double val, rval, displ; 232 double val, rval, displ; 235 233 236 if ( gamma == 0.0 ) return mean; 234 if ( gamma == 0.0 ) return mean; 237 val = std::atan(-mean/gamma); 235 val = std::atan(-mean/gamma); 238 rval = RandFlat::shoot(localEngine.get(),va 236 rval = RandFlat::shoot(localEngine.get(),val, CLHEP::halfpi); 239 displ = gamma*std::tan(rval); 237 displ = gamma*std::tan(rval); 240 238 241 return std::sqrt(mean*mean + mean*displ); 239 return std::sqrt(mean*mean + mean*displ); 242 } 240 } 243 241 244 double RandBreitWigner::fireM2(double mean, do 242 double RandBreitWigner::fireM2(double mean, double gamma, double cut ) 245 { 243 { 246 double rval, displ; 244 double rval, displ; 247 double lower, upper, tmp; 245 double lower, upper, tmp; 248 246 249 if ( gamma == 0.0 ) return mean; 247 if ( gamma == 0.0 ) return mean; 250 tmp = std::max(0.0,(mean-cut)); 248 tmp = std::max(0.0,(mean-cut)); 251 lower = std::atan( (tmp*tmp-mean*mean)/(mea 249 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) ); 252 upper = std::atan( ((mean+cut)*(mean+cut)-m 250 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) ); 253 rval = RandFlat::shoot(localEngine.get(),lo 251 rval = RandFlat::shoot(localEngine.get(),lower, upper); 254 displ = gamma*std::tan(rval); 252 displ = gamma*std::tan(rval); 255 253 256 return std::sqrt(std::max(0.0, mean*mean + 254 return std::sqrt(std::max(0.0, mean*mean + mean*displ)); 257 } 255 } 258 256 259 void RandBreitWigner::fireArray ( const int si 257 void RandBreitWigner::fireArray ( const int size, double* vect) 260 { 258 { 261 for( double* v = vect; v != vect + size; ++v 259 for( double* v = vect; v != vect + size; ++v ) 262 *v = fire(defaultA, defaultB ); 260 *v = fire(defaultA, defaultB ); 263 } 261 } 264 262 265 void RandBreitWigner::fireArray ( const int si 263 void RandBreitWigner::fireArray ( const int size, double* vect, 266 double a, do 264 double a, double b ) 267 { 265 { 268 for( double* v = vect; v != vect + size; ++v 266 for( double* v = vect; v != vect + size; ++v ) 269 *v = fire( a, b ); 267 *v = fire( a, b ); 270 } 268 } 271 269 272 void RandBreitWigner::fireArray ( const int si 270 void RandBreitWigner::fireArray ( const int size, double* vect, 273 double a, do 271 double a, double b, double c ) 274 { 272 { 275 for( double* v = vect; v != vect + size; ++v 273 for( double* v = vect; v != vect + size; ++v ) 276 *v = fire( a, b, c ); 274 *v = fire( a, b, c ); 277 } 275 } 278 276 279 277 280 std::ostream & RandBreitWigner::put ( std::ost 278 std::ostream & RandBreitWigner::put ( std::ostream & os ) const { 281 long pr=os.precision(20); << 279 int pr=os.precision(20); 282 std::vector<unsigned long> t(2); 280 std::vector<unsigned long> t(2); 283 os << " " << name() << "\n"; 281 os << " " << name() << "\n"; 284 os << "Uvec" << "\n"; 282 os << "Uvec" << "\n"; 285 t = DoubConv::dto2longs(defaultA); 283 t = DoubConv::dto2longs(defaultA); 286 os << defaultA << " " << t[0] << " " << t[1] 284 os << defaultA << " " << t[0] << " " << t[1] << "\n"; 287 t = DoubConv::dto2longs(defaultB); 285 t = DoubConv::dto2longs(defaultB); 288 os << defaultB << " " << t[0] << " " << t[1] 286 os << defaultB << " " << t[0] << " " << t[1] << "\n"; 289 os.precision(pr); 287 os.precision(pr); 290 return os; 288 return os; 291 } 289 } 292 290 293 std::istream & RandBreitWigner::get ( std::ist 291 std::istream & RandBreitWigner::get ( std::istream & is ) { 294 std::string inName; 292 std::string inName; 295 is >> inName; 293 is >> inName; 296 if (inName != name()) { 294 if (inName != name()) { 297 is.clear(std::ios::badbit | is.rdstate()); 295 is.clear(std::ios::badbit | is.rdstate()); 298 std::cerr << "Mismatch when expecting to r 296 std::cerr << "Mismatch when expecting to read state of a " 299 << name() << " distribution\n" 297 << name() << " distribution\n" 300 << "Name found was " << inName 298 << "Name found was " << inName 301 << "\nistream is left in the badbit st 299 << "\nistream is left in the badbit state\n"; 302 return is; 300 return is; 303 } 301 } 304 if (possibleKeywordInput(is, "Uvec", default 302 if (possibleKeywordInput(is, "Uvec", defaultA)) { 305 std::vector<unsigned long> t(2); 303 std::vector<unsigned long> t(2); 306 is >> defaultA >> t[0] >> t[1]; defaultA = 304 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 307 is >> defaultB >> t[0] >> t[1]; defaultB = 305 is >> defaultB >> t[0] >> t[1]; defaultB = DoubConv::longs2double(t); 308 return is; 306 return is; 309 } 307 } 310 // is >> defaultA encompassed by possibleKey 308 // is >> defaultA encompassed by possibleKeywordInput 311 is >> defaultB; 309 is >> defaultB; 312 return is; 310 return is; 313 } 311 } 314 312 315 313 316 } // namespace CLHEP 314 } // namespace CLHEP 317 315 318 316