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
G4BraggIonModel.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// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32//
33// File name: G4BraggIonModel
34//
35// Author: Vladimir Ivanchenko
36//
37// Creation date: 13.10.2004
38//
39// Modifications:
40// 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko)
41// 29-11-05 Do not use G4Alpha class (V.Ivantchenko)
42// 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
43// 25-04-06 Add stopping data from ASTAR (V.Ivanchenko)
44// 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko)
45// 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
46// CorrectionsAlongStep needed for ions(V.Ivanchenko)
47//
48
49// Class Description:
50//
51// Implementation of energy loss and delta-electron production by
52// slow charged heavy particles
53
54// -------------------------------------------------------------------
55//
56
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61#include "G4BraggIonModel.hh"
62#include "G4PhysicalConstants.hh"
63#include "G4SystemOfUnits.hh"
64#include "Randomize.hh"
65#include "G4Electron.hh"
66#include "G4ParticleChangeForLoss.hh"
67#include "G4LossTableManager.hh"
68#include "G4EmCorrections.hh"
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72using namespace std;
73
74G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p,
75 const G4String& nam)
76 : G4VEmModel(nam),
77 corr(0),
78 particle(0),
79 fParticleChange(0),
80 currentMaterial(0),
81 iMolecula(-1),
82 iASTAR(-1),
83 isIon(false),
84 isInitialised(false)
85{
86 SetHighEnergyLimit(2.0*MeV);
87
88 HeMass = 3.727417*GeV;
89 rateMassHe2p = HeMass/proton_mass_c2;
90 lowestKinEnergy = 1.0*keV/rateMassHe2p;
91 massFactor = 1000.*amu_c2/HeMass;
92 theZieglerFactor = eV*cm2*1.0e-15;
93 theElectron = G4Electron::Electron();
94 corrFactor = 1.0;
95 if(p) { SetParticle(p); }
96 else { SetParticle(theElectron); }
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101G4BraggIonModel::~G4BraggIonModel()
102{}
103
104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105
106void G4BraggIonModel::Initialise(const G4ParticleDefinition* p,
107 const G4DataVector&)
108{
109 if(p != particle) { SetParticle(p); }
110
111 corrFactor = chargeSquare;
112
113 // always false before the run
114 SetDeexcitationFlag(false);
115
116 if(!isInitialised) {
117 isInitialised = true;
118
119 G4String pname = particle->GetParticleName();
120 if(particle->GetParticleType() == "nucleus" &&
121 pname != "deuteron" && pname != "triton" &&
122 pname != "alpha+" && pname != "helium" &&
123 pname != "hydrogen") { isIon = true; }
124
125 corr = G4LossTableManager::Instance()->EmCorrections();
126
127 fParticleChange = GetParticleChangeForLoss();
128 }
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
134 const G4Material* mat,
135 G4double kineticEnergy)
136{
137 //G4cout << "G4BraggIonModel::GetChargeSquareRatio e= " << kineticEnergy << G4endl;
138 // this method is called only for ions
139 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
140 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
141 return corrFactor;
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145
146G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p,
147 const G4Material* mat,
148 G4double kineticEnergy)
149{
150 //G4cout << "G4BraggIonModel::GetParticleCharge e= " << kineticEnergy << G4endl;
151 // this method is called only for ions
152 return corr->GetParticleCharge(p,mat,kineticEnergy);
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156
157G4double G4BraggIonModel::ComputeCrossSectionPerElectron(
158 const G4ParticleDefinition* p,
159 G4double kineticEnergy,
160 G4double cutEnergy,
161 G4double maxKinEnergy)
162{
163 G4double cross = 0.0;
164 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
165 G4double maxEnergy = std::min(tmax,maxKinEnergy);
166 if(cutEnergy < tmax) {
167
168 G4double energy = kineticEnergy + mass;
169 G4double energy2 = energy*energy;
170 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
171 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
172
173 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
174 }
175 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
176 // << " tmax= " << tmax << " cross= " << cross << G4endl;
177
178 return cross;
179}
180
181//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182
183G4double G4BraggIonModel::ComputeCrossSectionPerAtom(
184 const G4ParticleDefinition* p,
185 G4double kineticEnergy,
186 G4double Z, G4double,
187 G4double cutEnergy,
188 G4double maxEnergy)
189{
190 G4double cross = Z*ComputeCrossSectionPerElectron
191 (p,kineticEnergy,cutEnergy,maxEnergy);
192 return cross;
193}
194
195//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196
197G4double G4BraggIonModel::CrossSectionPerVolume(
198 const G4Material* material,
199 const G4ParticleDefinition* p,
200 G4double kineticEnergy,
201 G4double cutEnergy,
202 G4double maxEnergy)
203{
204 G4double eDensity = material->GetElectronDensity();
205 G4double cross = eDensity*ComputeCrossSectionPerElectron
206 (p,kineticEnergy,cutEnergy,maxEnergy);
207 return cross;
208}
209
210//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211
212G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material,
213 const G4ParticleDefinition* p,
214 G4double kineticEnergy,
215 G4double cutEnergy)
216{
217 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
218 G4double tmin = min(cutEnergy, tmax);
219 G4double tkin = kineticEnergy/massRate;
220 G4double dedx = 0.0;
221
222 if(tkin < lowestKinEnergy) {
223 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
224 } else {
225 dedx = DEDX(material, tkin);
226 }
227
228 if (cutEnergy < tmax) {
229
230 G4double tau = kineticEnergy/mass;
231 G4double gam = tau + 1.0;
232 G4double bg2 = tau * (tau+2.0);
233 G4double beta2 = bg2/(gam*gam);
234 G4double x = tmin/tmax;
235
236 dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
237 * (material->GetElectronDensity())/beta2;
238 }
239
240 // now compute the total ionization loss
241
242 if (dedx < 0.0) dedx = 0.0 ;
243
244 dedx *= chargeSquare;
245
246 //G4cout << " tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = "
247 // << dedx*gram/(MeV*cm2*material->GetDensity())
248 // << " q2 = " << chargeSquare << G4endl;
249
250 return dedx;
251}
252
253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254
255void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
256 const G4DynamicParticle* dp,
257 G4double& eloss,
258 G4double&,
259 G4double /*length*/)
260{
261 // this method is called only for ions
262 const G4ParticleDefinition* p = dp->GetDefinition();
263 const G4Material* mat = couple->GetMaterial();
264 G4double preKinEnergy = dp->GetKineticEnergy();
265 G4double e = preKinEnergy - eloss*0.5;
266 if(e < 0.0) { e = preKinEnergy*0.5; }
267
268 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
269 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
270 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
271 eloss *= qfactor;
272
273 //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " << e
274 // << " qfactor= " << qfactor << " " << p->GetParticleName() <<G4endl;
275}
276
277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278
279void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
280 const G4MaterialCutsCouple*,
281 const G4DynamicParticle* dp,
282 G4double xmin,
283 G4double maxEnergy)
284{
285 G4double tmax = MaxSecondaryKinEnergy(dp);
286 G4double xmax = std::min(tmax, maxEnergy);
287 if(xmin >= xmax) { return; }
288
289 G4double kineticEnergy = dp->GetKineticEnergy();
290 G4double energy = kineticEnergy + mass;
291 G4double energy2 = energy*energy;
292 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
293 G4double grej = 1.0;
294 G4double deltaKinEnergy, f;
295
296 G4ThreeVector direction = dp->GetMomentumDirection();
297
298 // sampling follows ...
299 do {
300 G4double q = G4UniformRand();
301 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
302
303 f = 1.0 - beta2*deltaKinEnergy/tmax;
304
305 if(f > grej) {
306 G4cout << "G4BraggIonModel::SampleSecondary Warning! "
307 << "Majorant " << grej << " < "
308 << f << " for e= " << deltaKinEnergy
309 << G4endl;
310 }
311
312 } while( grej*G4UniformRand() >= f );
313
314 G4double deltaMomentum =
315 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
316 G4double totMomentum = energy*sqrt(beta2);
317 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
318 (deltaMomentum * totMomentum);
319 if(cost > 1.0) { cost = 1.0; }
320 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
321
322 G4double phi = twopi * G4UniformRand() ;
323
324 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
325 deltaDirection.rotateUz(direction);
326
327 // create G4DynamicParticle object for delta ray
328 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
329 deltaKinEnergy);
330
331 vdp->push_back(delta);
332
333 // Change kinematics of primary particle
334 kineticEnergy -= deltaKinEnergy;
335 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
336 finalP = finalP.unit();
337
338 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
339 fParticleChange->SetProposedMomentumDirection(finalP);
340}
341
342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
343
344G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
345 G4double kinEnergy)
346{
347 if(pd != particle) { SetParticle(pd); }
348 G4double tau = kinEnergy/mass;
349 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
350 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
351 return tmax;
352}
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355
356G4bool G4BraggIonModel::HasMaterial(const G4Material* material)
357{
358 G4String chFormula = material->GetChemicalFormula();
359 if("" == chFormula) { return false; }
360
361 // ICRU Report N49, 1993. Ziegler model for He.
362 const size_t numberOfMolecula = 11;
363 const G4String molName[numberOfMolecula] = {
364 "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate",
365 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate",
366 "Polysterene", "SiO_2", "NaI", "H_2O",
367 "Graphite" };
368 const G4int idxASTAR[numberOfMolecula] = {
369 17, 19, 33, 51,
370 52, 54,
371 56, 62, 43, 71,
372 13};
373
374 // Search for the material in the table
375 for (size_t i=0; i<numberOfMolecula; ++i) {
376 if (chFormula == molName[i]) {
377 iMolecula = -1;
378 iASTAR = idxASTAR[i];
379 return true;
380 }
381 }
382 return false ;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386
387G4double G4BraggIonModel::StoppingPower(const G4Material* material,
388 G4double kineticEnergy)
389{
390 G4double ionloss = 0.0 ;
391
392 if (iMolecula >= 0) {
393
394 // The data and the fit from:
395 // ICRU Report N49, 1993. Ziegler's model for alpha
396 // He energy in internal units of parametrisation formula (MeV)
397
398 G4double T = kineticEnergy*rateMassHe2p/MeV ;
399
400 static G4double a[11][5] = {
401 {9.43672, 0.54398, 84.341, 1.3705, 57.422},
402 {67.1503, 0.41409, 404.512, 148.97, 20.99},
403 {5.11203, 0.453, 36.718, 50.6, 28.058},
404 {61.793, 0.48445, 361.537, 57.889, 50.674},
405 {7.83464, 0.49804, 160.452, 3.192, 0.71922},
406 {19.729, 0.52153, 162.341, 58.35, 25.668},
407 {26.4648, 0.50112, 188.913, 30.079, 16.509},
408 {7.8655, 0.5205, 63.96, 51.32, 67.775},
409 {8.8965, 0.5148, 339.36, 1.7205, 0.70423},
410 {2.959, 0.53255, 34.247, 60.655, 15.153},
411 {3.80133, 0.41590, 12.9966, 117.83, 242.28} };
412
413 static G4double atomicWeight[11] = {
414 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
415 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
416
417 G4int i = iMolecula;
418
419 // Free electron gas model
420 if ( T < 0.001 ) {
421 G4double slow = a[i][0] ;
422 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
423 * a[i][2]*1000.0 ;
424 ionloss = slow*shigh / (slow + shigh) ;
425 ionloss *= sqrt(T*1000.0) ;
426
427 // Main parametrisation
428 } else {
429 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
430 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
431 ionloss = slow*shigh / (slow + shigh) ;
432 /*
433 G4cout << "## " << i << ". T= " << T << " slow= " << slow
434 << " a0= " << a[i][0] << " a1= " << a[i][1]
435 << " shigh= " << shigh
436 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
437 << G4endl;
438 */
439 }
440 if ( ionloss < 0.0) ionloss = 0.0 ;
441
442 // He effective charge
443 G4double aa = atomicWeight[iMolecula];
444 ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
445
446 // pure material (normally not the case for this function)
447 } else if(1 == (material->GetNumberOfElements())) {
448 G4double z = material->GetZ() ;
449 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
450 }
451
452 return ionloss;
453}
454
455//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456
457G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
458 G4double kineticEnergy) const
459{
460 G4double ionloss ;
461 G4int i = G4int(z)-1 ; // index of atom
462 if(i < 0) i = 0 ;
463 if(i > 91) i = 91 ;
464
465 // The data and the fit from:
466 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
467 // Proton kinetic energy for parametrisation (keV/amu)
468
469 // He energy in internal units of parametrisation formula (MeV)
470 G4double T = kineticEnergy*rateMassHe2p/MeV ;
471
472 static G4double a[92][5] = {
473 {0.35485, 0.6456, 6.01525, 20.8933, 4.3515
474 },{ 0.58, 0.59, 6.3, 130.0, 44.07
475 },{ 1.42, 0.49, 12.25, 32.0, 9.161
476 },{ 2.206, 0.51, 15.32, 0.25, 8.995 //Be Ziegler77
477 // },{ 2.1895, 0.47183,7.2362, 134.30, 197.96 //Be from ICRU
478 },{ 3.691, 0.4128, 18.48, 50.72, 9.0
479 },{ 3.83523, 0.42993,12.6125, 227.41, 188.97
480 },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768
481 },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382
482 },{ 1.533, 0.531, 40.44, 18.41, 2.718
483 },{ 2.303, 0.4861, 37.01, 37.96, 5.092
484 // Z= 11-20
485 },{ 9.894, 0.3081, 23.65, 0.384, 92.93
486 },{ 4.3, 0.47, 34.3, 3.3, 12.74
487 },{ 2.5, 0.625, 45.7, 0.1, 4.359
488 },{ 2.1, 0.65, 49.34, 1.788, 4.133
489 },{ 1.729, 0.6562, 53.41, 2.405, 3.845
490 },{ 1.402, 0.6791, 58.98, 3.528, 3.211
491 },{ 1.117, 0.7044, 69.69, 3.705, 2.156
492 },{ 2.291, 0.6284, 73.88, 4.478, 2.066
493 },{ 8.554, 0.3817, 83.61, 11.84, 1.875
494 },{ 6.297, 0.4622, 65.39, 10.14, 5.036
495 // Z= 21-30
496 },{ 5.307, 0.4918, 61.74, 12.4, 6.665
497 },{ 4.71, 0.5087, 65.28, 8.806, 5.948
498 },{ 6.151, 0.4524, 83.0, 18.31, 2.71
499 },{ 6.57, 0.4322, 84.76, 15.53, 2.779
500 },{ 5.738, 0.4492, 84.6, 14.18, 3.101
501 },{ 5.013, 0.4707, 85.8, 16.55, 3.211
502 },{ 4.32, 0.4947, 76.14, 10.85, 5.441
503 },{ 4.652, 0.4571, 80.73, 22.0, 4.952
504 },{ 3.114, 0.5236, 76.67, 7.62, 6.385
505 },{ 3.114, 0.5236, 76.67, 7.62, 7.502
506 // Z= 31-40
507 },{ 3.114, 0.5236, 76.67, 7.62, 8.514
508 },{ 5.746, 0.4662, 79.24, 1.185, 7.993
509 },{ 2.792, 0.6346, 106.1, 0.2986, 2.331
510 },{ 4.667, 0.5095, 124.3, 2.102, 1.667
511 },{ 2.44, 0.6346, 105.0, 0.83, 2.851
512 },{ 1.413, 0.7377, 147.9, 1.466, 1.016
513 },{ 11.72, 0.3826, 102.8, 9.231, 4.371
514 },{ 7.126, 0.4804, 119.3, 5.784, 2.454
515 },{ 11.61, 0.3955, 146.7, 7.031, 1.423
516 },{ 10.99, 0.41, 163.9, 7.1, 1.052
517 // Z= 41-50
518 },{ 9.241, 0.4275, 163.1, 7.954, 1.102
519 },{ 9.276, 0.418, 157.1, 8.038, 1.29
520 },{ 3.999, 0.6152, 97.6, 1.297, 5.792
521 },{ 4.306, 0.5658, 97.99, 5.514, 5.754
522 },{ 3.615, 0.6197, 86.26, 0.333, 8.689
523 },{ 5.8, 0.49, 147.2, 6.903, 1.289
524 },{ 5.6, 0.49, 130.0, 10.0, 2.844
525 },{ 3.55, 0.6068, 124.7, 1.112, 3.119
526 },{ 3.6, 0.62, 105.8, 0.1692, 6.026
527 },{ 5.4, 0.53, 103.1, 3.931, 7.767
528 // Z= 51-60
529 },{ 3.97, 0.6459, 131.8, 0.2233, 2.723
530 },{ 3.65, 0.64, 126.8, 0.6834, 3.411
531 },{ 3.118, 0.6519, 164.9, 1.208, 1.51
532 },{ 3.949, 0.6209, 200.5, 1.878, 0.9126
533 },{ 14.4, 0.3923, 152.5, 8.354, 2.597
534 },{ 10.99, 0.4599, 138.4, 4.811, 3.726
535 },{ 16.6, 0.3773, 224.1, 6.28, 0.9121
536 },{ 10.54, 0.4533, 159.3, 4.832, 2.529
537 },{ 10.33, 0.4502, 162.0, 5.132, 2.444
538 },{ 10.15, 0.4471, 165.6, 5.378, 2.328
539 // Z= 61-70
540 },{ 9.976, 0.4439, 168.0, 5.721, 2.258
541 },{ 9.804, 0.4408, 176.2, 5.675, 1.997
542 },{ 14.22, 0.363, 228.4, 7.024, 1.016
543 },{ 9.952, 0.4318, 233.5, 5.065, 0.9244
544 },{ 9.272, 0.4345, 210.0, 4.911, 1.258
545 },{ 10.13, 0.4146, 225.7, 5.525, 1.055
546 },{ 8.949, 0.4304, 213.3, 5.071, 1.221
547 },{ 11.94, 0.3783, 247.2, 6.655, 0.849
548 },{ 8.472, 0.4405, 195.5, 4.051, 1.604
549 },{ 8.301, 0.4399, 203.7, 3.667, 1.459
550 // Z= 71-80
551 },{ 6.567, 0.4858, 193.0, 2.65, 1.66
552 },{ 5.951, 0.5016, 196.1, 2.662, 1.589
553 },{ 7.495, 0.4523, 251.4, 3.433, 0.8619
554 },{ 6.335, 0.4825, 255.1, 2.834, 0.8228
555 },{ 4.314, 0.5558, 214.8, 2.354, 1.263
556 },{ 4.02, 0.5681, 219.9, 2.402, 1.191
557 },{ 3.836, 0.5765, 210.2, 2.742, 1.305
558 },{ 4.68, 0.5247, 244.7, 2.749, 0.8962
559 },{ 2.892, 0.6204, 208.6, 2.415, 1.416 //Au Z77
560 // },{ 3.223, 0.5883, 232.7, 2.954, 1.05 //Au ICRU
561 },{ 2.892, 0.6204, 208.6, 2.415, 1.416
562 // Z= 81-90
563 },{ 4.728, 0.5522, 217.0, 3.091, 1.386
564 },{ 6.18, 0.52, 170.0, 4.0, 3.224
565 },{ 9.0, 0.47, 198.0, 3.8, 2.032
566 },{ 2.324, 0.6997, 216.0, 1.599, 1.399
567 },{ 1.961, 0.7286, 223.0, 1.621, 1.296
568 },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507
569 },{ 10.31, 0.4613, 261.2, 4.738, 0.9899
570 },{ 7.962, 0.519, 235.7, 4.347, 1.313
571 },{ 6.227, 0.5645, 231.9, 3.961, 1.379
572 },{ 5.246, 0.5947, 228.6, 4.027, 1.432
573 // Z= 91-92
574 },{ 5.408, 0.5811, 235.7, 3.961, 1.358
575 },{ 5.218, 0.5828, 245.0, 3.838, 1.25}
576 };
577
578 // Free electron gas model
579 if ( T < 0.001 ) {
580 G4double slow = a[i][0] ;
581 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
582 * a[i][2]*1000.0 ;
583 ionloss = slow*shigh / (slow + shigh) ;
584 ionloss *= sqrt(T*1000.0) ;
585
586 // Main parametrisation
587 } else {
588 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
589 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
590 ionloss = slow*shigh / (slow + shigh) ;
591 /*
592 G4cout << "## " << i << ". T= " << T << " slow= " << slow
593 << " a0= " << a[i][0] << " a1= " << a[i][1]
594 << " shigh= " << shigh
595 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
596 << G4endl;
597 */
598 }
599 if ( ionloss < 0.0) { ionloss = 0.0; }
600
601 // He effective charge
602 ionloss /= HeEffChargeSquare(z, T);
603
604 return ionloss;
605}
606
607//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
608
609G4double G4BraggIonModel::DEDX(const G4Material* material,
610 G4double kineticEnergy)
611{
612 G4double eloss = 0.0;
613 // check DB
614 if(material != currentMaterial) {
615 currentMaterial = material;
616 iASTAR = -1;
617 iMolecula = -1;
618 if( !HasMaterial(material) ) { iASTAR = astar.GetIndex(material); }
619 }
620
621 const G4int numberOfElements = material->GetNumberOfElements();
622 const G4double* theAtomicNumDensityVector =
623 material->GetAtomicNumDensityVector();
624
625 if( iASTAR >= 0 ) {
626 G4double T = kineticEnergy*rateMassHe2p;
627 return astar.GetElectronicDEDX(iASTAR, T)*material->GetDensity()/
628 HeEffChargeSquare(astar.GetEffectiveZ(iASTAR), T/MeV);
629
630 } else if(iMolecula >= 0) {
631
632 eloss = StoppingPower(material, kineticEnergy)*
633 material->GetDensity()/amu;
634
635 // pure material
636 } else if(1 == numberOfElements) {
637
638 G4double z = material->GetZ();
639 eloss = ElectronicStoppingPower(z, kineticEnergy)
640 * (material->GetTotNbOfAtomsPerVolume());
641
642 // Brugg's rule calculation
643 } else {
644 const G4ElementVector* theElementVector =
645 material->GetElementVector() ;
646
647 // loop for the elements in the material
648 for (G4int i=0; i<numberOfElements; i++)
649 {
650 const G4Element* element = (*theElementVector)[i] ;
651 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
652 * theAtomicNumDensityVector[i];
653 }
654 }
655 return eloss*theZieglerFactor;
656}
657
658//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
659
660G4double G4BraggIonModel::HeEffChargeSquare(G4double z,
661 G4double kinEnergyHeInMeV) const
662{
663 // The aproximation of He effective charge from:
664 // J.F.Ziegler, J.P. Biersack, U. Littmark
665 // The Stopping and Range of Ions in Matter,
666 // Vol.1, Pergamon Press, 1985
667
668 static G4double c[6] = {0.2865, 0.1266, -0.001429,
669 0.02402,-0.01135, 0.001475};
670
671 G4double e = std::max(0.0,std::log(kinEnergyHeInMeV*massFactor));
672 G4double x = c[0] ;
673 G4double y = 1.0 ;
674 for (G4int i=1; i<6; i++) {
675 y *= e ;
676 x += y * c[i] ;
677 }
678
679 G4double w = 7.6 - e ;
680 w = 1.0 + (0.007 + 0.00005*z) * exp( -w*w ) ;
681 w = 4.0 * (1.0 - exp(-x)) * w * w ;
682
683 return w;
684}
685
686//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
687
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4ASTARStopping::GetEffectiveZ
G4double GetEffectiveZ(G4int idx)
Definition: G4ASTARStopping.hh:99
G4ASTARStopping::GetIndex
G4int GetIndex(const G4Material *)
Definition: G4ASTARStopping.cc:77
G4ASTARStopping::GetElectronicDEDX
G4double GetElectronicDEDX(G4int idx, G4double energy)
Definition: G4ASTARStopping.cc:93
G4BraggIonModel::MaxSecondaryEnergy
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4BraggIonModel.cc:344
G4BraggIonModel::ComputeCrossSectionPerElectron
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Definition: G4BraggIonModel.cc:157
G4BraggIonModel::GetChargeSquareRatio
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4BraggIonModel.cc:133
G4BraggIonModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4BraggIonModel.cc:106
G4BraggIonModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
Definition: G4BraggIonModel.cc:212
G4BraggIonModel::CorrectionsAlongStep
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4BraggIonModel.cc:255
G4BraggIonModel::GetParticleCharge
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy)
Definition: G4BraggIonModel.cc:146
G4BraggIonModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
Definition: G4BraggIonModel.cc:183
G4BraggIonModel::G4BraggIonModel
G4BraggIonModel(const G4ParticleDefinition *p=0, const G4String &nam="BraggIon")
Definition: G4BraggIonModel.cc:74
G4BraggIonModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4BraggIonModel.cc:279
G4BraggIonModel::~G4BraggIonModel
virtual ~G4BraggIonModel()
Definition: G4BraggIonModel.cc:101
G4BraggIonModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Definition: G4BraggIonModel.cc:197
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4EmCorrections::EffectiveChargeSquareRatio
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.hh:335
G4EmCorrections::GetParticleCharge
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.hh:327
G4EmCorrections::EffectiveChargeCorrection
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:742
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
G4LossTableManager::EmCorrections
G4EmCorrections * EmCorrections()
Definition: G4LossTableManager.cc:1078
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:179
G4Material::GetChemicalFormula
const G4String & GetChemicalFormula() const
Definition: G4Material.hh:178
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:208
G4Material::GetZ
G4double GetZ() const
Definition: G4Material.cc:604
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetAtomicNumDensityVector
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:215
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:216
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:136
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115
G4VEmFluctuationModel::SetParticleAndCharge
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
Definition: G4VEmFluctuationModel.cc:69
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4VEmModel::GetModelOfFluctuations
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:501
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
G4VEmModel::MaxSecondaryKinEnergy
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:399
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95