00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4DNAEmfietzoglouExcitationModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAChemistryManager.hh"
00033 #include "G4DNAMolecularMaterial.hh"
00034
00035
00036
00037 using namespace std;
00038
00039
00040
00041 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
00042 const G4String& nam)
00043 :G4VEmModel(nam),isInitialised(false)
00044 {
00045
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
00057
00058
00059
00060
00061
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
00076
00077 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
00078 {}
00079
00080
00081
00082 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* ,
00083 const G4DataVector& )
00084 {
00085
00086 if (verboseLevel > 3)
00087 G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
00088
00089
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
00116 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00117
00118 if (isInitialised) { return; }
00119 fParticleChangeForGamma = GetParticleChangeForGamma();
00120 isInitialised = true;
00121
00122 }
00123
00124
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
00136
00137 G4double sigma=0;
00138
00139 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00140
00141 if(waterDensity!= 0.0)
00142
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
00161 G4cout << "°°° G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
00162 }
00163
00164 }
00165
00166 return sigma*material->GetAtomicNumDensityVector()[1];
00167 }
00168
00169
00170
00171 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00172 const G4MaterialCutsCouple* ,
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
00202
00203 G4double G4DNAEmfietzoglouExcitationModel::PartialCrossSection(G4double t, G4int level)
00204 {
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
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
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
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