Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DNAMillerGreenExcitationModel.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: G4DNAMillerGreenExcitationModel.cc 95948 2016-03-03 10:40:33Z gcosmo $
27 // GEANT4 tag $Name: $
28 //
29 
30 #include "G4DNAMillerGreenExcitationModel.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4DNAChemistryManager.hh"
33 #include "G4DNAMolecularMaterial.hh"
34 #include "G4Exp.hh"
35 
36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37 
38 using namespace std;
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 
42 G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(const G4ParticleDefinition*,
43  const G4String& nam)
44  :G4VEmModel(nam),isInitialised(false)
45 {
46  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
47  fpMolWaterDensity = 0;
48 
49  nLevels=0;
50  kineticEnergyCorrection[0]=0.;
51  kineticEnergyCorrection[1]=0.;
52  kineticEnergyCorrection[2]=0.;
53  kineticEnergyCorrection[3]=0.;
54 
55  verboseLevel= 0;
56  // Verbosity scale:
57  // 0 = nothing
58  // 1 = warning for energy non-conservation
59  // 2 = details of energy budget
60  // 3 = calculation of cross sections, file openings, sampling of atoms
61  // 4 = entering in methods
62 
63  if( verboseLevel>0 )
64  {
65  G4cout << "Miller & Green excitation model is constructed " << G4endl;
66  }
67  fParticleChangeForGamma = 0;
68 
69  // Selection of stationary mode
70 
71  statCode = false;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel()
77 {}
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle,
82  const G4DataVector& /*cuts*/)
83 {
84 
85  if (verboseLevel > 3)
86  G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl;
87 
88  // Energy limits
89 
90  G4DNAGenericIonsManager *instance;
91  instance = G4DNAGenericIonsManager::Instance();
92  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
93  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
94  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
95  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
96  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
97 
98  G4String proton;
99  G4String hydrogen;
100  G4String alphaPlusPlus;
101  G4String alphaPlus;
102  G4String helium;
103 
104  // LIMITS AND CONSTANTS
105 
106  proton = protonDef->GetParticleName();
107  lowEnergyLimit[proton] = 10. * eV;
108  highEnergyLimit[proton] = 500. * keV;
109 
110  kineticEnergyCorrection[0] = 1.;
111  slaterEffectiveCharge[0][0] = 0.;
112  slaterEffectiveCharge[1][0] = 0.;
113  slaterEffectiveCharge[2][0] = 0.;
114  sCoefficient[0][0] = 0.;
115  sCoefficient[1][0] = 0.;
116  sCoefficient[2][0] = 0.;
117 
118  hydrogen = hydrogenDef->GetParticleName();
119  lowEnergyLimit[hydrogen] = 10. * eV;
120  highEnergyLimit[hydrogen] = 500. * keV;
121 
122  kineticEnergyCorrection[0] = 1.;
123  slaterEffectiveCharge[0][0] = 0.;
124  slaterEffectiveCharge[1][0] = 0.;
125  slaterEffectiveCharge[2][0] = 0.;
126  sCoefficient[0][0] = 0.;
127  sCoefficient[1][0] = 0.;
128  sCoefficient[2][0] = 0.;
129 
130  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
131  lowEnergyLimit[alphaPlusPlus] = 1. * keV;
132  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
133 
134  kineticEnergyCorrection[1] = 0.9382723/3.727417;
135  slaterEffectiveCharge[0][1]=0.;
136  slaterEffectiveCharge[1][1]=0.;
137  slaterEffectiveCharge[2][1]=0.;
138  sCoefficient[0][1]=0.;
139  sCoefficient[1][1]=0.;
140  sCoefficient[2][1]=0.;
141 
142  alphaPlus = alphaPlusDef->GetParticleName();
143  lowEnergyLimit[alphaPlus] = 1. * keV;
144  highEnergyLimit[alphaPlus] = 400. * MeV;
145 
146  kineticEnergyCorrection[2] = 0.9382723/3.727417;
147  slaterEffectiveCharge[0][2]=2.0;
148 
149  // Following values provided by M. Dingfelder
150  slaterEffectiveCharge[1][2]=2.00;
151  slaterEffectiveCharge[2][2]=2.00;
152  //
153  sCoefficient[0][2]=0.7;
154  sCoefficient[1][2]=0.15;
155  sCoefficient[2][2]=0.15;
156 
157  helium = heliumDef->GetParticleName();
158  lowEnergyLimit[helium] = 1. * keV;
159  highEnergyLimit[helium] = 400. * MeV;
160 
161  kineticEnergyCorrection[3] = 0.9382723/3.727417;
162  slaterEffectiveCharge[0][3]=1.7;
163  slaterEffectiveCharge[1][3]=1.15;
164  slaterEffectiveCharge[2][3]=1.15;
165  sCoefficient[0][3]=0.5;
166  sCoefficient[1][3]=0.25;
167  sCoefficient[2][3]=0.25;
168 
169  //
170 
171  if (particle==protonDef)
172  {
173  SetLowEnergyLimit(lowEnergyLimit[proton]);
174  SetHighEnergyLimit(highEnergyLimit[proton]);
175  }
176 
177  if (particle==hydrogenDef)
178  {
179  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
180  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
181  }
182 
183  if (particle==alphaPlusPlusDef)
184  {
185  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
186  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
187  }
188 
189  if (particle==alphaPlusDef)
190  {
191  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
192  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
193  }
194 
195  if (particle==heliumDef)
196  {
197  SetLowEnergyLimit(lowEnergyLimit[helium]);
198  SetHighEnergyLimit(highEnergyLimit[helium]);
199  }
200 
201  //
202 
203  nLevels = waterExcitation.NumberOfLevels();
204 
205  //
206  if( verboseLevel>0 )
207  {
208  G4cout << "Miller & Green excitation model is initialized " << G4endl
209  << "Energy range: "
210  << LowEnergyLimit() / eV << " eV - "
211  << HighEnergyLimit() / keV << " keV for "
212  << particle->GetParticleName()
213  << G4endl;
214  }
215 
216  // Initialize water density pointer
217  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
218 
219  if (isInitialised) { return; }
220  fParticleChangeForGamma = GetParticleChangeForGamma();
221  isInitialised = true;
222 
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 
227 G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material,
228  const G4ParticleDefinition* particleDefinition,
229  G4double k,
230  G4double,
231  G4double)
232 {
233  if (verboseLevel > 3)
234  G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl;
235 
236  // Calculate total cross section for model
237 
238  G4DNAGenericIonsManager *instance;
239  instance = G4DNAGenericIonsManager::Instance();
240 
241  if (
242  particleDefinition != G4Proton::ProtonDefinition()
243  &&
244  particleDefinition != instance->GetIon("hydrogen")
245  &&
246  particleDefinition != instance->GetIon("alpha++")
247  &&
248  particleDefinition != instance->GetIon("alpha+")
249  &&
250  particleDefinition != instance->GetIon("helium")
251  )
252 
253  return 0;
254 
255  G4double lowLim = 0;
256  G4double highLim = 0;
257  G4double crossSection = 0.;
258 
259  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
260 
261  if(waterDensity!= 0.0)
262  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
263  {
264 // G4cout << "WATER DENSITY = " << waterDensity*G4Material::GetMaterial("G4_WATER")->GetMassOfMolecule()/(g/cm3)
265 // << G4endl;
266  const G4String& particleName = particleDefinition->GetParticleName();
267 
268  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
269  pos1 = lowEnergyLimit.find(particleName);
270 
271  if (pos1 != lowEnergyLimit.end())
272  {
273  lowLim = pos1->second;
274  }
275 
276  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
277  pos2 = highEnergyLimit.find(particleName);
278 
279  if (pos2 != highEnergyLimit.end())
280  {
281  highLim = pos2->second;
282  }
283 
284  if (k >= lowLim && k < highLim)
285  {
286  crossSection = Sum(k,particleDefinition);
287 
288  // add ONE or TWO electron-water excitation for alpha+ and helium
289  /*
290  if ( particleDefinition == instance->GetIon("alpha+")
291  ||
292  particleDefinition == instance->GetIon("helium")
293  )
294  {
295 
296  G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
297  excitationXS->Initialise(G4Electron::ElectronDefinition());
298 
299  G4double sigmaExcitation=0;
300  G4double tmp =0.;
301 
302  if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation =
303  excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp)
304  /material->GetAtomicNumDensityVector()[1];
305 
306  if ( particleDefinition == instance->GetIon("alpha+") )
307  crossSection = crossSection + sigmaExcitation ;
308 
309  if ( particleDefinition == instance->GetIon("helium") )
310  crossSection = crossSection + 2*sigmaExcitation ;
311 
312  delete excitationXS;
313 
314  // Alternative excitation model
315 
316  G4DNABornExcitationModel * excitationXS = new G4DNABornExcitationModel();
317  excitationXS->Initialise(G4Electron::ElectronDefinition());
318 
319  G4double sigmaExcitation=0;
320  G4double tmp=0;
321 
322  if (k*0.511/3728 > 9*eV && k*0.511/3728 < 1*MeV ) sigmaExcitation =
323  excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp)
324  /material->GetAtomicNumDensityVector()[1];
325 
326  if ( particleDefinition == instance->GetIon("alpha+") )
327  crossSection = crossSection + sigmaExcitation ;
328 
329  if ( particleDefinition == instance->GetIon("helium") )
330  crossSection = crossSection + 2*sigmaExcitation ;
331 
332  delete excitationXS;
333 
334  }
335 */
336 
337  }
338 
339  if (verboseLevel > 2)
340  {
341  G4cout << "__________________________________" << G4endl;
342  G4cout << "G4DNAMillerGreenExcitationModel - XS INFO START" << G4endl;
343  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
344  G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
345  G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) << G4endl;
346  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
347  G4cout << "G4DNAMillerGreenExcitationModel - XS INFO END" << G4endl;
348  }
349  }
350  else
351  {
352  if (verboseLevel > 2)
353  {
354  G4cout << "MillerGreenExcitationModel : WARNING Water density is NULL" << G4endl;
355  }
356  }
357 
358  return crossSection*waterDensity;
359  // return crossSection*material->GetAtomicNumDensityVector()[1];
360 
361 }
362 
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
364 
365 void G4DNAMillerGreenExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
366  const G4MaterialCutsCouple* /*couple*/,
367  const G4DynamicParticle* aDynamicParticle,
368  G4double,
369  G4double)
370 {
371 
372  if (verboseLevel > 3)
373  G4cout << "Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" << G4endl;
374 
375  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
376 
377  G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition());
378 
379  // G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
380 
381  // Dingfelder's excitation levels
382  const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
383  G4double excitationEnergy = excitation[level];
384 
385  G4double newEnergy = 0.;
386 
387  if (!statCode) newEnergy = particleEnergy0 - excitationEnergy;
388 
389  else newEnergy = particleEnergy0;
390 
391  if (newEnergy>0)
392  {
393  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
394  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
395  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
396 
397  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
398  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
399  level,
400  theIncomingTrack);
401 
402  }
403 
404 }
405 
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407 
408 G4double G4DNAMillerGreenExcitationModel::GetPartialCrossSection(const G4Material*,
409  G4int level,
410  const G4ParticleDefinition* particleDefinition,
411  G4double kineticEnergy)
412 {
413  return PartialCrossSection(kineticEnergy, level, particleDefinition);
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
417 
418 G4double G4DNAMillerGreenExcitationModel::PartialCrossSection(G4double k, G4int excitationLevel,
419  const G4ParticleDefinition* particleDefinition)
420 {
421  // ( ( z * aj ) ^ omegaj ) * ( t - ej ) ^ nu
422  // sigma(t) = zEff^2 * sigma0 * --------------------------------------------
423  // jj ^ ( omegaj + nu ) + t ^ ( omegaj + nu )
424  //
425  // where t is the kinetic energy corrected by Helium mass over proton mass for Helium ions
426  //
427  // zEff is:
428  // 1 for protons
429  // 2 for alpha++
430  // and 2 - c1 S_1s - c2 S_2s - c3 S_2p for alpha+ and He
431  //
432  // Dingfelder et al., RPC 59, 255-275, 2000 from Miller and Green (1973)
433  // Formula (34) and Table 2
434 
435  const G4double sigma0(1.E+8 * barn);
436  const G4double nu(1.);
437  const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
438  const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
439  const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
440 
441  // Dingfelder's excitation levels
442  const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
443 
444  G4int particleTypeIndex = 0;
445  G4DNAGenericIonsManager* instance;
446  instance = G4DNAGenericIonsManager::Instance();
447 
448  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
449  if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0;
450  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
451  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
452  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3;
453 
454  G4double tCorrected;
455  tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
456 
457  // SI - added protection
458  if (tCorrected < Eliq[excitationLevel]) return 0;
459  //
460 
461  G4int z = 10;
462 
463  G4double numerator;
464  numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
465  std::pow(tCorrected - Eliq[excitationLevel], nu);
466 
467  // H case : see S. Uehara et al. IJRB 77, 2, 139-154 (2001) - section 3.3
468 
469  if (particleDefinition == instance->GetIon("hydrogen"))
470  numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
471  std::pow(tCorrected - Eliq[excitationLevel], nu);
472 
473 
474  G4double power;
475  power = omegaj[excitationLevel] + nu;
476 
477  G4double denominator;
478  denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
479 
480  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
481 
482  zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
483  sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
484  sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
485 
486  if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.;
487 
488  G4double cross = sigma0 * zEff * zEff * numerator / denominator;
489 
490 
491  return cross;
492 }
493 
494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
495 
496 G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle)
497 {
498  G4int i = nLevels;
499  G4double value = 0.;
500  std::deque<double> values;
501 
502  G4DNAGenericIonsManager *instance;
503  instance = G4DNAGenericIonsManager::Instance();
504 
505  if ( particle == instance->GetIon("alpha++") ||
506  particle == G4Proton::ProtonDefinition()||
507  particle == instance->GetIon("hydrogen") ||
508  particle == instance->GetIon("alpha+") ||
509  particle == instance->GetIon("helium")
510  )
511  {
512  while (i > 0)
513  {
514  i--;
515  G4double partial = PartialCrossSection(k,i,particle);
516  values.push_front(partial);
517  value += partial;
518  }
519 
520  value *= G4UniformRand();
521 
522  i = nLevels;
523 
524  while (i > 0)
525  {
526  i--;
527  if (values[i] > value) return i;
528  value -= values[i];
529  }
530  }
531 
532  /*
533  // add ONE or TWO electron-water excitation for alpha+ and helium
534 
535  if ( particle == instance->GetIon("alpha+")
536  ||
537  particle == instance->GetIon("helium")
538  )
539  {
540  while (i>0)
541  {
542  i--;
543 
544  G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
545  excitationXS->Initialise(G4Electron::ElectronDefinition());
546 
547  G4double sigmaExcitation=0;
548 
549  if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i);
550 
551  G4double partial = PartialCrossSection(k,i,particle);
552 
553  if (particle == instance->GetIon("alpha+")) partial = PartialCrossSection(k,i,particle) + sigmaExcitation;
554  if (particle == instance->GetIon("helium")) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation;
555 
556  values.push_front(partial);
557  value += partial;
558  delete excitationXS;
559  }
560 
561  value*=G4UniformRand();
562 
563  i=5;
564  while (i>0)
565  {
566  i--;
567 
568  if (values[i]>value) return i;
569 
570  value-=values[i];
571  }
572  }
573 */
574 
575  return 0;
576 }
577 
578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
579 
580 G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle)
581 {
582  G4double totalCrossSection = 0.;
583 
584  for (G4int i=0; i<nLevels; i++)
585  {
586  totalCrossSection += PartialCrossSection(k,i,particle);
587  }
588  return totalCrossSection;
589 }
590 
591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
592 
593 G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t,
594  G4double energyTransferred,
595  G4double _slaterEffectiveCharge,
596  G4double shellNumber)
597 {
598  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
599  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
600 
601  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
602  G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
603 
604  return value;
605 }
606 
607 
608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
609 
610 G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t,
611  G4double energyTransferred,
612  G4double _slaterEffectiveCharge,
613  G4double shellNumber)
614 {
615  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
616  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
617 
618  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
619  G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
620 
621  return value;
622 
623 }
624 
625 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
626 
627 G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t,
628  G4double energyTransferred,
629  G4double _slaterEffectiveCharge,
630  G4double shellNumber)
631 {
632  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
633  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
634 
635  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
636  G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
637 
638  return value;
639 }
640 
641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
642 
643 G4double G4DNAMillerGreenExcitationModel::R(G4double t,
644  G4double energyTransferred,
645  G4double _slaterEffectiveCharge,
646  G4double shellNumber)
647 {
648  // tElectron = m_electron / m_alpha * t
649  // Dingfelder, in Chattanooga 2005 proceedings, p 4
650 
651  G4double tElectron = 0.511/3728. * t;
652 
653  // The following is provided by M. Dingfelder
654  G4double H = 2.*13.60569172 * eV;
655  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber);
656 
657  return value;
658 }
659 
G4Exp.hh
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4DNAMillerGreenExcitationModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNAMillerGreenExcitationModel.cc:227
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel
G4DNAMillerGreenExcitationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAMillerGreenExcitationModel")
Definition: G4DNAMillerGreenExcitationModel.cc:42
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DNAMolecularMaterial.hh
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
G4DNAMillerGreenExcitationModel::GetPartialCrossSection
virtual G4double GetPartialCrossSection(const G4Material *, G4int, const G4ParticleDefinition *, G4double)
Definition: G4DNAMillerGreenExcitationModel.cc:408
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4DNAMolecularMaterial::GetNumMolPerVolTableFor
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
Definition: G4DNAMolecularMaterial.cc:437
G4DNAMillerGreenExcitationModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNAMillerGreenExcitationModel.hh:77
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
cm
static constexpr double cm
Definition: G4SIunits.hh:119
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4DNAGenericIonsManager::Instance
static G4DNAGenericIonsManager * Instance(void)
Definition: G4DNAGenericIonsManager.cc:36
G4DNAChemistryManager.hh
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:117
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:75
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:572
eExcitedMolecule
Definition: G4DNAChemistryManager.hh:71
G4ParticleChangeForGamma::GetCurrentTrack
const G4Track * GetCurrentTrack() const
Definition: G4ParticleChangeForGamma.hh:154
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4DNAMillerGreenExcitationModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNAMillerGreenExcitationModel.cc:365
G4double
double G4double
Definition: G4Types.hh:76
G4DNAMillerGreenExcitationModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4DNAMillerGreenExcitationModel.cc:81
G4SystemOfUnits.hh
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:731
barn
static constexpr double barn
Definition: G4SIunits.hh:105
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4DNAMillerGreenExcitationModel.hh
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
G4DNAGenericIonsManager::GetIon
G4ParticleDefinition * GetIon(const G4String &name)
Definition: G4DNAGenericIonsManager.cc:46