Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4DNADingfelderChargeDecreaseModel.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$
27//
28
31#include "G4SystemOfUnits.hh"
33
34//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35
36using namespace std;
37
38//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39
41 const G4String& nam)
42 :G4VEmModel(nam),isInitialised(false)
43{
44 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
45 fpMolWaterDensity = 0;
46 numberOfPartialCrossSections[0] = 0;
47 numberOfPartialCrossSections[1] = 0;
48 numberOfPartialCrossSections[2] = 0;
49
50 verboseLevel= 0;
51 // Verbosity scale:
52 // 0 = nothing
53 // 1 = warning for energy non-conservation
54 // 2 = details of energy budget
55 // 3 = calculation of cross sections, file openings, sampling of atoms
56 // 4 = entering in methods
57
58 if( verboseLevel>0 )
59 {
60 G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
61 }
63}
64
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
68{}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
73 const G4DataVector& /*cuts*/)
74{
75
76 if (verboseLevel > 3)
77 G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()" << G4endl;
78
79 // Energy limits
80
84 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
85 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
86
87 G4String proton;
88 G4String alphaPlusPlus;
89 G4String alphaPlus;
90
91 // LIMITS
92
93 proton = protonDef->GetParticleName();
94 lowEnergyLimit[proton] = 100. * eV;
95 highEnergyLimit[proton] = 100. * MeV;
96
97 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
98 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
99 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
100
101 alphaPlus = alphaPlusDef->GetParticleName();
102 lowEnergyLimit[alphaPlus] = 1. * keV;
103 highEnergyLimit[alphaPlus] = 400. * MeV;
104
105 //
106
107 if (particle==protonDef)
108 {
109 SetLowEnergyLimit(lowEnergyLimit[proton]);
110 SetHighEnergyLimit(highEnergyLimit[proton]);
111 }
112
113 if (particle==alphaPlusPlusDef)
114 {
115 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
116 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
117 }
118
119 if (particle==alphaPlusDef)
120 {
121 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
122 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
123 }
124
125 // Final state
126
127 //PROTON
128 f0[0][0]=1.;
129 a0[0][0]=-0.180;
130 a1[0][0]=-3.600;
131 b0[0][0]=-18.22;
132 b1[0][0]=-1.997;
133 c0[0][0]=0.215;
134 d0[0][0]=3.550;
135 x0[0][0]=3.450;
136 x1[0][0]=5.251;
137
138 numberOfPartialCrossSections[0] = 1;
139
140 //ALPHA++
141 f0[0][1]=1.; a0[0][1]=0.95;
142 a1[0][1]=-2.75;
143 b0[0][1]=-23.00;
144 c0[0][1]=0.215;
145 d0[0][1]=2.95;
146 x0[0][1]=3.50;
147
148 f0[1][1]=1.;
149 a0[1][1]=0.95;
150 a1[1][1]=-2.75;
151 b0[1][1]=-23.73;
152 c0[1][1]=0.250;
153 d0[1][1]=3.55;
154 x0[1][1]=3.72;
155
156 x1[0][1]=-1.;
157 b1[0][1]=-1.;
158
159 x1[1][1]=-1.;
160 b1[1][1]=-1.;
161
162 numberOfPartialCrossSections[1] = 2;
163
164 // ALPHA+
165 f0[0][2]=1.;
166 a0[0][2]=0.65;
167 a1[0][2]=-2.75;
168 b0[0][2]=-21.81;
169 c0[0][2]=0.232;
170 d0[0][2]=2.95;
171 x0[0][2]=3.53;
172
173 x1[0][2]=-1.;
174 b1[0][2]=-1.;
175
176 numberOfPartialCrossSections[2] = 1;
177
178 //
179
180 if( verboseLevel>0 )
181 {
182 G4cout << "Dingfelder charge decrease model is initialized " << G4endl
183 << "Energy range: "
184 << LowEnergyLimit() / keV << " keV - "
185 << HighEnergyLimit() / MeV << " MeV for "
186 << particle->GetParticleName()
187 << G4endl;
188 }
189
190 // Initialize water density pointer
192
193 if (isInitialised) { return; }
195 isInitialised = true;
196
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200
202 const G4ParticleDefinition* particleDefinition,
203 G4double k,
204 G4double,
205 G4double)
206{
207 if (verboseLevel > 3)
208 G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel" << G4endl;
209
210 // Calculate total cross section for model
211
212 G4DNAGenericIonsManager *instance;
214
215 if (
216 particleDefinition != G4Proton::ProtonDefinition()
217 &&
218 particleDefinition != instance->GetIon("alpha++")
219 &&
220 particleDefinition != instance->GetIon("alpha+")
221 )
222
223 return 0;
224
225 G4double lowLim = 0;
226 G4double highLim = 0;
227 G4double crossSection = 0.;
228
229 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
230
231 if(waterDensity!= 0.0)
232 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
233 {
234 const G4String& particleName = particleDefinition->GetParticleName();
235
236 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
237 pos1 = lowEnergyLimit.find(particleName);
238
239 if (pos1 != lowEnergyLimit.end())
240 {
241 lowLim = pos1->second;
242 }
243
244 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
245 pos2 = highEnergyLimit.find(particleName);
246
247 if (pos2 != highEnergyLimit.end())
248 {
249 highLim = pos2->second;
250 }
251
252 if (k >= lowLim && k < highLim)
253 {
254 crossSection = Sum(k,particleDefinition);
255 }
256
257 if (verboseLevel > 2)
258 {
259 G4cout << "_______________________________________" << G4endl;
260 G4cout << "°°° G4DNADingfelderChargeDecreaeModel" << G4endl;
261 G4cout << "---> Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
262 G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
263 G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*
264 waterDensity/(1./cm) << G4endl;
265// material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
266 }
267 }
268
269 return crossSection*waterDensity;
270// return crossSection*material->GetAtomicNumDensityVector()[1];
271
272}
273
274//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275
276void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
277 const G4MaterialCutsCouple* /*couple*/,
278 const G4DynamicParticle* aDynamicParticle,
279 G4double,
280 G4double)
281{
282 if (verboseLevel > 3)
283 G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel" << G4endl;
284
285 G4double inK = aDynamicParticle->GetKineticEnergy();
286
287 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
288
289 G4double particleMass = definition->GetPDGMass();
290
291 G4int finalStateIndex = RandomSelect(inK,definition);
292
293 G4int n = NumberOfFinalStates(definition, finalStateIndex);
294 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
295 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
296
297 G4double outK = 0.;
298 if (definition==G4Proton::Proton())
299 outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) - waterBindingEnergy + outgoingParticleBindingEnergy;
300 else
301 outK = inK - n*(inK*electron_mass_c2/particleMass) - waterBindingEnergy + outgoingParticleBindingEnergy;
302
303 if (outK<0)
304 {
305 G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
306 FatalException,"Final kinetic energy is negative.");
307 }
308
311
312 G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
313 aDynamicParticle->GetMomentumDirection(),
314 outK) ;
315 fvect->push_back(dp);
316
317}
318
319//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320
321G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
322 G4int finalStateIndex )
323
324{
325 if (particleDefinition == G4Proton::Proton()) return 1;
326
329
330 if (particleDefinition == instance->GetIon("alpha++") )
331 {
332 if (finalStateIndex==0) return 1;
333 return 2;
334 }
335
336 if (particleDefinition == instance->GetIon("alpha+") ) return 1;
337
338 return 0;
339}
340
341//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342
343G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
344 G4int finalStateIndex)
345{
347
348 if (particleDefinition == G4Proton::Proton()) return instance->GetIon("hydrogen");
349
350 if (particleDefinition == instance->GetIon("alpha++") )
351 {
352 if (finalStateIndex == 0) return instance->GetIon("alpha+");
353 return instance->GetIon("helium");
354 }
355
356 if (particleDefinition == instance->GetIon("alpha+") ) return instance->GetIon("helium");
357
358 return 0;
359}
360
361//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
362
363G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
364 G4int finalStateIndex )
365{
366 // Ionization energy of first water shell
367 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
368 // W + 10.79 eV -> W+ + e-
369
371
372 if(particleDefinition == G4Proton::Proton()) return 10.79*eV;
373
374 if (particleDefinition == instance->GetIon("alpha++") )
375 {
376 // Binding energy for W+ -> W++ + e- 10.79 eV
377 // Binding energy for W -> W+ + e- 10.79 eV
378 //
379 // Ionization energy of first water shell
380 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
381
382 if (finalStateIndex == 0) return 10.79*eV;
383
384 return 10.79*2*eV;
385 }
386
387 if (particleDefinition == instance->GetIon("alpha+") )
388 {
389 // Binding energy for W+ -> W++ + e- 10.79 eV
390 // Binding energy for W -> W+ + e- 10.79 eV
391 //
392 // Ionization energy of first water shell
393 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
394
395 return 10.79*eV;
396 }
397
398 return 0;
399}
400
401//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402
403G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
404 G4int finalStateIndex)
405{
407
408 if(particleDefinition == G4Proton::Proton()) return 13.6*eV;
409
410 if (particleDefinition == instance->GetIon("alpha++") )
411 {
412 // Binding energy for He+ -> He++ + e- 54.509 eV
413 // Binding energy for He -> He+ + e- 24.587 eV
414
415 if (finalStateIndex==0) return 54.509*eV;
416
417 return (54.509 + 24.587)*eV;
418 }
419
420 if (particleDefinition == instance->GetIon("alpha+") )
421 {
422 // Binding energy for He+ -> He++ + e- 54.509 eV
423 // Binding energy for He -> He+ + e- 24.587 eV
424
425 return 24.587*eV;
426 }
427
428 return 0;
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432
433G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k, G4int index,
434 const G4ParticleDefinition* particleDefinition)
435{
436 G4int particleTypeIndex = 0;
437 G4DNAGenericIonsManager* instance;
439
440 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
441 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
442 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
443
444 //
445 // sigma(T) = f0 10 ^ y(log10(T/eV))
446 //
447 // / a0 x + b0 x < x0
448 // |
449 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
450 // |
451 // \ a1 x + b1 x >= x1
452 //
453 //
454 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
455 //
456 // 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)
457 //
458 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
459 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
460 //
461
462 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
463 {
464 //
465 // 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)
466 //
467 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
468 //
469 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
470 //
471
472 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
473 / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
474 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
475 + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex]
476 - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
477 }
478
479 G4double x(std::log10(k/eV));
480 G4double y;
481
482 if (x<x0[index][particleTypeIndex])
483 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
484 else if (x<x1[index][particleTypeIndex])
485 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex]
486 * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
487 else
488 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
489
490 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
491
492}
493
494G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
495 const G4ParticleDefinition* particleDefinition)
496{
497 G4int particleTypeIndex = 0;
498 G4DNAGenericIonsManager *instance;
500
501 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex = 0;
502 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex = 1;
503 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex = 2;
504
505 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
506 G4double* values(new G4double[n]);
507 G4double value(0);
508 G4int i = n;
509
510 while (i>0)
511 {
512 i--;
513 values[i]=PartialCrossSection(k, i, particleDefinition);
514 value+=values[i];
515 }
516
517 value*=G4UniformRand();
518
519 i=n;
520 while (i>0)
521 {
522 i--;
523
524 if (values[i]>value)
525 break;
526
527 value-=values[i];
528 }
529
530 delete[] values;
531
532 return i;
533}
534
535//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
536
537G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
538{
539 G4int particleTypeIndex = 0;
540 G4DNAGenericIonsManager* instance;
542
543 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
544 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
545 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
546
547 G4double totalCrossSection = 0.;
548
549 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
550 {
551 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
552 }
553 return totalCrossSection;
554}
555
@ FatalException
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeDecreaseModel")
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
size_t GetIndex() const
Definition: G4Material.hh:261
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
static G4Proton * Proton()
Definition: G4Proton.cc:93
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41