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