StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RandBreitWigner.cc
1 /***************************************************************************
2  *
3  * $Id: RandBreitWigner.cc,v 1.5 2016/01/22 17:10:49 smirnovd Exp $
4  *
5  * Author: Gabriele Cosmo - Created: 5th September 1995
6  * modified for SCL bl
7  ***************************************************************************
8  *
9  * Description:
10  * RandBreitWigner.cc,v 1.3 1997/08/12 00:38:37
11  * -----------------------------------------------------------------------
12  * HEP Random
13  * --- RandBreitWigner ---
14  * class implementation file
15  * -----------------------------------------------------------------------
16  * This file is part of Geant4 (simulation toolkit for HEP).
17  *
18  ***************************************************************************
19  *
20  * $Log: RandBreitWigner.cc,v $
21  * Revision 1.5 2016/01/22 17:10:49 smirnovd
22  * StarClassLibrary: Removed deprecated storage class specifier 'register'
23  *
24  * This keyword is deprecated since C++11 and serves no purpose
25  *
26  * "
27  * The register specifier is only allowed for objects declared at block scope and
28  * in function parameter lists. It indicates automatic storage duration, which is
29  * the default for these kinds of declarations. Additionally, the presence of this
30  * keyword may be used as a hint for the optimizer to store the value of this
31  * variable in a CPU register.
32  * "
33  *
34  * Revision 1.4 2003/09/02 17:59:34 perev
35  * gcc 3.2 updates + WarnOff
36  *
37  * Revision 1.3 1999/12/21 15:13:53 ullrich
38  * Modified to cope with new compiler version on Sun (CC5.0).
39  *
40  * Revision 1.2 1999/12/07 23:43:04 ullrich
41  * Modified to get rid of warnings on Linux.
42  *
43  * Revision 1.1 1999/01/30 03:58:59 fisyak
44  * Root Version of StarClassLibrary
45  *
46  * Revision 1.1 1999/01/23 00:29:03 ullrich
47  * Initial Revision
48  *
49  **************************************************************************/
50 #include "RandBreitWigner.h"
51 #include <algorithm> // for max(), tu
52 #if !defined(ST_NO_NAMESPACES)
53 using std::max;
54 #endif
55 
56 RandBreitWigner::~RandBreitWigner() {
57  if ( deleteEngine ) delete localEngine;
58 }
59 
60 HepDouble RandBreitWigner::operator()() {
61  return fire();
62 }
63 
64 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma)
65 {
66  HepDouble rval, displ;
67 
68  rval = 2.0*HepRandom::getTheGenerator()->flat()-1.0;
69  displ = 0.5*gamma*tan(rval*M_PI_2);
70 
71  return mean + displ;
72 }
73 
74 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma, HepDouble cut)
75 {
76  HepDouble val, rval, displ;
77 
78  if ( gamma == 0.0 ) return mean;
79  val = atan(2.0*cut/gamma);
80  rval = 2.0*HepRandom::getTheGenerator()->flat()-1.0;
81  displ = 0.5*gamma*tan(rval*val);
82 
83  return mean + displ;
84 }
85 
86 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma )
87 {
88  HepDouble val, rval, displ;
89 
90  if ( gamma == 0.0 ) return mean;
91  val = atan(-mean/gamma);
92  rval = RandFlat::shoot(val, M_PI_2);
93  displ = gamma*tan(rval);
94 
95  return ::sqrt(mean*mean + mean*displ);
96 }
97 
98 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma, HepDouble cut )
99 {
100  HepDouble rval, displ;
101  HepDouble lower, upper, tmp;
102 
103  if ( gamma == 0.0 ) return mean;
104  tmp = max(0.0,(mean-cut));
105  lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
106  upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
107  rval = RandFlat::shoot(lower, upper);
108  displ = gamma*tan(rval);
109 
110  return ::sqrt(max(0.0, mean*mean + mean*displ));
111 }
112 
113 
114 void RandBreitWigner::shootArray ( const HepInt size, HepDouble* vect,
115  HepDouble a, HepDouble b,
116  HepDouble c )
117 {
118  HepInt i;
119 
120  for (i=0; i<size; ++i)
121  vect[i] = shoot( a, b, c );
122 }
123 
124 void
125 #ifndef ST_NO_TEMPLATE_DEF_ARGS
126 RandBreitWigner::shootArray ( vector<HepDouble>& vec,
127  HepDouble a, HepDouble b,
128  HepDouble c )
129 #else
130 RandBreitWigner::shootArray ( vector<HepDouble, allocator<HepDouble> >& vec,
131  HepDouble a, HepDouble b,
132  HepDouble c )
133 #endif
134 {
135  for (unsigned int i=0; i<vec.size(); ++i)
136  vec[i] = shoot( a, b, c );
137 }
138 
139 //----------------
140 
141 HepDouble RandBreitWigner::shoot(HepRandomEngine* anEngine,
142  HepDouble mean, HepDouble gamma)
143 {
144  HepDouble rval, displ;
145 
146  rval = 2.0*anEngine->flat()-1.0;
147  displ = 0.5*gamma*tan(rval*M_PI_2);
148 
149  return mean + displ;
150 }
151 
152 HepDouble RandBreitWigner::shoot(HepRandomEngine* anEngine,
153  HepDouble mean, HepDouble gamma, HepDouble cut )
154 {
155  HepDouble val, rval, displ;
156 
157  if ( gamma == 0.0 ) return mean;
158  val = atan(2.0*cut/gamma);
159  rval = 2.0*anEngine->flat()-1.0;
160  displ = 0.5*gamma*tan(rval*val);
161 
162  return mean + displ;
163 }
164 
165 HepDouble RandBreitWigner::shootM2(HepRandomEngine* anEngine,
166  HepDouble mean, HepDouble gamma )
167 {
168  HepDouble val, rval, displ;
169 
170  if ( gamma == 0.0 ) return mean;
171  val = atan(-mean/gamma);
172  rval = RandFlat::shoot(anEngine,val, M_PI_2);
173  displ = gamma*tan(rval);
174 
175  return ::sqrt(mean*mean + mean*displ);
176 }
177 
178 HepDouble RandBreitWigner::shootM2(HepRandomEngine* anEngine,
179  HepDouble mean, HepDouble gamma, HepDouble cut )
180 {
181  HepDouble rval, displ;
182  HepDouble lower, upper, tmp;
183 
184  if ( gamma == 0.0 ) return mean;
185  tmp = max(0.0,(mean-cut));
186  lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
187  upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
188  rval = RandFlat::shoot(anEngine, lower, upper);
189  displ = gamma*tan(rval);
190 
191  return ::sqrt( max(0.0, mean*mean+mean*displ) );
192 }
193 
194 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
195  const HepInt size, HepDouble* vect,
196  HepDouble a, HepDouble b,
197  HepDouble c )
198 {
199  HepInt i;
200 
201  for (i=0; i<size; ++i)
202  vect[i] = shoot( anEngine, a, b, c );
203 }
204 
205 void
206 #ifndef ST_NO_TEMPLATE_DEF_ARGS
207 RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
208  vector<HepDouble>& vec,
209  HepDouble a, HepDouble b,
210  HepDouble c )
211 #else
212 RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
213  vector<HepDouble,allocator<HepDouble> >& vec,
214  HepDouble a, HepDouble b,
215  HepDouble c )
216 #endif
217 {
218  for (unsigned int i=0; i<vec.size(); ++i)
219  vec[i] = shoot( anEngine, a, b, c );
220 }
221 //----------------
222 
223 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma)
224 {
225  HepDouble rval, displ;
226 
227  rval = 2.0*localEngine->flat()-1.0;
228  displ = 0.5*gamma*tan(rval*M_PI_2);
229 
230  return mean + displ;
231 }
232 
233 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma, HepDouble cut)
234 {
235  HepDouble val, rval, displ;
236 
237  if ( gamma == 0.0 ) return mean;
238  val = atan(2.0*cut/gamma);
239  rval = 2.0*localEngine->flat()-1.0;
240  displ = 0.5*gamma*tan(rval*val);
241 
242  return mean + displ;
243 }
244 
245 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma )
246 {
247  HepDouble val, rval, displ;
248 
249  if ( gamma == 0.0 ) return mean;
250  val = atan(-mean/gamma);
251  rval = RandFlat::shoot(localEngine,val, M_PI_2);
252  displ = gamma*tan(rval);
253 
254  return ::sqrt(mean*mean + mean*displ);
255 }
256 
257 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma, HepDouble cut )
258 {
259  HepDouble rval, displ;
260  HepDouble lower, upper, tmp;
261 
262  if ( gamma == 0.0 ) return mean;
263  tmp = max(0.0,(mean-cut));
264  lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
265  upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
266  rval = RandFlat::shoot(localEngine,lower, upper);
267  displ = gamma*tan(rval);
268 
269  return ::sqrt(max(0.0, mean*mean + mean*displ));
270 }
271 
272 void RandBreitWigner::fireArray ( const HepInt size, HepDouble* vect,
273  HepDouble a, HepDouble b,
274  HepDouble c )
275 {
276  HepInt i;
277 
278  for (i=0; i<size; ++i)
279  vect[i] = fire( a, b, c );
280 }
281 
282 void
283 #ifndef ST_NO_TEMPLATE_DEF_ARGS
284 RandBreitWigner::fireArray ( vector<HepDouble>& vec,
285  HepDouble a, HepDouble b,
286  HepDouble c )
287 #else
288 RandBreitWigner::fireArray ( vector<HepDouble, allocator<HepDouble> >& vec,
289  HepDouble a, HepDouble b,
290  HepDouble c )
291 #endif
292 {
293  for (unsigned int i=0; i<vec.size(); ++i)
294  vec[i] = fire( a, b, c );
295 }
Generated by   doxygen 1.8.5