G4BEChargedChannel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // Implementation of the HETC88 code into Geant4.
00027 // Evaporation and De-excitation parts
00028 // T. Lampen, Helsinki Institute of Physics, May-2000
00029 //
00030 // 20120608  M. Kelsey -- Change vars 's','m','m2' to avoid name collisions
00031 
00032 #include "G4BEChargedChannel.hh"
00033 #include "G4SystemOfUnits.hh"
00034 
00035 G4BEChargedChannel::G4BEChargedChannel()
00036 {
00037     verboseLevel = 0;
00038 }
00039 
00040 
00041 G4BEChargedChannel::~G4BEChargedChannel()
00042 {
00043 }
00044 
00045 
00046 void G4BEChargedChannel::calculateProbability()
00047 {
00048   G4int residualZ = nucleusZ - particleZ;
00049   G4int residualA = nucleusA - particleA;
00050 
00051 //    Check if nucleus is too small, if this evaporation channel
00052 //    leads to an impossible residual nucleus or if there is no enough
00053 //    energy.
00054   if ( nucleusA < 2.0 * particleA || 
00055        nucleusZ < 2.0 * particleZ ||
00056        residualA <= residualZ || 
00057        excitationEnergy - getThresh() - correction < 0 )
00058     {
00059       if ( verboseLevel >= 6 )
00060         G4cout << "G4BEChargedChannel : calculateProbability for " << getName() << " = 0 " << G4endl;
00061       emissionProbability = 0;
00062       return;
00063     }
00064 
00065   // In HETC88 s-s0 was used in std::exp( s ), in which s0 was either 50 or
00066   // max(s_i), where i goes over all channels.
00067 
00068   G4double levelParam = getLevelDensityParameter();
00069   G4double slevel = 2 * std::sqrt( levelParam  * ( excitationEnergy - getThresh() - correction ) );
00070   G4double constant = A / 2 * ( 2 * spin + 1 ) * ( 1 + coulombFactor() );
00071   G4double eye1     = ( slevel*slevel - 3 * slevel + 3 ) / ( 4 * levelParam*levelParam ) * std::exp( slevel );
00072 
00073   emissionProbability = constant * std::pow( G4double(residualA), 0.6666666 ) * eye1;
00074 
00075   if ( verboseLevel >= 6 )
00076     G4cout << "G4BEChargedChannel : calculateProbability for " << getName() << G4endl
00077            << "                    res A = " << residualA << G4endl 
00078            << "                    res Z = " << residualZ << G4endl 
00079            << "                 c*(c_i+1) = "<< constant << G4endl
00080            << "                  qmfactor = "<< qmFactor() << G4endl
00081            << "             coulombfactor = "<< coulombFactor() << G4endl
00082            << "                        E = " << excitationEnergy << G4endl
00083            << "               correction = " << correction << G4endl
00084            << "                     eye1 = " << eye1 << G4endl
00085            << "               levelParam = " << levelParam << G4endl
00086            << "                   thresh = " << getThresh() << G4endl
00087            << "                        s = " << s << G4endl
00088            << "              probability = " << emissionProbability << G4endl;
00089 
00090   return;
00091 }
00092 
00093 
00094 G4double G4BEChargedChannel::sampleKineticEnergy()
00095 {
00096   G4double levelParam;
00097   levelParam = getLevelDensityParameter();
00098   
00099   const G4double xMax  = excitationEnergy - getThresh() - correction; // maximum number
00100   const G4double xProb = ( - 1 + std::sqrt ( 1 + 4 * levelParam * xMax ) ) / ( 2 * levelParam ); // most probable value
00101   const G4double maxProb = xProb * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - xProb ) ) ); // maximum value of P(x)
00102 
00103   // Sample x according to density function P(x) with rejection method
00104   G4double r1;
00105   G4double r2;
00106   G4int koe=0;
00107   do
00108     {
00109       r1 = G4UniformRand() * xMax;
00110       r2 = G4UniformRand() * maxProb;
00111       koe++;
00112     }
00113   while (  r1 * std::exp ( 2 * std::sqrt ( levelParam * ( xMax - r1 ) ) )  < r2 );
00114 
00115 //  G4cout << "Q ch " << koe << G4endl;
00116   G4double kineticEnergy = r1 + getCoulomb(); // add coulomb potential;
00117 
00118   if ( verboseLevel >= 10 )
00119     G4cout << " G4BENeutronChannel : sampleKineticEnergy() " << G4endl
00120            << "       kinetic n e = " << kineticEnergy << G4endl
00121            << "        levelParam = " << levelParam << G4endl
00122            << "             thresh= " << getThresh() << G4endl;
00123 
00124   return kineticEnergy;
00125 }
00126 
00127 
00128 G4double G4BEChargedChannel::coulombFactorForProton()
00129 {
00130   //  Coefficient c_p:s for empirical cross section formula are
00131   //  defined with the proton constant.  See Dostrovsky, Phys. Rev.,
00132   //  vol. 116, 1959.
00133   G4double t[7] = { 0.08 , 0 , -0.06 , -0.1 , -0.1 , -0.1 , -0.1 };
00134   G4int Z = nucleusZ - particleZ;
00135 
00136   if ( Z >= 70.0 ) return t[6];
00137   if ( Z <= 10.0 ) return t[0];
00138   
00139   // Linear interpolation
00140   G4int   n = G4int( 0.1 * Z + 1.0 ); 
00141   G4float x = ( 10 * n - Z ) * 0.1; 
00142   G4double ret_val =  x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00143 
00144   return ret_val;
00145 }
00146 
00147 
00148 G4double G4BEChargedChannel::qmFactorForProton()
00149 {
00150   //  Coefficient k_p for empirical cross section formula are defined
00151   //  with the proton constant.  See Dostrovsky, Phys. Rev., vol. 116,
00152   //  1959
00153   G4double t[7] = {  0.36, 0.51, 0.60, 0.66, 0.68, 0.69, 0.69 };
00154   G4int Z = nucleusZ - particleZ;
00155 
00156   if ( Z >= 70.0 ) return t[6];
00157   if ( Z <= 10.0 ) return t[0];
00158 
00159   // Linear interpolation
00160   G4int   n = G4int( 0.1 * Z + 1.0 ); 
00161   G4float x = ( 10 * n - Z ) * 0.1; 
00162   return x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00163 }
00164 
00165 
00166 G4double G4BEChargedChannel::qmFactorForAlpha()
00167 {
00168 //  Coefficient k_alpha for empirical cross section formula presented
00169 //  in Dostrovsky, Phys. Rev., vol. 116, 1959
00170 
00171   G4double t[7] = {  0.77, 0.81, 0.85, 0.89, 0.93, 0.97,  1.00 };
00172   G4int Z = nucleusZ - particleZ;
00173 
00174   if ( Z >= 70.0 ) return t[6];
00175   if ( Z <= 10.0 ) return t[0];
00176 
00177   // Linear interpolation
00178   G4int   n = G4int( 0.1 * Z + 1.0 ); 
00179   G4float x = ( 10 * n - Z ) * 0.1; 
00180   return x * t[n - 2] + ( 1.0 - x ) * t[n-1];
00181 }

Generated on Mon May 27 17:47:43 2013 for Geant4 by  doxygen 1.4.7