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