Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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