| Geant4 Cross Reference |
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 //
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 constructor (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 27-01-03 Make models region aware (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko)
44 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
45 // 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI)
46 // 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
47 // 16-06-05 Fix problem of chemical formula (V.Ivantchenko)
48 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
49 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
50 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
51 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
52
53 // Class Description:
54 //
55 // Implementation of energy loss and delta-electron production by
56 // slow charged heavy particles
57
58 // -------------------------------------------------------------------
59 //
60
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
82
83 G4ICRU90StoppingData* G4BraggModel::fICRU90 = nullptr;
84 G4PSTARStopping* G4BraggModel::fPSTAR = nullptr;
85
86 namespace
87 {
88 G4Mutex ionMutex = G4MUTEX_INITIALIZER;
89 }
90
91 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
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-15;
98 theElectron = G4Electron::Electron();
99 expStopPower125 = 0.0;
100
101 corr = G4LossTableManager::Instance()->EmCorrections();
102 if(nullptr != p) { SetParticle(p); }
103 else { SetParticle(theElectron); }
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108 G4BraggModel::~G4BraggModel()
109 {
110 if(isFirst) {
111 delete fPSTAR;
112 fPSTAR = nullptr;
113 }
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117
118 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
119 const G4DataVector&)
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()->UseICRU90Data()) {
133 fICRU90 = G4NistManager::Instance()->GetICRU90StoppingData();
134 }
135 }
136 l.unlock();
137 }
138 if(isFirst) {
139 if(nullptr != fICRU90) { fICRU90->Initialise(); }
140 fPSTAR->Initialise();
141 }
142
143 if(nullptr == fParticleChange) {
144
145 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
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........oooOO0OOooo........oooOO0OOooo......
159
160 void G4BraggModel::SetParticle(const G4ParticleDefinition* p)
161 {
162 particle = p;
163 mass = particle->GetPDGMass();
164 spin = particle->GetPDGSpin();
165 G4double q = particle->GetPDGCharge()/CLHEP::eplus;
166 chargeSquare = q*q;
167 massRate = mass/CLHEP::proton_mass_c2;
168 ratio = CLHEP::electron_mass_c2/mass;
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172
173 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
174 const G4Material* mat,
175 G4double kinEnergy)
176 {
177 // this method is called only for ions
178 chargeSquare = corr->EffectiveChargeSquareRatio(p, mat, kinEnergy);
179 return chargeSquare;
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183
184 G4double G4BraggModel::MinEnergyCut(const G4ParticleDefinition*,
185 const G4MaterialCutsCouple* couple)
186 {
187 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
193 const G4Material* mat,
194 G4double kineticEnergy)
195 {
196 // this method is called only for ions, so no check if it is an ion
197 return corr->GetParticleCharge(p,mat,kineticEnergy);
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201
202 G4double G4BraggModel::ComputeCrossSectionPerElectron(
203 const G4ParticleDefinition* p,
204 G4double kineticEnergy,
205 G4double cut,
206 G4double maxKinEnergy)
207 {
208 G4double cross = 0.0;
209 const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
210 const G4double maxEnergy = std::min(tmax, maxKinEnergy);
211 const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate);
212 if(cutEnergy < maxEnergy) {
213
214 const G4double energy = kineticEnergy + mass;
215 const G4double energy2 = energy*energy;
216 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
217 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
218 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
219
220 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
221
222 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
223 }
224 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
225 // << " tmax= " << tmax << " cross= " << cross << G4endl;
226 return cross;
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230
231 G4double
232 G4BraggModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
233 G4double kineticEnergy,
234 G4double Z, G4double,
235 G4double cutEnergy,
236 G4double maxEnergy)
237 {
238 return
239 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
240 }
241
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243
244 G4double G4BraggModel::CrossSectionPerVolume(const G4Material* material,
245 const G4ParticleDefinition* p,
246 G4double kineticEnergy,
247 G4double cutEnergy,
248 G4double maxEnergy)
249 {
250 return material->GetElectronDensity()
251 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
252 }
253
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
257 const G4ParticleDefinition* p,
258 G4double kinEnergy,
259 G4double cut)
260 {
261 const G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
262 const G4double tkin = kinEnergy/massRate;
263 const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate);
264 G4double dedx = 0.0;
265
266 // tkin is the scaled energy to proton
267 if (tkin < lowestKinEnergy) {
268 dedx = DEDX(material, lowestKinEnergy)*std::sqrt(tkin/lowestKinEnergy);
269 } else {
270 dedx = DEDX(material, tkin);
271
272 // correction for low cut value using Bethe-Bloch main term
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.)/(tau * (tau + 2.0)) + 1.0 - x) *
278 CLHEP::twopi_mc2_rcl2 * material->GetElectronDensity();
279 }
280 }
281 dedx = std::max(dedx, 0.0) * chargeSquare;
282
283 //G4cout << "E(MeV)= " << tkin/MeV << " dedx= " << dedx
284 // << " " << material->GetName() << G4endl;
285 return dedx;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290 void G4BraggModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
291 const G4MaterialCutsCouple* couple,
292 const G4DynamicParticle* dp,
293 G4double minEnergy,
294 G4double maxEnergy)
295 {
296 const G4double tmax = MaxSecondaryKinEnergy(dp);
297 const G4double xmax = std::min(tmax, maxEnergy);
298 const G4double xmin = std::max(lowestKinEnergy*massRate, std::min(minEnergy, xmax));
299 if(xmin >= xmax) { return; }
300
301 G4double kineticEnergy = dp->GetKineticEnergy();
302 const G4double energy = kineticEnergy + mass;
303 const G4double energy2 = energy*energy;
304 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
305 const G4double grej = 1.0;
306 G4double deltaKinEnergy, f;
307
308 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
309 G4double rndm[2];
310
311 // sampling follows ...
312 do {
313 rndmEngineMod->flatArray(2, rndm);
314 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
315
316 f = 1.0 - beta2*deltaKinEnergy/tmax;
317
318 if(f > grej) {
319 G4cout << "G4BraggModel::SampleSecondary Warning! "
320 << "Majorant " << grej << " < "
321 << f << " for e= " << deltaKinEnergy
322 << G4endl;
323 }
324
325 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
326 } while( grej*rndm[1] >= f );
327
328 G4ThreeVector deltaDirection;
329
330 if(UseAngularGeneratorFlag()) {
331 const G4Material* mat = couple->GetMaterial();
332 G4int Z = SelectRandomAtomNumber(mat);
333
334 deltaDirection =
335 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
336
337 } else {
338
339 G4double deltaMomentum =
340 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
341 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
342 (deltaMomentum * dp->GetTotalMomentum());
343 if(cost > 1.0) { cost = 1.0; }
344 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
345
346 G4double phi = twopi*rndmEngineMod->flat();
347
348 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
349 deltaDirection.rotateUz(dp->GetMomentumDirection());
350 }
351
352 // create G4DynamicParticle object for delta ray
353 auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
354
355 // Change kinematics of primary particle
356 kineticEnergy -= deltaKinEnergy;
357 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
358 finalP = finalP.unit();
359
360 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
361 fParticleChange->SetProposedMomentumDirection(finalP);
362
363 vdp->push_back(delta);
364 }
365
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367
368 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
369 G4double kinEnergy)
370 {
371 if(pd != particle) { SetParticle(pd); }
372 G4double tau = kinEnergy/mass;
373 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
374 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
375 return tmax;
376 }
377
378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379
380 void G4BraggModel::HasMaterial(const G4Material* mat)
381 {
382 const G4String& chFormula = mat->GetChemicalFormula();
383 if(chFormula.empty()) { return; }
384
385 // ICRU Report N49, 1993. Power's model for H
386 static const G4int numberOfMolecula = 11;
387 static const G4String molName[numberOfMolecula] = {
388 "Al_2O_3", "CO_2", "CH_4",
389 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
390 "C_3H_8", "SiO_2", "H_2O",
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........oooOO0OOooo........oooOO0OOooo......
404
405 G4double G4BraggModel::StoppingPower(const G4Material* material,
406 G4double kineticEnergy)
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 for protons.
414 // Proton kinetic energy for parametrisation (keV/amu)
415
416 G4double T = kineticEnergy/(keV*protonMassAMU) ;
417
418 static const G4float a[11][5] = {
419 {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f, 2.153E-2f},
420 {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f, 7.966E-3f},
421 {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f, 8.153E-3f},
422 {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f, 9.383E-3f},
423 {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f, 3.512E-2f},
424 {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f, 1.273E-1f},
425 {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f, 3.775E-2f},
426 {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f, 1.007E-2f},
427 {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f, 7.904E-3f},
428 {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f, 8.572E-3f},
429 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f} };
430
431 static const G4float atomicWeight[11] = {
432 101.96128f, 44.0098f, 16.0426f, 28.0536f, 42.0804f,
433 104.1512f, 44.665f, 60.0843f, 18.0152f, 18.0152f, 12.0f};
434
435 if ( T < 10.0 ) {
436 ionloss = ((G4double)(a[iMolecula][0])) * std::sqrt(T) ;
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.45);
444 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
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/G4Log(10.);
451
452 if (T < 100.0) {
453 ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10);
454 }
455 else if (T < 700.0) {
456 ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10);
457 }
458 else if (T < 10000.0) {
459 ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10);
460 }
461 }
462 ionloss /= (G4double)atomicWeight[iMolecula];
463
464 // pure material (normally not the case for this function)
465 } else if(1 == (material->GetNumberOfElements())) {
466 G4double z = material->GetZ() ;
467 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
468 }
469
470 return ionloss;
471 }
472
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474
475 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
476 G4double kineticEnergy) const
477 {
478 G4double ionloss ;
479 G4int i = std::min(std::max(G4lrint(z)-1,0),91); // index of atom
480
481 // The data and the fit from:
482 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
483 // Proton kinetic energy for parametrisation (keV/amu)
484
485 G4double T = kineticEnergy/(keV*protonMassAMU) ;
486
487 static const G4float a[92][5] = {
488 {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f, 1.159E-1f},
489 {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f, 5.225E-2f},
490 {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f, 4.578E-2f},
491 {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f, 3.475E-2f},
492 {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f, 2.855E-2f},
493 {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f, 1.638E-2f},
494 {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f, 2.513E-2f},
495 {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f, 2.230E-2f},
496 {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f, 1.816E-2f},
497 {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f, 1.568E-2f},
498 // Z= 11-20
499 {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f, 1.472E-2f},
500 {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f, 1.397E-2f},
501 {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f, 2.023E-2f},
502 {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f, 1.419E-2f},
503 {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f, 1.267E-2f},
504 {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f, 1.211E-2f},
505 {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f, 1.183E-2f},
506 {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f, 1.123E-2f},
507 {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f, 1.129E-2f},
508 {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f, 1.112E-2f},
509 // Z= 21-30
510 {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f, 9.995E-3f},
511 {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f, 8.930E-3f},
512 {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f, 9.117E-3f},
513 {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f, 8.413E-3f},
514 {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f, 8.829E-3f},
515 {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f, 7.782E-3f},
516 {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f, 7.361E-3f},
517 {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f, 8.763E-3f},
518 {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f, 2.242E-2f},
519 {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f, 6.809E-3f},
520 // Z= 31-40
521 {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f, 6.725E-3f},
522 {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f, 9.689E-3f},
523 {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f, 6.447E-3f},
524 {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f, 7.684E-3f},
525 {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f, 6.509E-3f},
526 {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f, 2.880E-3f},
527 {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f, 6.087E-3f},
528 {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f, 6.003E-3f},
529 {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f, 5.889E-3f},
530 {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f, 5.765E-3f},
531 // Z= 41-50
532 {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f, 5.587E-3f},
533 {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f, 5.376E-3f},
534 {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f, 5.315E-3f},
535 {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f, 5.151E-3f},
536 {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f, 4.919E-3f},
537 {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f, 4.758E-3f},
538 // {5.623f, 6.354f, 7160.0f, 337.6f, 0.013940f}, // Ag Ziegler77
539 {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f, 1.676E-2f}, // Ag ICRU49
540 {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f, 4.626E-3f},
541 {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f, 4.540E-3f},
542 {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f, 4.474E-3f},
543 // Z= 51-60
544 {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f, 9.074E-3f},
545 {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f, 4.402E-3f},
546 {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f, 4.376E-3f},
547 {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f, 6.206E-3f},
548 {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f, 4.447E-3f},
549 {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f, 4.511E-3f},
550 {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f, 4.540E-3f},
551 {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f, 4.420E-3f},
552 {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f, 4.298E-3f},
553 {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f, 4.182E-3f},
554 // Z= 61-70
555 {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f, 4.058E-3f},
556 {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f, 3.976E-3f},
557 {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f, 3.877E-3f},
558 {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f, 3.863E-3f},
559 {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f, 3.725E-3f},
560 {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f, 3.632E-3f},
561 {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f, 3.498E-3f},
562 {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f, 3.405E-3f},
563 {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f, 3.342E-3f},
564 {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f, 3.292E-3f},
565 // Z= 71-80
566 {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f, 3.243E-3f},
567 {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f, 3.195E-3f},
568 {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f, 3.186E-3f},
569 {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f, 3.406E-3f},
570 {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f, 3.122E-3f},
571 {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f, 3.082E-3f},
572 {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f, 2.965E-3f},
573 {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f, 2.871E-3f},
574 // {4.856f, 5.460f, 18320.0f, 438.5f, 0.002542f}, //Ziegler77
575 {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f, 2.077E-2f}, //ICRU49
576 {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f, 2.882E-3f},
577 // Z= 81-90
578 {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f, 2.913E-3f},
579 {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f, 2.871E-3f},
580 {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f, 2.660E-3f},
581 {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f, 2.812E-3f},
582 {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f, 2.776E-3f},
583 {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f, 2.748E-3f},
584 {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f, 2.737E-3f},
585 {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f, 2.727E-3f},
586 {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f, 2.697E-3f},
587 {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f, 2.641E-3f},
588 // Z= 91-92
589 {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f, 2.603E-3f},
590 {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f, 2.673E-3f}
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 ) * x2/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........oooOO0OOooo........oooOO0OOooo......
621
622 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
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(baseMaterial) : -1;
634
635 if(iICRU90 < 0) {
636 iPSTAR = fPSTAR->GetIndex(baseMaterial);
637 if(iPSTAR < 0) { HasMaterial(baseMaterial); }
638 }
639 //G4cout << "%%% " <<material->GetName() << " iMolecula= "
640 // << iMolecula << " iPSTAR= " << iPSTAR
641 // << " iICRU90= " << iICRU90<< G4endl;
642 }
643
644 // ICRU90 parameterisation
645 if(iICRU90 >= 0) {
646 return fICRU90->GetElectronicDEDXforProton(iICRU90, kineticEnergy)
647 *material->GetDensity();
648 }
649 // PSTAR parameterisation
650 if( iPSTAR >= 0 ) {
651 return fPSTAR->GetElectronicDEDX(iPSTAR, kineticEnergy)
652 *material->GetDensity();
653
654 }
655 const std::size_t numberOfElements = material->GetNumberOfElements();
656 const G4double* theAtomicNumDensityVector =
657 material->GetAtomicNumDensityVector();
658
659
660 if(iMolecula >= 0) {
661 eloss = StoppingPower(baseMaterial, kineticEnergy)*
662 material->GetDensity()/amu;
663
664 // Pure material ICRU49 paralmeterisation
665 } else if(1 == numberOfElements) {
666
667 G4double z = material->GetZ();
668 eloss = ElectronicStoppingPower(z, kineticEnergy)
669 * (material->GetTotNbOfAtomsPerVolume());
670
671
672 // Experimental data exist only for kinetic energy 125 keV
673 } else if( MolecIsInZiegler1988(material) ) {
674
675 // Loop over elements - calculation based on Bragg's rule
676 G4double eloss125 = 0.0 ;
677 const G4ElementVector* theElementVector =
678 material->GetElementVector();
679
680 // Loop for the elements in the material
681 for (std::size_t i=0; i<numberOfElements; ++i) {
682 const G4Element* element = (*theElementVector)[i] ;
683 G4double z = element->GetZ() ;
684 eloss += ElectronicStoppingPower(z,kineticEnergy)
685 * theAtomicNumDensityVector[i] ;
686 eloss125 += ElectronicStoppingPower(z,125.0*keV)
687 * theAtomicNumDensityVector[i] ;
688 }
689
690 // Chemical factor is taken into account
691 if (eloss125 > 0.0) {
692 eloss *= ChemicalFactor(kineticEnergy, eloss125);
693 }
694
695 // Brugg's rule calculation
696 } else {
697 const G4ElementVector* theElementVector =
698 material->GetElementVector() ;
699
700 // loop for the elements in the material
701 for (std::size_t i=0; i<numberOfElements; ++i)
702 {
703 const G4Element* element = (*theElementVector)[i] ;
704 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
705 * theAtomicNumDensityVector[i];
706 }
707 }
708 return eloss*theZieglerFactor;
709 }
710
711 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
712
713 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
714 {
715 // The list of molecules from
716 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
717 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
718
719 G4String myFormula = G4String(" ") ;
720 const G4String chFormula = material->GetChemicalFormula() ;
721 if (myFormula == chFormula ) { return false; }
722
723 // There are no evidence for difference of stopping power depended on
724 // phase of the compound except for water. The stopping power of the
725 // water in gas phase can be predicted using Bragg's rule.
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 == chFormula) return false ;
732
733
734 // The coffecient from Table.4 of Ziegler & Manoyan
735 static const G4float HeEff = 2.8735f;
736
737 static const std::size_t numberOfMolecula = 53;
738 static const G4String nameOfMol[numberOfMolecula] = {
739 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
740 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
741 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
742 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
743 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
744 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
745 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
746 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
747 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
748 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
749 "C_3H_6S", "C_4H_4S", "C_7H_8"
750 };
751
752 static const G4float expStopping[numberOfMolecula] = {
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[numberOfMolecula] = {
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[numberOfMolecula] = {
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; ++i) {
795 if(chFormula == nameOfMol[i]) {
796 expStopPower125 = ((G4double)expStopping[i])
797 * (material->GetTotNbOfAtomsPerVolume()) /
798 ((G4double)(expCharge[i] * numberOfAtomsPerMolecula[i]));
799 return true;
800 }
801 }
802 return false;
803 }
804
805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806
807 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
808 G4double eloss125) const
809 {
810 // Approximation of Chemical Factor according to
811 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
812 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
813
814 static const G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2;
815 static const G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2;
816 static const G4double beta25 = std::sqrt(1.0 - 1.0/(gamma25*gamma25));
817 static const G4double beta125 = std::sqrt(1.0 - 1.0/(gamma125*gamma125));
818 static const G4double f12525 = 1.0 + G4Exp( 1.48*(beta125/beta25 - 7.0) );
819
820 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2;
821 G4double beta = std::sqrt(1.0 - 1.0/(gamma*gamma));
822
823 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) * f12525/
824 (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) );
825
826 return factor ;
827 }
828
829 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
830
831