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