| Geant4 Cross Reference |
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandBreitWigner ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 // This file is part of Geant4 (simulation toolkit for HEP).
9
10 // =======================================================================
11 // Gabriele Cosmo - Created: 5th September 1995
12 // - Added methods to shoot arrays: 28th July 1997
13 // J.Marraffino - Added default arguments as attributes and
14 // operator() with arguments: 16th Feb 1998
15 // M Fischler - put and get to/from streams 12/10/04
16 // M Fischler - put/get to/from streams uses pairs of ulongs when
17 // + storing doubles avoid problems with precision
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 {return "RandBreitWigner";}
33 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;}
34
35 RandBreitWigner::~RandBreitWigner() {
36 }
37
38 double RandBreitWigner::operator()() {
39 return fire( defaultA, defaultB );
40 }
41
42 double RandBreitWigner::operator()( double a, double b ) {
43 return fire( a, b );
44 }
45
46 double RandBreitWigner::operator()( double a, double b, double c ) {
47 return fire( a, b, c );
48 }
49
50 double RandBreitWigner::shoot(double mean, double gamma)
51 {
52 double rval, displ;
53
54 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0;
55 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi);
56
57 return mean + displ;
58 }
59
60 double RandBreitWigner::shoot(double mean, double gamma, double cut)
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()-1.0;
67 displ = 0.5*gamma*std::tan(rval*val);
68
69 return mean + displ;
70 }
71
72 double RandBreitWigner::shootM2(double mean, double gamma )
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, double gamma, double cut )
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)/(mean*gamma) );
92 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
93 rval = RandFlat::shoot(lower, upper);
94 displ = gamma*std::tan(rval);
95
96 return std::sqrt(std::max(0.0, mean*mean + mean*displ));
97 }
98
99 void RandBreitWigner::shootArray ( const int size, double* vect )
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 size, double* vect,
106 double a, double b )
107 {
108 for( double* v = vect; v != vect + size; ++v )
109 *v = shoot( a, b );
110 }
111
112 void RandBreitWigner::shootArray ( const int size, double* vect,
113 double a, double b,
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* anEngine,
123 double mean, double gamma)
124 {
125 double rval, displ;
126
127 rval = 2.0*anEngine->flat()-1.0;
128 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi);
129
130 return mean + displ;
131 }
132
133 double RandBreitWigner::shoot(HepRandomEngine* anEngine,
134 double mean, double gamma, double cut )
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(HepRandomEngine* anEngine,
147 double mean, double gamma )
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::halfpi);
154 displ = gamma*std::tan(rval);
155
156 return std::sqrt(mean*mean + mean*displ);
157 }
158
159 double RandBreitWigner::shootM2(HepRandomEngine* anEngine,
160 double mean, double gamma, double cut )
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)/(mean*gamma) );
168 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
169 rval = RandFlat::shoot(anEngine, lower, upper);
170 displ = gamma*std::tan(rval);
171
172 return std::sqrt( std::max(0.0, mean*mean+mean*displ) );
173 }
174
175 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
176 const int size, double* vect )
177 {
178 for( double* v = vect; v != vect + size; ++v )
179 *v = shoot( anEngine, 1.0, 0.2 );
180 }
181
182 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
183 const int size, double* vect,
184 double a, double b )
185 {
186 for( double* v = vect; v != vect + size; ++v )
187 *v = shoot( anEngine, a, b );
188 }
189
190 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
191 const int size, double* vect,
192 double a, double b, double c )
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, double gamma)
206 {
207 double rval, displ;
208
209 rval = 2.0*localEngine->flat()-1.0;
210 displ = 0.5*gamma*std::tan(rval*CLHEP::halfpi);
211
212 return mean + displ;
213 }
214
215 double RandBreitWigner::fire(double mean, double gamma, double cut)
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, double gamma )
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(),val, CLHEP::halfpi);
239 displ = gamma*std::tan(rval);
240
241 return std::sqrt(mean*mean + mean*displ);
242 }
243
244 double RandBreitWigner::fireM2(double mean, double gamma, double cut )
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)/(mean*gamma) );
252 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
253 rval = RandFlat::shoot(localEngine.get(),lower, upper);
254 displ = gamma*std::tan(rval);
255
256 return std::sqrt(std::max(0.0, mean*mean + mean*displ));
257 }
258
259 void RandBreitWigner::fireArray ( const int size, double* vect)
260 {
261 for( double* v = vect; v != vect + size; ++v )
262 *v = fire(defaultA, defaultB );
263 }
264
265 void RandBreitWigner::fireArray ( const int size, double* vect,
266 double a, double b )
267 {
268 for( double* v = vect; v != vect + size; ++v )
269 *v = fire( a, b );
270 }
271
272 void RandBreitWigner::fireArray ( const int size, double* vect,
273 double a, double b, double c )
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::ostream & os ) const {
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] << "\n";
287 t = DoubConv::dto2longs(defaultB);
288 os << defaultB << " " << t[0] << " " << t[1] << "\n";
289 os.precision(pr);
290 return os;
291 }
292
293 std::istream & RandBreitWigner::get ( std::istream & is ) {
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 read state of a "
299 << name() << " distribution\n"
300 << "Name found was " << inName
301 << "\nistream is left in the badbit state\n";
302 return is;
303 }
304 if (possibleKeywordInput(is, "Uvec", defaultA)) {
305 std::vector<unsigned long> t(2);
306 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
307 is >> defaultB >> t[0] >> t[1]; defaultB = DoubConv::longs2double(t);
308 return is;
309 }
310 // is >> defaultA encompassed by possibleKeywordInput
311 is >> defaultB;
312 return is;
313 }
314
315
316 } // namespace CLHEP
317
318