Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNARuddIonisationModel.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// GEANT4 tag $Name: $
28//
29
32#include "G4SystemOfUnits.hh"
34#include "G4LossTableManager.hh"
37
38//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39
40using namespace std;
41
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43
45 const G4String& nam)
46 :G4VEmModel(nam),isInitialised(false)
47{
48 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
49 fpWaterDensity = 0;
50
51 slaterEffectiveCharge[0]=0.;
52 slaterEffectiveCharge[1]=0.;
53 slaterEffectiveCharge[2]=0.;
54 sCoefficient[0]=0.;
55 sCoefficient[1]=0.;
56 sCoefficient[2]=0.;
57
58 lowEnergyLimitForZ1 = 0 * eV;
59 lowEnergyLimitForZ2 = 0 * eV;
60 lowEnergyLimitOfModelForZ1 = 100 * eV;
61 lowEnergyLimitOfModelForZ2 = 1 * keV;
62 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
63 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
64
65 verboseLevel= 0;
66 // Verbosity scale:
67 // 0 = nothing
68 // 1 = warning for energy non-conservation
69 // 2 = details of energy budget
70 // 3 = calculation of cross sections, file openings, sampling of atoms
71 // 4 = entering in methods
72
73 if( verboseLevel>0 )
74 {
75 G4cout << "Rudd ionisation model is constructed " << G4endl;
76 }
77
78 //Mark this model as "applicable" for atomic deexcitation
80 fAtomDeexcitation = 0;
82}
83
84//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
87{
88 // Cross section
89
90 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
91 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
92 {
93 G4DNACrossSectionDataSet* table = pos->second;
94 delete table;
95 }
96
97
98 // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
99 // Coverity however will signal this as an error
100 //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
101
102}
103
104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105
107 const G4DataVector& /*cuts*/)
108{
109
110 if (verboseLevel > 3)
111 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
112
113 // Energy limits
114
115 G4String fileProton("dna/sigma_ionisation_p_rudd");
116 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
117 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
118 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
119 G4String fileHelium("dna/sigma_ionisation_he_rudd");
120
121 G4DNAGenericIonsManager *instance;
124 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
125 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
126 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
127 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
128
129 G4String proton;
130 G4String hydrogen;
131 G4String alphaPlusPlus;
132 G4String alphaPlus;
133 G4String helium;
134
135 G4double scaleFactor = 1 * m*m;
136
137 // LIMITS AND DATA
138
139 // ********************************************************
140
141 proton = protonDef->GetParticleName();
142 tableFile[proton] = fileProton;
143
144 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
145 highEnergyLimit[proton] = 500. * keV;
146
147 // Cross section
148
150 eV,
151 scaleFactor );
152 tableProton->LoadData(fileProton);
153 tableData[proton] = tableProton;
154
155 // ********************************************************
156
157 hydrogen = hydrogenDef->GetParticleName();
158 tableFile[hydrogen] = fileHydrogen;
159
160 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
161 highEnergyLimit[hydrogen] = 100. * MeV;
162
163 // Cross section
164
166 eV,
167 scaleFactor );
168 tableHydrogen->LoadData(fileHydrogen);
169
170 tableData[hydrogen] = tableHydrogen;
171
172 // ********************************************************
173
174 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
175 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
176
177 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
178 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
179
180 // Cross section
181
183 eV,
184 scaleFactor );
185 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
186
187 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
188
189 // ********************************************************
190
191 alphaPlus = alphaPlusDef->GetParticleName();
192 tableFile[alphaPlus] = fileAlphaPlus;
193
194 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
195 highEnergyLimit[alphaPlus] = 400. * MeV;
196
197 // Cross section
198
200 eV,
201 scaleFactor );
202 tableAlphaPlus->LoadData(fileAlphaPlus);
203 tableData[alphaPlus] = tableAlphaPlus;
204
205 // ********************************************************
206
207 helium = heliumDef->GetParticleName();
208 tableFile[helium] = fileHelium;
209
210 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
211 highEnergyLimit[helium] = 400. * MeV;
212
213 // Cross section
214
216 eV,
217 scaleFactor );
218 tableHelium->LoadData(fileHelium);
219 tableData[helium] = tableHelium;
220
221 //
222
223 if (particle==protonDef)
224 {
225 SetLowEnergyLimit(lowEnergyLimit[proton]);
226 SetHighEnergyLimit(highEnergyLimit[proton]);
227 }
228
229 if (particle==hydrogenDef)
230 {
231 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
232 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
233 }
234
235 if (particle==heliumDef)
236 {
237 SetLowEnergyLimit(lowEnergyLimit[helium]);
238 SetHighEnergyLimit(highEnergyLimit[helium]);
239 }
240
241 if (particle==alphaPlusDef)
242 {
243 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
244 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
245 }
246
247 if (particle==alphaPlusPlusDef)
248 {
249 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
250 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
251 }
252
253 if( verboseLevel>0 )
254 {
255 G4cout << "Rudd ionisation model is initialized " << G4endl
256 << "Energy range: "
257 << LowEnergyLimit() / eV << " eV - "
258 << HighEnergyLimit() / keV << " keV for "
259 << particle->GetParticleName()
260 << G4endl;
261 }
262
263 // Initialize water density pointer
265
266 //
267
268 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
269
270 if (isInitialised) { return; }
272 isInitialised = true;
273
274}
275
276//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277
279 const G4ParticleDefinition* particleDefinition,
280 G4double k,
281 G4double,
282 G4double)
283{
284 if (verboseLevel > 3)
285 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel" << G4endl;
286
287 // Calculate total cross section for model
288
289 G4DNAGenericIonsManager *instance;
291
292 if (
293 particleDefinition != G4Proton::ProtonDefinition()
294 &&
295 particleDefinition != instance->GetIon("hydrogen")
296 &&
297 particleDefinition != instance->GetIon("alpha++")
298 &&
299 particleDefinition != instance->GetIon("alpha+")
300 &&
301 particleDefinition != instance->GetIon("helium")
302 )
303
304 return 0;
305
306 G4double lowLim = 0;
307
308 if ( particleDefinition == G4Proton::ProtonDefinition()
309 || particleDefinition == instance->GetIon("hydrogen")
310 )
311
312 lowLim = lowEnergyLimitOfModelForZ1;
313
314 if ( particleDefinition == instance->GetIon("alpha++")
315 || particleDefinition == instance->GetIon("alpha+")
316 || particleDefinition == instance->GetIon("helium")
317 )
318
319 lowLim = lowEnergyLimitOfModelForZ2;
320
321 G4double highLim = 0;
322 G4double sigma=0;
323
324 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
325
326 if(waterDensity!= 0.0)
327 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
328 {
329 const G4String& particleName = particleDefinition->GetParticleName();
330
331 // SI - the following is useless since lowLim is already defined
332/*
333 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
334 pos1 = lowEnergyLimit.find(particleName);
335
336 if (pos1 != lowEnergyLimit.end())
337 {
338 lowLim = pos1->second;
339 }
340*/
341
342 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
343 pos2 = highEnergyLimit.find(particleName);
344
345 if (pos2 != highEnergyLimit.end())
346 {
347 highLim = pos2->second;
348 }
349
350 if (k <= highLim)
351 {
352 //SI : XS must not be zero otherwise sampling of secondaries method ignored
353
354 if (k < lowLim) k = lowLim;
355
356 //
357
358 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
359 pos = tableData.find(particleName);
360
361 if (pos != tableData.end())
362 {
363 G4DNACrossSectionDataSet* table = pos->second;
364 if (table != 0)
365 {
366 sigma = table->FindValue(k);
367 }
368 }
369 else
370 {
371 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
372 FatalException,"Model not applicable to particle type.");
373 }
374
375 } // if (k >= lowLim && k < highLim)
376
377 if (verboseLevel > 2)
378 {
379 G4cout << "__________________________________" << G4endl;
380 G4cout << "°°° G4DNARuddIonisationModel - XS INFO START" << G4endl;
381 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
382 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
383 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
384 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
385 G4cout << "°°° G4DNARuddIonisationModel - XS INFO END" << G4endl;
386 }
387
388 } // if (waterMaterial)
389 else
390 {
391 if (verboseLevel > 2)
392 {
393 G4cout << "Warning : RuddIonisationModel: WATER DENSITY IS NULL" << G4endl;
394 }
395 }
396
397 return sigma*waterDensity;
398 // return sigma*material->GetAtomicNumDensityVector()[1];
399
400}
401
402//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
403
404void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
405 const G4MaterialCutsCouple* /*couple*/,
406 const G4DynamicParticle* particle,
407 G4double,
408 G4double)
409{
410 if (verboseLevel > 3)
411 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel" << G4endl;
412
413 G4double lowLim = 0;
414 G4double highLim = 0;
415
416 G4DNAGenericIonsManager *instance;
418
419 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
420 || particle->GetDefinition() == instance->GetIon("hydrogen")
421 )
422
423 lowLim = killBelowEnergyForZ1;
424
425 if ( particle->GetDefinition() == instance->GetIon("alpha++")
426 || particle->GetDefinition() == instance->GetIon("alpha+")
427 || particle->GetDefinition() == instance->GetIon("helium")
428 )
429
430 lowLim = killBelowEnergyForZ2;
431
432 G4double k = particle->GetKineticEnergy();
433
434 const G4String& particleName = particle->GetDefinition()->GetParticleName();
435
436 // SI - the following is useless since lowLim is already defined
437 /*
438 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
439 pos1 = lowEnergyLimit.find(particleName);
440
441 if (pos1 != lowEnergyLimit.end())
442 {
443 lowLim = pos1->second;
444 }
445 */
446
447 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
448 pos2 = highEnergyLimit.find(particleName);
449
450 if (pos2 != highEnergyLimit.end())
451 {
452 highLim = pos2->second;
453 }
454
455 if (k >= lowLim && k <= highLim)
456 {
457 G4ParticleDefinition* definition = particle->GetDefinition();
458 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
459 /*
460 G4double particleMass = definition->GetPDGMass();
461 G4double totalEnergy = k + particleMass;
462 G4double pSquare = k*(totalEnergy+particleMass);
463 G4double totalMomentum = std::sqrt(pSquare);
464 */
465
466 G4int ionizationShell = RandomSelect(k,particleName);
467
468 // sample deexcitation
469 // here we assume that H_{2}O electronic levels are the same of Oxigen.
470 // this can be considered true with a rough 10% error in energy on K-shell,
471
472 G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
473 G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
474
475 G4double bindingEnergy = 0;
476 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
477
478 if(fAtomDeexcitation) {
479 G4int Z = 8;
481
482 if (ionizationShell <5 && ionizationShell >1)
483 {
484 as = G4AtomicShellEnumerator(4-ionizationShell);
485 }
486 else if (ionizationShell <2)
487 {
489 }
490
491
492
493 // DEBUG
494 // if (ionizationShell == 4) {
495 //
496 // G4cout << "Z: " << Z << " as: " << as
497 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
498 // G4cout << "Press <Enter> key to continue..." << G4endl;
499 // G4cin.ignore();
500 // }
501
502 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
503 secNumberInit = fvect->size();
504 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
505 secNumberFinal = fvect->size();
506 }
507
508
509 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
510
511 G4double cosTheta = 0.;
512 G4double phi = 0.;
513 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
514
515 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
516 G4double dirX = sinTheta*std::cos(phi);
517 G4double dirY = sinTheta*std::sin(phi);
518 G4double dirZ = cosTheta;
519 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
520 deltaDirection.rotateUz(primaryDirection);
521
522 // Ignored for ions on electrons
523 /*
524 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
525
526 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
527 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
528 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
529 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
530 finalPx /= finalMomentum;
531 finalPy /= finalMomentum;
532 finalPz /= finalMomentum;
533
534 G4ThreeVector direction;
535 direction.set(finalPx,finalPy,finalPz);
536
537 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
538 */
540
541 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
542 G4double deexSecEnergy = 0;
543 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
544
545 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
546
547 }
548
550 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
551
552 // debug
553 // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
554 // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
555 // = bindingEnergy-deexSecEnergy
556 // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
557
558 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
559 fvect->push_back(dp);
560
561 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
563 ionizationShell,
564 theIncomingTrack);
565 }
566
567 // SI - not useful since low energy of model is 0 eV
568
569 if (k < lowLim)
570 {
574 }
575
576}
577
578//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
579
580G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
581 G4double k,
582 G4int shell)
583{
584 G4double maximumKineticEnergyTransfer = 0.;
585
586 G4DNAGenericIonsManager *instance;
588
589 if (particleDefinition == G4Proton::ProtonDefinition()
590 || particleDefinition == instance->GetIon("hydrogen"))
591 {
592 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
593 }
594
595 else if (particleDefinition == instance->GetIon("helium")
596 || particleDefinition == instance->GetIon("alpha+")
597 || particleDefinition == instance->GetIon("alpha++"))
598 {
599 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
600 }
601
602 G4double crossSectionMaximum = 0.;
603
604 for(G4double value=waterStructure.IonisationEnergy(shell); value<=5.*waterStructure.IonisationEnergy(shell) && k>=value ; value+=0.1*eV)
605 {
606 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
607 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
608 }
609
610
611 G4double secElecKinetic = 0.;
612
613 do
614 {
615 secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
616 } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
617 k,
618 secElecKinetic+waterStructure.IonisationEnergy(shell),
619 shell));
620
621 return(secElecKinetic);
622}
623
624//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
625
626
627void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
628 G4double k,
629 G4double secKinetic,
630 G4double & cosTheta,
631 G4double & phi )
632{
633 G4DNAGenericIonsManager *instance;
635
636 G4double maxSecKinetic = 0.;
637
638 if (particleDefinition == G4Proton::ProtonDefinition()
639 || particleDefinition == instance->GetIon("hydrogen"))
640 {
641 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
642 }
643
644 else if (particleDefinition == instance->GetIon("helium")
645 || particleDefinition == instance->GetIon("alpha+")
646 || particleDefinition == instance->GetIon("alpha++"))
647 {
648 maxSecKinetic = 4.* (0.511 / 3728) * k;
649 }
650
651 phi = twopi * G4UniformRand();
652
653 //cosTheta = std::sqrt(secKinetic / maxSecKinetic);
654
655 // Restriction below 100 eV from Emfietzoglou (2000)
656
657 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
658 else cosTheta = (2.*G4UniformRand())-1.;
659
660}
661
662//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
663
664
665G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
666 G4double k,
667 G4double energyTransfer,
668 G4int ionizationLevelIndex)
669{
670 // Shells ids are 0 1 2 3 4 (4 is k shell)
671 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
672 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
673 //
674 // ds S F1(nu) + w * F2(nu)
675 // ---- = G(k) * ---- -------------------------------------------
676 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
677 //
678 // w is the secondary electron kinetic Energy in eV
679 //
680 // All the other parameters can be found in Rudd's Papers
681 //
682 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
683 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
684 //
685
686 const G4int j=ionizationLevelIndex;
687
688 G4double A1 ;
689 G4double B1 ;
690 G4double C1 ;
691 G4double D1 ;
692 G4double E1 ;
693 G4double A2 ;
694 G4double B2 ;
695 G4double C2 ;
696 G4double D2 ;
697 G4double alphaConst ;
698
699 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
700 // The following values are provided by M. dingfelder (priv. comm)
701 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
702
703 if (j == 4)
704 {
705 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
706 A1 = 1.25;
707 B1 = 0.5;
708 C1 = 1.00;
709 D1 = 1.00;
710 E1 = 3.00;
711 A2 = 1.10;
712 B2 = 1.30;
713 C2 = 1.00;
714 D2 = 0.00;
715 alphaConst = 0.66;
716 }
717 else
718 {
719 //Data For Liquid Water from Dingfelder (Protons in Water)
720 A1 = 1.02;
721 B1 = 82.0;
722 C1 = 0.45;
723 D1 = -0.80;
724 E1 = 0.38;
725 A2 = 1.07;
726 // Value provided by M. Dingfelder (priv. comm)
727 B2 = 11.6;
728 //
729 C2 = 0.60;
730 D2 = 0.04;
731 alphaConst = 0.64;
732 }
733
734 const G4double n = 2.;
735 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
736
737 G4DNAGenericIonsManager* instance;
739
740 G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
741 if (wBig<0) return 0.;
742
743 G4double w = wBig / Bj[ionizationLevelIndex];
744 // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
745 if (j==4) w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
746
747 G4double Ry = 13.6*eV;
748
749 G4double tau = 0.;
750
751 G4bool isProtonOrHydrogen = false;
752 G4bool isHelium = false;
753
754 if (particleDefinition == G4Proton::ProtonDefinition()
755 || particleDefinition == instance->GetIon("hydrogen"))
756 {
757 isProtonOrHydrogen = true;
758 tau = (electron_mass_c2/proton_mass_c2) * k ;
759 }
760
761 else if ( particleDefinition == instance->GetIon("helium")
762 || particleDefinition == instance->GetIon("alpha+")
763 || particleDefinition == instance->GetIon("alpha++"))
764 {
765 isHelium = true;
766 tau = (0.511/3728.) * k ;
767 }
768
769 G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
770 if (j==4) S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
771
772 G4double v2 = tau / Bj[ionizationLevelIndex];
773 if (j==4) v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
774
775 G4double v = std::sqrt(v2);
776 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
777 if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
778
779 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
780 G4double L2 = C2*std::pow(v,(D2));
781 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
782 G4double H2 = (A2/v2) + (B2/(v2*v2));
783
784 G4double F1 = L1+H1;
785 G4double F2 = (L2*H2)/(L2+H2);
786
787 G4double sigma = CorrectionFactor(particleDefinition, k)
788 * Gj[j] * (S/Bj[ionizationLevelIndex])
789 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
790
791 if (j==4) sigma = CorrectionFactor(particleDefinition, k)
792 * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
793 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
794
795 if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4))
796 {
797 // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
798 sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
799 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
800 }
801// if ( particleDefinition == G4Proton::ProtonDefinition()
802// || particleDefinition == instance->GetIon("hydrogen")
803// )
804 if(isProtonOrHydrogen)
805 {
806 return(sigma);
807 }
808
809 if (particleDefinition == instance->GetIon("alpha++") )
810 {
811 slaterEffectiveCharge[0]=0.;
812 slaterEffectiveCharge[1]=0.;
813 slaterEffectiveCharge[2]=0.;
814 sCoefficient[0]=0.;
815 sCoefficient[1]=0.;
816 sCoefficient[2]=0.;
817 }
818
819 else if (particleDefinition == instance->GetIon("alpha+") )
820 {
821 slaterEffectiveCharge[0]=2.0;
822 // The following values are provided by M. Dingfelder (priv. comm)
823 slaterEffectiveCharge[1]=2.0;
824 slaterEffectiveCharge[2]=2.0;
825 //
826 sCoefficient[0]=0.7;
827 sCoefficient[1]=0.15;
828 sCoefficient[2]=0.15;
829 }
830
831 else if (particleDefinition == instance->GetIon("helium") )
832 {
833 slaterEffectiveCharge[0]=1.7;
834 slaterEffectiveCharge[1]=1.15;
835 slaterEffectiveCharge[2]=1.15;
836 sCoefficient[0]=0.5;
837 sCoefficient[1]=0.25;
838 sCoefficient[2]=0.25;
839 }
840
841// if ( particleDefinition == instance->GetIon("helium")
842// || particleDefinition == instance->GetIon("alpha+")
843// || particleDefinition == instance->GetIon("alpha++")
844// )
845 if(isHelium)
846 {
847 sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
848
849 if (j==4) sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
850 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
851
852 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
853
854 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
855 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
856 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
857
858 return zEff * zEff * sigma ;
859 }
860
861 return 0;
862}
863
864//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
865
866G4double G4DNARuddIonisationModel::S_1s(G4double t,
867 G4double energyTransferred,
868 G4double slaterEffectiveChg,
869 G4double shellNumber)
870{
871 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
872 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
873
874 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
875 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
876
877 return value;
878}
879
880//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
881
882G4double G4DNARuddIonisationModel::S_2s(G4double t,
883 G4double energyTransferred,
884 G4double slaterEffectiveChg,
885 G4double shellNumber)
886{
887 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
888 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
889
890 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
891 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
892
893 return value;
894
895}
896
897//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
898
899G4double G4DNARuddIonisationModel::S_2p(G4double t,
900 G4double energyTransferred,
901 G4double slaterEffectiveChg,
902 G4double shellNumber)
903{
904 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
905 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
906
907 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
908 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
909
910 return value;
911}
912
913//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
914
915G4double G4DNARuddIonisationModel::R(G4double t,
916 G4double energyTransferred,
917 G4double slaterEffectiveChg,
918 G4double shellNumber)
919{
920 // tElectron = m_electron / m_alpha * t
921 // Dingfelder, in Chattanooga 2005 proceedings, p 4
922
923 G4double tElectron = 0.511/3728. * t;
924 // The following values are provided by M. Dingfelder (priv. comm)
925 G4double H = 2.*13.60569172 * eV;
926 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
927
928 return value;
929}
930
931//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
932
933G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
934{
935 G4DNAGenericIonsManager *instance;
937
938 if (particleDefinition == G4Proton::Proton())
939 {
940 return(1.);
941 }
942 else
943 if (particleDefinition == instance->GetIon("hydrogen"))
944 {
945 G4double value = (std::log10(k/eV)-4.2)/0.5;
946 // The following values are provided by M. Dingfelder (priv. comm)
947 return((0.6/(1+std::exp(value))) + 0.9);
948 }
949 else
950 {
951 return(1.);
952 }
953}
954
955//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
956
957G4int G4DNARuddIonisationModel::RandomSelect(G4double k, const G4String& particle )
958{
959
960 // BEGIN PART 1/2 OF ELECTRON CORRECTION
961
962 // add ONE or TWO electron-water ionisation for alpha+ and helium
963
964 G4int level = 0;
965
966 // Retrieve data table corresponding to the current particle type
967
968 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
969 pos = tableData.find(particle);
970
971 if (pos != tableData.end())
972 {
973 G4DNACrossSectionDataSet* table = pos->second;
974
975 if (table != 0)
976 {
977 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
978
979 const size_t n(table->NumberOfComponents());
980 size_t i(n);
981 G4double value = 0.;
982
983 while (i>0)
984 {
985 i--;
986 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
987 value += valuesBuffer[i];
988 }
989
990 value *= G4UniformRand();
991
992 i = n;
993
994 while (i > 0)
995 {
996 i--;
997
998
999 if (valuesBuffer[i] > value)
1000 {
1001 delete[] valuesBuffer;
1002 return i;
1003 }
1004 value -= valuesBuffer[i];
1005 }
1006
1007 if (valuesBuffer) delete[] valuesBuffer;
1008
1009 }
1010 }
1011 else
1012 {
1013 G4Exception("G4DNARuddIonisationModel::RandomSelect","em0002",
1014 FatalException,"Model not applicable to particle type.");
1015 }
1016
1017 return level;
1018}
1019
1020//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1021
1022G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track )
1023{
1024 G4double sigma = 0.;
1025
1026 const G4DynamicParticle* particle = track.GetDynamicParticle();
1027 G4double k = particle->GetKineticEnergy();
1028
1029 G4double lowLim = 0;
1030 G4double highLim = 0;
1031
1032 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1033
1034 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1035 pos1 = lowEnergyLimit.find(particleName);
1036
1037 if (pos1 != lowEnergyLimit.end())
1038 {
1039 lowLim = pos1->second;
1040 }
1041
1042 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1043 pos2 = highEnergyLimit.find(particleName);
1044
1045 if (pos2 != highEnergyLimit.end())
1046 {
1047 highLim = pos2->second;
1048 }
1049
1050 if (k >= lowLim && k <= highLim)
1051 {
1052 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1053 pos = tableData.find(particleName);
1054
1055 if (pos != tableData.end())
1056 {
1057 G4DNACrossSectionDataSet* table = pos->second;
1058 if (table != 0)
1059 {
1060 sigma = table->FindValue(k);
1061 }
1062 }
1063 else
1064 {
1065 G4Exception("G4DNARuddIonisationModel::PartialCrossSection","em0002",
1066 FatalException,"Model not applicable to particle type.");
1067 }
1068 }
1069
1070 return sigma;
1071}
1072
1073//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1074
1075G4double G4DNARuddIonisationModel::Sum(G4double /* energy */, const G4String& /* particle */)
1076{
1077 return 0;
1078}
1079
G4AtomicShellEnumerator
@ eIonizedMolecule
@ FatalException
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define C1
#define G4UniformRand()
Definition: Randomize.hh:53
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4bool LoadData(const G4String &argFileName)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNARuddIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
size_t GetIndex() const
Definition: G4Material.hh:261
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
static G4Proton * Proton()
Definition: G4Proton.cc:93
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4double pi