Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNARelativisticIonisationModel.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// $Id: G4DNARelativisticIonisationModel.cc $
27//
28// Created on 2016/05/12
29//
30// Authors: D Sakata, S. Incerti
31//
32// This class perform ionisation for electron transportation in gold,
33// based on Relativistic Binary Encounter Bethe-Vriens(RBEBV) model.
34// See following reference paper,
35// M. Guerra et al, J. Phys. B: At. Mol. Opt. Phys. 48, 185202 (2015)
36// =======================================================================
37// Limitation of secondaries by GEANT4 atomic de-excitation:
38// The cross section and energy of secondary production is based on
39// EADL database. If there are no tabele for several orbitals, this class
40// will not provide secondaries for the orbitals.
41// For gold(Au), this class provide secondaries for inner 18 orbitals
42// but don't provide for outer 3 orbitals due to EADL databese limitation.
43// =======================================================================
44
45#include "G4DNARelativisticIonisationModel.hh"
46#include "G4SystemOfUnits.hh"
47#include "G4AtomicShell.hh"
48#include "G4UAtomicDeexcitation.hh"
49#include "G4LossTableManager.hh"
50#include "G4Gamma.hh"
51#include "G4RandomDirection.hh"
52
53#include "G4DNAMolecularMaterial.hh"
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57using namespace std;
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61G4DNARelativisticIonisationModel::G4DNARelativisticIonisationModel(
62 const G4ParticleDefinition*,
63 const G4String& nam) :
64G4VEmModel(nam), isInitialised(false),statCode(false),fasterCode(true)
65{
66 fHighEnergyLimit = 0;
67 fLowEnergyLimit = 0;
68
69 verboseLevel = 0;
70
71 SetDeexcitationFlag(true);
72 fAtomDeexcitation = 0;
73 fMaterialDensity = 0;
74 fParticleDefinition = 0;
75 fParticleChangeForGamma = 0;
76
77 if (verboseLevel > 0)
78 {
79 G4cout << "Relativistic Ionisation Model is constructed " << G4endl;
80 }
81}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85G4DNARelativisticIonisationModel::~G4DNARelativisticIonisationModel()
86{
87 // Cross section
88}
89
90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91
92void G4DNARelativisticIonisationModel::Initialise(const G4ParticleDefinition* particle,
93 const G4DataVector& /*cuts*/)
94{
95
96 if (verboseLevel > 3)
97 {
98 G4cout <<
99 "Calling G4DNARelativisticIonisationModel::Initialise()"
100 << G4endl;
101 }
102
103
104 if(fParticleDefinition != 0 && fParticleDefinition != particle)
105 {
106 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0001",
107 FatalException,"Model already initialized for another particle type.");
108 }
109
110 fParticleDefinition = particle;
111 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
112 if(particle == electronDef)
113 {
114 fLowEnergyLimit = 10 * eV;
115 fHighEnergyLimit = 1.0 * GeV;
116
117 std::ostringstream eFullFileNameZ;
118
119 const char *path = G4FindDataDir("G4LEDATA");
120 if (!path)
121 {
122 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0006",
123 FatalException,"G4LEDATA environment variable not set.");
124 return;
125 }
126
127
128 G4ProductionCutsTable *coupletable
129 = G4ProductionCutsTable::GetProductionCutsTable();
130 G4int Ncouple = (G4int)coupletable ->GetTableSize();
131 for(G4int i=0; i<Ncouple; ++i)
132 {
133 const G4MaterialCutsCouple* couple
134 = coupletable->GetMaterialCutsCouple(i);
135 const G4Material * material = couple ->GetMaterial();
136 {
137 // Protection: only for single element
138 if(material->GetNumberOfElements()>1) continue;
139
140 G4int Z = material->GetZ();
141 // Protection: only for GOLD
142 if(Z!=79) continue;
143
144 iState [Z].clear();
145 iShell [Z].clear();
146 iSubShell [Z].clear();
147 Nelectrons[Z].clear();
148 Ebinding [Z].clear();
149 Ekinetic [Z].clear();
150 LoadAtomicStates(Z,path);
151
152 /////////////Load cumulated DCS////////////////
153 eVecEZ.clear();
154 eVecEjeEZ.clear();
155 eProbaShellMapZ.clear();
156 eDiffCrossSectionDataZ.clear();
157
158 eFullFileNameZ.str("");
159 eFullFileNameZ.clear(stringstream::goodbit);
160
161 eFullFileNameZ
162 << path
163 << "/dna/sigmadiff_cumulated_ionisation_e_RBEBV_Z"
164 << Z << ".dat";
165 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
166 if (!eDiffCrossSectionZ)
167 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0003",
168 FatalException,
169 "Missing data file for cumulated DCS");
170
171 eVecEZ[Z].push_back(0.);
172 while(!eDiffCrossSectionZ.eof())
173 {
174 G4double tDummy;
175 G4double eDummy;
176 eDiffCrossSectionZ>>tDummy>>eDummy;
177 if (tDummy != eVecEZ[Z].back())
178 {
179 eVecEZ[Z].push_back(tDummy);
180 eVecEjeEZ[Z][tDummy].push_back(0.);
181 }
182
183 for(G4int istate=0;istate<(G4int)iState[Z].size();istate++)
184 {
185 eDiffCrossSectionZ>>
186 eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy];
187 eEjectedEnergyDataZ[Z][istate][tDummy]
188 [eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy]]
189 = eDummy;
190 eProbaShellMapZ[Z][istate][tDummy].push_back(
191 eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy]);
192 }
193
194 if (eDummy != eVecEjeEZ[Z][tDummy].back()){
195 eVecEjeEZ[Z][tDummy].push_back(eDummy);
196 }
197 }
198 }
199 }
200 }
201 else
202 {
203 G4cout<<
204 "Error : No particle Definition is found in G4DNARelativisticIonisationModel"
205 <<G4endl;
206 return;
207 }
208
209 if( verboseLevel>0 )
210 {
211 G4cout << "Relativistic Ionisation model is initialized " << G4endl
212 << "Energy range: "
213 << LowEnergyLimit() / eV << " eV - "
214 << HighEnergyLimit() / keV << " keV for "
215 << particle->GetParticleName()
216 << G4endl;
217 }
218
219 // Initialise gold density pointer
220 fMaterialDensity = G4DNAMolecularMaterial::Instance()
221 ->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_Au"));
222
223 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
224 fParticleChangeForGamma = GetParticleChangeForGamma();
225
226 if (isInitialised){return;}
227 isInitialised = true;
228}
229
230//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231
232G4double G4DNARelativisticIonisationModel::CrossSectionPerVolume(
233 const G4Material* material,
234 const G4ParticleDefinition* particleDefinition,
235 G4double ekin,
236 G4double,
237 G4double)
238{
239 if (verboseLevel > 3)
240 {
241 G4cout <<
242 "Calling CrossSectionPerVolume() of G4DNARelativisticIonisationModel"
243 << G4endl;
244 }
245
246 if(particleDefinition != fParticleDefinition) return 0;
247
248 // Calculate total cross section for model
249 G4double sigma=0;
250
251 if(material->GetNumberOfElements()>1) return 0.; // Protection for Molecules
252 G4double atomicNDensity = material->GetAtomicNumDensityVector()[0];
253 G4double z = material->GetZ();
254
255 if(atomicNDensity!= 0.0)
256 {
257 if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
258 {
259 sigma = GetTotalCrossSection(material,particleDefinition,ekin);
260 }
261
262 if (verboseLevel > 2)
263 {
264 G4cout << "__________________________________" << G4endl;
265 G4cout << "=== G4DNARelativisticIonisationModel - XS INFO START" <<G4endl;
266 G4cout << "=== Kinetic energy (eV)=" << ekin/eV << " particle : "
267 << particleDefinition->GetParticleName() << G4endl;
268 G4cout << "=== Cross section per atom for Z="<<z<<" is (cm^2)"
269 << sigma/cm/cm << G4endl;
270 G4cout << "=== Cross section per atom for Z="<<z<<" is (cm^-1)="
271 << sigma*atomicNDensity/(1./cm) << G4endl;
272 G4cout << "=== G4DNARelativisticIonisationModel - XS INFO END" << G4endl;
273 }
274 }
275 return sigma*atomicNDensity;
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279
280void G4DNARelativisticIonisationModel::SampleSecondaries(
281 std::vector<G4DynamicParticle*>* fvect,
282 const G4MaterialCutsCouple* couple,
283 const G4DynamicParticle* particle,
284 G4double,G4double)
285{
286 if (verboseLevel > 3)
287 {
288 G4cout <<
289 "Calling SampleSecondaries() of G4DNARelativisticIonisationModel"
290 << G4endl;
291 }
292
293
294 G4ParticleDefinition* particleDef = particle->GetDefinition();
295 G4double k = particle->GetKineticEnergy();
296 G4double ejectedE = 0.*eV;
297
298 if(fLowEnergyLimit <= k && k<fHighEnergyLimit)
299 {
300 G4ThreeVector primaryDir = particle ->GetMomentumDirection();
301
302 G4double particleMass = particleDef->GetPDGMass();
303 G4double totalEnergy = k+particleMass;
304 G4double pSquare = k*(totalEnergy+particleMass);
305 G4double totalMomentum = std::sqrt(pSquare);
306
307 const G4Material *material = couple->GetMaterial();
308 G4int z = material->GetZ();
309 G4int level = RandomSelect(material,particleDef,k);
310
311 if(k<Ebinding[z].at(level)) return;
312
313 G4int NumSecParticlesInit =0;
314 G4int NumSecParticlesFinal=0;
315
316 if(fAtomDeexcitation){
317 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(level);
318 const G4AtomicShell *shell = fAtomDeexcitation->GetAtomicShell(z,as);
319 NumSecParticlesInit = (G4int)fvect->size();
320 fAtomDeexcitation->GenerateParticles(fvect,shell,z,0,0);
321 NumSecParticlesFinal = (G4int)fvect->size();
322 }
323
324 ejectedE
325 = GetEjectedElectronEnergy (material,particleDef,k,level);
326 G4ThreeVector ejectedDir
327 = GetEjectedElectronDirection(particleDef,k,ejectedE);
328 ejectedDir.rotateUz(primaryDir);
329
330 G4double scatteredE = k - Ebinding[z].at(level) - ejectedE;
331
332 if(particleDef == G4Electron::ElectronDefinition()){
333 G4double secondaryTotMomentum
334 = std::sqrt(ejectedE*(ejectedE+2*CLHEP::electron_mass_c2));
335 G4double finalMomentumX
336 = totalMomentum*primaryDir.x()- secondaryTotMomentum*ejectedDir.x();
337 G4double finalMomentumY
338 = totalMomentum*primaryDir.y()- secondaryTotMomentum*ejectedDir.y();
339 G4double finalMomentumZ
340 = totalMomentum*primaryDir.z()- secondaryTotMomentum*ejectedDir.z();
341
342 G4ThreeVector scatteredDir(finalMomentumX,finalMomentumY,finalMomentumZ);
343 fParticleChangeForGamma->ProposeMomentumDirection(scatteredDir.unit());
344
345 }
346 else
347 {
348 fParticleChangeForGamma->ProposeMomentumDirection(primaryDir);
349 }
350
351 //G4double deexSecEnergy=0.;
352 G4double restEproduction = Ebinding[z].at(level);
353 for(G4int iparticle=NumSecParticlesInit;
354 iparticle<NumSecParticlesFinal;iparticle++)
355 {
356 //deexSecEnergy = deexSecEnergy + (*fvect)[iparticle]->GetKineticEnergy();
357 G4double Edeex = (*fvect)[iparticle]->GetKineticEnergy();
358 if(restEproduction>=Edeex){
359 restEproduction -= Edeex;
360 }
361 else{
362 delete (*fvect)[iparticle];
363 (*fvect)[iparticle]=0;
364 }
365 }
366 if(restEproduction < 0.0){
367 G4Exception("G4DNARelativisticIonisationModel::SampleSecondaries()",
368 "em0008",FatalException,"Negative local energy deposit");
369 }
370
371 if(!statCode)
372 {
373 if(scatteredE>0){
374 fParticleChangeForGamma->SetProposedKineticEnergy (scatteredE);
375 fParticleChangeForGamma->ProposeLocalEnergyDeposit(restEproduction);
376 //fParticleChangeForGamma
377 //->ProposeLocalEnergyDeposit(k-scatteredE-ejectedE-deexSecEnergy);
378 }
379 }
380 else
381 {
382 fParticleChangeForGamma->SetProposedKineticEnergy (k);
383 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredE);
384 }
385
386 if(ejectedE>0){
387 G4DynamicParticle* ejectedelectron
388 = new G4DynamicParticle(G4Electron::Electron(),ejectedDir,ejectedE);
389 fvect->push_back(ejectedelectron);
390 }
391 }
392}
393
394void G4DNARelativisticIonisationModel::LoadAtomicStates(
395 G4int z,const char* path)
396{
397
398 if (verboseLevel > 3)
399 {
400 G4cout <<
401 "Calling LoadAtomicStates() of G4DNARelativisticIonisationModel"
402 << G4endl;
403 }
404 const char *datadir = path;
405 if(!datadir)
406 {
407 datadir = G4FindDataDir("G4LEDATA");
408 if(!datadir)
409 {
410 G4Exception("G4DNARelativisticIonisationModel::LoadAtomicStates()",
411 "em0002",FatalException,"Enviroment variable G4LEDATA not defined");
412
413 return;
414 }
415 }
416 std::ostringstream targetfile;
417 targetfile << datadir <<"/dna/atomicstate_Z"<< z <<".dat";
418 std::ifstream fin(targetfile.str().c_str());
419 if(!fin)
420 {
421 G4cout<< " Error : "<< targetfile.str() <<" is not found "<<G4endl;
422 G4Exception("G4DNARelativisticIonisationModel::LoadAtomicStates()","em0002",
423 FatalException,"There is no target file");
424 return;
425 }
426
427 G4String buff0,buff1,buff2,buff3,buff4,buff5,buff6;
428 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
429 G4int iline=0;
430 while(true){
431 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
432 if(!fin.eof())
433 {
434 iState [z].push_back(stoi(buff0));
435 iShell [z].push_back(stoi(buff1));
436 iSubShell [z].push_back(stoi(buff2));
437 Nelectrons[z].push_back(stoi(buff3));
438 Ebinding [z].push_back(stod(buff4));
439 if(stod(buff5)==0.)
440 {// if there is no kinetic energy in the file, kinetic energy
441 // for Bhor atomic model will be calculated: !!! I's not realistic!!!
442 G4double radius = std::pow(iShell[z].at(iline),2)
443 *std::pow(CLHEP::hbar_Planck,2)*(4*CLHEP::pi*CLHEP::epsilon0)
444 /CLHEP::electron_mass_c2;
445 G4double momentum = iShell[z].at(iline)*CLHEP::hbar_Planck/radius;
446 Ekinetic[z].push_back(std::pow(momentum,2)/(2*CLHEP::electron_mass_c2));
447 }
448 else
449 {
450 Ekinetic [z].push_back(stod(buff5));
451 }
452 iline++;
453 }
454 else
455 {
456 break;
457 }
458 }
459}
460
461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
462
463G4double G4DNARelativisticIonisationModel::GetTotalCrossSection(
464 const G4Material* material,
465 const G4ParticleDefinition* particle,
466 G4double kineticEnergy)
467{
468 G4double value=0;
469 G4int z = material->GetZ();
470 if(z!=79){ return 0.;}
471 else {
472 std::size_t N=iState[z].size();
473 for(G4int i=0; i<(G4int)N; ++i){
474 value = value+GetPartialCrossSection(material,i,particle,kineticEnergy);
475 }
476 return value;
477 }
478}
479
480//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
481
482G4double G4DNARelativisticIonisationModel::GetPartialCrossSection(
483 const G4Material* material,
484 G4int level,
485 const G4ParticleDefinition* particle,
486 G4double kineticEnergy)
487{
488 G4double value = 0;
489 G4double constRy =13.6057E-6;//MeV
490
491 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
492 G4int z = material->GetZ();
493 if(particle==electronDef){
494
495 G4double t = kineticEnergy /Ebinding[z].at(level);
496 G4double tdash = kineticEnergy /CLHEP::electron_mass_c2;
497 G4double udash = Ekinetic[z].at(level)/CLHEP::electron_mass_c2;
498 G4double bdash = Ebinding[z].at(level)/CLHEP::electron_mass_c2;
499 G4double beta_t2 = 1.-1./std::pow(1.+tdash,2);
500 G4double beta_u2 = 1.-1./std::pow(1.+udash,2);
501 G4double beta_b2 = 1.-1./std::pow(1.+bdash,2);
502 G4double alpha = std::sqrt(2*constRy/CLHEP::electron_mass_c2);
503 G4double phi = std::cos(std::sqrt(std::pow(alpha,2)
504 /(beta_t2+beta_b2))*G4Log(beta_t2/beta_b2));
505 G4double constS = 4*CLHEP::pi*std::pow(CLHEP::Bohr_radius,2)
506 *Nelectrons[z].at(level)*std::pow(alpha,4);
507
508 if(Ebinding[z].at(level)<=kineticEnergy)
509 {
510 value =constS/((beta_t2+(beta_u2+beta_b2)/iShell[z].at(level))*2.*bdash)
511 *(1./2.*(G4Log(beta_t2/(1.-beta_t2))-beta_t2-G4Log(2.*bdash))
512 *(1.-1./std::pow(t,2.))
513 +1.-1./t-G4Log(t)/(t+1.)*(1.+2.*tdash)/(std::pow(1.+tdash/2.,2.))
514 *phi+std::pow(bdash,2)/(std::pow(1+tdash/2.,2))*(t-1)/2.);
515 }
516
517 }
518 return value;
519}
520
521//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
522
523G4double G4DNARelativisticIonisationModel::GetDifferentialCrossSection(
524 const G4Material* material,
525 const G4ParticleDefinition* particle,
526 G4double kineticEnergy,
527 G4double secondaryEnergy,
528 G4int level)
529{
530 G4double value=0.;
531 G4double constRy =13.6057E-6;//MeV
532
533 G4int z = material->GetZ();
534
535 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
536 if(particle==electronDef){
537 G4double w = secondaryEnergy /Ebinding[z].at(level);
538 G4double t = kineticEnergy /Ebinding[z].at(level);
539 G4double tdash = kineticEnergy /CLHEP::electron_mass_c2;
540 G4double udash = Ekinetic[z].at(level)/CLHEP::electron_mass_c2;
541 G4double bdash = Ebinding[z].at(level)/CLHEP::electron_mass_c2;
542 G4double beta_t2 = 1.-1./std::pow(1.+tdash,2);
543 G4double beta_u2 = 1.-1./std::pow(1.+udash,2);
544 G4double beta_b2 = 1.-1./std::pow(1.+bdash,2);
545 G4double alpha = std::sqrt(2*constRy/CLHEP::electron_mass_c2);
546 G4double phi = std::cos(std::sqrt(std::pow(alpha,2)/(beta_t2+beta_b2))
547 *G4Log(beta_t2/beta_b2));
548 G4double constS = 4*CLHEP::pi*std::pow(CLHEP::Bohr_radius,2)
549 *Nelectrons[z].at(level)*std::pow(alpha,4);
550
551 if(secondaryEnergy<=((kineticEnergy-Ebinding[z].at(level))/2.))
552 {
553 value = constS/((beta_t2+(beta_u2+beta_b2)/iShell[z].at(level))*2.*bdash)
554 *(-phi/(t+1.)*(1./std::pow(w+1.,1.)+1./std::pow(t-w,1.))
555 *(1.+2*tdash)/std::pow(1.+tdash/2.,2.)
556 +1./std::pow(w+1.,2.)+1./std::pow(t-w,2.)
557 +std::pow(bdash,2)/std::pow(1+tdash/2.,2)
558 +(1./std::pow(w+1.,3.)+1./std::pow(t-w,3.))
559 *(G4Log(beta_t2/(1.-beta_t2))-beta_t2-G4Log(2*bdash)));
560 }
561 }
562 return value;
563}
564
565//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
566
567G4int G4DNARelativisticIonisationModel::RandomSelect(
568 const G4Material* material,
569 const G4ParticleDefinition* particle,
570 G4double kineticEnergy)
571{
572 G4double value = 0.;
573 G4int z = material->GetZ();
574 std::size_t numberOfShell = iShell[z].size();
575 auto valuesBuffer = new G4double[numberOfShell];
576 const G4int n = (G4int)iShell[z].size();
577 G4int i(n);
578
579 while (i > 0)
580 {
581 i--;
582 if((fLowEnergyLimit<=kineticEnergy)&&(kineticEnergy<fHighEnergyLimit))
583 {
584 valuesBuffer[i]=GetPartialCrossSection(material,i,particle,kineticEnergy);
585 }
586 value += valuesBuffer[i];
587 }
588
589 value *= G4UniformRand();
590 i = n;
591
592 while (i > 0)
593 {
594 i--;
595
596 if (valuesBuffer[i] > value)
597 {
598 delete[] valuesBuffer;
599 return i;
600 }
601 value -= valuesBuffer[i];
602 }
603
604 if (valuesBuffer) delete[] valuesBuffer;
605
606 return 9999;
607}
608
609//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
610
611
612G4double G4DNARelativisticIonisationModel::GetEjectedElectronEnergy(
613 const G4Material* material,
614 const G4ParticleDefinition* particle,
615 G4double energy, G4int ishell)
616{
617 G4double secondaryEnergy=0;
618
619 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
620 G4int z = material->GetZ();
621 if(!fasterCode){ // for 2D rejection method
622 if(particle==electronDef){
623 G4double maximumsecondaryEnergy = (energy-Ebinding[z].at(ishell))/2.;
624 if(maximumsecondaryEnergy<0.) return 0.;
625 G4double maximumCrossSection=-999.;
626
627 maximumCrossSection
628 = GetDifferentialCrossSection(material,particle,energy,0.,ishell);
629 do{
630 secondaryEnergy = G4UniformRand()* maximumsecondaryEnergy;
631 }while(G4UniformRand()*maximumCrossSection >
632 GetDifferentialCrossSection(
633 material,particle,energy,secondaryEnergy,ishell));
634 }
635 }
636 else { // for cumulative method using cumulated DCS file
637
638 G4double valueE1 =0.;
639 G4double valueE2 =0.;
640 G4double valueXS21=0.;
641 G4double valueXS22=0.;
642 G4double valueXS11=0.;
643 G4double valueXS12=0.;
644 G4double ejeE21 =0.;
645 G4double ejeE22 =0.;
646 G4double ejeE11 =0.;
647 G4double ejeE12 =0.;
648 G4double random = G4UniformRand();
649
650 if (particle == G4Electron::ElectronDefinition())
651 {
652 if((eVecEZ[z].at(0)<=energy)&&(energy<eVecEZ[z].back()))
653 {
654 std::vector<G4double>::iterator k2
655 = std::upper_bound(eVecEZ[z].begin(),eVecEZ[z].end(), energy);
656 std::vector<G4double>::iterator k1 = k2-1;
657
658 if ( random < eProbaShellMapZ[z][ishell][(*k1)].back()
659 && random < eProbaShellMapZ[z][ishell][(*k2)].back() )
660 {
661 std::vector<G4double>::iterator xs12 =
662 std::upper_bound(eProbaShellMapZ[z][ishell][(*k1)].begin(),
663 eProbaShellMapZ[z][ishell][(*k1)].end(), random);
664 std::vector<G4double>::iterator xs11 = xs12-1;
665
666 std::vector<G4double>::iterator xs22 =
667 std::upper_bound(eProbaShellMapZ[z][ishell][(*k2)].begin(),
668 eProbaShellMapZ[z][ishell][(*k2)].end(), random);
669 std::vector<G4double>::iterator xs21 = xs22-1;
670
671 valueE1 =*k1;
672 valueE2 =*k2;
673 valueXS21 =*xs21;
674 valueXS22 =*xs22;
675 valueXS12 =*xs12;
676 valueXS11 =*xs11;
677
678 ejeE11 = eEjectedEnergyDataZ[z][ishell][valueE1][valueXS11];
679 ejeE12 = eEjectedEnergyDataZ[z][ishell][valueE1][valueXS12];
680 ejeE21 = eEjectedEnergyDataZ[z][ishell][valueE2][valueXS21];
681 ejeE22 = eEjectedEnergyDataZ[z][ishell][valueE2][valueXS22];
682
683 secondaryEnergy = QuadInterpolator( valueXS11, valueXS12,
684 valueXS21, valueXS22,
685 ejeE11 , ejeE12 ,
686 ejeE21 , ejeE22 ,
687 valueE1, valueE2,
688 energy, random );
689 }
690 }
691 }
692 }
693
694 if(secondaryEnergy<0) secondaryEnergy=0;
695 return secondaryEnergy;
696}
697
698//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
699
700G4ThreeVector G4DNARelativisticIonisationModel::GetEjectedElectronDirection(
701 const G4ParticleDefinition* ,
702 G4double energy,G4double secondaryenergy)
703{
704 G4double phi = 2*CLHEP::pi*G4UniformRand();
705 G4double sintheta = std::sqrt((1.-secondaryenergy/energy)
706 / (1.+secondaryenergy/(2*CLHEP::electron_mass_c2)));
707
708 G4double dirX = sintheta*std::cos(phi);
709 G4double dirY = sintheta*std::sin(phi);
710 G4double dirZ = std::sqrt(1.-sintheta*sintheta);
711
712 G4ThreeVector vec(dirX,dirY,dirZ);
713 return vec;
714}
715
716//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
717
718G4double G4DNARelativisticIonisationModel::Interpolate( G4double e1,
719 G4double e2,
720 G4double e,
721 G4double xs1,
722 G4double xs2)
723{
724
725 G4double value = 0.;
726
727 if((xs1!=0)&&(e1!=0)){
728 // Log-log interpolation by default
729 G4double a = (std::log10(xs2)-std::log10(xs1))
730 / (std::log10(e2)-std::log10(e1));
731 G4double b = std::log10(xs2) - a*std::log10(e2);
732 G4double sigma = a*std::log10(e) + b;
733 value = (std::pow(10.,sigma));
734 }
735 else{
736 // Lin-Lin interpolation
737 G4double d1 = xs1;
738 G4double d2 = xs2;
739 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
740 }
741
742 return value;
743}
744
745//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
746
747G4double G4DNARelativisticIonisationModel::QuadInterpolator(
748 G4double e11, G4double e12,
749 G4double e21, G4double e22,
750 G4double xs11, G4double xs12,
751 G4double xs21, G4double xs22,
752 G4double t1, G4double t2,
753 G4double t, G4double e)
754{
755 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
756 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
757 G4double value
758 = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
759 return value;
760}
761
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition: G4FindDataDir.cc:93
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
G4DNAMolecularMaterial::GetNumMolPerVolTableFor
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
Definition: G4DNAMolecularMaterial.cc:418
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:80
G4DNARelativisticIonisationModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
Definition: G4DNARelativisticIonisationModel.cc:92
G4DNARelativisticIonisationModel::LoadAtomicStates
virtual void LoadAtomicStates(G4int z, const char *path)
Definition: G4DNARelativisticIonisationModel.cc:394
G4DNARelativisticIonisationModel::GetTotalCrossSection
virtual G4double GetTotalCrossSection(const G4Material *material, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4DNARelativisticIonisationModel.cc:463
G4DNARelativisticIonisationModel::G4DNARelativisticIonisationModel
G4DNARelativisticIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARelativisticIonisationModel")
Definition: G4DNARelativisticIonisationModel.cc:61
G4DNARelativisticIonisationModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNARelativisticIonisationModel.cc:232
G4DNARelativisticIonisationModel::GetPartialCrossSection
virtual G4double GetPartialCrossSection(const G4Material *material, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4DNARelativisticIonisationModel.cc:482
G4DNARelativisticIonisationModel::GetDifferentialCrossSection
virtual G4double GetDifferentialCrossSection(const G4Material *material, const G4ParticleDefinition *particle, G4double kineticEnergy, G4double secondaryEnergy, G4int level)
Definition: G4DNARelativisticIonisationModel.cc:523
G4DNARelativisticIonisationModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNARelativisticIonisationModel.cc:280
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:88
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:389
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4Material::GetZ
G4double GetZ() const
Definition: G4Material.cc:745
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
G4Material::GetAtomicNumDensityVector
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:211
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:691
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:106
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
Definition: G4ParticleChangeForGamma.hh:120
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:246
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.cc:234
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:802
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
N
#define N
Definition: crc32.c:56
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158