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