Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNACPA100IonisationModel.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// CPA100 ionisation model class for electrons
27//
28// Based on the work of M. Terrissol and M. C. Bordage
29//
30// Users are requested to cite the following papers:
31// - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32// - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33// M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34//
35// Authors of this class:
36// M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37//
38// 15.01.2014: creation
39//
40
41#include "G4DNACPA100IonisationModel.hh"
42#include "G4PhysicalConstants.hh"
43#include "G4SystemOfUnits.hh"
44#include "G4UAtomicDeexcitation.hh"
45#include "G4LossTableManager.hh"
46#include "G4DNAChemistryManager.hh"
47#include "G4DNAMolecularMaterial.hh"
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51using namespace std;
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55G4DNACPA100IonisationModel::G4DNACPA100IonisationModel(const G4ParticleDefinition*,
56 const G4String& nam)
57:G4VEmModel(nam),isInitialised(false)
58{
59 verboseLevel= 0;
60 // Verbosity scale:
61 // 0 = nothing
62 // 1 = warning for energy non-conservation
63 // 2 = details of energy budget
64 // 3 = calculation of cross sections, file openings, sampling of atoms
65 // 4 = entering in methods
66
67 if( verboseLevel>0 )
68 {
69 G4cout << "CPA100 ionisation model is constructed " << G4endl;
70 }
71
72 SetLowEnergyLimit(11*eV);
73 SetHighEnergyLimit(255955*eV);
74
75 // Mark this model as "applicable" for atomic deexcitation
76 SetDeexcitationFlag(true);
77 fAtomDeexcitation = 0;
78 fParticleChangeForGamma = 0;
79 fpMolWaterDensity = 0;
80
81 // Selection of computation method
82
83 // useDcs = true if usage of dcs for sampling of secondaries
84 // useDcs = false if usage of composition sampling (DEFAULT)
85
86 useDcs = true;
87
88 // if useDcs is true, one has the following choice
89 // fasterCode = true for usage of cumulated dcs (DEFAULT)
90 // fasterCode = false for usage of non-cumulated dcs
91
92 fasterCode = true;
93
94 // Selection of stationary mode
95
96 statCode = false;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101G4DNACPA100IonisationModel::~G4DNACPA100IonisationModel()
102{
103 // Cross section
104
105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
107 {
108 G4DNACrossSectionDataSet* table = pos->second;
109 delete table;
110 }
111
112 // Final state
113
114 eVecm.clear();
115
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
120void G4DNACPA100IonisationModel::Initialise(const G4ParticleDefinition* particle,
121 const G4DataVector& /*cuts*/)
122{
123
124 if (verboseLevel > 3)
125 G4cout << "Calling G4DNACPA100IonisationModel::Initialise()" << G4endl;
126
127 // Energy limits
128
129 // The following file is proved by M. Terrissol et al. (sigion3)
130
131 G4String fileElectron("dna/sigma_ionisation_e_cpa100_form_rel");
132
133 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
134
135 G4String electron;
136
137 G4double scaleFactor = 1.e-20 * m*m;
138
139 const char *path = G4FindDataDir("G4LEDATA");
140
141 // *** ELECTRON
142
143 electron = electronDef->GetParticleName();
144
145 tableFile[electron] = fileElectron;
146
147 // Cross section
148
149 G4DNACrossSectionDataSet* tableE =
150 new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
151
152 //G4DNACrossSectionDataSet* tableE =
153 // new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation, eV,scaleFactor );
154
155 tableE->LoadData(fileElectron);
156
157 tableData[electron] = tableE;
158
159 // Final state
160
161 // ******************************
162
163 if (useDcs)
164 {
165
166 std::ostringstream eFullFileName;
167
168 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat";
169
170 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_cpa100_rel.dat";
171
172 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
173
174 if (!eDiffCrossSection)
175 {
176 if (fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
177 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat");
178
179 if (!fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
180 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_cpa100_rel.dat");
181 }
182
183 // Clear the arrays for re-initialization case (MT mode)
184 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
185
186 eTdummyVec.clear();
187 eVecm.clear();
188 eProbaShellMap->clear();
189 eDiffCrossSectionData->clear();
190 eNrjTransfData->clear();
191
192 //
193
194 eTdummyVec.push_back(0.);
195 while(!eDiffCrossSection.eof())
196 {
197 G4double tDummy;
198 G4double eDummy;
199 eDiffCrossSection>>tDummy>>eDummy;
200 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
201 for (G4int j=0; j<5; j++)
202 {
203 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
204
205 if (fasterCode)
206 {
207 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
208 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
209 }
210
211 // SI - only if eof is not reached
212 if (!eDiffCrossSection.eof() && !fasterCode)
213 eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
214
215 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
216
217 }
218 }
219
220 //
221
222 } // end of if (useDcs)
223
224 // ******************************
225
226 //
227
228 if( verboseLevel>0 )
229 {
230 G4cout << "CPA100 ionisation model is initialized " << G4endl
231 << "Energy range: "
232 << LowEnergyLimit() / eV << " eV - "
233 << HighEnergyLimit() / keV << " keV for "
234 << particle->GetParticleName()
235 << G4endl;
236 }
237
238 // Initialize water density pointer
239 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
240
241 // AD
242 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
243
244 if (isInitialised) return;
245 fParticleChangeForGamma = GetParticleChangeForGamma();
246 isInitialised = true;
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250
251G4double G4DNACPA100IonisationModel::CrossSectionPerVolume( const G4Material* material,
252 const G4ParticleDefinition* particleDefinition,
253 G4double ekin,
254 G4double,
255 G4double)
256{
257
258 if (verboseLevel > 3)
259 G4cout << "Calling CrossSectionPerVolume() of G4DNACPA100IonisationModel" << G4endl;
260
261 if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
262
263 // Calculate total cross section for model
264
265 G4double sigma=0;
266
267 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
268
269 const G4String& particleName = particleDefinition->GetParticleName();
270
271 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
272 {
273 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
274 pos = tableData.find(particleName);
275
276 if (pos != tableData.end())
277 {
278 G4DNACrossSectionDataSet* table = pos->second;
279 if (table != 0) sigma = table->FindValue(ekin);
280 }
281 else
282 {
283 G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume","em0002",
284 FatalException,"Model not applicable to particle type.");
285 }
286 }
287
288 if (verboseLevel > 2)
289 {
290 G4cout << "__________________________________" << G4endl;
291 G4cout << "G4DNACPA100IonisationModel - XS INFO START" << G4endl;
292 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
293 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
294 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
295 G4cout << "G4DNACPA100IonisationModel - XS INFO END" << G4endl;
296 }
297
298 return sigma*waterDensity;
299}
300
301//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302
303void G4DNACPA100IonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
304 const G4MaterialCutsCouple* ,//must be set!
305 const G4DynamicParticle* particle,
306 G4double,
307 G4double)
308{
309 if (verboseLevel > 3)
310 G4cout << "Calling SampleSecondaries() of G4DNACPA100IonisationModel" << G4endl;
311
312 G4double k = particle->GetKineticEnergy();
313
314 const G4String& particleName = particle->GetDefinition()->GetParticleName();
315
316 if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
317 {
318 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
319 G4double particleMass = particle->GetDefinition()->GetPDGMass();
320 G4double totalEnergy = k + particleMass;
321 G4double pSquare = k * (totalEnergy + particleMass);
322 G4double totalMomentum = std::sqrt(pSquare);
323
324 G4int ionizationShell = -1;
325
326 ionizationShell = RandomSelect(k,particleName);
327
328 //SI: PROTECTION FOR G4LOGLOGINTERPOLATION ON UPPER VALUE
329 if (k<waterStructure.IonisationEnergy(ionizationShell)) { return; }
330
331 G4double bindingEnergy = 0;
332 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
333
334 G4double secondaryKinetic=-1000*eV;
335
336 if (useDcs && !fasterCode)
337 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
338
339 if (useDcs && fasterCode)
340 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
341
342 if (!useDcs)
343 secondaryKinetic = RandomizeEjectedElectronEnergyFromCompositionSampling(particle->GetDefinition(),k,ionizationShell);
344
345 // Quick test
346 /*
347 FILE* myFile;
348 myFile=fopen("nrj.txt","a");
349 fprintf(myFile,"%e\n", secondaryKinetic/eV );
350 fclose(myFile);
351 */
352
353 G4double cosTheta = 0.;
354 G4double phi = 0.;
355 RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
356
357 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
358 G4double dirX = sinTheta*std::cos(phi);
359 G4double dirY = sinTheta*std::sin(phi);
360 G4double dirZ = cosTheta;
361 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
362 deltaDirection.rotateUz(primaryDirection);
363
364 // SI - For atom. deexc. tagging - 23/05/2017
365 if (secondaryKinetic>0)
366 {
367 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
368 fvect->push_back(dp);
369 }
370 //
371
372 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
373 {
374 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
375
376 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
377 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
378 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
379 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
380 finalPx /= finalMomentum;
381 finalPy /= finalMomentum;
382 finalPz /= finalMomentum;
383
384 G4ThreeVector direction;
385 direction.set(finalPx,finalPy,finalPz);
386
387 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
388 }
389
390 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
391
392 // SI - For atom. deexc. tagging - 23/05/2017
393
394 // AM: sample deexcitation
395 // here we assume that H_{2}O electronic levels are the same of Oxigen.
396 // this can be considered true with a rough 10% error in energy on K-shell,
397
398 std::size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
399 std::size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
400
401 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
402
403 // SI: only atomic deexcitation from K shell is considered
404 if(fAtomDeexcitation && ionizationShell == 4)
405 {
406 G4int Z = 8;
407 const G4AtomicShell* shell =
408 fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
409 secNumberInit = fvect->size();
410 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
411 secNumberFinal = fvect->size();
412
413 if(secNumberFinal > secNumberInit)
414 {
415 for (std::size_t i=secNumberInit; i<secNumberFinal; ++i)
416 {
417 //Check if there is enough residual energy
418 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
419 {
420 //Ok, this is a valid secondary: keep it
421 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
422 }
423 else
424 {
425 //Invalid secondary: not enough energy to create it!
426 //Keep its energy in the local deposit
427 delete (*fvect)[i];
428 (*fvect)[i]=0;
429 }
430 }
431 }
432
433 }
434
435 //This should never happen
436 if(bindingEnergy < 0.0)
437 G4Exception("G4DNACPA100IonisatioModel1::SampleSecondaries()",
438 "em2050",FatalException,"Negative local energy deposit");
439
440 //bindingEnergy has been decreased
441 //by the amount of energy taken away by deexc. products
442
443 if (!statCode)
444 {
445 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
446 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
447 }
448 else
449 {
450 fParticleChangeForGamma->SetProposedKineticEnergy(k);
451 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
452 }
453
454 // TEST //////////////////////////
455 // if (secondaryKinetic<0) abort();
456 // if (scatteredEnergy<0) abort();
457 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
458 // if (k-scatteredEnergy<0) abort();
459 /////////////////////////////////
460
461 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
462 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
463 ionizationShell,
464 theIncomingTrack);
465 }
466
467}
468
469//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
470
471G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
472 G4double k, G4int shell)
473{
474 // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
475
476 if (particleDefinition == G4Electron::ElectronDefinition())
477 {
478 G4double maximumEnergyTransfer=0.;
479 if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
480 else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
481
482 // SI : original method
483 /*
484 G4double crossSectionMaximum = 0.;
485 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
486 {
487 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
488 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
489 }
490 */
491
492 // SI : alternative method
493
494 G4double crossSectionMaximum = 0.;
495
496 G4double minEnergy = waterStructure.IonisationEnergy(shell);
497 G4double maxEnergy = maximumEnergyTransfer;
498
499 // nEnergySteps can be optimized - 100 by default
500 G4int nEnergySteps = 50;
501
502 // *** METHOD 1
503 // FOR SLOW COMPUTATION ONLY
504 /*
505 G4double value(minEnergy);
506 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
507 G4int step(nEnergySteps);
508 while (step>0)
509 {
510 step--;
511 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
512 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
513 value*=stpEnergy;
514 }
515 */
516
517 // *** METHOD 2 : Faster method for CPA100 only since DCS is monotonously decreasing
518 // FOR SLOW COMPUTATION ONLY
519
520 G4double value(minEnergy);
521 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
522 G4int step(nEnergySteps);
523 G4double differentialCrossSection = 0.;
524 while (step>0)
525 {
526 step--;
527 differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
528 if(differentialCrossSection >0)
529 {
530 crossSectionMaximum=differentialCrossSection;
531 break;
532 }
533 value*=stpEnergy;
534 }
535
536 //
537
538 G4double secondaryElectronKineticEnergy=0.;
539 do
540 {
541 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
542 } while(G4UniformRand()*crossSectionMaximum >
543 DifferentialCrossSection(particleDefinition, k/eV,
544 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
545
546 return secondaryElectronKineticEnergy;
547
548 }
549
550 return 0;
551}
552
553//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
554
555void G4DNACPA100IonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition*,
556 G4double k,
557 G4double secKinetic,
558 G4double & cosTheta,
559 G4double & phi )
560{
561
562 phi = twopi * G4UniformRand();
563 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
564 cosTheta = std::sqrt(1.-sin2O);
565
566}
567
568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
569
570G4double G4DNACPA100IonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
571 G4double k,
572 G4double energyTransfer,
573 G4int ionizationLevelIndex)
574{
575 G4double sigma = 0.;
576
577 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
578 {
579 G4double valueT1 = 0;
580 G4double valueT2 = 0;
581 G4double valueE21 = 0;
582 G4double valueE22 = 0;
583 G4double valueE12 = 0;
584 G4double valueE11 = 0;
585
586 G4double xs11 = 0;
587 G4double xs12 = 0;
588 G4double xs21 = 0;
589 G4double xs22 = 0;
590
591 if (particleDefinition == G4Electron::ElectronDefinition())
592 {
593 // Protection against out of boundary access
594 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
595 //
596
597 // k should be in eV and energy transfer eV also
598
599 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
600
601 std::vector<G4double>::iterator t1 = t2-1;
602
603 // SI : the following condition avoids situations where energyTransfer >last vector element
604
605 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
606 {
607 std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
608 std::vector<G4double>::iterator e11 = e12-1;
609
610 std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
611 std::vector<G4double>::iterator e21 = e22-1;
612
613 valueT1 =*t1;
614 valueT2 =*t2;
615 valueE21 =*e21;
616 valueE22 =*e22;
617 valueE12 =*e12;
618 valueE11 =*e11;
619
620 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
621 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
622 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
623 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
624
625 }
626
627 }
628
629 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
630 if (xsProduct != 0.)
631 {
632 sigma = QuadInterpolator( valueE11, valueE12,
633 valueE21, valueE22,
634 xs11, xs12,
635 xs21, xs22,
636 valueT1, valueT2,
637 k, energyTransfer);
638 }
639
640 }
641 return sigma;
642}
643
644//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
645
646G4double G4DNACPA100IonisationModel::Interpolate( G4double e1,
647 G4double e2,
648 G4double e,
649 G4double xs1,
650 G4double xs2)
651{
652
653 G4double value = 0.;
654
655 // Log-log interpolation by default
656
657 if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !fasterCode && useDcs)
658 {
659 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
660 G4double b = std::log10(xs2) - a*std::log10(e2);
661 G4double sigma = a*std::log10(e) + b;
662 value = (std::pow(10.,sigma));
663 }
664
665 // Switch to lin-lin interpolation
666 /*
667 if ((e2-e1)!=0)
668 {
669 G4double d1 = xs1;
670 G4double d2 = xs2;
671 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
672 }
673 */
674
675 // Switch to log-lin interpolation for faster code
676
677 if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 && fasterCode && useDcs )
678 {
679 G4double d1 = std::log10(xs1);
680 G4double d2 = std::log10(xs2);
681 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
682 }
683
684 // Switch to lin-lin interpolation for faster code
685 // in case one of xs1 or xs2 (=cum proba) value is zero
686
687 if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) && fasterCode && useDcs )
688 {
689 G4double d1 = xs1;
690 G4double d2 = xs2;
691 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
692 }
693
694 /*
695 G4cout
696 << e1 << " "
697 << e2 << " "
698 << e << " "
699 << xs1 << " "
700 << xs2 << " "
701 << value
702 << G4endl;
703 */
704
705 return value;
706}
707
708//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
709
710G4double G4DNACPA100IonisationModel::QuadInterpolator(G4double e11, G4double e12,
711 G4double e21, G4double e22,
712 G4double xs11, G4double xs12,
713 G4double xs21, G4double xs22,
714 G4double t1, G4double t2,
715 G4double t, G4double e)
716{
717 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
718 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
719 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
720
721 return value;
722}
723
724//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
725
726G4int G4DNACPA100IonisationModel::RandomSelect(G4double k, const G4String& particle )
727{
728 G4int level = 0;
729
730 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
731 pos = tableData.find(particle);
732
733 if (pos != tableData.end())
734 {
735 G4DNACrossSectionDataSet* table = pos->second;
736
737 if (table != 0)
738 {
739 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
740 const G4int n = (G4int)table->NumberOfComponents();
741 G4int i(n);
742 G4double value = 0.;
743
744 //Verification
745 /*
746 G4double tmp=200*keV;
747 G4cout << table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) << G4endl;
748 G4cout << table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) << G4endl;
749 G4cout << table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) << G4endl;
750 G4cout << table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m) << G4endl;
751 G4cout << table->GetComponent(4)->FindValue(tmp)/(1e-20*m*m) << G4endl;
752 G4cout <<
753 table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) +
754 table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) +
755 table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) +
756 table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m)
757 << G4endl;
758 abort();
759 */
760 //
761 //Dump
762 //
763 /*
764 G4double minEnergy = 10.985 * eV;
765 G4double maxEnergy = 255955. * eV;
766 G4int nEnergySteps = 1000;
767 G4double energy(minEnergy);
768 G4double stpEnergy(std::pow(maxEnergy/energy, 1./static_cast<G4double>(nEnergySteps-1)));
769 G4int step(nEnergySteps);
770 system ("rm -rf ionisation-cpa100.out");
771 FILE* myFile=fopen("ionisation-cpa100.out","a");
772 while (step>0)
773 {
774 step--;
775 fprintf (myFile,"%16.9le %16.9le %16.9le %16.9le %16.9le %16.9le %16.9le \n",
776 energy/eV,
777 table->GetComponent(0)->FindValue(energy)/(1e-20*m*m),
778 table->GetComponent(1)->FindValue(energy)/(1e-20*m*m),
779 table->GetComponent(2)->FindValue(energy)/(1e-20*m*m),
780 table->GetComponent(3)->FindValue(energy)/(1e-20*m*m),
781 table->GetComponent(4)->FindValue(energy)/(1e-20*m*m),
782 table->GetComponent(0)->FindValue(energy)/(1e-20*m*m)+
783 table->GetComponent(1)->FindValue(energy)/(1e-20*m*m)+
784 table->GetComponent(2)->FindValue(energy)/(1e-20*m*m)+
785 table->GetComponent(3)->FindValue(energy)/(1e-20*m*m)+
786 table->GetComponent(4)->FindValue(energy)/(1e-20*m*m)
787 );
788 energy*=stpEnergy;
789 }
790 fclose (myFile);
791 abort();
792 */
793 //
794 // end of dump
795 //
796 // Test of diff XS
797 // G4double nrj1 = .26827E+04; // in eV
798 // G4double nrj2 = .57991E+03; // in eV
799 // Shells run from 0 to 4
800 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 0)/(1e-20*m*m) << G4endl;
801 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 1)/(1e-20*m*m) << G4endl;
802 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 2)/(1e-20*m*m) << G4endl;
803 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 3)/(1e-20*m*m) << G4endl;
804 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 4)/(1e-20*m*m) << G4endl;
805 // abort();
806 //
807
808 while (i>0)
809 {
810 i--;
811 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
812 value += valuesBuffer[i];
813 }
814
815 value *= G4UniformRand();
816
817 i = n;
818
819 while (i > 0)
820 {
821 i--;
822 if (valuesBuffer[i] > value)
823 {
824 delete[] valuesBuffer;
825
826 return i;
827 }
828 value -= valuesBuffer[i];
829 }
830
831 if (valuesBuffer) delete[] valuesBuffer;
832
833 }
834 }
835 else
836 {
837 G4Exception("G4DNACPA100IonisationModel::RandomSelect","em0002",
838 FatalException,"Model not applicable to particle type.");
839 }
840
841 return level;
842}
843
844//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
845
846G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs
847(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
848{
849 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
850
851 G4double secondaryElectronKineticEnergy = 0.;
852
853 secondaryElectronKineticEnergy=
854 RandomTransferedEnergy(particleDefinition, k/eV, shell)*eV-waterStructure.IonisationEnergy(shell);
855
856 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
857 if (secondaryElectronKineticEnergy<0.) return 0.;
858 //
859
860 return secondaryElectronKineticEnergy;
861}
862
863//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
864
865G4double G4DNACPA100IonisationModel::RandomTransferedEnergy
866(G4ParticleDefinition* particleDefinition,G4double k, G4int ionizationLevelIndex)
867{
868
869 G4double random = G4UniformRand();
870
871 G4double nrj = 0.;
872
873 G4double valueK1 = 0;
874 G4double valueK2 = 0;
875 G4double valuePROB21 = 0;
876 G4double valuePROB22 = 0;
877 G4double valuePROB12 = 0;
878 G4double valuePROB11 = 0;
879
880 G4double nrjTransf11 = 0;
881 G4double nrjTransf12 = 0;
882 G4double nrjTransf21 = 0;
883 G4double nrjTransf22 = 0;
884
885 if (particleDefinition == G4Electron::ElectronDefinition())
886 {
887
888 // Protection against out of boundary access
889 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
890 //
891
892 // k should be in eV
893
894 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
895
896 std::vector<G4double>::iterator k1 = k2-1;
897
898 /*
899 G4cout << "----> k=" << k
900 << " " << *k1
901 << " " << *k2
902 << " " << random
903 << " " << ionizationLevelIndex
904 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
905 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
906 << G4endl;
907 */
908
909 // SI : the following condition avoids situations where random >last vector element
910
911 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
912 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
913
914 {
915
916 std::vector<G4double>::iterator prob12 =
917 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
918 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
919
920 std::vector<G4double>::iterator prob11 = prob12-1;
921
922
923 std::vector<G4double>::iterator prob22 =
924 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
925 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
926
927 std::vector<G4double>::iterator prob21 = prob22-1;
928
929 valueK1 =*k1;
930 valueK2 =*k2;
931 valuePROB21 =*prob21;
932 valuePROB22 =*prob22;
933 valuePROB12 =*prob12;
934 valuePROB11 =*prob11;
935
936
937 /*
938 G4cout << " " << random << " " << valuePROB11 << " "
939 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
940 */
941
942 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
943 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
944 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
945 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
946
947 /*
948 G4cout << " " << ionizationLevelIndex << " "
949 << random << " " <<valueK1 << " " << valueK2 << G4endl;
950
951 G4cout << " " << random << " " << nrjTransf11 << " "
952 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
953 */
954
955 // Avoids cases where cum xs is 1 in last bin only and zero everywhere else
956
957 if (valuePROB11==0. && valuePROB12==1.)
958 {
959 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
960 G4double valueNrjTransf = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
961 return valueNrjTransf;
962 }
963
964 }
965
966 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
967
968 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
969 {
970
971 std::vector<G4double>::iterator prob22 =
972
973 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
974 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
975
976 std::vector<G4double>::iterator prob21 = prob22-1;
977
978 valueK1 =*k1;
979 valueK2 =*k2;
980 valuePROB21 =*prob21;
981 valuePROB22 =*prob22;
982
983 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
984
985 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
986 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
987
988 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
989
990 // zero is explicitly set
991
992 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
993
994 /*
995 G4cout << " " << ionizationLevelIndex << " "
996 << random << " " <<valueK1 << " " << valueK2 << G4endl;
997
998 G4cout << " " << random << " " << nrjTransf11 << " "
999 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1000
1001 G4cout << "ici" << " " << value << G4endl;
1002 */
1003
1004 return value;
1005 }
1006
1007 }
1008
1009 // End electron case
1010
1011 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1012
1013 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1014
1015 if (nrjTransfProduct != 0.)
1016 {
1017 nrj = QuadInterpolator( valuePROB11, valuePROB12,
1018 valuePROB21, valuePROB22,
1019 nrjTransf11, nrjTransf12,
1020 nrjTransf21, nrjTransf22,
1021 valueK1, valueK2,
1022 k, random);
1023 }
1024
1025 //G4cout << nrj << endl;
1026
1027 return nrj ;
1028}
1029
1030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1031
1032G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCompositionSampling
1033(G4ParticleDefinition*, G4double tt, G4int shell)
1034{
1035 //G4cout << "*** Rejection method for " << " " << particleDefinition->GetParticleName() << G4endl;
1036
1037 // ***** METHOD 1 ***** (sequential)
1038 /*
1039
1040 // ww is KINETIC ENERGY OF SECONDARY ELECTRON
1041 G4double un=1.;
1042 G4double deux=2.;
1043
1044 G4double bb = waterStructure.IonisationEnergy(shell);
1045 G4double uu = waterStructure.UEnergy(shell);
1046
1047 if (tt<=bb) return 0.;
1048
1049 G4double t = tt/bb;
1050 G4double u = uu/bb;
1051 G4double tp1 = t + un;
1052 G4double tu1 = t + u + un;
1053 G4double tm1 = t - un;
1054 G4double tp12 = tp1 * tp1;
1055 G4double dlt = std::log(t);
1056
1057 G4double a1 = t * tm1 / tu1 / tp12;
1058 G4double a2 = tm1 / tu1 / t / tp1 / deux;
1059 G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
1060 G4double ato = a1 + a2 + a3;
1061
1062 // 15
1063
1064 G4double r1 =G4UniformRand();
1065 G4double r2 =G4UniformRand();
1066 G4double r3 =G4UniformRand();
1067
1068 while (r1<=a1/ato)
1069 {
1070 G4double fx1=r2*tm1/tp1;
1071 G4double wx1=un/(un-fx1)-un;
1072 G4double gx1=(t-wx1)/t;
1073 if(r3 <= gx1) return wx1*bb;
1074
1075 r1 =G4UniformRand();
1076 r2 =G4UniformRand();
1077 r3 =G4UniformRand();
1078
1079 }
1080
1081 // 20
1082
1083 while (r1<=(a1+a2)/ato)
1084 {
1085 G4double fx2=tp1+r2*tm1;
1086 G4double wx2=t-t*tp1/fx2;
1087 G4double gx2=deux*(un-(t-wx2)/tp1);
1088 if(r3 <= gx2) return wx2*bb;
1089
1090 // REPEAT 15
1091 r1 =G4UniformRand();
1092 r2 =G4UniformRand();
1093 r3 =G4UniformRand();
1094
1095 while (r1<=a1/ato)
1096 {
1097 G4double fx1=r2*tm1/tp1;
1098 G4double wx1=un/(un-fx1)-un;
1099 G4double gx1=(t-wx1)/t;
1100 if(r3 <= gx1) return wx1*bb;
1101 r1 =G4UniformRand();
1102 r2 =G4UniformRand();
1103 r3 =G4UniformRand();
1104 }
1105 // END 15
1106
1107 }
1108
1109 // 30
1110
1111 G4double wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
1112 G4double gg3=(wx3+un)/(t-wx3);
1113 G4double gx3=(un+gg3*gg3*gg3)/deux;
1114
1115 while (r3>gx3)
1116 {
1117
1118 // 15
1119
1120 r1 =G4UniformRand();
1121 r2 =G4UniformRand();
1122 r3 =G4UniformRand();
1123
1124 while (r1<=a1/ato)
1125 {
1126 G4double fx1=r2*tm1/tp1;
1127 G4double wx1=un/(un-fx1)-un;
1128 G4double gx1=(t-wx1)/t;
1129 if(r3 <= gx1) return wx1*bb;
1130
1131 r1 =G4UniformRand();
1132 r2 =G4UniformRand();
1133 r3 =G4UniformRand();
1134
1135 }
1136
1137 // 20
1138
1139 while (r1<=(a1+a2)/ato)
1140 {
1141 G4double fx2=tp1+r2*tm1;
1142 G4double wx2=t-t*tp1/fx2;
1143 G4double gx2=deux*(un-(t-wx2)/tp1);
1144 if(r3 <= gx2)return wx2*bb;
1145
1146 // REPEAT 15
1147 r1 =G4UniformRand();
1148 r2 =G4UniformRand();
1149 r3 =G4UniformRand();
1150
1151 while (r1<=a1/ato)
1152 {
1153 G4double fx1=r2*tm1/tp1;
1154 G4double wx1=un/(un-fx1)-un;
1155 G4double gx1=(t-wx1)/t;
1156 if(r3 <= gx1) return wx1*bb;
1157
1158 r1 =G4UniformRand();
1159 r2 =G4UniformRand();
1160 r3 =G4UniformRand();
1161 }
1162 //
1163
1164 }
1165
1166 wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
1167 gg3=(wx3+un)/(t-wx3);
1168 gx3=(un+gg3*gg3*gg3)/deux;
1169
1170 }
1171
1172 //
1173
1174 return wx3*bb;
1175 */
1176
1177 // ***** METHOD by M. C. Bordage ***** (optimized)
1178
1179 G4double un=1.;
1180 G4double deux=2.;
1181
1182 G4double bb = waterStructure.IonisationEnergy(shell);
1183 G4double uu = waterStructure.UEnergy(shell);
1184
1185 if (tt<=bb) return 0.;
1186
1187 G4double t = tt/bb;
1188 G4double u = uu/bb;
1189 G4double tp1 = t + un;
1190 G4double tu1 = t + u + un;
1191 G4double tm1 = t - un;
1192 G4double tp12 = tp1 * tp1;
1193 G4double dlt = std::log(t);
1194
1195 G4double a1 = t * tm1 / tu1 / tp12;
1196 G4double a2 = tm1 / tu1 / t / tp1 / deux;
1197 G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
1198 G4double ato = a1 + a2 + a3;
1199
1200 G4double A1 = a1/ato;
1201 G4double A2 = (a1+a2)/ato;
1202 G4int F = 0;
1203 G4double fx=0;
1204 G4double gx=0;
1205 G4double gg=0;
1206 G4double wx=0;
1207
1208 G4double r1=0;
1209 G4double r2=0;
1210 G4double r3=0;
1211
1212 //
1213
1214 do
1215 {
1216 r1 =G4UniformRand();
1217 r2 =G4UniformRand();
1218 r3 =G4UniformRand();
1219
1220 if (r1>A2)
1221 F=3;
1222 else if ((r1>A1) && (r1< A2))
1223 F=2;
1224 else
1225 F=1;
1226
1227 switch (F)
1228 {
1229 case 1:
1230 {
1231 fx=r2*tm1/tp1;
1232 wx=un/(un-fx)-un;
1233 gx=(t-wx)/t;
1234 break;
1235 }
1236
1237 case 2:
1238 {
1239 fx=tp1+r2*tm1;
1240 wx=t-t*tp1/fx;
1241 gx=deux*(un-(t-wx)/tp1);
1242 break;
1243 }
1244
1245 case 3:
1246 {
1247 fx=un-r2*(tp12-deux*deux)/tp12;
1248 wx=sqrt(un/fx)-un;
1249 gg=(wx+un)/(t-wx);
1250 gx=(un+gg*gg*gg)/deux;
1251 break;
1252 }
1253 } // switch
1254
1255 } while (r3>gx);
1256
1257 return wx*bb;
1258
1259}
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
eIonizedMolecule
@ eIonizedMolecule
Definition: G4DNAChemistryManager.hh:70
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
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::set
void set(double x, double y, double z)
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
G4DNACPA100IonisationModel::DifferentialCrossSection
G4double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
Definition: G4DNACPA100IonisationModel.cc:570
G4DNACPA100IonisationModel::G4DNACPA100IonisationModel
G4DNACPA100IonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNACPA100IonisationModel")
Definition: G4DNACPA100IonisationModel.cc:55
G4DNACPA100IonisationModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNACPA100IonisationModel.cc:251
G4DNACPA100IonisationModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNACPA100IonisationModel.hh:94
G4DNACPA100IonisationModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
Definition: G4DNACPA100IonisationModel.cc:120
G4DNACPA100IonisationModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNACPA100IonisationModel.cc:303
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:132
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:607
G4DNACrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4DNACrossSectionDataSet.cc:483
G4DNACrossSectionDataSet::NumberOfComponents
virtual size_t NumberOfComponents(void) const
Definition: G4DNACrossSectionDataSet.hh:74
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:68
G4DNACrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4DNACrossSectionDataSet.cc:83
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
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
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:255
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
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
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:746
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::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:802
G4VParticleChange::GetCurrentTrack
const G4Track * GetCurrentTrack() const
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)