| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4PAIySection.cc,v 1.4 2009/07/26 15:51:01 vnivanch Exp $
>> 28 // GEANT4 tag $Name: geant4-09-03-patch-02 $
>> 29 //
27 // 30 //
28 // G4PAIySection.cc -- class implementation fi 31 // G4PAIySection.cc -- class implementation file
29 // 32 //
30 // GEANT 4 class implementation file 33 // GEANT 4 class implementation file
31 // 34 //
32 // For information related to this code, pleas 35 // For information related to this code, please, contact
33 // the Geant4 Collaboration. 36 // the Geant4 Collaboration.
34 // 37 //
35 // R&D: Vladimir.Grichine@cern.ch 38 // R&D: Vladimir.Grichine@cern.ch
36 // 39 //
37 // History: 40 // History:
38 // 41 //
39 // 01.10.07, V.Ivanchenko create using V.Grich 42 // 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class
40 // 26.07.09, V.Ivanchenko added protection for 43 // 26.07.09, V.Ivanchenko added protection for mumerical exceptions for
41 // low-density materials 44 // low-density materials
42 // 21.11.10 V. Grichine bug fixed in Initialis <<
43 // material. Warning: the table is <<
44 // 23.06.13 V.Grichine arrays->G4DataVectors <<
45 // 45 //
46 46
47 #include "G4PAIySection.hh" 47 #include "G4PAIySection.hh"
48 48
49 #include "globals.hh" 49 #include "globals.hh"
50 #include "G4PhysicalConstants.hh" <<
51 #include "G4SystemOfUnits.hh" <<
52 #include "G4ios.hh" 50 #include "G4ios.hh"
53 #include "G4Poisson.hh" 51 #include "G4Poisson.hh"
54 #include "G4Material.hh" 52 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh" 53 #include "G4MaterialCutsCouple.hh"
56 #include "G4SandiaTable.hh" 54 #include "G4SandiaTable.hh"
57 #include "G4Exp.hh" <<
58 #include "G4Log.hh" <<
59 55
60 using namespace std; 56 using namespace std;
61 57
62 // Local class constants 58 // Local class constants
63 59
64 const G4double G4PAIySection::fDelta = 0.005; << 60 const G4double G4PAIySection::fDelta = 0.005 ; // energy shift from interval border
65 const G4double G4PAIySection::fError = 0.005; << 61 const G4double G4PAIySection::fError = 0.005 ; // error in lin-log approximation
66 62
67 const G4int G4PAIySection::fMaxSplineSize = 50 << 63 const G4int G4PAIySection::fMaxSplineSize = 500 ; // Max size of output spline
68 << 64 // arrays
69 65
70 ////////////////////////////////////////////// 66 //////////////////////////////////////////////////////////////////
71 // 67 //
72 // Constructor 68 // Constructor
73 // 69 //
74 70
75 G4PAIySection::G4PAIySection() 71 G4PAIySection::G4PAIySection()
76 { << 72 {}
77 fSandia = nullptr; <<
78 fDensity = fElectronDensity = fNormalization <<
79 fIntervalNumber = fSplineNumber = 0; <<
80 fVerbose = 0; <<
81 <<
82 betaBohr = fine_structure_const; <<
83 G4double cofBetaBohr = 4.0; <<
84 G4double betaBohr2 = fine_structure_const*fi <<
85 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; <<
86 <<
87 fSplineEnergy = G4DataVector(fMaxSp <<
88 fRePartDielectricConst = G4DataVector(fMaxSp <<
89 fImPartDielectricConst = G4DataVector(fMaxSp <<
90 fIntegralTerm = G4DataVector(fMaxSp <<
91 fDifPAIySection = G4DataVector(fMaxSp <<
92 fdNdxCerenkov = G4DataVector(fMaxSp <<
93 fdNdxPlasmon = G4DataVector(fMaxSp <<
94 fIntegralPAIySection = G4DataVector(fMaxSp <<
95 fIntegralPAIdEdx = G4DataVector(fMaxSp <<
96 fIntegralCerenkov = G4DataVector(fMaxSp <<
97 fIntegralPlasmon = G4DataVector(fMaxSp <<
98 <<
99 for( G4int i = 0; i < 500; ++i ) <<
100 { <<
101 for( G4int j = 0; j < 112; ++j ) { fPAItab <<
102 } <<
103 } <<
104 73
105 ////////////////////////////////////////////// 74 ////////////////////////////////////////////////////////////////////////////
106 // 75 //
107 // << 76 // Destructor
108 77
109 G4double G4PAIySection::GetLorentzFactor(G4int << 78 G4PAIySection::~G4PAIySection()
110 { << 79 {}
111 return fLorentzFactor[j]; <<
112 } <<
113 80
114 ////////////////////////////////////////////// 81 ////////////////////////////////////////////////////////////////////////
115 // 82 //
116 // Constructor with beta*gamma square value ca << 83 // Test Constructor with beta*gamma square value
117 84
118 void G4PAIySection::Initialize( const G4Materi 85 void G4PAIySection::Initialize( const G4Material* material,
119 G4double maxEn << 86 G4double maxEnergyTransfer,
120 G4double betaG << 87 G4double betaGammaSq)
121 G4SandiaTable* <<
122 { 88 {
123 if(fVerbose > 0) << 89 G4int i, j, numberOfElements ;
124 { << 90
125 G4cout<<G4endl; << 91 fDensity = material->GetDensity();
126 G4cout<<"G4PAIySection::Initialize(...,G4S << 92 fElectronDensity = material->GetElectronDensity() ;
127 G4cout<<G4endl; << 93 numberOfElements = material->GetNumberOfElements() ;
128 } <<
129 G4int i, j; <<
130 94
131 fSandia = sandia; << 95 fSandia = material->GetSandiaTable();
132 fIntervalNumber = sandia->GetMaxInterval(); <<
133 fDensity = material->GetDensity(); <<
134 fElectronDensity = material->GetElectronDens <<
135 96
136 // fIntervalNumber--; << 97 fIntervalNumber = fSandia->GetMaxInterval();
137 98
138 if( fVerbose > 0 ) << 99 fIntervalNumber--;
139 { <<
140 G4cout<<"fDensity = "<<fDensity<<"\t"<<fEl <<
141 <<fIntervalNumber<< " (beta*gamma)^2 <<
142 } <<
143 fEnergyInterval = G4DataVector(fIntervalNumb <<
144 fA1 = G4DataVector(fIntervalNumb <<
145 fA2 = G4DataVector(fIntervalNumb <<
146 fA3 = G4DataVector(fIntervalNumb <<
147 fA4 = G4DataVector(fIntervalNumb <<
148 100
149 for( i = 1; i <= fIntervalNumber; ++i ) << 101 for(i=1;i<=fIntervalNumber;i++)
150 { << 102 {
151 if ( sandia->GetSandiaMatTablePAI(i-1,0) < << 103 G4double e = fSandia->GetSandiaMatTablePAI(i,0);
152 { << 104 if(e >= maxEnergyTransfer || i > fIntervalNumber)
153 fIntervalNumber--; << 105 {
154 continue; << 106 fEnergyInterval[i] = maxEnergyTransfer ;
155 } << 107 fIntervalNumber = i ;
156 if( ( sandia->GetSandiaMatTablePAI(i-1,0) << 108 break;
157 || i >= fIntervalNumber ) << 109 }
158 { << 110 fEnergyInterval[i] = e;
159 fEnergyInterval[i] = maxEnergyTransfer; << 111 fA1[i] = fSandia->GetSandiaMatTablePAI(i,1);
160 fIntervalNumber = i; << 112 fA2[i] = fSandia->GetSandiaMatTablePAI(i,2);
161 break; << 113 fA3[i] = fSandia->GetSandiaMatTablePAI(i,3);
162 } << 114 fA4[i] = fSandia->GetSandiaMatTablePAI(i,4);
163 fEnergyInterval[i] = sandia->GetSandiaMatT << 115
164 fA1[i] = sandia->GetSandiaMatT << 116 }
165 fA2[i] = sandia->GetSandiaMatT << 117 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
166 fA3[i] = sandia->GetSandiaMatT << 118 {
167 fA4[i] = sandia->GetSandiaMatT << 119 fIntervalNumber++;
168 << 120 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
169 if( fVerbose > 0 ) { << 121 fA1[fIntervalNumber] = fA1[fIntervalNumber-1] ;
170 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< << 122 fA2[fIntervalNumber] = fA2[fIntervalNumber-1] ;
171 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4en << 123 fA3[fIntervalNumber] = fA3[fIntervalNumber-1] ;
172 } << 124 fA4[fIntervalNumber] = fA4[fIntervalNumber-1] ;
173 } << 125 }
174 if( fVerbose > 0 ) { << 126
175 G4cout<<"last i = "<<i<<"; "<<"fIntervalNu << 127 // Now checking, if two borders are too close together
176 <<fIntervalNumber<<G4endl; << 128 for(i=1;i<fIntervalNumber;i++)
177 } <<
178 if( fEnergyInterval[fIntervalNumber] != maxE <<
179 { <<
180 fIntervalNumber++; <<
181 fEnergyInterval[fIntervalNumber] = maxEn <<
182 } <<
183 if( fVerbose > 0 ) <<
184 { <<
185 for( i = 1; i <= fIntervalNumber; ++i ) <<
186 { <<
187 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< <<
188 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; <<
189 } <<
190 } <<
191 if( fVerbose > 0 ) { <<
192 G4cout<<"Now checking, if two borders are <<
193 } <<
194 for( i = 1; i < fIntervalNumber; ++i ) <<
195 { <<
196 if( fEnergyInterval[i+1]-fEnergyInterval[i <<
197 1.5*fDelta*(fEnergyInterval[i+1]+fEne <<
198 else <<
199 { <<
200 for( j = i; j < fIntervalNumber; j++ ) <<
201 { 129 {
202 fEnergyInterval[j] = fEnergyInte << 130 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
203 fA1[j] = fA1[j+1]; << 131 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
204 fA2[j] = fA2[j+1]; << 132 if(fEnergyInterval[i+1]-fEnergyInterval[i] <
205 fA3[j] = fA3[j+1]; << 133 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
206 fA4[j] = fA4[j+1]; << 134 {
>> 135 for(j=i;j<fIntervalNumber;j++)
>> 136 {
>> 137 fEnergyInterval[j] = fEnergyInterval[j+1] ;
>> 138 fA1[j] = fA1[j+1] ;
>> 139 fA2[j] = fA2[j+1] ;
>> 140 fA3[j] = fA3[j+1] ;
>> 141 fA4[j] = fA4[j+1] ;
>> 142 }
>> 143 fIntervalNumber-- ;
>> 144 i-- ;
>> 145 }
207 } 146 }
208 fIntervalNumber--; <<
209 } <<
210 } <<
211 if( fVerbose > 0 ) <<
212 { <<
213 for( i = 1; i <= fIntervalNumber; ++i ) <<
214 { <<
215 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< <<
216 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; <<
217 } <<
218 } <<
219 // Preparation of fSplineEnergy array corres <<
220 <<
221 ComputeLowEnergyCof(material); <<
222 <<
223 G4double betaGammaSqRef = <<
224 fLorentzFactor[fRefGammaNumber]*fLorentzFa <<
225 147
226 NormShift(betaGammaSqRef); << 148 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
227 SplainPAI(betaGammaSqRef); <<
228 149
229 // Preparation of integral PAI cross section << 150 G4double betaGammaSqRef =
230 << 151 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
231 for( i = 1; i <= fSplineNumber; ++i ) <<
232 { <<
233 fDifPAIySection[i] = DifPAIySection(i,bet <<
234 <<
235 if( fVerbose > 0 ) G4cout<<i<<"; dNdxPAI <<
236 } <<
237 IntegralPAIySection(); <<
238 } <<
239 <<
240 ////////////////////////////////////////////// <<
241 // <<
242 // Compute low energy cof. It reduces PAI xsc <<
243 // <<
244 <<
245 void G4PAIySection::ComputeLowEnergyCof(const <<
246 { <<
247 G4int i, numberOfElements = (G4int)material- <<
248 G4double sumZ = 0., sumCof = 0.; <<
249 152
250 static const G4double p0 = 1.20923e+00; << 153 NormShift(betaGammaSqRef) ;
251 static const G4double p1 = 3.53256e-01; << 154 SplainPAI(betaGammaSqRef) ;
252 static const G4double p2 = -1.45052e-03; << 155
253 << 156 // Preparation of integral PAI cross section for input betaGammaSq
254 G4double* thisMaterialZ = new G4double[num <<
255 G4double* thisMaterialCof = new G4double[num <<
256 157
257 for( i = 0; i < numberOfElements; ++i ) << 158 for(i = 1 ; i <= fSplineNumber ; i++)
258 { << 159 {
259 thisMaterialZ[i] = material->GetElement(i) << 160 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
260 sumZ += thisMaterialZ[i]; << 161 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
261 thisMaterialCof[i] = p0+p1*thisMaterialZ[i << 162 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
262 } << 163 }
263 for( i = 0; i < numberOfElements; ++i ) << 164 IntegralPAIySection() ;
264 { << 165 IntegralCerenkov() ;
265 sumCof += thisMaterialCof[i]*thisMaterialZ << 166 IntegralPlasmon() ;
266 } <<
267 fLowEnergyCof = sumCof; <<
268 delete [] thisMaterialZ; <<
269 delete [] thisMaterialCof; <<
270 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof <<
271 } 167 }
272 168
273 ////////////////////////////////////////////// 169 /////////////////////////////////////////////////////////////////////////
274 // 170 //
275 // General control function for class G4PAIySe 171 // General control function for class G4PAIySection
276 // 172 //
277 173
278 void G4PAIySection::InitPAI() 174 void G4PAIySection::InitPAI()
279 { 175 {
280 G4int i; << 176 G4int i ;
281 G4double betaGammaSq = fLorentzFactor[fRefG 177 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
282 fLorentzFactor[fRefG 178 fLorentzFactor[fRefGammaNumber] - 1;
283 179
284 // Preparation of integral PAI cross sectio 180 // Preparation of integral PAI cross section for reference gamma
285 181
286 NormShift(betaGammaSq); << 182 NormShift(betaGammaSq) ;
287 SplainPAI(betaGammaSq); << 183 SplainPAI(betaGammaSq) ;
288 184
289 IntegralPAIySection(); << 185 IntegralPAIySection() ;
290 IntegralCerenkov(); << 186 IntegralCerenkov() ;
291 IntegralPlasmon(); << 187 IntegralPlasmon() ;
292 188
293 for( i = 0; i<= fSplineNumber; ++i) << 189 for(i = 0 ; i<=fSplineNumber ; i++)
294 { 190 {
295 fPAItable[i][fRefGammaNumber] = fIntegral << 191 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i] ;
296 << 192 if(i != 0)
297 if(i != 0) fPAItable[i][0] = fSplineEner << 193 {
>> 194 fPAItable[i][0] = fSplineEnergy[i] ;
>> 195 }
298 } 196 }
299 fPAItable[0][0] = fSplineNumber; << 197 fPAItable[0][0] = fSplineNumber ;
300 198
301 for( G4int j = 1; j < 112; ++j) // fo << 199 for(G4int j = 1 ; j < 112 ; j++) // for other gammas
302 { 200 {
303 if( j == fRefGammaNumber ) continue; << 201 if( j == fRefGammaNumber ) continue ;
304 202
305 betaGammaSq = fLorentzFactor[j]*fLorentz << 203 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1 ;
306 204
307 for(i = 1; i <= fSplineNumber; ++i) << 205 for(i = 1 ; i <= fSplineNumber ; i++)
308 { 206 {
309 fDifPAIySection[i] = DifPAIySection(i 207 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
310 fdNdxCerenkov[i] = PAIdNdxCerenkov( 208 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
311 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 209 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
312 } 210 }
313 IntegralPAIySection(); << 211 IntegralPAIySection() ;
314 IntegralCerenkov(); << 212 IntegralCerenkov() ;
315 IntegralPlasmon(); << 213 IntegralPlasmon() ;
316 214
317 for(i = 0; i <= fSplineNumber; ++i) << 215 for(i = 0 ; i <= fSplineNumber ; i++)
318 { 216 {
319 fPAItable[i][j] = fIntegralPAIySection << 217 fPAItable[i][j] = fIntegralPAIySection[i] ;
320 } 218 }
321 } 219 }
>> 220
322 } 221 }
323 222
324 ////////////////////////////////////////////// 223 ///////////////////////////////////////////////////////////////////////
325 // 224 //
326 // Shifting from borders to intervals Creation 225 // Shifting from borders to intervals Creation of first energy points
327 // 226 //
328 227
329 void G4PAIySection::NormShift(G4double betaGam 228 void G4PAIySection::NormShift(G4double betaGammaSq)
330 { 229 {
331 G4int i, j; << 230 G4int i, j ;
332 231
333 for( i = 1; i <= fIntervalNumber-1; ++i) << 232 for( i = 1 ; i <= fIntervalNumber-1 ; i++ )
334 { 233 {
335 for( j = 1; j <= 2; ++j) << 234 for( j = 1 ; j <= 2 ; j++ )
336 { 235 {
337 fSplineNumber = (i-1)*2 + j; << 236 fSplineNumber = (i-1)*2 + j ;
338 237
339 if( j == 1 ) fSplineEnergy[fSplineNumber 238 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
340 else fSplineEnergy[fSplineNumber 239 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
341 // G4cout<<"cn = "<<fSplineNumber<<"; 240 // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
342 // <<fSplineEnergy[fSplineNumber]<<G4en 241 // <<fSplineEnergy[fSplineNumber]<<G4endl;
343 } 242 }
344 } 243 }
345 fIntegralTerm[1]=RutherfordIntegral(1,fEnerg 244 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
346 245
347 j = 1; << 246 j = 1 ;
348 247
349 for(i=2;i<=fSplineNumber;++i) << 248 for(i=2;i<=fSplineNumber;i++)
350 { 249 {
351 if(fSplineEnergy[i]<fEnergyInterval[j+1]) 250 if(fSplineEnergy[i]<fEnergyInterval[j+1])
352 { 251 {
353 fIntegralTerm[i] = fIntegralTerm[i-1] 252 fIntegralTerm[i] = fIntegralTerm[i-1] +
354 RutherfordIntegral << 253 RutherfordIntegral(j,fSplineEnergy[i-1],
355 << 254 fSplineEnergy[i] ) ;
356 } 255 }
357 else 256 else
358 { 257 {
359 G4double x = RutherfordIntegral(j,fSpli 258 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
360 fEn << 259 fEnergyInterval[j+1] ) ;
361 j++; 260 j++;
362 fIntegralTerm[i] = fIntegralTerm[i-1] 261 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
363 RutherfordIntegral << 262 RutherfordIntegral(j,fEnergyInterval[j],
364 << 263 fSplineEnergy[i] ) ;
365 } 264 }
366 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t" 265 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
367 } 266 }
368 static const G4double nfactor = << 267 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
369 2*pi*pi*hbarc*hbarc*fine_structure_const/e << 268 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber] ;
370 fNormalizationCof = nfactor*fElectronDensity <<
371 269
372 // G4cout<<"fNormalizationCof = "<<fNormaliz << 270 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl ;
373 271
374 // Calculation of PAI differrential cross-se << 272 // Calculation of PAI differrential cross-section (1/(keV*cm))
375 // in the energy points near borders of ener << 273 // in the energy points near borders of energy intervals
376 274
377 for(G4int k=1; k<=fIntervalNumber-1; ++k) << 275 for(G4int k=1;k<=fIntervalNumber-1;k++)
378 { 276 {
379 for(j=1; j<=2; ++j) << 277 for(j=1;j<=2;j++)
380 { 278 {
381 i = (k-1)*2 + j; << 279 i = (k-1)*2 + j ;
382 fImPartDielectricConst[i] = fNormaliz 280 fImPartDielectricConst[i] = fNormalizationCof*
383 ImPartDie << 281 ImPartDielectricConst(k,fSplineEnergy[i]);
384 fRePartDielectricConst[i] = fNormaliz 282 fRePartDielectricConst[i] = fNormalizationCof*
385 RePartDie << 283 RePartDielectricConst(fSplineEnergy[i]);
386 fIntegralTerm[i] *= fNormalizationCof 284 fIntegralTerm[i] *= fNormalizationCof;
387 285
388 fDifPAIySection[i] = DifPAIySection(i 286 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
389 fdNdxCerenkov[i] = PAIdNdxCerenkov( 287 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
390 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 288 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
391 } 289 }
392 } 290 }
393 291
394 } // end of NormShift 292 } // end of NormShift
395 293
396 ////////////////////////////////////////////// 294 /////////////////////////////////////////////////////////////////////////
397 // 295 //
398 // Creation of new energy points as geometrica 296 // Creation of new energy points as geometrical mean of existing
399 // one, calculation PAI_cs for them, while the 297 // one, calculation PAI_cs for them, while the error of logarithmic
400 // linear approximation would be smaller than 298 // linear approximation would be smaller than 'fError'
401 299
402 void G4PAIySection::SplainPAI(G4double betaGam 300 void G4PAIySection::SplainPAI(G4double betaGammaSq)
403 { 301 {
404 G4int k = 1; << 302 G4int k = 1 ;
405 G4int i = 1; << 303 G4int i = 1 ;
406 304
407 while ( (i < fSplineNumber) && (fSplineNumb 305 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
408 { 306 {
409 if(fSplineEnergy[i+1] > fEnergyInterval[ 307 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
410 { 308 {
411 k++; // Here next energy point is << 309 k++ ; // Here next energy point is in next energy interval
412 ++i; << 310 i++;
413 continue; 311 continue;
414 } 312 }
415 // Shifting of arrayes for inserting the << 313 // Shifting of arrayes for inserting the geometrical
416 // average of 'i' and 'i+1' energy point << 314 // average of 'i' and 'i+1' energy points to 'i+1' place
417 fSplineNumber++; 315 fSplineNumber++;
418 316
419 for(G4int j = fSplineNumber; j >= i+2; j << 317 for(G4int j = fSplineNumber; j >= i+2 ; j-- )
420 { 318 {
421 fSplineEnergy[j] = fSplineEn 319 fSplineEnergy[j] = fSplineEnergy[j-1];
422 fImPartDielectricConst[j] = fImPartDi 320 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
423 fRePartDielectricConst[j] = fRePartDi << 321 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
424 fIntegralTerm[j] = fIntegral << 322 fIntegralTerm[j] = fIntegralTerm[j-1];
425 323
426 fDifPAIySection[j] = fDifPAIySection[ << 324 fDifPAIySection[j] = fDifPAIySection[j-1];
427 fdNdxCerenkov[j] = fdNdxCerenkov[j- 325 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
428 fdNdxPlasmon[j] = fdNdxPlasmon[j-1 326 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
429 } 327 }
430 G4double x1 = fSplineEnergy[i]; 328 G4double x1 = fSplineEnergy[i];
431 G4double x2 = fSplineEnergy[i+1]; 329 G4double x2 = fSplineEnergy[i+1];
432 G4double yy1 = fDifPAIySection[i]; 330 G4double yy1 = fDifPAIySection[i];
433 G4double y2 = fDifPAIySection[i+1]; 331 G4double y2 = fDifPAIySection[i+1];
434 332
435 G4double en1 = sqrt(x1*x2); 333 G4double en1 = sqrt(x1*x2);
436 fSplineEnergy[i+1] = en1; 334 fSplineEnergy[i+1] = en1;
437 335
438 // Calculation of logarithmic linear app << 336 // Calculation of logarithmic linear approximation
439 // in this (enr) energy point, which num << 337 // in this (enr) energy point, which number is 'i+1' now
440 338
441 G4double a = log10(y2/yy1)/log10(x2/x1); 339 G4double a = log10(y2/yy1)/log10(x2/x1);
442 G4double b = log10(yy1) - a*log10(x1); 340 G4double b = log10(yy1) - a*log10(x1);
443 G4double y = a*log10(en1) + b; << 341 G4double y = a*log10(en1) + b ;
444 y = pow(10.,y); 342 y = pow(10.,y);
445 343
446 // Calculation of the PAI dif. cross-sec << 344 // Calculation of the PAI dif. cross-section at this point
447 345
448 fImPartDielectricConst[i+1] = fNormaliza 346 fImPartDielectricConst[i+1] = fNormalizationCof*
449 ImPartDiel << 347 ImPartDielectricConst(k,fSplineEnergy[i+1]);
450 fRePartDielectricConst[i+1] = fNormaliza 348 fRePartDielectricConst[i+1] = fNormalizationCof*
451 RePartDiel << 349 RePartDielectricConst(fSplineEnergy[i+1]);
452 fIntegralTerm[i+1] = fIntegralTerm[i] + 350 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
453 RutherfordIntegral( << 351 RutherfordIntegral(k,fSplineEnergy[i],
454 352 fSplineEnergy[i+1]);
455 353
456 fDifPAIySection[i+1] = DifPAIySection(i+ 354 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
457 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i 355 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
458 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+ 356 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
459 357
460 // Condition for next divisi << 358 // Condition for next division of this segment or to pass
461 // to higher energies << 359 // to higher energies
462 360
463 G4double x = 2*(fDifPAIySection[i+1] - y 361 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
464 362
465 G4double delta = 2.*(fSplineEnergy[i+1]- <<
466 /(fSplineEnergy[i+1]+fSplineEnergy[i]) <<
467 <<
468 if( x < 0 ) 363 if( x < 0 )
469 { 364 {
470 x = -x; << 365 x = -x ;
471 } 366 }
472 if( x > fError && fSplineNumber < fMaxSp << 367 if( x > fError && fSplineNumber < fMaxSplineSize-1 )
473 { 368 {
474 continue; // next division << 369 continue; // next division
475 } 370 }
476 i += 2; // pass to next segment 371 i += 2; // pass to next segment
477 372
478 // Loop checking, 03-Aug-2015, Vladimir <<
479 } // close 'while' 373 } // close 'while'
480 374
481 } // end of SplainPAI 375 } // end of SplainPAI
482 376
483 377
484 ////////////////////////////////////////////// 378 ////////////////////////////////////////////////////////////////////
485 // 379 //
486 // Integration over electrons that could be co 380 // Integration over electrons that could be considered
487 // quasi-free at energy transfer of interest 381 // quasi-free at energy transfer of interest
488 382
489 G4double G4PAIySection::RutherfordIntegral( G4 383 G4double G4PAIySection::RutherfordIntegral( G4int k,
490 G4 << 384 G4double x1,
491 << 385 G4double x2 )
492 { 386 {
493 G4double c1, c2, c3; << 387 G4double c1, c2, c3 ;
494 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<< << 388 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
495 G4double x12 = x1*x2; << 389 c1 = (x2 - x1)/x1/x2 ;
496 c1 = (x2 - x1)/x12; << 390 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
497 c2 = (x2 - x1)*(x2 + x1)/(x12*x12); << 391 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
498 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12 <<
499 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "< 392 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
500 393
501 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3 << 394 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3 ;
502 395
503 } // end of RutherfordIntegral 396 } // end of RutherfordIntegral
504 397
505 398
506 ////////////////////////////////////////////// 399 /////////////////////////////////////////////////////////////////
507 // 400 //
508 // Imaginary part of dielectric constant 401 // Imaginary part of dielectric constant
509 // (G4int k - interval number, G4double en1 - 402 // (G4int k - interval number, G4double en1 - energy point)
510 403
511 G4double G4PAIySection::ImPartDielectricConst( << 404 G4double G4PAIySection::ImPartDielectricConst( G4int k ,
>> 405 G4double energy1 )
512 { 406 {
513 G4double energy2,energy3,energy4,result; 407 G4double energy2,energy3,energy4,result;
514 408
515 energy2 = energy1*energy1; 409 energy2 = energy1*energy1;
516 energy3 = energy2*energy1; 410 energy3 = energy2*energy1;
517 energy4 = energy3*energy1; 411 energy4 = energy3*energy1;
518 412
519 result = fA1[k]/energy1+fA2[k]/energy2+fA3[ << 413 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4 ;
520 result *=hbarc/energy1; << 414 result *=hbarc/energy1 ;
521 415
522 return result; << 416 return result ;
523 417
524 } // end of ImPartDielectricConst 418 } // end of ImPartDielectricConst
525 419
526 420
527 ////////////////////////////////////////////// 421 //////////////////////////////////////////////////////////////////////////////
528 // 422 //
529 // Real part of dielectric constant minus unit 423 // Real part of dielectric constant minus unit: epsilon_1 - 1
530 // (G4double enb - energy point) 424 // (G4double enb - energy point)
531 // 425 //
532 426
533 G4double G4PAIySection::RePartDielectricConst( 427 G4double G4PAIySection::RePartDielectricConst(G4double enb)
534 { 428 {
535 G4double x0, x02, x03, x04, x05, x1, x2, xx 429 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
536 c1, c2, c3, cof1, cof2, xln1, xln2 << 430 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
537 431
538 x0 = enb; << 432 x0 = enb ;
539 result = 0; << 433 result = 0 ;
540 434
541 for(G4int i=1;i<=fIntervalNumber-1;++i) << 435 for(G4int i=1;i<=fIntervalNumber-1;i++)
542 { 436 {
543 x1 = fEnergyInterval[i]; << 437 x1 = fEnergyInterval[i] ;
544 x2 = fEnergyInterval[i+1]; << 438 x2 = fEnergyInterval[i+1] ;
545 xx1 = x1 - x0; << 439 xx1 = x1 - x0 ;
546 xx2 = x2 - x0; << 440 xx2 = x2 - x0 ;
547 xx12 = xx2/xx1; << 441 xx12 = xx2/xx1 ;
548 442
549 if(xx12<0.) << 443 if(xx12<0)
550 { 444 {
551 xx12 = -xx12; << 445 xx12 = -xx12;
552 } 446 }
553 xln1 = log(x2/x1); << 447 xln1 = log(x2/x1) ;
554 xln2 = log(xx12); << 448 xln2 = log(xx12) ;
555 xln3 = log((x2 + x0)/(x1 + x0)); << 449 xln3 = log((x2 + x0)/(x1 + x0)) ;
556 x02 = x0*x0; << 450 x02 = x0*x0 ;
557 x03 = x02*x0; << 451 x03 = x02*x0 ;
558 x04 = x03*x0; << 452 x04 = x03*x0 ;
559 x05 = x04*x0; 453 x05 = x04*x0;
560 G4double x12 = x1*x2; << 454 c1 = (x2 - x1)/x1/x2 ;
561 c1 = (x2 - x1)/x12; << 455 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
562 c2 = (x2 - x1)*(x2 +x1)/(x12*x12); << 456 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
563 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/( << 457
564 << 458 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 ;
565 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 << 459 result -= (fA2[i]/x02 + fA4[i]/x04)*c1 ;
566 result -= (fA2[i]/x02 + fA4[i]/x04)*c1; << 460 result -= fA3[i]*c2/2/x02 ;
567 result -= fA3[i]*c2/2/x02; << 461 result -= fA4[i]*c3/3/x02 ;
568 result -= fA4[i]*c3/3/x02; <<
569 462
570 cof1 = fA1[i]/x02 + fA3[i]/x04; << 463 cof1 = fA1[i]/x02 + fA3[i]/x04 ;
571 cof2 = fA2[i]/x03 + fA4[i]/x05; << 464 cof2 = fA2[i]/x03 + fA4[i]/x05 ;
572 465
573 result += 0.5*(cof1 +cof2)*xln2; << 466 result += 0.5*(cof1 +cof2)*xln2 ;
574 result += 0.5*(cof1 - cof2)*xln3; << 467 result += 0.5*(cof1 - cof2)*xln3 ;
575 } 468 }
576 result *= 2*hbarc/pi; << 469 result *= 2*hbarc/pi ;
577 470
578 return result; << 471 return result ;
579 472
580 } // end of RePartDielectricConst 473 } // end of RePartDielectricConst
581 474
582 ////////////////////////////////////////////// 475 //////////////////////////////////////////////////////////////////////
583 // 476 //
584 // PAI differential cross-section in terms of 477 // PAI differential cross-section in terms of
585 // simplified Allison's equation 478 // simplified Allison's equation
586 // 479 //
587 480
588 G4double G4PAIySection::DifPAIySection( G4int 481 G4double G4PAIySection::DifPAIySection( G4int i ,
589 G4doub 482 G4double betaGammaSq )
590 { 483 {
591 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7 << 484 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
592 be2 = betaGammaSq/(1 + betaGammaSq); << 485 //G4double beta, be4 ;
593 beta = std::sqrt(be2); << 486 G4double be4 ;
594 cof = 1; << 487 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
595 x1 = log(2*electron_mass_c2/fSplineEnergy[i << 488 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
>> 489 be2 = betaGammaSq/(1 + betaGammaSq) ;
>> 490 be4 = be2*be2 ;
>> 491 // beta = sqrt(be2) ;
>> 492 cof = 1 ;
>> 493 x1 = log(2*electron_mass_c2/fSplineEnergy[i]) ;
596 494
597 if( betaGammaSq < 0.01 ) x2 = log(be2); << 495 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
598 else 496 else
599 { 497 {
600 x2 = -log( (1/betaGammaSq - fRePartDielec 498 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
601 (1/betaGammaSq - fRePartDielec << 499 (1/betaGammaSq - fRePartDielectricConst[i]) +
602 fImPartDielectricConst[i]*fImP << 500 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2 ;
603 } 501 }
604 if( fImPartDielectricConst[i] == 0.0 ||beta 502 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
605 { 503 {
606 x6=0; << 504 x6=0 ;
607 } 505 }
608 else 506 else
609 { 507 {
610 x3 = -fRePartDielectricConst[i] + 1/betaG << 508 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq ;
611 x5 = -1 - fRePartDielectricConst[i] + 509 x5 = -1 - fRePartDielectricConst[i] +
612 be2*((1 +fRePartDielectricConst[i])* 510 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
613 fImPartDielectricConst[i]*fImPartDie << 511 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
614 512
615 x7 = std::atan2(fImPartDielectricConst[i] << 513 x7 = atan2(fImPartDielectricConst[i],x3) ;
616 x6 = x5 * x7; << 514 x6 = x5 * x7 ;
617 } 515 }
618 x4 = ((x1 + x2)*fImPartDielectricConst[i] + << 516 // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
>> 517
>> 518 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc ;
>> 519 // if( x4 < 0.0 ) x4 = 0.0 ;
619 x8 = (1 + fRePartDielectricConst[i])*(1 + f 520 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
620 fImPartDielectricConst[i]*fImPartDiele << 521 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
621 <<
622 result = (x4 + cof*fIntegralTerm[i]/fSpline <<
623 result = std::max(result, 1.0e-8); <<
624 result *= fine_structure_const/(be2*pi); <<
625 // low energy correction <<
626 522
627 G4double lowCof = fLowEnergyCof; // 6.0 ; / << 523 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]) ;
628 << 524 if(result < 1.0e-8) result = 1.0e-8 ;
629 result *= (1 - std::exp(-beta/(betaBohr*low << 525 result *= fine_structure_const/be2/pi ;
>> 526 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
>> 527 // result *= (1-exp(-be2/betaBohr2)) ;
>> 528 result *= (1-exp(-be4/betaBohr4)) ;
>> 529 // if(fDensity >= 0.1)
630 if(x8 > 0.) 530 if(x8 > 0.)
631 { 531 {
632 result /= x8; << 532 result /= x8 ;
633 } 533 }
634 return result; << 534 return result ;
635 535
636 } // end of DifPAIySection 536 } // end of DifPAIySection
637 537
638 ////////////////////////////////////////////// 538 //////////////////////////////////////////////////////////////////////////
639 // 539 //
640 // Calculation od dN/dx of collisions with cre 540 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
641 541
642 G4double G4PAIySection::PAIdNdxCerenkov( G4int << 542 G4double G4PAIySection::PAIdNdxCerenkov( G4int i ,
>> 543 G4double betaGammaSq )
643 { 544 {
644 G4double logarithm, x3, x5, argument, modul << 545 G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ;
645 G4double be2, be4; << 546 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
>> 547
>> 548 cof = 1.0 ;
>> 549 cofBetaBohr = 4.0 ;
>> 550 betaBohr2 = fine_structure_const*fine_structure_const ;
>> 551 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
646 552
647 be2 = betaGammaSq/(1 + betaGammaSq); << 553 be2 = betaGammaSq/(1 + betaGammaSq) ;
648 be4 = be2*be2; << 554 be4 = be2*be2 ;
649 555
650 if( betaGammaSq < 0.01 ) logarithm = log(1. << 556 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
651 else 557 else
652 { 558 {
653 logarithm = -std::log( (1/betaGammaSq - f << 559 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
654 (1/betaGammaSq - fRePa << 560 (1/betaGammaSq - fRePartDielectricConst[i]) +
655 fImPartDielectricConst << 561 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5 ;
656 logarithm += std::log(1+1.0/betaGammaSq); << 562 logarithm += log(1+1.0/betaGammaSq) ;
657 } 563 }
658 564
659 if( fImPartDielectricConst[i] == 0.0 || bet 565 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
660 { 566 {
661 argument = 0.0; << 567 argument = 0.0 ;
662 } 568 }
663 else 569 else
664 { 570 {
665 x3 = -fRePartDielectricConst[i] + 1.0/bet << 571 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq ;
666 x5 = -1.0 - fRePartDielectricConst[i] + 572 x5 = -1.0 - fRePartDielectricConst[i] +
667 be2*((1.0 +fRePartDielectricConst[i] 573 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
668 fImPartDielectricConst[i]*fImPartDie << 574 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
669 if( x3 == 0.0 ) argument = 0.5*pi; 575 if( x3 == 0.0 ) argument = 0.5*pi;
670 else argument = std::atan2(fIm << 576 else argument = atan2(fImPartDielectricConst[i],x3) ;
671 argument *= x5 ; << 577 argument *= x5 ;
672 } 578 }
673 dNdxC = ( logarithm*fImPartDielectricConst[ << 579 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc ;
674 580
675 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; << 581 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
676 582
677 dNdxC *= fine_structure_const/be2/pi; << 583 dNdxC *= fine_structure_const/be2/pi ;
678 584
679 dNdxC *= (1 - std::exp(-be4/betaBohr4)); << 585 dNdxC *= (1-exp(-be4/betaBohr4)) ;
680 586
>> 587 // if(fDensity >= 0.1)
>> 588 // {
681 modul2 = (1.0 + fRePartDielectricConst[i])* 589 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
682 fImPartDielectricConst[i]* << 590 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
683 if(modul2 > 0.) 591 if(modul2 > 0.)
684 { 592 {
685 dNdxC /= modul2; << 593 dNdxC /= modul2 ;
686 } 594 }
687 return dNdxC; << 595 return dNdxC ;
688 596
689 } // end of PAIdNdxCerenkov 597 } // end of PAIdNdxCerenkov
690 598
691 ////////////////////////////////////////////// 599 //////////////////////////////////////////////////////////////////////////
692 // 600 //
693 // Calculation od dN/dx of collisions with cre 601 // Calculation od dN/dx of collisions with creation of longitudinal EM
694 // excitations (plasmons, delta-electrons) 602 // excitations (plasmons, delta-electrons)
695 603
696 G4double G4PAIySection::PAIdNdxPlasmon( G4int << 604 G4double G4PAIySection::PAIdNdxPlasmon( G4int i ,
>> 605 G4double betaGammaSq )
697 { 606 {
698 G4double cof, resonance, modul2, dNdxP; << 607 G4double cof, resonance, modul2, dNdxP ;
699 G4double be2, be4; << 608 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
700 609
701 cof = 1; << 610 cof = 1 ;
>> 611 cofBetaBohr = 4.0 ;
>> 612 betaBohr2 = fine_structure_const*fine_structure_const ;
>> 613 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
702 614
703 be2 = betaGammaSq/(1 + betaGammaSq); << 615 be2 = betaGammaSq/(1 + betaGammaSq) ;
704 be4 = be2*be2; << 616 be4 = be2*be2 ;
705 617
706 resonance = std::log(2*electron_mass_c2*be2 << 618 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]) ;
707 resonance *= fImPartDielectricConst[i]/hbar << 619 resonance *= fImPartDielectricConst[i]/hbarc ;
708 620
709 dNdxP = ( resonance + cof*fIntegralTerm[i]/ <<
710 621
711 dNdxP = std::max(dNdxP, 1.0e-8); << 622 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ) ;
712 623
713 dNdxP *= fine_structure_const/be2/pi; << 624 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
714 dNdxP *= (1 - std::exp(-be4/betaBohr4)); <<
715 625
>> 626 dNdxP *= fine_structure_const/be2/pi ;
>> 627 dNdxP *= (1-exp(-be4/betaBohr4)) ;
>> 628
>> 629 // if( fDensity >= 0.1 )
>> 630 // {
716 modul2 = (1 + fRePartDielectricConst[i])*(1 631 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
717 fImPartDielectricConst[i]*fImPartDielectr << 632 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
718 if(modul2 > 0.) 633 if(modul2 > 0.)
719 { 634 {
720 dNdxP /= modul2; << 635 dNdxP /= modul2 ;
721 } 636 }
722 return dNdxP; << 637 return dNdxP ;
723 638
724 } // end of PAIdNdxPlasmon 639 } // end of PAIdNdxPlasmon
725 640
726 ////////////////////////////////////////////// 641 ////////////////////////////////////////////////////////////////////////
727 // 642 //
728 // Calculation of the PAI integral cross-secti 643 // Calculation of the PAI integral cross-section
729 // fIntegralPAIySection[1] = specific primary 644 // fIntegralPAIySection[1] = specific primary ionisation, 1/cm
730 // and fIntegralPAIySection[0] = mean energy l 645 // and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm
731 646
732 void G4PAIySection::IntegralPAIySection() 647 void G4PAIySection::IntegralPAIySection()
733 { 648 {
734 fIntegralPAIySection[fSplineNumber] = 0; << 649 fIntegralPAIySection[fSplineNumber] = 0 ;
735 fIntegralPAIdEdx[fSplineNumber] = 0; << 650 fIntegralPAIdEdx[fSplineNumber] = 0 ;
736 fIntegralPAIySection[0] = 0; << 651 fIntegralPAIySection[0] = 0 ;
737 G4int k = fIntervalNumber -1; << 652 G4int k = fIntervalNumber -1 ;
738 653
739 for(G4int i = fSplineNumber-1; i >= 1; i--) << 654 for(G4int i = fSplineNumber-1 ; i >= 1 ; i--)
740 { 655 {
741 if(fSplineEnergy[i] >= fEnergyInterval[k]) 656 if(fSplineEnergy[i] >= fEnergyInterval[k])
742 { 657 {
743 fIntegralPAIySection[i] = fIntegralPAIyS << 658 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i) ;
744 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i << 659 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i) ;
745 } 660 }
746 else 661 else
747 { 662 {
748 fIntegralPAIySection[i] = fIntegralPAIyS 663 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
749 SumOverBord << 664 SumOverBorder(i+1,fEnergyInterval[k]) ;
750 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i 665 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
751 SumOverBord << 666 SumOverBorderdEdx(i+1,fEnergyInterval[k]) ;
752 k--; << 667 k-- ;
753 } 668 }
754 } 669 }
755 } // end of IntegralPAIySection 670 } // end of IntegralPAIySection
756 671
757 ////////////////////////////////////////////// 672 ////////////////////////////////////////////////////////////////////////
758 // 673 //
759 // Calculation of the PAI Cerenkov integral cr 674 // Calculation of the PAI Cerenkov integral cross-section
760 // fIntegralCrenkov[1] = specific Crenkov ioni 675 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
761 // and fIntegralCerenkov[0] = mean Cerenkov lo 676 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
762 677
763 void G4PAIySection::IntegralCerenkov() 678 void G4PAIySection::IntegralCerenkov()
764 { 679 {
765 G4int i, k; << 680 G4int i, k ;
766 fIntegralCerenkov[fSplineNumber] = 0; << 681 fIntegralCerenkov[fSplineNumber] = 0 ;
767 fIntegralCerenkov[0] = 0; << 682 fIntegralCerenkov[0] = 0 ;
768 k = fIntervalNumber -1; << 683 k = fIntervalNumber -1 ;
769 684
770 for( i = fSplineNumber-1; i >= 1; i-- ) << 685 for( i = fSplineNumber-1 ; i >= 1 ; i-- )
771 { 686 {
772 if(fSplineEnergy[i] >= fEnergyInterval[k 687 if(fSplineEnergy[i] >= fEnergyInterval[k])
773 { 688 {
774 fIntegralCerenkov[i] = fIntegralCerenk << 689 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i) ;
775 // G4cout<<"int: i = "<<i<<"; sumC = " << 690 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
776 } 691 }
777 else 692 else
778 { 693 {
779 fIntegralCerenkov[i] = fIntegralCerenk 694 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
780 SumOverBord << 695 SumOverBordCerenkov(i+1,fEnergyInterval[k]) ;
781 k--; << 696 k-- ;
782 // G4cout<<"bord: i = "<<i<<"; sumC = << 697 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
783 } 698 }
784 } 699 }
785 700
786 } // end of IntegralCerenkov 701 } // end of IntegralCerenkov
787 702
788 ////////////////////////////////////////////// 703 ////////////////////////////////////////////////////////////////////////
789 // 704 //
790 // Calculation of the PAI Plasmon integral cro 705 // Calculation of the PAI Plasmon integral cross-section
791 // fIntegralPlasmon[1] = splasmon primary ioni 706 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
792 // and fIntegralPlasmon[0] = mean plasmon loss 707 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
793 708
794 void G4PAIySection::IntegralPlasmon() 709 void G4PAIySection::IntegralPlasmon()
795 { 710 {
796 fIntegralPlasmon[fSplineNumber] = 0; << 711 fIntegralPlasmon[fSplineNumber] = 0 ;
797 fIntegralPlasmon[0] = 0; << 712 fIntegralPlasmon[0] = 0 ;
798 G4int k = fIntervalNumber -1; << 713 G4int k = fIntervalNumber -1 ;
799 for(G4int i=fSplineNumber-1;i>=1;i--) 714 for(G4int i=fSplineNumber-1;i>=1;i--)
800 { 715 {
801 if(fSplineEnergy[i] >= fEnergyInterval[k 716 if(fSplineEnergy[i] >= fEnergyInterval[k])
802 { 717 {
803 fIntegralPlasmon[i] = fIntegralPlasmon << 718 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i) ;
804 } 719 }
805 else 720 else
806 { 721 {
807 fIntegralPlasmon[i] = fIntegralPlasmon 722 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
808 SumOverBord << 723 SumOverBordPlasmon(i+1,fEnergyInterval[k]) ;
809 k--; << 724 k-- ;
810 } 725 }
811 } 726 }
>> 727
812 } // end of IntegralPlasmon 728 } // end of IntegralPlasmon
813 729
814 ////////////////////////////////////////////// 730 //////////////////////////////////////////////////////////////////////
815 // 731 //
816 // Calculation the PAI integral cross-section 732 // Calculation the PAI integral cross-section inside
817 // of interval of continuous values of photo-i 733 // of interval of continuous values of photo-ionisation
818 // cross-section. Parameter 'i' is the number 734 // cross-section. Parameter 'i' is the number of interval.
819 735
820 G4double G4PAIySection::SumOverInterval( G4int 736 G4double G4PAIySection::SumOverInterval( G4int i )
821 { 737 {
822 G4double x0,x1,y0,yy1,a,b,c,result; << 738 G4double x0,x1,y0,yy1,a,b,c,result ;
823 <<
824 x0 = fSplineEnergy[i]; <<
825 x1 = fSplineEnergy[i+1]; <<
826 <<
827 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) <<
828 739
829 y0 = fDifPAIySection[i]; << 740 x0 = fSplineEnergy[i] ;
>> 741 x1 = fSplineEnergy[i+1] ;
>> 742 y0 = fDifPAIySection[i] ;
830 yy1 = fDifPAIySection[i+1]; 743 yy1 = fDifPAIySection[i+1];
831 //G4cout << "## x0= " << x0 << " x1= " << x <<
832 c = x1/x0; 744 c = x1/x0;
833 //G4cout << "c= " << c << " y0= " << y0 << << 745 a = log10(yy1/y0)/log10(c) ;
834 a = log10(yy1/y0)/log10(c); << 746 // b = log10(y0) - a*log10(x0) ;
835 //G4cout << "a= " << a << G4endl; << 747 b = y0/pow(x0,a) ;
836 << 748 a += 1 ;
837 b = 0.0; <<
838 if(a < 20.) b = y0/pow(x0,a); <<
839 <<
840 a += 1; <<
841 if(a == 0) 749 if(a == 0)
842 { 750 {
843 result = b*log(x1/x0); << 751 result = b*log(x1/x0) ;
844 } 752 }
845 else 753 else
846 { 754 {
847 result = y0*(x1*pow(c,a-1) - x0)/a; << 755 result = y0*(x1*pow(c,a-1) - x0)/a ;
848 } 756 }
849 a++; 757 a++;
850 if(a == 0) 758 if(a == 0)
851 { 759 {
852 fIntegralPAIySection[0] += b*log(x1/x0); << 760 fIntegralPAIySection[0] += b*log(x1/x0) ;
853 } 761 }
854 else 762 else
855 { 763 {
856 fIntegralPAIySection[0] += y0*(x1*x1*pow << 764 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
857 } 765 }
858 return result; << 766 return result ;
859 767
860 } // end of SumOverInterval 768 } // end of SumOverInterval
861 769
862 ///////////////////////////////// 770 /////////////////////////////////
863 771
864 G4double G4PAIySection::SumOverIntervaldEdx( G 772 G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
865 { 773 {
866 G4double x0,x1,y0,yy1,a,b,c,result; << 774 G4double x0,x1,y0,yy1,a,b,c,result ;
867 <<
868 x0 = fSplineEnergy[i]; <<
869 x1 = fSplineEnergy[i+1]; <<
870 775
871 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 776 x0 = fSplineEnergy[i] ;
872 << 777 x1 = fSplineEnergy[i+1] ;
873 y0 = fDifPAIySection[i]; << 778 y0 = fDifPAIySection[i] ;
874 yy1 = fDifPAIySection[i+1]; 779 yy1 = fDifPAIySection[i+1];
875 c = x1/x0; 780 c = x1/x0;
876 a = log10(yy1/y0)/log10(c); << 781 a = log10(yy1/y0)/log10(c) ;
877 << 782 // b = log10(y0) - a*log10(x0) ;
878 b = 0.0; << 783 b = y0/pow(x0,a) ;
879 if(a < 20.) b = y0/pow(x0,a); << 784 a += 2 ;
880 <<
881 a += 2; <<
882 if(a == 0) 785 if(a == 0)
883 { 786 {
884 result = b*log(x1/x0); << 787 result = b*log(x1/x0) ;
885 } 788 }
886 else 789 else
887 { 790 {
888 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a; << 791 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
889 } 792 }
890 return result; << 793 return result ;
891 794
892 } // end of SumOverInterval 795 } // end of SumOverInterval
893 796
894 ////////////////////////////////////////////// 797 //////////////////////////////////////////////////////////////////////
895 // 798 //
896 // Calculation the PAI Cerenkov integral cross 799 // Calculation the PAI Cerenkov integral cross-section inside
897 // of interval of continuous values of photo-i 800 // of interval of continuous values of photo-ionisation Cerenkov
898 // cross-section. Parameter 'i' is the number 801 // cross-section. Parameter 'i' is the number of interval.
899 802
900 G4double G4PAIySection::SumOverInterCerenkov( 803 G4double G4PAIySection::SumOverInterCerenkov( G4int i )
901 { 804 {
902 G4double x0,x1,y0,yy1,a,c,result; << 805 G4double x0,x1,y0,yy1,a,c,result ;
903 <<
904 x0 = fSplineEnergy[i]; <<
905 x1 = fSplineEnergy[i+1]; <<
906 806
907 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 807 x0 = fSplineEnergy[i] ;
908 << 808 x1 = fSplineEnergy[i+1] ;
909 y0 = fdNdxCerenkov[i]; << 809 y0 = fdNdxCerenkov[i] ;
910 yy1 = fdNdxCerenkov[i+1]; 810 yy1 = fdNdxCerenkov[i+1];
911 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<" 811 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
912 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4en 812 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
913 813
914 c = x1/x0; 814 c = x1/x0;
915 a = log10(yy1/y0)/log10(c); << 815 a = log10(yy1/y0)/log10(c) ;
916 G4double b = 0.0; 816 G4double b = 0.0;
917 if(a < 20.) b = y0/pow(x0,a); << 817 if(a < 20.) b = y0/pow(x0,a) ;
918 818
919 a += 1.0; << 819 a += 1.0 ;
920 if(a == 0) result = b*log(c); << 820 if(a == 0) result = b*log(c) ;
921 else result = y0*(x1*pow(c,a-1) - x0) << 821 else result = y0*(x1*pow(c,a-1) - x0)/a ;
922 a += 1.0; << 822 a += 1.0 ;
923 823
924 if( a == 0 ) fIntegralCerenkov[0] += b*log( << 824 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0) ;
925 else fIntegralCerenkov[0] += y0*(x1 << 825 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
926 // G4cout<<"a = "<<a<<"; b = "<<b<<"; resu 826 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
927 return result; << 827 return result ;
928 828
929 } // end of SumOverInterCerenkov 829 } // end of SumOverInterCerenkov
930 830
931 ////////////////////////////////////////////// 831 //////////////////////////////////////////////////////////////////////
932 // 832 //
933 // Calculation the PAI Plasmon integral cross- 833 // Calculation the PAI Plasmon integral cross-section inside
934 // of interval of continuous values of photo-i 834 // of interval of continuous values of photo-ionisation Plasmon
935 // cross-section. Parameter 'i' is the number 835 // cross-section. Parameter 'i' is the number of interval.
936 836
937 G4double G4PAIySection::SumOverInterPlasmon( G 837 G4double G4PAIySection::SumOverInterPlasmon( G4int i )
938 { 838 {
939 G4double x0,x1,y0,yy1,a,c,result; << 839 G4double x0,x1,y0,yy1,a,c,result ;
940 <<
941 x0 = fSplineEnergy[i]; <<
942 x1 = fSplineEnergy[i+1]; <<
943 <<
944 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) <<
945 840
946 y0 = fdNdxPlasmon[i]; << 841 x0 = fSplineEnergy[i] ;
>> 842 x1 = fSplineEnergy[i+1] ;
>> 843 y0 = fdNdxPlasmon[i] ;
947 yy1 = fdNdxPlasmon[i+1]; 844 yy1 = fdNdxPlasmon[i+1];
948 c = x1/x0; << 845 c =x1/x0;
949 a = log10(yy1/y0)/log10(c); << 846 a = log10(yy1/y0)/log10(c) ;
950 847
951 G4double b = 0.0; 848 G4double b = 0.0;
952 if(a < 20.) b = y0/pow(x0,a); << 849 if(a < 20.) b = y0/pow(x0,a) ;
953 850
954 a += 1.0; << 851 a += 1.0 ;
955 if(a == 0) result = b*log(x1/x0); << 852 if(a == 0) result = b*log(x1/x0) ;
956 else result = y0*(x1*pow(c,a-1) - x0) << 853 else result = y0*(x1*pow(c,a-1) - x0)/a ;
957 a += 1.0; << 854 a += 1.0 ;
958 855
959 if( a == 0 ) fIntegralPlasmon[0] += b*log(x << 856 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0) ;
960 else fIntegralPlasmon[0] += y0*(x1* << 857 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
961 858
962 return result; << 859 return result ;
963 860
964 } // end of SumOverInterPlasmon 861 } // end of SumOverInterPlasmon
965 862
966 ////////////////////////////////////////////// 863 ///////////////////////////////////////////////////////////////////////////////
967 // 864 //
968 // Integration of PAI cross-section for the ca 865 // Integration of PAI cross-section for the case of
969 // passing across border between intervals 866 // passing across border between intervals
970 867
971 G4double G4PAIySection::SumOverBorder( G4int 868 G4double G4PAIySection::SumOverBorder( G4int i ,
972 G4doubl 869 G4double en0 )
973 { 870 {
974 G4double x0,x1,y0,yy1,a,d,e0,result; << 871 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
975 872
976 e0 = en0; << 873 e0 = en0 ;
977 x0 = fSplineEnergy[i]; << 874 x0 = fSplineEnergy[i] ;
978 x1 = fSplineEnergy[i+1]; << 875 x1 = fSplineEnergy[i+1] ;
979 y0 = fDifPAIySection[i]; << 876 y0 = fDifPAIySection[i] ;
980 yy1 = fDifPAIySection[i+1]; << 877 yy1 = fDifPAIySection[i+1] ;
981 878
>> 879 c = x1/x0;
982 d = e0/x0; 880 d = e0/x0;
983 a = log10(yy1/y0)/log10(x1/x0); << 881 a = log10(yy1/y0)/log10(x1/x0) ;
984 882
985 G4double b = 0.0; 883 G4double b = 0.0;
986 if(a < 20.) b = y0/pow(x0,a); << 884 if(a < 20.) b = y0/pow(x0,a) ;
987 885
988 a += 1; << 886 a += 1 ;
989 if(a == 0) 887 if(a == 0)
990 { 888 {
991 result = b*log(x0/e0); << 889 result = b*log(x0/e0) ;
992 } 890 }
993 else 891 else
994 { 892 {
995 result = y0*(x0 - e0*pow(d,a-1))/a; << 893 result = y0*(x0 - e0*pow(d,a-1))/a ;
996 } 894 }
997 a++; << 895 a++ ;
998 if(a == 0) 896 if(a == 0)
999 { 897 {
1000 fIntegralPAIySection[0] += b*log(x0/e0) << 898 fIntegralPAIySection[0] += b*log(x0/e0) ;
1001 } 899 }
1002 else 900 else
1003 { 901 {
1004 fIntegralPAIySection[0] += y0*(x0*x0 - << 902 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1005 } 903 }
1006 x0 = fSplineEnergy[i - 1]; << 904 x0 = fSplineEnergy[i - 1] ;
1007 x1 = fSplineEnergy[i - 2]; << 905 x1 = fSplineEnergy[i - 2] ;
1008 y0 = fDifPAIySection[i - 1]; << 906 y0 = fDifPAIySection[i - 1] ;
1009 yy1 = fDifPAIySection[i - 2]; << 907 yy1 = fDifPAIySection[i - 2] ;
1010 908
1011 //c = x1/x0; << 909 c = x1/x0;
1012 d = e0/x0; 910 d = e0/x0;
1013 a = log10(yy1/y0)/log10(x1/x0); << 911 a = log10(yy1/y0)/log10(x1/x0) ;
1014 << 912 // b0 = log10(y0) - a*log10(x0) ;
1015 b = 0.0; << 913 b = y0/pow(x0,a) ;
1016 if(a < 20.) b = y0/pow(x0,a); << 914 a += 1 ;
1017 <<
1018 a += 1; <<
1019 if(a == 0) 915 if(a == 0)
1020 { 916 {
1021 result += b*log(e0/x0); << 917 result += b*log(e0/x0) ;
1022 } 918 }
1023 else 919 else
1024 { 920 {
1025 result += y0*(e0*pow(d,a-1) - x0)/a; << 921 result += y0*(e0*pow(d,a-1) - x0)/a ;
1026 } 922 }
1027 a++; << 923 a++ ;
1028 if(a == 0) 924 if(a == 0)
1029 { 925 {
1030 fIntegralPAIySection[0] += b*log(e0/x0) << 926 fIntegralPAIySection[0] += b*log(e0/x0) ;
1031 } 927 }
1032 else 928 else
1033 { 929 {
1034 fIntegralPAIySection[0] += y0*(e0*e0*po << 930 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1035 } 931 }
1036 return result; << 932 return result ;
1037 933
1038 } 934 }
1039 935
1040 ///////////////////////////////////////////// 936 ///////////////////////////////////////////////////////////////////////
1041 937
1042 G4double G4PAIySection::SumOverBorderdEdx( G4 938 G4double G4PAIySection::SumOverBorderdEdx( G4int i ,
1043 G4doub 939 G4double en0 )
1044 { 940 {
1045 G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result; << 941 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
1046 942
1047 e0 = en0; << 943 e0 = en0 ;
1048 x0 = fSplineEnergy[i]; << 944 x0 = fSplineEnergy[i] ;
1049 x1 = fSplineEnergy[i+1]; << 945 x1 = fSplineEnergy[i+1] ;
1050 y0 = fDifPAIySection[i]; << 946 y0 = fDifPAIySection[i] ;
1051 yy1 = fDifPAIySection[i+1]; << 947 yy1 = fDifPAIySection[i+1] ;
1052 948
>> 949 c = x1/x0;
1053 d = e0/x0; 950 d = e0/x0;
1054 a = log10(yy1/y0)/log10(x1/x0); << 951 a = log10(yy1/y0)/log10(x1/x0) ;
1055 952
1056 G4double b = 0.0; 953 G4double b = 0.0;
1057 if(a < 20.) b = y0/pow(x0,a); << 954 if(a < 20.) b = y0/pow(x0,a) ;
1058 955
1059 a += 2; << 956 a += 2 ;
1060 if(a == 0) 957 if(a == 0)
1061 { 958 {
1062 result = b*log(x0/e0); << 959 result = b*log(x0/e0) ;
1063 } 960 }
1064 else 961 else
1065 { 962 {
1066 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/ << 963 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1067 } 964 }
1068 x0 = fSplineEnergy[i - 1]; << 965 x0 = fSplineEnergy[i - 1] ;
1069 x1 = fSplineEnergy[i - 2]; << 966 x1 = fSplineEnergy[i - 2] ;
1070 y0 = fDifPAIySection[i - 1]; << 967 y0 = fDifPAIySection[i - 1] ;
1071 yy1 = fDifPAIySection[i - 2]; << 968 yy1 = fDifPAIySection[i - 2] ;
1072 969
>> 970 c = x1/x0;
1073 d = e0/x0; 971 d = e0/x0;
1074 a = log10(yy1/y0)/log10(x1/x0); << 972 a = log10(yy1/y0)/log10(x1/x0) ;
1075 973
1076 b = 0.0; << 974 if(a < 20.) b = y0/pow(x0,a) ;
1077 if(a < 20.) b = y0/pow(x0,a); <<
1078 975
1079 a += 2; << 976 a += 2 ;
1080 if(a == 0) 977 if(a == 0)
1081 { 978 {
1082 result += b*log(e0/x0); << 979 result += b*log(e0/x0) ;
1083 } 980 }
1084 else 981 else
1085 { 982 {
1086 result += y0*(e0*e0*pow(d,a-2) - x0*x0) << 983 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1087 } 984 }
1088 return result; << 985 return result ;
>> 986
1089 } 987 }
1090 988
1091 ///////////////////////////////////////////// 989 ///////////////////////////////////////////////////////////////////////////////
1092 // 990 //
1093 // Integration of Cerenkov cross-section for 991 // Integration of Cerenkov cross-section for the case of
1094 // passing across border between intervals 992 // passing across border between intervals
1095 993
1096 G4double G4PAIySection::SumOverBordCerenkov( 994 G4double G4PAIySection::SumOverBordCerenkov( G4int i ,
1097 995 G4double en0 )
1098 { 996 {
1099 G4double x0,x1,y0,yy1,a,e0,c,d,result; << 997 G4double x0,x1,y0,yy1,a,e0,c,d,result ;
1100 998
1101 e0 = en0; << 999 e0 = en0 ;
1102 x0 = fSplineEnergy[i]; << 1000 x0 = fSplineEnergy[i] ;
1103 x1 = fSplineEnergy[i+1]; << 1001 x1 = fSplineEnergy[i+1] ;
1104 y0 = fdNdxCerenkov[i]; << 1002 y0 = fdNdxCerenkov[i] ;
1105 yy1 = fdNdxCerenkov[i+1]; << 1003 yy1 = fdNdxCerenkov[i+1] ;
1106 1004
1107 // G4cout<<G4endl; 1005 // G4cout<<G4endl;
1108 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "< 1006 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1109 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G 1007 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1110 c = x1/x0; << 1008 c = x1/x0 ;
1111 d = e0/x0; << 1009 d = e0/x0 ;
1112 a = log10(yy1/y0)/log10(c); << 1010 a = log10(yy1/y0)/log10(c) ;
1113 1011
1114 G4double b = 0.0; 1012 G4double b = 0.0;
1115 if(a < 20.) b = y0/pow(x0,a); << 1013 if(a < 20.) b = y0/pow(x0,a) ;
1116 1014
1117 a += 1.0; << 1015 a += 1.0 ;
1118 if( a == 0 ) result = b*log(x0/e0); << 1016 if( a == 0 ) result = b*log(x0/e0) ;
1119 else result = y0*(x0 - e0*pow(d,a- << 1017 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1120 a += 1.0; << 1018 a += 1.0 ;
1121 1019
1122 if( a == 0 ) fIntegralCerenkov[0] += b*log << 1020 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0) ;
1123 else fIntegralCerenkov[0] += y0*(x << 1021 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1124 1022
1125 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul 1023 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1126 1024
1127 x0 = fSplineEnergy[i - 1]; << 1025 x0 = fSplineEnergy[i - 1] ;
1128 x1 = fSplineEnergy[i - 2]; << 1026 x1 = fSplineEnergy[i - 2] ;
1129 y0 = fdNdxCerenkov[i - 1]; << 1027 y0 = fdNdxCerenkov[i - 1] ;
1130 yy1 = fdNdxCerenkov[i - 2]; << 1028 yy1 = fdNdxCerenkov[i - 2] ;
1131 1029
1132 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 1030 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1133 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 1031 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1134 1032
1135 c = x1/x0; << 1033 c = x1/x0 ;
1136 d = e0/x0; << 1034 d = e0/x0 ;
1137 a = log10(yy1/y0)/log10(x1/x0); << 1035 a = log10(yy1/y0)/log10(x1/x0) ;
1138 1036
1139 // G4cout << "a= " << a << G4endl; 1037 // G4cout << "a= " << a << G4endl;
>> 1038 if(a < 20.) b = y0/pow(x0,a) ;
>> 1039
1140 if(a > 20.0) b = 0.0; 1040 if(a > 20.0) b = 0.0;
1141 else b = y0/pow(x0,a); << 1041 else b = y0/pow(x0,a); // pow(10.,b0) ;
1142 1042
1143 //G4cout << "b= " << b << G4endl; 1043 //G4cout << "b= " << b << G4endl;
1144 1044
1145 a += 1.0; << 1045 a += 1.0 ;
1146 if( a == 0 ) result += b*log(e0/x0); << 1046 if( a == 0 ) result += b*log(e0/x0) ;
1147 else result += y0*(e0*pow(d,a-1) - << 1047 else result += y0*(e0*pow(d,a-1) - x0 )/a ;
1148 a += 1.0; << 1048 a += 1.0 ;
1149 //G4cout << "result= " << result << G4endl 1049 //G4cout << "result= " << result << G4endl;
1150 1050
1151 if( a == 0 ) fIntegralCerenkov[0] += b*l << 1051 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0) ;
1152 else fIntegralCerenkov[0] += y0* << 1052 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1153 1053
1154 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul 1054 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1155 1055
1156 return result; << 1056 return result ;
>> 1057
1157 } 1058 }
1158 1059
1159 ///////////////////////////////////////////// 1060 ///////////////////////////////////////////////////////////////////////////////
1160 // 1061 //
1161 // Integration of Plasmon cross-section for t 1062 // Integration of Plasmon cross-section for the case of
1162 // passing across border between intervals 1063 // passing across border between intervals
1163 1064
1164 G4double G4PAIySection::SumOverBordPlasmon( G 1065 G4double G4PAIySection::SumOverBordPlasmon( G4int i ,
1165 1066 G4double en0 )
1166 { 1067 {
1167 G4double x0,x1,y0,yy1,a,c,d,e0,result; << 1068 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
1168 1069
1169 e0 = en0; << 1070 e0 = en0 ;
1170 x0 = fSplineEnergy[i]; << 1071 x0 = fSplineEnergy[i] ;
1171 x1 = fSplineEnergy[i+1]; << 1072 x1 = fSplineEnergy[i+1] ;
1172 y0 = fdNdxPlasmon[i]; << 1073 y0 = fdNdxPlasmon[i] ;
1173 yy1 = fdNdxPlasmon[i+1]; << 1074 yy1 = fdNdxPlasmon[i+1] ;
1174 << 1075
1175 c = x1/x0; << 1076 c = x1/x0 ;
1176 d = e0/x0; << 1077 d = e0/x0 ;
1177 a = log10(yy1/y0)/log10(c); << 1078 a = log10(yy1/y0)/log10(c) ;
1178 1079
1179 G4double b = 0.0; 1080 G4double b = 0.0;
1180 if(a < 20.) b = y0/pow(x0,a); << 1081 if(a < 20.) b = y0/pow(x0,a) ;
1181 1082
1182 a += 1.0; << 1083 a += 1.0 ;
1183 if( a == 0 ) result = b*log(x0/e0); << 1084 if( a == 0 ) result = b*log(x0/e0) ;
1184 else result = y0*(x0 - e0*pow(d,a- << 1085 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1185 a += 1.0; << 1086 a += 1.0 ;
1186 1087
1187 if( a == 0 ) fIntegralPlasmon[0] += b*log( << 1088 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0) ;
1188 else fIntegralPlasmon[0] += y0*(x0 << 1089 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1189 1090
1190 x0 = fSplineEnergy[i - 1]; << 1091 x0 = fSplineEnergy[i - 1] ;
1191 x1 = fSplineEnergy[i - 2]; << 1092 x1 = fSplineEnergy[i - 2] ;
1192 y0 = fdNdxPlasmon[i - 1]; << 1093 y0 = fdNdxPlasmon[i - 1] ;
1193 yy1 = fdNdxPlasmon[i - 2]; << 1094 yy1 = fdNdxPlasmon[i - 2] ;
1194 << 1095
1195 c = x1/x0; << 1096 c = x1/x0 ;
1196 d = e0/x0; << 1097 d = e0/x0 ;
1197 a = log10(yy1/y0)/log10(c); << 1098 a = log10(yy1/y0)/log10(c) ;
1198 1099
1199 if(a < 20.) b = y0/pow(x0,a); << 1100 if(a < 20.) b = y0/pow(x0,a) ;
1200 1101
1201 a += 1.0; << 1102 a += 1.0 ;
1202 if( a == 0 ) result += b*log(e0/x0); << 1103 if( a == 0 ) result += b*log(e0/x0) ;
1203 else result += y0*(e0*pow(d,a-1) - << 1104 else result += y0*(e0*pow(d,a-1) - x0)/a ;
1204 a += 1.0; << 1105 a += 1.0 ;
1205 1106
1206 if( a == 0 ) fIntegralPlasmon[0] += b*lo << 1107 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0) ;
1207 else fIntegralPlasmon[0] += y0*( << 1108 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1208 1109
1209 return result; << 1110 return result ;
1210 1111
1211 } 1112 }
1212 1113
1213 ///////////////////////////////////////////// 1114 /////////////////////////////////////////////////////////////////////////
1214 // 1115 //
1215 // 1116 //
1216 1117
1217 G4double G4PAIySection::GetStepEnergyLoss( G4 1118 G4double G4PAIySection::GetStepEnergyLoss( G4double step )
1218 { 1119 {
1219 G4int iTransfer ; << 1120 G4int iTransfer ;
1220 G4long numOfCollisions; << 1121 G4long numOfCollisions ;
1221 G4double loss = 0.0; << 1122 G4double loss = 0.0 ;
1222 G4double meanNumber, position; << 1123 G4double meanNumber, position ;
1223 1124
1224 // G4cout<<" G4PAIySection::GetStepEnergyLo << 1125 // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl ;
1225 1126
1226 1127
1227 1128
1228 meanNumber = fIntegralPAIySection[1]*step; << 1129 meanNumber = fIntegralPAIySection[1]*step ;
1229 numOfCollisions = G4Poisson(meanNumber); << 1130 numOfCollisions = G4Poisson(meanNumber) ;
1230 1131
1231 // G4cout<<"numOfCollisions = "<<numOfCol << 1132 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1232 1133
1233 while(numOfCollisions) 1134 while(numOfCollisions)
1234 { 1135 {
1235 position = fIntegralPAIySection[1]*G4Unif << 1136 position = fIntegralPAIySection[1]*G4UniformRand() ;
1236 1137
1237 for( iTransfer=1; iTransfer<=fSplineNumbe << 1138 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1238 { 1139 {
1239 if( position >= fIntegralPAIySection[ << 1140 if( position >= fIntegralPAIySection[iTransfer] ) break ;
1240 } 1141 }
1241 loss += fSplineEnergy[iTransfer] ; << 1142 loss += fSplineEnergy[iTransfer] ;
1242 numOfCollisions--; << 1143 numOfCollisions-- ;
1243 // Loop checking, 03-Aug-2015, Vladimir I <<
1244 } 1144 }
1245 // G4cout<<"PAI energy loss = "<<loss/keV<< << 1145 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl ;
1246 1146
1247 return loss; << 1147 return loss ;
1248 } 1148 }
1249 1149
1250 ///////////////////////////////////////////// 1150 /////////////////////////////////////////////////////////////////////////
1251 // 1151 //
1252 // 1152 //
1253 1153
1254 G4double G4PAIySection::GetStepCerenkovLoss( 1154 G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
1255 { 1155 {
1256 G4int iTransfer ; << 1156 G4int iTransfer ;
1257 G4long numOfCollisions; << 1157 G4long numOfCollisions ;
1258 G4double loss = 0.0; << 1158 G4double loss = 0.0 ;
1259 G4double meanNumber, position; << 1159 G4double meanNumber, position ;
1260 1160
1261 // G4cout<<" G4PAIySection::GetStepCreLosnk << 1161 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ;
1262 1162
1263 1163
1264 1164
1265 meanNumber = fIntegralCerenkov[1]*step; << 1165 meanNumber = fIntegralCerenkov[1]*step ;
1266 numOfCollisions = G4Poisson(meanNumber); << 1166 numOfCollisions = G4Poisson(meanNumber) ;
1267 1167
1268 // G4cout<<"numOfCollisions = "<<numOfCol << 1168 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1269 1169
1270 while(numOfCollisions) 1170 while(numOfCollisions)
1271 { 1171 {
1272 position = fIntegralCerenkov[1]*G4Uniform << 1172 position = fIntegralCerenkov[1]*G4UniformRand() ;
1273 1173
1274 for( iTransfer=1; iTransfer<=fSplineNumbe << 1174 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1275 { 1175 {
1276 if( position >= fIntegralCerenkov[iTr << 1176 if( position >= fIntegralCerenkov[iTransfer] ) break ;
1277 } 1177 }
1278 loss += fSplineEnergy[iTransfer] ; << 1178 loss += fSplineEnergy[iTransfer] ;
1279 numOfCollisions--; << 1179 numOfCollisions-- ;
1280 // Loop checking, 03-Aug-2015, Vladimir I <<
1281 } 1180 }
1282 // G4cout<<"PAI Cerenkov loss = "<<loss/keV << 1181 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl ;
1283 1182
1284 return loss; << 1183 return loss ;
1285 } 1184 }
1286 1185
1287 ///////////////////////////////////////////// 1186 /////////////////////////////////////////////////////////////////////////
1288 // 1187 //
1289 // 1188 //
1290 1189
1291 G4double G4PAIySection::GetStepPlasmonLoss( G 1190 G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
1292 { 1191 {
1293 G4int iTransfer ; << 1192 G4int iTransfer ;
1294 G4long numOfCollisions; << 1193 G4long numOfCollisions ;
1295 G4double loss = 0.0; << 1194 G4double loss = 0.0 ;
1296 G4double meanNumber, position; << 1195 G4double meanNumber, position ;
1297 1196
1298 // G4cout<<" G4PAIySection::GetStepCreLosnk << 1197 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ;
1299 1198
1300 1199
1301 1200
1302 meanNumber = fIntegralPlasmon[1]*step; << 1201 meanNumber = fIntegralPlasmon[1]*step ;
1303 numOfCollisions = G4Poisson(meanNumber); << 1202 numOfCollisions = G4Poisson(meanNumber) ;
1304 1203
1305 // G4cout<<"numOfCollisions = "<<numOfCol << 1204 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1306 1205
1307 while(numOfCollisions) 1206 while(numOfCollisions)
1308 { 1207 {
1309 position = fIntegralPlasmon[1]*G4UniformR << 1208 position = fIntegralPlasmon[1]*G4UniformRand() ;
1310 1209
1311 for( iTransfer=1; iTransfer<=fSplineNumbe << 1210 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1312 { 1211 {
1313 if( position >= fIntegralPlasmon[iTra << 1212 if( position >= fIntegralPlasmon[iTransfer] ) break ;
1314 } 1213 }
1315 loss += fSplineEnergy[iTransfer] ; << 1214 loss += fSplineEnergy[iTransfer] ;
1316 numOfCollisions--; << 1215 numOfCollisions-- ;
1317 // Loop checking, 03-Aug-2015, Vladimir I <<
1318 } 1216 }
1319 // G4cout<<"PAI Plasmon loss = "<<loss/keV< << 1217 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl ;
1320 1218
1321 return loss; << 1219 return loss ;
1322 } 1220 }
1323 1221
1324 ///////////////////////////////////////////// <<
1325 // <<
1326 1222
1327 void G4PAIySection::CallError(G4int i, const <<
1328 { <<
1329 G4String head = "G4PAIySection::" + methodN <<
1330 G4ExceptionDescription ed; <<
1331 ed << "Wrong index " << i << " fSplineNumbe <<
1332 G4Exception(head,"pai001",FatalException,ed <<
1333 } <<
1334 1223
1335 ///////////////////////////////////////////// 1224 /////////////////////////////////////////////////////////////////////////////
1336 // 1225 //
1337 // Init array of Lorentz factors 1226 // Init array of Lorentz factors
1338 // 1227 //
1339 1228
1340 G4int G4PAIySection::fNumberOfGammas = 111; << 1229 G4int G4PAIySection::fNumberOfGammas = 111 ;
1341 1230
1342 const G4double G4PAIySection::fLorentzFactor[ 1231 const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1
1343 { 1232 {
1344 0.0, 1233 0.0,
1345 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.1 1234 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
1346 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.2 1235 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
1347 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.4 1236 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
1348 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.9 1237 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
1349 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.7 1238 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
1350 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.2 1239 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
1351 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.2 1240 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
1352 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.2 1241 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
1353 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.3 1242 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
1354 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.2 1243 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
1355 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.9 1244 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
1356 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.4 1245 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
1357 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.7 1246 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
1358 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.2 1247 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
1359 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.8 1248 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
1360 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.8 1249 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
1361 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.4 1250 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
1362 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.5 1251 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
1363 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.2 1252 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
1364 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.3 1253 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
1365 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.3 1254 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
1366 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.2 1255 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
1367 1.065799e+05 1256 1.065799e+05
1368 }; << 1257 } ;
1369 1258
1370 ///////////////////////////////////////////// 1259 ///////////////////////////////////////////////////////////////////////
1371 // 1260 //
1372 // The number of gamma for creation of splin 1261 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
1373 // 1262 //
1374 1263
1375 const G4int G4PAIySection::fRefGammaNumber = << 1264 const
>> 1265 G4int G4PAIySection::fRefGammaNumber = 29 ;
1376 1266
>> 1267
1377 // 1268 //
1378 // end of G4PAIySection implementation file 1269 // end of G4PAIySection implementation file
1379 // 1270 //
1380 ///////////////////////////////////////////// 1271 ////////////////////////////////////////////////////////////////////////////
1381 1272
1382 1273