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