50 #include "RandBreitWigner.h"
52 #if !defined(ST_NO_NAMESPACES)
56 RandBreitWigner::~RandBreitWigner() {
57 if ( deleteEngine )
delete localEngine;
60 HepDouble RandBreitWigner::operator()() {
64 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma)
66 HepDouble rval, displ;
68 rval = 2.0*HepRandom::getTheGenerator()->flat()-1.0;
69 displ = 0.5*gamma*tan(rval*M_PI_2);
74 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma, HepDouble cut)
76 HepDouble val, rval, displ;
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);
86 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma )
88 HepDouble val, rval, displ;
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);
95 return ::sqrt(mean*mean + mean*displ);
98 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma, HepDouble cut )
100 HepDouble rval, displ;
101 HepDouble lower, upper, tmp;
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);
110 return ::sqrt(max(0.0, mean*mean + mean*displ));
114 void RandBreitWigner::shootArray (
const HepInt size, HepDouble* vect,
115 HepDouble a, HepDouble b,
120 for (i=0; i<size; ++i)
121 vect[i] = shoot( a, b, c );
125 #ifndef ST_NO_TEMPLATE_DEF_ARGS
126 RandBreitWigner::shootArray ( vector<HepDouble>& vec,
127 HepDouble a, HepDouble b,
130 RandBreitWigner::shootArray ( vector<HepDouble, allocator<HepDouble> >& vec,
131 HepDouble a, HepDouble b,
135 for (
unsigned int i=0; i<vec.size(); ++i)
136 vec[i] = shoot( a, b, c );
142 HepDouble mean, HepDouble gamma)
144 HepDouble rval, displ;
146 rval = 2.0*anEngine->flat()-1.0;
147 displ = 0.5*gamma*tan(rval*M_PI_2);
153 HepDouble mean, HepDouble gamma, HepDouble cut )
155 HepDouble val, rval, displ;
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);
166 HepDouble mean, HepDouble gamma )
168 HepDouble val, rval, displ;
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);
175 return ::sqrt(mean*mean + mean*displ);
179 HepDouble mean, HepDouble gamma, HepDouble cut )
181 HepDouble rval, displ;
182 HepDouble lower, upper, tmp;
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);
191 return ::sqrt( max(0.0, mean*mean+mean*displ) );
195 const HepInt size, HepDouble* vect,
196 HepDouble a, HepDouble b,
201 for (i=0; i<size; ++i)
202 vect[i] = shoot( anEngine, a, b, c );
206 #ifndef ST_NO_TEMPLATE_DEF_ARGS
208 vector<HepDouble>& vec,
209 HepDouble a, HepDouble b,
213 vector<HepDouble,allocator<HepDouble> >& vec,
214 HepDouble a, HepDouble b,
218 for (
unsigned int i=0; i<vec.size(); ++i)
219 vec[i] = shoot( anEngine, a, b, c );
223 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma)
225 HepDouble rval, displ;
227 rval = 2.0*localEngine->flat()-1.0;
228 displ = 0.5*gamma*tan(rval*M_PI_2);
233 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma, HepDouble cut)
235 HepDouble val, rval, displ;
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);
245 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma )
247 HepDouble val, rval, displ;
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);
254 return ::sqrt(mean*mean + mean*displ);
257 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma, HepDouble cut )
259 HepDouble rval, displ;
260 HepDouble lower, upper, tmp;
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);
269 return ::sqrt(max(0.0, mean*mean + mean*displ));
272 void RandBreitWigner::fireArray (
const HepInt size, HepDouble* vect,
273 HepDouble a, HepDouble b,
278 for (i=0; i<size; ++i)
279 vect[i] = fire( a, b, c );
283 #ifndef ST_NO_TEMPLATE_DEF_ARGS
284 RandBreitWigner::fireArray ( vector<HepDouble>& vec,
285 HepDouble a, HepDouble b,
288 RandBreitWigner::fireArray ( vector<HepDouble, allocator<HepDouble> >& vec,
289 HepDouble a, HepDouble b,
293 for (
unsigned int i=0; i<vec.size(); ++i)
294 vec[i] = fire( a, b, c );