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