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