Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNARuddIonisationExtendedModel.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: G4DNARuddIonisationExtendedModel.cc 70171 2013-05-24 13:34:18Z gcosmo $
27 // GEANT4 tag $Name: $
28 //
29 
30 
31 // Modified by Z. Francis to handle HZE && inverse rudd function sampling 26-10-2010
32 
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4UAtomicDeexcitation.hh"
37 #include "G4LossTableManager.hh"
38 #include "G4DNAChemistryManager.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  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
52  fpWaterDensity = 0;
53 
54  slaterEffectiveCharge[0]=0.;
55  slaterEffectiveCharge[1]=0.;
56  slaterEffectiveCharge[2]=0.;
57  sCoefficient[0]=0.;
58  sCoefficient[1]=0.;
59  sCoefficient[2]=0.;
60 
61  lowEnergyLimitForA[1] = 0 * eV;
62  lowEnergyLimitForA[2] = 0 * eV;
63  lowEnergyLimitForA[3] = 0 * eV;
64  lowEnergyLimitOfModelForA[1] = 100 * eV;
65  lowEnergyLimitOfModelForA[4] = 1 * keV;
66  lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
67  killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
68  killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
69  killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
70 
71 
72  verboseLevel= 0;
73  // Verbosity scale:
74  // 0 = nothing
75  // 1 = warning for energy non-conservation
76  // 2 = details of energy budget
77  // 3 = calculation of cross sections, file openings, sampling of atoms
78  // 4 = entering in methods
79 
80  if( verboseLevel>0 )
81  {
82  G4cout << "Rudd ionisation model is constructed " << G4endl;
83  }
84 
85  //Mark this model as "applicable" for atomic deexcitation
86  SetDeexcitationFlag(true);
87  fAtomDeexcitation = 0;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94 {
95  // Cross section
96 
97  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
98  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99  {
100  G4DNACrossSectionDataSet* table = pos->second;
101  delete table;
102  }
103 
104  // The following removal is forbidden G4VEnergyLossModel takes care of deletion
105  // however coverity will signal this as an error
106  //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107 
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
113  const G4DataVector& /*cuts*/)
114 {
115 
116  if (verboseLevel > 3)
117  G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
118 
119  // Energy limits
120 
121  G4String fileProton("dna/sigma_ionisation_p_rudd");
122  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
123  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
124  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
125  G4String fileHelium("dna/sigma_ionisation_he_rudd");
126  G4String fileCarbon("dna/sigma_ionisation_c_rudd");
127  G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
128  G4String fileOxygen("dna/sigma_ionisation_o_rudd");
129  G4String fileIron("dna/sigma_ionisation_fe_rudd");
130 
131  G4DNAGenericIonsManager *instance;
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  G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
139  G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
140  G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
141  G4ParticleDefinition* ironDef = instance->GetIon("iron");
142 
144  G4String hydrogen;
145  G4String alphaPlusPlus;
146  G4String alphaPlus;
147  G4String helium;
148  G4String carbon;
149  G4String nitrogen;
150  G4String oxygen;
151  G4String iron;
152 
153  G4double scaleFactor = 1 * m*m;
154 
155  // LIMITS AND DATA
156 
157  proton = protonDef->GetParticleName();
158  tableFile[proton] = fileProton;
159  lowEnergyLimit[proton] = lowEnergyLimitForA[1];
160  highEnergyLimit[proton] = 500. * keV;
161 
162  // Cross section
163 
165  eV,
166  scaleFactor );
167  tableProton->LoadData(fileProton);
168  tableData[proton] = tableProton;
169 
170  // **********************************************************************************************
171 
172  hydrogen = hydrogenDef->GetParticleName();
173  tableFile[hydrogen] = fileHydrogen;
174 
175  lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
176  highEnergyLimit[hydrogen] = 100. * MeV;
177 
178  // Cross section
179 
181  eV,
182  scaleFactor );
183  tableHydrogen->LoadData(fileHydrogen);
184 
185  tableData[hydrogen] = tableHydrogen;
186 
187  // **********************************************************************************************
188 
189  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
190  tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
191 
192  lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
193  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
194 
195  // Cross section
196 
198  eV,
199  scaleFactor );
200  tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
201 
202  tableData[alphaPlusPlus] = tableAlphaPlusPlus;
203 
204  // **********************************************************************************************
205 
206  alphaPlus = alphaPlusDef->GetParticleName();
207  tableFile[alphaPlus] = fileAlphaPlus;
208 
209  lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
210  highEnergyLimit[alphaPlus] = 400. * MeV;
211 
212  // Cross section
213 
215  eV,
216  scaleFactor );
217  tableAlphaPlus->LoadData(fileAlphaPlus);
218  tableData[alphaPlus] = tableAlphaPlus;
219 
220  // **********************************************************************************************
221 
222  helium = heliumDef->GetParticleName();
223  tableFile[helium] = fileHelium;
224 
225  lowEnergyLimit[helium] = lowEnergyLimitForA[4];
226  highEnergyLimit[helium] = 400. * MeV;
227 
228  // Cross section
229 
231  eV,
232  scaleFactor );
233  tableHelium->LoadData(fileHelium);
234  tableData[helium] = tableHelium;
235 
236  // **********************************************************************************************
237 
238  carbon = carbonDef->GetParticleName();
239  tableFile[carbon] = fileCarbon;
240 
241  lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
242  highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
243 
244  // Cross section
245 
247  eV,
248  scaleFactor );
249  tableCarbon->LoadData(fileCarbon);
250  tableData[carbon] = tableCarbon;
251 
252  // **********************************************************************************************
253 
254  oxygen = oxygenDef->GetParticleName();
255  tableFile[oxygen] = fileOxygen;
256 
257  lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
258  highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
259 
260  // Cross section
261 
263  eV,
264  scaleFactor );
265  tableOxygen->LoadData(fileOxygen);
266  tableData[oxygen] = tableOxygen;
267 
268  // **********************************************************************************************
269 
270  nitrogen = nitrogenDef->GetParticleName();
271  tableFile[nitrogen] = fileNitrogen;
272 
273  lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
274  highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
275 
276  // Cross section
277 
279  eV,
280  scaleFactor );
281  tableNitrogen->LoadData(fileNitrogen);
282  tableData[nitrogen] = tableNitrogen;
283 
284  // **********************************************************************************************
285 
286  iron = ironDef->GetParticleName();
287  tableFile[iron] = fileIron;
288 
289  lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
290  highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
291 
292  // Cross section
293 
295  eV,
296  scaleFactor );
297  tableIron->LoadData(fileIron);
298  tableData[iron] = tableIron;
299 
300  // ZF Following lines can be replaced by:
301  SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
302  SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
303  // at least for HZE
304 
305  /*
306  if (particle==protonDef)
307  {
308  SetLowEnergyLimit(lowEnergyLimit[proton]);
309  SetHighEnergyLimit(highEnergyLimit[proton]);
310  }
311 
312  if (particle==hydrogenDef)
313  {
314  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
315  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
316  }
317 
318  if (particle==heliumDef)
319  {
320  SetLowEnergyLimit(lowEnergyLimit[helium]);
321  SetHighEnergyLimit(highEnergyLimit[helium]);
322  }
323 
324  if (particle==alphaPlusDef)
325  {
326  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
327  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
328  }
329 
330  if (particle==alphaPlusPlusDef)
331  {
332  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
333  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
334  }
335 
336  if (particle==carbonDef)
337  {
338  SetLowEnergyLimit(lowEnergyLimit[carbon]);
339  SetHighEnergyLimit(highEnergyLimit[carbon]);
340  }
341 
342  if (particle==oxygenDef)
343  {
344  SetLowEnergyLimit(lowEnergyLimit[oxygen]);
345  SetHighEnergyLimit(highEnergyLimit[oxygen]);
346  }*/
347 
348  //----------------------------------------------------------------------
349 
350  if( verboseLevel>0 )
351  {
352  G4cout << "Rudd ionisation model is initialized " << G4endl
353  << "Energy range: "
354  << LowEnergyLimit() / eV << " eV - "
355  << HighEnergyLimit() / keV << " keV for "
356  << particle->GetParticleName()
357  << G4endl;
358  }
359 
360  // Initialize water density pointer
362 
363  //
364 
365  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
366 
367  if (isInitialised) { return; }
369  isInitialised = true;
370 }
371 
372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373 
375  const G4ParticleDefinition* particleDefinition,
376  G4double k,
377  G4double,
378  G4double)
379 {
380  if (verboseLevel > 3)
381  G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
382 
383  // Calculate total cross section for model
384 
385  G4DNAGenericIonsManager *instance;
387 
388  if (
389  particleDefinition != G4Proton::ProtonDefinition()
390  &&
391  particleDefinition != instance->GetIon("hydrogen")
392  &&
393  particleDefinition != instance->GetIon("alpha++")
394  &&
395  particleDefinition != instance->GetIon("alpha+")
396  &&
397  particleDefinition != instance->GetIon("helium")
398  &&
399  particleDefinition != instance->GetIon("carbon")
400  &&
401  particleDefinition != instance->GetIon("nitrogen")
402  &&
403  particleDefinition != instance->GetIon("oxygen")
404  &&
405  particleDefinition != instance->GetIon("iron")
406  )
407 
408  return 0;
409 
410  G4double lowLim = 0;
411 
412  if ( particleDefinition == G4Proton::ProtonDefinition()
413  || particleDefinition == instance->GetIon("hydrogen")
414  )
415 
416  lowLim = lowEnergyLimitOfModelForA[1];
417 
418  else if ( particleDefinition == instance->GetIon("alpha++")
419  || particleDefinition == instance->GetIon("alpha+")
420  || particleDefinition == instance->GetIon("helium")
421  )
422 
423  lowLim = lowEnergyLimitOfModelForA[4];
424 
425  else lowLim = lowEnergyLimitOfModelForA[5];
426 
427  G4double highLim = 0;
428  G4double sigma=0;
429 
430 
431  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
432 
433  if(waterDensity!= 0.0)
434 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
435  {
436  const G4String& particleName = particleDefinition->GetParticleName();
437 
438  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439  pos2 = highEnergyLimit.find(particleName);
440 
441  if (pos2 != highEnergyLimit.end())
442  {
443  highLim = pos2->second;
444  }
445 
446  if (k <= highLim)
447  {
448 
449  //SI : XS must not be zero otherwise sampling of secondaries method ignored
450 
451  if (k < lowLim) k = lowLim;
452 
453  //
454 
455  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
456  pos = tableData.find(particleName);
457 
458  if (pos != tableData.end())
459  {
460  G4DNACrossSectionDataSet* table = pos->second;
461  if (table != 0)
462  {
463  sigma = table->FindValue(k);
464  }
465  }
466  else
467  {
468  G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
469  FatalException,"Model not applicable to particle type.");
470  }
471 
472  } // if (k >= lowLim && k < highLim)
473 
474  if (verboseLevel > 2)
475  {
476  G4cout << "__________________________________" << G4endl;
477  G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
478  G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
479  G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
480  G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
481  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
482  G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
483 
484  }
485 
486  } // if (waterMaterial)
487 
488  return sigma*waterDensity;
489 // return sigma*material->GetAtomicNumDensityVector()[1];
490 
491 }
492 
493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494 
495 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
496  const G4MaterialCutsCouple* /*couple*/,
497  const G4DynamicParticle* particle,
498  G4double,
499  G4double)
500 {
501  if (verboseLevel > 3)
502  G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
503 
504  G4double lowLim = 0;
505  G4double highLim = 0;
506 
507  // ZF: the following line summarizes the commented part
508  if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
509  else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
510 
511  /*
512  if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
513 
514 
515  if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
516  || particle->GetDefinition() == instance->GetIon("hydrogen")
517  )
518 
519  lowLim = killBelowEnergyForA[1];
520 
521  if ( particle->GetDefinition() == instance->GetIon("alpha++")
522  || particle->GetDefinition() == instance->GetIon("alpha+")
523  || particle->GetDefinition() == instance->GetIon("helium")
524  )
525 
526  lowLim = killBelowEnergyForA[4];*/
527 
528  G4double k = particle->GetKineticEnergy();
529 
530  const G4String& particleName = particle->GetDefinition()->GetParticleName();
531 
532  // SI - the following is useless since lowLim is already defined
533  /*
534  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
535  pos1 = lowEnergyLimit.find(particleName);
536 
537  if (pos1 != lowEnergyLimit.end())
538  {
539  lowLim = pos1->second;
540  }
541  */
542 
543  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
544  pos2 = highEnergyLimit.find(particleName);
545 
546  if (pos2 != highEnergyLimit.end())highLim = pos2->second;
547 
548  if (k >= lowLim && k < highLim)
549  {
550  G4ParticleDefinition* definition = particle->GetDefinition();
551  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
552  /*
553  G4double particleMass = definition->GetPDGMass();
554  G4double totalEnergy = k + particleMass;
555  G4double pSquare = k*(totalEnergy+particleMass);
556  G4double totalMomentum = std::sqrt(pSquare);
557  */
558 
559  G4int ionizationShell = RandomSelect(k,particleName);
560 
561  // sample deexcitation
562  // here we assume that H_{2}O electronic levels are the same of Oxigen.
563  // this can be considered true with a rough 10% error in energy on K-shell,
564 
565  G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
566  G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
568  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
569 
570  if(fAtomDeexcitation) {
571  G4int Z = 8;
573 
574  if (ionizationShell <5 && ionizationShell >1)
575  {
576  as = G4AtomicShellEnumerator(4-ionizationShell);
577  }
578  else if (ionizationShell <2)
579  {
580  as = G4AtomicShellEnumerator(3);
581  }
582 
583 
584  // DEBUG
585  // if (ionizationShell == 4) {
586  //
587  // G4cout << "Z: " << Z << " as: " << as
588  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
589  // G4cout << "Press <Enter> key to continue..." << G4endl;
590  // G4cin.ignore();
591  // }
592 
593 
594  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
595  secNumberInit = fvect->size();
596  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
597  secNumberFinal = fvect->size();
598  }
599 
600 
601  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
602 
603  G4double cosTheta = 0.;
604  G4double phi = 0.;
605  RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
606 
607  G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
608  G4double dirX = sinTheta*std::cos(phi);
609  G4double dirY = sinTheta*std::sin(phi);
610  G4double dirZ = cosTheta;
611  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
612  deltaDirection.rotateUz(primaryDirection);
613 
614  // Ignored for ions on electrons
615  /*
616  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
617 
618  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
619  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
620  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
621  G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
622  finalPx /= finalMomentum;
623  finalPy /= finalMomentum;
624  finalPz /= finalMomentum;
625 
626  G4ThreeVector direction;
627  direction.set(finalPx,finalPy,finalPz);
628 
629  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
630  */
631 
633  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
634  G4double deexSecEnergy = 0;
635  for (G4int j=secNumberInit; j < secNumberFinal; j++) {
636 
637  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
638 
639  }
640 
642  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
643 
644  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
645  fvect->push_back(dp);
646 
647  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
649  ionizationShell,
650  theIncomingTrack);
651  }
652 
653  // SI - not useful since low energy of model is 0 eV
654 
655  if (k < lowLim)
656  {
660  }
661 
662 }
663 
664 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665 
666 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
667  G4double k,
668  G4int shell)
669 {
670  //-- Fast sampling method -----
671  G4double proposed_energy;
672  G4double random1;
673  G4double value_sampling;
674  G4double max1;
675 
676  do
677  {
678  proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
679 
680  max1=0.;
681 
682  for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
683  max1=RejectionFunction(particleDefinition, k, en, shell);
684 
685  random1 = G4UniformRand()*max1;
686 
687  value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
688 
689  } while(random1 > value_sampling);
690 
691  return(proposed_energy);
692 }
693 
694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
695 
696 
697 void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
698  G4double k,
699  G4double secKinetic,
700  G4double & cosTheta,
701  G4double & phi,
702  G4int shell )
703 {
704  G4double maxSecKinetic = 0.;
705  G4double maximumEnergyTransfer = 0.;
706 
707  /* if (particleDefinition == G4Proton::ProtonDefinition()
708  || particleDefinition == instance->GetIon("hydrogen"))
709  {
710  if(k/MeV < 100.)maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
711  else {
712  G4double beta2 = 1.-(1.+k*);
713  maxSecKinetic =
714  }
715  }
716 
717  if (particleDefinition == instance->GetIon("helium")
718  || particleDefinition == instance->GetIon("alpha+")
719  || particleDefinition == instance->GetIon("alpha++"))
720  {
721  maxSecKinetic = 4.* (0.511 / 3728) * k;
722  }
723 
724  if (particleDefinition == G4Carbon::Carbon())
725  {
726  maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k / 12.;
727  }*/
728 
729  // ZF. generalized & relativistic version
730 
731  if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
732  {
733  maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
734  maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
735  }
736  else
737  {
738  G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
739  G4double en_per_nucleon = k/approx_nuc_number;
740  G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
741  G4double gamma = 1./sqrt(1.-beta2);
742  maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
743  maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
744  }
745 
746  maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
747 
748  phi = twopi * G4UniformRand();
749 
750  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
751  else cosTheta = (2.*G4UniformRand())-1.;
752 
753 }
754 
755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
756 
757 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
758  G4double k,
759  G4double proposed_ws,
760  G4int ionizationLevelIndex)
761 {
762  const G4int j=ionizationLevelIndex;
763  G4double Bj_energy, alphaConst;
764  G4double Ry = 13.6*eV;
765  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
766 
767  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
768 
769  // Following values provided by M. Dingfelder (priv. comm)
770  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
771 
772  if (j == 4)
773  {
774  alphaConst = 0.66;
775  //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
776  Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
777  //---
778  }
779  else
780  {
781  alphaConst = 0.64;
782  Bj_energy = Bj[ionizationLevelIndex];
783  }
784 
785  G4double energyTransfer = proposed_ws + Bj_energy;
786  proposed_ws/=Bj_energy;
787  G4DNAGenericIonsManager *instance;
789  G4double tau = 0.;
790  G4double A_ion = 0.;
791  tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
792  A_ion = particleDefinition->GetAtomicMass();
793 
794  G4double v2;
795  G4double beta2;
796 
797  if((tau/MeV)<5.447761194e-2)
798  {
799  v2 = tau / Bj_energy;
800  beta2 = 2.*tau / electron_mass_c2;
801  }
802  // Relativistic
803  else
804  {
805  v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
806  beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
807  }
808 
809  G4double v = std::sqrt(v2);
810  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
811  G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
812  rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
813  //* (S/Bj_energy) ; Not needed anymore
814 
815  G4bool isHelium = false;
816 
817  if ( particleDefinition == G4Proton::ProtonDefinition()
818  || particleDefinition == instance->GetIon("hydrogen")
819  )
820  {
821  return(rejection_term);
822  }
823 
824  else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
825  {
826  G4double Z = particleDefinition->GetAtomicNumber();
827  G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
828  G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
829  rejection_term*=Zeffion*Zeffion;
830  }
831 
832  else if (particleDefinition == instance->GetIon("alpha++") )
833  {
834  isHelium = true;
835  slaterEffectiveCharge[0]=0.;
836  slaterEffectiveCharge[1]=0.;
837  slaterEffectiveCharge[2]=0.;
838  sCoefficient[0]=0.;
839  sCoefficient[1]=0.;
840  sCoefficient[2]=0.;
841  }
842 
843  else if (particleDefinition == instance->GetIon("alpha+") )
844  {
845  isHelium = true;
846  slaterEffectiveCharge[0]=2.0;
847  // The following values are provided by M. Dingfelder (priv. comm)
848  slaterEffectiveCharge[1]=2.0;
849  slaterEffectiveCharge[2]=2.0;
850  //
851  sCoefficient[0]=0.7;
852  sCoefficient[1]=0.15;
853  sCoefficient[2]=0.15;
854  }
855 
856  else if (particleDefinition == instance->GetIon("helium") )
857  {
858  isHelium = true;
859  slaterEffectiveCharge[0]=1.7;
860  slaterEffectiveCharge[1]=1.15;
861  slaterEffectiveCharge[2]=1.15;
862  sCoefficient[0]=0.5;
863  sCoefficient[1]=0.25;
864  sCoefficient[2]=0.25;
865  }
866 
867 // if ( particleDefinition == instance->GetIon("helium")
868 // || particleDefinition == instance->GetIon("alpha+")
869 // || particleDefinition == instance->GetIon("alpha++")
870 // )
871  if (isHelium)
872  {
873 
874  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
875 
876  zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
877  sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
878  sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
879 
880  rejection_term*= zEff * zEff;
881  }
882 
883  return (rejection_term);
884 }
885 
886 
887 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
888 
889 
890 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
891  G4double k,
892  G4int ionizationLevelIndex)
893 {
894 
895  const G4int j=ionizationLevelIndex;
896 
897  G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
898  //G4double alphaConst ;
899  G4double Bj_energy;
900 
901  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
902  // Following values provided by M. Dingfelder (priv. comm)
903 
904  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
905 
906  if (j == 4)
907  {
908  //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
909  A1 = 1.25;
910  B1 = 0.5;
911  C1 = 1.00;
912  D1 = 1.00;
913  E1 = 3.00;
914  A2 = 1.10;
915  B2 = 1.30;
916  C2 = 1.00;
917  D2 = 0.00;
918  //alphaConst = 0.66;
919  //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
920  Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
921  //---
922  }
923  else
924  {
925  //Data For Liquid Water from Dingfelder (Protons in Water)
926  A1 = 1.02;
927  B1 = 82.0;
928  C1 = 0.45;
929  D1 = -0.80;
930  E1 = 0.38;
931  A2 = 1.07;
932  //B2 = 14.6; From Ding Paper
933  // Value provided by M. Dingfelder (priv. comm)
934  B2 = 11.6;
935  //
936  C2 = 0.60;
937  D2 = 0.04;
938  //alphaConst = 0.64;
939 
940  Bj_energy = Bj[ionizationLevelIndex];
941  }
942 
943  G4double tau = 0.;
944  G4double A_ion = 0.;
945  tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
946 
947  A_ion = particle->GetAtomicMass();
948 
949  G4double v2;
950  G4double beta2;
951  if((tau/MeV)<5.447761194e-2)
952  {
953  v2 = tau / Bj_energy;
954  beta2 = 2.*tau / electron_mass_c2;
955  }
956  // Relativistic
957  else
958  {
959  v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
960  beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
961  }
962 
963  G4double v = std::sqrt(v2);
964  //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
965  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
966  G4double L2 = C2*std::pow(v,(D2));
967  G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
968  G4double H2 = (A2/v2) + (B2/(v2*v2));
969  G4double F1 = L1+H1;
970  G4double F2 = (L2*H2)/(L2+H2);
971 
972  // ZF. generalized & relativistic version
973  G4double maximumEnergy;
974 
975  //---- maximum kinetic energy , non relativistic ------
976  if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
977  {
978  maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
979  }
980  //---- relativistic -----------------------------------
981  else
982  {
983  G4double gamma = 1./sqrt(1.-beta2);
984  maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
985  (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
986  }
987 
988  //either it is transfered energy or secondary electron energy ...
989  //maximumEnergy-=Bj_energy;
990 
991  //-----------------------------------------------------
992  G4double wmax = maximumEnergy/Bj_energy;
993  G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
994  c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
995  G4double randVal = G4UniformRand();
996  G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
997  proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
998  // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
999  proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1000  proposed_ws*=Bj_energy;
1001 
1002  return(proposed_ws);
1003 
1004 }
1005 
1006 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1007 
1008 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
1009  G4double energyTransferred,
1010  G4double slaterEffectiveChg,
1011  G4double shellNumber)
1012 {
1013  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
1014  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
1015 
1016  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1017  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1018 
1019  return value;
1020 }
1021 
1022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1023 
1024 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
1025  G4double energyTransferred,
1026  G4double slaterEffectiveChg,
1027  G4double shellNumber)
1028 {
1029  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
1030  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
1031 
1032  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1033  G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1034 
1035  return value;
1036 
1037 }
1038 
1039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1040 
1041 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
1042  G4double energyTransferred,
1043  G4double slaterEffectiveChg,
1044  G4double shellNumber)
1045 {
1046  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
1047  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
1048 
1049  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1050  G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1051 
1052  return value;
1053 }
1054 
1055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1056 
1057 G4double G4DNARuddIonisationExtendedModel::R(G4double t,
1058  G4double energyTransferred,
1059  G4double slaterEffectiveChg,
1060  G4double shellNumber)
1061 {
1062  // tElectron = m_electron / m_alpha * t
1063  // Dingfelder, in Chattanooga 2005 proceedings, p 4
1064 
1065  G4double tElectron = 0.511/3728. * t;
1066  // The following values are provided by M. Dingfelder (priv. comm)
1067  G4double H = 2.*13.60569172 * eV;
1068  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1069 
1070  return value;
1071 }
1072 
1073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1074 
1075 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
1076 {
1077  // ZF Shortened
1078  G4DNAGenericIonsManager *instance;
1079  instance = G4DNAGenericIonsManager::Instance();
1080 
1081  if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
1082  {
1083  G4double value = (std::log10(k/eV)-4.2)/0.5;
1084  // The following values are provided by M. Dingfelder (priv. comm)
1085  return((0.6/(1+std::exp(value))) + 0.9);
1086  }
1087  else
1088  {
1089  return(1.);
1090  }
1091 }
1092 
1093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1094 
1095 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
1096 {
1097 
1098  G4int level = 0;
1099 
1100  // Retrieve data table corresponding to the current particle type
1101 
1102  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1103  pos = tableData.find(particle);
1104 
1105  if (pos != tableData.end())
1106  {
1107  G4DNACrossSectionDataSet* table = pos->second;
1108 
1109  if (table != 0)
1110  {
1111  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1112 
1113  const size_t n(table->NumberOfComponents());
1114  size_t i(n);
1115  G4double value = 0.;
1116 
1117  while (i>0)
1118  {
1119  i--;
1120  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1121 
1122  value += valuesBuffer[i];
1123  }
1124 
1125  value *= G4UniformRand();
1126 
1127  i = n;
1128 
1129  while (i > 0)
1130  {
1131  i--;
1132 
1133  if (valuesBuffer[i] > value)
1134  {
1135  delete[] valuesBuffer;
1136  return i;
1137  }
1138  value -= valuesBuffer[i];
1139  }
1140 
1141  if (valuesBuffer) delete[] valuesBuffer;
1142 
1143  }
1144  }
1145  else
1146  {
1147  G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
1148  FatalException,"Model not applicable to particle type.");
1149  }
1150 
1151  return level;
1152 }
1153 
1154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1155 
1156 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
1157 {
1158  G4double sigma = 0.;
1159 
1160  const G4DynamicParticle* particle = track.GetDynamicParticle();
1161  G4double k = particle->GetKineticEnergy();
1162 
1163  G4double lowLim = 0;
1164  G4double highLim = 0;
1165 
1166  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1167 
1168  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1169  pos1 = lowEnergyLimit.find(particleName);
1170 
1171  if (pos1 != lowEnergyLimit.end())
1172  {
1173  lowLim = pos1->second;
1174  }
1175 
1176  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1177  pos2 = highEnergyLimit.find(particleName);
1178 
1179  if (pos2 != highEnergyLimit.end())
1180  {
1181  highLim = pos2->second;
1182  }
1183 
1184  if (k >= lowLim && k <= highLim)
1185  {
1186  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1187  pos = tableData.find(particleName);
1188 
1189  if (pos != tableData.end())
1190  {
1191  G4DNACrossSectionDataSet* table = pos->second;
1192  if (table != 0)
1193  {
1194  sigma = table->FindValue(k);
1195  }
1196  }
1197  else
1198  {
1199  G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
1200  FatalException,"Model not applicable to particle type.");
1201  }
1202  }
1203 
1204  return sigma;
1205 }
1206 
1207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1208 
1209 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1210 {
1211  return 0;
1212 }
1213 
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:578
size_t GetIndex() const
Definition: G4Material.hh:260
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
string material
Definition: eplot.py:19
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4DNAGenericIonsManager * Instance(void)
const G4int n
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4int GetAtomicMass() const
virtual size_t NumberOfComponents(void) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define C1
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
const G4Track * GetCurrentTrack() const
const XML_Char int const XML_Char * value
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:739
G4double GetPDGCharge() const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4double bindingEnergy(G4int A, G4int Z)
G4AtomicShellEnumerator
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121
G4ParticleDefinition * GetIon(const G4String &name)