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