Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4ICRU49NuclearStoppingModel.cc 106442 2017-10-10 08:00:30Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4ICRU49NuclearStoppingModel
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 09.04.2008 from G4MuMscModel
38 //
39 // Modifications:
40 //
41 //
42 // Class Description:
43 //
44 // Implementation of the model of ICRU'49 nuclear stopping
45 
46 // -------------------------------------------------------------------
47 //
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "Randomize.hh"
56 #include "G4LossTableManager.hh"
58 #include "G4ElementVector.hh"
59 #include "G4ProductionCutsTable.hh"
60 #include "G4Step.hh"
61 #include "G4Pow.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
66 
67 #ifdef G4MULTITHREADED
68 G4Mutex G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex = G4MUTEX_INITIALIZER;
69 #endif
70 
72  : G4VEmModel(nam)
73 {
74  theZieglerFactor = eV*cm2*1.0e-15;
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
82 {}
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
87  const G4DataVector&)
88 {}
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
93 {
94  if(0.0 == Z23[1]) {
95 #ifdef G4MULTITHREADED
96  G4MUTEXLOCK(&G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex);
97 #endif
98  if(0.0 == Z23[2]) {
99  for(G4int i=1; i<100; ++i) {
100  Z23[i] = g4calc->powZ(i, 0.23);
101  }
102  }
103 #ifdef G4MULTITHREADED
104  G4MUTEXUNLOCK(&G4ICRU49NuclearStoppingModel::ICRU49NuclearMutex);
105 #endif
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 void
113  std::vector<G4DynamicParticle*>*,
114  const G4MaterialCutsCouple*,
115  const G4DynamicParticle*,
117 {}
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 G4double
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
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......
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Z13(G4int Z) const
Definition: G4Pow.hh:132
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double keV
Definition: G4SIunits.hh:216
const char * p
Definition: xmltok.h:285
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
G4double GetPDGCharge() const
G4double GetPDGMass() const
static constexpr double proton_mass_c2
G4double GetN() const
Definition: G4Element.hh:135
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:237
double G4double
Definition: G4Types.hh:76
static constexpr double amu_c2
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
static constexpr double cm2
Definition: G4SIunits.hh:120
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) final
G4bool lossFlucFlag
Definition: G4VEmModel.hh:436
double energy
Definition: plottest35.C:25
static constexpr double eV
Definition: G4SIunits.hh:215
Float_t mat
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double) final
ThreeVector shoot(const G4int Ap, const G4int Af)
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
static constexpr double eplus
Definition: G4SIunits.hh:199
int G4lrint(double ad)
Definition: templates.hh:151
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX) final
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4double GetZ() const
Definition: G4Element.hh:131
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:217
G4ICRU49NuclearStoppingModel(const G4String &nam="ICRU49NucStopping")
G4double NuclearStoppingPower(G4double kineticEnergy, G4double Z1, G4double Z2, G4double A1, G4double A2)
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::mutex G4Mutex
Definition: G4Threading.hh:84