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"
23 #include "CLHEP/Random/RandBreitWigner.h"
24 #include "CLHEP/Units/PhysicalConstants.h"
25 #include "CLHEP/Random/DoubConv.hh"
26 #include <algorithm> // for min() and max()
27 #include <cmath>
28 
29 namespace CLHEP {
30 
31 std::string RandBreitWigner::name() const {return "RandBreitWigner";}
32 HepRandomEngine & RandBreitWigner::engine() {return *localEngine;}
33 
35 }
36 
38  return fire( defaultA, defaultB );
39 }
40 
41 double RandBreitWigner::operator()( double a, double b ) {
42  return fire( a, b );
43 }
44 
45 double RandBreitWigner::operator()( double a, double b, double c ) {
46  return fire( a, b, c );
47 }
48 
49 double 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 
59 double 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 
71 double 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 
83 double 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 
98 void 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 
104 void 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 
111 void 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 
204 double 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 
214 double 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 
231 double 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 
243 double 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 
258 void RandBreitWigner::fireArray ( const int size, double* vect)
259 {
260  for( double* v = vect; v != vect + size; ++v )
261  *v = fire(defaultA, defaultB );
262 }
263 
264 void 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 
271 void 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 
279 std::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 
299 std::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 
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
static void shootArray(const int size, double *vect)
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
static double shoot(double a=1.0, double b=0.2)
std::string name() const
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
virtual double flat()=0
HepRandomEngine & engine()
void fireArray(const int size, double *vect)
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
std::ostream & put(std::ostream &os) const
static double shoot()
Definition: RandFlat.cc:60
static double shootM2(double a=1.0, double b=0.2)
std::istream & get(std::istream &is)