StRoot: StarClassLibrary/RandBreitWigner.cc Source File

00001 /***************************************************************************
00002  *
00003  * $Id: RandBreitWigner.cc,v 1.4 2003/09/02 17:59:34 perev Exp $
00004  *
00005  * Author: Gabriele Cosmo - Created: 5th September 1995
00006  *         modified for SCL bl
00007  ***************************************************************************
00008  *
00009  * Description:
00010  *               RandBreitWigner.cc,v 1.3 1997/08/12 00:38:37
00011  * -----------------------------------------------------------------------
00012  *                             HEP Random
00013  *                       --- RandBreitWigner ---
00014  *                      class implementation file
00015  * -----------------------------------------------------------------------
00016  * This file is part of Geant4 (simulation toolkit for HEP).
00017  *
00018  ***************************************************************************
00019  *
00020  * $Log: RandBreitWigner.cc,v $
00021  * Revision 1.4  2003/09/02 17:59:34  perev
00022  * gcc 3.2 updates + WarnOff
00023  *
00024  * Revision 1.3  1999/12/21 15:13:53  ullrich
00025  * Modified to cope with new compiler version on Sun (CC5.0).
00026  *
00027  * Revision 1.2  1999/12/07 23:43:04  ullrich
00028  * Modified to get rid of warnings on Linux.
00029  *
00030  * Revision 1.1  1999/01/30 03:58:59  fisyak
00031  * Root Version of StarClassLibrary
00032  *
00033  * Revision 1.1  1999/01/23 00:29:03  ullrich
00034  * Initial Revision
00035  *
00036  **************************************************************************/
00037 #include "RandBreitWigner.h"
00038 #include <algorithm>        // for max(), tu
00039 #if !defined(ST_NO_NAMESPACES)
00040 using std::max;
00041 #endif
00042 
00043 RandBreitWigner::~RandBreitWigner() {
00044   if ( deleteEngine ) delete localEngine;
00045 }
00046 
00047 HepDouble RandBreitWigner::operator()() {
00048    return fire();
00049 }
00050 
00051 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma)
00052 {
00053    HepDouble rval, displ;
00054 
00055    rval = 2.0*HepRandom::getTheGenerator()->flat()-1.0;
00056    displ = 0.5*gamma*tan(rval*M_PI_2);
00057 
00058    return mean + displ;
00059 }
00060 
00061 HepDouble RandBreitWigner::shoot(HepDouble mean, HepDouble gamma, HepDouble cut)
00062 {
00063    HepDouble val, rval, displ;
00064 
00065    if ( gamma == 0.0 ) return mean;
00066    val = atan(2.0*cut/gamma);
00067    rval = 2.0*HepRandom::getTheGenerator()->flat()-1.0;
00068    displ = 0.5*gamma*tan(rval*val);
00069 
00070    return mean + displ;
00071 }
00072 
00073 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma )
00074 {
00075    HepDouble val, rval, displ;
00076 
00077    if ( gamma == 0.0 ) return mean;
00078    val = atan(-mean/gamma);
00079    rval = RandFlat::shoot(val, M_PI_2);
00080    displ = gamma*tan(rval);
00081 
00082    return ::sqrt(mean*mean + mean*displ);
00083 }
00084 
00085 HepDouble RandBreitWigner::shootM2(HepDouble mean, HepDouble gamma, HepDouble cut )
00086 {
00087    HepDouble rval, displ;
00088    HepDouble lower, upper, tmp;
00089 
00090    if ( gamma == 0.0 ) return mean;
00091    tmp = max(0.0,(mean-cut));
00092    lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
00093    upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
00094    rval = RandFlat::shoot(lower, upper);
00095    displ = gamma*tan(rval);
00096 
00097    return ::sqrt(max(0.0, mean*mean + mean*displ));
00098 }
00099 
00100 
00101 void RandBreitWigner::shootArray ( const HepInt size, HepDouble* vect,
00102                                    HepDouble a, HepDouble b,
00103                                    HepDouble c )
00104 {
00105    register HepInt i;
00106 
00107    for (i=0; i<size; ++i)
00108      vect[i] = shoot( a, b, c );
00109 }
00110 
00111 void
00112 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00113 RandBreitWigner::shootArray ( vector<HepDouble>& vec,
00114                                    HepDouble a, HepDouble b,
00115                                    HepDouble c )
00116 #else
00117 RandBreitWigner::shootArray ( vector<HepDouble, allocator<HepDouble> >& vec,
00118                                    HepDouble a, HepDouble b,
00119                                    HepDouble c )
00120 #endif
00121 {
00122    for (unsigned int i=0; i<vec.size(); ++i)
00123      vec[i] = shoot( a, b, c );
00124 }
00125 
00126 //----------------
00127 
00128 HepDouble RandBreitWigner::shoot(HepRandomEngine* anEngine,
00129                                  HepDouble mean, HepDouble gamma)
00130 {
00131    HepDouble rval, displ;
00132 
00133    rval = 2.0*anEngine->flat()-1.0;
00134    displ = 0.5*gamma*tan(rval*M_PI_2);
00135 
00136    return mean + displ;
00137 }
00138 
00139 HepDouble RandBreitWigner::shoot(HepRandomEngine* anEngine,
00140                                  HepDouble mean, HepDouble gamma, HepDouble cut )
00141 {
00142    HepDouble val, rval, displ;
00143 
00144    if ( gamma == 0.0 ) return mean;
00145    val = atan(2.0*cut/gamma);
00146    rval = 2.0*anEngine->flat()-1.0;
00147    displ = 0.5*gamma*tan(rval*val);
00148 
00149    return mean + displ;
00150 }
00151 
00152 HepDouble RandBreitWigner::shootM2(HepRandomEngine* anEngine,
00153                                    HepDouble mean, HepDouble gamma )
00154 {
00155    HepDouble val, rval, displ;
00156 
00157    if ( gamma == 0.0 ) return mean;
00158    val = atan(-mean/gamma);
00159    rval = RandFlat::shoot(anEngine,val, M_PI_2);
00160    displ = gamma*tan(rval);
00161 
00162    return ::sqrt(mean*mean + mean*displ);
00163 }
00164 
00165 HepDouble RandBreitWigner::shootM2(HepRandomEngine* anEngine,
00166                                    HepDouble mean, HepDouble gamma, HepDouble cut )
00167 {
00168    HepDouble rval, displ;
00169    HepDouble lower, upper, tmp;
00170 
00171    if ( gamma == 0.0 ) return mean;
00172    tmp = max(0.0,(mean-cut));
00173    lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
00174    upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
00175    rval = RandFlat::shoot(anEngine, lower, upper);
00176    displ = gamma*tan(rval);
00177 
00178    return ::sqrt( max(0.0, mean*mean+mean*displ) );
00179 }
00180 
00181 void RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
00182                                    const HepInt size, HepDouble* vect,
00183                                    HepDouble a, HepDouble b,
00184                                    HepDouble c )
00185 {
00186    register HepInt i;
00187 
00188    for (i=0; i<size; ++i)
00189      vect[i] = shoot( anEngine, a, b, c );
00190 }
00191 
00192 void
00193 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00194 RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
00195                                    vector<HepDouble>& vec,
00196                                    HepDouble a, HepDouble b,
00197                                    HepDouble c )
00198 #else
00199 RandBreitWigner::shootArray ( HepRandomEngine* anEngine,
00200                                    vector<HepDouble,allocator<HepDouble> >& vec,
00201                                    HepDouble a, HepDouble b,
00202                                    HepDouble c )
00203 #endif
00204 {
00205    for (unsigned int i=0; i<vec.size(); ++i)
00206      vec[i] = shoot( anEngine, a, b, c );
00207 }
00208 //----------------
00209 
00210 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma)
00211 {
00212    HepDouble rval, displ;
00213 
00214    rval = 2.0*localEngine->flat()-1.0;
00215    displ = 0.5*gamma*tan(rval*M_PI_2);
00216 
00217    return mean + displ;
00218 }
00219 
00220 HepDouble RandBreitWigner::fire(HepDouble mean, HepDouble gamma, HepDouble cut)
00221 {
00222    HepDouble val, rval, displ;
00223 
00224    if ( gamma == 0.0 ) return mean;
00225    val = atan(2.0*cut/gamma);
00226    rval = 2.0*localEngine->flat()-1.0;
00227    displ = 0.5*gamma*tan(rval*val);
00228 
00229    return mean + displ;
00230 }
00231 
00232 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma )
00233 {
00234    HepDouble val, rval, displ;
00235 
00236    if ( gamma == 0.0 ) return mean;
00237    val = atan(-mean/gamma);
00238    rval = RandFlat::shoot(localEngine,val, M_PI_2);
00239    displ = gamma*tan(rval);
00240 
00241    return ::sqrt(mean*mean + mean*displ);
00242 }
00243 
00244 HepDouble RandBreitWigner::fireM2(HepDouble mean, HepDouble gamma, HepDouble cut )
00245 {
00246    HepDouble rval, displ;
00247    HepDouble lower, upper, tmp;
00248 
00249    if ( gamma == 0.0 ) return mean;
00250    tmp = max(0.0,(mean-cut));
00251    lower = atan( (tmp*tmp-mean*mean)/(mean*gamma) );
00252    upper = atan( ((mean+cut)*(mean+cut)-mean*mean)/(mean*gamma) );
00253    rval = RandFlat::shoot(localEngine,lower, upper);
00254    displ = gamma*tan(rval);
00255 
00256    return ::sqrt(max(0.0, mean*mean + mean*displ));
00257 }
00258 
00259 void RandBreitWigner::fireArray ( const HepInt size, HepDouble* vect,
00260                                   HepDouble a, HepDouble b,
00261                                   HepDouble c )
00262 {
00263    register HepInt i;
00264 
00265    for (i=0; i<size; ++i)
00266      vect[i] = fire( a, b, c );
00267 }
00268 
00269 void
00270 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00271 RandBreitWigner::fireArray ( vector<HepDouble>& vec,
00272                                   HepDouble a, HepDouble b,
00273                                   HepDouble c )
00274 #else
00275 RandBreitWigner::fireArray ( vector<HepDouble, allocator<HepDouble> >& vec,
00276                              HepDouble a, HepDouble b,
00277                              HepDouble c )
00278 #endif
00279 {
00280    for (unsigned int i=0; i<vec.size(); ++i)
00281      vec[i] = fire( a, b, c );
00282 }