Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
30 #include "G4DNARuddIonisationModel.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4UAtomicDeexcitation.hh"
34 #include "G4LossTableManager.hh"
35 #include "G4DNAChemistryManager.hh"
36 #include "G4DNAMolecularMaterial.hh"
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
40 using namespace std;
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 
44 G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
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
79  SetDeexcitationFlag(true);
80  fAtomDeexcitation = 0;
81  fParticleChangeForGamma = 0;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
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 
106 void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
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;
122  instance = G4DNAGenericIonsManager::Instance();
123  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
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 
149  G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
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 
165  G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
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 
182  G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
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 
199  G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
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 
215  G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
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
264  fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
265 
266  //
267 
268  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
269 
270  if (isInitialised) { return; }
271  fParticleChangeForGamma = GetParticleChangeForGamma();
272  isInitialised = true;
273 
274 }
275 
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277 
278 G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
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;
290  instance = 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 
404 void 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;
417  instance = 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;
480  G4AtomicShellEnumerator as = fKShell;
481 
482  if (ionizationShell <5 && ionizationShell >1)
483  {
484  as = G4AtomicShellEnumerator(4-ionizationShell);
485  }
486  else if (ionizationShell <2)
487  {
488  as = G4AtomicShellEnumerator(3);
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  */
539  fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
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 
549  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
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();
562  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
563  ionizationShell,
564  theIncomingTrack);
565  }
566 
567  // SI - not useful since low energy of model is 0 eV
568 
569  if (k < lowLim)
570  {
571  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
572  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
573  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
574  }
575 
576 }
577 
578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
579 
580 G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
581  G4double k,
582  G4int shell)
583 {
584  G4double maximumKineticEnergyTransfer = 0.;
585 
586  G4DNAGenericIonsManager *instance;
587  instance = 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 
627 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
628  G4double k,
629  G4double secKinetic,
630  G4double & cosTheta,
631  G4double & phi )
632 {
633  G4DNAGenericIonsManager *instance;
634  instance = 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 
665 G4double 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;
738  instance = 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 
866 G4double 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 
882 G4double 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 
899 G4double 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 
915 G4double 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 
933 G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
934 {
935  G4DNAGenericIonsManager *instance;
936  instance = 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 
957 G4int 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 
1022 G4double 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 
1075 G4double G4DNARuddIonisationModel::Sum(G4double /* energy */, const G4String& /* particle */)
1076 {
1077  return 0;
1078 }
1079 
Generated on Sat May 25 2013 14:33:31 for Geant4 by   doxygen 1.8.4