CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandBreitWigner.cc
Go to the documentation of this file.
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