Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNADingfelderChargeIncreaseModel.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
47 numberOfPartialCrossSections[0]=0;
48 numberOfPartialCrossSections[1]=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 increase 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 G4DNADingfelderChargeIncreaseModel::Initialise()" << G4endl;
78
79 // Energy limits
80
83 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
84 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
85 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
86
87 G4String hydrogen;
88 G4String alphaPlus;
89 G4String helium;
90
91 // LIMITS
92
93 hydrogen = hydrogenDef->GetParticleName();
94 lowEnergyLimit[hydrogen] = 100. * eV;
95 highEnergyLimit[hydrogen] = 100. * MeV;
96
97 alphaPlus = alphaPlusDef->GetParticleName();
98 lowEnergyLimit[alphaPlus] = 1. * keV;
99 highEnergyLimit[alphaPlus] = 400. * MeV;
100
101 helium = heliumDef->GetParticleName();
102 lowEnergyLimit[helium] = 1. * keV;
103 highEnergyLimit[helium] = 400. * MeV;
104
105 //
106
107 if (particle==hydrogenDef)
108 {
109 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
110 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
111 }
112
113 if (particle==alphaPlusDef)
114 {
115 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
116 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
117 }
118
119 if (particle==heliumDef)
120 {
121 SetLowEnergyLimit(lowEnergyLimit[helium]);
122 SetHighEnergyLimit(highEnergyLimit[helium]);
123 }
124
125 // Final state
126
127 //ALPHA+
128
129 f0[0][0]=1.;
130 a0[0][0]=2.25;
131 a1[0][0]=-0.75;
132 b0[0][0]=-32.10;
133 c0[0][0]=0.600;
134 d0[0][0]=2.40;
135 x0[0][0]=4.60;
136
137 x1[0][0]=-1.;
138 b1[0][0]=-1.;
139
140 numberOfPartialCrossSections[0]=1;
141
142 //HELIUM
143
144 f0[0][1]=1.;
145 a0[0][1]=2.25;
146 a1[0][1]=-0.75;
147 b0[0][1]=-30.93;
148 c0[0][1]=0.590;
149 d0[0][1]=2.35;
150 x0[0][1]=4.29;
151
152 f0[1][1]=1.;
153 a0[1][1]=2.25;
154 a1[1][1]=-0.75;
155 b0[1][1]=-32.61;
156 c0[1][1]=0.435;
157 d0[1][1]=2.70;
158 x0[1][1]=4.45;
159
160 x1[0][1]=-1.;
161 b1[0][1]=-1.;
162
163 x1[1][1]=-1.;
164 b1[1][1]=-1.;
165
166 numberOfPartialCrossSections[1]=2;
167
168 //
169
170 if( verboseLevel>0 )
171 {
172 G4cout << "Dingfelder charge increase model is initialized " << G4endl
173 << "Energy range: "
174 << LowEnergyLimit() / keV << " keV - "
175 << HighEnergyLimit() / MeV << " MeV for "
176 << particle->GetParticleName()
177 << G4endl;
178 }
179
180 // Initialize water density pointer
182
183 if (isInitialised) { return; }
185 isInitialised = true;
186}
187
188//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189
191 const G4ParticleDefinition* particleDefinition,
192 G4double k,
193 G4double,
194 G4double)
195{
196 if (verboseLevel > 3)
197 G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" << G4endl;
198
199 // Calculate total cross section for model
200
201 G4DNAGenericIonsManager *instance;
203
204 if (
205 particleDefinition != instance->GetIon("hydrogen")
206 &&
207 particleDefinition != instance->GetIon("alpha+")
208 &&
209 particleDefinition != instance->GetIon("helium")
210 )
211
212 return 0;
213
214 G4double lowLim = 0;
215 G4double highLim = 0;
216 G4double totalCrossSection = 0.;
217
218 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
219
220 if(waterDensity!= 0.0)
221 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
222 {
223 const G4String& particleName = particleDefinition->GetParticleName();
224
225 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
226 pos1 = lowEnergyLimit.find(particleName);
227
228 if (pos1 != lowEnergyLimit.end())
229 {
230 lowLim = pos1->second;
231 }
232
233 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
234 pos2 = highEnergyLimit.find(particleName);
235
236 if (pos2 != highEnergyLimit.end())
237 {
238 highLim = pos2->second;
239 }
240
241 if (k >= lowLim && k < highLim)
242 {
243 //HYDROGEN
244 if (particleDefinition == instance->GetIon("hydrogen"))
245 {
246 const G4double aa = 2.835;
247 const G4double bb = 0.310;
248 const G4double cc = 2.100;
249 const G4double dd = 0.760;
250 const G4double fac = 1.0e-18;
251 const G4double rr = 13.606 * eV;
252
253 G4double t = k / (proton_mass_c2/electron_mass_c2);
254 G4double x = t / rr;
255 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
256 G4double sigmal = temp * cc * (std::pow(x,dd));
257 G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
258 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
259 }
260 else
261 {
262 totalCrossSection = Sum(k,particleDefinition);
263 }
264 }
265
266 if (verboseLevel > 2)
267 {
268 G4cout << "__________________________________" << G4endl;
269 G4cout << "°°° G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl;
270 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl;
271 G4cout << "°°° Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl;
272 G4cout << "°°° Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl;
273 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
274 G4cout << "°°° G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl;
275 }
276
277 }
278
279 return totalCrossSection*waterDensity;
280// return totalCrossSection*material->GetAtomicNumDensityVector()[1];
281
282}
283
284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285
286void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
287 const G4MaterialCutsCouple* /*couple*/,
288 const G4DynamicParticle* aDynamicParticle,
289 G4double,
290 G4double)
291{
292 if (verboseLevel > 3)
293 G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" << G4endl;
294
297
298 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
299
300 G4double particleMass = definition->GetPDGMass();
301
302 G4double inK = aDynamicParticle->GetKineticEnergy();
303
304 G4int finalStateIndex = RandomSelect(inK,definition);
305
306 G4int n = NumberOfFinalStates(definition,finalStateIndex);
307
308 G4double outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
309
310 G4DNAGenericIonsManager* instance;
312
313 G4double electronK;
314 if (definition == instance->GetIon("hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
315 else electronK = inK*electron_mass_c2/(particleMass);
316
317 if (outK<0)
318 {
319 G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004",
320 FatalException,"Final kinetic energy is negative.");
321 }
322
324
325 (OutgoingParticleDefinition(definition,finalStateIndex),
326 aDynamicParticle->GetMomentumDirection(),
327 outK) ;
328
329 fvect->push_back(dp);
330
331 n = n - 1;
332
333 while (n>0)
334 {
335 n--;
336 fvect->push_back(new G4DynamicParticle
337 (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) );
338 }
339
340}
341
342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
343
344G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
345 G4int finalStateIndex )
346
347{
348 G4DNAGenericIonsManager* instance;
350 if (particleDefinition == instance->GetIon("hydrogen")) return 2;
351 if (particleDefinition == instance->GetIon("alpha+")) return 2;
352
353 if (particleDefinition == instance->GetIon("helium"))
354 { if (finalStateIndex==0) return 2;
355 return 3;
356 }
357 return 0;
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361
362G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
363 G4int finalStateIndex)
364{
366 if (particleDefinition == instance->GetIon("hydrogen")) return G4Proton::Proton();
367 if (particleDefinition == instance->GetIon("alpha+")) return instance->GetIon("alpha++");
368
369 if (particleDefinition == instance->GetIon("helium"))
370 {
371 if (finalStateIndex==0) return instance->GetIon("alpha+");
372 return instance->GetIon("alpha++");
373 }
374 return 0;
375}
376
377//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378
379G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
380 G4int finalStateIndex )
381{
383 if(particleDefinition == instance->GetIon("hydrogen")) return 13.6*eV;
384
385 if(particleDefinition == instance->GetIon("alpha+"))
386 {
387 // Binding energy for He+ -> He++ + e- 54.509 eV
388 // Binding energy for He -> He+ + e- 24.587 eV
389 return 54.509*eV;
390 }
391
392 if(particleDefinition == instance->GetIon("helium"))
393 {
394 // Binding energy for He+ -> He++ + e- 54.509 eV
395 // Binding energy for He -> He+ + e- 24.587 eV
396
397 if (finalStateIndex==0) return 24.587*eV;
398 return (54.509 + 24.587)*eV;
399 }
400
401 return 0;
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
405
406//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407
408G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(G4double k, G4int index,
409 const G4ParticleDefinition* particleDefinition)
410{
411
412 G4int particleTypeIndex = 0;
413 G4DNAGenericIonsManager *instance;
415
416 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
417 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
418
419 //
420 // sigma(T) = f0 10 ^ y(log10(T/eV))
421 //
422 // / a0 x + b0 x < x0
423 // |
424 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
425 // |
426 // \ a1 x + b1 x >= x1
427 //
428 //
429 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
430 //
431 // 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)
432 //
433 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
434 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
435 //
436
437 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
438 {
439 //
440 // 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)
441 //
442 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
443 //
444 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
445 //
446
447 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
448 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex] - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
449 }
450
451 G4double x(std::log10(k/eV));
452 G4double y;
453
454 if (x<x0[index][particleTypeIndex])
455 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
456 else if (x<x1[index][particleTypeIndex])
457 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
458 else
459 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
460
461 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
462
463}
464
465//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466
467G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(G4double k,
468 const G4ParticleDefinition* particleDefinition)
469{
470
471 G4int particleTypeIndex = 0;
472 G4DNAGenericIonsManager *instance;
474
475 if (particleDefinition == instance->GetIon("hydrogen")) return 0;
476 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
477 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
478
479 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
480 G4double* values(new G4double[n]);
481 G4double value = 0;
482 G4int i = n;
483
484 while (i>0)
485 {
486 i--;
487 values[i]=PartialCrossSection(k, i, particleDefinition);
488 value+=values[i];
489 }
490
491 value*=G4UniformRand();
492
493 i=n;
494 while (i>0)
495 {
496 i--;
497
498 if (values[i]>value)
499 break;
500
501 value-=values[i];
502 }
503
504 delete[] values;
505
506 return i;
507}
508
509//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
510
511G4double G4DNADingfelderChargeIncreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
512{
513 G4int particleTypeIndex = 0;
514 G4DNAGenericIonsManager *instance;
516
517 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
518 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
519
520 G4double totalCrossSection = 0.;
521
522 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
523 {
524 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
525 }
526 return totalCrossSection;
527}
@ 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 void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeIncreaseModel")
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
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
static G4Electron * Electron()
Definition: G4Electron.cc:94
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 * 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