| 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 // $Id: G4InitXscPAI.cc,v 1.9 2006-06-29 19:53:00 gunter Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 // 29 //
28 // G4InitXscPAI.cc -- class implementation fil 30 // G4InitXscPAI.cc -- class implementation file
29 // 31 //
30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file
31 // 33 //
32 // For information related to this code, pleas 34 // For information related to this code, please, contact
33 // the Geant4 Collaboration. 35 // the Geant4 Collaboration.
34 // 36 //
35 // R&D: Vladimir.Grichine@cern.ch 37 // R&D: Vladimir.Grichine@cern.ch
36 // 38 //
37 // History: 39 // History:
38 // 40 //
39 41
40 42
41 43
42 #include "G4InitXscPAI.hh" 44 #include "G4InitXscPAI.hh"
43 45
44 #include "globals.hh" 46 #include "globals.hh"
45 #include "G4PhysicalConstants.hh" <<
46 #include "G4SystemOfUnits.hh" <<
47 #include "G4ios.hh" 47 #include "G4ios.hh"
48 #include "G4Poisson.hh" 48 #include "G4Poisson.hh"
49 #include "G4Integrator.hh" 49 #include "G4Integrator.hh"
50 #include "G4Material.hh" 50 #include "G4Material.hh"
51 #include "G4MaterialCutsCouple.hh" 51 #include "G4MaterialCutsCouple.hh"
52 #include "G4SandiaTable.hh" 52 #include "G4SandiaTable.hh"
53 53
54 54
55 55
56 // Local class constants 56 // Local class constants
57 57
58 const G4double G4InitXscPAI::fDelta = 0 58 const G4double G4InitXscPAI::fDelta = 0.005 ; // energy shift from interval border
59 const G4int G4InitXscPAI::fPAIbin = 1 59 const G4int G4InitXscPAI::fPAIbin = 100 ; // size of energy transfer vectors
60 const G4double G4InitXscPAI::fSolidDensity = 0 60 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border
61 61
62 ////////////////////////////////////////////// 62 //////////////////////////////////////////////////////////////////
63 // 63 //
64 // Constructor 64 // Constructor
65 // 65 //
66 66
67 using namespace std; 67 using namespace std;
68 68
69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCu 69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC)
70 : fPAIxscVector(nullptr), << 70 : fPAIxscVector(NULL),
71 fPAIdEdxVector(nullptr), << 71 fPAIdEdxVector(NULL),
72 fPAIphotonVector(nullptr), << 72 fPAIphotonVector(NULL),
73 fPAIelectronVector(nullptr), << 73 fPAIelectronVector(NULL),
74 fChCosSqVector(nullptr), << 74 fChCosSqVector(NULL),
75 fChWidthVector(nullptr) << 75 fChWidthVector(NULL)
76 { 76 {
77 G4int i, j, matIndex; 77 G4int i, j, matIndex;
78 78
79 fDensity = matCC->GetMaterial()->Get 79 fDensity = matCC->GetMaterial()->GetDensity();
80 fElectronDensity = matCC->GetMaterial()->Get 80 fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
81 matIndex = (G4int)matCC->GetMaterial << 81 matIndex = matCC->GetMaterial()->GetIndex();
82 82
83 fSandia = new G4SandiaTable(matInde 83 fSandia = new G4SandiaTable(matIndex);
84 fIntervalNumber = fSandia->GetMaxInterval() 84 fIntervalNumber = fSandia->GetMaxInterval()-1;
85 85
86 fMatSandiaMatrix = new G4OrderedTable(); 86 fMatSandiaMatrix = new G4OrderedTable();
87 87
88 for (i = 0; i < fIntervalNumber; ++i) << 88 for (i = 0; i < fIntervalNumber; i++)
89 { 89 {
90 fMatSandiaMatrix->push_back(new G4DataVect 90 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
91 } 91 }
92 for (i = 0; i < fIntervalNumber; ++i) << 92 for (i = 0; i < fIntervalNumber; i++)
93 { 93 {
94 (*(*fMatSandiaMatrix)[i])[0] = fSandia->Ge 94 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
95 95
96 for(j = 1; j < 5 ; ++j) << 96 for(j = 1; j < 5 ; j++)
97 { 97 {
98 (*(*fMatSandiaMatrix)[i])[j] = fSandia-> 98 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
99 } 99 }
100 } 100 }
101 KillCloseIntervals(); 101 KillCloseIntervals();
102 Normalisation(); 102 Normalisation();
103 fBetaGammaSq = fTmax = 0.0; 103 fBetaGammaSq = fTmax = 0.0;
104 fIntervalTmax = fCurrentInterval = 0; 104 fIntervalTmax = fCurrentInterval = 0;
105 } 105 }
106 106
107 107
108 108
109 109
110 ////////////////////////////////////////////// 110 ////////////////////////////////////////////////////////////////////////////
111 // 111 //
112 // Destructor 112 // Destructor
113 113
114 G4InitXscPAI::~G4InitXscPAI() 114 G4InitXscPAI::~G4InitXscPAI()
115 { 115 {
116 delete fPAIxscVector; << 116 if(fPAIxscVector) delete fPAIxscVector;
117 delete fPAIdEdxVector; << 117 if(fPAIdEdxVector) delete fPAIdEdxVector;
118 delete fPAIphotonVector; << 118 if(fPAIphotonVector) delete fPAIphotonVector;
119 delete fPAIelectronVector; << 119 if(fPAIelectronVector) delete fPAIelectronVector;
120 delete fChCosSqVector; << 120 if(fChCosSqVector) delete fChCosSqVector;
121 delete fChWidthVector; << 121 if(fChWidthVector) delete fChWidthVector;
122 delete fSandia; 122 delete fSandia;
123 delete fMatSandiaMatrix; 123 delete fMatSandiaMatrix;
124 } 124 }
125 125
126 ////////////////////////////////////////////// 126 ////////////////////////////////////////////////////////////////////////
127 // 127 //
128 // Kill close intervals, recalculate fInterval 128 // Kill close intervals, recalculate fIntervalNumber
129 129
130 void G4InitXscPAI::KillCloseIntervals() 130 void G4InitXscPAI::KillCloseIntervals()
131 { 131 {
132 G4int i, j, k; 132 G4int i, j, k;
133 G4double energy1, energy2; 133 G4double energy1, energy2;
134 134
135 for( i = 0 ; i < fIntervalNumber - 1 ; i++ ) 135 for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
136 { 136 {
137 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 137 energy1 = (*(*fMatSandiaMatrix)[i])[0];
138 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 138 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
139 139
140 if( energy2 - energy1 > 1.5*fDelta*(energy 140 if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ;
141 else 141 else
142 { 142 {
143 for(j = i; j < fIntervalNumber-1; j++) 143 for(j = i; j < fIntervalNumber-1; j++)
144 { 144 {
145 for( k = 0; k < 5; k++ ) 145 for( k = 0; k < 5; k++ )
146 { 146 {
147 (*(*fMatSandiaMatrix)[j])[k] = (*(*f 147 (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
148 } 148 }
149 } 149 }
150 fIntervalNumber-- ; 150 fIntervalNumber-- ;
151 i-- ; 151 i-- ;
152 } 152 }
153 } 153 }
154 154
155 } 155 }
156 156
157 ////////////////////////////////////////////// 157 ////////////////////////////////////////////////////////////////////////
158 // 158 //
159 // Kill close intervals, recalculate fInterval 159 // Kill close intervals, recalculate fIntervalNumber
160 160
161 void G4InitXscPAI::Normalisation() 161 void G4InitXscPAI::Normalisation()
162 { 162 {
163 G4int i, j; 163 G4int i, j;
164 G4double energy1, energy2, /*delta,*/ cof; / 164 G4double energy1, energy2, /*delta,*/ cof; // , shift;
165 165
166 energy1 = (*(*fMatSandiaMatrix)[fIntervalNum 166 energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
167 energy2 = 2.*(*(*fMatSandiaMatrix)[fInterval 167 energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
168 168
169 169
170 cof = RutherfordIntegral(fIntervalNumber-1,e 170 cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
171 171
172 for( i = fIntervalNumber-2; i >= 0; i-- ) 172 for( i = fIntervalNumber-2; i >= 0; i-- )
173 { 173 {
174 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 174 energy1 = (*(*fMatSandiaMatrix)[i])[0];
175 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 175 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
176 176
177 cof += RutherfordIntegral(i,energy1,energy 177 cof += RutherfordIntegral(i,energy1,energy2);
178 // G4cout<<"norm. cof = "<<cof<<G4endl; 178 // G4cout<<"norm. cof = "<<cof<<G4endl;
179 } 179 }
180 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fin 180 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
181 fNormalizationCof *= fElectronDensity; 181 fNormalizationCof *= fElectronDensity;
182 //delta = fNormalizationCof - cof; 182 //delta = fNormalizationCof - cof;
183 fNormalizationCof /= cof; 183 fNormalizationCof /= cof;
184 // G4cout<<"G4InitXscPAI::fNormalizationCof 184 // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
185 // <<"; at delta ="<<delta<<G4endl ; 185 // <<"; at delta ="<<delta<<G4endl ;
186 186
187 for (i = 0; i < fIntervalNumber; i++) // ren << 187 for (G4int i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
188 { 188 {
189 for(j = 1; j < 5 ; j++) 189 for(j = 1; j < 5 ; j++)
190 { 190 {
191 (*(*fMatSandiaMatrix)[i])[j] *= fNormali 191 (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
192 } 192 }
193 } 193 }
194 /* 194 /*
195 if(delta > 0) // shift the first energy inte 195 if(delta > 0) // shift the first energy interval
196 { 196 {
197 for(i=1;i<100;i++) 197 for(i=1;i<100;i++)
198 { 198 {
199 energy1 = (1.-i/100.)*(*(*fMatSandiaMatr 199 energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
200 energy2 = (*(*fMatSandiaMatrix)[0])[0]; 200 energy2 = (*(*fMatSandiaMatrix)[0])[0];
201 shift = RutherfordIntegral(0,energy1,e 201 shift = RutherfordIntegral(0,energy1,energy2);
202 G4cout<<shift<<"\t"; 202 G4cout<<shift<<"\t";
203 if(shift >= delta) break; 203 if(shift >= delta) break;
204 } 204 }
205 (*(*fMatSandiaMatrix)[0])[0] = energy1; 205 (*(*fMatSandiaMatrix)[0])[0] = energy1;
206 cof += shift; 206 cof += shift;
207 } 207 }
208 else if(delta < 0) 208 else if(delta < 0)
209 { 209 {
210 for(i=1;i<100;i++) 210 for(i=1;i<100;i++)
211 { 211 {
212 energy1 = (*(*fMatSandiaMatrix)[0])[0]; 212 energy1 = (*(*fMatSandiaMatrix)[0])[0];
213 energy2 = (*(*fMatSandiaMatrix)[0])[0] + 213 energy2 = (*(*fMatSandiaMatrix)[0])[0] +
214 ( (*(*fMatSandiaMatrix)[0])[0] - 214 ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
215 shift = RutherfordIntegral(0,energy1,e 215 shift = RutherfordIntegral(0,energy1,energy2);
216 if( shift >= std::abs(delta) ) break; 216 if( shift >= std::abs(delta) ) break;
217 } 217 }
218 (*(*fMatSandiaMatrix)[0])[0] = energy2; 218 (*(*fMatSandiaMatrix)[0])[0] = energy2;
219 cof -= shift; 219 cof -= shift;
220 } 220 }
221 G4cout<<G4cout<<"G4InitXscPAI::fNormalizatio 221 G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
222 <<"; at delta ="<<delta<<" and i = " 222 <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ;
223 */ 223 */
224 } 224 }
225 225
>> 226
>> 227
>> 228
>> 229
226 ////////////////////////////////////////////// 230 ////////////////////////////////////////////////////////////////////
227 // 231 //
228 // Integration over electrons that could be co 232 // Integration over electrons that could be considered
229 // quasi-free at energy transfer of interest 233 // quasi-free at energy transfer of interest
230 234
231 G4double G4InitXscPAI::RutherfordIntegral( G4i 235 G4double G4InitXscPAI::RutherfordIntegral( G4int k,
232 G4double x1, 236 G4double x1,
233 G4double x2 ) 237 G4double x2 )
234 { 238 {
235 G4double c1, c2, c3, a1, a2, a3, a4 ; 239 G4double c1, c2, c3, a1, a2, a3, a4 ;
236 240
237 a1 = (*(*fMatSandiaMatrix)[k])[1]; 241 a1 = (*(*fMatSandiaMatrix)[k])[1];
238 a2 = (*(*fMatSandiaMatrix)[k])[2]; 242 a2 = (*(*fMatSandiaMatrix)[k])[2];
239 a3 = (*(*fMatSandiaMatrix)[k])[3]; 243 a3 = (*(*fMatSandiaMatrix)[k])[3];
240 a4 = (*(*fMatSandiaMatrix)[k])[4]; 244 a4 = (*(*fMatSandiaMatrix)[k])[4];
241 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<< 245 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
242 c1 = (x2 - x1)/x1/x2 ; 246 c1 = (x2 - x1)/x1/x2 ;
243 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ; 247 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
244 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x 248 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
245 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "< 249 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
246 250
247 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a 251 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
248 252
249 } // end of RutherfordIntegral 253 } // end of RutherfordIntegral
250 254
251 ////////////////////////////////////////////// 255 ///////////////////////////////////////////////////////////////
252 // 256 //
253 // Integrate photo-absorption cross-section f 257 // Integrate photo-absorption cross-section from I1 up to omega
254 258
255 G4double G4InitXscPAI::IntegralTerm(G4double o 259 G4double G4InitXscPAI::IntegralTerm(G4double omega)
256 { 260 {
257 G4int i; 261 G4int i;
258 G4double energy1, energy2, result = 0.; 262 G4double energy1, energy2, result = 0.;
259 263
260 for( i = 0; i <= fIntervalTmax; i++ ) 264 for( i = 0; i <= fIntervalTmax; i++ )
261 { 265 {
262 if(i == fIntervalTmax) 266 if(i == fIntervalTmax)
263 { 267 {
264 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 268 energy1 = (*(*fMatSandiaMatrix)[i])[0];
265 result += RutherfordIntegral(i,energy1,o 269 result += RutherfordIntegral(i,energy1,omega);
266 } 270 }
267 else 271 else
268 { 272 {
269 if( omega <= (*(*fMatSandiaMatrix)[i+1]) 273 if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
270 { 274 {
271 energy1 = (*(*fMatSandiaMatrix)[i])[0] 275 energy1 = (*(*fMatSandiaMatrix)[i])[0];
272 result += RutherfordIntegral(i,energy1 276 result += RutherfordIntegral(i,energy1,omega);
273 break; 277 break;
274 } 278 }
275 else 279 else
276 { 280 {
277 energy1 = (*(*fMatSandiaMatrix)[i])[0] 281 energy1 = (*(*fMatSandiaMatrix)[i])[0];
278 energy2 = (*(*fMatSandiaMatrix)[i+1])[ 282 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
279 result += RutherfordIntegral(i,energy1 283 result += RutherfordIntegral(i,energy1,energy2);
280 } 284 }
281 } 285 }
282 // G4cout<<"IntegralTerm<<"("<<omega<<")"< 286 // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
283 } 287 }
284 return result; 288 return result;
285 } 289 }
286 290
287 291
288 ////////////////////////////////////////////// 292 ////////////////////////////////////////////////////////////////
289 // 293 //
290 // Imaginary part of dielectric constant 294 // Imaginary part of dielectric constant
291 // (G4int k - interval number, G4double en1 - 295 // (G4int k - interval number, G4double en1 - energy point)
292 296
293 G4double G4InitXscPAI::ImPartDielectricConst( 297 G4double G4InitXscPAI::ImPartDielectricConst( G4int k ,
294 G4double energy1 298 G4double energy1 )
295 { 299 {
296 G4double energy2,energy3,energy4,a1,a2,a3,a 300 G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
297 301
298 a1 = (*(*fMatSandiaMatrix)[k])[1]; 302 a1 = (*(*fMatSandiaMatrix)[k])[1];
299 a2 = (*(*fMatSandiaMatrix)[k])[2]; 303 a2 = (*(*fMatSandiaMatrix)[k])[2];
300 a3 = (*(*fMatSandiaMatrix)[k])[3]; 304 a3 = (*(*fMatSandiaMatrix)[k])[3];
301 a4 = (*(*fMatSandiaMatrix)[k])[4]; 305 a4 = (*(*fMatSandiaMatrix)[k])[4];
302 306
303 energy2 = energy1*energy1; 307 energy2 = energy1*energy1;
304 energy3 = energy2*energy1; 308 energy3 = energy2*energy1;
305 energy4 = energy3*energy1; 309 energy4 = energy3*energy1;
306 310
307 result = a1/energy1+a2/energy2+a3/energy3+ 311 result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
308 result *= hbarc/energy1 ; 312 result *= hbarc/energy1 ;
309 313
310 return result ; 314 return result ;
311 315
312 } // end of ImPartDielectricConst 316 } // end of ImPartDielectricConst
313 317
314 ////////////////////////////////////////////// 318 ////////////////////////////////////////////////////////////////
315 // 319 //
316 // Modulus squared of dielectric constant 320 // Modulus squared of dielectric constant
317 // (G4int k - interval number, G4double omega 321 // (G4int k - interval number, G4double omega - energy point)
318 322
319 G4double G4InitXscPAI::ModuleSqDielectricConst 323 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int k ,
320 G4double omega ) 324 G4double omega )
321 { 325 {
322 G4double eIm2, eRe2, result; 326 G4double eIm2, eRe2, result;
323 327
324 result = ImPartDielectricConst(k,omega); 328 result = ImPartDielectricConst(k,omega);
325 eIm2 = result*result; 329 eIm2 = result*result;
326 330
327 result = RePartDielectricConst(omega); 331 result = RePartDielectricConst(omega);
328 eRe2 = result*result; 332 eRe2 = result*result;
329 333
330 result = eIm2 + eRe2; 334 result = eIm2 + eRe2;
331 335
332 return result ; 336 return result ;
333 } 337 }
334 338
335 339
336 ////////////////////////////////////////////// 340 //////////////////////////////////////////////////////////////////////////////
337 // 341 //
338 // Real part of dielectric constant minus unit 342 // Real part of dielectric constant minus unit: epsilon_1 - 1
339 // (G4double enb - energy point) 343 // (G4double enb - energy point)
340 // 344 //
341 345
342 G4double G4InitXscPAI::RePartDielectricConst(G 346 G4double G4InitXscPAI::RePartDielectricConst(G4double enb)
343 { 347 {
344 G4int i; 348 G4int i;
345 G4double x0, x02, x03, x04, x05, x1, x2, a1 349 G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
346 c1, c2, c3, cof1, cof2, xln1, xln2 350 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
347 351
348 x0 = enb ; 352 x0 = enb ;
349 result = 0 ; 353 result = 0 ;
350 354
351 for( i = 0; i < fIntervalNumber-1; i++) 355 for( i = 0; i < fIntervalNumber-1; i++)
352 { 356 {
353 x1 = (*(*fMatSandiaMatrix)[i])[0]; 357 x1 = (*(*fMatSandiaMatrix)[i])[0];
354 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ; 358 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
355 359
356 a1 = (*(*fMatSandiaMatrix)[i])[1]; 360 a1 = (*(*fMatSandiaMatrix)[i])[1];
357 a2 = (*(*fMatSandiaMatrix)[i])[2]; 361 a2 = (*(*fMatSandiaMatrix)[i])[2];
358 a3 = (*(*fMatSandiaMatrix)[i])[3]; 362 a3 = (*(*fMatSandiaMatrix)[i])[3];
359 a4 = (*(*fMatSandiaMatrix)[i])[4]; 363 a4 = (*(*fMatSandiaMatrix)[i])[4];
360 364
361 if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta 365 if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta )
362 { 366 {
363 if(x0 >= x1) x0 = x1*(1+fDelta); 367 if(x0 >= x1) x0 = x1*(1+fDelta);
364 else x0 = x1*(1-fDelta); 368 else x0 = x1*(1-fDelta);
365 } 369 }
366 if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta 370 if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta )
367 { 371 {
368 if(x0 >= x2) x0 = x2*(1+fDelta); 372 if(x0 >= x2) x0 = x2*(1+fDelta);
369 else x0 = x2*(1-fDelta); 373 else x0 = x2*(1-fDelta);
370 } 374 }
371 xx1 = x1 - x0 ; 375 xx1 = x1 - x0 ;
372 xx2 = x2 - x0 ; 376 xx2 = x2 - x0 ;
373 xx12 = xx2/xx1 ; 377 xx12 = xx2/xx1 ;
374 378
375 if( xx12 < 0 ) xx12 = -xx12; 379 if( xx12 < 0 ) xx12 = -xx12;
376 380
377 xln1 = log(x2/x1) ; 381 xln1 = log(x2/x1) ;
378 xln2 = log(xx12) ; 382 xln2 = log(xx12) ;
379 xln3 = log((x2 + x0)/(x1 + x0)) ; 383 xln3 = log((x2 + x0)/(x1 + x0)) ;
380 384
381 x02 = x0*x0 ; 385 x02 = x0*x0 ;
382 x03 = x02*x0 ; 386 x03 = x02*x0 ;
383 x04 = x03*x0 ; 387 x04 = x03*x0 ;
384 x05 = x04*x0; 388 x05 = x04*x0;
385 389
386 c1 = (x2 - x1)/x1/x2 ; 390 c1 = (x2 - x1)/x1/x2 ;
387 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ; 391 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
388 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x 392 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
389 393
390 result -= (a1/x02 + a3/x04)*xln1 ; 394 result -= (a1/x02 + a3/x04)*xln1 ;
391 result -= (a2/x02 + a4/x04)*c1 ; 395 result -= (a2/x02 + a4/x04)*c1 ;
392 result -= a3*c2/2/x02 ; 396 result -= a3*c2/2/x02 ;
393 result -= a4*c3/3/x02 ; 397 result -= a4*c3/3/x02 ;
394 398
395 cof1 = a1/x02 + a3/x04 ; 399 cof1 = a1/x02 + a3/x04 ;
396 cof2 = a2/x03 + a4/x05 ; 400 cof2 = a2/x03 + a4/x05 ;
397 401
398 result += 0.5*(cof1 +cof2)*xln2 ; 402 result += 0.5*(cof1 +cof2)*xln2 ;
399 result += 0.5*(cof1 - cof2)*xln3 ; 403 result += 0.5*(cof1 - cof2)*xln3 ;
400 } 404 }
401 result *= 2*hbarc/pi ; 405 result *= 2*hbarc/pi ;
402 406
403 return result ; 407 return result ;
404 408
405 } // end of RePartDielectricConst 409 } // end of RePartDielectricConst
406 410
407 ////////////////////////////////////////////// 411 //////////////////////////////////////////////////////////////////////
408 // 412 //
409 // PAI differential cross-section in terms of 413 // PAI differential cross-section in terms of
410 // simplified Allison's equation 414 // simplified Allison's equation
411 // 415 //
412 416
413 G4double G4InitXscPAI::DifPAIxSection( G4doubl 417 G4double G4InitXscPAI::DifPAIxSection( G4double omega )
414 { 418 {
415 G4int i = fCurrentInterval; 419 G4int i = fCurrentInterval;
416 G4double betaGammaSq = fBetaGammaSq; 420 G4double betaGammaSq = fBetaGammaSq;
417 G4double integralTerm = IntegralTerm(omega); 421 G4double integralTerm = IntegralTerm(omega);
418 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,res 422 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
419 G4double epsilonRe = RePartDielectricConst(o 423 G4double epsilonRe = RePartDielectricConst(omega);
420 G4double epsilonIm = ImPartDielectricConst(i 424 G4double epsilonIm = ImPartDielectricConst(i,omega);
421 G4double be4 ; 425 G4double be4 ;
422 static const G4double betaBohr2 = fine_struc << 426 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
423 static const G4double betaBohr4 = betaBohr2* << 427 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
424 be2 = betaGammaSq/(1 + betaGammaSq) ; 428 be2 = betaGammaSq/(1 + betaGammaSq) ;
425 be4 = be2*be2 ; 429 be4 = be2*be2 ;
426 430
427 cof = 1 ; 431 cof = 1 ;
428 x1 = log(2*electron_mass_c2/omega) ; 432 x1 = log(2*electron_mass_c2/omega) ;
429 433
430 if( betaGammaSq < 0.01 ) x2 = log(be2) ; 434 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
431 else 435 else
432 { 436 {
433 x2 = -log( (1/betaGammaSq - epsilonRe)* 437 x2 = -log( (1/betaGammaSq - epsilonRe)*
434 (1/betaGammaSq - epsilonRe) + 438 (1/betaGammaSq - epsilonRe) +
435 epsilonIm*epsilonIm )/2 ; 439 epsilonIm*epsilonIm )/2 ;
436 } 440 }
437 if( epsilonIm == 0.0 || betaGammaSq < 0.01 441 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
438 { 442 {
439 x6=0 ; 443 x6=0 ;
440 } 444 }
441 else 445 else
442 { 446 {
443 x3 = -epsilonRe + 1/betaGammaSq ; 447 x3 = -epsilonRe + 1/betaGammaSq ;
444 x5 = -1 - epsilonRe + be2*((1 +epsilonRe) 448 x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
445 epsilonIm*epsilonIm) ; 449 epsilonIm*epsilonIm) ;
446 450
447 x7 = atan2(epsilonIm,x3) ; 451 x7 = atan2(epsilonIm,x3) ;
448 x6 = x5 * x7 ; 452 x6 = x5 * x7 ;
449 } 453 }
450 // if(fImPartDielectricConst[i] == 0) x6 = 454 // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
451 455
452 x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ; 456 x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
453 // if( x4 < 0.0 ) x4 = 0.0 ; 457 // if( x4 < 0.0 ) x4 = 0.0 ;
454 x8 = (1 + epsilonRe)*(1 + epsilonRe) + 458 x8 = (1 + epsilonRe)*(1 + epsilonRe) +
455 epsilonIm*epsilonIm; 459 epsilonIm*epsilonIm;
456 460
457 result = (x4 + cof*integralTerm/omega/omega 461 result = (x4 + cof*integralTerm/omega/omega) ;
458 if(result < 1.0e-8) result = 1.0e-8 ; 462 if(result < 1.0e-8) result = 1.0e-8 ;
459 result *= fine_structure_const/be2/pi ; 463 result *= fine_structure_const/be2/pi ;
460 // result *= (1-exp(-beta/betaBohr))*(1-e 464 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
461 // result *= (1-exp(-be2/betaBohr2)) ; 465 // result *= (1-exp(-be2/betaBohr2)) ;
462 result *= (1-exp(-be4/betaBohr4)) ; 466 result *= (1-exp(-be4/betaBohr4)) ;
463 if(fDensity >= fSolidDensity) 467 if(fDensity >= fSolidDensity)
464 { 468 {
465 result /= x8 ; 469 result /= x8 ;
466 } 470 }
467 return result ; 471 return result ;
468 472
469 } // end of DifPAIxSection 473 } // end of DifPAIxSection
470 474
471 ////////////////////////////////////////////// 475 //////////////////////////////////////////////////////////////////////
472 // 476 //
473 // Differential PAI dEdx(omega)=omega*dNdx(ome 477 // Differential PAI dEdx(omega)=omega*dNdx(omega)
474 // 478 //
475 479
476 G4double G4InitXscPAI::DifPAIdEdx( G4double om 480 G4double G4InitXscPAI::DifPAIdEdx( G4double omega )
477 { 481 {
478 G4double dEdx = omega*DifPAIxSection(omega); 482 G4double dEdx = omega*DifPAIxSection(omega);
479 return dEdx; 483 return dEdx;
480 } 484 }
481 485
482 ////////////////////////////////////////////// 486 //////////////////////////////////////////////////////////////////////////
483 // 487 //
484 // Calculation od dN/dx of collisions with cre 488 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
485 489
486 G4double G4InitXscPAI::PAIdNdxCherenkov( G4dou 490 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega )
487 { 491 {
488 G4int i = fCurrentInterval; 492 G4int i = fCurrentInterval;
489 G4double betaGammaSq = fBetaGammaSq; 493 G4double betaGammaSq = fBetaGammaSq;
490 G4double epsilonRe = RePartDielectricConst(o 494 G4double epsilonRe = RePartDielectricConst(omega);
491 G4double epsilonIm = ImPartDielectricConst(i 495 G4double epsilonIm = ImPartDielectricConst(i,omega);
492 496
493 G4double /*cof,*/ logarithm, x3, x5, argumen 497 G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ;
494 G4double be2, be4; << 498 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
495 499
496 //cof = 1.0 ; 500 //cof = 1.0 ;
497 static const G4double cofBetaBohr = 4.0 ; << 501 cofBetaBohr = 4.0 ;
498 static const G4double betaBohr2 = fine_str << 502 betaBohr2 = fine_structure_const*fine_structure_const ;
499 static const G4double betaBohr4 = betaBohr << 503 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
500 504
501 be2 = betaGammaSq/(1 + betaGammaSq) ; 505 be2 = betaGammaSq/(1 + betaGammaSq) ;
502 be4 = be2*be2 ; 506 be4 = be2*be2 ;
503 507
504 if( betaGammaSq < 0.01 ) logarithm = log(1. 508 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
505 else 509 else
506 { 510 {
507 logarithm = -log( (1/betaGammaSq - epsil 511 logarithm = -log( (1/betaGammaSq - epsilonRe)*
508 (1/betaGammaSq - epsilonRe) 512 (1/betaGammaSq - epsilonRe) +
509 epsilonIm*epsilonIm )*0.5 ; 513 epsilonIm*epsilonIm )*0.5 ;
510 logarithm += log(1+1.0/betaGammaSq) ; 514 logarithm += log(1+1.0/betaGammaSq) ;
511 } 515 }
512 516
513 if( epsilonIm == 0.0 || betaGammaSq < 0.01 517 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
514 { 518 {
515 argument = 0.0 ; 519 argument = 0.0 ;
516 } 520 }
517 else 521 else
518 { 522 {
519 x3 = -epsilonRe + 1.0/betaGammaSq ; 523 x3 = -epsilonRe + 1.0/betaGammaSq ;
520 x5 = -1.0 - epsilonRe + 524 x5 = -1.0 - epsilonRe +
521 be2*((1.0 +epsilonRe)*(1.0 + epsilon 525 be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
522 epsilonIm*epsilonIm) ; 526 epsilonIm*epsilonIm) ;
523 if( x3 == 0.0 ) argument = 0.5*pi; 527 if( x3 == 0.0 ) argument = 0.5*pi;
524 else argument = atan2(epsilonI 528 else argument = atan2(epsilonIm,x3) ;
525 argument *= x5 ; 529 argument *= x5 ;
526 } 530 }
527 dNdxC = ( logarithm*epsilonIm + argument )/ 531 dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
528 532
529 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ; 533 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
530 534
531 dNdxC *= fine_structure_const/be2/pi ; 535 dNdxC *= fine_structure_const/be2/pi ;
532 536
533 dNdxC *= (1-exp(-be4/betaBohr4)) ; 537 dNdxC *= (1-exp(-be4/betaBohr4)) ;
534 538
535 if(fDensity >= fSolidDensity) 539 if(fDensity >= fSolidDensity)
536 { 540 {
537 modul2 = (1.0 + epsilonRe)*(1.0 + epsilo 541 modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
538 epsilonIm*epsilonIm; 542 epsilonIm*epsilonIm;
539 dNdxC /= modul2 ; 543 dNdxC /= modul2 ;
540 } 544 }
541 return dNdxC ; 545 return dNdxC ;
542 546
543 } // end of PAIdNdxCerenkov 547 } // end of PAIdNdxCerenkov
544 548
545 ////////////////////////////////////////////// 549 //////////////////////////////////////////////////////////////////////////
546 // 550 //
547 // Calculation od dN/dx of collisions with cre 551 // Calculation od dN/dx of collisions with creation of longitudinal EM
548 // excitations (plasmons, delta-electrons) 552 // excitations (plasmons, delta-electrons)
549 553
550 G4double G4InitXscPAI::PAIdNdxPlasmon( G4doubl 554 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega )
551 { 555 {
552 G4int i = fCurrentInterval; 556 G4int i = fCurrentInterval;
553 G4double betaGammaSq = fBetaGammaSq; 557 G4double betaGammaSq = fBetaGammaSq;
554 G4double integralTerm = IntegralTerm(omega); 558 G4double integralTerm = IntegralTerm(omega);
555 G4double epsilonRe = RePartDielectricConst(o 559 G4double epsilonRe = RePartDielectricConst(omega);
556 G4double epsilonIm = ImPartDielectricConst(i 560 G4double epsilonIm = ImPartDielectricConst(i,omega);
557 561
558 G4double cof, resonance, modul2, dNdxP ; 562 G4double cof, resonance, modul2, dNdxP ;
559 G4double be2, be4; << 563 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
560 564
561 cof = 1 ; 565 cof = 1 ;
562 static const G4double cofBetaBohr = 4.0 ; << 566 cofBetaBohr = 4.0 ;
563 static const G4double betaBohr2 = fine_st << 567 betaBohr2 = fine_structure_const*fine_structure_const ;
564 static const G4double betaBohr4 = betaBoh << 568 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
565 569
566 be2 = betaGammaSq/(1 + betaGammaSq) ; 570 be2 = betaGammaSq/(1 + betaGammaSq) ;
567 be4 = be2*be2 ; 571 be4 = be2*be2 ;
568 572
569 resonance = log(2*electron_mass_c2*be2/ome 573 resonance = log(2*electron_mass_c2*be2/omega) ;
570 resonance *= epsilonIm/hbarc ; 574 resonance *= epsilonIm/hbarc ;
571 575
572 576
573 dNdxP = ( resonance + cof*integralTerm/omeg 577 dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
574 578
575 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ; 579 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
576 580
577 dNdxP *= fine_structure_const/be2/pi ; 581 dNdxP *= fine_structure_const/be2/pi ;
578 dNdxP *= (1-exp(-be4/betaBohr4)) ; 582 dNdxP *= (1-exp(-be4/betaBohr4)) ;
579 583
580 if( fDensity >= fSolidDensity ) 584 if( fDensity >= fSolidDensity )
581 { 585 {
582 modul2 = (1 + epsilonRe)*(1 + epsilonRe) 586 modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
583 epsilonIm*epsilonIm; 587 epsilonIm*epsilonIm;
584 dNdxP /= modul2 ; 588 dNdxP /= modul2 ;
585 } 589 }
586 return dNdxP ; 590 return dNdxP ;
587 591
588 } // end of PAIdNdxPlasmon 592 } // end of PAIdNdxPlasmon
589 593
590 ////////////////////////////////////////////// 594 ////////////////////////////////////////////////////////////////////////
591 // 595 //
592 // Calculation of the PAI integral cross-secti 596 // Calculation of the PAI integral cross-section
593 // = specific primary ionisation, 1/cm 597 // = specific primary ionisation, 1/cm
594 // 598 //
595 599
596 void G4InitXscPAI::IntegralPAIxSection(G4doubl 600 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax)
597 { 601 {
598 G4int i,k,i1,i2; 602 G4int i,k,i1,i2;
599 G4double energy1, energy2, result = 0.; 603 G4double energy1, energy2, result = 0.;
600 604
601 fBetaGammaSq = bg2; 605 fBetaGammaSq = bg2;
602 fTmax = Tmax; 606 fTmax = Tmax;
603 607
604 delete fPAIxscVector; << 608 if(fPAIxscVector) delete fPAIxscVector;
605 609
606 fPAIxscVector = new G4PhysicsLogVector( (*(* 610 fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
607 fPAIxscVector->PutValue(fPAIbin-1,result); 611 fPAIxscVector->PutValue(fPAIbin-1,result);
608 612
609 for( i = fIntervalNumber - 1; i >= 0; i-- ) 613 for( i = fIntervalNumber - 1; i >= 0; i-- )
610 { 614 {
611 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 615 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
612 } 616 }
613 if (i < 0) i = 0; // Tmax should be more tha 617 if (i < 0) i = 0; // Tmax should be more than
614 // first ionisation potent 618 // first ionisation potential
615 fIntervalTmax = i; 619 fIntervalTmax = i;
616 620
617 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 621 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
618 622
619 for( k = fPAIbin - 2; k >= 0; k-- ) 623 for( k = fPAIbin - 2; k >= 0; k-- )
620 { 624 {
621 energy1 = fPAIxscVector->GetLowEdgeEnergy( 625 energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
622 energy2 = fPAIxscVector->GetLowEdgeEnergy( 626 energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
623 627
624 for( i = fIntervalTmax; i >= 0; i-- ) 628 for( i = fIntervalTmax; i >= 0; i-- )
625 { 629 {
626 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 630 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
627 } 631 }
628 if(i < 0) i = 0; 632 if(i < 0) i = 0;
629 i2 = i; 633 i2 = i;
630 634
631 for( i = fIntervalTmax; i >= 0; i-- ) 635 for( i = fIntervalTmax; i >= 0; i-- )
632 { 636 {
633 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 637 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
634 } 638 }
635 if(i < 0) i = 0; 639 if(i < 0) i = 0;
636 i1 = i; 640 i1 = i;
637 641
638 if( i1 == i2 ) 642 if( i1 == i2 )
639 { 643 {
640 fCurrentInterval = i1; 644 fCurrentInterval = i1;
641 result += integral.Legendre10(this,&G4In 645 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
642 energy1,en 646 energy1,energy2);
643 fPAIxscVector->PutValue(k,result); 647 fPAIxscVector->PutValue(k,result);
644 } 648 }
645 else 649 else
646 { 650 {
647 for( i = i2; i >= i1; i-- ) 651 for( i = i2; i >= i1; i-- )
648 { 652 {
649 fCurrentInterval = i; 653 fCurrentInterval = i;
650 654
651 if( i==i2 ) result += integral. 655 if( i==i2 ) result += integral.Legendre10(this,
652 &G4InitXscPAI::DifP 656 &G4InitXscPAI::DifPAIxSection,
653 (*(*fMatSandiaMatri 657 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
654 658
655 else if( i == i1 ) result += integral.Legend 659 else if( i == i1 ) result += integral.Legendre10(this,
656 &G4InitXscPAI::DifP 660 &G4InitXscPAI::DifPAIxSection,energy1,
657 (*(*fMatSandiaMatri 661 (*(*fMatSandiaMatrix)[i+1])[0]);
658 662
659 else result += integral. 663 else result += integral.Legendre10(this,
660 &G4InitXscPAI::DifP 664 &G4InitXscPAI::DifPAIxSection,
661 (*(*fMatSandiaMatrix)[i 665 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
662 } 666 }
663 fPAIxscVector->PutValue(k,result); 667 fPAIxscVector->PutValue(k,result);
664 } 668 }
665 // G4cout<<k<<"\t"<<result<<G4endl; 669 // G4cout<<k<<"\t"<<result<<G4endl;
666 } 670 }
667 return ; 671 return ;
668 } 672 }
669 673
670 674
671 ////////////////////////////////////////////// 675 ////////////////////////////////////////////////////////////////////////
672 // 676 //
673 // Calculation of the PAI integral dEdx 677 // Calculation of the PAI integral dEdx
674 // = mean energy loss per unit length, keV/cm 678 // = mean energy loss per unit length, keV/cm
675 // 679 //
676 680
677 void G4InitXscPAI::IntegralPAIdEdx(G4double bg 681 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax)
678 { 682 {
679 G4int i,k,i1,i2; 683 G4int i,k,i1,i2;
680 G4double energy1, energy2, result = 0.; 684 G4double energy1, energy2, result = 0.;
681 685
682 fBetaGammaSq = bg2; 686 fBetaGammaSq = bg2;
683 fTmax = Tmax; 687 fTmax = Tmax;
684 688
685 delete fPAIdEdxVector; << 689 if(fPAIdEdxVector) delete fPAIdEdxVector;
686 690
687 fPAIdEdxVector = new G4PhysicsLogVector( (*( 691 fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
688 fPAIdEdxVector->PutValue(fPAIbin-1,result); 692 fPAIdEdxVector->PutValue(fPAIbin-1,result);
689 693
690 for( i = fIntervalNumber - 1; i >= 0; i-- ) 694 for( i = fIntervalNumber - 1; i >= 0; i-- )
691 { 695 {
692 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 696 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
693 } 697 }
694 if (i < 0) i = 0; // Tmax should be more tha 698 if (i < 0) i = 0; // Tmax should be more than
695 // first ionisation potent 699 // first ionisation potential
696 fIntervalTmax = i; 700 fIntervalTmax = i;
697 701
698 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 702 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
699 703
700 for( k = fPAIbin - 2; k >= 0; k-- ) 704 for( k = fPAIbin - 2; k >= 0; k-- )
701 { 705 {
702 energy1 = fPAIdEdxVector->GetLowEdgeEnergy 706 energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
703 energy2 = fPAIdEdxVector->GetLowEdgeEnergy 707 energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
704 708
705 for( i = fIntervalTmax; i >= 0; i-- ) 709 for( i = fIntervalTmax; i >= 0; i-- )
706 { 710 {
707 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 711 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
708 } 712 }
709 if(i < 0) i = 0; 713 if(i < 0) i = 0;
710 i2 = i; 714 i2 = i;
711 715
712 for( i = fIntervalTmax; i >= 0; i-- ) 716 for( i = fIntervalTmax; i >= 0; i-- )
713 { 717 {
714 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 718 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
715 } 719 }
716 if(i < 0) i = 0; 720 if(i < 0) i = 0;
717 i1 = i; 721 i1 = i;
718 722
719 if( i1 == i2 ) 723 if( i1 == i2 )
720 { 724 {
721 fCurrentInterval = i1; 725 fCurrentInterval = i1;
722 result += integral.Legendre10(this,&G4In 726 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
723 energy1,en 727 energy1,energy2);
724 fPAIdEdxVector->PutValue(k,result); 728 fPAIdEdxVector->PutValue(k,result);
725 } 729 }
726 else 730 else
727 { 731 {
728 for( i = i2; i >= i1; i-- ) 732 for( i = i2; i >= i1; i-- )
729 { 733 {
730 fCurrentInterval = i; 734 fCurrentInterval = i;
731 735
732 if( i==i2 ) result += integral. 736 if( i==i2 ) result += integral.Legendre10(this,
733 &G4InitXscPAI::DifP 737 &G4InitXscPAI::DifPAIdEdx,
734 (*(*fMatSandiaMatri 738 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
735 739
736 else if( i == i1 ) result += integral.Legend 740 else if( i == i1 ) result += integral.Legendre10(this,
737 &G4InitXscPAI::DifP 741 &G4InitXscPAI::DifPAIdEdx,energy1,
738 (*(*fMatSandiaMatri 742 (*(*fMatSandiaMatrix)[i+1])[0]);
739 743
740 else result += integral. 744 else result += integral.Legendre10(this,
741 &G4InitXscPAI::DifP 745 &G4InitXscPAI::DifPAIdEdx,
742 (*(*fMatSandiaMatrix)[i 746 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
743 } 747 }
744 fPAIdEdxVector->PutValue(k,result); 748 fPAIdEdxVector->PutValue(k,result);
745 } 749 }
746 // G4cout<<k<<"\t"<<result<<G4endl; 750 // G4cout<<k<<"\t"<<result<<G4endl;
747 } 751 }
748 return ; 752 return ;
749 } 753 }
750 754
751 ////////////////////////////////////////////// 755 ////////////////////////////////////////////////////////////////////////
752 // 756 //
753 // Calculation of the PAI Cerenkov integral cr 757 // Calculation of the PAI Cerenkov integral cross-section
754 // fIntegralCrenkov[1] = specific Crenkov ioni 758 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
755 // and fIntegralCerenkov[0] = mean Cerenkov lo 759 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
756 760
757 void G4InitXscPAI::IntegralCherenkov(G4double 761 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax)
758 { 762 {
759 G4int i,k,i1,i2; 763 G4int i,k,i1,i2;
760 G4double energy1, energy2, beta2, module2, c 764 G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
761 765
762 fBetaGammaSq = bg2; 766 fBetaGammaSq = bg2;
763 fTmax = Tmax; 767 fTmax = Tmax;
764 beta2 = bg2/(1+bg2); 768 beta2 = bg2/(1+bg2);
765 769
766 delete fPAIphotonVector; << 770 if(fPAIphotonVector) delete fPAIphotonVector;
767 delete fChCosSqVector; << 771 if(fChCosSqVector) delete fChCosSqVector;
768 delete fChWidthVector; << 772 if(fChWidthVector) delete fChWidthVector;
769 773
770 fPAIphotonVector = new G4PhysicsLogVector( ( 774 fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
771 fChCosSqVector = new G4PhysicsLogVector( (*( 775 fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
772 fChWidthVector = new G4PhysicsLogVector( (*( 776 fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
773 777
774 fPAIphotonVector->PutValue(fPAIbin-1,result) 778 fPAIphotonVector->PutValue(fPAIbin-1,result);
775 fChCosSqVector->PutValue(fPAIbin-1,1.); 779 fChCosSqVector->PutValue(fPAIbin-1,1.);
776 fChWidthVector->PutValue(fPAIbin-1,1e-7); 780 fChWidthVector->PutValue(fPAIbin-1,1e-7);
777 781
778 for( i = fIntervalNumber - 1; i >= 0; i-- ) 782 for( i = fIntervalNumber - 1; i >= 0; i-- )
779 { 783 {
780 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 784 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
781 } 785 }
782 if (i < 0) i = 0; // Tmax should be more tha 786 if (i < 0) i = 0; // Tmax should be more than
783 // first ionisation potent 787 // first ionisation potential
784 fIntervalTmax = i; 788 fIntervalTmax = i;
785 789
786 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 790 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
787 791
788 for( k = fPAIbin - 2; k >= 0; k-- ) 792 for( k = fPAIbin - 2; k >= 0; k-- )
789 { 793 {
790 energy1 = fPAIphotonVector->GetLowEdgeEner 794 energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
791 energy2 = fPAIphotonVector->GetLowEdgeEner 795 energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
792 796
793 for( i = fIntervalTmax; i >= 0; i-- ) 797 for( i = fIntervalTmax; i >= 0; i-- )
794 { 798 {
795 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 799 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
796 } 800 }
797 if(i < 0) i = 0; 801 if(i < 0) i = 0;
798 i2 = i; 802 i2 = i;
799 803
800 for( i = fIntervalTmax; i >= 0; i-- ) 804 for( i = fIntervalTmax; i >= 0; i-- )
801 { 805 {
802 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 806 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
803 } 807 }
804 if(i < 0) i = 0; 808 if(i < 0) i = 0;
805 i1 = i; 809 i1 = i;
806 810
807 module2 = ModuleSqDielectricConst(i1,energ 811 module2 = ModuleSqDielectricConst(i1,energy1);
808 cos2 = RePartDielectricConst(energy1)/m 812 cos2 = RePartDielectricConst(energy1)/module2/beta2;
809 width = ImPartDielectricConst(i1,energy1 813 width = ImPartDielectricConst(i1,energy1)/module2/beta2;
810 814
811 fChCosSqVector->PutValue(k,cos2); 815 fChCosSqVector->PutValue(k,cos2);
812 fChWidthVector->PutValue(k,width); 816 fChWidthVector->PutValue(k,width);
813 817
814 if( i1 == i2 ) 818 if( i1 == i2 )
815 { 819 {
816 fCurrentInterval = i1; 820 fCurrentInterval = i1;
817 result += integral.Legendre10(this,&G4In 821 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
818 energy1,en 822 energy1,energy2);
819 fPAIphotonVector->PutValue(k,result); 823 fPAIphotonVector->PutValue(k,result);
820 824
821 } 825 }
822 else 826 else
823 { 827 {
824 for( i = i2; i >= i1; i-- ) 828 for( i = i2; i >= i1; i-- )
825 { 829 {
826 fCurrentInterval = i; 830 fCurrentInterval = i;
827 831
828 if( i==i2 ) result += integral. 832 if( i==i2 ) result += integral.Legendre10(this,
829 &G4InitXscPAI::PAId 833 &G4InitXscPAI::PAIdNdxCherenkov,
830 (*(*fMatSandiaMatri 834 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
831 835
832 else if( i == i1 ) result += integral.Legend 836 else if( i == i1 ) result += integral.Legendre10(this,
833 &G4InitXscPAI::PAId 837 &G4InitXscPAI::PAIdNdxCherenkov,energy1,
834 (*(*fMatSandiaMatri 838 (*(*fMatSandiaMatrix)[i+1])[0]);
835 839
836 else result += integral. 840 else result += integral.Legendre10(this,
837 &G4InitXscPAI::PAId 841 &G4InitXscPAI::PAIdNdxCherenkov,
838 (*(*fMatSandiaMatrix)[i 842 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
839 } 843 }
840 fPAIphotonVector->PutValue(k,result); 844 fPAIphotonVector->PutValue(k,result);
841 } 845 }
842 // G4cout<<k<<"\t"<<result<<G4endl; 846 // G4cout<<k<<"\t"<<result<<G4endl;
843 } 847 }
844 return; 848 return;
845 } // end of IntegralCerenkov 849 } // end of IntegralCerenkov
846 850
847 ////////////////////////////////////////////// 851 ////////////////////////////////////////////////////////////////////////
848 // 852 //
849 // Calculation of the PAI Plasmon integral cro 853 // Calculation of the PAI Plasmon integral cross-section
850 // fIntegralPlasmon[1] = splasmon primary ioni 854 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
851 // and fIntegralPlasmon[0] = mean plasmon loss 855 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
852 856
853 void G4InitXscPAI::IntegralPlasmon(G4double bg 857 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax)
854 { 858 {
855 G4int i,k,i1,i2; 859 G4int i,k,i1,i2;
856 G4double energy1, energy2, result = 0.; 860 G4double energy1, energy2, result = 0.;
857 861
858 fBetaGammaSq = bg2; 862 fBetaGammaSq = bg2;
859 fTmax = Tmax; 863 fTmax = Tmax;
860 864
861 delete fPAIelectronVector; << 865 if(fPAIelectronVector) delete fPAIelectronVector;
862 866
863 fPAIelectronVector = new G4PhysicsLogVector( 867 fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
864 fPAIelectronVector->PutValue(fPAIbin-1,resul 868 fPAIelectronVector->PutValue(fPAIbin-1,result);
865 869
866 for( i = fIntervalNumber - 1; i >= 0; i-- ) 870 for( i = fIntervalNumber - 1; i >= 0; i-- )
867 { 871 {
868 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 872 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
869 } 873 }
870 if (i < 0) i = 0; // Tmax should be more tha 874 if (i < 0) i = 0; // Tmax should be more than
871 // first ionisation potent 875 // first ionisation potential
872 fIntervalTmax = i; 876 fIntervalTmax = i;
873 877
874 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 878 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
875 879
876 for( k = fPAIbin - 2; k >= 0; k-- ) 880 for( k = fPAIbin - 2; k >= 0; k-- )
877 { 881 {
878 energy1 = fPAIelectronVector->GetLowEdgeEn 882 energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
879 energy2 = fPAIelectronVector->GetLowEdgeEn 883 energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
880 884
881 for( i = fIntervalTmax; i >= 0; i-- ) 885 for( i = fIntervalTmax; i >= 0; i-- )
882 { 886 {
883 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 887 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
884 } 888 }
885 if(i < 0) i = 0; 889 if(i < 0) i = 0;
886 i2 = i; 890 i2 = i;
887 891
888 for( i = fIntervalTmax; i >= 0; i-- ) 892 for( i = fIntervalTmax; i >= 0; i-- )
889 { 893 {
890 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 894 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
891 } 895 }
892 if(i < 0) i = 0; 896 if(i < 0) i = 0;
893 i1 = i; 897 i1 = i;
894 898
895 if( i1 == i2 ) 899 if( i1 == i2 )
896 { 900 {
897 fCurrentInterval = i1; 901 fCurrentInterval = i1;
898 result += integral.Legendre10(this,&G4In 902 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
899 energy1,en 903 energy1,energy2);
900 fPAIelectronVector->PutValue(k,result); 904 fPAIelectronVector->PutValue(k,result);
901 } 905 }
902 else 906 else
903 { 907 {
904 for( i = i2; i >= i1; i-- ) 908 for( i = i2; i >= i1; i-- )
905 { 909 {
906 fCurrentInterval = i; 910 fCurrentInterval = i;
907 911
908 if( i==i2 ) result += integral. 912 if( i==i2 ) result += integral.Legendre10(this,
909 &G4InitXscPAI::PAId 913 &G4InitXscPAI::PAIdNdxPlasmon,
910 (*(*fMatSandiaMatri 914 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
911 915
912 else if( i == i1 ) result += integral.Legend 916 else if( i == i1 ) result += integral.Legendre10(this,
913 &G4InitXscPAI::PAId 917 &G4InitXscPAI::PAIdNdxPlasmon,energy1,
914 (*(*fMatSandiaMatri 918 (*(*fMatSandiaMatrix)[i+1])[0]);
915 919
916 else result += integral. 920 else result += integral.Legendre10(this,
917 &G4InitXscPAI::PAId 921 &G4InitXscPAI::PAIdNdxPlasmon,
918 (*(*fMatSandiaMatrix)[i 922 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
919 } 923 }
920 fPAIelectronVector->PutValue(k,result); 924 fPAIelectronVector->PutValue(k,result);
921 } 925 }
922 // G4cout<<k<<"\t"<<result<<G4endl; 926 // G4cout<<k<<"\t"<<result<<G4endl;
923 } 927 }
924 return; 928 return;
925 } // end of IntegralPlasmon 929 } // end of IntegralPlasmon
926 930
927 931
928 ////////////////////////////////////////////// 932 /////////////////////////////////////////////////////////////////////////
929 // 933 //
930 // 934 //
931 935
932 G4double G4InitXscPAI::GetPhotonLambda( G4doub 936 G4double G4InitXscPAI::GetPhotonLambda( G4double omega )
933 { 937 {
934 G4int i ; 938 G4int i ;
935 G4double omega2, omega3, omega4, a1, a2, a3, 939 G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
936 940
937 omega2 = omega*omega ; 941 omega2 = omega*omega ;
938 omega3 = omega2*omega ; 942 omega3 = omega2*omega ;
939 omega4 = omega2*omega2 ; 943 omega4 = omega2*omega2 ;
940 944
941 for(i = 0; i < fIntervalNumber;i++) 945 for(i = 0; i < fIntervalNumber;i++)
942 { 946 {
943 if( omega < (*(*fMatSandiaMatrix)[i])[0] ) 947 if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
944 } 948 }
945 if( i == 0 ) 949 if( i == 0 )
946 { 950 {
947 G4cout<<"Warning: energy in G4InitXscPAI:: 951 G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
948 } 952 }
949 else i-- ; 953 else i-- ;
950 954
951 a1 = (*(*fMatSandiaMatrix)[i])[1]; 955 a1 = (*(*fMatSandiaMatrix)[i])[1];
952 a2 = (*(*fMatSandiaMatrix)[i])[2]; 956 a2 = (*(*fMatSandiaMatrix)[i])[2];
953 a3 = (*(*fMatSandiaMatrix)[i])[3]; 957 a3 = (*(*fMatSandiaMatrix)[i])[3];
954 a4 = (*(*fMatSandiaMatrix)[i])[4]; 958 a4 = (*(*fMatSandiaMatrix)[i])[4];
955 959
956 lambda = 1./(a1/omega + a2/omega2 + a3/omega 960 lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
957 961
958 return lambda ; 962 return lambda ;
959 } 963 }
960 964
961 ////////////////////////////////////////////// 965 /////////////////////////////////////////////////////////////////////////
962 // 966 //
963 // 967 //
964 968
965 ////////////////////////////////////////////// 969 /////////////////////////////////////////////////////////////////////////
966 // 970 //
967 // 971 //
968 972
969 G4double G4InitXscPAI::GetStepEnergyLoss( G4do 973 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step )
970 { 974 {
971 G4double loss = 0.0 ; 975 G4double loss = 0.0 ;
972 loss *= step; 976 loss *= step;
973 977
974 return loss ; 978 return loss ;
975 } 979 }
976 980
977 ////////////////////////////////////////////// 981 /////////////////////////////////////////////////////////////////////////
978 // 982 //
979 // 983 //
980 984
981 G4double G4InitXscPAI::GetStepCerenkovLoss( G4 985 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step )
982 { 986 {
983 G4double loss = 0.0 ; 987 G4double loss = 0.0 ;
984 loss *= step; 988 loss *= step;
985 989
986 return loss ; 990 return loss ;
987 } 991 }
988 992
989 ////////////////////////////////////////////// 993 /////////////////////////////////////////////////////////////////////////
990 // 994 //
991 // 995 //
992 996
993 G4double G4InitXscPAI::GetStepPlasmonLoss( G4d 997 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step )
994 { 998 {
995 999
996 1000
997 G4double loss = 0.0 ; 1001 G4double loss = 0.0 ;
998 loss *= step; 1002 loss *= step;
999 return loss ; 1003 return loss ;
1000 } 1004 }
1001 1005
1002 1006
1003 // 1007 //
1004 // end of G4InitXscPAI implementation file 1008 // end of G4InitXscPAI implementation file
1005 // 1009 //
1006 ///////////////////////////////////////////// 1010 ////////////////////////////////////////////////////////////////////////////
1007 1011
1008 1012