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