G4DNAEmfietzoglouExcitationModel.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 // $Id$
00027 //
00028 
00029 #include "G4DNAEmfietzoglouExcitationModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAChemistryManager.hh"
00033 #include "G4DNAMolecularMaterial.hh"
00034 
00035 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00036 
00037 using namespace std;
00038 
00039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00040 
00041 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
00042                                                                    const G4String& nam)
00043     :G4VEmModel(nam),isInitialised(false)
00044 {
00045     //  nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00046     fpWaterDensity = 0;
00047 
00048     lowEnergyLimit = 8.23 * eV;
00049     highEnergyLimit = 10 * MeV;
00050     SetLowEnergyLimit(lowEnergyLimit);
00051     SetHighEnergyLimit(highEnergyLimit);
00052 
00053     nLevels = waterExcitation.NumberOfLevels();
00054 
00055     verboseLevel= 0;
00056     // Verbosity scale:
00057     // 0 = nothing
00058     // 1 = warning for energy non-conservation
00059     // 2 = details of energy budget
00060     // 3 = calculation of cross sections, file openings, sampling of atoms
00061     // 4 = entering in methods
00062 
00063     if (verboseLevel > 3)
00064         if( verboseLevel>0 )
00065         {
00066             G4cout << "Emfietzoglou Excitation model is constructed " << G4endl
00067                    << "Energy range: "
00068                    << lowEnergyLimit / eV << " eV - "
00069                    << highEnergyLimit / MeV << " MeV"
00070                    << G4endl;
00071         }
00072     fParticleChangeForGamma = 0;
00073 }
00074 
00075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00076 
00077 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
00078 {}
00079 
00080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00081 
00082 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
00083                                                   const G4DataVector& /*cuts*/)
00084 {
00085 
00086     if (verboseLevel > 3)
00087         G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
00088 
00089     // Energy limits
00090 
00091     if (LowEnergyLimit() < lowEnergyLimit)
00092     {
00093         G4cout << "G4DNAEmfietzoglouExcitationModel: low energy limit increased from " <<
00094                   LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00095         SetLowEnergyLimit(lowEnergyLimit);
00096     }
00097 
00098     if (HighEnergyLimit() > highEnergyLimit)
00099     {
00100         G4cout << "G4DNAEmfietzoglouExcitationModel: high energy limit decreased from " <<
00101                   HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
00102         SetHighEnergyLimit(highEnergyLimit);
00103     }
00104 
00105     //
00106     if( verboseLevel>0 )
00107     {
00108         G4cout << "Emfietzoglou Excitation model is initialized " << G4endl
00109                << "Energy range: "
00110                << LowEnergyLimit() / eV << " eV - "
00111                << HighEnergyLimit() / MeV << " MeV"
00112                << G4endl;
00113     }
00114 
00115     // Initialize water density pointer
00116     fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00117 
00118     if (isInitialised) { return; }
00119     fParticleChangeForGamma = GetParticleChangeForGamma();
00120     isInitialised = true;
00121 
00122 }
00123 
00124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00125 
00126 G4double G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume(const G4Material* material,
00127                                                                  const G4ParticleDefinition* particleDefinition,
00128                                                                  G4double ekin,
00129                                                                  G4double,
00130                                                                  G4double)
00131 {
00132     if (verboseLevel > 3)
00133         G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
00134 
00135     // Calculate total cross section for model
00136 
00137     G4double sigma=0;
00138 
00139     G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00140 
00141     if(waterDensity!= 0.0)
00142   //  if (material == nistwater || material->GetBaseMaterial() == nistwater)
00143     {
00144 
00145         if (particleDefinition == G4Electron::ElectronDefinition())
00146         {
00147             if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
00148             {
00149                 sigma = Sum(ekin);
00150             }
00151         }
00152 
00153         if (verboseLevel > 2)
00154         {
00155             G4cout << "__________________________________" << G4endl;
00156             G4cout << "°°° G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
00157             G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00158             G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00159             G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00160             //      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
00161             G4cout << "°°° G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
00162         }
00163 
00164     }
00165 
00166     return sigma*material->GetAtomicNumDensityVector()[1];
00167 }
00168 
00169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00170 
00171 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
00172                                                          const G4MaterialCutsCouple* /*couple*/,
00173                                                          const G4DynamicParticle* aDynamicElectron,
00174                                                          G4double,
00175                                                          G4double)
00176 {
00177 
00178     if (verboseLevel > 3)
00179         G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
00180 
00181     G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00182 
00183     G4int level = RandomSelect(electronEnergy0);
00184 
00185     G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
00186     G4double newEnergy = electronEnergy0 - excitationEnergy;
00187 
00188     if (electronEnergy0 < highEnergyLimit)
00189     {
00190             fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
00191             fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00192             fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00193 
00194         const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00195         G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
00196                                                                level,
00197                                                                theIncomingTrack);
00198     }
00199 }
00200 
00201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00202 
00203 G4double G4DNAEmfietzoglouExcitationModel::PartialCrossSection(G4double t, G4int level)
00204 {
00205     //                 Aj                        T
00206     // Sigma(T) = ------------- (Bj /  T) ln(Cj ---) [1 - Bj / T]^Pj
00207     //             2 pi alpha0                   R
00208     //
00209     // Sigma is the macroscopic cross section = N sigma, where N = number of target particles per unit volume
00210     // and sigma is the microscopic cross section
00211     // T      is the incoming electron kinetic energy
00212     // alpha0 is the Bohr Radius (Bohr_radius)
00213     // Aj, Bj, Cj & Pj are parameters that can be found in Emfietzoglou's papers
00214     //
00215     // From Phys. Med. Biol. 48 (2003) 2355-2371, D.Emfietzoglou,
00216     // Monte Carlo Simulation of the energy loss of low energy electrons in liquid Water
00217     //
00218     // Scaling for macroscopic cross section: number of water moleculs per unit volume
00219     // const G4double sigma0 = (10. / 3.343e22) * cm2;
00220 
00221     const G4double density = 3.34192e+19 * mm3;
00222 
00223     const G4double aj[]={0.0205, 0.0209, 0.0130, 0.0026, 0.0025};
00224     const G4double cj[]={4.9801, 3.3850, 2.8095, 1.9242, 3.4624};
00225     const G4double pj[]={0.4757, 0.3483, 0.4443, 0.3429, 0.4379};
00226     const G4double r = 13.6 * eV;
00227 
00228     G4double sigma = 0.;
00229 
00230     G4double exc = waterExcitation.ExcitationEnergy(level);
00231 
00232     if (t >= exc)
00233     {
00234         G4double excitationSigma = ( aj[level] / (2.*pi*Bohr_radius))
00235                 * (exc / t)
00236                 * std::log(cj[level]*(t/r))
00237                 * std::pow((1.- (exc/t)), pj[level]);
00238         sigma = excitationSigma / density;
00239     }
00240 
00241     return sigma;
00242 }
00243 
00244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00245 
00246 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k)
00247 {
00248     G4int i = nLevels;
00249     G4double value = 0.;
00250     std::deque<double> values;
00251 
00252     while (i > 0)
00253     {
00254         i--;
00255         G4double partial = PartialCrossSection(k,i);
00256         values.push_front(partial);
00257         value += partial;
00258     }
00259 
00260     value *= G4UniformRand();
00261     
00262     i = nLevels;
00263 
00264     while (i > 0)
00265     {
00266         i--;
00267         if (values[i] > value) return i;
00268         value -= values[i];
00269     }
00270     
00271     return 0;
00272 }
00273 
00274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00275 
00276 G4double G4DNAEmfietzoglouExcitationModel::Sum(G4double k)
00277 {
00278     G4double totalCrossSection = 0.;
00279 
00280     for (G4int i=0; i<nLevels; i++)
00281     {
00282         totalCrossSection += PartialCrossSection(k,i);
00283     }
00284     return totalCrossSection;
00285 }
00286 

Generated on Mon May 27 17:48:02 2013 for Geant4 by  doxygen 1.4.7