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