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