Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAQuinnPlasmonExcitationModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// Created on 2016/04/08
27//
28// Authors: D. Sakata, S. Incerti
29//
30// This class perform transmission term of volume plasmon excitation,
31// based on Quinn Model, see Phys. Rev. vol 126, number 4 (1962)
32
34#include "G4SystemOfUnits.hh"
35#include "G4RandomDirection.hh"
36
37//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38
39using namespace std;
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
45 const G4String& nam):
46G4VEmModel(nam)
47{
48 fpMaterialDensity = nullptr;
49 fLowEnergyLimit = 10 * eV;
50 fHighEnergyLimit = 1.0 * GeV;
51
52 for(int & i : nValenceElectron) i=0;
53
54 verboseLevel = 0;
55
56 if (verboseLevel > 0)
57 {
58 G4cout << "Quinn plasmon excitation model is constructed " << G4endl;
59 }
61 statCode = false;
62}
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
67= default;
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
72 (const G4ParticleDefinition* particle,
73 const G4DataVector& /*cuts*/)
74{
75 for(int & i : nValenceElectron) i=0;
76
77 if (verboseLevel > 3)
78 {
79 G4cout <<
80 "Calling G4DNAQuinnPlasmonExcitationModel::Initialise()"
81 << G4endl;
82 }
83
84 if(particle == G4Electron::ElectronDefinition())
85 {
86 fLowEnergyLimit = 10 * eV;
87 fHighEnergyLimit = 1.0 * GeV;
88 }
89 else
90 {
91 G4Exception("G4DNAQuinnPlasmonExcitationModel::Initialise","em0001",
92 FatalException,"Not defined for other particles than electrons.");
93 return;
94 }
95
96 // Get Number of valence electrons
97 G4ProductionCutsTable* theCoupleTable =
99
100 auto numOfCouples = (G4int)theCoupleTable->GetTableSize();
101
102 for(G4int i=0;i<numOfCouples;i++){
103
104 const G4MaterialCutsCouple* couple =
105 theCoupleTable->GetMaterialCutsCouple(i);
106
107 const G4Material* material = couple->GetMaterial();
108
109 const G4ElementVector* theElementVector =material->GetElementVector();
110
111 std::size_t nelm = material->GetNumberOfElements();
112 if (nelm==1) // Protection: only for single element
113 {
114 G4int z = G4lrint((*theElementVector)[0]->GetZ());
115 if(z<=100)
116 {
117 nValenceElectron[z] = GetNValenceElectron(z);
118 }
119 else
120 {
121 G4Exception("G4DNAQuinnPlasmonExcitationModel::Initialise","em0002",
122 FatalException,"The model is not applied for z>100");
123 }
124 }
125 //for(G4int j=0;j<nelm;j++){
126 // G4int z=G4lrint((*theElementVector)[j]->GetZ());
127 // if(z<=100){nValenceElectron[z] = GetNValenceElectron(z);}
128 //}
129 }
130
131 if( verboseLevel>0 )
132 {
133 G4cout << "Quinn plasmon excitation model is initialized " << G4endl
134 << "Energy range: "
135 << LowEnergyLimit() / eV << " eV - "
136 << HighEnergyLimit() / keV << " keV for "
137 << particle->GetParticleName()
138 << G4endl;
139 }
140
141 if (isInitialised){return;}
143 isInitialised = true;
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
149 (const G4Material* material,
150 const G4ParticleDefinition* particleDefinition,
151 G4double ekin,
152 G4double,
153 G4double)
154{
155 if (verboseLevel > 3)
156 {
157 G4cout <<
158 "Calling CrossSectionPerVolume() of G4DNAQuinnPlasmonExcitationModel"
159 << G4endl;
160 }
161
162 // Protection: only for single element
163 if(material->GetNumberOfElements()>1) return 0.;
164 G4double z = material->GetZ();
165
166 // Protection: only for Gold
167 if (z!=79){return 0.;}
168
169
170 G4double sigma = 0;
171 G4double atomicNDensity = material->GetAtomicNumDensityVector()[0];
172
173 if(atomicNDensity!= 0.0)
174 {
175 if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
176 {
177 sigma = GetCrossSection(material,particleDefinition,ekin);
178 }
179
180 if (verboseLevel > 2)
181 {
182 G4cout<<"__________________________________" << G4endl;
183 G4cout<<"=== G4DNAQuinnPlasmonExcitationModel - XS INFO START"<<G4endl;
184 G4cout<<"=== Kinetic energy (eV)=" << ekin/eV << " particle : "
185 <<particleDefinition->GetParticleName() << G4endl;
186 G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^2)"
187 <<sigma/cm/cm << G4endl;
188 G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^-1)="
189 <<sigma*atomicNDensity/(1./cm) << G4endl;
190 G4cout<<"=== G4DNAQuinnPlasmonExcitationModel - XS INFO END" << G4endl;
191 }
192 }
193
194 return sigma*atomicNDensity;
195}
196
197//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198
200 (std::vector<G4DynamicParticle*>* /*fvect*/,
201 const G4MaterialCutsCouple* couple,
202 const G4DynamicParticle* aDynamicParticle,
204{
205
206 if (verboseLevel > 3)
207 {
208 G4cout <<
209 "Calling SampleSecondaries() of G4DNAQuinnPlasmonExcitationModel"
210 << G4endl;
211 }
212
213 const G4Material *material = couple->GetMaterial();
214
215 G4ParticleDefinition* particle = aDynamicParticle->GetDefinition();
216
217 G4double k = aDynamicParticle->GetKineticEnergy();
218
219 if(particle == G4Electron::ElectronDefinition())
220 {
221 G4double e = 1.;
222 G4int z = material->GetZ();
223 G4int Nve = 0;
224
225 //TODO: have to be change to realistic!!
226 if(z<100) Nve = nValenceElectron[z];
227
228 G4double A = material->GetA()/g/mole;
229 G4double Dens = material->GetDensity()/g*cm*cm*cm;
230 G4double veDens = Dens*CLHEP::Avogadro*Nve/A;
231
232 G4double omega_p = std::sqrt(veDens*std::pow(e,2)/
233 (CLHEP::epsilon0/(1./cm)*CLHEP::electron_mass_c2
234 /(CLHEP::c_squared/cm/cm)));
235
236 G4double excitationEnergy = CLHEP::hbar_Planck*omega_p;
237 G4double newEnergy = k - excitationEnergy;
238
239
240 if (newEnergy > 0)
241 {
243 ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
244
246
247 if(!statCode)
248 {
250 }
251 else
252 {
254
255 }
256 }
257 }
258}
259
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261
263 (const G4Material* material,
264 const G4ParticleDefinition* particle,
265 G4double kineticEnergy)
266{
267 G4double value=0;
268
269 if(particle == G4Electron::ElectronDefinition())
270 {
271 G4double e = 1.;
272 G4int z = material->GetZ();
273 G4int Nve = 0;
274 if(z<100) Nve = nValenceElectron[z];
275 G4double A = material->GetA()/g/mole;
276 G4double Dens = material->GetDensity()/g*cm*cm*cm;
277 G4double veDens = Dens*CLHEP::Avogadro*Nve/A;
278
279 G4double omega_p = std::sqrt(veDens*std::pow(e,2)
280 /(CLHEP::epsilon0/(1./cm)*CLHEP::electron_mass_c2/
281 (CLHEP::c_squared/cm/cm)));
282
283 G4double fEnergy = std::pow(CLHEP::h_Planck,2)/(8*CLHEP::electron_mass_c2)*
284 std::pow(3*veDens/CLHEP::pi,2./3.)/e
285 *(CLHEP::c_squared/cm/cm);
286
287 G4double p0 = sqrt(2*CLHEP::electron_mass_c2
288 /(CLHEP::c_squared/cm/cm)*fEnergy);
289
290 G4double p = sqrt(2*CLHEP::electron_mass_c2
291 /(CLHEP::c_squared/cm/cm)*kineticEnergy);
292
293 G4double mfp = 2*CLHEP::Bohr_radius/cm*kineticEnergy
294 /(CLHEP::hbar_Planck*omega_p)/
295 (G4Log((std::pow(std::pow(p0,2)
296 +2*CLHEP::electron_mass_c2/
297 (CLHEP::c_squared/cm/cm)*omega_p
298 *CLHEP::hbar_Planck,1./2.)-p0)
299 /(p-std::pow(std::pow(p,2)-2*CLHEP::electron_mass_c2/
300 (CLHEP::c_squared/cm/cm)*omega_p
301 *CLHEP::hbar_Planck,1./2.))));
302
303 G4double excitationEnergy = CLHEP::hbar_Planck*omega_p;
304
305 if((0<mfp)&&(0<veDens)&&(excitationEnergy<kineticEnergy)){
306 value = 1./(veDens*mfp);
307 }
308 }
309 return value*cm*cm;
310}
311
312//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
313
314G4int G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron(G4int z)
315{
316
317 G4int Nve=0;
318
319 // Current limitation to gold
320 if (z!=79){return 0.;}
321
322 if (verboseLevel > 3)
323 {
324 G4cout <<
325 "Calling GetNValenceElectron() of G4DNAQuinnPlasmonExcitationModel"
326 << G4endl;
327 }
328
329 const char *datadir=nullptr;
330
331 if(datadir == nullptr)
332 {
333 datadir = G4FindDataDir("G4LEDATA");
334 if(datadir == nullptr)
335 {
336 G4Exception("G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron()"
337 ,"em0002",FatalException,
338 "Enviroment variable G4LEDATA not defined");
339 return 0;
340 }
341 }
342
343 std::ostringstream targetfile;
344 targetfile.str("");
345 targetfile.clear(stringstream::goodbit);
346 targetfile << datadir <<"/dna/atomicstate_Z"<< z <<".dat";
347 std::ifstream fin(targetfile.str().c_str());
348
349 if(!fin)
350 {
351 G4cout<< " Error : "<< targetfile.str() <<" is not found "<<endl;
352 G4Exception("G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron()"
353 ,"em0003",FatalException,
354 "There is no target file");
355 return 0;
356 }
357
358 string buff0,buff1,buff2,buff3,buff4,buff5,buff6;
359 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
360
361 while(true){
362 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
363 if(!fin.eof())
364 {
365 Nve = stoi(buff3);
366 }
367 else
368 {
369 break;
370 }
371 }
372 return Nve;
373}
374
std::vector< const G4Element * > G4ElementVector
const char * G4FindDataDir(const char *)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Log(G4double x)
Definition G4Log.hh:227
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
const G4double A[17]
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector())) override
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
virtual G4double GetCrossSection(const G4Material *material, const G4ParticleDefinition *, G4double kineticEnergy)
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4DNAQuinnPlasmonExcitationModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="DNAQuinnPlasmonExcitationModel")
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
const G4Material * GetMaterial() const
G4double GetDensity() const
const G4ElementVector * GetElementVector() const
G4double GetZ() const
const G4double * GetAtomicNumDensityVector() const
G4double GetA() const
std::size_t GetNumberOfElements() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
const G4String & GetParticleName() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
int G4lrint(double ad)
Definition templates.hh:134