22 #include "CLHEP/Random/RandBreitWigner.h"
23 #include "CLHEP/Units/PhysicalConstants.h"
24 #include "CLHEP/Random/DoubConv.h"
33 RandBreitWigner::~RandBreitWigner() {
36 double RandBreitWigner::operator()() {
37 return fire( defaultA, defaultB );
40 double RandBreitWigner::operator()(
double a,
double b ) {
44 double RandBreitWigner::operator()(
double a,
double b,
double c ) {
45 return fire( a, b, c );
52 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0;
60 double val, rval, displ;
62 if ( gamma == 0.0 )
return mean;
63 val = std::atan(2.0*cut/gamma);
64 rval = 2.0*HepRandom::getTheEngine()->flat()-1.0;
65 displ = 0.5*gamma*std::tan(rval*val);
70 double RandBreitWigner::shootM2(
double mean,
double gamma )
72 double val, rval, displ;
74 if ( gamma == 0.0 )
return mean;
75 val = std::atan(-mean/gamma);
77 displ = gamma*std::tan(rval);
79 return std::sqrt(mean*mean + mean*displ);
82 double RandBreitWigner::shootM2(
double mean,
double gamma,
double cut )
85 double lower, upper, tmp;
87 if ( gamma == 0.0 )
return mean;
89 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) );
90 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
92 displ = gamma*std::tan(rval);
94 return std::sqrt(
std::max(0.0, mean*mean + mean*displ));
97 void RandBreitWigner::shootArray (
const int size,
double* vect )
99 for(
double* v = vect; v != vect + size; ++v )
100 *v =
shoot( 1.0, 0.2 );
103 void RandBreitWigner::shootArray (
const int size,
double* vect,
106 for(
double* v = vect; v != vect + size; ++v )
110 void RandBreitWigner::shootArray (
const int size,
double* vect,
114 for(
double* v = vect; v != vect + size; ++v )
115 *v =
shoot( a, b, c );
121 double mean,
double gamma)
125 rval = 2.0*anEngine->flat()-1.0;
132 double mean,
double gamma,
double cut )
134 double val, rval, displ;
136 if ( gamma == 0.0 )
return mean;
137 val = std::atan(2.0*cut/gamma);
138 rval = 2.0*anEngine->flat()-1.0;
139 displ = 0.5*gamma*std::tan(rval*val);
144 double RandBreitWigner::shootM2(HepRandomEngine* anEngine,
145 double mean,
double gamma )
147 double val, rval, displ;
149 if ( gamma == 0.0 )
return mean;
150 val = std::atan(-mean/gamma);
152 displ = gamma*std::tan(rval);
154 return std::sqrt(mean*mean + mean*displ);
157 double RandBreitWigner::shootM2(HepRandomEngine* anEngine,
158 double mean,
double gamma,
double cut )
161 double lower, upper, tmp;
163 if ( gamma == 0.0 )
return mean;
165 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) );
166 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
168 displ = gamma*std::tan(rval);
170 return std::sqrt(
std::max(0.0, mean*mean+mean*displ) );
173 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
174 const int size,
double* vect )
176 for(
double* v = vect; v != vect + size; ++v )
177 *v =
shoot( anEngine, 1.0, 0.2 );
180 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
181 const int size,
double* vect,
184 for(
double* v = vect; v != vect + size; ++v )
185 *v =
shoot( anEngine, a, b );
188 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
189 const int size,
double* vect,
190 double a,
double b,
double c )
192 for(
double* v = vect; v != vect + size; ++v )
193 *v =
shoot( anEngine, a, b, c );
198 double RandBreitWigner::fire()
200 return fire( defaultA, defaultB );
203 double RandBreitWigner::fire(
double mean,
double gamma)
207 rval = 2.0*localEngine->flat()-1.0;
213 double RandBreitWigner::fire(
double mean,
double gamma,
double cut)
215 double val, rval, displ;
217 if ( gamma == 0.0 )
return mean;
218 val = std::atan(2.0*cut/gamma);
219 rval = 2.0*localEngine->flat()-1.0;
220 displ = 0.5*gamma*std::tan(rval*val);
225 double RandBreitWigner::fireM2()
227 return fireM2( defaultA, defaultB );
230 double RandBreitWigner::fireM2(
double mean,
double gamma )
232 double val, rval, displ;
234 if ( gamma == 0.0 )
return mean;
235 val = std::atan(-mean/gamma);
237 displ = gamma*std::tan(rval);
239 return std::sqrt(mean*mean + mean*displ);
242 double RandBreitWigner::fireM2(
double mean,
double gamma,
double cut )
245 double lower, upper, tmp;
247 if ( gamma == 0.0 )
return mean;
249 lower = std::atan( (tmp*tmp-mean*mean)/(mean*gamma) );
250 upper = std::atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
252 displ = gamma*std::tan(rval);
254 return std::sqrt(
std::max(0.0, mean*mean + mean*displ));
257 void RandBreitWigner::fireArray (
const int size,
double* vect)
259 for(
double* v = vect; v != vect + size; ++v )
260 *v = fire(defaultA, defaultB );
263 void RandBreitWigner::fireArray (
const int size,
double* vect,
266 for(
double* v = vect; v != vect + size; ++v )
270 void RandBreitWigner::fireArray (
const int size,
double* vect,
271 double a,
double b,
double c )
273 for(
double* v = vect; v != vect + size; ++v )
274 *v = fire( a, b, c );
278 std::ostream & RandBreitWigner::put ( std::ostream & os )
const {
279 int pr=os.precision(20);
280 std::vector<unsigned long> t(2);
281 os <<
" " <<
name() <<
"\n";
282 os <<
"Uvec" <<
"\n";
283 t = DoubConv::dto2longs(defaultA);
284 os << defaultA <<
" " << t[0] <<
" " << t[1] <<
"\n";
285 t = DoubConv::dto2longs(defaultB);
286 os << defaultB <<
" " << t[0] <<
" " << t[1] <<
"\n";
291 std::istream & RandBreitWigner::get ( std::istream & is ) {
294 if (inName !=
name()) {
295 is.clear(std::ios::badbit | is.rdstate());
296 std::cerr <<
"Mismatch when expecting to read state of a "
297 <<
name() <<
" distribution\n"
298 <<
"Name found was " << inName
299 <<
"\nistream is left in the badbit state\n";
302 if (possibleKeywordInput(is,
"Uvec", defaultA)) {
303 std::vector<unsigned long> t(2);
304 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
305 is >> defaultB >> t[0] >> t[1]; defaultB = DoubConv::longs2double(t);
ThreeVector shoot(const G4int Ap, const G4int Af)
static const double halfpi
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
T max(const T t1, const T t2)
brief Return the largest of the two arguments