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