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