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