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

Last change on this file since 1005 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 31.9 KB
RevLine 
[968]1// ********************************************************************
2// * License and Disclaimer *
3// * *
4// * The Geant4 software is copyright of the Copyright Holders of *
5// * the Geant4 Collaboration. It is provided under the terms and *
6// * conditions of the Geant4 Software License, included in the file *
7// * LICENSE and available at http://cern.ch/geant4/license . These *
8// * include a list of copyright holders. *
9// * *
10// * Neither the authors of this software system, nor their employing *
11// * institutes,nor the agencies providing financial support for this *
12// * work make any representation or warranty, express or implied, *
13// * regarding this software system or assume any liability for its *
14// * use. Please see the license in the file LICENSE and URL above *
15// * for the full disclaimer and the limitation of liability. *
16// * *
17// * This code implementation is the result of the scientific and *
18// * technical work of the GEANT4 collaboration. *
19// * By using, copying, modifying or distributing the software (or *
20// * any work based on the software) you agree to acknowledge its *
21// * use in resulting scientific publications, and indicate your *
22// * acceptance of all terms of the Geant4 Software license. *
23// ********************************************************************
24//
25// $Id: G4DNARuddIonisationModel.cc,v 1.4 2009/02/16 11:00:11 sincerti Exp $
26// GEANT4 tag $Name: geant4-09-02-ref-02 $
27//
28
29#include "G4DNARuddIonisationModel.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
32
33using namespace std;
34
35//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36
37G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
38 const G4String& nam)
39:G4VEmModel(nam),isInitialised(false)
40{
41 lowEnergyLimitForZ1 = 0 * eV;
42 lowEnergyLimitForZ2 = 0 * eV;
43 lowEnergyLimitOfModelForZ1 = 100 * eV;
44 lowEnergyLimitOfModelForZ2 = 1 * keV;
45 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
46 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
47
48 verboseLevel= 0;
49 // Verbosity scale:
50 // 0 = nothing
51 // 1 = warning for energy non-conservation
52 // 2 = details of energy budget
53 // 3 = calculation of cross sections, file openings, sampling of atoms
54 // 4 = entering in methods
55
56 G4cout << "Rudd ionisation model is constructed " << G4endl;
57}
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
62{
63 // Cross section
64
65 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
66 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
67 {
68 G4DNACrossSectionDataSet* table = pos->second;
69 delete table;
70 }
71
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
77 const G4DataVector& /*cuts*/)
78{
79
80 if (verboseLevel > 3)
81 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
82
83 // Energy limits
84
85 G4String fileProton("dna/sigma_ionisation_p_rudd");
86 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
87 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
88 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
89 G4String fileHelium("dna/sigma_ionisation_he_rudd");
90
91 G4DNAGenericIonsManager *instance;
92 instance = G4DNAGenericIonsManager::Instance();
93 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
94 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
95 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
96 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
97 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
98
99 G4String proton;
100 G4String hydrogen;
101 G4String alphaPlusPlus;
102 G4String alphaPlus;
103 G4String helium;
104
105 G4double scaleFactor = 1 * m*m;
106
107 if (protonDef != 0)
108 {
109 proton = protonDef->GetParticleName();
110 tableFile[proton] = fileProton;
111
112 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
113 highEnergyLimit[proton] = 500. * keV;
114
115 // Cross section
116
117 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
118 eV,
119 scaleFactor );
120 tableProton->LoadData(fileProton);
121 tableData[proton] = tableProton;
122 }
123 else
124 {
125 G4Exception("G4DNARuddIonisationModel::Initialise: proton is not defined");
126 }
127
128 if (hydrogenDef != 0)
129 {
130 hydrogen = hydrogenDef->GetParticleName();
131 tableFile[hydrogen] = fileHydrogen;
132
133 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
134 highEnergyLimit[hydrogen] = 100. * MeV;
135
136 // Cross section
137
138 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
139 eV,
140 scaleFactor );
141 tableHydrogen->LoadData(fileHydrogen);
142
143 tableData[hydrogen] = tableHydrogen;
144 }
145 else
146 {
147 G4Exception("G4DNARuddIonisationModel::Initialise: hydrogen is not defined");
148 }
149
150 if (alphaPlusPlusDef != 0)
151 {
152 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
153 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
154
155 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
156 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
157
158 // Cross section
159
160 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
161 eV,
162 scaleFactor );
163 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
164
165 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
166 }
167 else
168 {
169 G4Exception("G4DNARuddIonisationModel::Initialise: alphaPlusPlus is not defined");
170 }
171
172 if (alphaPlusDef != 0)
173 {
174 alphaPlus = alphaPlusDef->GetParticleName();
175 tableFile[alphaPlus] = fileAlphaPlus;
176
177 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
178 highEnergyLimit[alphaPlus] = 10. * MeV;
179
180 // Cross section
181
182 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
183 eV,
184 scaleFactor );
185 tableAlphaPlus->LoadData(fileAlphaPlus);
186 tableData[alphaPlus] = tableAlphaPlus;
187 }
188 else
189 {
190 G4Exception("G4DNARuddIonisationModel::Initialise: alphaPlus is not defined");
191 }
192
193 if (heliumDef != 0)
194 {
195 helium = heliumDef->GetParticleName();
196 tableFile[helium] = fileHelium;
197
198 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
199 highEnergyLimit[helium] = 10. * MeV;
200
201 // Cross section
202
203 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
204 eV,
205 scaleFactor );
206 tableHelium->LoadData(fileHelium);
207 tableData[helium] = tableHelium;
208 }
209 else
210 {
211 G4Exception("G4DNARuddIonisationModel::Initialise: helium is not defined");
212 }
213
214 if (particle==protonDef)
215 {
216 SetLowEnergyLimit(lowEnergyLimit[proton]);
217 SetHighEnergyLimit(highEnergyLimit[proton]);
218 }
219
220 if (particle==hydrogenDef)
221 {
222 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
223 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
224 }
225
226 if (particle==heliumDef)
227 {
228 SetLowEnergyLimit(lowEnergyLimit[helium]);
229 SetHighEnergyLimit(highEnergyLimit[helium]);
230 }
231
232 if (particle==alphaPlusDef)
233 {
234 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
235 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
236 }
237
238 if (particle==alphaPlusPlusDef)
239 {
240 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
241 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
242 }
243
244 G4cout << "Rudd ionisation model is initialized " << G4endl
245 << "Energy range: "
246 << LowEnergyLimit() / eV << " eV - "
247 << HighEnergyLimit() / keV << " keV for "
248 << particle->GetParticleName()
249 << G4endl;
250
251 //
252
253 if(!isInitialised)
254 {
255 isInitialised = true;
256
257 if(pParticleChange)
258 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
259 else
260 fParticleChangeForGamma = new G4ParticleChangeForGamma();
261 }
262
263 // InitialiseElementSelectors(particle,cuts);
264
265 // Test if water material
266
267 flagMaterialIsWater= false;
268 densityWater = 0;
269
270 const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
271
272 if(theCoupleTable)
273 {
274 G4int numOfCouples = theCoupleTable->GetTableSize();
275
276 if(numOfCouples>0)
277 {
278 for (G4int i=0; i<numOfCouples; i++)
279 {
280 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
281 const G4Material* material = couple->GetMaterial();
282
283 size_t j = material->GetNumberOfElements();
284 while (j>0)
285 {
286 j--;
287 const G4Element* element(material->GetElement(j));
288 if (element->GetZ() == 8.)
289 {
290 G4double density = material->GetAtomicNumDensityVector()[j];
291 if (density > 0.)
292 {
293 flagMaterialIsWater = true;
294 densityWater = density;
295
296 if (verboseLevel > 3)
297 G4cout << "Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
298 }
299 }
300 }
301
302 }
303 } // if(numOfCouples>0)
304
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
564 }
565
566 // SI - not useful since low energy of model is 0 eV
567
568 if (k < lowLim)
569 {
570 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
571 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
572 }
573
574}
575
576//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
577
578G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
579 G4double k,
580 G4int shell)
581{
582 G4double maximumKineticEnergyTransfer = 0.;
583
584 G4DNAGenericIonsManager *instance;
585 instance = G4DNAGenericIonsManager::Instance();
586
587 if (particleDefinition == G4Proton::ProtonDefinition()
588 || particleDefinition == instance->GetIon("hydrogen"))
589 {
590 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
591 }
592
593 if (particleDefinition == instance->GetIon("helium")
594 || particleDefinition == instance->GetIon("alpha+")
595 || particleDefinition == instance->GetIon("alpha++"))
596 {
597 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
598 }
599
600 G4double crossSectionMaximum = 0.;
601
602 for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
603 {
604 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
605 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
606 }
607
608 G4double secElecKinetic = 0.;
609
610 do
611 {
612 secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
613 } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
614 k,
615 secElecKinetic+waterStructure.IonisationEnergy(shell),
616 shell));
617
618 return(secElecKinetic);
619}
620
621//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
622
623
624void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
625 G4double k,
626 G4double secKinetic,
627 G4double & cosTheta,
628 G4double & phi )
629{
630 G4DNAGenericIonsManager *instance;
631 instance = G4DNAGenericIonsManager::Instance();
632
633 G4double maxSecKinetic = 0.;
634
635 if (particleDefinition == G4Proton::ProtonDefinition()
636 || particleDefinition == instance->GetIon("hydrogen"))
637 {
638 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
639 }
640
641 if (particleDefinition == instance->GetIon("helium")
642 || particleDefinition == instance->GetIon("alpha+")
643 || particleDefinition == instance->GetIon("alpha++"))
644 {
645 maxSecKinetic = 4.* (0.511 / 3728) * k;
646 }
647
648 phi = twopi * G4UniformRand();
649 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
650}
651
652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
653
654
655G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
656 G4double k,
657 G4double energyTransfer,
658 G4int ionizationLevelIndex)
659{
660 // Shells ids are 0 1 2 3 4 (4 is k shell)
661 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
662 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
663 //
664 // ds S F1(nu) + w * F2(nu)
665 // ---- = G(k) * ---- -------------------------------------------
666 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
667 //
668 // w is the secondary electron kinetic Energy in eV
669 //
670 // All the other parameters can be found in Rudd's Papers
671 //
672 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
673 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
674 //
675
676 const G4int j=ionizationLevelIndex;
677
678 G4double A1 ;
679 G4double B1 ;
680 G4double C1 ;
681 G4double D1 ;
682 G4double E1 ;
683 G4double A2 ;
684 G4double B2 ;
685 G4double C2 ;
686 G4double D2 ;
687 G4double alphaConst ;
688
689 if (j == 4)
690 {
691 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
692 A1 = 1.25;
693 B1 = 0.5;
694 C1 = 1.00;
695 D1 = 1.00;
696 E1 = 3.00;
697 A2 = 1.10;
698 B2 = 1.30;
699 C2 = 1.00;
700 D2 = 0.00;
701 alphaConst = 0.66;
702 }
703 else
704 {
705 //Data For Liquid Water from Dingfelder (Protons in Water)
706 A1 = 1.02;
707 B1 = 82.0;
708 C1 = 0.45;
709 D1 = -0.80;
710 E1 = 0.38;
711 A2 = 1.07;
712 B2 = 14.6;
713 C2 = 0.60;
714 D2 = 0.04;
715 alphaConst = 0.64;
716 }
717
718 const G4double n = 2.;
719 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
720
721 G4DNAGenericIonsManager* instance;
722 instance = G4DNAGenericIonsManager::Instance();
723
724 G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
725 G4double w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
726 G4double Ry = 13.6*eV;
727
728 G4double tau = 0.;
729
730 if (particleDefinition == G4Proton::ProtonDefinition()
731 || particleDefinition == instance->GetIon("hydrogen"))
732 {
733 tau = (electron_mass_c2/proton_mass_c2) * k ;
734 }
735
736 if ( particleDefinition == instance->GetIon("helium")
737 || particleDefinition == instance->GetIon("alpha+")
738 || particleDefinition == instance->GetIon("alpha++"))
739 {
740 tau = (0.511/3728.) * k ;
741 }
742
743 G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
744 G4double v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
745 G4double v = std::sqrt(v2);
746 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
747
748 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
749 G4double L2 = C2*std::pow(v,(D2));
750 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
751 G4double H2 = (A2/v2) + (B2/(v2*v2));
752
753 G4double F1 = L1+H1;
754 G4double F2 = (L2*H2)/(L2+H2);
755
756 G4double sigma = CorrectionFactor(particleDefinition, k/eV)
757 * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
758 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
759
760 if ( particleDefinition == G4Proton::ProtonDefinition()
761 || particleDefinition == instance->GetIon("hydrogen")
762 )
763 {
764 return(sigma);
765 }
766
767 if (particleDefinition == instance->GetIon("alpha++") )
768 {
769 slaterEffectiveCharge[0]=0.;
770 slaterEffectiveCharge[1]=0.;
771 slaterEffectiveCharge[2]=0.;
772 sCoefficient[0]=0.;
773 sCoefficient[1]=0.;
774 sCoefficient[2]=0.;
775 }
776
777 if (particleDefinition == instance->GetIon("alpha+") )
778 {
779 slaterEffectiveCharge[0]=2.0;
780 slaterEffectiveCharge[1]=1.15;
781 slaterEffectiveCharge[2]=1.15;
782 sCoefficient[0]=0.7;
783 sCoefficient[1]=0.15;
784 sCoefficient[2]=0.15;
785 }
786
787 if (particleDefinition == instance->GetIon("helium") )
788 {
789 slaterEffectiveCharge[0]=1.7;
790 slaterEffectiveCharge[1]=1.15;
791 slaterEffectiveCharge[2]=1.15;
792 sCoefficient[0]=0.5;
793 sCoefficient[1]=0.25;
794 sCoefficient[2]=0.25;
795 }
796
797 if ( particleDefinition == instance->GetIon("helium")
798 || particleDefinition == instance->GetIon("alpha+")
799 || particleDefinition == instance->GetIon("alpha++")
800 )
801 {
802 sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
803
804 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
805
806 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
807 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
808 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
809
810 return zEff * zEff * sigma ;
811 }
812
813 return 0;
814}
815
816//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
817
818G4double G4DNARuddIonisationModel::S_1s(G4double t,
819 G4double energyTransferred,
820 G4double slaterEffectiveChg,
821 G4double shellNumber)
822{
823 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
824 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
825
826 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
827 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
828
829 return value;
830}
831
832//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
833
834G4double G4DNARuddIonisationModel::S_2s(G4double t,
835 G4double energyTransferred,
836 G4double slaterEffectiveChg,
837 G4double shellNumber)
838{
839 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
840 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
841
842 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
843 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
844
845 return value;
846
847}
848
849//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
850
851G4double G4DNARuddIonisationModel::S_2p(G4double t,
852 G4double energyTransferred,
853 G4double slaterEffectiveChg,
854 G4double shellNumber)
855{
856 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
857 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
858
859 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
860 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
861
862 return value;
863}
864
865//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
866
867G4double G4DNARuddIonisationModel::R(G4double t,
868 G4double energyTransferred,
869 G4double slaterEffectiveChg,
870 G4double shellNumber)
871{
872 // tElectron = m_electron / m_alpha * t
873 // Dingfelder, in Chattanooga 2005 proceedings, p 4
874
875 G4double tElectron = 0.511/3728. * t;
876 G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
877
878 return value;
879}
880
881//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
882
883G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
884{
885 G4DNAGenericIonsManager *instance;
886 instance = G4DNAGenericIonsManager::Instance();
887
888 if (particleDefinition == G4Proton::Proton())
889 {
890 return(1.);
891 }
892 else
893 if (particleDefinition == instance->GetIon("hydrogen"))
894 {
895 G4double value = (std::log(k/eV)-4.2)/0.5;
896 return((0.8/(1+std::exp(value))) + 0.9);
897 }
898 else
899 {
900 return(1.);
901 }
902}
903
904//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
905
906G4int G4DNARuddIonisationModel::RandomSelect(G4double k, const G4String& particle )
907{
908
909 // BEGIN PART 1/2 OF ELECTRON CORRECTION
910
911 // add ONE or TWO electron-water excitation for alpha+ and helium
912
913 G4DNAGenericIonsManager *instance;
914 instance = G4DNAGenericIonsManager::Instance();
915 G4double kElectron(0);
916 G4double electronComponent(0);
917 G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
918
919 if ( particle == instance->GetIon("alpha+")->GetParticleName()
920 ||
921 particle == instance->GetIon("helium")->GetParticleName()
922 )
923 {
924 electronDataset->LoadData("dna/sigma_ionisation_e_born");
925
926 kElectron = k * 0.511/3728;
927
928 electronComponent = electronDataset->FindValue(kElectron);
929
930 }
931
932 delete electronDataset;
933
934 // END PART 1/2 OF ELECTRON CORRECTION
935
936 G4int level = 0;
937
938 // Retrieve data table corresponding to the current particle type
939
940 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
941 pos = tableData.find(particle);
942
943 if (pos != tableData.end())
944 {
945 G4DNACrossSectionDataSet* table = pos->second;
946
947 if (table != 0)
948 {
949 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
950
951 const size_t n(table->NumberOfComponents());
952 size_t i(n);
953 G4double value = 0.;
954
955 while (i>0)
956 {
957 i--;
958 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
959
960 // BEGIN PART 2/2 OF ELECTRON CORRECTION
961
962 if (particle == instance->GetIon("alpha+")->GetParticleName())
963 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronComponent; }
964
965 if (particle == instance->GetIon("helium")->GetParticleName())
966 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronComponent; }
967
968 // BEGIN PART 2/2 OF ELECTRON CORRECTION
969
970 value += valuesBuffer[i];
971 }
972
973 value *= G4UniformRand();
974
975 i = n;
976
977 while (i > 0)
978 {
979 i--;
980
981 if (valuesBuffer[i] > value)
982 {
983 delete[] valuesBuffer;
984 return i;
985 }
986 value -= valuesBuffer[i];
987 }
988
989 if (valuesBuffer) delete[] valuesBuffer;
990
991 }
992 }
993 else
994 {
995 G4Exception("G4DNARuddIonisationModel::RandomSelect: attempting to calculate cross section for wrong particle");
996 }
997
998 return level;
999}
1000
1001//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1002
1003G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track )
1004{
1005 G4double sigma = 0.;
1006
1007 const G4DynamicParticle* particle = track.GetDynamicParticle();
1008 G4double k = particle->GetKineticEnergy();
1009
1010 G4double lowLim = 0;
1011 G4double highLim = 0;
1012
1013 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1014
1015 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1016 pos1 = lowEnergyLimit.find(particleName);
1017
1018 if (pos1 != lowEnergyLimit.end())
1019 {
1020 lowLim = pos1->second;
1021 }
1022
1023 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1024 pos2 = highEnergyLimit.find(particleName);
1025
1026 if (pos2 != highEnergyLimit.end())
1027 {
1028 highLim = pos2->second;
1029 }
1030
1031 if (k >= lowLim && k <= highLim)
1032 {
1033 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1034 pos = tableData.find(particleName);
1035
1036 if (pos != tableData.end())
1037 {
1038 G4DNACrossSectionDataSet* table = pos->second;
1039 if (table != 0)
1040 {
1041 sigma = table->FindValue(k);
1042 }
1043 }
1044 else
1045 {
1046 G4Exception("G4DNARuddIonisationModel::PartialCrossSection: attempting to calculate cross section for wrong particle");
1047 }
1048 }
1049
1050 return sigma;
1051}
1052
1053//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1054
1055G4double G4DNARuddIonisationModel::Sum(G4double /* energy */, const G4String& /* particle */)
1056{
1057 return 0;
1058}
1059
Note: See TracBrowser for help on using the repository browser.