Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNABornIonisationModel.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$
27//
28
29#include "G4DNABornIonisationModel.hh"
30#include "G4PhysicalConstants.hh"
31#include "G4SystemOfUnits.hh"
32#include "G4UAtomicDeexcitation.hh"
33#include "G4LossTableManager.hh"
34#include "G4DNAChemistryManager.hh"
35#include "G4DNAMolecularMaterial.hh"
36
37//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38
39using namespace std;
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43G4DNABornIonisationModel::G4DNABornIonisationModel(const G4ParticleDefinition*,
44 const G4String& nam)
45 :G4VEmModel(nam),isInitialised(false)
46{
47 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
48
49 verboseLevel= 0;
50 // Verbosity scale:
51 // 0 = nothing
52 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file openings, sampling of atoms
55 // 4 = entering in methods
56
57 if( verboseLevel>0 )
58 {
59 G4cout << "Born ionisation model is constructed " << G4endl;
60 }
61
62 //Mark this model as "applicable" for atomic deexcitation
63 SetDeexcitationFlag(true);
64 fAtomDeexcitation = 0;
65 fParticleChangeForGamma = 0;
66 fpMolWaterDensity = 0;
67}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71G4DNABornIonisationModel::~G4DNABornIonisationModel()
72{
73 // Cross section
74
75 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
76 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
77 {
78 G4DNACrossSectionDataSet* table = pos->second;
79 delete table;
80 }
81
82 // Final state
83
84 eVecm.clear();
85 pVecm.clear();
86
87}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91void G4DNABornIonisationModel::Initialise(const G4ParticleDefinition* particle,
92 const G4DataVector& /*cuts*/)
93{
94
95 if (verboseLevel > 3)
96 G4cout << "Calling G4DNABornIonisationModel::Initialise()" << G4endl;
97
98 // Energy limits
99
100 G4String fileElectron("dna/sigma_ionisation_e_born");
101 G4String fileProton("dna/sigma_ionisation_p_born");
102
103 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
104 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
105
106 G4String electron;
107 G4String proton;
108
109 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
110
111 char *path = getenv("G4LEDATA");
112
113 // *** ELECTRON
114
115 electron = electronDef->GetParticleName();
116
117 tableFile[electron] = fileElectron;
118
119 lowEnergyLimit[electron] = 11. * eV;
120 highEnergyLimit[electron] = 1. * MeV;
121
122 // Cross section
123
124 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
125 tableE->LoadData(fileElectron);
126
127 tableData[electron] = tableE;
128
129 // Final state
130
131 std::ostringstream eFullFileName;
132 eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
133 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
134
135 if (!eDiffCrossSection)
136 {
137 G4Exception("G4DNABornIonisationModel::Initialise","em0003",
138 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
139 }
140
141 eTdummyVec.push_back(0.);
142 while(!eDiffCrossSection.eof())
143 {
144 double tDummy;
145 double eDummy;
146 eDiffCrossSection>>tDummy>>eDummy;
147 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
148 for (int j=0; j<5; j++)
149 {
150 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
151
152 // SI - only if eof is not reached !
153 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
154
155 eVecm[tDummy].push_back(eDummy);
156
157 }
158 }
159
160 // *** PROTON
161
162 proton = protonDef->GetParticleName();
163
164 tableFile[proton] = fileProton;
165
166 lowEnergyLimit[proton] = 500. * keV;
167 highEnergyLimit[proton] = 100. * MeV;
168
169 // Cross section
170
171 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
172 tableP->LoadData(fileProton);
173
174 tableData[proton] = tableP;
175
176 // Final state
177
178 std::ostringstream pFullFileName;
179 pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
180 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
181
182 if (!pDiffCrossSection)
183 {
184 G4Exception("G4DNABornIonisationModel::Initialise","em0003",
185 FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
186 }
187
188 pTdummyVec.push_back(0.);
189 while(!pDiffCrossSection.eof())
190 {
191 double tDummy;
192 double eDummy;
193 pDiffCrossSection>>tDummy>>eDummy;
194 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
195 for (int j=0; j<5; j++)
196 {
197 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
198
199 // SI - only if eof is not reached !
200 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
201
202 pVecm[tDummy].push_back(eDummy);
203 }
204 }
205
206 //
207
208 if (particle==electronDef)
209 {
210 SetLowEnergyLimit(lowEnergyLimit[electron]);
211 SetHighEnergyLimit(highEnergyLimit[electron]);
212 }
213
214 if (particle==protonDef)
215 {
216 SetLowEnergyLimit(lowEnergyLimit[proton]);
217 SetHighEnergyLimit(highEnergyLimit[proton]);
218 }
219
220 if( verboseLevel>0 )
221 {
222 G4cout << "Born ionisation model is initialized " << G4endl
223 << "Energy range: "
224 << LowEnergyLimit() / eV << " eV - "
225 << HighEnergyLimit() / keV << " keV for "
226 << particle->GetParticleName()
227 << G4endl;
228 }
229
230 // Initialize water density pointer
231 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
232
233 //
234 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
235
236 if (isInitialised) { return; }
237 fParticleChangeForGamma = GetParticleChangeForGamma();
238 isInitialised = true;
239}
240
241//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242
243G4double G4DNABornIonisationModel::CrossSectionPerVolume(const G4Material* material,
244 const G4ParticleDefinition* particleDefinition,
245 G4double ekin,
246 G4double,
247 G4double)
248{
249 if (verboseLevel > 3)
250 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel" << G4endl;
251
252 if (
253 particleDefinition != G4Proton::ProtonDefinition()
254 &&
255 particleDefinition != G4Electron::ElectronDefinition()
256 )
257
258 return 0;
259
260 // Calculate total cross section for model
261
262 G4double lowLim = 0;
263 G4double highLim = 0;
264 G4double sigma=0;
265
266 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
267
268 if(waterDensity!= 0.0)
269 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
270 {
271 const G4String& particleName = particleDefinition->GetParticleName();
272
273 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
274 pos1 = lowEnergyLimit.find(particleName);
275 if (pos1 != lowEnergyLimit.end())
276 {
277 lowLim = pos1->second;
278 }
279
280 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
281 pos2 = highEnergyLimit.find(particleName);
282 if (pos2 != highEnergyLimit.end())
283 {
284 highLim = pos2->second;
285 }
286
287 if (ekin >= lowLim && ekin < highLim)
288 {
289 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
290 pos = tableData.find(particleName);
291
292 if (pos != tableData.end())
293 {
294 G4DNACrossSectionDataSet* table = pos->second;
295 if (table != 0)
296 {
297 sigma = table->FindValue(ekin);
298 }
299 }
300 else
301 {
302 G4Exception("G4DNABornIonisationModel::CrossSectionPerVolume","em0002",
303 FatalException,"Model not applicable to particle type.");
304 }
305 }
306
307 if (verboseLevel > 2)
308 {
309 G4cout << "__________________________________" << G4endl;
310 G4cout << "°°° G4DNABornIonisationModel - XS INFO START" << G4endl;
311 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
312 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
313 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
314 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
315 G4cout << "°°° G4DNABornIonisationModel - XS INFO END" << G4endl;
316 }
317
318 } // if (waterMaterial)
319
320 // return sigma*material->GetAtomicNumDensityVector()[1];
321 return sigma*waterDensity;
322}
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326void G4DNABornIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
327 const G4MaterialCutsCouple* ,//must be set!
328 const G4DynamicParticle* particle,
329 G4double,
330 G4double)
331{
332
333 if (verboseLevel > 3)
334 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel" << G4endl;
335
336 G4double lowLim = 0;
337 G4double highLim = 0;
338
339 G4double k = particle->GetKineticEnergy();
340
341 const G4String& particleName = particle->GetDefinition()->GetParticleName();
342
343 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
344 pos1 = lowEnergyLimit.find(particleName);
345
346 if (pos1 != lowEnergyLimit.end())
347 {
348 lowLim = pos1->second;
349 }
350
351 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
352 pos2 = highEnergyLimit.find(particleName);
353
354 if (pos2 != highEnergyLimit.end())
355 {
356 highLim = pos2->second;
357 }
358
359 if (k >= lowLim && k < highLim)
360 {
361 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
362 G4double particleMass = particle->GetDefinition()->GetPDGMass();
363 G4double totalEnergy = k + particleMass;
364 G4double pSquare = k * (totalEnergy + particleMass);
365 G4double totalMomentum = std::sqrt(pSquare);
366
367 G4int ionizationShell = RandomSelect(k,particleName);
368
369 // sample deexcitation
370 // here we assume that H_{2}O electronic levels are the same of Oxigen.
371 // this can be considered true with a rough 10% error in energy on K-shell,
372
373 G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
374 G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
375
376 G4double bindingEnergy = 0;
377 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
378
379 if(fAtomDeexcitation) {
380 G4int Z = 8;
381 G4AtomicShellEnumerator as = fKShell;
382
383 if (ionizationShell <5 && ionizationShell >1)
384 {
385 as = G4AtomicShellEnumerator(4-ionizationShell);
386 }
387 else if (ionizationShell <2)
388 {
389 as = G4AtomicShellEnumerator(3);
390 }
391
392 // FOR DEBUG ONLY
393 // if (ionizationShell == 4) {
394 //
395 // G4cout << "Z: " << Z << " as: " << as
396 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
397 // G4cout << "Press <Enter> key to continue..." << G4endl;
398 // G4cin.ignore();
399 // }
400
401 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
402 secNumberInit = fvect->size();
403 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
404 secNumberFinal = fvect->size();
405 }
406
407 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
408
409 G4double cosTheta = 0.;
410 G4double phi = 0.;
411 RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
412
413 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
414 G4double dirX = sinTheta*std::cos(phi);
415 G4double dirY = sinTheta*std::sin(phi);
416 G4double dirZ = cosTheta;
417 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
418 deltaDirection.rotateUz(primaryDirection);
419
420 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
421 {
422 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
423
424 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
425 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
426 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
427 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
428 finalPx /= finalMomentum;
429 finalPy /= finalMomentum;
430 finalPz /= finalMomentum;
431
432 G4ThreeVector direction;
433 direction.set(finalPx,finalPy,finalPz);
434
435 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
436 }
437
438 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
439
440 // note thta secondaryKinetic is the nergy of the delta ray, not of all secondaries.
441 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
442 G4double deexSecEnergy = 0;
443 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
444
445 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
446
447 }
448
449 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
450 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
451
452 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
453 fvect->push_back(dp);
454
455
456 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
457 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
458 ionizationShell,
459 theIncomingTrack);
460 }
461
462}
463
464//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
465
466G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
467 G4double k, G4int shell)
468{
469 if (particleDefinition == G4Electron::ElectronDefinition())
470 {
471 G4double maximumEnergyTransfer=0.;
472 if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
473 else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
474
475 // SI : original method
476 /*
477 G4double crossSectionMaximum = 0.;
478 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
479 {
480 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
481 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
482 }
483*/
484
485
486 // SI : alternative method
487
488 G4double crossSectionMaximum = 0.;
489
490 G4double minEnergy = waterStructure.IonisationEnergy(shell);
491 G4double maxEnergy = maximumEnergyTransfer;
492 G4int nEnergySteps = 50;
493
494 G4double value(minEnergy);
495 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
496 G4int step(nEnergySteps);
497 while (step>0)
498 {
499 step--;
500 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
501 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
502 value*=stpEnergy;
503 }
504 //
505
506 G4double secondaryElectronKineticEnergy=0.;
507 do
508 {
509 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
510 } while(G4UniformRand()*crossSectionMaximum >
511 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
512
513 return secondaryElectronKineticEnergy;
514
515 }
516
517 else if (particleDefinition == G4Proton::ProtonDefinition())
518 {
519 G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
520
521 G4double crossSectionMaximum = 0.;
522 for (G4double value = waterStructure.IonisationEnergy(shell);
523 value<=4.*waterStructure.IonisationEnergy(shell) ;
524 value+=0.1*eV)
525 {
526 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
527 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
528 }
529
530 G4double secondaryElectronKineticEnergy = 0.;
531 do
532 {
533 secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
534 } while(G4UniformRand()*crossSectionMaximum >=
535 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
536
537 return secondaryElectronKineticEnergy;
538 }
539
540 return 0;
541}
542
543//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
544
545void G4DNABornIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
546 G4double k,
547 G4double secKinetic,
548 G4double & cosTheta,
549 G4double & phi )
550{
551 if (particleDefinition == G4Electron::ElectronDefinition())
552 {
553 phi = twopi * G4UniformRand();
554 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
555 else if (secKinetic <= 200.*eV)
556 {
557 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
558 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
559 }
560 else
561 {
562 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
563 cosTheta = std::sqrt(1.-sin2O);
564 }
565 }
566
567 else if (particleDefinition == G4Proton::ProtonDefinition())
568 {
569 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
570 phi = twopi * G4UniformRand();
571
572 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
573
574 // Restriction below 100 eV from Emfietzoglou (2000)
575
576 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
577 else cosTheta = (2.*G4UniformRand())-1.;
578
579 }
580}
581
582//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
583
584double G4DNABornIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
585 G4double k,
586 G4double energyTransfer,
587 G4int ionizationLevelIndex)
588{
589 G4double sigma = 0.;
590
591 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
592 {
593 G4double valueT1 = 0;
594 G4double valueT2 = 0;
595 G4double valueE21 = 0;
596 G4double valueE22 = 0;
597 G4double valueE12 = 0;
598 G4double valueE11 = 0;
599
600 G4double xs11 = 0;
601 G4double xs12 = 0;
602 G4double xs21 = 0;
603 G4double xs22 = 0;
604
605 if (particleDefinition == G4Electron::ElectronDefinition())
606 {
607 // k should be in eV and energy transfer eV also
608
609 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
610
611 std::vector<double>::iterator t1 = t2-1;
612
613 // SI : the following condition avoids situations where energyTransfer >last vector element
614 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
615 {
616 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
617 std::vector<double>::iterator e11 = e12-1;
618
619 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
620 std::vector<double>::iterator e21 = e22-1;
621
622 valueT1 =*t1;
623 valueT2 =*t2;
624 valueE21 =*e21;
625 valueE22 =*e22;
626 valueE12 =*e12;
627 valueE11 =*e11;
628
629 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
630 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
631 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
632 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
633 }
634
635 }
636
637 if (particleDefinition == G4Proton::ProtonDefinition())
638 {
639 // k should be in eV and energy transfer eV also
640 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
641 std::vector<double>::iterator t1 = t2-1;
642
643 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
644 std::vector<double>::iterator e11 = e12-1;
645
646 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
647 std::vector<double>::iterator e21 = e22-1;
648
649 valueT1 =*t1;
650 valueT2 =*t2;
651 valueE21 =*e21;
652 valueE22 =*e22;
653 valueE12 =*e12;
654 valueE11 =*e11;
655
656 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
657 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
658 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
659 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
660
661 }
662
663 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
664 if (xsProduct != 0.)
665 {
666 sigma = QuadInterpolator( valueE11, valueE12,
667 valueE21, valueE22,
668 xs11, xs12,
669 xs21, xs22,
670 valueT1, valueT2,
671 k, energyTransfer);
672 }
673
674 }
675
676 return sigma;
677}
678
679//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
680
681G4double G4DNABornIonisationModel::LogLogInterpolate(G4double e1,
682 G4double e2,
683 G4double e,
684 G4double xs1,
685 G4double xs2)
686{
687 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
688 G4double b = std::log10(xs2) - a*std::log10(e2);
689 G4double sigma = a*std::log10(e) + b;
690 G4double value = (std::pow(10.,sigma));
691 return value;
692}
693
694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
695
696G4double G4DNABornIonisationModel::QuadInterpolator(G4double e11, G4double e12,
697 G4double e21, G4double e22,
698 G4double xs11, G4double xs12,
699 G4double xs21, G4double xs22,
700 G4double t1, G4double t2,
701 G4double t, G4double e)
702{
703 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
704 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
705 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
706 return value;
707}
708
709//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
710
711G4int G4DNABornIonisationModel::RandomSelect(G4double k, const G4String& particle )
712{
713 G4int level = 0;
714
715 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
716 pos = tableData.find(particle);
717
718 if (pos != tableData.end())
719 {
720 G4DNACrossSectionDataSet* table = pos->second;
721
722 if (table != 0)
723 {
724 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
725 const size_t n(table->NumberOfComponents());
726 size_t i(n);
727 G4double value = 0.;
728
729 while (i>0)
730 {
731 i--;
732 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
733 value += valuesBuffer[i];
734 }
735
736 value *= G4UniformRand();
737
738 i = n;
739
740 while (i > 0)
741 {
742 i--;
743
744 if (valuesBuffer[i] > value)
745 {
746 delete[] valuesBuffer;
747 return i;
748 }
749 value -= valuesBuffer[i];
750 }
751
752 if (valuesBuffer) delete[] valuesBuffer;
753
754 }
755 }
756 else
757 {
758 G4Exception("G4DNABornIonisationModel::RandomSelect","em0002",
759 FatalException,"Model not applicable to particle type.");
760 }
761
762 return level;
763}
764
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:52
eIonizedMolecule
@ eIonizedMolecule
Definition: G4DNAChemistryManager.hh:54
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
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:72
G4DNABornIonisationModel::G4DNABornIonisationModel
G4DNABornIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
Definition: G4DNABornIonisationModel.cc:43
G4DNABornIonisationModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNABornIonisationModel.cc:243
G4DNABornIonisationModel::DifferentialCrossSection
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
Definition: G4DNABornIonisationModel.cc:584
G4DNABornIonisationModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNABornIonisationModel.hh:76
G4DNABornIonisationModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
Definition: G4DNABornIonisationModel.cc:91
G4DNABornIonisationModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNABornIonisationModel.cc:326
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:64
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:144
G4DNACrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4DNACrossSectionDataSet.cc:500
G4DNACrossSectionDataSet::NumberOfComponents
virtual size_t NumberOfComponents(void) const
Definition: G4DNACrossSectionDataSet.hh:75
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:69
G4DNACrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4DNACrossSectionDataSet.cc:84
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:64
G4DNAMolecularMaterial::GetNumMolPerVolTableFor
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
Definition: G4DNAMolecularMaterial.cc:367
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.cc:1106
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:261
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
G4ParticleChangeForGamma::GetCurrentTrack
const G4Track * GetCurrentTrack() const
Definition: G4ParticleChangeForGamma.hh:154
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
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.hh:284
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41