include/G4HepEm/G4HepEmInteractionUtils.icc

0001
0002 #include "G4HepEmInteractionUtils.hh"
0003
0004 #include "G4HepEmRandomEngine.hh"
0005 #include "G4HepEmConstants.hh"
0006 #include "G4HepEmMath.hh"
0007
0008 #include <cmath>
0009
0010 G4HepEmHostDevice
0011 double SampleCostModifiedTsai(const double thePrimEkin, G4HepEmRandomEngine* rnge) {
0012   // sample photon direction (modified Tsai sampling):
0013   const double uMax = 2.0*(1.0 + thePrimEkin*kInvElectronMassC2);
0014   double rndm3[3];
0015   double u;
0016   do {
0017     rnge->flatArray(3, rndm3);
0018     const double uu = -G4HepEmLog(rndm3[0]*rndm3[1]);
0019     u = (0.25 > rndm3[2]) ? uu*1.6 : uu*0.533333333;
0020   } while (u > uMax);
0021   // cost = 1.0 - 2.0*u*u/(uMax*uMax);
0022   return 1.0 - 2.0*u*u/(uMax*uMax);
0023 }
0024
0025 // times = 1.0 for Brem and -1.0 for Pair production
0026 // densityCor = 0.0  for Pair production
0027 G4HepEmHostDevice
0028 void EvaluateLPMFunctions(double& funcXiS, double& funcGS, double& funcPhiS, const double egamma,
0029      const double etotal, const double elpm, const double z23, const double ilVarS1,
0030      const double ilVarS1Cond, const double densityCor, const double times) {
0031   const double     sqrt2 = 1.414213562373095;
0032   const double redegamma = egamma / etotal;
0033   const double varSprime = std::sqrt( 0.125 * redegamma * elpm / ( times*( 1.0 - redegamma ) * etotal ) );
0034   const double     varS1 = z23 / ( 184.15 * 184.15 );
0035   const double condition = sqrt2*varS1;
0036   double funcXiSprime = 2.0;
0037   if (varSprime > 1.0) {
0038     funcXiSprime = 1.0;
0039   } else if (varSprime > condition) {
0040     const double funcHSprime = G4HepEmLog(varSprime)*ilVarS1Cond;
0041     funcXiSprime = 1.0 + funcHSprime - 0.08*(1.0-funcHSprime)*funcHSprime*(2.0-funcHSprime)*ilVarS1Cond;
0042   }
0043   funcXiS = funcXiSprime;
0044   const double    varS = varSprime / std::sqrt( funcXiSprime );
0045   // - include dielectric suppression effect into s according to Migdal (only in case of Brem !)
0046   double varShat = varS;
0047   if (densityCor != 0.0) {
0048     varShat *= ( 1.0 + densityCor / (egamma*egamma) );
0049     funcXiS = 2.0;
0050     if (varShat > 1.0) {
0051       funcXiS = 1.0;
0052     } else if (varShat > varS1) {
0053       funcXiS = 1.0 + G4HepEmLog ( varShat ) * ilVarS1;
0054     }
0055   }
0056   // avluate the LPM G(s) and Phi(s) function (approximations) at s = s-hat
0057   const double lpmSLimit =  2.0;
0058   const double lpmISDelt = 20.0;
0059   if (varShat < lpmSLimit) {
0060     double  val = varShat*lpmISDelt;
0061     int    ilow = (int)val;
0062     val        -= ilow;
0063     ilow       *= 2;
0064     funcGS      = ( kFuncLPM[ilow+2] - kFuncLPM[ilow]   ) * val + kFuncLPM[ilow];
0065     funcPhiS    = ( kFuncLPM[ilow+3] - kFuncLPM[ilow+1] ) * val + kFuncLPM[ilow+1];
0066   } else {
0067     double ss = 1.0/(varShat*varShat);
0068     ss *= ss;
0069     funcGS   = 1.0-0.0230655*ss;
0070     funcPhiS = 1.0-0.01190476*ss;
0071   }
0072   //MAKE SURE SUPPRESSION IS SMALLER THAN 1: due to Migdal's approximation on xi
0073   if (funcXiS*funcPhiS > 1.0 || varShat > 0.57) {
0074     funcXiS = 1.0/funcPhiS;
0075   }
0076 }