| 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 // 27 //
28 // 28 //
29 // G4PAIxSection.cc -- class implementation fi 29 // G4PAIxSection.cc -- class implementation file
30 // 30 //
31 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
32 // 32 //
33 // For information related to this code, pleas 33 // For information related to this code, please, contact
34 // the Geant4 Collaboration. 34 // the Geant4 Collaboration.
35 // 35 //
36 // R&D: Vladimir.Grichine@cern.ch 36 // R&D: Vladimir.Grichine@cern.ch
37 // 37 //
38 // History: 38 // History:
39 // 39 //
40 // 13.05.03 V. Grichine, bug fixed for maxEner 40 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
41 // 28.05.01 V.Ivanchenko minor changes to prov 41 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
42 // 17.05.01 V. Grichine, low energy extension 42 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
43 // 20.11.98 adapted to a new Material/SandiaTa 43 // 20.11.98 adapted to a new Material/SandiaTable interface, mma
44 // 11.06.97 V. Grichine, 1st version 44 // 11.06.97 V. Grichine, 1st version
45 // 45 //
46 46
47 #include "G4PAIxSection.hh" 47 #include "G4PAIxSection.hh"
48 48
49 #include "globals.hh" 49 #include "globals.hh"
50 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
52 #include "G4ios.hh" 52 #include "G4ios.hh"
53 #include "G4Poisson.hh" 53 #include "G4Poisson.hh"
54 #include "G4Material.hh" 54 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh" 55 #include "G4MaterialCutsCouple.hh"
56 #include "G4SandiaTable.hh" 56 #include "G4SandiaTable.hh"
57 57
58 using namespace std; 58 using namespace std;
59 59
60 /* ******************************************* 60 /* ******************************************************************
61 61
62 // Init array of Lorentz factors 62 // Init array of Lorentz factors
63 63
64 const G4double G4PAIxSection::fLorentzFactor[2 64 const G4double G4PAIxSection::fLorentzFactor[22] =
65 { 65 {
66 0.0 , 1.1 , 1.2 , 1.3 , 1 66 0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
67 2.5 , 3.0 , 4.0 , 7.0 , 10 67 2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
68 70.0 , 100.0 , 300.0 , 600.0 , 1000 68 70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
69 10000.0 , 50000.0 69 10000.0 , 50000.0
70 }; 70 };
71 71
72 const G4int G4PAIxSection:: 72 const G4int G4PAIxSection::
73 fRefGammaNumber = 29; // The number of 73 fRefGammaNumber = 29; // The number of gamma for creation of
74 // spline (9) 74 // spline (9)
75 75
76 ********************************************** 76 ***************************************************************** */
77 77
78 // Local class constants 78 // Local class constants
79 79
80 const G4double G4PAIxSection::fDelta = 0.005; 80 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
81 const G4double G4PAIxSection::fError = 0.005; 81 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
82 82
83 const G4int G4PAIxSection::fMaxSplineSize = 10 83 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
84 84 // arrays
85 ////////////////////////////////////////////// 85 //////////////////////////////////////////////////////////////////
86 // 86 //
87 // Constructor 87 // Constructor
88 // 88 //
89 89
90 G4PAIxSection::G4PAIxSection() 90 G4PAIxSection::G4PAIxSection()
91 { 91 {
92 fSandia = nullptr; 92 fSandia = nullptr;
93 fMatSandiaMatrix = nullptr; 93 fMatSandiaMatrix = nullptr;
94 fDensity = fElectronDensity = fNormalization 94 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
95 fIntervalNumber = fSplineNumber = 0; 95 fIntervalNumber = fSplineNumber = 0;
96 fVerbose = 0; 96 fVerbose = 0;
97 97
98 fSplineEnergy = G4DataVector(fMaxSp 98 fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
99 fRePartDielectricConst = G4DataVector(fMaxSp 99 fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
100 fImPartDielectricConst = G4DataVector(fMaxSp 100 fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
101 fIntegralTerm = G4DataVector(fMaxSp 101 fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
102 fDifPAIxSection = G4DataVector(fMaxSp 102 fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
103 fdNdxCerenkov = G4DataVector(fMaxSp 103 fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
104 fdNdxPlasmon = G4DataVector(fMaxSp 104 fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
105 fdNdxMM = G4DataVector(fMaxSp 105 fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
106 fdNdxResonance = G4DataVector(fMaxSp 106 fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
107 fIntegralPAIxSection = G4DataVector(fMaxSp 107 fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
108 fIntegralPAIdEdx = G4DataVector(fMaxSp 108 fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
109 fIntegralCerenkov = G4DataVector(fMaxSp 109 fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
110 fIntegralPlasmon = G4DataVector(fMaxSp 110 fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
111 fIntegralMM = G4DataVector(fMaxSp 111 fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
112 fIntegralResonance = G4DataVector(fMaxSp 112 fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
113 113
114 fMaterialIndex = 0; 114 fMaterialIndex = 0;
115 115
116 for( G4int i = 0; i < 500; ++i ) 116 for( G4int i = 0; i < 500; ++i )
117 { 117 {
118 for( G4int j = 0; j < 112; ++j ) fPAItabl 118 for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
119 } 119 }
120 } 120 }
121 121
122 ////////////////////////////////////////////// 122 //////////////////////////////////////////////////////////////////
123 // 123 //
124 // Constructor 124 // Constructor
125 // 125 //
126 126
127 G4PAIxSection::G4PAIxSection(G4MaterialCutsCou 127 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
128 { 128 {
129 fDensity = matCC->GetMaterial()->GetDe 129 fDensity = matCC->GetMaterial()->GetDensity();
130 G4int matIndex = (G4int)matCC->GetMaterial() 130 G4int matIndex = (G4int)matCC->GetMaterial()->GetIndex();
131 fMaterialIndex = matIndex; 131 fMaterialIndex = matIndex;
132 132
133 const G4MaterialTable* theMaterialTable = G4 133 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
134 fSandia = (*theMaterialTable)[matIndex]->Get 134 fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
135 135
136 fVerbose = 0; 136 fVerbose = 0;
137 137
138 G4int i, j; 138 G4int i, j;
139 fMatSandiaMatrix = new G4OrderedTable(); 139 fMatSandiaMatrix = new G4OrderedTable();
140 140
141 for (i = 0; i < fSandia->GetMaxInterval()-1; 141 for (i = 0; i < fSandia->GetMaxInterval()-1; ++i)
142 { 142 {
143 fMatSandiaMatrix->push_back(new G4DataVec 143 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
144 } 144 }
145 for (i = 0; i < fSandia->GetMaxInterval()-1; 145 for (i = 0; i < fSandia->GetMaxInterval()-1; ++i)
146 { 146 {
147 (*(*fMatSandiaMatrix)[i])[0] = fSandia->Ge 147 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
148 148
149 for(j = 1; j < 5; ++j) 149 for(j = 1; j < 5; ++j)
150 { 150 {
151 (*(*fMatSandiaMatrix)[i])[j] = fSandia-> 151 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
152 } 152 }
153 } 153 }
154 ComputeLowEnergyCof(); 154 ComputeLowEnergyCof();
155 // fEnergyInterval = fA1 = fA2 = fA3 = fA4 155 // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
156 } 156 }
157 157
158 ////////////////////////////////////////////// 158 ////////////////////////////////////////////////////////////////
159 159
160 G4PAIxSection::G4PAIxSection(G4int materialInd 160 G4PAIxSection::G4PAIxSection(G4int materialIndex,
161 G4double maxEnerg 161 G4double maxEnergyTransfer)
162 { 162 {
163 fSandia = nullptr; 163 fSandia = nullptr;
164 fMatSandiaMatrix = nullptr; 164 fMatSandiaMatrix = nullptr;
165 fVerbose = 0; 165 fVerbose = 0;
166 const G4MaterialTable* theMaterialTable = G4 166 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
167 G4int i, j; 167 G4int i, j;
168 168
169 fMaterialIndex = materialIndex; 169 fMaterialIndex = materialIndex;
170 fDensity = (*theMaterialTable 170 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
171 fElectronDensity = (*theMaterialTable 171 fElectronDensity = (*theMaterialTable)[materialIndex]->
172 GetElectronDensit 172 GetElectronDensity();
173 fIntervalNumber = (*theMaterialTable 173 fIntervalNumber = (*theMaterialTable)[materialIndex]->
174 GetSandiaTable()- 174 GetSandiaTable()->GetMatNbOfIntervals();
175 fIntervalNumber--; 175 fIntervalNumber--;
176 // G4cout<<fDensity<<"\t"<<fElectronDensity< 176 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
177 177
178 fEnergyInterval = G4DataVector(fIntervalNumb 178 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
179 fA1 = G4DataVector(fIntervalNumb 179 fA1 = G4DataVector(fIntervalNumber+2,0.0);
180 fA2 = G4DataVector(fIntervalNumb 180 fA2 = G4DataVector(fIntervalNumber+2,0.0);
181 fA3 = G4DataVector(fIntervalNumb 181 fA3 = G4DataVector(fIntervalNumber+2,0.0);
182 fA4 = G4DataVector(fIntervalNumb 182 fA4 = G4DataVector(fIntervalNumber+2,0.0);
183 183
184 for(i = 1; i <= fIntervalNumber; i++ ) 184 for(i = 1; i <= fIntervalNumber; i++ )
185 { 185 {
186 if(((*theMaterialTable)[materialIndex]-> 186 if(((*theMaterialTable)[materialIndex]->
187 GetSandiaTable()->GetSandiaCofForMaterial( 187 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
188 i > fIntervalNumber 188 i > fIntervalNumber )
189 { 189 {
190 fEnergyInterval[i] = maxEnergyTransf 190 fEnergyInterval[i] = maxEnergyTransfer;
191 fIntervalNumber = i; 191 fIntervalNumber = i;
192 break; 192 break;
193 } 193 }
194 fEnergyInterval[i] = (*theMaterialTab 194 fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
195 GetSandiaTable() 195 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
196 fA1[i] = (*theMaterialTab 196 fA1[i] = (*theMaterialTable)[materialIndex]->
197 GetSandiaTable() 197 GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
198 fA2[i] = (*theMaterialTab 198 fA2[i] = (*theMaterialTable)[materialIndex]->
199 GetSandiaTable() 199 GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
200 fA3[i] = (*theMaterialTab 200 fA3[i] = (*theMaterialTable)[materialIndex]->
201 GetSandiaTable() 201 GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
202 fA4[i] = (*theMaterialTab 202 fA4[i] = (*theMaterialTable)[materialIndex]->
203 GetSandiaTable() 203 GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
204 // G4cout<<i<<"\t"<<fEnergyInterval[i 204 // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
205 // <<fA 205 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
206 } 206 }
207 if(fEnergyInterval[fIntervalNumber] != maxEn 207 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
208 { 208 {
209 fIntervalNumber++; 209 fIntervalNumber++;
210 fEnergyInterval[fIntervalNumber] = ma 210 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
211 } 211 }
212 212
213 // Now checking, if two borders are too clos 213 // Now checking, if two borders are too close together
214 214
215 for(i=1;i<fIntervalNumber;i++) 215 for(i=1;i<fIntervalNumber;i++)
216 { 216 {
217 if(fEnergyInterval[i+1]-fEnergyInterva 217 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
218 1.5*fDelta*(fEnergyInterval[i+1]+fE 218 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
219 { 219 {
220 continue; 220 continue;
221 } 221 }
222 else 222 else
223 { 223 {
224 for(j=i;j<fIntervalNumber;j++) 224 for(j=i;j<fIntervalNumber;j++)
225 { 225 {
226 fEnergyInterval[j] = fEnergyInterv 226 fEnergyInterval[j] = fEnergyInterval[j+1];
227 fA1[j] = fA1[j+1]; 227 fA1[j] = fA1[j+1];
228 fA2[j] = fA2[j+1]; 228 fA2[j] = fA2[j+1];
229 fA3[j] = fA3[j+1]; 229 fA3[j] = fA3[j+1];
230 fA4[j] = fA4[j+1]; 230 fA4[j] = fA4[j+1];
231 } 231 }
232 fIntervalNumber--; 232 fIntervalNumber--;
233 i--; 233 i--;
234 } 234 }
235 } 235 }
236 236
237 237
238 /* ********************************* 238 /* *********************************
239 239
240 fSplineEnergy = new G4double[fM 240 fSplineEnergy = new G4double[fMaxSplineSize];
241 fRePartDielectricConst = new G4double[fM 241 fRePartDielectricConst = new G4double[fMaxSplineSize];
242 fImPartDielectricConst = new G4double[fM 242 fImPartDielectricConst = new G4double[fMaxSplineSize];
243 fIntegralTerm = new G4double[fM 243 fIntegralTerm = new G4double[fMaxSplineSize];
244 fDifPAIxSection = new G4double[fM 244 fDifPAIxSection = new G4double[fMaxSplineSize];
245 fIntegralPAIxSection = new G4double[fM 245 fIntegralPAIxSection = new G4double[fMaxSplineSize];
246 246
247 for(i=0;i<fMaxSplineSize;i++) 247 for(i=0;i<fMaxSplineSize;i++)
248 { 248 {
249 fSplineEnergy[i] = 0.0; 249 fSplineEnergy[i] = 0.0;
250 fRePartDielectricConst[i] = 0.0; 250 fRePartDielectricConst[i] = 0.0;
251 fImPartDielectricConst[i] = 0.0; 251 fImPartDielectricConst[i] = 0.0;
252 fIntegralTerm[i] = 0.0; 252 fIntegralTerm[i] = 0.0;
253 fDifPAIxSection[i] = 0.0; 253 fDifPAIxSection[i] = 0.0;
254 fIntegralPAIxSection[i] = 0.0; 254 fIntegralPAIxSection[i] = 0.0;
255 } 255 }
256 **************************************** 256 ************************************************** */
257 ComputeLowEnergyCof(); 257 ComputeLowEnergyCof();
258 InitPAI(); // create arrays allocated a 258 InitPAI(); // create arrays allocated above
259 /* 259 /*
260 delete[] fEnergyInterval; 260 delete[] fEnergyInterval;
261 delete[] fA1; 261 delete[] fA1;
262 delete[] fA2; 262 delete[] fA2;
263 delete[] fA3; 263 delete[] fA3;
264 delete[] fA4; 264 delete[] fA4;
265 */ 265 */
266 } 266 }
267 267
268 ////////////////////////////////////////////// 268 ////////////////////////////////////////////////////////////////////////
269 // 269 //
270 // Constructor called from G4PAIPhotonModel !! 270 // Constructor called from G4PAIPhotonModel !!!
271 271
272 G4PAIxSection::G4PAIxSection( G4int materialIn 272 G4PAIxSection::G4PAIxSection( G4int materialIndex,
273 G4double maxEner 273 G4double maxEnergyTransfer,
274 G4double betaGam 274 G4double betaGammaSq,
275 G4double** photo 275 G4double** photoAbsCof,
276 G4int intNumber 276 G4int intNumber )
277 { 277 {
278 fSandia = nullptr; 278 fSandia = nullptr;
279 fDensity = fElectronDensity = fNormalization 279 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
280 fIntervalNumber = fSplineNumber = 0; 280 fIntervalNumber = fSplineNumber = 0;
281 fVerbose = 0; 281 fVerbose = 0;
282 282
283 fSplineEnergy = G4DataVector(500,0. 283 fSplineEnergy = G4DataVector(500,0.0);
284 fRePartDielectricConst = G4DataVector(500,0. 284 fRePartDielectricConst = G4DataVector(500,0.0);
285 fImPartDielectricConst = G4DataVector(500,0. 285 fImPartDielectricConst = G4DataVector(500,0.0);
286 fIntegralTerm = G4DataVector(500,0. 286 fIntegralTerm = G4DataVector(500,0.0);
287 fDifPAIxSection = G4DataVector(500,0. 287 fDifPAIxSection = G4DataVector(500,0.0);
288 fdNdxCerenkov = G4DataVector(500,0. 288 fdNdxCerenkov = G4DataVector(500,0.0);
289 fdNdxPlasmon = G4DataVector(500,0. 289 fdNdxPlasmon = G4DataVector(500,0.0);
290 fdNdxMM = G4DataVector(500,0. 290 fdNdxMM = G4DataVector(500,0.0);
291 fdNdxResonance = G4DataVector(500,0. 291 fdNdxResonance = G4DataVector(500,0.0);
292 fIntegralPAIxSection = G4DataVector(500,0. 292 fIntegralPAIxSection = G4DataVector(500,0.0);
293 fIntegralPAIdEdx = G4DataVector(500,0. 293 fIntegralPAIdEdx = G4DataVector(500,0.0);
294 fIntegralCerenkov = G4DataVector(500,0. 294 fIntegralCerenkov = G4DataVector(500,0.0);
295 fIntegralPlasmon = G4DataVector(500,0. 295 fIntegralPlasmon = G4DataVector(500,0.0);
296 fIntegralMM = G4DataVector(500,0. 296 fIntegralMM = G4DataVector(500,0.0);
297 fIntegralResonance = G4DataVector(500,0. 297 fIntegralResonance = G4DataVector(500,0.0);
298 298
299 for( G4int i = 0; i < 500; ++i ) 299 for( G4int i = 0; i < 500; ++i )
300 { 300 {
301 for( G4int j = 0; j < 112; ++j ) fPAItabl 301 for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
302 } 302 }
303 303
304 fSandia = nullptr; 304 fSandia = nullptr;
305 fMatSandiaMatrix = nullptr; 305 fMatSandiaMatrix = nullptr;
306 const G4MaterialTable* theMaterialTable = G4 306 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
307 G4int i, j; 307 G4int i, j;
308 308
309 fMaterialIndex = materialIndex; 309 fMaterialIndex = materialIndex;
310 fDensity = (*theMaterialTable)[mater 310 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
311 fElectronDensity = (*theMaterialTable)[mater 311 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
312 312
313 fIntervalNumber = intNumber; 313 fIntervalNumber = intNumber;
314 fIntervalNumber--; 314 fIntervalNumber--;
315 // G4cout<<fDensity<<"\t"<<fElectronDensit 315 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
316 316
317 fEnergyInterval = G4DataVector(fIntervalNumb 317 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
318 fA1 = G4DataVector(fIntervalNumb 318 fA1 = G4DataVector(fIntervalNumber+2,0.0);
319 fA2 = G4DataVector(fIntervalNumb 319 fA2 = G4DataVector(fIntervalNumber+2,0.0);
320 fA3 = G4DataVector(fIntervalNumb 320 fA3 = G4DataVector(fIntervalNumber+2,0.0);
321 fA4 = G4DataVector(fIntervalNumb 321 fA4 = G4DataVector(fIntervalNumber+2,0.0);
322 322
323 323
324 /* 324 /*
325 fEnergyInterval = new G4double[fInterval 325 fEnergyInterval = new G4double[fIntervalNumber+2];
326 fA1 = new G4double[fInterval 326 fA1 = new G4double[fIntervalNumber+2];
327 fA2 = new G4double[fInterval 327 fA2 = new G4double[fIntervalNumber+2];
328 fA3 = new G4double[fInterval 328 fA3 = new G4double[fIntervalNumber+2];
329 fA4 = new G4double[fInterval 329 fA4 = new G4double[fIntervalNumber+2];
330 */ 330 */
331 for( i = 1; i <= fIntervalNumber; i++ ) 331 for( i = 1; i <= fIntervalNumber; i++ )
332 { 332 {
333 if( ( photoAbsCof[i-1][0] >= maxEnerg 333 if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
334 i > fIntervalNumber ) 334 i > fIntervalNumber )
335 { 335 {
336 fEnergyInterval[i] = maxEnergyTran 336 fEnergyInterval[i] = maxEnergyTransfer;
337 fIntervalNumber = i; 337 fIntervalNumber = i;
338 break; 338 break;
339 } 339 }
340 fEnergyInterval[i] = photoAbsCof[i-1] 340 fEnergyInterval[i] = photoAbsCof[i-1][0];
341 fA1[i] = photoAbsCof[i-1] 341 fA1[i] = photoAbsCof[i-1][1];
342 fA2[i] = photoAbsCof[i-1] 342 fA2[i] = photoAbsCof[i-1][2];
343 fA3[i] = photoAbsCof[i-1] 343 fA3[i] = photoAbsCof[i-1][3];
344 fA4[i] = photoAbsCof[i-1] 344 fA4[i] = photoAbsCof[i-1][4];
345 // G4cout<<i<<"\t"<<fEnergyInterval[i 345 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
346 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"< 346 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
347 } 347 }
348 // G4cout<<"i last = "<<i<<"; "<<"fInter 348 // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
349 349
350 if(fEnergyInterval[fIntervalNumber] != maxEn 350 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
351 { 351 {
352 fIntervalNumber++; 352 fIntervalNumber++;
353 fEnergyInterval[fIntervalNumber] = ma 353 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
354 } 354 }
355 // G4cout<<"after check of max energy tr 355 // G4cout<<"after check of max energy transfer"<<G4endl;
356 356
357 for( i = 1; i <= fIntervalNumber; i++ ) 357 for( i = 1; i <= fIntervalNumber; i++ )
358 { 358 {
359 // G4cout<<i<<"\t"<<fEnergyInterval[i] 359 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
360 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e 360 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
361 } 361 }
362 // Now checking, if two borders are too 362 // Now checking, if two borders are too close together
363 363
364 for( i = 1; i < fIntervalNumber; i++ ) 364 for( i = 1; i < fIntervalNumber; i++ )
365 { 365 {
366 if(fEnergyInterval[i+1]-fEnergyInterva 366 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
367 1.5*fDelta*(fEnergyInterval[i+1]+fE 367 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
368 { 368 {
369 continue; 369 continue;
370 } 370 }
371 else 371 else
372 { 372 {
373 for(j=i;j<fIntervalNumber;j++) 373 for(j=i;j<fIntervalNumber;j++)
374 { 374 {
375 fEnergyInterval[j] = fEnergyInterv 375 fEnergyInterval[j] = fEnergyInterval[j+1];
376 fA1[j] = fA1[j+1]; 376 fA1[j] = fA1[j+1];
377 fA2[j] = fA2[j+1]; 377 fA2[j] = fA2[j+1];
378 fA3[j] = fA3[j+1]; 378 fA3[j] = fA3[j+1];
379 fA4[j] = fA4[j+1]; 379 fA4[j] = fA4[j+1];
380 } 380 }
381 fIntervalNumber--; 381 fIntervalNumber--;
382 i--; 382 i--;
383 } 383 }
384 } 384 }
385 // G4cout<<"after check of close borders"<<G 385 // G4cout<<"after check of close borders"<<G4endl;
386 386
387 for( i = 1; i <= fIntervalNumber; i++ ) 387 for( i = 1; i <= fIntervalNumber; i++ )
388 { 388 {
389 // G4cout<<i<<"\t"<<fEnergyInterval[i] 389 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
390 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4en 390 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
391 } 391 }
392 392
393 // Preparation of fSplineEnergy array corres 393 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
394 394
395 ComputeLowEnergyCof(); 395 ComputeLowEnergyCof();
396 G4double betaGammaSqRef = 396 G4double betaGammaSqRef =
397 fLorentzFactor[fRefGammaNumber]*fLorentzFa 397 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
398 398
399 NormShift(betaGammaSqRef); 399 NormShift(betaGammaSqRef);
400 SplainPAI(betaGammaSqRef); 400 SplainPAI(betaGammaSqRef);
401 401
402 // Preparation of integral PAI cross section 402 // Preparation of integral PAI cross section for input betaGammaSq
403 403
404 for(i = 1; i <= fSplineNumber; i++) 404 for(i = 1; i <= fSplineNumber; i++)
405 { 405 {
406 fdNdxCerenkov[i] = PAIdNdxCerenkov( 406 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
407 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 407 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
408 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 408 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
409 fdNdxResonance[i] = PAIdNdxResonance 409 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
410 fDifPAIxSection[i] = DifPAIxSection(i 410 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
411 411
412 // G4cout<<i<<"; dNdxC = "<<fdNdxCere 412 // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
413 // <<"; dNdxPAI = "<<fDifPAIxSecti 413 // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
414 } 414 }
415 IntegralCerenkov(); 415 IntegralCerenkov();
416 IntegralMM(); 416 IntegralMM();
417 IntegralPlasmon(); 417 IntegralPlasmon();
418 IntegralResonance(); 418 IntegralResonance();
419 IntegralPAIxSection(); 419 IntegralPAIxSection();
420 /* 420 /*
421 delete[] fEnergyInterval; 421 delete[] fEnergyInterval;
422 delete[] fA1; 422 delete[] fA1;
423 delete[] fA2; 423 delete[] fA2;
424 delete[] fA3; 424 delete[] fA3;
425 delete[] fA4; 425 delete[] fA4;
426 */ 426 */
427 } 427 }
428 428
429 ////////////////////////////////////////////// 429 ////////////////////////////////////////////////////////////////////////
430 // 430 //
431 // Test Constructor with beta*gamma square val 431 // Test Constructor with beta*gamma square value
432 432
433 G4PAIxSection::G4PAIxSection( G4int materialIn 433 G4PAIxSection::G4PAIxSection( G4int materialIndex,
434 G4double maxEner 434 G4double maxEnergyTransfer,
435 G4double betaGam 435 G4double betaGammaSq )
436 { 436 {
437 fSandia = nullptr; 437 fSandia = nullptr;
438 fMatSandiaMatrix = nullptr; 438 fMatSandiaMatrix = nullptr;
439 fVerbose = 0; 439 fVerbose = 0;
440 const G4MaterialTable* theMaterialTable = G4 440 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
441 441
442 G4int i, j, numberOfElements; 442 G4int i, j, numberOfElements;
443 443
444 fMaterialIndex = materialIndex; 444 fMaterialIndex = materialIndex;
445 fDensity = (*theMaterialTable)[mater 445 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
446 fElectronDensity = (*theMaterialTable)[mater 446 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
447 numberOfElements = (G4int)(*theMaterialTable 447 numberOfElements = (G4int)(*theMaterialTable)[materialIndex]->GetNumberOfElements();
448 448
449 G4int* thisMaterialZ = new G4int[numberOfEle 449 G4int* thisMaterialZ = new G4int[numberOfElements];
450 450
451 for( i = 0; i < numberOfElements; ++i ) 451 for( i = 0; i < numberOfElements; ++i )
452 { 452 {
453 thisMaterialZ[i] = (G4int)(*theMateri 453 thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
454 GetEleme 454 GetElement(i)->GetZ();
455 } 455 }
456 // fSandia = new G4SandiaTable(materialIndex 456 // fSandia = new G4SandiaTable(materialIndex);
457 fSandia = (*theMaterialTable)[materialIndex] 457 fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
458 G4SandiaTable thisMaterialSandiaTable(ma 458 G4SandiaTable thisMaterialSandiaTable(materialIndex);
459 fIntervalNumber = thisMaterialSandiaTable.Sa 459 fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
460 460 numberOfElements);
461 fIntervalNumber = thisMaterialSandiaTable.Sa 461 fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
462 ( thisMaterialZ , 462 ( thisMaterialZ ,
463 (*theMaterialTable)[mate 463 (*theMaterialTable)[materialIndex]->GetFractionVector() ,
464 numberOfElements, 464 numberOfElements,fIntervalNumber);
465 465
466 fIntervalNumber--; 466 fIntervalNumber--;
467 467
468 fEnergyInterval = G4DataVector(fIntervalNumb 468 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
469 fA1 = G4DataVector(fIntervalNumb 469 fA1 = G4DataVector(fIntervalNumber+2,0.0);
470 fA2 = G4DataVector(fIntervalNumb 470 fA2 = G4DataVector(fIntervalNumber+2,0.0);
471 fA3 = G4DataVector(fIntervalNumb 471 fA3 = G4DataVector(fIntervalNumber+2,0.0);
472 fA4 = G4DataVector(fIntervalNumb 472 fA4 = G4DataVector(fIntervalNumber+2,0.0);
473 473
474 /* 474 /*
475 fEnergyInterval = new G4double[fInterval 475 fEnergyInterval = new G4double[fIntervalNumber+2];
476 fA1 = new G4double[fInterval 476 fA1 = new G4double[fIntervalNumber+2];
477 fA2 = new G4double[fInterval 477 fA2 = new G4double[fIntervalNumber+2];
478 fA3 = new G4double[fInterval 478 fA3 = new G4double[fIntervalNumber+2];
479 fA4 = new G4double[fInterval 479 fA4 = new G4double[fIntervalNumber+2];
480 */ 480 */
481 for( i = 1; i <= fIntervalNumber; i++ ) 481 for( i = 1; i <= fIntervalNumber; i++ )
482 { 482 {
483 if((thisMaterialSandiaTable.GetPhotoAbsorpCo 483 if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
484 i > fIntervalNumber) 484 i > fIntervalNumber)
485 { 485 {
486 fEnergyInterval[i] = maxEnergyTran 486 fEnergyInterval[i] = maxEnergyTransfer;
487 fIntervalNumber = i; 487 fIntervalNumber = i;
488 break; 488 break;
489 } 489 }
490 fEnergyInterval[i] = thisMaterialSandiaTabl 490 fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
491 fA1[i] = thisMaterialSandiaTabl 491 fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
492 fA2[i] = thisMaterialSandiaTabl 492 fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
493 fA3[i] = thisMaterialSandiaTabl 493 fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
494 fA4[i] = thisMaterialSandiaTabl 494 fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
495 495
496 } 496 }
497 if(fEnergyInterval[fIntervalNumber] != maxEn 497 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
498 { 498 {
499 fIntervalNumber++; 499 fIntervalNumber++;
500 fEnergyInterval[fIntervalNumber] = ma 500 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
501 fA1[fIntervalNumber] = fA1[fIntervalN 501 fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
502 fA2[fIntervalNumber] = fA2[fIntervalN 502 fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
503 fA3[fIntervalNumber] = fA3[fIntervalN 503 fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
504 fA4[fIntervalNumber] = fA4[fIntervalN 504 fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
505 } 505 }
506 for(i=1;i<=fIntervalNumber;i++) 506 for(i=1;i<=fIntervalNumber;i++)
507 { 507 {
508 // G4cout<<fEnergyInterval[i]<<"\t"<<f 508 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
509 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e 509 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
510 } 510 }
511 // Now checking, if two borders are too clos 511 // Now checking, if two borders are too close together
512 512
513 for( i = 1; i < fIntervalNumber; i++ ) 513 for( i = 1; i < fIntervalNumber; i++ )
514 { 514 {
515 if(fEnergyInterval[i+1]-fEnergyInterva 515 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
516 1.5*fDelta*(fEnergyInterval[i+1]+fE 516 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
517 { 517 {
518 continue; 518 continue;
519 } 519 }
520 else 520 else
521 { 521 {
522 for( j = i; j < fIntervalNumber; j++ 522 for( j = i; j < fIntervalNumber; j++ )
523 { 523 {
524 fEnergyInterval[j] = fEnergyInterv 524 fEnergyInterval[j] = fEnergyInterval[j+1];
525 fA1[j] = fA1[j+1]; 525 fA1[j] = fA1[j+1];
526 fA2[j] = fA2[j+1]; 526 fA2[j] = fA2[j+1];
527 fA3[j] = fA3[j+1]; 527 fA3[j] = fA3[j+1];
528 fA4[j] = fA4[j+1]; 528 fA4[j] = fA4[j+1];
529 } 529 }
530 fIntervalNumber--; 530 fIntervalNumber--;
531 i--; 531 i--;
532 } 532 }
533 } 533 }
534 534
535 /* ********************************* 535 /* *********************************
536 fSplineEnergy = new G4double[fM 536 fSplineEnergy = new G4double[fMaxSplineSize];
537 fRePartDielectricConst = new G4double[fM 537 fRePartDielectricConst = new G4double[fMaxSplineSize];
538 fImPartDielectricConst = new G4double[fM 538 fImPartDielectricConst = new G4double[fMaxSplineSize];
539 fIntegralTerm = new G4double[fM 539 fIntegralTerm = new G4double[fMaxSplineSize];
540 fDifPAIxSection = new G4double[fM 540 fDifPAIxSection = new G4double[fMaxSplineSize];
541 fIntegralPAIxSection = new G4double[fM 541 fIntegralPAIxSection = new G4double[fMaxSplineSize];
542 542
543 for(i=0;i<fMaxSplineSize;i++) 543 for(i=0;i<fMaxSplineSize;i++)
544 { 544 {
545 fSplineEnergy[i] = 0.0; 545 fSplineEnergy[i] = 0.0;
546 fRePartDielectricConst[i] = 0.0; 546 fRePartDielectricConst[i] = 0.0;
547 fImPartDielectricConst[i] = 0.0; 547 fImPartDielectricConst[i] = 0.0;
548 fIntegralTerm[i] = 0.0; 548 fIntegralTerm[i] = 0.0;
549 fDifPAIxSection[i] = 0.0; 549 fDifPAIxSection[i] = 0.0;
550 fIntegralPAIxSection[i] = 0.0; 550 fIntegralPAIxSection[i] = 0.0;
551 } 551 }
552 */ //////////////////////// 552 */ ////////////////////////
553 553
554 // Preparation of fSplineEnergy array co 554 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
555 555
556 ComputeLowEnergyCof(); 556 ComputeLowEnergyCof();
557 G4double betaGammaSqRef = 557 G4double betaGammaSqRef =
558 fLorentzFactor[fRefGammaNumber]*fLorentzFa 558 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
559 559
560 NormShift(betaGammaSqRef); 560 NormShift(betaGammaSqRef);
561 SplainPAI(betaGammaSqRef); 561 SplainPAI(betaGammaSqRef);
562 562
563 // Preparation of integral PAI cross section 563 // Preparation of integral PAI cross section for input betaGammaSq
564 564
565 for(i = 1; i <= fSplineNumber; i++) 565 for(i = 1; i <= fSplineNumber; i++)
566 { 566 {
567 fDifPAIxSection[i] = DifPAIxSection(i 567 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
568 fdNdxCerenkov[i] = PAIdNdxCerenkov( 568 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
569 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 569 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
570 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 570 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
571 fdNdxResonance[i] = PAIdNdxResonance 571 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
572 } 572 }
573 IntegralPAIxSection(); 573 IntegralPAIxSection();
574 IntegralCerenkov(); 574 IntegralCerenkov();
575 IntegralMM(); 575 IntegralMM();
576 IntegralPlasmon(); 576 IntegralPlasmon();
577 IntegralResonance(); 577 IntegralResonance();
578 } 578 }
579 579
580 ////////////////////////////////////////////// 580 ////////////////////////////////////////////////////////////////////////////
581 // 581 //
582 // Destructor 582 // Destructor
583 583
584 G4PAIxSection::~G4PAIxSection() 584 G4PAIxSection::~G4PAIxSection()
585 { 585 {
586 /* ************************ 586 /* ************************
587 delete[] fSplineEnergy ; 587 delete[] fSplineEnergy ;
588 delete[] fRePartDielectricConst; 588 delete[] fRePartDielectricConst;
589 delete[] fImPartDielectricConst; 589 delete[] fImPartDielectricConst;
590 delete[] fIntegralTerm ; 590 delete[] fIntegralTerm ;
591 delete[] fDifPAIxSection ; 591 delete[] fDifPAIxSection ;
592 delete[] fIntegralPAIxSection ; 592 delete[] fIntegralPAIxSection ;
593 */ //////////////////////// 593 */ ////////////////////////
594 delete fMatSandiaMatrix; 594 delete fMatSandiaMatrix;
595 } 595 }
596 596
597 G4double G4PAIxSection::GetLorentzFactor(G4int 597 G4double G4PAIxSection::GetLorentzFactor(G4int j) const
598 { 598 {
599 return fLorentzFactor[j]; 599 return fLorentzFactor[j];
600 } 600 }
601 601
602 ////////////////////////////////////////////// 602 ////////////////////////////////////////////////////////////////////////
603 // 603 //
604 // Constructor with beta*gamma square value ca 604 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data
605 605
606 void G4PAIxSection::Initialize( const G4Materi 606 void G4PAIxSection::Initialize( const G4Material* material,
607 G4double maxEn 607 G4double maxEnergyTransfer,
608 G4double betaG 608 G4double betaGammaSq,
609 G4SandiaTable* 609 G4SandiaTable* sandia)
610 { 610 {
611 if(fVerbose > 0) 611 if(fVerbose > 0)
612 { 612 {
613 G4cout<<G4endl; 613 G4cout<<G4endl;
614 G4cout<<"G4PAIxSection::Initialize(...,G4S 614 G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
615 G4cout<<G4endl; 615 G4cout<<G4endl;
616 } 616 }
617 G4int i, j; 617 G4int i, j;
618 618
619 fSandia = sandia; 619 fSandia = sandia;
620 fIntervalNumber = sandia->GetMaxInterval(); 620 fIntervalNumber = sandia->GetMaxInterval();
621 fDensity = material->GetDensity(); 621 fDensity = material->GetDensity();
622 fElectronDensity = material->GetElectronDens 622 fElectronDensity = material->GetElectronDensity();
623 623
624 // fIntervalNumber--; 624 // fIntervalNumber--;
625 625
626 if( fVerbose > 0 ) 626 if( fVerbose > 0 )
627 { 627 {
628 G4cout<<"fDensity = "<<fDensity<<"\t"<<fEl 628 G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
629 } 629 }
630 fEnergyInterval = G4DataVector(fIntervalNumb 630 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
631 fA1 = G4DataVector(fIntervalNumb 631 fA1 = G4DataVector(fIntervalNumber+2,0.0);
632 fA2 = G4DataVector(fIntervalNumb 632 fA2 = G4DataVector(fIntervalNumber+2,0.0);
633 fA3 = G4DataVector(fIntervalNumb 633 fA3 = G4DataVector(fIntervalNumber+2,0.0);
634 fA4 = G4DataVector(fIntervalNumb 634 fA4 = G4DataVector(fIntervalNumber+2,0.0);
635 635
636 for( i = 1; i <= fIntervalNumber; i++ ) 636 for( i = 1; i <= fIntervalNumber; i++ )
637 { 637 {
638 if ( sandia->GetSandiaMatTablePAI(i-1,0) < 638 if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
639 { 639 {
640 fIntervalNumber--; 640 fIntervalNumber--;
641 continue; 641 continue;
642 } 642 }
643 if( ( sandia->GetSandiaMatTablePAI(i-1,0) 643 if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
644 { 644 {
645 fEnergyInterval[i] = maxEnergyTransfer; 645 fEnergyInterval[i] = maxEnergyTransfer;
646 fIntervalNumber = i; 646 fIntervalNumber = i;
647 break; 647 break;
648 } 648 }
649 fEnergyInterval[i] = sandia->GetSandiaMatT 649 fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
650 fA1[i] = sandia->GetSandiaMatT 650 fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
651 fA2[i] = sandia->GetSandiaMatT 651 fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
652 fA3[i] = sandia->GetSandiaMatT 652 fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
653 fA4[i] = sandia->GetSandiaMatT 653 fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
654 654
655 if( fVerbose > 0 ) 655 if( fVerbose > 0 )
656 { 656 {
657 G4cout<<i<<"\t"<<fEnergyInterval[i]/ke 657 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
658 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e 658 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
659 } 659 }
660 } 660 }
661 if( fVerbose > 0 ) G4cout<<"last i = "<<i<<" 661 if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
662 662
663 if( fEnergyInterval[fIntervalNumber] != maxE 663 if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
664 { 664 {
665 fIntervalNumber++; 665 fIntervalNumber++;
666 fEnergyInterval[fIntervalNumber] = maxEn 666 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
667 } 667 }
668 if( fVerbose > 0 ) 668 if( fVerbose > 0 )
669 { 669 {
670 for( i = 1; i <= fIntervalNumber; i++ ) 670 for( i = 1; i <= fIntervalNumber; i++ )
671 { 671 {
672 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< 672 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
673 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 673 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
674 } 674 }
675 } 675 }
676 if( fVerbose > 0 ) G4cout<<"Now checking, 676 if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
677 677
678 for( i = 1; i < fIntervalNumber; i++ ) 678 for( i = 1; i < fIntervalNumber; i++ )
679 { 679 {
680 if( fEnergyInterval[i+1]-fEnergyInterval[i 680 if( fEnergyInterval[i+1]-fEnergyInterval[i] >
681 1.5*fDelta*(fEnergyInterval[i+1]+fEne 681 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
682 else 682 else
683 { 683 {
684 if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEn 684 if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
685 685
686 for( j = i; j < fIntervalNumber; j++ ) 686 for( j = i; j < fIntervalNumber; j++ )
687 { 687 {
688 fEnergyInterval[j] = fEnergyInte 688 fEnergyInterval[j] = fEnergyInterval[j+1];
689 fA1[j] = fA1[j+1]; 689 fA1[j] = fA1[j+1];
690 fA2[j] = fA2[j+1]; 690 fA2[j] = fA2[j+1];
691 fA3[j] = fA3[j+1]; 691 fA3[j] = fA3[j+1];
692 fA4[j] = fA4[j+1]; 692 fA4[j] = fA4[j+1];
693 } 693 }
694 fIntervalNumber--; 694 fIntervalNumber--;
695 i--; 695 i--;
696 } 696 }
697 } 697 }
698 if( fVerbose > 0 ) 698 if( fVerbose > 0 )
699 { 699 {
700 for( i = 1; i <= fIntervalNumber; i++ ) 700 for( i = 1; i <= fIntervalNumber; i++ )
701 { 701 {
702 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< 702 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
703 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 703 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
704 } 704 }
705 } 705 }
706 // Preparation of fSplineEnergy array corres 706 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
707 707
708 ComputeLowEnergyCof(material); 708 ComputeLowEnergyCof(material);
709 709
710 G4double betaGammaSqRef = 710 G4double betaGammaSqRef =
711 fLorentzFactor[fRefGammaNumber]*fLorentzFa 711 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
712 712
713 NormShift(betaGammaSqRef); 713 NormShift(betaGammaSqRef);
714 SplainPAI(betaGammaSqRef); 714 SplainPAI(betaGammaSqRef);
715 715
716 // Preparation of integral PAI cross section 716 // Preparation of integral PAI cross section for input betaGammaSq
717 717
718 for( i = 1; i <= fSplineNumber; i++ ) 718 for( i = 1; i <= fSplineNumber; i++ )
719 { 719 {
720 fDifPAIxSection[i] = DifPAIxSection(i,bet 720 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
721 721
722 722
723 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,be 723 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
724 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 724 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
725 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,bet 725 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
726 fdNdxResonance[i] = PAIdNdxResonance(i,b 726 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
727 } 727 }
728 IntegralPAIxSection(); 728 IntegralPAIxSection();
729 IntegralCerenkov(); 729 IntegralCerenkov();
730 IntegralMM(); 730 IntegralMM();
731 IntegralPlasmon(); 731 IntegralPlasmon();
732 IntegralResonance(); 732 IntegralResonance();
733 733
734 for( i = 1; i <= fSplineNumber; i++ ) 734 for( i = 1; i <= fSplineNumber; i++ )
735 { 735 {
736 if(fVerbose>0) G4cout<<i<<"; w = "<<fSplin 736 if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
737 } 737 }
738 } 738 }
739 739
740 740
741 ////////////////////////////////////////////// 741 /////////////////////////////////////////////////////////////////////////
742 // 742 //
743 // Compute low energy cof. It reduces PAI xsc 743 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
744 // 744 //
745 745
746 void G4PAIxSection::ComputeLowEnergyCof(const 746 void G4PAIxSection::ComputeLowEnergyCof(const G4Material* material)
747 { 747 {
748 G4int i, numberOfElements = (G4int)material- 748 G4int i, numberOfElements = (G4int)material->GetNumberOfElements();
749 G4double sumZ = 0., sumCof = 0.; 749 G4double sumZ = 0., sumCof = 0.;
750 750
751 static const G4double p0 = 1.20923e+00; 751 static const G4double p0 = 1.20923e+00;
752 static const G4double p1 = 3.53256e-01; 752 static const G4double p1 = 3.53256e-01;
753 static const G4double p2 = -1.45052e-03; 753 static const G4double p2 = -1.45052e-03;
754 754
755 G4double* thisMaterialZ = new G4double[num 755 G4double* thisMaterialZ = new G4double[numberOfElements];
756 G4double* thisMaterialCof = new G4double[num 756 G4double* thisMaterialCof = new G4double[numberOfElements];
757 757
758 for( i = 0; i < numberOfElements; ++i ) 758 for( i = 0; i < numberOfElements; ++i )
759 { 759 {
760 thisMaterialZ[i] = material->GetElement(i) 760 thisMaterialZ[i] = material->GetElement(i)->GetZ();
761 sumZ += thisMaterialZ[i]; 761 sumZ += thisMaterialZ[i];
762 thisMaterialCof[i] = p0+p1*thisMaterialZ[i 762 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
763 } 763 }
764 for( i = 0; i < numberOfElements; ++i ) 764 for( i = 0; i < numberOfElements; ++i )
765 { 765 {
766 sumCof += thisMaterialCof[i]*thisMaterialZ 766 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
767 } 767 }
768 fLowEnergyCof = sumCof; 768 fLowEnergyCof = sumCof;
769 delete [] thisMaterialZ; 769 delete [] thisMaterialZ;
770 delete [] thisMaterialCof; 770 delete [] thisMaterialCof;
771 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof 771 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
772 } 772 }
773 773
774 ////////////////////////////////////////////// 774 /////////////////////////////////////////////////////////////////////////
775 // 775 //
776 // Compute low energy cof. It reduces PAI xsc 776 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
777 // 777 //
778 778
779 void G4PAIxSection::ComputeLowEnergyCof() 779 void G4PAIxSection::ComputeLowEnergyCof()
780 { 780 {
781 const G4MaterialTable* theMaterialTable = G4 781 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
782 G4int i, numberOfElements = (G4int)(*theMate 782 G4int i, numberOfElements = (G4int)(*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
783 G4double sumZ = 0., sumCof = 0.; 783 G4double sumZ = 0., sumCof = 0.;
784 784
785 const G4double p0 = 1.20923e+00; 785 const G4double p0 = 1.20923e+00;
786 const G4double p1 = 3.53256e-01; 786 const G4double p1 = 3.53256e-01;
787 const G4double p2 = -1.45052e-03; 787 const G4double p2 = -1.45052e-03;
788 788
789 G4double* thisMaterialZ = new G4double[num 789 G4double* thisMaterialZ = new G4double[numberOfElements];
790 G4double* thisMaterialCof = new G4double[num 790 G4double* thisMaterialCof = new G4double[numberOfElements];
791 791
792 for( i = 0; i < numberOfElements; ++i ) 792 for( i = 0; i < numberOfElements; ++i )
793 { 793 {
794 thisMaterialZ[i] = (*theMaterialTable)[fMa 794 thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
795 sumZ += thisMaterialZ[i]; 795 sumZ += thisMaterialZ[i];
796 thisMaterialCof[i] = p0+p1*thisMaterialZ[i 796 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
797 } 797 }
798 for( i = 0; i < numberOfElements; ++i ) 798 for( i = 0; i < numberOfElements; ++i )
799 { 799 {
800 sumCof += thisMaterialCof[i]*thisMaterialZ 800 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
801 } 801 }
802 fLowEnergyCof = sumCof; 802 fLowEnergyCof = sumCof;
803 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof 803 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
804 delete [] thisMaterialZ; 804 delete [] thisMaterialZ;
805 delete [] thisMaterialCof; 805 delete [] thisMaterialCof;
806 } 806 }
807 807
808 ////////////////////////////////////////////// 808 /////////////////////////////////////////////////////////////////////////
809 // 809 //
810 // General control function for class G4PAIxSe 810 // General control function for class G4PAIxSection
811 // 811 //
812 812
813 void G4PAIxSection::InitPAI() 813 void G4PAIxSection::InitPAI()
814 { 814 {
815 G4int i; 815 G4int i;
816 G4double betaGammaSq = fLorentzFactor[fRefG 816 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
817 fLorentzFactor[fRefG 817 fLorentzFactor[fRefGammaNumber] - 1;
818 818
819 // Preparation of integral PAI cross sectio 819 // Preparation of integral PAI cross section for reference gamma
820 820
821 NormShift(betaGammaSq); 821 NormShift(betaGammaSq);
822 SplainPAI(betaGammaSq); 822 SplainPAI(betaGammaSq);
823 823
824 IntegralPAIxSection(); 824 IntegralPAIxSection();
825 IntegralCerenkov(); 825 IntegralCerenkov();
826 IntegralMM(); 826 IntegralMM();
827 IntegralPlasmon(); 827 IntegralPlasmon();
828 IntegralResonance(); 828 IntegralResonance();
829 829
830 for(i = 0; i<= fSplineNumber; i++) 830 for(i = 0; i<= fSplineNumber; i++)
831 { 831 {
832 fPAItable[i][fRefGammaNumber] = fIntegra 832 fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
833 if(i != 0) 833 if(i != 0)
834 { 834 {
835 fPAItable[i][0] = fSplineEnergy[i]; 835 fPAItable[i][0] = fSplineEnergy[i];
836 } 836 }
837 } 837 }
838 fPAItable[0][0] = fSplineNumber; 838 fPAItable[0][0] = fSplineNumber;
839 839
840 for(G4int j = 1; j < 112; j++) // for 840 for(G4int j = 1; j < 112; j++) // for other gammas
841 { 841 {
842 if( j == fRefGammaNumber ) continue; 842 if( j == fRefGammaNumber ) continue;
843 843
844 betaGammaSq = fLorentzFactor[j]*fLorentz 844 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
845 845
846 for(i = 1; i <= fSplineNumber; i++) 846 for(i = 1; i <= fSplineNumber; i++)
847 { 847 {
848 fDifPAIxSection[i] = DifPAIxSection(i 848 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
849 fdNdxCerenkov[i] = PAIdNdxCerenkov( 849 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
850 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 850 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
851 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 851 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
852 fdNdxResonance[i] = PAIdNdxResonance 852 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
853 } 853 }
854 IntegralPAIxSection(); 854 IntegralPAIxSection();
855 IntegralCerenkov(); 855 IntegralCerenkov();
856 IntegralMM(); 856 IntegralMM();
857 IntegralPlasmon(); 857 IntegralPlasmon();
858 IntegralResonance(); 858 IntegralResonance();
859 859
860 for(i = 0; i <= fSplineNumber; i++) 860 for(i = 0; i <= fSplineNumber; i++)
861 { 861 {
862 fPAItable[i][j] = fIntegralPAIxSectio 862 fPAItable[i][j] = fIntegralPAIxSection[i];
863 } 863 }
864 } 864 }
865 865
866 } 866 }
867 867
868 ////////////////////////////////////////////// 868 ///////////////////////////////////////////////////////////////////////
869 // 869 //
870 // Shifting from borders to intervals Creation 870 // Shifting from borders to intervals Creation of first energy points
871 // 871 //
872 872
873 void G4PAIxSection::NormShift(G4double betaGam 873 void G4PAIxSection::NormShift(G4double betaGammaSq)
874 { 874 {
875 G4int i, j; 875 G4int i, j;
876 876
877 if(fVerbose>0) G4cout<<" G4PAIxSection: 877 if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
878 878
879 879
880 for( i = 1; i <= fIntervalNumber-1; i++ ) 880 for( i = 1; i <= fIntervalNumber-1; i++ )
881 { 881 {
882 for( j = 1; j <= 2; j++ ) 882 for( j = 1; j <= 2; j++ )
883 { 883 {
884 fSplineNumber = (i-1)*2 + j; 884 fSplineNumber = (i-1)*2 + j;
885 885
886 if( j == 1 ) fSplineEnergy[fSplineNumber 886 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
887 else fSplineEnergy[fSplineNumber 887 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
888 if(fVerbose>0) G4cout<<"cn = "<<fSplineN 888 if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
889 } 889 }
890 } 890 }
891 fIntegralTerm[1]=RutherfordIntegral(1,fEnerg 891 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
892 892
893 j = 1; 893 j = 1;
894 894
895 for( i = 2; i <= fSplineNumber; i++ ) 895 for( i = 2; i <= fSplineNumber; i++ )
896 { 896 {
897 if( fSplineEnergy[i]<fEnergyInterval[j+1] 897 if( fSplineEnergy[i]<fEnergyInterval[j+1] )
898 { 898 {
899 fIntegralTerm[i] = fIntegralTerm[i-1] 899 fIntegralTerm[i] = fIntegralTerm[i-1] +
900 RutherfordIntegral 900 RutherfordIntegral(j,fSplineEnergy[i-1],
901 901 fSplineEnergy[i] );
902 } 902 }
903 else 903 else
904 { 904 {
905 G4double x = RutherfordIntegral(j,fSpli 905 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
906 fEn 906 fEnergyInterval[j+1] );
907 j++; 907 j++;
908 fIntegralTerm[i] = fIntegralTerm[i-1] 908 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
909 RutherfordIntegral 909 RutherfordIntegral(j,fEnergyInterval[j],
910 910 fSplineEnergy[i] );
911 } 911 }
912 if(fVerbose>0) G4cout<<i<<" Shift: w = "< 912 if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
913 } 913 }
914 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine 914 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
915 fNormalizationCof *= fElectronDensity/fInteg 915 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
916 916
917 // G4cout<<"fNormalizationCof = "<<fNormaliz 917 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
918 918
919 // Calculation of PAI differrential 919 // Calculation of PAI differrential cross-section (1/(keV*cm))
920 // in the energy points near borders 920 // in the energy points near borders of energy intervals
921 921
922 for(G4int k = 1; k <= fIntervalNumber-1; k+ 922 for(G4int k = 1; k <= fIntervalNumber-1; k++ )
923 { 923 {
924 for( j = 1; j <= 2; j++ ) 924 for( j = 1; j <= 2; j++ )
925 { 925 {
926 i = (k-1)*2 + j; 926 i = (k-1)*2 + j;
927 fImPartDielectricConst[i] = fNormaliz 927 fImPartDielectricConst[i] = fNormalizationCof*
928 ImPartDie 928 ImPartDielectricConst(k,fSplineEnergy[i]);
929 fRePartDielectricConst[i] = fNormaliz 929 fRePartDielectricConst[i] = fNormalizationCof*
930 RePartDie 930 RePartDielectricConst(fSplineEnergy[i]);
931 fIntegralTerm[i] *= fNormalizationCof 931 fIntegralTerm[i] *= fNormalizationCof;
932 932
933 fDifPAIxSection[i] = DifPAIxSection(i 933 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
934 fdNdxCerenkov[i] = PAIdNdxCerenkov( 934 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
935 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 935 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
936 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 936 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
937 fdNdxResonance[i] = PAIdNdxResonan 937 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
938 if(fVerbose>0) G4cout<<i<<" Shift: w = "< 938 if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
939 } 939 }
940 } 940 }
941 941
942 } // end of NormShift 942 } // end of NormShift
943 943
944 ////////////////////////////////////////////// 944 /////////////////////////////////////////////////////////////////////////
945 // 945 //
946 // Creation of new energy points as geometrica 946 // Creation of new energy points as geometrical mean of existing
947 // one, calculation PAI_cs for them, while the 947 // one, calculation PAI_cs for them, while the error of logarithmic
948 // linear approximation would be smaller than 948 // linear approximation would be smaller than 'fError'
949 949
950 void G4PAIxSection::SplainPAI(G4double betaGam 950 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
951 { 951 {
952 G4int j, k = 1, i = 1; 952 G4int j, k = 1, i = 1;
953 953
954 if(fVerbose>0) G4cout<<" G 954 if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
955 955
956 while ( (i < fSplineNumber) && (fSplineNumbe 956 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
957 { 957 {
958 // if( std::abs(fSplineEnergy[i+1]-fEnerg 958 // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
959 if( fSplineEnergy[i+1] > fEnergyInterval[ 959 if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
960 { 960 {
961 k++; // Here next energy point is 961 k++; // Here next energy point is in next energy interval
962 i++; 962 i++;
963 if(fVerbose>0) G4cout<<" 963 if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
964 continue; 964 continue;
965 } 965 }
966 if(fVerbose>0) G4cout<<" out if: i 966 if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
967 967
968 // Shifting of arrayes 968 // Shifting of arrayes for inserting the geometrical
969 // average of 'i' and ' 969 // average of 'i' and 'i+1' energy points to 'i+1' place
970 fSplineNumber++; 970 fSplineNumber++;
971 971
972 for( j = fSplineNumber; j >= i+2; j-- ) 972 for( j = fSplineNumber; j >= i+2; j-- )
973 { 973 {
974 fSplineEnergy[j] = fSplineEn 974 fSplineEnergy[j] = fSplineEnergy[j-1];
975 fImPartDielectricConst[j] = fImPartDi 975 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
976 fRePartDielectricConst[j] = fRePartDi 976 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
977 fIntegralTerm[j] = fIntegral 977 fIntegralTerm[j] = fIntegralTerm[j-1];
978 978
979 fDifPAIxSection[j] = fDifPAIxSection[ 979 fDifPAIxSection[j] = fDifPAIxSection[j-1];
980 fdNdxCerenkov[j] = fdNdxCerenkov[j- 980 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
981 fdNdxMM[j] = fdNdxMM[j-1]; 981 fdNdxMM[j] = fdNdxMM[j-1];
982 fdNdxPlasmon[j] = fdNdxPlasmon[j-1 982 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
983 fdNdxResonance[j] = fdNdxResonance[j 983 fdNdxResonance[j] = fdNdxResonance[j-1];
984 } 984 }
985 G4double x1 = fSplineEnergy[i]; 985 G4double x1 = fSplineEnergy[i];
986 G4double x2 = fSplineEnergy[i+1]; 986 G4double x2 = fSplineEnergy[i+1];
987 G4double yy1 = fDifPAIxSection[i]; 987 G4double yy1 = fDifPAIxSection[i];
988 G4double y2 = fDifPAIxSection[i+1]; 988 G4double y2 = fDifPAIxSection[i+1];
989 989
990 if(fVerbose>0) G4cout<<"Spline: x1 = "<< 990 if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
991 991
992 992
993 G4double en1 = sqrt(x1*x2); 993 G4double en1 = sqrt(x1*x2);
994 // G4double en1 = 0.5*(x1 + x2); 994 // G4double en1 = 0.5*(x1 + x2);
995 995
996 996
997 fSplineEnergy[i+1] = en1; 997 fSplineEnergy[i+1] = en1;
998 998
999 // Calculation of logarithmic 999 // Calculation of logarithmic linear approximation
1000 // in this (enr) energy poin 1000 // in this (enr) energy point, which number is 'i+1' now
1001 1001
1002 G4double a = log10(y2/yy1)/log10(x2/x1) 1002 G4double a = log10(y2/yy1)/log10(x2/x1);
1003 G4double b = log10(yy1) - a*log10(x1); 1003 G4double b = log10(yy1) - a*log10(x1);
1004 G4double y = a*log10(en1) + b; 1004 G4double y = a*log10(en1) + b;
1005 1005
1006 y = pow(10.,y); 1006 y = pow(10.,y);
1007 1007
1008 // Calculation of the PAI di 1008 // Calculation of the PAI dif. cross-section at this point
1009 1009
1010 fImPartDielectricConst[i+1] = fNormaliz 1010 fImPartDielectricConst[i+1] = fNormalizationCof*
1011 ImPartDie 1011 ImPartDielectricConst(k,fSplineEnergy[i+1]);
1012 fRePartDielectricConst[i+1] = fNormaliz 1012 fRePartDielectricConst[i+1] = fNormalizationCof*
1013 RePartDie 1013 RePartDielectricConst(fSplineEnergy[i+1]);
1014 fIntegralTerm[i+1] = fIntegralTerm[i] + 1014 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1015 RutherfordIntegral 1015 RutherfordIntegral(k,fSplineEnergy[i],
1016 1016 fSplineEnergy[i+1]);
1017 1017
1018 fDifPAIxSection[i+1] = DifPAIxSection(i 1018 fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1019 fdNdxCerenkov[i+1] = PAIdNdxCerenkov( 1019 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1020 fdNdxMM[i+1] = PAIdNdxMM(i+1,be 1020 fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1021 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i 1021 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1022 fdNdxResonance[i+1] = PAIdNdxResonance 1022 fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1023 1023
1024 // Condition for next divis 1024 // Condition for next division of this segment or to pass
1025 1025
1026 if(fVerbose>0) G4cout<<"Spline, a = "<<a< 1026 if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1027 1027
1028 // to higher energies 1028 // to higher energies
1029 1029
1030 G4double x = 2*(fDifPAIxSection[i+1] - 1030 G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1031 1031
1032 G4double delta = 2.*(fSplineEnergy[i+1] 1032 G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1033 1033
1034 if( x < 0 ) 1034 if( x < 0 )
1035 { 1035 {
1036 x = -x; 1036 x = -x;
1037 } 1037 }
1038 if( x > fError && fSplineNumber < fMaxS 1038 if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1039 { 1039 {
1040 continue; // next division 1040 continue; // next division
1041 } 1041 }
1042 i += 2; // pass to next segment 1042 i += 2; // pass to next segment
1043 1043
1044 // Loop checking, 03-Aug-2015, Vladimir 1044 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1045 } // close 'while' 1045 } // close 'while'
1046 1046
1047 } // end of SplainPAI 1047 } // end of SplainPAI
1048 1048
1049 1049
1050 ///////////////////////////////////////////// 1050 ////////////////////////////////////////////////////////////////////
1051 // 1051 //
1052 // Integration over electrons that could be c 1052 // Integration over electrons that could be considered
1053 // quasi-free at energy transfer of interest 1053 // quasi-free at energy transfer of interest
1054 1054
1055 G4double G4PAIxSection::RutherfordIntegral( G 1055 G4double G4PAIxSection::RutherfordIntegral( G4int k,
1056 G 1056 G4double x1,
1057 1057 G4double x2 )
1058 { 1058 {
1059 G4double c1, c2, c3; 1059 G4double c1, c2, c3;
1060 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "< 1060 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
1061 c1 = (x2 - x1)/x1/x2; 1061 c1 = (x2 - x1)/x1/x2;
1062 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2; 1062 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
1063 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/ 1063 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1064 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = " 1064 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
1065 1065
1066 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA 1066 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
1067 1067
1068 } // end of RutherfordIntegral 1068 } // end of RutherfordIntegral
1069 1069
1070 1070
1071 ///////////////////////////////////////////// 1071 /////////////////////////////////////////////////////////////////
1072 // 1072 //
1073 // Imaginary part of dielectric constant 1073 // Imaginary part of dielectric constant
1074 // (G4int k - interval number, G4double en1 - 1074 // (G4int k - interval number, G4double en1 - energy point)
1075 1075
1076 G4double G4PAIxSection::ImPartDielectricConst 1076 G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
1077 1077 G4double energy1 )
1078 { 1078 {
1079 G4double energy2,energy3,energy4,result; 1079 G4double energy2,energy3,energy4,result;
1080 1080
1081 energy2 = energy1*energy1; 1081 energy2 = energy1*energy1;
1082 energy3 = energy2*energy1; 1082 energy3 = energy2*energy1;
1083 energy4 = energy3*energy1; 1083 energy4 = energy3*energy1;
1084 1084
1085 result = fA1[k]/energy1+fA2[k]/energy2+fA3 1085 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1086 result *=hbarc/energy1; 1086 result *=hbarc/energy1;
1087 1087
1088 return result; 1088 return result;
1089 1089
1090 } // end of ImPartDielectricConst 1090 } // end of ImPartDielectricConst
1091 1091
1092 ///////////////////////////////////////////// 1092 /////////////////////////////////////////////////////////////////
1093 // 1093 //
1094 // Returns lambda of photon with energy1 in c 1094 // Returns lambda of photon with energy1 in current material
1095 1095
1096 G4double G4PAIxSection::GetPhotonRange( G4dou 1096 G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
1097 { 1097 {
1098 G4int i; 1098 G4int i;
1099 G4double energy2, energy3, energy4, result, 1099 G4double energy2, energy3, energy4, result, lambda;
1100 1100
1101 energy2 = energy1*energy1; 1101 energy2 = energy1*energy1;
1102 energy3 = energy2*energy1; 1102 energy3 = energy2*energy1;
1103 energy4 = energy3*energy1; 1103 energy4 = energy3*energy1;
1104 1104
>> 1105 // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
>> 1106 // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
>> 1107 // result *= fDensity;
>> 1108
1105 for( i = 1; i <= fIntervalNumber; i++ ) 1109 for( i = 1; i <= fIntervalNumber; i++ )
1106 { 1110 {
1107 if( energy1 < fEnergyInterval[i]) break; 1111 if( energy1 < fEnergyInterval[i]) break;
1108 } 1112 }
1109 i--; 1113 i--;
1110 if(i == 0) i = 1; 1114 if(i == 0) i = 1;
1111 1115
1112 result = fA1[i]/energy1+fA2[i]/energy2+fA3[ 1116 result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1113 1117
1114 if( result > DBL_MIN ) lambda = 1./result; 1118 if( result > DBL_MIN ) lambda = 1./result;
1115 else lambda = DBL_MAX; 1119 else lambda = DBL_MAX;
1116 1120
1117 return lambda; 1121 return lambda;
1118 } 1122 }
1119 1123
1120 ///////////////////////////////////////////// 1124 /////////////////////////////////////////////////////////////////
1121 // 1125 //
1122 // Return lambda of electron with energy1 in 1126 // Return lambda of electron with energy1 in current material
1123 // parametrisation from NIM A554(2005)474-493 1127 // parametrisation from NIM A554(2005)474-493
1124 1128
1125 G4double G4PAIxSection::GetElectronRange( G4d 1129 G4double G4PAIxSection::GetElectronRange( G4double energy )
1126 { 1130 {
1127 G4double range; 1131 G4double range;
1128 /* 1132 /*
1129 const G4MaterialTable* theMaterialTable = G 1133 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1130 1134
1131 G4double Z = (*theMaterialTable)[fMaterialI 1135 G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
1132 G4double A = (*theMaterialTable)[fMaterialI 1136 G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
1133 1137
1134 energy /= keV; // energy in keV in parametr 1138 energy /= keV; // energy in keV in parametrised formula
1135 1139
1136 if (energy < 10.) 1140 if (energy < 10.)
1137 { 1141 {
1138 range = 3.872e-3*A/Z; 1142 range = 3.872e-3*A/Z;
1139 range *= pow(energy,1.492); 1143 range *= pow(energy,1.492);
1140 } 1144 }
1141 else 1145 else
1142 { 1146 {
1143 range = 6.97e-3*pow(energy,1.6); 1147 range = 6.97e-3*pow(energy,1.6);
1144 } 1148 }
1145 */ 1149 */
1146 // Blum&Rolandi Particle Detection with Dri 1150 // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
1147 1151
1148 G4double cofA = 5.37e-4*g/cm2/keV; 1152 G4double cofA = 5.37e-4*g/cm2/keV;
1149 G4double cofB = 0.9815; 1153 G4double cofB = 0.9815;
1150 G4double cofC = 3.123e-3/keV; 1154 G4double cofC = 3.123e-3/keV;
1151 // energy /= keV; 1155 // energy /= keV;
1152 1156
1153 range = cofA*energy*( 1 - cofB/(1 + cofC*en 1157 range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
1154 1158
1155 // range *= g/cm2; 1159 // range *= g/cm2;
1156 range /= fDensity; 1160 range /= fDensity;
1157 1161
1158 return range; 1162 return range;
1159 } 1163 }
1160 1164
1161 ///////////////////////////////////////////// 1165 //////////////////////////////////////////////////////////////////////////////
1162 // 1166 //
1163 // Real part of dielectric constant minus uni 1167 // Real part of dielectric constant minus unit: epsilon_1 - 1
1164 // (G4double enb - energy point) 1168 // (G4double enb - energy point)
1165 // 1169 //
1166 1170
1167 G4double G4PAIxSection::RePartDielectricConst 1171 G4double G4PAIxSection::RePartDielectricConst(G4double enb)
1168 { 1172 {
1169 G4double x0, x02, x03, x04, x05, x1, x2, x 1173 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
1170 c1, c2, c3, cof1, cof2, xln1, xln 1174 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
1171 1175
1172 x0 = enb; 1176 x0 = enb;
1173 result = 0; 1177 result = 0;
1174 1178
1175 for(G4int i=1;i<=fIntervalNumber-1;i++) 1179 for(G4int i=1;i<=fIntervalNumber-1;i++)
1176 { 1180 {
1177 x1 = fEnergyInterval[i]; 1181 x1 = fEnergyInterval[i];
1178 x2 = fEnergyInterval[i+1]; 1182 x2 = fEnergyInterval[i+1];
1179 xx1 = x1 - x0; 1183 xx1 = x1 - x0;
1180 xx2 = x2 - x0; 1184 xx2 = x2 - x0;
1181 xx12 = xx2/xx1; 1185 xx12 = xx2/xx1;
1182 1186
1183 if(xx12<0) 1187 if(xx12<0)
1184 { 1188 {
1185 xx12 = -xx12; 1189 xx12 = -xx12;
1186 } 1190 }
1187 xln1 = log(x2/x1); 1191 xln1 = log(x2/x1);
1188 xln2 = log(xx12); 1192 xln2 = log(xx12);
1189 xln3 = log((x2 + x0)/(x1 + x0)); 1193 xln3 = log((x2 + x0)/(x1 + x0));
1190 x02 = x0*x0; 1194 x02 = x0*x0;
1191 x03 = x02*x0; 1195 x03 = x02*x0;
1192 x04 = x03*x0; 1196 x04 = x03*x0;
1193 x05 = x04*x0; 1197 x05 = x04*x0;
1194 c1 = (x2 - x1)/x1/x2; 1198 c1 = (x2 - x1)/x1/x2;
1195 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2; 1199 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
1196 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/ 1200 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1197 1201
1198 result -= (fA1[i]/x02 + fA3[i]/x04)*xln 1202 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
1199 result -= (fA2[i]/x02 + fA4[i]/x04)*c1; 1203 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
1200 result -= fA3[i]*c2/2/x02; 1204 result -= fA3[i]*c2/2/x02;
1201 result -= fA4[i]*c3/3/x02; 1205 result -= fA4[i]*c3/3/x02;
1202 1206
1203 cof1 = fA1[i]/x02 + fA3[i]/x04; 1207 cof1 = fA1[i]/x02 + fA3[i]/x04;
1204 cof2 = fA2[i]/x03 + fA4[i]/x05; 1208 cof2 = fA2[i]/x03 + fA4[i]/x05;
1205 1209
1206 result += 0.5*(cof1 +cof2)*xln2; 1210 result += 0.5*(cof1 +cof2)*xln2;
1207 result += 0.5*(cof1 - cof2)*xln3; 1211 result += 0.5*(cof1 - cof2)*xln3;
1208 } 1212 }
1209 result *= 2*hbarc/pi; 1213 result *= 2*hbarc/pi;
1210 1214
1211 return result; 1215 return result;
1212 1216
1213 } // end of RePartDielectricConst 1217 } // end of RePartDielectricConst
1214 1218
1215 ///////////////////////////////////////////// 1219 //////////////////////////////////////////////////////////////////////
1216 // 1220 //
1217 // PAI differential cross-section in terms of 1221 // PAI differential cross-section in terms of
1218 // simplified Allison's equation 1222 // simplified Allison's equation
1219 // 1223 //
1220 1224
1221 G4double G4PAIxSection::DifPAIxSection( G4int << 1225 G4double G4PAIxSection::DifPAIxSection( G4int i ,
>> 1226 G4double betaGammaSq )
1222 { 1227 {
1223 G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,resul 1228 G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
1224 1229
1225 G4double betaBohr = fine_structure_const; 1230 G4double betaBohr = fine_structure_const;
>> 1231 // G4double betaBohr2 = fine_structure_const*fine_structure_const;
>> 1232 // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
>> 1233
1226 G4double be2 = betaGammaSq/(1 + betaGamma 1234 G4double be2 = betaGammaSq/(1 + betaGammaSq);
1227 G4double beta = std::sqrt(be2); << 1235 G4double beta = sqrt(be2);
>> 1236 // G4double be3 = beta*be2;
1228 1237
1229 cof = 1.; 1238 cof = 1.;
1230 x1 = std::log(2*electron_mass_c2/fSplineE << 1239 x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
1231 1240
1232 if( betaGammaSq < 0.01 ) x2 = std::log(be2 << 1241 if( betaGammaSq < 0.01 ) x2 = log(be2);
1233 else 1242 else
1234 { 1243 {
1235 x2 = -log( (1/betaGammaSq - fRePartDiele 1244 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1236 (1/betaGammaSq - fRePartDiele 1245 (1/betaGammaSq - fRePartDielectricConst[i]) +
1237 fImPartDielectricConst[i]*fIm 1246 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
1238 } 1247 }
1239 if( fImPartDielectricConst[i] == 0.0 ||bet 1248 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
1240 { 1249 {
1241 x6 = 0.; 1250 x6 = 0.;
1242 } 1251 }
1243 else 1252 else
1244 { 1253 {
1245 x3 = -fRePartDielectricConst[i] + 1/beta 1254 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
1246 x5 = -1 - fRePartDielectricConst[i] + 1255 x5 = -1 - fRePartDielectricConst[i] +
1247 be2*((1 +fRePartDielectricConst[i]) 1256 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1248 fImPartDielectricConst[i]*fImPartDi 1257 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1249 1258
1250 x7 = atan2(fImPartDielectricConst[i],x3) 1259 x7 = atan2(fImPartDielectricConst[i],x3);
1251 x6 = x5 * x7; 1260 x6 = x5 * x7;
1252 } 1261 }
>> 1262 // if(fImPartDielectricConst[i] == 0) x6 = 0.;
1253 1263
1254 x4 = ((x1 + x2)*fImPartDielectricConst[i] 1264 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
1255 1265
>> 1266 // if( x4 < 0.0 ) x4 = 0.0;
>> 1267
1256 x8 = (1 + fRePartDielectricConst[i])*(1 + 1268 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1257 fImPartDielectricConst[i]*fImPartDiel 1269 fImPartDielectricConst[i]*fImPartDielectricConst[i];
1258 1270
1259 result = (x4 + cof*fIntegralTerm[i]/fSplin 1271 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
1260 1272
1261 if( result < 1.0e-8 ) result = 1.0e-8; 1273 if( result < 1.0e-8 ) result = 1.0e-8;
1262 1274
1263 result *= fine_structure_const/be2/pi; 1275 result *= fine_structure_const/be2/pi;
1264 1276
1265 // low energy correction 1277 // low energy correction
1266 1278
1267 G4double lowCof = fLowEnergyCof; // 6.0 ; 1279 G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
1268 1280
1269 result *= (1 - std::exp(-beta/betaBohr/low << 1281 result *= (1 - exp(-beta/betaBohr/lowCof));
1270 if(x8 >= 0.0) << 1282
>> 1283
>> 1284 // result *= (1 - exp(-be2/betaBohr2/lowCof));
>> 1285
>> 1286 // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
>> 1287
>> 1288 // result *= (1 - exp(-be4/betaBohr4/lowCof));
>> 1289
>> 1290 if(fDensity >= 0.1)
1271 { 1291 {
1272 result /= x8; << 1292 result /= x8;
1273 } 1293 }
1274 return result; 1294 return result;
1275 1295
1276 } // end of DifPAIxSection 1296 } // end of DifPAIxSection
1277 1297
1278 ///////////////////////////////////////////// 1298 //////////////////////////////////////////////////////////////////////////
1279 // 1299 //
1280 // Calculation od dN/dx of collisions with cr 1300 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
1281 1301
1282 G4double G4PAIxSection::PAIdNdxCerenkov( G4in 1302 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
1283 G4do 1303 G4double betaGammaSq )
1284 { 1304 {
1285 G4double logarithm, x3, x5, argument, modu 1305 G4double logarithm, x3, x5, argument, modul2, dNdxC;
1286 G4double be2, betaBohr2, cofBetaBohr; 1306 G4double be2, betaBohr2, cofBetaBohr;
1287 1307
1288 cofBetaBohr = 4.0; 1308 cofBetaBohr = 4.0;
1289 betaBohr2 = fine_structure_const*fine_stru << 1309 betaBohr2 = fine_structure_const*fine_structure_const;
1290 G4double betaBohr4 = betaBohr2*betaBohr2*c << 1310 G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1291 1311
1292 be2 = betaGammaSq/(1 + betaGammaSq); 1312 be2 = betaGammaSq/(1 + betaGammaSq);
1293 G4double be4 = be2*be2; 1313 G4double be4 = be2*be2;
1294 1314
1295 if( betaGammaSq < 0.01 ) logarithm = std:: << 1315 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1296 else 1316 else
1297 { 1317 {
1298 logarithm = -log( (1/betaGammaSq - fReP 1318 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1299 (1/betaGammaSq - fReP 1319 (1/betaGammaSq - fRePartDielectricConst[i]) +
1300 fImPartDielectricCons 1320 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1301 logarithm += log(1+1.0/betaGammaSq); 1321 logarithm += log(1+1.0/betaGammaSq);
1302 } 1322 }
1303 1323
1304 if( fImPartDielectricConst[i] == 0.0 || be 1324 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1305 { 1325 {
1306 argument = 0.0; 1326 argument = 0.0;
1307 } 1327 }
1308 else 1328 else
1309 { 1329 {
1310 x3 = -fRePartDielectricConst[i] + 1.0/be 1330 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1311 x5 = -1.0 - fRePartDielectricConst[i] + 1331 x5 = -1.0 - fRePartDielectricConst[i] +
1312 be2*((1.0 +fRePartDielectricConst[i 1332 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1313 fImPartDielectricConst[i]*fImPartDi 1333 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1314 if( x3 == 0.0 ) argument = 0.5*pi; 1334 if( x3 == 0.0 ) argument = 0.5*pi;
1315 else argument = std::atan2(fI << 1335 else argument = atan2(fImPartDielectricConst[i],x3);
1316 argument *= x5 ; 1336 argument *= x5 ;
1317 } 1337 }
1318 dNdxC = ( logarithm*fImPartDielectricConst 1338 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1319 1339
1320 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1340 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1321 1341
1322 dNdxC *= fine_structure_const/be2/pi; 1342 dNdxC *= fine_structure_const/be2/pi;
1323 1343
1324 dNdxC *= (1-std::exp(-be4/betaBohr4)); << 1344 dNdxC *= (1-exp(-be4/betaBohr4));
1325 1345
1326 modul2 = (1.0 + fRePartDielectricConst[i]) << 1346 if(fDensity >= 0.1)
1327 fImPartDielectricConst[i]*fImPartDielect <<
1328 if(modul2 >= 0.0) <<
1329 { 1347 {
1330 dNdxC /= modul2; << 1348 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
>> 1349 fImPartDielectricConst[i]*fImPartDielectricConst[i];
>> 1350 dNdxC /= modul2;
1331 } 1351 }
1332 return dNdxC; 1352 return dNdxC;
1333 1353
1334 } // end of PAIdNdxCerenkov 1354 } // end of PAIdNdxCerenkov
1335 1355
1336 ///////////////////////////////////////////// 1356 //////////////////////////////////////////////////////////////////////////
1337 // 1357 //
1338 // Calculation od dN/dx of collisions of MM w 1358 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1339 1359
1340 G4double G4PAIxSection::PAIdNdxMM( G4int i 1360 G4double G4PAIxSection::PAIdNdxMM( G4int i ,
1341 G4do 1361 G4double betaGammaSq )
1342 { 1362 {
1343 G4double logarithm, x3, x5, argument, dNdx 1363 G4double logarithm, x3, x5, argument, dNdxC;
1344 G4double be2, be4, betaBohr2,betaBohr4,cof 1364 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
1345 1365
1346 cofBetaBohr = 4.0; 1366 cofBetaBohr = 4.0;
1347 betaBohr2 = fine_structure_const*fine_st 1367 betaBohr2 = fine_structure_const*fine_structure_const;
1348 betaBohr4 = betaBohr2*betaBohr2*cofBetaB 1368 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1349 1369
1350 be2 = betaGammaSq/(1 + betaGammaSq); 1370 be2 = betaGammaSq/(1 + betaGammaSq);
1351 be4 = be2*be2; 1371 be4 = be2*be2;
1352 1372
1353 if( betaGammaSq < 0.01 ) logarithm = log(1 1373 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1354 else 1374 else
1355 { 1375 {
1356 logarithm = -log( (1/betaGammaSq - fReP 1376 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1357 (1/betaGammaSq - fReP 1377 (1/betaGammaSq - fRePartDielectricConst[i]) +
1358 fImPartDielectricCons 1378 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1359 logarithm += log(1+1.0/betaGammaSq); 1379 logarithm += log(1+1.0/betaGammaSq);
1360 } 1380 }
1361 1381
1362 if( fImPartDielectricConst[i] == 0.0 || be 1382 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1363 { 1383 {
1364 argument = 0.0; 1384 argument = 0.0;
1365 } 1385 }
1366 else 1386 else
1367 { 1387 {
1368 x3 = -fRePartDielectricConst[i] + 1.0/be 1388 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1369 x5 = be2*( 1.0 + fRePartDielectricConst[ 1389 x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1370 if( x3 == 0.0 ) argument = 0.5*pi; 1390 if( x3 == 0.0 ) argument = 0.5*pi;
1371 else argument = atan2(fImPart 1391 else argument = atan2(fImPartDielectricConst[i],x3);
1372 argument *= x5 ; 1392 argument *= x5 ;
1373 } 1393 }
1374 dNdxC = ( logarithm*fImPartDielectricConst 1394 dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1375 1395
1376 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1396 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1377 1397
1378 dNdxC *= fine_structure_const/be2/pi; 1398 dNdxC *= fine_structure_const/be2/pi;
1379 1399
1380 dNdxC *= (1-std::exp(-be4/betaBohr4)); << 1400 dNdxC *= (1-exp(-be4/betaBohr4));
1381 return dNdxC; 1401 return dNdxC;
1382 1402
1383 } // end of PAIdNdxMM 1403 } // end of PAIdNdxMM
1384 1404
1385 ///////////////////////////////////////////// 1405 //////////////////////////////////////////////////////////////////////////
1386 // 1406 //
1387 // Calculation od dN/dx of collisions with cr 1407 // Calculation od dN/dx of collisions with creation of longitudinal EM
1388 // excitations (plasmons, delta-electrons) 1408 // excitations (plasmons, delta-electrons)
1389 1409
1390 G4double G4PAIxSection::PAIdNdxPlasmon( G4int 1410 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
1391 G4dou 1411 G4double betaGammaSq )
1392 { 1412 {
1393 G4double resonance, modul2, dNdxP, cof = 1 1413 G4double resonance, modul2, dNdxP, cof = 1.;
1394 G4double be2, betaBohr; 1414 G4double be2, betaBohr;
1395 1415
1396 betaBohr = fine_structure_const; 1416 betaBohr = fine_structure_const;
1397 be2 = betaGammaSq/(1 + betaGammaSq); 1417 be2 = betaGammaSq/(1 + betaGammaSq);
1398 1418
1399 G4double beta = std::sqrt(be2); << 1419 G4double beta = sqrt(be2);
1400 1420
1401 resonance = std::log(2*electron_mass_c2*be << 1421 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1402 resonance *= fImPartDielectricConst[i]/hba 1422 resonance *= fImPartDielectricConst[i]/hbarc;
1403 1423
>> 1424
1404 dNdxP = ( resonance + cof*fIntegralTerm[i] 1425 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1405 1426
1406 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1427 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1407 1428
1408 dNdxP *= fine_structure_const/be2/pi; 1429 dNdxP *= fine_structure_const/be2/pi;
1409 1430
1410 dNdxP *= (1 - std::exp(-beta/betaBohr/fLo << 1431 dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof));
1411 1432
1412 modul2 = (1 + fRePartDielectricConst[i])*( << 1433 // dNdxP *= (1-exp(-be4/betaBohr4));
1413 fImPartDielectricConst[i]*fImPartDielect << 1434
1414 if( modul2 >= 0.0 ) << 1435 if( fDensity >= 0.1 )
1415 { 1436 {
>> 1437 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
>> 1438 fImPartDielectricConst[i]*fImPartDielectricConst[i];
1416 dNdxP /= modul2; 1439 dNdxP /= modul2;
1417 } 1440 }
1418 return dNdxP; 1441 return dNdxP;
1419 1442
1420 } // end of PAIdNdxPlasmon 1443 } // end of PAIdNdxPlasmon
1421 1444
1422 ///////////////////////////////////////////// 1445 //////////////////////////////////////////////////////////////////////////
1423 // 1446 //
1424 // Calculation od dN/dx of collisions with cr 1447 // Calculation od dN/dx of collisions with creation of longitudinal EM
1425 // resonance excitations (plasmons, delta-ele 1448 // resonance excitations (plasmons, delta-electrons)
1426 1449
1427 G4double G4PAIxSection::PAIdNdxResonance( G4i 1450 G4double G4PAIxSection::PAIdNdxResonance( G4int i ,
1428 G4dou 1451 G4double betaGammaSq )
1429 { 1452 {
1430 G4double resonance, modul2, dNdxP; 1453 G4double resonance, modul2, dNdxP;
1431 G4double be2, be4, betaBohr2, betaBohr4, c 1454 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1432 1455
1433 cofBetaBohr = 4.0; 1456 cofBetaBohr = 4.0;
1434 betaBohr2 = fine_structure_const*fine_st 1457 betaBohr2 = fine_structure_const*fine_structure_const;
1435 betaBohr4 = betaBohr2*betaBohr2*cofBetaB 1458 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1436 1459
1437 be2 = betaGammaSq/(1 + betaGammaSq); 1460 be2 = betaGammaSq/(1 + betaGammaSq);
1438 be4 = be2*be2; 1461 be4 = be2*be2;
1439 1462
1440 resonance = log(2*electron_mass_c2*be2/fSp 1463 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1441 resonance *= fImPartDielectricConst[i]/hba 1464 resonance *= fImPartDielectricConst[i]/hbarc;
1442 1465
>> 1466
1443 dNdxP = resonance; 1467 dNdxP = resonance;
1444 1468
1445 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1469 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1446 1470
1447 dNdxP *= fine_structure_const/be2/pi; 1471 dNdxP *= fine_structure_const/be2/pi;
1448 dNdxP *= (1 - std::exp(-be4/betaBohr4)); << 1472 dNdxP *= (1-exp(-be4/betaBohr4));
1449 1473
1450 modul2 = (1 + fRePartDielectricConst[i])*( << 1474 if( fDensity >= 0.1 )
1451 fImPartDielectricConst[i]*fImPartDielect <<
1452 if( modul2 >= 0.0 ) <<
1453 { 1475 {
>> 1476 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
>> 1477 fImPartDielectricConst[i]*fImPartDielectricConst[i];
1454 dNdxP /= modul2; 1478 dNdxP /= modul2;
1455 } 1479 }
1456 return dNdxP; 1480 return dNdxP;
1457 1481
1458 } // end of PAIdNdxResonance 1482 } // end of PAIdNdxResonance
1459 1483
1460 ///////////////////////////////////////////// 1484 ////////////////////////////////////////////////////////////////////////
1461 // 1485 //
1462 // Calculation of the PAI integral cross-sect 1486 // Calculation of the PAI integral cross-section
1463 // fIntegralPAIxSection[1] = specific primary 1487 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1464 // and fIntegralPAIxSection[0] = mean energy 1488 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
1465 1489
1466 void G4PAIxSection::IntegralPAIxSection() 1490 void G4PAIxSection::IntegralPAIxSection()
1467 { 1491 {
1468 fIntegralPAIxSection[fSplineNumber] = 0; 1492 fIntegralPAIxSection[fSplineNumber] = 0;
1469 fIntegralPAIdEdx[fSplineNumber] = 0; 1493 fIntegralPAIdEdx[fSplineNumber] = 0;
1470 fIntegralPAIxSection[0] = 0; 1494 fIntegralPAIxSection[0] = 0;
1471 G4int i, k = fIntervalNumber -1; 1495 G4int i, k = fIntervalNumber -1;
1472 1496
1473 for( i = fSplineNumber-1; i >= 1; i--) 1497 for( i = fSplineNumber-1; i >= 1; i--)
1474 { 1498 {
1475 if(fSplineEnergy[i] >= fEnergyInterval[k] 1499 if(fSplineEnergy[i] >= fEnergyInterval[k])
1476 { 1500 {
1477 fIntegralPAIxSection[i] = fIntegralPAIx 1501 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1478 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[ 1502 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1479 } 1503 }
1480 else 1504 else
1481 { 1505 {
1482 fIntegralPAIxSection[i] = fIntegralPAIx 1506 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1483 SumOverBor 1507 SumOverBorder(i+1,fEnergyInterval[k]);
1484 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[ 1508 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1485 SumOverBor 1509 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1486 k--; 1510 k--;
1487 } 1511 }
1488 if(fVerbose>0) G4cout<<"i = "<<i<<"; k = 1512 if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1489 } 1513 }
1490 } // end of IntegralPAIxSection 1514 } // end of IntegralPAIxSection
1491 1515
1492 ///////////////////////////////////////////// 1516 ////////////////////////////////////////////////////////////////////////
1493 // 1517 //
1494 // Calculation of the PAI Cerenkov integral c 1518 // Calculation of the PAI Cerenkov integral cross-section
1495 // fIntegralCrenkov[1] = specific Crenkov ion 1519 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1496 // and fIntegralCerenkov[0] = mean Cerenkov l 1520 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
1497 1521
1498 void G4PAIxSection::IntegralCerenkov() 1522 void G4PAIxSection::IntegralCerenkov()
1499 { 1523 {
1500 G4int i, k; 1524 G4int i, k;
1501 fIntegralCerenkov[fSplineNumber] = 0; 1525 fIntegralCerenkov[fSplineNumber] = 0;
1502 fIntegralCerenkov[0] = 0; 1526 fIntegralCerenkov[0] = 0;
1503 k = fIntervalNumber -1; 1527 k = fIntervalNumber -1;
1504 1528
1505 for( i = fSplineNumber-1; i >= 1; i-- ) 1529 for( i = fSplineNumber-1; i >= 1; i-- )
1506 { 1530 {
1507 if(fSplineEnergy[i] >= fEnergyInterval[ 1531 if(fSplineEnergy[i] >= fEnergyInterval[k])
1508 { 1532 {
1509 fIntegralCerenkov[i] = fIntegralCeren 1533 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1510 // G4cout<<"int: i = "<<i<<"; sumC = 1534 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1511 } 1535 }
1512 else 1536 else
1513 { 1537 {
1514 fIntegralCerenkov[i] = fIntegralCeren 1538 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1515 SumOverBor 1539 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1516 k--; 1540 k--;
1517 // G4cout<<"bord: i = "<<i<<"; sumC = 1541 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1518 } 1542 }
1519 } 1543 }
1520 1544
1521 } // end of IntegralCerenkov 1545 } // end of IntegralCerenkov
1522 1546
1523 ///////////////////////////////////////////// 1547 ////////////////////////////////////////////////////////////////////////
1524 // 1548 //
1525 // Calculation of the PAI MM-Cerenkov integra 1549 // Calculation of the PAI MM-Cerenkov integral cross-section
1526 // fIntegralMM[1] = specific MM-Cerenkov ioni 1550 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1527 // and fIntegralMM[0] = mean MM-Cerenkov loss 1551 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm
1528 1552
1529 void G4PAIxSection::IntegralMM() 1553 void G4PAIxSection::IntegralMM()
1530 { 1554 {
1531 G4int i, k; 1555 G4int i, k;
1532 fIntegralMM[fSplineNumber] = 0; 1556 fIntegralMM[fSplineNumber] = 0;
1533 fIntegralMM[0] = 0; 1557 fIntegralMM[0] = 0;
1534 k = fIntervalNumber -1; 1558 k = fIntervalNumber -1;
1535 1559
1536 for( i = fSplineNumber-1; i >= 1; i-- ) 1560 for( i = fSplineNumber-1; i >= 1; i-- )
1537 { 1561 {
1538 if(fSplineEnergy[i] >= fEnergyInterval[ 1562 if(fSplineEnergy[i] >= fEnergyInterval[k])
1539 { 1563 {
1540 fIntegralMM[i] = fIntegralMM[i+1] + S 1564 fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1541 // G4cout<<"int: i = "<<i<<"; sumC = 1565 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1542 } 1566 }
1543 else 1567 else
1544 { 1568 {
1545 fIntegralMM[i] = fIntegralMM[i+1] + 1569 fIntegralMM[i] = fIntegralMM[i+1] +
1546 SumOverBor 1570 SumOverBordMM(i+1,fEnergyInterval[k]);
1547 k--; 1571 k--;
1548 // G4cout<<"bord: i = "<<i<<"; sumC = 1572 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1549 } 1573 }
1550 } 1574 }
1551 1575
1552 } // end of IntegralMM 1576 } // end of IntegralMM
1553 1577
1554 ///////////////////////////////////////////// 1578 ////////////////////////////////////////////////////////////////////////
1555 // 1579 //
1556 // Calculation of the PAI Plasmon integral cr 1580 // Calculation of the PAI Plasmon integral cross-section
1557 // fIntegralPlasmon[1] = splasmon primary ion 1581 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1558 // and fIntegralPlasmon[0] = mean plasmon los 1582 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1559 1583
1560 void G4PAIxSection::IntegralPlasmon() 1584 void G4PAIxSection::IntegralPlasmon()
1561 { 1585 {
1562 fIntegralPlasmon[fSplineNumber] = 0; 1586 fIntegralPlasmon[fSplineNumber] = 0;
1563 fIntegralPlasmon[0] = 0; 1587 fIntegralPlasmon[0] = 0;
1564 G4int k = fIntervalNumber -1; 1588 G4int k = fIntervalNumber -1;
1565 for(G4int i=fSplineNumber-1;i>=1;i--) 1589 for(G4int i=fSplineNumber-1;i>=1;i--)
1566 { 1590 {
1567 if(fSplineEnergy[i] >= fEnergyInterval[ 1591 if(fSplineEnergy[i] >= fEnergyInterval[k])
1568 { 1592 {
1569 fIntegralPlasmon[i] = fIntegralPlasmo 1593 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1570 } 1594 }
1571 else 1595 else
1572 { 1596 {
1573 fIntegralPlasmon[i] = fIntegralPlasmo 1597 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1574 SumOverBor 1598 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1575 k--; 1599 k--;
1576 } 1600 }
1577 } 1601 }
1578 1602
1579 } // end of IntegralPlasmon 1603 } // end of IntegralPlasmon
1580 1604
1581 ///////////////////////////////////////////// 1605 ////////////////////////////////////////////////////////////////////////
1582 // 1606 //
1583 // Calculation of the PAI resonance integral 1607 // Calculation of the PAI resonance integral cross-section
1584 // fIntegralResonance[1] = resonance primary 1608 // fIntegralResonance[1] = resonance primary ionisation, 1/cm
1585 // and fIntegralResonance[0] = mean resonance 1609 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm
1586 1610
1587 void G4PAIxSection::IntegralResonance() 1611 void G4PAIxSection::IntegralResonance()
1588 { 1612 {
1589 fIntegralResonance[fSplineNumber] = 0; 1613 fIntegralResonance[fSplineNumber] = 0;
1590 fIntegralResonance[0] = 0; 1614 fIntegralResonance[0] = 0;
1591 G4int k = fIntervalNumber -1; 1615 G4int k = fIntervalNumber -1;
1592 for(G4int i=fSplineNumber-1;i>=1;i--) 1616 for(G4int i=fSplineNumber-1;i>=1;i--)
1593 { 1617 {
1594 if(fSplineEnergy[i] >= fEnergyInterval[ 1618 if(fSplineEnergy[i] >= fEnergyInterval[k])
1595 { 1619 {
1596 fIntegralResonance[i] = fIntegralReso 1620 fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1597 } 1621 }
1598 else 1622 else
1599 { 1623 {
1600 fIntegralResonance[i] = fIntegralReso 1624 fIntegralResonance[i] = fIntegralResonance[i+1] +
1601 SumOverBor 1625 SumOverBordResonance(i+1,fEnergyInterval[k]);
1602 k--; 1626 k--;
1603 } 1627 }
1604 } 1628 }
1605 1629
1606 } // end of IntegralResonance 1630 } // end of IntegralResonance
1607 1631
1608 ///////////////////////////////////////////// 1632 //////////////////////////////////////////////////////////////////////
1609 // 1633 //
1610 // Calculation the PAI integral cross-section 1634 // Calculation the PAI integral cross-section inside
1611 // of interval of continuous values of photo- 1635 // of interval of continuous values of photo-ionisation
1612 // cross-section. Parameter 'i' is the numbe 1636 // cross-section. Parameter 'i' is the number of interval.
1613 1637
1614 G4double G4PAIxSection::SumOverInterval( G4in 1638 G4double G4PAIxSection::SumOverInterval( G4int i )
1615 { 1639 {
1616 G4double x0,x1,y0,yy1,a,b,c,result; 1640 G4double x0,x1,y0,yy1,a,b,c,result;
1617 1641
1618 x0 = fSplineEnergy[i]; 1642 x0 = fSplineEnergy[i];
1619 x1 = fSplineEnergy[i+1]; 1643 x1 = fSplineEnergy[i+1];
1620 if(fVerbose>0) G4cout<<"SumOverInterval i= 1644 if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl;
1621 1645
1622 if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/( 1646 if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1623 1647
1624 y0 = fDifPAIxSection[i]; 1648 y0 = fDifPAIxSection[i];
1625 yy1 = fDifPAIxSection[i+1]; 1649 yy1 = fDifPAIxSection[i+1];
1626 1650
1627 if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 1651 if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1628 1652
1629 c = x1/x0; 1653 c = x1/x0;
1630 a = log10(yy1/y0)/log10(c); 1654 a = log10(yy1/y0)/log10(c);
1631 1655
1632 if(fVerbose>0) G4cout<<"SumOverInterval, a 1656 if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1633 1657
1634 b = 0.0; << 1658 // b = log10(y0) - a*log10(x0);
1635 if(a < 20.) b = y0/pow(x0,a); << 1659 b = y0/pow(x0,a);
1636 <<
1637 a += 1.; 1660 a += 1.;
1638 if( std::abs(a) < 1.e-6 ) 1661 if( std::abs(a) < 1.e-6 )
1639 { 1662 {
1640 result = b*log(x1/x0); 1663 result = b*log(x1/x0);
1641 } 1664 }
1642 else 1665 else
1643 { 1666 {
1644 result = y0*(x1*pow(c,a-1) - x0)/a; 1667 result = y0*(x1*pow(c,a-1) - x0)/a;
1645 } 1668 }
1646 a += 1.; 1669 a += 1.;
1647 if( std::abs(a) < 1.e-6 ) 1670 if( std::abs(a) < 1.e-6 )
1648 { 1671 {
1649 fIntegralPAIxSection[0] += b*log(x1/x0); << 1672 fIntegralPAIxSection[0] += b*log(x1/x0);
1650 } 1673 }
1651 else 1674 else
1652 { 1675 {
1653 fIntegralPAIxSection[0] += y0*(x1*x1*po 1676 fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1654 } 1677 }
1655 if(fVerbose>0) G4cout<<"SumOverInterval, r 1678 if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1656 return result; 1679 return result;
1657 1680
1658 } // end of SumOverInterval 1681 } // end of SumOverInterval
1659 1682
1660 ///////////////////////////////// 1683 /////////////////////////////////
1661 1684
1662 G4double G4PAIxSection::SumOverIntervaldEdx( 1685 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1663 { 1686 {
1664 G4double x0,x1,y0,yy1,a,b,c,result; 1687 G4double x0,x1,y0,yy1,a,b,c,result;
1665 1688
1666 x0 = fSplineEnergy[i]; 1689 x0 = fSplineEnergy[i];
1667 x1 = fSplineEnergy[i+1]; 1690 x1 = fSplineEnergy[i+1];
1668 1691
1669 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1692 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1670 1693
1671 y0 = fDifPAIxSection[i]; 1694 y0 = fDifPAIxSection[i];
1672 yy1 = fDifPAIxSection[i+1]; 1695 yy1 = fDifPAIxSection[i+1];
1673 c = x1/x0; 1696 c = x1/x0;
1674 a = log10(yy1/y0)/log10(c); 1697 a = log10(yy1/y0)/log10(c);
1675 << 1698 // b = log10(y0) - a*log10(x0);
1676 b = 0.0; << 1699 b = y0/pow(x0,a);
1677 if(a < 20.) b = y0/pow(x0,a); <<
1678 <<
1679 a += 2; 1700 a += 2;
1680 if(a == 0) 1701 if(a == 0)
1681 { 1702 {
1682 result = b*log(x1/x0); 1703 result = b*log(x1/x0);
1683 } 1704 }
1684 else 1705 else
1685 { 1706 {
1686 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a 1707 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1687 } 1708 }
1688 return result; 1709 return result;
1689 1710
1690 } // end of SumOverInterval 1711 } // end of SumOverInterval
1691 1712
1692 ///////////////////////////////////////////// 1713 //////////////////////////////////////////////////////////////////////
1693 // 1714 //
1694 // Calculation the PAI Cerenkov integral cros 1715 // Calculation the PAI Cerenkov integral cross-section inside
1695 // of interval of continuous values of photo- 1716 // of interval of continuous values of photo-ionisation Cerenkov
1696 // cross-section. Parameter 'i' is the numbe 1717 // cross-section. Parameter 'i' is the number of interval.
1697 1718
1698 G4double G4PAIxSection::SumOverInterCerenkov( 1719 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1699 { 1720 {
1700 G4double x0,x1,y0,yy1,a,b,c,result; 1721 G4double x0,x1,y0,yy1,a,b,c,result;
1701 1722
1702 x0 = fSplineEnergy[i]; 1723 x0 = fSplineEnergy[i];
1703 x1 = fSplineEnergy[i+1]; 1724 x1 = fSplineEnergy[i+1];
1704 1725
1705 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1726 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1706 1727
1707 y0 = fdNdxCerenkov[i]; 1728 y0 = fdNdxCerenkov[i];
1708 yy1 = fdNdxCerenkov[i+1]; 1729 yy1 = fdNdxCerenkov[i+1];
1709 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<< 1730 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1710 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4e 1731 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1711 1732
1712 c = x1/x0; 1733 c = x1/x0;
1713 a = log10(yy1/y0)/log10(c); 1734 a = log10(yy1/y0)/log10(c);
1714 << 1735 b = y0/pow(x0,a);
1715 if(a > 20.0) b = 0.0; <<
1716 else b = y0/pow(x0,a); <<
1717 1736
1718 a += 1.0; 1737 a += 1.0;
1719 if(a == 0) result = b*log(c); 1738 if(a == 0) result = b*log(c);
1720 else result = y0*(x1*pow(c,a-1) - x0 1739 else result = y0*(x1*pow(c,a-1) - x0)/a;
1721 a += 1.0; 1740 a += 1.0;
1722 1741
1723 if( a == 0 ) fIntegralCerenkov[0] += b*log 1742 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1724 else fIntegralCerenkov[0] += y0*(x 1743 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1725 // G4cout<<"a = "<<a<<"; b = "<<b<<"; res 1744 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1726 return result; 1745 return result;
1727 1746
1728 } // end of SumOverInterCerenkov 1747 } // end of SumOverInterCerenkov
1729 1748
1730 ///////////////////////////////////////////// 1749 //////////////////////////////////////////////////////////////////////
1731 // 1750 //
1732 // Calculation the PAI MM-Cerenkov integral c 1751 // Calculation the PAI MM-Cerenkov integral cross-section inside
1733 // of interval of continuous values of photo- 1752 // of interval of continuous values of photo-ionisation Cerenkov
1734 // cross-section. Parameter 'i' is the numbe 1753 // cross-section. Parameter 'i' is the number of interval.
1735 1754
1736 G4double G4PAIxSection::SumOverInterMM( G4int 1755 G4double G4PAIxSection::SumOverInterMM( G4int i )
1737 { 1756 {
1738 G4double x0,x1,y0,yy1,a,b,c,result; 1757 G4double x0,x1,y0,yy1,a,b,c,result;
1739 1758
1740 x0 = fSplineEnergy[i]; 1759 x0 = fSplineEnergy[i];
1741 x1 = fSplineEnergy[i+1]; 1760 x1 = fSplineEnergy[i+1];
1742 1761
1743 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1762 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1744 1763
1745 y0 = fdNdxMM[i]; 1764 y0 = fdNdxMM[i];
1746 yy1 = fdNdxMM[i+1]; 1765 yy1 = fdNdxMM[i+1];
1747 //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<" 1766 //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1748 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4e 1767 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1749 1768
1750 c = x1/x0; 1769 c = x1/x0;
1751 //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/ 1770 //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl;
1752 a = log10(yy1/y0)/log10(c); 1771 a = log10(yy1/y0)/log10(c);
1753 << 1772 if(a > 10.0) return 0.;
1754 b = 0.0; << 1773 b = y0/pow(x0,a);
1755 if(a < 20.) b = y0/pow(x0,a); <<
1756 1774
1757 a += 1.0; 1775 a += 1.0;
1758 if(a == 0) result = b*log(c); 1776 if(a == 0) result = b*log(c);
1759 else result = y0*(x1*pow(c,a-1) - x0 1777 else result = y0*(x1*pow(c,a-1) - x0)/a;
1760 a += 1.0; 1778 a += 1.0;
1761 1779
1762 if( a == 0 ) fIntegralMM[0] += b*log(c); 1780 if( a == 0 ) fIntegralMM[0] += b*log(c);
1763 else fIntegralMM[0] += y0*(x1*x1*p 1781 else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1764 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul 1782 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1765 return result; 1783 return result;
1766 1784
1767 } // end of SumOverInterMM 1785 } // end of SumOverInterMM
1768 1786
1769 ///////////////////////////////////////////// 1787 //////////////////////////////////////////////////////////////////////
1770 // 1788 //
1771 // Calculation the PAI Plasmon integral cross 1789 // Calculation the PAI Plasmon integral cross-section inside
1772 // of interval of continuous values of photo- 1790 // of interval of continuous values of photo-ionisation Plasmon
1773 // cross-section. Parameter 'i' is the numbe 1791 // cross-section. Parameter 'i' is the number of interval.
1774 1792
1775 G4double G4PAIxSection::SumOverInterPlasmon( 1793 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1776 { 1794 {
1777 G4double x0,x1,y0,yy1,a,b,c,result; 1795 G4double x0,x1,y0,yy1,a,b,c,result;
1778 1796
1779 x0 = fSplineEnergy[i]; 1797 x0 = fSplineEnergy[i];
1780 x1 = fSplineEnergy[i+1]; 1798 x1 = fSplineEnergy[i+1];
1781 1799
1782 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1800 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1783 1801
1784 y0 = fdNdxPlasmon[i]; 1802 y0 = fdNdxPlasmon[i];
1785 yy1 = fdNdxPlasmon[i+1]; 1803 yy1 = fdNdxPlasmon[i+1];
1786 c = x1/x0; << 1804 c =x1/x0;
1787 a = log10(yy1/y0)/log10(c); 1805 a = log10(yy1/y0)/log10(c);
1788 << 1806 if(a > 10.0) return 0.;
1789 b = 0.0; << 1807 // b = log10(y0) - a*log10(x0);
1790 if(a < 20.) b = y0/pow(x0,a); << 1808 b = y0/pow(x0,a);
1791 1809
1792 a += 1.0; 1810 a += 1.0;
1793 if(a == 0) result = b*log(x1/x0); 1811 if(a == 0) result = b*log(x1/x0);
1794 else result = y0*(x1*pow(c,a-1) - x0 1812 else result = y0*(x1*pow(c,a-1) - x0)/a;
1795 a += 1.0; 1813 a += 1.0;
1796 1814
1797 if( a == 0 ) fIntegralPlasmon[0] += b*log( 1815 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1798 else fIntegralPlasmon[0] += y0*(x1 1816 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1799 1817
1800 return result; 1818 return result;
1801 1819
1802 } // end of SumOverInterPlasmon 1820 } // end of SumOverInterPlasmon
1803 1821
1804 ///////////////////////////////////////////// 1822 //////////////////////////////////////////////////////////////////////
1805 // 1823 //
1806 // Calculation the PAI resonance integral cro 1824 // Calculation the PAI resonance integral cross-section inside
1807 // of interval of continuous values of photo- 1825 // of interval of continuous values of photo-ionisation resonance
1808 // cross-section. Parameter 'i' is the numbe 1826 // cross-section. Parameter 'i' is the number of interval.
1809 1827
1810 G4double G4PAIxSection::SumOverInterResonance 1828 G4double G4PAIxSection::SumOverInterResonance( G4int i )
1811 { 1829 {
1812 G4double x0,x1,y0,yy1,a,b,c,result; 1830 G4double x0,x1,y0,yy1,a,b,c,result;
1813 1831
1814 x0 = fSplineEnergy[i]; 1832 x0 = fSplineEnergy[i];
1815 x1 = fSplineEnergy[i+1]; 1833 x1 = fSplineEnergy[i+1];
1816 1834
1817 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1835 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1818 1836
1819 y0 = fdNdxResonance[i]; 1837 y0 = fdNdxResonance[i];
1820 yy1 = fdNdxResonance[i+1]; 1838 yy1 = fdNdxResonance[i+1];
1821 c =x1/x0; 1839 c =x1/x0;
1822 a = log10(yy1/y0)/log10(c); 1840 a = log10(yy1/y0)/log10(c);
1823 << 1841 if(a > 10.0) return 0.;
1824 b = 0.0; << 1842 // b = log10(y0) - a*log10(x0);
1825 if(a < 20.) b = y0/pow(x0,a); << 1843 b = y0/pow(x0,a);
1826 1844
1827 a += 1.0; 1845 a += 1.0;
1828 if(a == 0) result = b*log(x1/x0); 1846 if(a == 0) result = b*log(x1/x0);
1829 else result = y0*(x1*pow(c,a-1) - x0 1847 else result = y0*(x1*pow(c,a-1) - x0)/a;
1830 a += 1.0; 1848 a += 1.0;
1831 1849
1832 if( a == 0 ) fIntegralResonance[0] += b*lo 1850 if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1833 else fIntegralResonance[0] += y0*( 1851 else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1834 1852
1835 return result; 1853 return result;
1836 1854
1837 } // end of SumOverInterResonance 1855 } // end of SumOverInterResonance
1838 1856
1839 ///////////////////////////////////////////// 1857 ///////////////////////////////////////////////////////////////////////////////
1840 // 1858 //
1841 // Integration of PAI cross-section for the c 1859 // Integration of PAI cross-section for the case of
1842 // passing across border between intervals 1860 // passing across border between intervals
1843 1861
1844 G4double G4PAIxSection::SumOverBorder( G4int 1862 G4double G4PAIxSection::SumOverBorder( G4int i ,
1845 G4doub 1863 G4double en0 )
1846 { 1864 {
1847 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result 1865 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1848 1866
1849 e0 = en0; 1867 e0 = en0;
1850 x0 = fSplineEnergy[i]; 1868 x0 = fSplineEnergy[i];
1851 x1 = fSplineEnergy[i+1]; 1869 x1 = fSplineEnergy[i+1];
1852 y0 = fDifPAIxSection[i]; 1870 y0 = fDifPAIxSection[i];
1853 yy1 = fDifPAIxSection[i+1]; 1871 yy1 = fDifPAIxSection[i+1];
1854 1872
1855 //c = x1/x0; 1873 //c = x1/x0;
1856 d = e0/x0; 1874 d = e0/x0;
1857 a = log10(yy1/y0)/log10(x1/x0); 1875 a = log10(yy1/y0)/log10(x1/x0);
>> 1876 if(a > 10.0) return 0.;
1858 1877
1859 if(fVerbose>0) G4cout<<"SumOverBorder, a = 1878 if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1860 1879
1861 b = 0.0; << 1880 // b0 = log10(y0) - a*log10(x0);
1862 if(a < 20.) b = y0/pow(x0,a); << 1881 b = y0/pow(x0,a); // pow(10.,b);
1863 1882
1864 a += 1.; 1883 a += 1.;
1865 if( std::abs(a) < 1.e-6 ) 1884 if( std::abs(a) < 1.e-6 )
1866 { 1885 {
1867 result = b*log(x0/e0); 1886 result = b*log(x0/e0);
1868 } 1887 }
1869 else 1888 else
1870 { 1889 {
1871 result = y0*(x0 - e0*pow(d,a-1))/a; 1890 result = y0*(x0 - e0*pow(d,a-1))/a;
1872 } 1891 }
1873 a += 1.; 1892 a += 1.;
1874 if( std::abs(a) < 1.e-6 ) 1893 if( std::abs(a) < 1.e-6 )
1875 { 1894 {
1876 fIntegralPAIxSection[0] += b*log(x0/e0) 1895 fIntegralPAIxSection[0] += b*log(x0/e0);
1877 } 1896 }
1878 else 1897 else
1879 { 1898 {
1880 fIntegralPAIxSection[0] += y0*(x0*x0 - 1899 fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1881 } 1900 }
1882 x0 = fSplineEnergy[i - 1]; 1901 x0 = fSplineEnergy[i - 1];
1883 x1 = fSplineEnergy[i - 2]; 1902 x1 = fSplineEnergy[i - 2];
1884 y0 = fDifPAIxSection[i - 1]; 1903 y0 = fDifPAIxSection[i - 1];
1885 yy1 = fDifPAIxSection[i - 2]; 1904 yy1 = fDifPAIxSection[i - 2];
1886 1905
>> 1906 //c = x1/x0;
1887 d = e0/x0; 1907 d = e0/x0;
1888 a = log10(yy1/y0)/log10(x1/x0); 1908 a = log10(yy1/y0)/log10(x1/x0);
1889 << 1909 // b0 = log10(y0) - a*log10(x0);
1890 b = 0.0; << 1910 b = y0/pow(x0,a);
1891 if(a < 20.) b = y0/pow(x0,a); <<
1892 <<
1893 a += 1.; 1911 a += 1.;
1894 if( std::abs(a) < 1.e-6 ) 1912 if( std::abs(a) < 1.e-6 )
1895 { 1913 {
1896 result += b*log(e0/x0); 1914 result += b*log(e0/x0);
1897 } 1915 }
1898 else 1916 else
1899 { 1917 {
1900 result += y0*(e0*pow(d,a-1) - x0)/a; 1918 result += y0*(e0*pow(d,a-1) - x0)/a;
1901 } 1919 }
1902 a += 1.; 1920 a += 1.;
1903 if( std::abs(a) < 1.e-6 ) 1921 if( std::abs(a) < 1.e-6 )
1904 { 1922 {
1905 fIntegralPAIxSection[0] += b*log(e0/x0) 1923 fIntegralPAIxSection[0] += b*log(e0/x0);
1906 } 1924 }
1907 else 1925 else
1908 { 1926 {
1909 fIntegralPAIxSection[0] += y0*(e0*e0*po 1927 fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1910 } 1928 }
1911 return result; 1929 return result;
1912 1930
1913 } 1931 }
1914 1932
1915 ///////////////////////////////////////////// 1933 ///////////////////////////////////////////////////////////////////////
1916 1934
1917 G4double G4PAIxSection::SumOverBorderdEdx( G4 << 1935 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
>> 1936 G4double en0 )
1918 { 1937 {
1919 G4double x0,x1,y0,yy1,a,b,d,e0,result; << 1938 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1920 1939
1921 e0 = en0; 1940 e0 = en0;
1922 x0 = fSplineEnergy[i]; 1941 x0 = fSplineEnergy[i];
1923 x1 = fSplineEnergy[i+1]; 1942 x1 = fSplineEnergy[i+1];
1924 y0 = fDifPAIxSection[i]; 1943 y0 = fDifPAIxSection[i];
1925 yy1 = fDifPAIxSection[i+1]; 1944 yy1 = fDifPAIxSection[i+1];
1926 1945
>> 1946 //c = x1/x0;
1927 d = e0/x0; 1947 d = e0/x0;
1928 a = log10(yy1/y0)/log10(x1/x0); 1948 a = log10(yy1/y0)/log10(x1/x0);
1929 << 1949 if(a > 10.0) return 0.;
1930 b = 0.0; << 1950 // b0 = log10(y0) - a*log10(x0);
1931 if(a < 20.) b = y0/pow(x0,a); << 1951 b = y0/pow(x0,a); // pow(10.,b);
1932 1952
1933 a += 2; 1953 a += 2;
1934 if(a == 0) 1954 if(a == 0)
1935 { 1955 {
1936 result = b*log(x0/e0); 1956 result = b*log(x0/e0);
1937 } 1957 }
1938 else 1958 else
1939 { 1959 {
1940 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/ 1960 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1941 } 1961 }
1942 x0 = fSplineEnergy[i - 1]; 1962 x0 = fSplineEnergy[i - 1];
1943 x1 = fSplineEnergy[i - 2]; 1963 x1 = fSplineEnergy[i - 2];
1944 y0 = fDifPAIxSection[i - 1]; 1964 y0 = fDifPAIxSection[i - 1];
1945 yy1 = fDifPAIxSection[i - 2]; 1965 yy1 = fDifPAIxSection[i - 2];
1946 1966
1947 // c = x1/x0; 1967 // c = x1/x0;
1948 d = e0/x0; 1968 d = e0/x0;
1949 a = log10(yy1/y0)/log10(x1/x0); 1969 a = log10(yy1/y0)/log10(x1/x0);
1950 << 1970 // b0 = log10(y0) - a*log10(x0);
1951 b = 0.0; << 1971 b = y0/pow(x0,a);
1952 if(a < 20.) b = y0/pow(x0,a); <<
1953 <<
1954 a += 2; 1972 a += 2;
1955 if(a == 0) 1973 if(a == 0)
1956 { 1974 {
1957 result += b*log(e0/x0); 1975 result += b*log(e0/x0);
1958 } 1976 }
1959 else 1977 else
1960 { 1978 {
1961 result += y0*(e0*e0*pow(d,a-2) - x0*x0) 1979 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1962 } 1980 }
1963 return result; 1981 return result;
1964 1982
1965 } 1983 }
1966 1984
1967 ///////////////////////////////////////////// 1985 ///////////////////////////////////////////////////////////////////////////////
1968 // 1986 //
1969 // Integration of Cerenkov cross-section for 1987 // Integration of Cerenkov cross-section for the case of
1970 // passing across border between intervals 1988 // passing across border between intervals
1971 1989
1972 G4double G4PAIxSection::SumOverBordCerenkov( << 1990 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
>> 1991 G4double en0 )
1973 { 1992 {
1974 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 1993 G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1975 1994
1976 e0 = en0; 1995 e0 = en0;
1977 x0 = fSplineEnergy[i]; 1996 x0 = fSplineEnergy[i];
1978 x1 = fSplineEnergy[i+1]; 1997 x1 = fSplineEnergy[i+1];
1979 y0 = fdNdxCerenkov[i]; 1998 y0 = fdNdxCerenkov[i];
1980 yy1 = fdNdxCerenkov[i+1]; 1999 yy1 = fdNdxCerenkov[i+1];
1981 2000
1982 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "< << 2001 // G4cout<<G4endl;
1983 //<<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl << 2002 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
>> 2003 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1984 c = x1/x0; 2004 c = x1/x0;
1985 d = e0/x0; 2005 d = e0/x0;
1986 a = log10(yy1/y0)/log10(c); 2006 a = log10(yy1/y0)/log10(c);
1987 //G4cout << " a= " << a << " c=" << c < << 2007 if(a > 10.0) return 0.;
1988 << 2008 // b0 = log10(y0) - a*log10(x0);
1989 b = 0.0; << 2009 b = y0/pow(x0,a); // pow(10.,b0);
1990 if(a < 20.) b = y0/pow(x0,a); <<
1991 2010
1992 a += 1.0; 2011 a += 1.0;
1993 if( a == 0 ) result = b*log(x0/e0); 2012 if( a == 0 ) result = b*log(x0/e0);
1994 else result = y0*(x0 - e0*pow(d,a- 2013 else result = y0*(x0 - e0*pow(d,a-1))/a;
1995 a += 1.0; 2014 a += 1.0;
1996 2015
1997 if( a == 0 ) fIntegralCerenkov[0] += b*log 2016 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1998 else fIntegralCerenkov[0] += y0*(x 2017 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1999 2018
>> 2019 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
>> 2020
2000 x0 = fSplineEnergy[i - 1]; 2021 x0 = fSplineEnergy[i - 1];
2001 x1 = fSplineEnergy[i - 2]; 2022 x1 = fSplineEnergy[i - 2];
2002 y0 = fdNdxCerenkov[i - 1]; 2023 y0 = fdNdxCerenkov[i - 1];
2003 yy1 = fdNdxCerenkov[i - 2]; 2024 yy1 = fdNdxCerenkov[i - 2];
2004 2025
2005 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2026 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2006 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 2027 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2007 2028
2008 c = x1/x0; 2029 c = x1/x0;
2009 d = e0/x0; 2030 d = e0/x0;
2010 a = log10(yy1/y0)/log10(c); << 2031 a = log10(yy1/y0)/log10(x1/x0);
2011 << 2032 // b0 = log10(y0) - a*log10(x0);
2012 b = 0.0; << 2033 b = y0/pow(x0,a); // pow(10.,b0);
2013 if(a < 20.) b = y0/pow(x0,a); <<
2014 2034
2015 a += 1.0; 2035 a += 1.0;
2016 if( a == 0 ) result += b*log(e0/x0); 2036 if( a == 0 ) result += b*log(e0/x0);
2017 else result += y0*(e0*pow(d,a-1) - 2037 else result += y0*(e0*pow(d,a-1) - x0 )/a;
2018 a += 1.0; 2038 a += 1.0;
2019 2039
2020 if( a == 0 ) fIntegralCerenkov[0] += b*log << 2040 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2021 else fIntegralCerenkov[0] += y0*(e << 2041 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
>> 2042
>> 2043 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
>> 2044 // <<b<<"; result = "<<result<<G4endl;
2022 2045
2023 //G4cout<<" a="<< a <<" b="<< b <<" res <<
2024 return result; 2046 return result;
>> 2047
2025 } 2048 }
2026 2049
2027 ///////////////////////////////////////////// 2050 ///////////////////////////////////////////////////////////////////////////////
2028 // 2051 //
2029 // Integration of MM-Cerenkov cross-section f 2052 // Integration of MM-Cerenkov cross-section for the case of
2030 // passing across border between intervals 2053 // passing across border between intervals
2031 2054
2032 G4double G4PAIxSection::SumOverBordMM( G4int << 2055 G4double G4PAIxSection::SumOverBordMM( G4int i ,
>> 2056 G4double en0 )
2033 { 2057 {
2034 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 2058 G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2035 2059
2036 e0 = en0; 2060 e0 = en0;
2037 x0 = fSplineEnergy[i]; 2061 x0 = fSplineEnergy[i];
2038 x1 = fSplineEnergy[i+1]; 2062 x1 = fSplineEnergy[i+1];
2039 y0 = fdNdxMM[i]; 2063 y0 = fdNdxMM[i];
2040 yy1 = fdNdxMM[i+1]; 2064 yy1 = fdNdxMM[i+1];
2041 2065
>> 2066 // G4cout<<G4endl;
2042 // G4cout<<"SumBordC, i = "<<i<<"; en0 = 2067 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2043 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G 2068 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2044 c = x1/x0; 2069 c = x1/x0;
2045 d = e0/x0; 2070 d = e0/x0;
2046 a = log10(yy1/y0)/log10(c); 2071 a = log10(yy1/y0)/log10(c);
2047 << 2072 if(a > 10.0) return 0.;
2048 if(a > 20.0) b = 0.0; << 2073 // b0 = log10(y0) - a*log10(x0);
2049 else b = y0/pow(x0,a); << 2074 b = y0/pow(x0,a); // pow(10.,b0);
2050 2075
2051 a += 1.0; 2076 a += 1.0;
2052 if( a == 0 ) result = b*log(x0/e0); 2077 if( a == 0 ) result = b*log(x0/e0);
2053 else result = y0*(x0 - e0*pow(d,a- 2078 else result = y0*(x0 - e0*pow(d,a-1))/a;
2054 a += 1.0; 2079 a += 1.0;
2055 2080
2056 if( a == 0 ) fIntegralMM[0] += b*log(x0/e0 2081 if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2057 else fIntegralMM[0] += y0*(x0*x0 - 2082 else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2058 2083
2059 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b << 2084 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2060 2085
2061 x0 = fSplineEnergy[i - 1]; 2086 x0 = fSplineEnergy[i - 1];
2062 x1 = fSplineEnergy[i - 2]; 2087 x1 = fSplineEnergy[i - 2];
2063 y0 = fdNdxMM[i - 1]; 2088 y0 = fdNdxMM[i - 1];
2064 yy1 = fdNdxMM[i - 2]; 2089 yy1 = fdNdxMM[i - 2];
2065 2090
2066 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2091 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2067 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 2092 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2068 2093
2069 c = x1/x0; 2094 c = x1/x0;
2070 d = e0/x0; 2095 d = e0/x0;
2071 a = log10(yy1/y0)/log10(x1/x0); 2096 a = log10(yy1/y0)/log10(x1/x0);
2072 << 2097 // b0 = log10(y0) - a*log10(x0);
2073 if(a > 20.0) b = 0.0; << 2098 b = y0/pow(x0,a); // pow(10.,b0);
2074 else b = y0/pow(x0,a); <<
2075 2099
2076 a += 1.0; 2100 a += 1.0;
2077 if( a == 0 ) result += b*log(e0/x0); 2101 if( a == 0 ) result += b*log(e0/x0);
2078 else result += y0*(e0*pow(d,a-1) - 2102 else result += y0*(e0*pow(d,a-1) - x0 )/a;
2079 a += 1.0; 2103 a += 1.0;
2080 2104
2081 if( a == 0 ) fIntegralMM[0] += b*log(e0/x0 << 2105 if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2082 else fIntegralMM[0] += y0*(e0*e0*p << 2106 else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2083 2107
2084 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b 2108 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2085 // <<b<<"; result = "<<result<<G4endl; 2109 // <<b<<"; result = "<<result<<G4endl;
2086 2110
2087 return result; 2111 return result;
2088 2112
2089 } 2113 }
2090 2114
2091 ///////////////////////////////////////////// 2115 ///////////////////////////////////////////////////////////////////////////////
2092 // 2116 //
2093 // Integration of Plasmon cross-section for t 2117 // Integration of Plasmon cross-section for the case of
2094 // passing across border between intervals 2118 // passing across border between intervals
2095 2119
2096 G4double G4PAIxSection::SumOverBordPlasmon( G 2120 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
2097 2121 G4double en0 )
2098 { 2122 {
2099 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2123 G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2100 2124
2101 e0 = en0; 2125 e0 = en0;
2102 x0 = fSplineEnergy[i]; 2126 x0 = fSplineEnergy[i];
2103 x1 = fSplineEnergy[i+1]; 2127 x1 = fSplineEnergy[i+1];
2104 y0 = fdNdxPlasmon[i]; 2128 y0 = fdNdxPlasmon[i];
2105 yy1 = fdNdxPlasmon[i+1]; 2129 yy1 = fdNdxPlasmon[i+1];
2106 2130
2107 c = x1/x0; 2131 c = x1/x0;
2108 d = e0/x0; 2132 d = e0/x0;
2109 a = log10(yy1/y0)/log10(c); 2133 a = log10(yy1/y0)/log10(c);
2110 << 2134 if(a > 10.0) return 0.;
2111 if(a > 20.0) b = 0.0; << 2135 // b0 = log10(y0) - a*log10(x0);
2112 else b = y0/pow(x0,a); << 2136 b = y0/pow(x0,a); //pow(10.,b);
2113 2137
2114 a += 1.0; 2138 a += 1.0;
2115 if( a == 0 ) result = b*log(x0/e0); 2139 if( a == 0 ) result = b*log(x0/e0);
2116 else result = y0*(x0 - e0*pow(d,a- 2140 else result = y0*(x0 - e0*pow(d,a-1))/a;
2117 a += 1.0; 2141 a += 1.0;
2118 2142
2119 if( a == 0 ) fIntegralPlasmon[0] += b*log( 2143 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2120 else fIntegralPlasmon[0] += y0*(x0 2144 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2121 2145
2122 x0 = fSplineEnergy[i - 1]; 2146 x0 = fSplineEnergy[i - 1];
2123 x1 = fSplineEnergy[i - 2]; 2147 x1 = fSplineEnergy[i - 2];
2124 y0 = fdNdxPlasmon[i - 1]; 2148 y0 = fdNdxPlasmon[i - 1];
2125 yy1 = fdNdxPlasmon[i - 2]; 2149 yy1 = fdNdxPlasmon[i - 2];
2126 2150
2127 c = x1/x0; 2151 c = x1/x0;
2128 d = e0/x0; 2152 d = e0/x0;
2129 a = log10(yy1/y0)/log10(c); 2153 a = log10(yy1/y0)/log10(c);
2130 << 2154 // b0 = log10(y0) - a*log10(x0);
2131 if(a > 20.0) b = 0.0; << 2155 b = y0/pow(x0,a);// pow(10.,b0);
2132 else b = y0/pow(x0,a); <<
2133 2156
2134 a += 1.0; 2157 a += 1.0;
2135 if( a == 0 ) result += b*log(e0/x0); 2158 if( a == 0 ) result += b*log(e0/x0);
2136 else result += y0*(e0*pow(d,a-1) - 2159 else result += y0*(e0*pow(d,a-1) - x0)/a;
2137 a += 1.0; 2160 a += 1.0;
2138 2161
2139 if( a == 0 ) fIntegralPlasmon[0] += b*lo 2162 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2140 else fIntegralPlasmon[0] += y0*( 2163 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2141 2164
2142 return result; 2165 return result;
>> 2166
2143 } 2167 }
2144 2168
2145 ///////////////////////////////////////////// 2169 ///////////////////////////////////////////////////////////////////////////////
2146 // 2170 //
2147 // Integration of resonance cross-section for 2171 // Integration of resonance cross-section for the case of
2148 // passing across border between intervals 2172 // passing across border between intervals
2149 2173
2150 G4double G4PAIxSection::SumOverBordResonance( 2174 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
2151 2175 G4double en0 )
2152 { 2176 {
2153 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2177 G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2154 2178
2155 e0 = en0; 2179 e0 = en0;
2156 x0 = fSplineEnergy[i]; 2180 x0 = fSplineEnergy[i];
2157 x1 = fSplineEnergy[i+1]; 2181 x1 = fSplineEnergy[i+1];
2158 y0 = fdNdxResonance[i]; 2182 y0 = fdNdxResonance[i];
2159 yy1 = fdNdxResonance[i+1]; 2183 yy1 = fdNdxResonance[i+1];
2160 2184
2161 c = x1/x0; 2185 c = x1/x0;
2162 d = e0/x0; 2186 d = e0/x0;
2163 a = log10(yy1/y0)/log10(c); 2187 a = log10(yy1/y0)/log10(c);
2164 << 2188 if(a > 10.0) return 0.;
2165 if(a > 20.0) b = 0.0; << 2189 // b0 = log10(y0) - a*log10(x0);
2166 else b = y0/pow(x0,a); << 2190 b = y0/pow(x0,a); //pow(10.,b);
2167 2191
2168 a += 1.0; 2192 a += 1.0;
2169 if( a == 0 ) result = b*log(x0/e0); 2193 if( a == 0 ) result = b*log(x0/e0);
2170 else result = y0*(x0 - e0*pow(d,a- 2194 else result = y0*(x0 - e0*pow(d,a-1))/a;
2171 a += 1.0; 2195 a += 1.0;
2172 2196
2173 if( a == 0 ) fIntegralResonance[0] += b*lo 2197 if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2174 else fIntegralResonance[0] += y0*( 2198 else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2175 2199
2176 x0 = fSplineEnergy[i - 1]; 2200 x0 = fSplineEnergy[i - 1];
2177 x1 = fSplineEnergy[i - 2]; 2201 x1 = fSplineEnergy[i - 2];
2178 y0 = fdNdxResonance[i - 1]; 2202 y0 = fdNdxResonance[i - 1];
2179 yy1 = fdNdxResonance[i - 2]; 2203 yy1 = fdNdxResonance[i - 2];
2180 2204
2181 c = x1/x0; 2205 c = x1/x0;
2182 d = e0/x0; 2206 d = e0/x0;
2183 a = log10(yy1/y0)/log10(c); 2207 a = log10(yy1/y0)/log10(c);
2184 << 2208 // b0 = log10(y0) - a*log10(x0);
2185 if(a > 20.0) b = 0.0; << 2209 b = y0/pow(x0,a);// pow(10.,b0);
2186 else b = y0/pow(x0,a); <<
2187 2210
2188 a += 1.0; 2211 a += 1.0;
2189 if( a == 0 ) result += b*log(e0/x0); 2212 if( a == 0 ) result += b*log(e0/x0);
2190 else result += y0*(e0*pow(d,a-1) - 2213 else result += y0*(e0*pow(d,a-1) - x0)/a;
2191 a += 1.0; 2214 a += 1.0;
2192 2215
2193 if( a == 0 ) fIntegralResonance[0] += b*lo << 2216 if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2194 else fIntegralResonance[0] += y0*( << 2217 else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2195 2218
2196 return result; 2219 return result;
2197 2220
2198 } 2221 }
2199 2222
2200 ///////////////////////////////////////////// 2223 /////////////////////////////////////////////////////////////////////////
2201 // 2224 //
2202 // Returns random PAI-total energy loss over 2225 // Returns random PAI-total energy loss over step
2203 2226
2204 G4double G4PAIxSection::GetStepEnergyLoss( G4 2227 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
2205 { 2228 {
2206 G4long numOfCollisions; 2229 G4long numOfCollisions;
2207 G4double meanNumber, loss = 0.0; 2230 G4double meanNumber, loss = 0.0;
2208 2231
2209 // G4cout<<" G4PAIxSection::GetStepEnergyLo 2232 // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2210 2233
2211 meanNumber = fIntegralPAIxSection[1]*step; 2234 meanNumber = fIntegralPAIxSection[1]*step;
2212 numOfCollisions = G4Poisson(meanNumber); 2235 numOfCollisions = G4Poisson(meanNumber);
2213 2236
2214 // G4cout<<"numOfCollisions = "<<numOfColli << 2237 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2215 2238
2216 while(numOfCollisions) 2239 while(numOfCollisions)
2217 { 2240 {
2218 loss += GetEnergyTransfer(); 2241 loss += GetEnergyTransfer();
2219 numOfCollisions--; 2242 numOfCollisions--;
2220 // Loop checking, 03-Aug-2015, Vladimir I 2243 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2221 } 2244 }
2222 // G4cout<<"PAI energy loss = "<<loss/keV<< 2245 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2223 2246
2224 return loss; 2247 return loss;
2225 } 2248 }
2226 2249
2227 ///////////////////////////////////////////// 2250 /////////////////////////////////////////////////////////////////////////
2228 // 2251 //
2229 // Returns random PAI-total energy transfer i 2252 // Returns random PAI-total energy transfer in one collision
2230 2253
2231 G4double G4PAIxSection::GetEnergyTransfer() 2254 G4double G4PAIxSection::GetEnergyTransfer()
2232 { 2255 {
2233 G4int iTransfer ; 2256 G4int iTransfer ;
2234 2257
2235 G4double energyTransfer, position; 2258 G4double energyTransfer, position;
2236 2259
2237 position = fIntegralPAIxSection[1]*G4Unifor 2260 position = fIntegralPAIxSection[1]*G4UniformRand();
2238 2261
2239 for( iTransfer = 1; iTransfer <= fSplineNum 2262 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2240 { 2263 {
2241 if( position >= fIntegralPAIxSection[iTra << 2264 if( position >= fIntegralPAIxSection[iTransfer] ) break;
2242 } 2265 }
2243 if(iTransfer > fSplineNumber) iTransfer--; 2266 if(iTransfer > fSplineNumber) iTransfer--;
2244 2267
2245 energyTransfer = fSplineEnergy[iTransfer]; 2268 energyTransfer = fSplineEnergy[iTransfer];
2246 2269
2247 if(iTransfer > 1) 2270 if(iTransfer > 1)
2248 { 2271 {
2249 energyTransfer -= (fSplineEnergy[iTransfe 2272 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2250 } 2273 }
2251 return energyTransfer; 2274 return energyTransfer;
2252 } 2275 }
2253 2276
2254 ///////////////////////////////////////////// 2277 /////////////////////////////////////////////////////////////////////////
2255 // 2278 //
2256 // Returns random Cerenkov energy loss over s 2279 // Returns random Cerenkov energy loss over step
2257 2280
2258 G4double G4PAIxSection::GetStepCerenkovLoss( 2281 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
2259 { 2282 {
2260 G4long numOfCollisions; 2283 G4long numOfCollisions;
2261 G4double meanNumber, loss = 0.0; 2284 G4double meanNumber, loss = 0.0;
2262 2285
2263 // G4cout<<" G4PAIxSection::GetStepCerenkov 2286 // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2264 2287
2265 meanNumber = fIntegralCerenkov[1]*step; 2288 meanNumber = fIntegralCerenkov[1]*step;
2266 numOfCollisions = G4Poisson(meanNumber); 2289 numOfCollisions = G4Poisson(meanNumber);
2267 2290
2268 // G4cout<<"numOfCollisions = "<<numOfCol 2291 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2269 2292
2270 while(numOfCollisions) 2293 while(numOfCollisions)
2271 { 2294 {
2272 loss += GetCerenkovEnergyTransfer(); 2295 loss += GetCerenkovEnergyTransfer();
2273 numOfCollisions--; 2296 numOfCollisions--;
2274 // Loop checking, 03-Aug-2015, Vladimir I 2297 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2275 } 2298 }
2276 // G4cout<<"PAI Cerenkov loss = "<<loss/keV 2299 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2277 2300
2278 return loss; 2301 return loss;
2279 } 2302 }
2280 2303
2281 ///////////////////////////////////////////// 2304 /////////////////////////////////////////////////////////////////////////
2282 // 2305 //
2283 // Returns random MM-Cerenkov energy loss ove 2306 // Returns random MM-Cerenkov energy loss over step
2284 2307
2285 G4double G4PAIxSection::GetStepMMLoss( G4doub 2308 G4double G4PAIxSection::GetStepMMLoss( G4double step )
2286 { 2309 {
2287 G4long numOfCollisions; 2310 G4long numOfCollisions;
2288 G4double meanNumber, loss = 0.0; 2311 G4double meanNumber, loss = 0.0;
2289 2312
2290 // G4cout<<" G4PAIxSection::GetStepMMLoss " 2313 // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2291 2314
2292 meanNumber = fIntegralMM[1]*step; 2315 meanNumber = fIntegralMM[1]*step;
2293 numOfCollisions = G4Poisson(meanNumber); 2316 numOfCollisions = G4Poisson(meanNumber);
2294 2317
2295 // G4cout<<"numOfCollisions = "<<numOfCol 2318 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2296 2319
2297 while(numOfCollisions) 2320 while(numOfCollisions)
2298 { 2321 {
2299 loss += GetMMEnergyTransfer(); 2322 loss += GetMMEnergyTransfer();
2300 numOfCollisions--; 2323 numOfCollisions--;
2301 // Loop checking, 03-Aug-2015, Vladimir I 2324 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2302 } 2325 }
2303 // G4cout<<"PAI MM-Cerenkov loss = "<<loss/ 2326 // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2304 2327
2305 return loss; 2328 return loss;
2306 } 2329 }
2307 2330
2308 ///////////////////////////////////////////// 2331 /////////////////////////////////////////////////////////////////////////
2309 // 2332 //
2310 // Returns Cerenkov energy transfer in one co 2333 // Returns Cerenkov energy transfer in one collision
2311 2334
2312 G4double G4PAIxSection::GetCerenkovEnergyTran 2335 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
2313 { 2336 {
2314 G4int iTransfer ; 2337 G4int iTransfer ;
2315 2338
2316 G4double energyTransfer, position; 2339 G4double energyTransfer, position;
2317 2340
2318 position = fIntegralCerenkov[1]*G4UniformRa 2341 position = fIntegralCerenkov[1]*G4UniformRand();
2319 2342
2320 for( iTransfer = 1; iTransfer <= fSplineNum 2343 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2321 { 2344 {
2322 if( position >= fIntegralCerenkov[iTr 2345 if( position >= fIntegralCerenkov[iTransfer] ) break;
2323 } 2346 }
2324 if(iTransfer > fSplineNumber) iTransfer--; 2347 if(iTransfer > fSplineNumber) iTransfer--;
2325 2348
2326 energyTransfer = fSplineEnergy[iTransfer]; 2349 energyTransfer = fSplineEnergy[iTransfer];
2327 2350
2328 if(iTransfer > 1) 2351 if(iTransfer > 1)
2329 { 2352 {
2330 energyTransfer -= (fSplineEnergy[iTransfe 2353 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2331 } 2354 }
2332 return energyTransfer; 2355 return energyTransfer;
2333 } 2356 }
2334 2357
2335 ///////////////////////////////////////////// 2358 /////////////////////////////////////////////////////////////////////////
2336 // 2359 //
2337 // Returns MM-Cerenkov energy transfer in one 2360 // Returns MM-Cerenkov energy transfer in one collision
2338 2361
2339 G4double G4PAIxSection::GetMMEnergyTransfer() 2362 G4double G4PAIxSection::GetMMEnergyTransfer()
2340 { 2363 {
2341 G4int iTransfer ; 2364 G4int iTransfer ;
2342 2365
2343 G4double energyTransfer, position; 2366 G4double energyTransfer, position;
2344 2367
2345 position = fIntegralMM[1]*G4UniformRand(); 2368 position = fIntegralMM[1]*G4UniformRand();
2346 2369
2347 for( iTransfer = 1; iTransfer <= fSplineNum 2370 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2348 { 2371 {
2349 if( position >= fIntegralMM[iTransfer] ) << 2372 if( position >= fIntegralMM[iTransfer] ) break;
2350 } 2373 }
2351 if(iTransfer > fSplineNumber) iTransfer--; 2374 if(iTransfer > fSplineNumber) iTransfer--;
2352 2375
2353 energyTransfer = fSplineEnergy[iTransfer]; 2376 energyTransfer = fSplineEnergy[iTransfer];
2354 2377
2355 if(iTransfer > 1) 2378 if(iTransfer > 1)
2356 { 2379 {
2357 energyTransfer -= (fSplineEnergy[iTransfe 2380 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2358 } 2381 }
2359 return energyTransfer; 2382 return energyTransfer;
2360 } 2383 }
2361 2384
2362 ///////////////////////////////////////////// 2385 /////////////////////////////////////////////////////////////////////////
2363 // 2386 //
2364 // Returns random plasmon energy loss over st 2387 // Returns random plasmon energy loss over step
2365 2388
2366 G4double G4PAIxSection::GetStepPlasmonLoss( G 2389 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
2367 { 2390 {
2368 G4long numOfCollisions; 2391 G4long numOfCollisions;
2369 G4double meanNumber, loss = 0.0; 2392 G4double meanNumber, loss = 0.0;
2370 2393
2371 // G4cout<<" G4PAIxSection::GetStepPlasmonL 2394 // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2372 2395
2373 meanNumber = fIntegralPlasmon[1]*step; 2396 meanNumber = fIntegralPlasmon[1]*step;
2374 numOfCollisions = G4Poisson(meanNumber); 2397 numOfCollisions = G4Poisson(meanNumber);
2375 2398
2376 // G4cout<<"numOfCollisions = "<<numOfCol 2399 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2377 2400
2378 while(numOfCollisions) 2401 while(numOfCollisions)
2379 { 2402 {
2380 loss += GetPlasmonEnergyTransfer(); 2403 loss += GetPlasmonEnergyTransfer();
2381 numOfCollisions--; 2404 numOfCollisions--;
2382 // Loop checking, 03-Aug-2015, Vladimir I 2405 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2383 } 2406 }
2384 // G4cout<<"PAI Plasmon loss = "<<loss/keV< 2407 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2385 2408
2386 return loss; 2409 return loss;
2387 } 2410 }
2388 2411
2389 ///////////////////////////////////////////// 2412 /////////////////////////////////////////////////////////////////////////
2390 // 2413 //
2391 // Returns plasmon energy transfer in one col 2414 // Returns plasmon energy transfer in one collision
2392 2415
2393 G4double G4PAIxSection::GetPlasmonEnergyTrans 2416 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
2394 { 2417 {
2395 G4int iTransfer ; 2418 G4int iTransfer ;
2396 2419
2397 G4double energyTransfer, position; 2420 G4double energyTransfer, position;
2398 2421
2399 position = fIntegralPlasmon[1]*G4UniformRan 2422 position = fIntegralPlasmon[1]*G4UniformRand();
2400 2423
2401 for( iTransfer = 1; iTransfer <= fSplineNum 2424 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2402 { 2425 {
2403 if( position >= fIntegralPlasmon[iTransfe << 2426 if( position >= fIntegralPlasmon[iTransfer] ) break;
2404 } 2427 }
2405 if(iTransfer > fSplineNumber) iTransfer--; 2428 if(iTransfer > fSplineNumber) iTransfer--;
2406 2429
2407 energyTransfer = fSplineEnergy[iTransfer]; 2430 energyTransfer = fSplineEnergy[iTransfer];
2408 2431
2409 if(iTransfer > 1) 2432 if(iTransfer > 1)
2410 { 2433 {
2411 energyTransfer -= (fSplineEnergy[iTransfe 2434 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2412 } 2435 }
2413 return energyTransfer; 2436 return energyTransfer;
2414 } 2437 }
2415 2438
2416 ///////////////////////////////////////////// 2439 /////////////////////////////////////////////////////////////////////////
2417 // 2440 //
2418 // Returns random resonance energy loss over 2441 // Returns random resonance energy loss over step
2419 2442
2420 G4double G4PAIxSection::GetStepResonanceLoss( 2443 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
2421 { 2444 {
2422 G4long numOfCollisions; 2445 G4long numOfCollisions;
2423 G4double meanNumber, loss = 0.0; 2446 G4double meanNumber, loss = 0.0;
2424 2447
2425 // G4cout<<" G4PAIxSection::GetStepCreLosnk 2448 // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2426 2449
2427 meanNumber = fIntegralResonance[1]*step; 2450 meanNumber = fIntegralResonance[1]*step;
2428 numOfCollisions = G4Poisson(meanNumber); 2451 numOfCollisions = G4Poisson(meanNumber);
2429 2452
2430 // G4cout<<"numOfCollisions = "<<numOfCol 2453 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2431 2454
2432 while(numOfCollisions) 2455 while(numOfCollisions)
2433 { 2456 {
2434 loss += GetResonanceEnergyTransfer(); 2457 loss += GetResonanceEnergyTransfer();
2435 numOfCollisions--; 2458 numOfCollisions--;
2436 // Loop checking, 03-Aug-2015, Vladimir I 2459 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2437 } 2460 }
2438 // G4cout<<"PAI resonance loss = "<<loss/ke 2461 // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2439 2462
2440 return loss; 2463 return loss;
2441 } 2464 }
2442 2465
2443 2466
2444 ///////////////////////////////////////////// 2467 /////////////////////////////////////////////////////////////////////////
2445 // 2468 //
2446 // Returns resonance energy transfer in one c 2469 // Returns resonance energy transfer in one collision
2447 2470
2448 G4double G4PAIxSection::GetResonanceEnergyTra 2471 G4double G4PAIxSection::GetResonanceEnergyTransfer()
2449 { 2472 {
2450 G4int iTransfer ; 2473 G4int iTransfer ;
2451 2474
2452 G4double energyTransfer, position; 2475 G4double energyTransfer, position;
2453 2476
2454 position = fIntegralResonance[1]*G4UniformR 2477 position = fIntegralResonance[1]*G4UniformRand();
2455 2478
2456 for( iTransfer = 1; iTransfer <= fSplineNum 2479 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2457 { 2480 {
2458 if( position >= fIntegralResonance[iTrans << 2481 if( position >= fIntegralResonance[iTransfer] ) break;
2459 } 2482 }
2460 if(iTransfer > fSplineNumber) iTransfer--; 2483 if(iTransfer > fSplineNumber) iTransfer--;
2461 2484
2462 energyTransfer = fSplineEnergy[iTransfer]; 2485 energyTransfer = fSplineEnergy[iTransfer];
2463 2486
2464 if(iTransfer > 1) 2487 if(iTransfer > 1)
2465 { 2488 {
2466 energyTransfer -= (fSplineEnergy[iTransfe 2489 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2467 } 2490 }
2468 return energyTransfer; 2491 return energyTransfer;
2469 } 2492 }
2470 2493
2471 2494
2472 ///////////////////////////////////////////// 2495 /////////////////////////////////////////////////////////////////////////
2473 // 2496 //
2474 // Returns Rutherford energy transfer in one 2497 // Returns Rutherford energy transfer in one collision
2475 2498
2476 G4double G4PAIxSection::GetRutherfordEnergyTr 2499 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
2477 { 2500 {
2478 G4int iTransfer ; 2501 G4int iTransfer ;
2479 2502
2480 G4double energyTransfer, position; 2503 G4double energyTransfer, position;
2481 2504
2482 position = (fIntegralPlasmon[1]-fIntegralRe 2505 position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2483 2506
2484 for( iTransfer = 1; iTransfer <= fSplineNum 2507 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2485 { 2508 {
2486 if( position >= (fIntegralPlasmon[iTransf << 2509 if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2487 } 2510 }
2488 if(iTransfer > fSplineNumber) iTransfer--; 2511 if(iTransfer > fSplineNumber) iTransfer--;
2489 2512
2490 energyTransfer = fSplineEnergy[iTransfer]; 2513 energyTransfer = fSplineEnergy[iTransfer];
2491 2514
2492 if(iTransfer > 1) 2515 if(iTransfer > 1)
2493 { 2516 {
2494 energyTransfer -= (fSplineEnergy[iTransfe 2517 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2495 } 2518 }
2496 return energyTransfer; 2519 return energyTransfer;
2497 } 2520 }
2498 2521
2499 ///////////////////////////////////////////// 2522 /////////////////////////////////////////////////////////////////////////////
2500 // 2523 //
2501 2524
2502 void G4PAIxSection::CallError(G4int i, const 2525 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2503 { 2526 {
2504 G4String head = "G4PAIxSection::" + methodN 2527 G4String head = "G4PAIxSection::" + methodName + "()";
2505 G4ExceptionDescription ed; 2528 G4ExceptionDescription ed;
2506 ed << "Wrong index " << i << " fSplineNumbe 2529 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2507 G4Exception(head,"pai001",FatalException,ed 2530 G4Exception(head,"pai001",FatalException,ed);
2508 } 2531 }
2509 2532
2510 ///////////////////////////////////////////// 2533 /////////////////////////////////////////////////////////////////////////////
2511 // 2534 //
2512 // Init array of Lorentz factors 2535 // Init array of Lorentz factors
2513 // 2536 //
2514 2537
2515 G4int G4PAIxSection::fNumberOfGammas = 111; 2538 G4int G4PAIxSection::fNumberOfGammas = 111;
2516 2539
2517 const G4double G4PAIxSection::fLorentzFactor[ 2540 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2518 { 2541 {
2519 0.0, 2542 0.0,
2520 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.1 2543 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2521 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.2 2544 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2522 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.4 2545 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2523 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.9 2546 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2524 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.7 2547 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2525 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.2 2548 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2526 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.2 2549 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2527 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.2 2550 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2528 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.3 2551 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2529 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.2 2552 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2530 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.9 2553 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2531 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.4 2554 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2532 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.7 2555 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2533 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.2 2556 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2534 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.8 2557 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2535 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.8 2558 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2536 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.4 2559 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2537 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.5 2560 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2538 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.2 2561 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2539 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.3 2562 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2540 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.3 2563 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2541 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.2 2564 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2542 1.065799e+05 2565 1.065799e+05
2543 }; 2566 };
2544 2567
2545 ///////////////////////////////////////////// 2568 ///////////////////////////////////////////////////////////////////////
2546 // 2569 //
2547 // The number of gamma for creation of splin 2570 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2548 // 2571 //
2549 2572
2550 const 2573 const
2551 G4int G4PAIxSection::fRefGammaNumber = 29; 2574 G4int G4PAIxSection::fRefGammaNumber = 29;
2552 2575
2553 2576
2554 // 2577 //
2555 // end of G4PAIxSection implementation file 2578 // end of G4PAIxSection implementation file
2556 // 2579 //
2557 ///////////////////////////////////////////// 2580 ////////////////////////////////////////////////////////////////////////////
2558 2581
2559 2582