Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4ICRU49NuclearStoppingModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 "G4AutoLock.hh"
 61 #include "G4Pow.hh"
 62 
 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64 
 65 G4double G4ICRU49NuclearStoppingModel::Z23[] = {0.0};
 66 
 67 namespace
 68 {
 69   G4Mutex ICRU49NuclearMutex = G4MUTEX_INITIALIZER;
 70 }
 71 
 72 G4ICRU49NuclearStoppingModel::G4ICRU49NuclearStoppingModel(const G4String& nam) 
 73   : G4VEmModel(nam)
 74 {
 75   theZieglerFactor = eV*cm2*1.0e-15;
 76   g4calc = G4Pow::GetInstance();
 77   InitialiseArray();
 78 }
 79 
 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 81 
 82 G4ICRU49NuclearStoppingModel::~G4ICRU49NuclearStoppingModel() = default;
 83 
 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 85 
 86 void G4ICRU49NuclearStoppingModel::Initialise(const G4ParticleDefinition*, 
 87                 const G4DataVector&)
 88 {}
 89 
 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 91 
 92 void G4ICRU49NuclearStoppingModel::InitialiseArray()
 93 {
 94   if(0.0 != Z23[1]) { return; }
 95   G4AutoLock l(&ICRU49NuclearMutex);
 96   if(0.0 == Z23[1]) {
 97     for(G4int i=2; i<100; ++i) { 
 98       Z23[i] = g4calc->powZ(i, 0.23);
 99     }
100     Z23[1] = 1.0;
101   }
102   l.unlock();
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 
107 void 
108 G4ICRU49NuclearStoppingModel::SampleSecondaries(
109                          std::vector<G4DynamicParticle*>*,
110        const G4MaterialCutsCouple*,
111        const G4DynamicParticle*,
112        G4double, G4double)
113 {}
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 G4double 
118 G4ICRU49NuclearStoppingModel::ComputeDEDXPerVolume(
119                          const G4Material* mat,
120        const G4ParticleDefinition* p,
121        G4double kinEnergy,
122        G4double)
123 {
124   G4double nloss = 0.0;
125   if(kinEnergy <= 0.0) { return nloss; }
126 
127   // projectile
128   G4double mass1 = p->GetPDGMass();
129   G4double z1 = std::abs(p->GetPDGCharge()/eplus);
130 
131   if(kinEnergy*proton_mass_c2/mass1 > z1*z1*MeV) { return nloss; }
132 
133   // Projectile nucleus
134   mass1 /= amu_c2;
135 
136   //  loop for the elements in the material
137   std::size_t numberOfElements = mat->GetNumberOfElements();
138   const G4ElementVector* theElementVector = mat->GetElementVector();
139   const G4double* atomDensity  = mat->GetAtomicNumDensityVector();
140  
141   for (std::size_t iel=0; iel<numberOfElements; ++iel) {
142     const G4Element* element = (*theElementVector)[iel] ;
143     G4double z2 = element->GetZ();
144     G4double mass2 = element->GetN();
145     nloss += (NuclearStoppingPower(kinEnergy, z1, z2, mass1, mass2))
146            * atomDensity[iel];
147   }
148   nloss *= theZieglerFactor;
149   //G4cout << "    nloss= " << nloss << G4endl;
150   return nloss;
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
155 G4double 
156 G4ICRU49NuclearStoppingModel::NuclearStoppingPower(G4double kineticEnergy,
157                G4double z1, G4double z2,
158                G4double mass1, G4double mass2)
159 {
160   G4double energy = kineticEnergy/keV ;  // energy in keV
161   G4double nloss = 0.0;
162   G4double z12 = z1*z2;
163   G4int iz1 = std::min(99, G4lrint(z1));
164   G4int iz2 = std::min(99, G4lrint(z2));
165   
166   G4double rm;
167   if(z1 > 1.5) { 
168     rm = (mass1 + mass2)*(Z23[iz1] + Z23[iz2]); 
169   } else { 
170     rm = (mass1 + mass2)*g4calc->Z13(G4lrint(z2)); 
171   }
172   G4double er = 32.536 * mass2 * energy / ( z12 * rm ) ;  // reduced energy
173   /*
174   G4cout << "   z1= " << iz1 << " z2= " << iz2 << " mass1= " << mass1
175    << " mass2= " << mass2 << " er= " << er << G4endl;
176   */
177   static const G4double nuca[104][2] = {
178   { 1.0E+8, 5.831E-8},
179   { 8.0E+7, 7.288E-8},
180   { 6.0E+7, 9.719E-8},
181   { 5.0E+7, 1.166E-7},
182   { 4.0E+7, 1.457E-7},
183   { 3.0E+7, 1.942E-7},
184   { 2.0E+7, 2.916E-7},
185   { 1.5E+7, 3.887E-7},
186 
187   { 1.0E+7, 5.833E-7},
188   { 8.0E+6, 7.287E-7},
189   { 6.0E+6, 9.712E-7},
190   { 5.0E+6, 1.166E-6},
191   { 4.0E+6, 1.457E-6},
192   { 3.0E+6, 1.941E-6},
193   { 2.0E+6, 2.911E-6},
194   { 1.5E+6, 3.878E-6},
195 
196   { 1.0E+6, 5.810E-6},
197   { 8.0E+5, 7.262E-6},
198   { 6.0E+5, 9.663E-6},
199   { 5.0E+5, 1.157E-5},
200   { 4.0E+5, 1.442E-5},
201   { 3.0E+5, 1.913E-5},
202   { 2.0E+5, 2.845E-5},
203   { 1.5E+5, 3.762E-5},
204 
205   { 1.0E+5, 5.554E-5},
206   { 8.0E+4, 6.866E-5},
207   { 6.0E+4, 9.020E-5},
208   { 5.0E+4, 1.070E-4},
209   { 4.0E+4, 1.319E-4},
210   { 3.0E+4, 1.722E-4},
211   { 2.0E+4, 2.499E-4},
212   { 1.5E+4, 3.248E-4},
213 
214   { 1.0E+4, 4.688E-4},
215   { 8.0E+3, 5.729E-4},
216   { 6.0E+3, 7.411E-4},
217   { 5.0E+3, 8.718E-4},
218   { 4.0E+3, 1.063E-3},
219   { 3.0E+3, 1.370E-3},
220   { 2.0E+3, 1.955E-3},
221   { 1.5E+3, 2.511E-3},
222 
223   { 1.0E+3, 3.563E-3},
224   { 8.0E+2, 4.314E-3},
225   { 6.0E+2, 5.511E-3},
226   { 5.0E+2, 6.430E-3},
227   { 4.0E+2, 7.756E-3},
228   { 3.0E+2, 9.855E-3},
229   { 2.0E+2, 1.375E-2},
230   { 1.5E+2, 1.736E-2},
231 
232   { 1.0E+2, 2.395E-2},
233   { 8.0E+1, 2.850E-2},
234   { 6.0E+1, 3.552E-2},
235   { 5.0E+1, 4.073E-2},
236   { 4.0E+1, 4.802E-2},
237   { 3.0E+1, 5.904E-2},
238   { 1.5E+1, 9.426E-2},
239 
240   { 1.0E+1, 1.210E-1},
241   { 8.0E+0, 1.377E-1},
242   { 6.0E+0, 1.611E-1},
243   { 5.0E+0, 1.768E-1},
244   { 4.0E+0, 1.968E-1},
245   { 3.0E+0, 2.235E-1},
246   { 2.0E+0, 2.613E-1},
247   { 1.5E+0, 2.871E-1},
248 
249   { 1.0E+0, 3.199E-1},
250   { 8.0E-1, 3.354E-1},
251   { 6.0E-1, 3.523E-1},
252   { 5.0E-1, 3.609E-1},
253   { 4.0E-1, 3.693E-1},
254   { 3.0E-1, 3.766E-1},
255   { 2.0E-1, 3.803E-1},
256   { 1.5E-1, 3.788E-1},
257 
258   { 1.0E-1, 3.711E-1},
259   { 8.0E-2, 3.644E-1},
260   { 6.0E-2, 3.530E-1},
261   { 5.0E-2, 3.444E-1},
262   { 4.0E-2, 3.323E-1},
263   { 3.0E-2, 3.144E-1},
264   { 2.0E-2, 2.854E-1},
265   { 1.5E-2, 2.629E-1},
266 
267   { 1.0E-2, 2.298E-1},
268   { 8.0E-3, 2.115E-1},
269   { 6.0E-3, 1.883E-1},
270   { 5.0E-3, 1.741E-1},
271   { 4.0E-3, 1.574E-1},
272   { 3.0E-3, 1.372E-1},
273   { 2.0E-3, 1.116E-1},
274   { 1.5E-3, 9.559E-2},
275 
276   { 1.0E-3, 7.601E-2},
277   { 8.0E-4, 6.668E-2},
278   { 6.0E-4, 5.605E-2},
279   { 5.0E-4, 5.008E-2},
280   { 4.0E-4, 4.352E-2},
281   { 3.0E-4, 3.617E-2},
282   { 2.0E-4, 2.768E-2},
283   { 1.5E-4, 2.279E-2},
284 
285   { 1.0E-4, 1.723E-2},
286   { 8.0E-5, 1.473E-2},
287   { 6.0E-5, 1.200E-2},
288   { 5.0E-5, 1.052E-2},
289   { 4.0E-5, 8.950E-3},
290   { 3.0E-5, 7.246E-3},
291   { 2.0E-5, 5.358E-3},
292   { 1.5E-5, 4.313E-3},
293   { 0.0, 3.166E-3}
294   };
295 
296   if (er >= nuca[0][0]) { nloss = nuca[0][1]; }
297   else {
298     // the table is inverse in energy
299     for (G4int i=102; i>=0; --i) {
300       G4double edi = nuca[i][0];
301       if (er <= edi) {
302         G4double edi1 = nuca[i+1][0];
303         G4double ai   = nuca[i][1];
304         G4double ai1  = nuca[i+1][1];
305         nloss = (ai - ai1)*(er - edi1)/(edi - edi1) + ai1;
306         break;
307       }
308     }
309   }
310 
311   // Stragling
312   if(lossFlucFlag) {
313     G4double sig = 4.0 * mass1 * mass2 / ((mass1 + mass2)*(mass1 + mass2)*
314             (4.0 + 0.197/(er*er) + 6.584/er));
315 
316     nloss *= G4RandGauss::shoot(1.0,sig);
317   }
318    
319   nloss *= 8.462 * z12 * mass1 / rm; // Return to [ev/(10^15 atoms/cm^2]
320 
321   nloss = std::max(nloss, 0.0);
322   return nloss;
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
326