71 : fPAIxscVector(nullptr),
72 fPAIdEdxVector(nullptr),
73 fPAIphotonVector(nullptr),
74 fPAIelectronVector(nullptr),
75 fChCosSqVector(nullptr),
76 fChWidthVector(nullptr)
97 for(j = 1; j < 5 ; j++)
138 energy1 = (*(*fMatSandiaMatrix)[i])[0];
139 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
141 if( energy2 - energy1 > 1.5*
fDelta*(energy1 + energy2) )
continue ;
144 for(j = i; j < fIntervalNumber-1; j++)
146 for( k = 0; k < 5; k++ )
175 energy1 = (*(*fMatSandiaMatrix)[i])[0];
176 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
190 for(j = 1; j < 5 ; j++)
238 a1 = (*(*fMatSandiaMatrix)[k])[1];
239 a2 = (*(*fMatSandiaMatrix)[k])[2];
240 a3 = (*(*fMatSandiaMatrix)[k])[3];
241 a4 = (*(*fMatSandiaMatrix)[k])[4];
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 ;
248 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
265 energy1 = (*(*fMatSandiaMatrix)[i])[0];
272 energy1 = (*(*fMatSandiaMatrix)[i])[0];
278 energy1 = (*(*fMatSandiaMatrix)[i])[0];
279 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
299 a1 = (*(*fMatSandiaMatrix)[k])[1];
300 a2 = (*(*fMatSandiaMatrix)[k])[2];
301 a3 = (*(*fMatSandiaMatrix)[k])[3];
302 a4 = (*(*fMatSandiaMatrix)[k])[4];
304 energy2 = energy1*energy1;
305 energy3 = energy2*energy1;
306 energy4 = energy3*energy1;
308 result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
309 result *=
hbarc/energy1 ;
331 result = eIm2 + eRe2;
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 ;
354 x1 = (*(*fMatSandiaMatrix)[i])[0];
355 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
357 a1 = (*(*fMatSandiaMatrix)[i])[1];
358 a2 = (*(*fMatSandiaMatrix)[i])[2];
359 a3 = (*(*fMatSandiaMatrix)[i])[3];
360 a4 = (*(*fMatSandiaMatrix)[i])[4];
362 if( std::abs(x0-x1) < 0.5*(x0+
x1)*
fDelta )
364 if(x0 >= x1) x0 = x1*(1+
fDelta);
367 if( std::abs(x0-x2) < 0.5*(x0+
x2)*fDelta )
369 if(x0 >= x2) x0 = x2*(1+
fDelta);
376 if( xx12 < 0 ) xx12 = -xx12;
380 xln3 = log((x2 + x0)/(x1 + x0)) ;
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 ;
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 ;
396 cof1 = a1/x02 + a3/x04 ;
397 cof2 = a2/x03 + a4/x05 ;
399 result += 0.5*(cof1 +cof2)*xln2 ;
400 result += 0.5*(cof1 - cof2)*xln3 ;
419 G4double be2,cof,
x1,
x2,x3,x4,x5,x6,x7,x8,
result ;
424 static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
425 be2 = betaGammaSq/(1 + betaGammaSq) ;
431 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
434 x2 = -log( (1/betaGammaSq - epsilonRe)*
435 (1/betaGammaSq - epsilonRe) +
436 epsilonIm*epsilonIm )/2 ;
438 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
444 x3 = -epsilonRe + 1/betaGammaSq ;
445 x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
446 epsilonIm*epsilonIm) ;
448 x7 = atan2(epsilonIm,x3) ;
453 x4 = ((x1 +
x2)*epsilonIm + x6)/
hbarc ;
455 x8 = (1 + epsilonRe)*(1 + epsilonRe) +
458 result = (x4 + cof*integralTerm/omega/omega) ;
459 if(result < 1.0
e-8) result = 1.0e-8 ;
460 result *= fine_structure_const/be2/
pi ;
463 result *= (1-exp(-be4/betaBohr4)) ;
494 G4double logarithm, x3, x5, argument, modul2, dNdxC ;
498 static const G4double cofBetaBohr = 4.0 ;
500 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
502 be2 = betaGammaSq/(1 + betaGammaSq) ;
505 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ;
508 logarithm = -log( (1/betaGammaSq - epsilonRe)*
509 (1/betaGammaSq - epsilonRe) +
510 epsilonIm*epsilonIm )*0.5 ;
511 logarithm += log(1+1.0/betaGammaSq) ;
514 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
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) ;
528 dNdxC = ( logarithm*epsilonIm + argument )/
hbarc ;
530 if(dNdxC < 1.0
e-8) dNdxC = 1.0e-8 ;
532 dNdxC *= fine_structure_const/be2/
pi ;
534 dNdxC *= (1-exp(-be4/betaBohr4)) ;
538 modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
559 G4double cof, resonance, modul2, dNdxP ;
563 static const G4double cofBetaBohr = 4.0 ;
565 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
567 be2 = betaGammaSq/(1 + betaGammaSq) ;
571 resonance *= epsilonIm/
hbarc ;
574 dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
576 if( dNdxP < 1.0
e-8 ) dNdxP = 1.0e-8 ;
578 dNdxP *= fine_structure_const/be2/
pi ;
579 dNdxP *= (1-exp(-be4/betaBohr4)) ;
583 modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
620 for( k =
fPAIbin - 2; k >= 0; k-- )
648 for( i = i2; i >= i1; i-- )
652 if( i==i2 ) result += integral.
Legendre10(
this,
656 else if( i == i1 ) result += integral.
Legendre10(
this,
701 for( k =
fPAIbin - 2; k >= 0; k-- )
729 for( i = i2; i >= i1; i-- )
733 if( i==i2 ) result += integral.
Legendre10(
this,
737 else if( i == i1 ) result += integral.
Legendre10(
this,
789 for( k =
fPAIbin - 2; k >= 0; k-- )
825 for( i = i2; i >= i1; i-- )
829 if( i==i2 ) result += integral.
Legendre10(
this,
833 else if( i == i1 ) result += integral.
Legendre10(
this,
877 for( k =
fPAIbin - 2; k >= 0; k-- )
905 for( i = i2; i >= i1; i-- )
909 if( i==i2 ) result += integral.
Legendre10(
this,
913 else if( i == i1 ) result += integral.
Legendre10(
this,
938 omega2 = omega*omega ;
939 omega3 = omega2*omega ;
940 omega4 = omega2*omega2 ;
948 G4cout<<
"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<
G4endl;
952 a1 = (*(*fMatSandiaMatrix)[i])[1];
953 a2 = (*(*fMatSandiaMatrix)[i])[2];
954 a3 = (*(*fMatSandiaMatrix)[i])[3];
955 a4 = (*(*fMatSandiaMatrix)[i])[4];
957 lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
G4PhysicsLogVector * fPAIxscVector
G4PhysicsLogVector * fPAIelectronVector
G4double GetLowEdgeEnergy(size_t binNumber) const
G4OrderedTable * fMatSandiaMatrix
G4double ModuleSqDielectricConst(G4int intervalNumber, G4double energy)
static const G4double fSolidDensity
Float_t x1[n_points_granero]
G4double DifPAIdEdx(G4double omega)
static constexpr double hbarc
G4PhysicsLogVector * fPAIphotonVector
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double GetStepPlasmonLoss(G4double step)
G4int GetMaxInterval() const
void IntegralPAIdEdx(G4double bg2, G4double Tmax)
void IntegralPlasmon(G4double bg2, G4double Tmax)
void IntegralPAIxSection(G4double bg2, G4double Tmax)
static const G4int fPAIbin
static constexpr double g
G4PhysicsLogVector * fChWidthVector
static const G4double fDelta
G4PhysicsLogVector * fPAIdEdxVector
static constexpr double electron_mass_c2
G4PhysicsLogVector * fChCosSqVector
G4double PAIdNdxPlasmon(G4double omega)
G4double fNormalizationCof
G4double GetStepEnergyLoss(G4double step)
G4double DifPAIxSection(G4double omega)
G4double RePartDielectricConst(G4double energy)
G4double G4ParticleHPJENDLHEData::G4double result
G4double GetStepCerenkovLoss(G4double step)
G4InitXscPAI(const G4MaterialCutsCouple *matCC)
G4double fElectronDensity
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
void KillCloseIntervals()
G4GLOB_DLL std::ostream G4cout
const G4Material * GetMaterial() const
static constexpr double pi
G4double GetSandiaMatTable(G4int, G4int) const
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
G4double PAIdNdxCherenkov(G4double omega)
static constexpr double cm3
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
void IntegralCherenkov(G4double bg2, G4double Tmax)
G4double GetElectronDensity() const
G4double IntegralTerm(G4double omega)
void PutValue(size_t index, G4double theValue)
G4double GetPhotonLambda(G4double omega)
G4double GetDensity() const