source: trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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