/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Radiation.cc
Go to the documentation of this file.
00001 #include "GenDecay/Radiation.h"
00002 
00003 #include "GaudiKernel/ServiceHandle.h"
00004 #include "GaudiKernel/IRndmGenSvc.h"
00005 #include "GaudiKernel/IMessageSvc.h"
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/GaudiException.h"
00008 
00009 #include "CLHEP/Units/PhysicalConstants.h"
00010 #include "CLHEP/Units/SystemOfUnits.h"
00011 
00012 #include <cmath>
00013 #include <sstream>
00014 
00015 using namespace std;
00016 using namespace GenDecay;
00017 
00018 std::ostream& operator<<(std::ostream& o, const Radiation& r) 
00019 {
00020     return o << r.asString(); 
00021 }
00022 
00023 // base
00024 
00025 Radiation::Radiation(double energy) : m_energy(energy) 
00026 {
00027 }
00028 
00029 Radiation::~Radiation() 
00030 {
00031 }
00032 
00033 double Radiation::kineticEnergy() const
00034 {
00035     return m_energy; 
00036 }
00037 
00038 // alpha
00039 
00040 AlphaRadiation::AlphaRadiation(double energy, int parentA)
00041     : Radiation(energy), m_parentA(parentA)
00042 {
00043 }
00044 AlphaRadiation::~AlphaRadiation()
00045 {
00046 }
00047 
00048 std::string AlphaRadiation::asString() const
00049 {
00050     stringstream ss;
00051     ss << "alpha: A=" << m_parentA << ", KE=" << m_energy << ends;
00052     return ss.str().c_str();
00053 }
00054 
00055 double AlphaRadiation::kineticEnergy() const
00056 {
00057     // NNDC tables give alpha kinetic energy, not Qvalue
00058     //return (m_parentA-4.0)/m_parentA * m_energy;
00059     return m_energy;
00060 }
00061 
00062 int AlphaRadiation::pid() const
00063 {
00064     return 1000020040;
00065 }
00066 double AlphaRadiation::mass() const
00067 {
00068     return 3.727000*CLHEP::GeV;
00069 }
00070 
00071 // beta
00072 
00073 BetaRadiation::BetaRadiation(double energy, int parentZ)
00074     : Radiation(energy)
00075     , m_parentZ(parentZ)
00076     , m_daughterZ(0)
00077     , m_betaSign(0)
00078     , m_endpoint(0.0)
00079 {
00080     if (parentZ < 0) {          // beta+ decay 
00081         m_parentZ = -parentZ;
00082         m_daughterZ = m_parentZ - 1;
00083         m_betaSign = +1;
00084         m_endpoint = m_energy - 2.0*CLHEP::electron_mass_c2;
00085     }
00086     else {                      // # beta- decay
00087         m_parentZ = parentZ;
00088         m_daughterZ = m_parentZ + 1;
00089         m_betaSign = -1;
00090         m_endpoint = m_energy;
00091     }
00092 
00093     ServiceHandle<IRndmGenSvc> rgsh("RndmGenSvc","BetaRadiation");
00094     IRndmGenSvc *rgs = rgsh.operator->();
00095     if (m_rand.initialize(rgs, Rndm::Flat(0,1)).isFailure()) {
00096         throw GaudiException("Failed to initialize uniform random numbers",
00097                              "GenDecay::Radiation",StatusCode::FAILURE);
00098     }
00099     if (!m_rand) {
00100         throw GaudiException("Got null Rndm::Numbers object",
00101                              "GenDecay::Radiation",StatusCode::FAILURE);
00102     }
00103 
00104     this->normalize();
00105 }
00106 
00107 void BetaRadiation::normalize()
00108 {
00109     const int steps = 1000;
00110     double dx = m_endpoint/steps;
00111     double lo = dx/2.0;
00112 
00113     m_norm = 1.0;
00114     double sum = 0.0;
00115     for (int ind=0; ind<steps; ++ind) {
00116         sum += dx * this->dnde(ind*dx+lo);
00117     }
00118     m_norm = sum;
00119 
00120     double max = 0;
00121     for (int ind=0; ind<steps; ++ind) {
00122         double tmp = this->dnde(ind*dx+lo);
00123         if (tmp > max) max = tmp;
00124     }
00125     m_maximum = max;
00126 }
00127 
00128 BetaRadiation::~BetaRadiation()
00129 {
00130 }
00131 
00132 std::string BetaRadiation::asString() const
00133 {
00134     stringstream ss;
00135     ss << "beta"<< (m_betaSign < 0 ? '-' : '+') <<": Z=" 
00136        << m_parentZ << ", E_endpoint=" << m_energy << ends;
00137     return ss.str().c_str();
00138 }
00139 
00140 double BetaRadiation::kineticEnergy() const
00141 {
00142     // MC/rejection method to sample dN/dE spectrum
00143     while (true) {
00144         double T = m_rand() * m_endpoint;
00145         double P = m_rand() * m_maximum;
00146         if (P <= this->dnde(T)) return T;
00147     }
00148     return 0.0;
00149 }
00150 
00151 double BetaRadiation::dnde(double kineticEnergy) const
00152 {
00153     if (kineticEnergy > m_endpoint) return 0.0;
00154     if (kineticEnergy < 0.0) return 0.0;
00155     return this->fermi_function(kineticEnergy) * this->dnde_noff(kineticEnergy) / m_norm;
00156 }
00157 
00158 double BetaRadiation::dnde_noff(double kineticEnergy) const
00159 {
00160     double W = m_endpoint / CLHEP::electron_mass_c2 + 1.0;
00161     double E = kineticEnergy / CLHEP::electron_mass_c2 + 1.0;
00162     return sqrt(E*E-1.0) * (W-E)*(W-E) * E;
00163 }
00164 
00165 double BetaRadiation::fermi_function(double kineticEnergy) const
00166 {
00167     double E = kineticEnergy / CLHEP::electron_mass_c2 + 1.0;
00168     double P = sqrt(E*E-1.0);
00169     double U = -1*m_betaSign*m_daughterZ/137.0;
00170     double S = sqrt(1.0 - U*U) - 1.0;
00171     double Y = 2.0*M_PI*U*E/P;
00172     double A1 = U*U * E*E + P*P/4.0;
00173     double A2 = fabs(Y/(1.0-exp(-Y)));
00174     return pow(A1,S)*A2;
00175 }
00176 
00177 int BetaRadiation::pid() const
00178 {
00179     return -1*m_betaSign*11;
00180 }
00181 double BetaRadiation::mass() const
00182 {
00183     return CLHEP::electron_mass_c2;
00184 }
00185 
00186 GammaRadiation::GammaRadiation(double energy)
00187     : Radiation(energy)
00188 {
00189 }
00190 GammaRadiation::~GammaRadiation()
00191 {
00192 }
00193 
00194 std::string GammaRadiation::asString() const
00195 {
00196     stringstream ss;
00197     ss << "gamma: Energy=" << m_energy << ends;
00198     return ss.str().c_str();
00199 }
00200 
00201 int GammaRadiation::pid() const
00202 {
00203     return 22;
00204 }
00205 
00206 double GammaRadiation::mass() const
00207 {
00208     return 0;
00209 }
00210 
00211 ElectronCapture::ElectronCapture(double characteristic_energy)
00212     : Radiation(characteristic_energy)
00213 {
00214 }
00215 ElectronCapture::~ElectronCapture()
00216 {
00217 }
00218 
00219 std::string ElectronCapture::asString() const
00220 {
00221     stringstream ss;
00222     ss << "electroncapture: Energy=" << m_energy << ends;
00223     return ss.str().c_str();
00224 }
00225 
00226 int ElectronCapture::pid() const
00227 {
00228     return 0;                   // What to return?
00229 }
00230 
00231 double ElectronCapture::mass() const
00232 {
00233     return 0.0;
00234 }
00235 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:21:06 for GenDecay by doxygen 1.7.4