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