ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ICRU49NuclearStoppingModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ICRU49NuclearStoppingModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4ICRU49NuclearStoppingModel
33 //
34 // Author: V.Ivanchenko
35 //
36 // Creation date: 09.04.2008 from G4MuMscModel
37 //
38 // Modifications:
39 //
40 //
41 // Class Description:
42 //
43 // Implementation of the model of ICRU'49 nuclear stopping
44 
45 // -------------------------------------------------------------------
46 //
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 #include "G4ICRU49NuclearStoppingModel.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "Randomize.hh"
55 #include "G4LossTableManager.hh"
56 #include "G4ParticleChangeForLoss.hh"
57 #include "G4ElementVector.hh"
58 #include "G4ProductionCutsTable.hh"
59 #include "G4Step.hh"
60 #include "G4Pow.hh"
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
64 G4double G4ICRU49NuclearStoppingModel::Z23[] = {0.0};
65 
66 #ifdef G4MULTITHREADED
67 G4Mutex G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex = G4MUTEX_INITIALIZER;
68 #endif
69 
70 G4ICRU49NuclearStoppingModel::G4ICRU49NuclearStoppingModel(const G4String& nam)
71  : G4VEmModel(nam)
72 {
73  theZieglerFactor = eV*cm2*1.0e-15;
74  g4calc = G4Pow::GetInstance();
75  InitialiseArray();
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 
80 G4ICRU49NuclearStoppingModel::~G4ICRU49NuclearStoppingModel()
81 {}
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 void G4ICRU49NuclearStoppingModel::Initialise(const G4ParticleDefinition*,
86  const G4DataVector&)
87 {}
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void G4ICRU49NuclearStoppingModel::InitialiseArray()
92 {
93  if(0.0 == Z23[1]) {
94 #ifdef G4MULTITHREADED
95  G4MUTEXLOCK(&G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex);
96 #endif
97  if(0.0 == Z23[1]) {
98  for(G4int i=2; i<100; ++i) {
99  Z23[i] = g4calc->powZ(i, 0.23);
100  }
101  Z23[1] = 1.0;
102  }
103 #ifdef G4MULTITHREADED
104  G4MUTEXUNLOCK(&G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex);
105 #endif
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 void
112 G4ICRU49NuclearStoppingModel::SampleSecondaries(
113  std::vector<G4DynamicParticle*>*,
114  const G4MaterialCutsCouple*,
115  const G4DynamicParticle*,
116  G4double, G4double)
117 {}
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 G4double
122 G4ICRU49NuclearStoppingModel::ComputeDEDXPerVolume(
123  const G4Material* mat,
124  const G4ParticleDefinition* p,
125  G4double kinEnergy,
126  G4double)
127 {
128  G4double nloss = 0.0;
129  if(kinEnergy <= 0.0) { return nloss; }
130 
131  // projectile
132  G4double mass1 = p->GetPDGMass();
133  G4double z1 = std::abs(p->GetPDGCharge()/eplus);
134 
135  if(kinEnergy*proton_mass_c2/mass1 > z1*z1*MeV) { return nloss; }
136 
137  // Projectile nucleus
138  mass1 /= amu_c2;
139 
140  // loop for the elements in the material
141  G4int numberOfElements = mat->GetNumberOfElements();
142  const G4ElementVector* theElementVector = mat->GetElementVector();
143  const G4double* atomDensity = mat->GetAtomicNumDensityVector();
144 
145  for (G4int iel=0; iel<numberOfElements; ++iel) {
146  const G4Element* element = (*theElementVector)[iel] ;
147  G4double z2 = element->GetZ();
148  G4double mass2 = element->GetN();
149  nloss += (NuclearStoppingPower(kinEnergy, z1, z2, mass1, mass2))
150  * atomDensity[iel];
151  }
152  nloss *= theZieglerFactor;
153  //G4cout << " nloss= " << nloss << G4endl;
154  return nloss;
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
159 G4double
160 G4ICRU49NuclearStoppingModel::NuclearStoppingPower(G4double kineticEnergy,
161  G4double z1, G4double z2,
162  G4double mass1, G4double mass2)
163 {
164  G4double energy = kineticEnergy/keV ; // energy in keV
165  G4double nloss = 0.0;
166  G4double z12 = z1*z2;
167  G4int iz1 = std::min(99, G4lrint(z1));
168  G4int iz2 = std::min(99, G4lrint(z2));
169 
170  G4double rm;
171  if(z1 > 1.5) {
172  rm = (mass1 + mass2)*(Z23[iz1] + Z23[iz2]);
173  } else {
174  rm = (mass1 + mass2)*g4calc->Z13(G4lrint(z2));
175  }
176  G4double er = 32.536 * mass2 * energy / ( z12 * rm ) ; // reduced energy
177  /*
178  G4cout << " z1= " << iz1 << " z2= " << iz2 << " mass1= " << mass1
179  << " mass2= " << mass2 << " er= " << er << G4endl;
180  */
181  static const G4double nuca[104][2] = {
182  { 1.0E+8, 5.831E-8},
183  { 8.0E+7, 7.288E-8},
184  { 6.0E+7, 9.719E-8},
185  { 5.0E+7, 1.166E-7},
186  { 4.0E+7, 1.457E-7},
187  { 3.0E+7, 1.942E-7},
188  { 2.0E+7, 2.916E-7},
189  { 1.5E+7, 3.887E-7},
190 
191  { 1.0E+7, 5.833E-7},
192  { 8.0E+6, 7.287E-7},
193  { 6.0E+6, 9.712E-7},
194  { 5.0E+6, 1.166E-6},
195  { 4.0E+6, 1.457E-6},
196  { 3.0E+6, 1.941E-6},
197  { 2.0E+6, 2.911E-6},
198  { 1.5E+6, 3.878E-6},
199 
200  { 1.0E+6, 5.810E-6},
201  { 8.0E+5, 7.262E-6},
202  { 6.0E+5, 9.663E-6},
203  { 5.0E+5, 1.157E-5},
204  { 4.0E+5, 1.442E-5},
205  { 3.0E+5, 1.913E-5},
206  { 2.0E+5, 2.845E-5},
207  { 1.5E+5, 3.762E-5},
208 
209  { 1.0E+5, 5.554E-5},
210  { 8.0E+4, 6.866E-5},
211  { 6.0E+4, 9.020E-5},
212  { 5.0E+4, 1.070E-4},
213  { 4.0E+4, 1.319E-4},
214  { 3.0E+4, 1.722E-4},
215  { 2.0E+4, 2.499E-4},
216  { 1.5E+4, 3.248E-4},
217 
218  { 1.0E+4, 4.688E-4},
219  { 8.0E+3, 5.729E-4},
220  { 6.0E+3, 7.411E-4},
221  { 5.0E+3, 8.718E-4},
222  { 4.0E+3, 1.063E-3},
223  { 3.0E+3, 1.370E-3},
224  { 2.0E+3, 1.955E-3},
225  { 1.5E+3, 2.511E-3},
226 
227  { 1.0E+3, 3.563E-3},
228  { 8.0E+2, 4.314E-3},
229  { 6.0E+2, 5.511E-3},
230  { 5.0E+2, 6.430E-3},
231  { 4.0E+2, 7.756E-3},
232  { 3.0E+2, 9.855E-3},
233  { 2.0E+2, 1.375E-2},
234  { 1.5E+2, 1.736E-2},
235 
236  { 1.0E+2, 2.395E-2},
237  { 8.0E+1, 2.850E-2},
238  { 6.0E+1, 3.552E-2},
239  { 5.0E+1, 4.073E-2},
240  { 4.0E+1, 4.802E-2},
241  { 3.0E+1, 5.904E-2},
242  { 1.5E+1, 9.426E-2},
243 
244  { 1.0E+1, 1.210E-1},
245  { 8.0E+0, 1.377E-1},
246  { 6.0E+0, 1.611E-1},
247  { 5.0E+0, 1.768E-1},
248  { 4.0E+0, 1.968E-1},
249  { 3.0E+0, 2.235E-1},
250  { 2.0E+0, 2.613E-1},
251  { 1.5E+0, 2.871E-1},
252 
253  { 1.0E+0, 3.199E-1},
254  { 8.0E-1, 3.354E-1},
255  { 6.0E-1, 3.523E-1},
256  { 5.0E-1, 3.609E-1},
257  { 4.0E-1, 3.693E-1},
258  { 3.0E-1, 3.766E-1},
259  { 2.0E-1, 3.803E-1},
260  { 1.5E-1, 3.788E-1},
261 
262  { 1.0E-1, 3.711E-1},
263  { 8.0E-2, 3.644E-1},
264  { 6.0E-2, 3.530E-1},
265  { 5.0E-2, 3.444E-1},
266  { 4.0E-2, 3.323E-1},
267  { 3.0E-2, 3.144E-1},
268  { 2.0E-2, 2.854E-1},
269  { 1.5E-2, 2.629E-1},
270 
271  { 1.0E-2, 2.298E-1},
272  { 8.0E-3, 2.115E-1},
273  { 6.0E-3, 1.883E-1},
274  { 5.0E-3, 1.741E-1},
275  { 4.0E-3, 1.574E-1},
276  { 3.0E-3, 1.372E-1},
277  { 2.0E-3, 1.116E-1},
278  { 1.5E-3, 9.559E-2},
279 
280  { 1.0E-3, 7.601E-2},
281  { 8.0E-4, 6.668E-2},
282  { 6.0E-4, 5.605E-2},
283  { 5.0E-4, 5.008E-2},
284  { 4.0E-4, 4.352E-2},
285  { 3.0E-4, 3.617E-2},
286  { 2.0E-4, 2.768E-2},
287  { 1.5E-4, 2.279E-2},
288 
289  { 1.0E-4, 1.723E-2},
290  { 8.0E-5, 1.473E-2},
291  { 6.0E-5, 1.200E-2},
292  { 5.0E-5, 1.052E-2},
293  { 4.0E-5, 8.950E-3},
294  { 3.0E-5, 7.246E-3},
295  { 2.0E-5, 5.358E-3},
296  { 1.5E-5, 4.313E-3},
297  { 0.0, 3.166E-3}
298  };
299 
300  if (er >= nuca[0][0]) { nloss = nuca[0][1]; }
301  else {
302  // the table is inverse in energy
303  for (G4int i=102; i>=0; --i) {
304  G4double edi = nuca[i][0];
305  if (er <= edi) {
306  G4double edi1 = nuca[i+1][0];
307  G4double ai = nuca[i][1];
308  G4double ai1 = nuca[i+1][1];
309  nloss = (ai - ai1)*(er - edi1)/(edi - edi1) + ai1;
310  break;
311  }
312  }
313  }
314 
315  // Stragling
316  if(lossFlucFlag) {
317  G4double sig = 4.0 * mass1 * mass2 / ((mass1 + mass2)*(mass1 + mass2)*
318  (4.0 + 0.197/(er*er) + 6.584/er));
319 
320  nloss *= G4RandGauss::shoot(1.0,sig);
321  }
322 
323  nloss *= 8.462 * z12 * mass1 / rm; // Return to [ev/(10^15 atoms/cm^2]
324 
325  nloss = std::max(nloss, 0.0);
326  return nloss;
327 }
328 
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......