Geant4  10.02.p01
RandBreitWigner.cc
Go to the documentation of this file.
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 
22 #include "CLHEP/Random/RandBreitWigner.h"
23 #include "CLHEP/Units/PhysicalConstants.h"
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";}
31 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;}
32 
33 RandBreitWigner::~RandBreitWigner() {
34 }
35 
36 double RandBreitWigner::operator()() {
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 
120 double RandBreitWigner::shoot(HepRandomEngine* anEngine,
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 
131 double RandBreitWigner::shoot(HepRandomEngine* anEngine,
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 
144 double RandBreitWigner::shootM2(HepRandomEngine* anEngine,
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 
157 double RandBreitWigner::shootM2(HepRandomEngine* anEngine,
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 
173 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
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 
180 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
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 
188 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
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 
198 double RandBreitWigner::fire()
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 
225 double RandBreitWigner::fireM2()
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";
283  t = DoubConv::dto2longs(defaultA);
284  os << defaultA << " " << t[0] << " " << t[1] << "\n";
285  t = DoubConv::dto2longs(defaultB);
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()) {
295  is.clear(std::ios::badbit | is.rdstate());
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 
ThreeVector shoot(const G4int Ap, const G4int Af)
static const double halfpi
Definition: G4SIunits.hh:76
G4String name
Definition: TRTMaterials.hh:40
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
Definition: Evaluator.cc:358
G4double a
Definition: TRTMaterials.hh:39
T max(const T t1, const T t2)
brief Return the largest of the two arguments