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