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 
a
@ a
Definition: testCategories.cc:125
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandBreitWigner::shootM2
static double shootM2(double a=1.0, double b=0.2)
Definition: RandBreitWigner.cc:71
b
@ b
Definition: testCategories.cc:125
CLHEP::RandBreitWigner::shoot
static double shoot(double a=1.0, double b=0.2)
Definition: RandBreitWigner.cc:49
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandBreitWigner::shootArray
static void shootArray(const int size, double *vect)
Definition: RandBreitWigner.cc:98
CLHEP::RandBreitWigner::fire
double fire()
Definition: RandBreitWigner.cc:199
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandBreitWigner::fireArray
void fireArray(const int size, double *vect)
Definition: RandBreitWigner.cc:258
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandBreitWigner::fireM2
double fireM2()
Definition: RandBreitWigner.cc:226
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandBreitWigner::~RandBreitWigner
virtual ~RandBreitWigner()
Definition: RandBreitWigner.cc:34
CLHEP::RandBreitWigner::engine
HepRandomEngine & engine()
Definition: RandBreitWigner.cc:32
CLHEP::RandBreitWigner::operator()
double operator()()
Definition: RandBreitWigner.cc:37
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandBreitWigner::put
std::ostream & put(std::ostream &os) const
Definition: RandBreitWigner.cc:279
CLHEP::RandBreitWigner::name
std::string name() const
Definition: RandBreitWigner.cc:31
CLHEP::RandBreitWigner::get
std::istream & get(std::istream &is)
Definition: RandBreitWigner.cc:299
CLHEP::HepRandomEngine::flat
virtual double flat()=0
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:60