Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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