| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4BraggModel
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 03.01.2002
37 //
38 // Modifications:
39 //
40 // 04-12-02 Fix problem of G4DynamicParticle c
41 // 23-12-02 Change interface in order to move
42 // 27-01-03 Make models region aware (V.Ivanch
43 // 13-02-03 Add name (V.Ivanchenko)
44 // 04-06-03 Fix compilation warnings (V.Ivanch
45 // 12-09-04 Add lowestKinEnergy and change ord
46 // 11-04-05 Major optimisation of internal int
47 // 16-06-05 Fix problem of chemical formula (V
48 // 15-02-06 ComputeCrossSectionPerElectron, Co
49 // 25-04-06 Add stopping data from PSTAR (V.Iv
50 // 12-08-08 Added methods GetParticleCharge, G
51 // CorrectionsAlongStep needed for io
52
53 // Class Description:
54 //
55 // Implementation of energy loss and delta-ele
56 // slow charged heavy particles
57
58 // -------------------------------------------
59 //
60
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 #include "G4BraggModel.hh"
66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh"
68 #include "Randomize.hh"
69 #include "G4Electron.hh"
70 #include "G4ParticleChangeForLoss.hh"
71 #include "G4LossTableManager.hh"
72 #include "G4EmCorrections.hh"
73 #include "G4EmParameters.hh"
74 #include "G4DeltaAngle.hh"
75 #include "G4ICRU90StoppingData.hh"
76 #include "G4NistManager.hh"
77 #include "G4Log.hh"
78 #include "G4Exp.hh"
79 #include "G4AutoLock.hh"
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4ICRU90StoppingData* G4BraggModel::fICRU90 =
84 G4PSTARStopping* G4BraggModel::fPSTAR = nullpt
85
86 namespace
87 {
88 G4Mutex ionMutex = G4MUTEX_INITIALIZER;
89 }
90
91 G4BraggModel::G4BraggModel(const G4ParticleDef
92 : G4VEmModel(nam)
93 {
94 SetHighEnergyLimit(2.0*CLHEP::MeV);
95
96 lowestKinEnergy = 0.25*CLHEP::keV;
97 theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e
98 theElectron = G4Electron::Electron();
99 expStopPower125 = 0.0;
100
101 corr = G4LossTableManager::Instance()->EmCor
102 if(nullptr != p) { SetParticle(p); }
103 else { SetParticle(theElectron); }
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 G4BraggModel::~G4BraggModel()
109 {
110 if(isFirst) {
111 delete fPSTAR;
112 fPSTAR = nullptr;
113 }
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117
118 void G4BraggModel::Initialise(const G4Particle
119 const G4DataVect
120 {
121 if(p != particle) { SetParticle(p); }
122
123 // always false before the run
124 SetDeexcitationFlag(false);
125
126 // initialise data only once
127 if(nullptr == fPSTAR) {
128 G4AutoLock l(&ionMutex);
129 if(nullptr == fPSTAR) {
130 isFirst = true;
131 fPSTAR = new G4PSTARStopping();
132 if(G4EmParameters::Instance()->UseICRU90
133 fICRU90 = G4NistManager::Instance()->GetICRU
134 }
135 }
136 l.unlock();
137 }
138 if(isFirst) {
139 if(nullptr != fICRU90) { fICRU90->Initiali
140 fPSTAR->Initialise();
141 }
142
143 if(nullptr == fParticleChange) {
144
145 if(UseAngularGeneratorFlag() && !GetAngula
146 SetAngularDistribution(new G4DeltaAngle(
147 }
148 G4String pname = particle->GetParticleName
149 if(particle->GetParticleType() == "nucleus
150 pname != "deuteron" && pname != "triton
151 pname != "alpha+" && pname != "helium
152 pname != "hydrogen") { isIon = true; }
153
154 fParticleChange = GetParticleChangeForLoss
155 }
156 }
157
158 //....oooOO0OOooo........oooOO0OOooo........oo
159
160 void G4BraggModel::SetParticle(const G4Particl
161 {
162 particle = p;
163 mass = particle->GetPDGMass();
164 spin = particle->GetPDGSpin();
165 G4double q = particle->GetPDGCharge()/CLHEP:
166 chargeSquare = q*q;
167 massRate = mass/CLHEP::proton_mass_c2;
168 ratio = CLHEP::electron_mass_c2/mass;
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oo
172
173 G4double G4BraggModel::GetChargeSquareRatio(co
174 co
175 G4
176 {
177 // this method is called only for ions
178 chargeSquare = corr->EffectiveChargeSquareRa
179 return chargeSquare;
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oo
183
184 G4double G4BraggModel::MinEnergyCut(const G4Pa
185 const G4Ma
186 {
187 return couple->GetMaterial()->GetIonisation(
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oo
191
192 G4double G4BraggModel::GetParticleCharge(const
193 const
194 G4dou
195 {
196 // this method is called only for ions, so n
197 return corr->GetParticleCharge(p,mat,kinetic
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oo
201
202 G4double G4BraggModel::ComputeCrossSectionPerE
203 con
204
205
206
207 {
208 G4double cross = 0.0;
209 const G4double tmax = MaxSecondaryEnerg
210 const G4double maxEnergy = std::min(tmax, ma
211 const G4double cutEnergy = std::max(cut, low
212 if(cutEnergy < maxEnergy) {
213
214 const G4double energy = kineticEnergy + m
215 const G4double energy2 = energy*energy;
216 const G4double beta2 = kineticEnergy*(ki
217 cross = (maxEnergy - cutEnergy)/(cutEnergy
218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
219
220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy
221
222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar
223 }
224 // G4cout << "BR: e= " << kineticEnergy <<
225 // << " tmax= " << tmax << " cross=
226 return cross;
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oo
230
231 G4double
232 G4BraggModel::ComputeCrossSectionPerAtom(const
233 G4dou
234 G4dou
235 G4dou
236 G4dou
237 {
238 return
239 Z*ComputeCrossSectionPerElectron(p,kinetic
240 }
241
242 //....oooOO0OOooo........oooOO0OOooo........oo
243
244 G4double G4BraggModel::CrossSectionPerVolume(c
245 c
246 G
247 G
248 G
249 {
250 return material->GetElectronDensity()
251 *ComputeCrossSectionPerElectron(p,kineticE
252 }
253
254 //....oooOO0OOooo........oooOO0OOooo........oo
255
256 G4double G4BraggModel::ComputeDEDXPerVolume(co
257 co
258 G4
259 G4
260 {
261 const G4double tmax = MaxSecondaryEnergy(p,
262 const G4double tkin = kinEnergy/massRate;
263 const G4double cutEnergy = std::max(cut, low
264 G4double dedx = 0.0;
265
266 // tkin is the scaled energy to proton
267 if (tkin < lowestKinEnergy) {
268 dedx = DEDX(material, lowestKinEnergy)*std
269 } else {
270 dedx = DEDX(material, tkin);
271
272 // correction for low cut value using Beth
273 if (cutEnergy < tmax) {
274 const G4double tau = kinEnergy/mass;
275 const G4double x = cutEnergy/tmax;
276
277 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/
278 CLHEP::twopi_mc2_rcl2 * material->GetElectro
279 }
280 }
281 dedx = std::max(dedx, 0.0) * chargeSquare;
282
283 //G4cout << "E(MeV)= " << tkin/MeV << " dedx
284 // << " " << material->GetName() <<
285 return dedx;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 void G4BraggModel::SampleSecondaries(std::vect
291 const G4M
292 const G4D
293 G4double
294 G4double
295 {
296 const G4double tmax = MaxSecondaryKinEnergy(
297 const G4double xmax = std::min(tmax, maxEner
298 const G4double xmin = std::max(lowestKinEner
299 if(xmin >= xmax) { return; }
300
301 G4double kineticEnergy = dp->GetKineticEnerg
302 const G4double energy = kineticEnergy + mas
303 const G4double energy2 = energy*energy;
304 const G4double beta2 = kineticEnergy*(kine
305 const G4double grej = 1.0;
306 G4double deltaKinEnergy, f;
307
308 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
309 G4double rndm[2];
310
311 // sampling follows ...
312 do {
313 rndmEngineMod->flatArray(2, rndm);
314 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rn
315
316 f = 1.0 - beta2*deltaKinEnergy/tmax;
317
318 if(f > grej) {
319 G4cout << "G4BraggModel::SampleSeconda
320 << "Majorant " << grej << " < "
321 << f << " for e= " << deltaKinE
322 << G4endl;
323 }
324
325 // Loop checking, 03-Aug-2015, Vladimir Iv
326 } while( grej*rndm[1] >= f );
327
328 G4ThreeVector deltaDirection;
329
330 if(UseAngularGeneratorFlag()) {
331 const G4Material* mat = couple->GetMateri
332 G4int Z = SelectRandomAtomNumber(mat);
333
334 deltaDirection =
335 GetAngularDistribution()->SampleDirectio
336
337 } else {
338
339 G4double deltaMomentum =
340 std::sqrt(deltaKinEnergy * (deltaKinEner
341 G4double cost = deltaKinEnergy * (energy +
342 (deltaMomentum * dp->GetTotalMomentum())
343 if(cost > 1.0) { cost = 1.0; }
344 G4double sint = std::sqrt((1.0 - cost)*(1.
345
346 G4double phi = twopi*rndmEngineMod->flat()
347
348 deltaDirection.set(sint*std::cos(phi),sint
349 deltaDirection.rotateUz(dp->GetMomentumDir
350 }
351
352 // create G4DynamicParticle object for delta
353 auto delta = new G4DynamicParticle(theElectr
354
355 // Change kinematics of primary particle
356 kineticEnergy -= deltaKinEnergy;
357 G4ThreeVector finalP = dp->GetMomentum() - d
358 finalP = finalP.unit();
359
360 fParticleChange->SetProposedKineticEnergy(ki
361 fParticleChange->SetProposedMomentumDirectio
362
363 vdp->push_back(delta);
364 }
365
366 //....oooOO0OOooo........oooOO0OOooo........oo
367
368 G4double G4BraggModel::MaxSecondaryEnergy(cons
369 G4do
370 {
371 if(pd != particle) { SetParticle(pd); }
372 G4double tau = kinEnergy/mass;
373 G4double tmax = 2.0*electron_mass_c2*tau*(ta
374 (1. + 2.0*(tau + 1.)*ratio +
375 return tmax;
376 }
377
378 //....oooOO0OOooo........oooOO0OOooo........oo
379
380 void G4BraggModel::HasMaterial(const G4Materia
381 {
382 const G4String& chFormula = mat->GetChemical
383 if(chFormula.empty()) { return; }
384
385 // ICRU Report N49, 1993. Power's model for
386 static const G4int numberOfMolecula = 11;
387 static const G4String molName[numberOfMolecu
388 "Al_2O_3", "CO_2",
389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol
390 "C_3H_8", "SiO_2",
391 "H_2O-Gas", "Graphite" } ;
392
393 // Search for the material in the table
394 for (G4int i=0; i<numberOfMolecula; ++i) {
395 if (chFormula == molName[i]) {
396 iMolecula = i;
397 return;
398 }
399 }
400 return;
401 }
402
403 //....oooOO0OOooo........oooOO0OOooo........oo
404
405 G4double G4BraggModel::StoppingPower(const G4M
406 G4double
407 {
408 G4double ionloss = 0.0 ;
409
410 if (iMolecula >= 0) {
411
412 // The data and the fit from:
413 // ICRU Report N49, 1993. Ziegler's model
414 // Proton kinetic energy for parametrisati
415
416 G4double T = kineticEnergy/(keV*protonMass
417
418 static const G4float a[11][5] = {
419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f
420 {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f
421 {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f
422 {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f
423 {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f
424 {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f
425 {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f
426 {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f
427 {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f
428 {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f
429 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f
430
431 static const G4float atomicWeight[11] = {
432 101.96128f, 44.0098f, 16.0426f, 28.0536f,
433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18
434
435 if ( T < 10.0 ) {
436 ionloss = ((G4double)(a[iMolecula][0]))
437
438 } else if ( T < 10000.0 ) {
439 G4double x1 = (G4double)(a[iMolecula][1]
440 G4double x2 = (G4double)(a[iMolecula][2]
441 G4double x3 = (G4double)(a[iMolecula][3]
442 G4double x4 = (G4double)(a[iMolecula][4]
443 G4double slow = x1 * G4Exp(G4Log(T)* 0.
444 G4double shigh = G4Log( 1.0 + x3/T + x4
445 ionloss = slow*shigh / (slow + shigh) ;
446 }
447
448 ionloss = std::max(ionloss, 0.0);
449 if ( 10 == iMolecula ) {
450 static const G4double invLog10 = 1.0/G4L
451
452 if (T < 100.0) {
453 ionloss *= (1.0+0.023+0.0066*G4Log(T)*
454 }
455 else if (T < 700.0) {
456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99
457 }
458 else if (T < 10000.0) {
459 ionloss *=(1.0+0.089-0.0248*G4Log(700.
460 }
461 }
462 ionloss /= (G4double)atomicWeight[iMolecul
463
464 // pure material (normally not the case for
465 } else if(1 == (material->GetNumberOfElement
466 G4double z = material->GetZ() ;
467 ionloss = ElectronicStoppingPower( z, kine
468 }
469
470 return ionloss;
471 }
472
473 //....oooOO0OOooo........oooOO0OOooo........oo
474
475 G4double G4BraggModel::ElectronicStoppingPower
476
477 {
478 G4double ionloss ;
479 G4int i = std::min(std::max(G4lrint(z)-1,0),
480
481 // The data and the fit from:
482 // ICRU Report 49, 1993. Ziegler's type of p
483 // Proton kinetic energy for parametrisation
484
485 G4double T = kineticEnergy/(keV*protonMassAM
486
487 static const G4float a[92][5] = {
488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f
489 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f
490 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f
491 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f
492 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f
493 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f
494 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f
495 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f
496 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f
497 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f
498 // Z= 11-20
499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f
500 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f
501 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f
502 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f
503 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f
504 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f
505 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f
506 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f
507 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f
508 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f
509 // Z= 21-30
510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f
511 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f
512 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f
513 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f
514 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f
515 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f
516 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f
517 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f
518 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f
519 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f
520 // Z= 31-40
521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f
522 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f
523 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f
524 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f
525 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f
526 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f
527 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f
528 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f
529 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f
530 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f
531 // Z= 41-50
532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f
533 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f
534 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f
535 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f
536 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f
537 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f
538 // {5.623f, 6.354f, 7160.0f, 337.6f
539 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f
540 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f
541 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f
542 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f
543 // Z= 51-60
544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f
545 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f
546 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f
547 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f
548 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f
549 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f
550 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f
551 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f
552 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f
553 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f
554 // Z= 61-70
555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f
556 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f
557 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f
558 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f
559 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f
560 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f
561 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f
562 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f
563 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f
564 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f
565 // Z= 71-80
566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f
567 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f
568 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f
569 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f
570 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f
571 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f
572 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f
573 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f
574 // {4.856f, 5.460f, 18320.0f, 438.5
575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f
576 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f
577 // Z= 81-90
578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f
579 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f
580 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f
581 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f
582 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f
583 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f
584 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f
585 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f
586 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f
587 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f
588 // Z= 91-92
589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f
590 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f
591 };
592
593 G4double fac = 1.0 ;
594
595 // Carbon specific case for E < 40 keV
596 if ( T < 40.0 && 5 == i) {
597 fac = std::sqrt(T*0.025);
598 T = 40.0;
599
600 // Free electron gas model
601 } else if ( T < 10.0 ) {
602 fac = std::sqrt(T*0.1) ;
603 T = 10.0;
604 }
605
606 // Main parametrisation
607 G4double x1 = (G4double)(a[i][1]);
608 G4double x2 = (G4double)(a[i][2]);
609 G4double x3 = (G4double)(a[i][3]);
610 G4double x4 = (G4double)(a[i][4]);
611 G4double slow = x1 * G4Exp(G4Log(T) * 0.45)
612 G4double shigh = G4Log( 1.0 + x3/T + x4*T )
613 ionloss = slow*shigh*fac / (slow + shigh);
614
615 ionloss = std::max(ionloss, 0.0);
616
617 return ionloss;
618 }
619
620 //....oooOO0OOooo........oooOO0OOooo........oo
621
622 G4double G4BraggModel::DEDX(const G4Material*
623 {
624 G4double eloss = 0.0;
625
626 // check DB
627 if(material != currentMaterial) {
628 currentMaterial = material;
629 baseMaterial = material->GetBaseMaterial()
630 ? material->GetBaseMaterial() : material
631 iPSTAR = -1;
632 iMolecula = -1;
633 iICRU90 = fICRU90 ? fICRU90->GetIndex(base
634
635 if(iICRU90 < 0) {
636 iPSTAR = fPSTAR->GetIndex(baseMaterial);
637 if(iPSTAR < 0) { HasMaterial(baseMateria
638 }
639 //G4cout << "%%% " <<material->GetName() <
640 // << iMolecula << " iPSTAR= " << i
641 // << " iICRU90= " << iICRU90<< G4e
642 }
643
644 // ICRU90 parameterisation
645 if(iICRU90 >= 0) {
646 return fICRU90->GetElectronicDEDXforProton
647 *material->GetDensity();
648 }
649 // PSTAR parameterisation
650 if( iPSTAR >= 0 ) {
651 return fPSTAR->GetElectronicDEDX(iPSTAR, k
652 *material->GetDensity();
653
654 }
655 const std::size_t numberOfElements = materia
656 const G4double* theAtomicNumDensityVector =
657 material->Get
658
659
660 if(iMolecula >= 0) {
661 eloss = StoppingPower(baseMaterial, kineti
662 material->GetDensity
663
664 // Pure material ICRU49 paralmeterisation
665 } else if(1 == numberOfElements) {
666
667 G4double z = material->GetZ();
668 eloss = ElectronicStoppingPower(z, kinetic
669 * (material->Ge
670
671
672 // Experimental data exist only for kinetic
673 } else if( MolecIsInZiegler1988(material) )
674
675 // Loop over elements - calculation based
676 G4double eloss125 = 0.0 ;
677 const G4ElementVector* theElementVector =
678 material->GetElemen
679
680 // Loop for the elements in the material
681 for (std::size_t i=0; i<numberOfElements;
682 const G4Element* element = (*theElementV
683 G4double z = element->GetZ() ;
684 eloss += ElectronicStoppingPower(z,ki
685 * theAtomi
686 eloss125 += ElectronicStoppingPower(z,12
687 * theAtomi
688 }
689
690 // Chemical factor is taken into account
691 if (eloss125 > 0.0) {
692 eloss *= ChemicalFactor(kineticEnergy, e
693 }
694
695 // Brugg's rule calculation
696 } else {
697 const G4ElementVector* theElementVector =
698 material->GetElemen
699
700 // loop for the elements in the material
701 for (std::size_t i=0; i<numberOfElements;
702 {
703 const G4Element* element = (*theElementV
704 eloss += ElectronicStoppingPower(eleme
705 * theAtomic
706 }
707 }
708 return eloss*theZieglerFactor;
709 }
710
711 //....oooOO0OOooo........oooOO0OOooo........oo
712
713 G4bool G4BraggModel::MolecIsInZiegler1988(cons
714 {
715 // The list of molecules from
716 // J.F.Ziegler and J.M.Manoyan, The stopping
717 // Nucl. Inst. & Meth. in Phys. Res. B35 (19
718
719 G4String myFormula = G4String(" ") ;
720 const G4String chFormula = material->GetChem
721 if (myFormula == chFormula ) { return false;
722
723 // There are no evidence for difference of
724 // phase of the compound except for water.
725 // water in gas phase can be predicted usin
726 //
727 // No chemical factor for water-gas
728
729 myFormula = G4String("H_2O") ;
730 const G4State theState = material->GetState(
731 if( theState == kStateGas && myFormula == ch
732
733
734 // The coffecient from Table.4 of Ziegler &
735 static const G4float HeEff = 2.8735f;
736
737 static const std::size_t numberOfMolecula =
738 static const G4String nameOfMol[numberOfMole
739 "H_2O", "C_2H_4O", "C_3H_6O", "C_
740 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH
741 "C_6H_6", "C_4H_10", "C_4H_6", "C_
742 "CF_4", "C_6H_8", "C_6H_12", "C_
743 "C_8H_16", "C_5H_10", "C_5H_8", "C_
744 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_
745 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_
746 "SH_2", "CH_4", "CCLF_3", "CC
747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_
748 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_
749 "C_3H_6S", "C_4H_4S", "C_7H_8"
750 };
751
752 static const G4float expStopping[numberOfMol
753 66.1f, 190.4f, 258.7f, 42.2f, 141.5f,
754 210.9f, 279.6f, 198.8f, 31.0f, 267.5f,
755 122.8f, 311.4f, 260.3f, 328.9f, 391.3f,
756 206.6f, 374.0f, 422.0f, 432.0f, 398.0f,
757 554.0f, 353.0f, 326.0f, 74.6f, 220.5f,
758 197.4f, 362.0f, 170.0f, 330.5f, 211.3f,
759 262.3f, 349.6f, 51.3f, 187.0f, 236.9f,
760 121.9f, 35.8f, 247.0f, 292.6f, 268.0f,
761 262.3f, 49.0f, 398.9f, 444.0f, 22.91f,
762 68.0f, 155.0f, 84.0f, 74.2f, 254.7f,
763 306.8f, 324.4f, 420.0f
764 } ;
765
766 static const G4float expCharge[numberOfMolec
767 HeEff, HeEff, HeEff, 1.0f, HeEff,
768 HeEff, HeEff, HeEff, 1.0f, 1.0f,
769 1.0f, HeEff, HeEff, HeEff, HeEff,
770 HeEff, HeEff, HeEff, HeEff, HeEff,
771 HeEff, HeEff, HeEff, 1.0f, HeEff,
772 HeEff, HeEff, HeEff, HeEff, HeEff,
773 HeEff, HeEff, 1.0f, HeEff, HeEff,
774 HeEff, 1.0f, HeEff, HeEff, HeEff,
775 HeEff, 1.0f, HeEff, HeEff, 1.0f,
776 1.0f, 1.0f, 1.0f, 1.0f, HeEff,
777 HeEff, HeEff, HeEff
778 } ;
779
780 static const G4int numberOfAtomsPerMolecula[
781 3, 7, 10, 4, 6,
782 9, 12, 7, 4, 24,
783 12,14, 10, 13, 5,
784 5, 14, 18, 17, 17,
785 24,15, 13, 9, 8,
786 6, 14, 8, 8, 9,
787 10,15, 6, 7, 7,
788 3, 5, 5, 5, 5,
789 9, 3, 16, 14, 3,
790 9, 16, 11, 9, 10,
791 10, 9, 15};
792
793 // Search for the compaund in the table
794 for (std::size_t i=0; i<numberOfMolecula; ++
795 if(chFormula == nameOfMol[i]) {
796 expStopPower125 = ((G4double)expStopping
797 * (material->GetTotNbOfAtomsPerVolume(
798 ((G4double)(expCharge[i] * numberOfAto
799 return true;
800 }
801 }
802 return false;
803 }
804
805 //....oooOO0OOooo........oooOO0OOooo........oo
806
807 G4double G4BraggModel::ChemicalFactor(G4double
808 G4double
809 {
810 // Approximation of Chemical Factor accordin
811 // J.F.Ziegler and J.M.Manoyan, The stopping
812 // Nucl. Inst. & Meth. in Phys. Res. B35 (19
813
814 static const G4double gamma25 = 1.0 + 25.0*
815 static const G4double gamma125 = 1.0 + 125.0
816 static const G4double beta25 = std::sqrt(1
817 static const G4double beta125 = std::sqrt(1
818 static const G4double f12525 = 1.0 + G4Exp
819
820 G4double gamma = 1.0 + kineticEnergy/proton_
821 G4double beta = std::sqrt(1.0 - 1.0/(gamma*
822
823 G4double factor = 1.0 + (expStopPower125/elo
824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.
825
826 return factor ;
827 }
828
829 //....oooOO0OOooo........oooOO0OOooo........oo
830
831