source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNAMillerGreenExcitationModel.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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