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