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