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