source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc@ 1229

Last change on this file since 1229 was 1228, checked in by garnier, 16 years ago

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