CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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