CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandSkewNormal.cc,v 1.1 2011/05/27 20:36:57 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandSkewNormal --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // M Fischler and L Garren - Created: 26 May 2011 00012 // ======================================================================= 00013 00014 #include "CLHEP/Random/defs.h" 00015 #include "CLHEP/Random/RandSkewNormal.h" 00016 #include "CLHEP/Random/RandGaussT.h" 00017 #include "CLHEP/Random/DoubConv.hh" 00018 00019 namespace CLHEP { 00020 00021 std::string RandSkewNormal::name() const {return "RandSkewNormal";} 00022 HepRandomEngine & RandSkewNormal::engine() {return *localEngine;} 00023 00024 RandSkewNormal::~RandSkewNormal() {} 00025 00026 double RandSkewNormal::operator()() 00027 { 00028 return fire( shapeParameter ); 00029 } 00030 00031 double RandSkewNormal::operator()( double shape ) 00032 { 00033 return fire( shape ); 00034 } 00035 00036 //------------- 00037 00038 double RandSkewNormal::shoot(HepRandomEngine* anEngine) 00039 { 00040 // really dumb use of RandSkewNormal 00041 double k = 1; 00042 return gaussianSkewNormal( anEngine, k ); 00043 } 00044 00045 double RandSkewNormal::shoot(HepRandomEngine* anEngine, double shape) 00046 { 00047 return gaussianSkewNormal( anEngine, shape ); 00048 } 00049 00050 double RandSkewNormal::shoot() 00051 { 00052 // really dumb use of RandSkewNormal 00053 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00054 double k = 1; 00055 return gaussianSkewNormal( anEngine, k ); 00056 } 00057 00058 double RandSkewNormal::shoot(double shape) 00059 { 00060 HepRandomEngine* anEngine = HepRandom::getTheEngine(); 00061 return gaussianSkewNormal( anEngine, shape ); 00062 } 00063 00064 void RandSkewNormal::shootArray( const int size, double* vect, 00065 double shape ) 00066 { 00067 for( double* v = vect; v != vect+size; ++v ) 00068 *v = shoot(shape); 00069 } 00070 00071 void RandSkewNormal::shootArray(HepRandomEngine* anEngine, const int size, 00072 double* vect, double shape ) 00073 { 00074 for( double* v = vect; v != vect+size; ++v ) 00075 *v = shoot(anEngine, shape); 00076 } 00077 00078 //------------- 00079 00080 double RandSkewNormal::fire() { 00081 return gaussianSkewNormal( getLocalEngine(), getShapeParameter() ); 00082 } 00083 00084 double RandSkewNormal::fire(double shape) { 00085 return gaussianSkewNormal( getLocalEngine(), shape ); 00086 } 00087 00088 void RandSkewNormal::fireArray( const int size, double* vect) 00089 { 00090 for( double* v = vect; v != vect+size; ++v ) 00091 *v = fire( shapeParameter ); 00092 } 00093 00094 void RandSkewNormal::fireArray( const int size, double* vect, 00095 double shape ) 00096 { 00097 for( double* v = vect; v != vect+size; ++v ) 00098 *v = fire( shape ); 00099 } 00100 00101 //------------- 00102 00103 double RandSkewNormal::gaussianSkewNormal ( HepRandomEngine* e, double k) 00104 { 00105 // RandSkewNormal is an implementation of Azzalini's SN generator 00106 // http://azzalini.stat.unipd.it/SN/ 00107 // K. McFarlane, June 2010. 00108 // For location parameter m = 0., scale = 1. 00109 // To get a distribution with scale parameter b and location m: 00110 // r = m + b * RandSkewNormal.fire(k); 00111 double u[2] = {0.}; 00112 RandGaussT::shootArray(e, 2, u, 0, 1); 00113 double delta = k/std::sqrt(1. + k*k); 00114 double u1 = delta*u[0] + std::sqrt(1 - delta*delta)*u[1]; 00115 double r = u[0] >= 0 ? u1 : -u1; 00116 return r; 00117 } 00118 00119 //------------- 00120 00121 std::ostream & RandSkewNormal::put ( std::ostream & os ) const { 00122 int pr=os.precision(20); 00123 std::vector<unsigned long> t(2); 00124 os << " " << name() << "\n"; 00125 os << "Uvec" << "\n"; 00126 t = DoubConv::dto2longs(shapeParameter); 00127 os << shapeParameter << " " << t[0] << " " << t[1] << "\n"; 00128 os.precision(pr); 00129 return os; 00130 } 00131 00132 std::istream & RandSkewNormal::get ( std::istream & is ) { 00133 std::string inName; 00134 is >> inName; 00135 if (inName != name()) { 00136 is.clear(std::ios::badbit | is.rdstate()); 00137 std::cerr << "Mismatch when expecting to read state of a " 00138 << name() << " distribution\n" 00139 << "Name found was " << inName 00140 << "\nistream is left in the badbit state\n"; 00141 return is; 00142 } 00143 if (possibleKeywordInput(is, "Uvec", shapeParameter)) { 00144 std::vector<unsigned long> t(2); 00145 is >> shapeParameter >> t[0] >> t[1]; shapeParameter = DoubConv::longs2double(t); 00146 return is; 00147 } 00148 // is >> shapeParameter encompassed by possibleKeywordInput 00149 return is; 00150 } 00151 00152 00153 } // namespace CLHEP