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

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

geant4 tag 9.4
File size: 19.9 KB
RevLine 
[1058]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//
[1340]26// $Id: G4DNAMillerGreenExcitationModel.cc,v 1.11 2010/10/08 08:53:17 sincerti Exp $
[1347]27// GEANT4 tag $Name: geant4-09-04-ref-00 $
[1058]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
[1192]51 if( verboseLevel>0 )
52 {
53 G4cout << "Miller & Green excitation model is constructed " << G4endl;
54 }
[1058]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();
[1340]76 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
[1058]77 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
78 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
79 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
80
81 G4String proton;
[1340]82 G4String hydrogen;
[1058]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
[1340]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 }
[1058]125 if (alphaPlusPlusDef != 0)
126 {
127 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
128 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
[1340]129 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
[1058]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;
[1340]148 highEnergyLimit[alphaPlus] = 400. * MeV;
[1058]149
150 kineticEnergyCorrection[2] = 0.9382723/3.727417;
151 slaterEffectiveCharge[0][2]=2.0;
[1340]152
153// Following values provided by M. Dingfelder
154 slaterEffectiveCharge[1][2]=2.00;
155 slaterEffectiveCharge[2][2]=2.00;
156//
[1058]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;
[1340]170 highEnergyLimit[helium] = 400. * MeV;
[1058]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;
[1340]179
[1058]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
[1340]192 if (particle==hydrogenDef)
193 {
194 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
195 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
196 }
197
[1058]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 //
[1192]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 }
[1058]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 &&
[1340]264 particleDefinition != instance->GetIon("hydrogen")
265 &&
[1058]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
[1315]279 if (material->GetName() == "G4_WATER")
[1058]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
[1340]307/*
[1058]308 if ( particleDefinition == instance->GetIon("alpha+")
309 ||
310 particleDefinition == instance->GetIon("helium")
311 )
312 {
[1340]313
[1058]314 G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
[1340]315 excitationXS->Initialise(G4Electron::ElectronDefinition());
[1058]316
317 G4double sigmaExcitation=0;
318 G4double tmp =0.;
319
[1315]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];
[1058]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;
[1340]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
[1058]352 }
[1340]353*/
[1058]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;
[1315]361 G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
[1058]362 }
363
[1315]364 }
[1058]365
[1315]366 return crossSection*material->GetAtomicNumDensityVector()[1];
[1058]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
[1340]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
[1058]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
[1340]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
[1058]431 G4int particleTypeIndex = 0;
432 G4DNAGenericIonsManager* instance;
433 instance = G4DNAGenericIonsManager::Instance();
434
435 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
[1340]436 if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0;
[1058]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
[1340]445 if (tCorrected < Eliq[excitationLevel]) return 0;
[1058]446 //
447
448 G4int z = 10;
449
450 G4double numerator;
451 numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
[1340]452 std::pow(tCorrected - Eliq[excitationLevel], nu);
[1058]453
[1340]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
[1058]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
[1340]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.) );
[1058]472
[1340]473 if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.;
474
[1058]475 G4double cross = sigma0 * zEff * zEff * numerator / denominator;
476
[1340]477
[1058]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++") ||
[1340]493 particle == G4Proton::ProtonDefinition()||
494 particle == instance->GetIon("hydrogen") ||
495 particle == instance->GetIon("alpha+") ||
496 particle == instance->GetIon("helium")
497 )
[1058]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
[1340]519/*
[1058]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();
[1340]532 excitationXS->Initialise(G4Electron::ElectronDefinition());
[1058]533
534 G4double sigmaExcitation=0;
535
[1315]536 if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i);
[1058]537
538 G4double partial = PartialCrossSection(k,i,particle);
[1340]539
[1058]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;
[1340]542
[1058]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 }
[1340]560*/
[1058]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
[1340]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
[1058]644 return value;
645}
646
Note: See TracBrowser for help on using the repository browser.