| 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