source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNADingfelderChargeDecreaseModel.cc@ 1346

Last change on this file since 1346 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 17.3 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//
[1315]26// $Id: G4DNADingfelderChargeDecreaseModel.cc,v 1.9 2010/04/06 11:00:35 sincerti Exp $
[1337]27// GEANT4 tag $Name: geant4-09-04-beta-01 $
[1058]28//
29
30#include "G4DNADingfelderChargeDecreaseModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNADingfelderChargeDecreaseModel::G4DNADingfelderChargeDecreaseModel(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 << "Dingfelder charge decrease model is constructed " << G4endl;
54 }
[1058]55}
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59G4DNADingfelderChargeDecreaseModel::~G4DNADingfelderChargeDecreaseModel()
60{}
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64void G4DNADingfelderChargeDecreaseModel::Initialise(const G4ParticleDefinition* particle,
65 const G4DataVector& /*cuts*/)
66{
67
68 if (verboseLevel > 3)
69 G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()" << G4endl;
70
71 // Energy limits
72
73 G4DNAGenericIonsManager *instance;
74 instance = G4DNAGenericIonsManager::Instance();
75 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
76 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
77 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
78
79 G4String proton;
80 G4String alphaPlusPlus;
81 G4String alphaPlus;
82
83 if (protonDef != 0)
84 {
85 proton = protonDef->GetParticleName();
[1315]86 lowEnergyLimit[proton] = 100. * eV;
[1058]87 highEnergyLimit[proton] = 10. * MeV;
88 }
89 else
90 {
91 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: proton is not defined");
92 }
93
94 if (alphaPlusPlusDef != 0)
95 {
96 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
97 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
98 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
99 }
100 else
101 {
102 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: alphaPlusPlus is not defined");
103 }
104
105 if (alphaPlusDef != 0)
106 {
107 alphaPlus = alphaPlusDef->GetParticleName();
108 lowEnergyLimit[alphaPlus] = 1. * keV;
109 highEnergyLimit[alphaPlus] = 10. * MeV;
110 }
111 else
112 {
113 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: alphaPlus is not defined");
114 }
115
116 if (particle==protonDef)
117 {
118 SetLowEnergyLimit(lowEnergyLimit[proton]);
119 SetHighEnergyLimit(highEnergyLimit[proton]);
120 }
121
122 if (particle==alphaPlusPlusDef)
123 {
124 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
125 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
126 }
127
128 if (particle==alphaPlusDef)
129 {
130 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
131 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
132 }
133
134 // Final state
135
136 //PROTON
137 f0[0][0]=1.;
138 a0[0][0]=-0.180;
139 a1[0][0]=-3.600;
140 b0[0][0]=-18.22;
141 b1[0][0]=-1.997;
142 c0[0][0]=0.215;
143 d0[0][0]=3.550;
144 x0[0][0]=3.450;
145 x1[0][0]=5.251;
146
147 numberOfPartialCrossSections[0] = 1;
148
149 //ALPHA++
150 f0[0][1]=1.; a0[0][1]=0.95;
151 a1[0][1]=-2.75;
152 b0[0][1]=-23.00;
153 c0[0][1]=0.215;
154 d0[0][1]=2.95;
155 x0[0][1]=3.50;
156
157 f0[1][1]=1.;
158 a0[1][1]=0.95;
159 a1[1][1]=-2.75;
160 b0[1][1]=-23.73;
161 c0[1][1]=0.250;
162 d0[1][1]=3.55;
163 x0[1][1]=3.72;
164
165 x1[0][1]=-1.;
166 b1[0][1]=-1.;
167
168 x1[1][1]=-1.;
169 b1[1][1]=-1.;
170
171 numberOfPartialCrossSections[1] = 2;
172
173 // ALPHA+
174 f0[0][2]=1.;
175 a0[0][2]=0.65;
176 a1[0][2]=-2.75;
177 b0[0][2]=-21.81;
178 c0[0][2]=0.232;
179 d0[0][2]=2.95;
180 x0[0][2]=3.53;
181
182 x1[0][2]=-1.;
183 b1[0][2]=-1.;
184
185 numberOfPartialCrossSections[2] = 1;
186
187 //
188
[1192]189 if( verboseLevel>0 )
190 {
191 G4cout << "Dingfelder charge decrease model is initialized " << G4endl
192 << "Energy range: "
193 << LowEnergyLimit() / keV << " keV - "
194 << HighEnergyLimit() / MeV << " MeV for "
195 << particle->GetParticleName()
196 << G4endl;
197 }
198
[1058]199 if(!isInitialised)
200 {
201 isInitialised = true;
202
203 if(pParticleChange)
204 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
205 else
206 fParticleChangeForGamma = new G4ParticleChangeForGamma();
207 }
208
209 // InitialiseElementSelectors(particle,cuts);
210
211}
212
213//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214
[1315]215G4double G4DNADingfelderChargeDecreaseModel::CrossSectionPerVolume(const G4Material* material,
[1058]216 const G4ParticleDefinition* particleDefinition,
217 G4double k,
218 G4double,
219 G4double)
220{
221 if (verboseLevel > 3)
222 G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel" << G4endl;
223
224 // Calculate total cross section for model
225
226 G4DNAGenericIonsManager *instance;
227 instance = G4DNAGenericIonsManager::Instance();
228
229 if (
230 particleDefinition != G4Proton::ProtonDefinition()
231 &&
232 particleDefinition != instance->GetIon("alpha++")
233 &&
234 particleDefinition != instance->GetIon("alpha+")
235 )
236
237 return 0;
238
239 G4double lowLim = 0;
240 G4double highLim = 0;
241 G4double crossSection = 0.;
242
[1315]243 if (material->GetName() == "G4_WATER")
[1058]244 {
245 const G4String& particleName = particleDefinition->GetParticleName();
246
247 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
248 pos1 = lowEnergyLimit.find(particleName);
249
250 if (pos1 != lowEnergyLimit.end())
251 {
252 lowLim = pos1->second;
253 }
254
255 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
256 pos2 = highEnergyLimit.find(particleName);
257
258 if (pos2 != highEnergyLimit.end())
259 {
260 highLim = pos2->second;
261 }
262
263 if (k >= lowLim && k < highLim)
264 {
265 crossSection = Sum(k,particleDefinition);
266 }
267
268 if (verboseLevel > 3)
269 {
270 G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
271 G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
[1315]272 G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
[1058]273 }
274
[1315]275 }
[1058]276
[1315]277 return crossSection*material->GetAtomicNumDensityVector()[1];
[1058]278
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282
283void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
284 const G4MaterialCutsCouple* /*couple*/,
285 const G4DynamicParticle* aDynamicParticle,
286 G4double,
287 G4double)
288{
289 if (verboseLevel > 3)
290 G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel" << G4endl;
291
292 G4double inK = aDynamicParticle->GetKineticEnergy();
293
294 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
[1315]295
296 G4double particleMass = definition->GetPDGMass();
[1058]297
298 G4int finalStateIndex = RandomSelect(inK,definition);
299
300 G4int n = NumberOfFinalStates(definition, finalStateIndex);
301 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
302 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
303
304 G4double outK = 0.;
305 if (definition==G4Proton::Proton())
306 outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) - waterBindingEnergy + outgoingParticleBindingEnergy;
307 else
[1315]308 outK = inK - n*(inK*electron_mass_c2/particleMass) - waterBindingEnergy + outgoingParticleBindingEnergy;
[1058]309
310 if (outK<0)
311 {
312 G4String message;
313 message="G4DNADingfelderChargeDecreaseModel::SampleSecondaries: final kinetic energy is below 0! Process ";
314 G4Exception(message);
315 }
316
317 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
318 fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
319
320 G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
321 aDynamicParticle->GetMomentumDirection(),
322 outK) ;
323 fvect->push_back(dp);
324
325}
326
327//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328
329G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
330 G4int finalStateIndex )
331
332{
333 if (particleDefinition == G4Proton::Proton()) return 1;
334
335 G4DNAGenericIonsManager*instance;
336 instance = G4DNAGenericIonsManager::Instance();
337
338 if (particleDefinition == instance->GetIon("alpha++") )
339 {
340 if (finalStateIndex==0) return 1;
341 return 2;
342 }
343
344 if (particleDefinition == instance->GetIon("alpha+") ) return 1;
345
346 return 0;
347}
348
349//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
350
351G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
352 G4int finalStateIndex)
353{
354 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
355
356 if (particleDefinition == G4Proton::Proton()) return instance->GetIon("hydrogen");
357
358 if (particleDefinition == instance->GetIon("alpha++") )
359 {
360 if (finalStateIndex == 0) return instance->GetIon("alpha+");
361 return instance->GetIon("helium");
362 }
363
364 if (particleDefinition == instance->GetIon("alpha+") ) return instance->GetIon("helium");
365
366 return 0;
367}
368
369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370
371G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
372 G4int finalStateIndex )
373{
374 // Ionization energy of first water shell
375 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
376 // W + 10.79 eV -> W+ + e-
377
378 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
379
380 if(particleDefinition == G4Proton::Proton()) return 10.79*eV;
381
382 if (particleDefinition == instance->GetIon("alpha++") )
383 {
384 // Binding energy for W+ -> W++ + e- 10.79 eV
385 // Binding energy for W -> W+ + e- 10.79 eV
386 //
387 // Ionization energy of first water shell
388 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
389
390 if (finalStateIndex == 0) return 10.79*eV;
391
392 return 10.79*2*eV;
393 }
394
395 if (particleDefinition == instance->GetIon("alpha+") )
396 {
397 // Binding energy for W+ -> W++ + e- 10.79 eV
398 // Binding energy for W -> W+ + e- 10.79 eV
399 //
400 // Ionization energy of first water shell
401 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
402
403 return 10.79*eV;
404 }
405
406 return 0;
407}
408
409//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410
411G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
412 G4int finalStateIndex)
413{
414 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
415
416 if(particleDefinition == G4Proton::Proton()) return 13.6*eV;
417
418 if (particleDefinition == instance->GetIon("alpha++") )
419 {
420 // Binding energy for He+ -> He++ + e- 54.509 eV
421 // Binding energy for He -> He+ + e- 24.587 eV
422
423 if (finalStateIndex==0) return 54.509*eV;
424
425 return (54.509 + 24.587)*eV;
426 }
427
428 if (particleDefinition == instance->GetIon("alpha+") )
429 {
430 // Binding energy for He+ -> He++ + e- 54.509 eV
431 // Binding energy for He -> He+ + e- 24.587 eV
432
433 return 24.587*eV;
434 }
435
436 return 0;
437}
438
439//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440
441G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k, G4int index,
442 const G4ParticleDefinition* particleDefinition)
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("alpha++")) particleTypeIndex=1;
450 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
451
452 //
453 // sigma(T) = f0 10 ^ y(log10(T/eV))
454 //
455 // / a0 x + b0 x < x0
456 // |
457 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
458 // |
459 // \ a1 x + b1 x >= x1
460 //
461 //
462 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
463 //
464 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
465 //
466 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
467 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
468 //
469
470 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
471 {
472 //
473 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
474 //
475 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
476 //
477 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
478 //
479
480 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
481 / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
482 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
483 + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex]
484 - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
485 }
486
487 G4double x(std::log10(k/eV));
488 G4double y;
489
490 if (x<x0[index][particleTypeIndex])
491 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
492 else if (x<x1[index][particleTypeIndex])
493 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex]
494 * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
495 else
496 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
497
498 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
499
500}
501
502G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
503 const G4ParticleDefinition* particleDefinition)
504{
505 G4int particleTypeIndex = 0;
506 G4DNAGenericIonsManager *instance;
507 instance = G4DNAGenericIonsManager::Instance();
508
509 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex = 0;
510 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex = 1;
511 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex = 2;
512
513 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
514 G4double* values(new G4double[n]);
515 G4double value(0);
516 G4int i = n;
517
518 while (i>0)
519 {
520 i--;
521 values[i]=PartialCrossSection(k, i, particleDefinition);
522 value+=values[i];
523 }
524
525 value*=G4UniformRand();
526
527 i=n;
528 while (i>0)
529 {
530 i--;
531
532 if (values[i]>value)
533 break;
534
535 value-=values[i];
536 }
537
538 delete[] values;
539
540 return i;
541}
542
543//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
544
545G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
546{
547 G4int particleTypeIndex = 0;
548 G4DNAGenericIonsManager* instance;
549 instance = G4DNAGenericIonsManager::Instance();
550
551 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
552 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
553 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
554
555 G4double totalCrossSection = 0.;
556
557 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
558 {
559 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
560 }
561 return totalCrossSection;
562}
563
Note: See TracBrowser for help on using the repository browser.