include/G4HepEm/G4HepEmGammaInteractionConversion.icc

0001 #include "G4HepEmGammaInteractionConversion.hh"
0002
0003 #include "G4HepEmTLData.hh"
0004 #include "G4HepEmRandomEngine.hh"
0005 #include "G4HepEmData.hh"
0006 #include "G4HepEmMatCutData.hh"
0007 #include "G4HepEmMaterialData.hh"
0008 #include "G4HepEmElementData.hh"
0009 #include "G4HepEmGammaData.hh"
0010
0011 #include "G4HepEmConstants.hh"
0012 #include "G4HepEmInteractionUtils.hh"
0013 #include "G4HepEmRunUtils.hh"
0014
0015 #include "G4HepEmMath.hh"
0016
0017 #include <iostream>
0018
0019 void G4HepEmGammaInteractionConversion::Perform(G4HepEmTLData* tlData, struct G4HepEmData* hepEmData) {
0020   G4HepEmTrack* thePrimaryTrack = tlData->GetPrimaryGammaTrack()->GetTrack();
0021   const double       thePrimGmE = thePrimaryTrack->GetEKin();
0022   //
0023   // check kinematical limit: gamma energy(Eg) must be at least 2 e- rest mass
0024   // (but the model should be used at higher energies above 100 MeV)
0025   if (thePrimGmE < 2.*kElectronMassC2) {
0026     return;
0027   }
0028   //
0029   // Sample/compute kinetic energies of the e- and e+ pair
0030   const double theLogPrimGmE = thePrimaryTrack->GetLogEKin();
0031   const int        theMCIndx = thePrimaryTrack->GetMCIndex();
0032   double elKinEnergy;  // e- kinetic energy
0033   double posKinEnergy; // e+ kinetic energy
0034   SampleKinEnergies(hepEmData, thePrimGmE, theLogPrimGmE, theMCIndx, elKinEnergy, posKinEnergy, tlData->GetRNGEngine());
0035   //
0036   // Sample/compute secondary e-/e+ directions:
0037   // obtain 2 secondary electorn track (one with +1.0 charge for e+)
0038   G4HepEmTrack* theSecElTrack  = tlData->AddSecondaryElectronTrack()->GetTrack();
0039   G4HepEmTrack* theSecPosTrack = tlData->AddSecondaryElectronTrack()->GetTrack();
0040   SampleDirections(thePrimaryTrack->GetDirection(), theSecElTrack->GetDirection(), theSecPosTrack->GetDirection(),
0041                    elKinEnergy, posKinEnergy, tlData->GetRNGEngine());
0042   //
0043   // Set remaining properties of the secondary tracks
0044   theSecElTrack->SetEKin(elKinEnergy);
0045   theSecElTrack->SetParentID(thePrimaryTrack->GetID());
0046   theSecPosTrack->SetEKin(posKinEnergy);
0047   theSecPosTrack->SetParentID(thePrimaryTrack->GetID());
0048   theSecPosTrack->SetCharge(+1.0);
0049
0050   //
0051   // Kill the primary gamma track by setting its energy to zero.
0052   thePrimaryTrack->SetEKin(0.0);
0053 }
0054
0055
0056 void G4HepEmGammaInteractionConversion::SampleKinEnergies(struct G4HepEmData* hepEmData, double thePrimEkin,
0057                                                           double theLogEkin, int theMCIndx, double& eKinEnergy,
0058                                                           double& pKinEnergy, G4HepEmRandomEngine* rnge) {
0059   // get the material data
0060   const int               matIndx = (hepEmData->fTheMatCutData->fMatCutData[theMCIndx]).fHepEmMatIndex;
0061   const G4HepEmMatData&  theMData = hepEmData->fTheMaterialData->fMaterialData[matIndx];
0062   // sample target element
0063   const int  elemIndx = (theMData.fNumOfElement > 1)
0064                        ? SelectTargetAtom(hepEmData->fTheGammaData, matIndx, thePrimEkin, theLogEkin, rnge->flat())
0065                        : 0;
0066   const int      iZet = theMData.fElementVect[elemIndx];
0067   const double lpmEnr = kLPMconstant * theMData.fRadiationLength;
0068
0069   // get the corresponding element data
0070   const G4HepEmElemData& theElemData = hepEmData->fTheElementData->fElementData[G4HepEmMin(iZet, hepEmData->fTheElementData->fMaxZet)];
0071   //
0072   // == Sampling of the total energy, transferred to one of the e-/e+ pair in
0073   //    units of initial photon energy
0074   const double eps0 = kElectronMassC2/thePrimEkin;
0075   double eps = 0.0;
0076   if (thePrimEkin < 2.0) {
0077     // uniform sampling at low energies (flat DCS) between the kinematical limits
0078     // of eps0=mc^2/Eg <= eps <= 0.5 ( symmetric DCS around eps = 0.5)
0079     eps = eps0 + (0.5-eps0)*rnge->flat();
0080   } else {
0081     // use a gamma energy limit of 50.0 [MeV] to turn off Coulomb correction below
0082     const double deltaFactor = eps0*136./theElemData.fZet13;
0083     const double    deltaMin = 4.*deltaFactor;
0084     const double    deltaMax = (thePrimEkin < 50.0) ? theElemData.fDeltaMaxLow : theElemData.fDeltaMaxHigh;
0085     const double    logZ13   = 0.333333*theElemData.fLogZ;
0086     const double          FZ = (thePrimEkin < 50.0) ? 8.*logZ13 : 8.*(logZ13 + theElemData.fCoulomb);
0087     // compute the limits of eps
0088     const double        epsp = 0.5 - 0.5*std::sqrt(1. - deltaMin/deltaMax) ;
0089     const double      epsMin = G4HepEmMax(eps0, epsp);
0090     const double    epsRange = 0.5 - epsMin;
0091     //
0092     // sample the energy rate (eps) of the created electron (or positron)
0093     double F10, F20;
0094     ScreenFunction12(deltaMin, F10, F20);
0095     F10 -= FZ;
0096     F20 -= FZ;
0097     const double NormF1   = G4HepEmMax(F10 * epsRange * epsRange, 0.);
0098     const double NormF2   = G4HepEmMax(1.5 * F20                , 0.);
0099     const double NormCond = NormF1/(NormF1 + NormF2);
0100     // check if LPM correction is active ( active if gamma energy > 100 [GeV])
0101     eps = (thePrimEkin < 100000.0)
0102           ? SampleEnergyRateNoLPM  (NormCond, epsMin, epsRange, deltaFactor, 1./F10, 1./F20, FZ, rnge)
0103           : SampleEnergyRateWithLPM(NormCond, epsMin, epsRange, deltaFactor, 1./F10, 1./F20, FZ, rnge,
0104                                     thePrimEkin, lpmEnr, &theElemData);
0105   }
0106   //
0107   // select charges randomly and compute kinetic
0108   double eTotEnergy, pTotEnergy;
0109   if (rnge->flat() > 0.5) {
0110     eTotEnergy = (1.-eps)*thePrimEkin;
0111     pTotEnergy = eps*thePrimEkin;
0112   } else {
0113     pTotEnergy = (1.-eps)*thePrimEkin;
0114     eTotEnergy = eps*thePrimEkin;
0115   }
0116   //
0117   // compute kinetic energies of the e- and e+ pair.
0118   eKinEnergy = G4HepEmMax(0.,eTotEnergy - kElectronMassC2);
0119   pKinEnergy = G4HepEmMax(0.,pTotEnergy - kElectronMassC2);
0120 }
0121
0122
0123 void G4HepEmGammaInteractionConversion::SampleDirections(const double* orgGammaDir, double* secElDir,
0124                                                          double* secPosDir, const double secElEkin,
0125                                                          const double secPosEkin, G4HepEmRandomEngine* rnge) {
0126   // sample azimuthal angle (2Pi symmetric)
0127   const double  phi    = k2Pi*rnge->flat();
0128   const double cosPhi  = std::cos(phi);
0129   const double sinPhi  = std::sin(phi);
0130   // sample e- cos(theta) by using the (modified Tsai sampling:
0131   const double costEl  = SampleCostModifiedTsai(secElEkin, rnge);
0132   const double sintEl  = std::sqrt((1.0-costEl)*(1.0+costEl));
0133   secElDir[0] = sintEl * cosPhi;
0134   secElDir[1] = sintEl * sinPhi;
0135   secElDir[2] = costEl;
0136   // rotate to refernce frame (G4HepEmRunUtils function) to get it in lab. frame
0137   RotateToReferenceFrame(secElDir, orgGammaDir);
0138   // sample e+ cos(theta) by using the (modified Tsai sampling:
0139   const double costPos = SampleCostModifiedTsai(secPosEkin, rnge);
0140   const double sintPos  = std::sqrt((1.0-costPos)*(1.0+costPos));
0141   secPosDir[0] = -sintPos * cosPhi;
0142   secPosDir[1] = -sintPos * sinPhi;
0143   secPosDir[2] = costPos;
0144   // rotate to refernce frame (G4HepEmRunUtils function) to get it in lab. frame
0145   RotateToReferenceFrame(secPosDir, orgGammaDir);
0146 }
0147
0148
0149
0150 // should be called only for mat-cuts with more than one elements in their material
0151 int G4HepEmGammaInteractionConversion::SelectTargetAtom(const struct G4HepEmGammaData* gmData, const int imat,
0152                                                         const double ekin, const double lekin, const double urndn) {
0153   // start index for this material
0154   const int   indxStart = gmData->fElemSelectorConvStartIndexPerMat[imat];
0155   const double* theData = &(gmData->fElemSelectorConvData[indxStart]);
0156   const int     numData = gmData->fElemSelectorConvEgridSize;
0157   const int     numElem = theData[0]; // the very first element for each material
0158   const double    logE0 = gmData->fElemSelectorConvLogMinEkin;
0159   const double    invLD = gmData->fElemSelectorConvEILDelta;
0160   const double*   xdata = gmData->fElemSelectorConvEgrid;
0161   // make sure that $x \in  [x[0],x[ndata-1]]$
0162   const double   xv = G4HepEmMax(xdata[0], G4HepEmMin(xdata[numData-1], ekin));
0163   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0164   const int idxEkin = G4HepEmMax(0.0, G4HepEmMin((lekin-logE0)*invLD, numData-2.0));
0165   // linear interpolation
0166   const double   x1 = xdata[idxEkin];
0167   const double   x2 = xdata[idxEkin+1];
0168   const double   dl = x2-x1;
0169   const double    b = G4HepEmMax(0., G4HepEmMin(1., (xv - x1)/dl));
0170   // the real index position of the y-data: idxEkin x (numElem-1)+1 (+1 the very first #element)
0171   const int  indx0 = idxEkin*(numElem-1) + 1;
0172   const int  indx1 = indx0 + (numElem-1);
0173   int theElemIndex = 0;
0174   while (theElemIndex<numElem-1 && urndn > theData[indx0+theElemIndex]+b*(theData[indx1+theElemIndex]-theData[indx0+theElemIndex])) { ++theElemIndex; }
0175   return theElemIndex;
0176 }
0177
0178
0179 double G4HepEmGammaInteractionConversion::SampleEnergyRateNoLPM(
0180     const double normCond, const double epsMin, const double epsRange, const double deltaFactor,
0181     const double invF10, const double invF20, const double fz, G4HepEmRandomEngine* rnge) {
0182   double rndmv[3];
0183   double greject = 0.;
0184   double eps     = 0.;
0185   do {
0186     rnge->flatArray(3, rndmv);
0187     if (normCond > rndmv[0]) {
0188       eps = 0.5 - epsRange * std::pow(rndmv[1], 1./3.);//G4HepEmX13(rndmv[1]);
0189       const double delta = deltaFactor/(eps*(1.-eps));
0190       greject = (ScreenFunction1(delta)-fz)*invF10;
0191     } else {
0192       eps = epsMin + epsRange*rndmv[1];
0193       const double delta = deltaFactor/(eps*(1.-eps));
0194       greject = (ScreenFunction2(delta)-fz)*invF20;
0195     }
0196     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
0197   } while (greject < rndmv[2]);
0198   //  end of eps sampling
0199   return eps;
0200 }
0201
0202
0203 double G4HepEmGammaInteractionConversion::SampleEnergyRateWithLPM(
0204     const double normCond, const double epsMin, const double epsRange, const double deltaFactor,
0205     const double invF10, const double invF20, const double fz, G4HepEmRandomEngine* rnge,
0206     const double eGamma, const double lpmEnergy, const struct G4HepEmElemData* elemData) {
0207   const double         z23 = elemData->fZet23;
0208   const double     ilVarS1 = elemData->fILVarS1;
0209   const double ilVarS1Cond = elemData->fILVarS1Cond;
0210   double rndmv[3];
0211   double greject = 0.;
0212   double eps     = 0.;
0213   do {
0214     rnge->flatArray(3, rndmv);
0215     if (normCond > rndmv[0]) {
0216       eps = 0.5 - epsRange * std::pow(rndmv[1], 1./3.); //G4HepEmX13(rndmv[1]);
0217       const double delta = deltaFactor/(eps*(1.-eps));
0218       double funcXiS, funcGS, funcPhiS, phi1, phi2;
0219       ComputePhi12(delta, phi1, phi2);
0220       //  0.0 = no density effect correction (only in case of Brem.)
0221       // +1.0 => Brem:  s' = sqrt{ 0.125 E_lpm E_g / [ E_t ( E_t - E_g) ]  }
0222       // -1.0 => Pair:  s' = sqrt{ 0.125 E_lpm E_g / [ E_t ( E_g - E_t) ]  } with E_t = eps*E_g
0223       EvaluateLPMFunctions(funcXiS, funcGS, funcPhiS, eGamma, eps*eGamma, lpmEnergy, z23, ilVarS1, ilVarS1Cond, 0.0, -1.0);
0224       greject = funcXiS*((2.*funcPhiS+funcGS)*phi1-funcGS*phi2-funcPhiS*fz)*invF10;
0225     } else {
0226       eps = epsMin + epsRange*rndmv[1];
0227       const double delta = deltaFactor/(eps*(1.-eps));
0228       double funcXiS, funcGS, funcPhiS, phi1, phi2;
0229       ComputePhi12(delta, phi1, phi2);
0230       EvaluateLPMFunctions(funcXiS, funcGS, funcPhiS, eGamma, eps*eGamma, lpmEnergy, z23, ilVarS1, ilVarS1Cond, 0.0, -1.0);
0231       greject = funcXiS*( (funcPhiS+0.5*funcGS)*phi1 + 0.5*funcGS*phi2
0232                          -0.5*(funcGS+funcPhiS)*fz)*invF20;
0233     }
0234   } while (greject < rndmv[2]);
0235   //  end of eps sampling
0236   return eps;
0237 }
0238
0239
0240 void G4HepEmGammaInteractionConversion::ComputePhi12(const double delta, double &phi1, double &phi2) {
0241    if (delta > 1.4) {
0242      phi1 = 21.0190 - 4.145*G4HepEmLog(delta + 0.958);
0243      phi2 = phi1;
0244    } else {
0245      phi1 = 20.806 - delta*(3.190 - 0.5710*delta);
0246      phi2 = 20.234 - delta*(2.126 - 0.0903*delta);
0247    }
0248 }
0249
0250
0251 // Compute the value of the screening function 3*PHI1(delta) - PHI2(delta):
0252 double G4HepEmGammaInteractionConversion::ScreenFunction1(const double delta) {
0253  return (delta > 1.4) ? 42.038 - 8.29*G4HepEmLog(delta + 0.958)
0254                       : 42.184 - delta*(7.444 - 1.623*delta);
0255 }
0256
0257
0258 // Compute the value of the screening function 1.5*PHI1(delta) +0.5*PHI2(delta):
0259 double G4HepEmGammaInteractionConversion::ScreenFunction2(const double delta) {
0260  return (delta > 1.4) ? 42.038 - 8.29*G4HepEmLog(delta + 0.958)
0261                       : 41.326 - delta*(5.848 - 0.902*delta);
0262 }
0263
0264
0265 // Same as ScreenFunction1 and ScreenFunction2 but computes them at once
0266 void G4HepEmGammaInteractionConversion::ScreenFunction12(const double delta, double &f1, double &f2) {
0267  if (delta > 1.4) {
0268    f1 = 42.038 - 8.29*G4HepEmLog(delta + 0.958);
0269    f2 = f1;
0270  } else {
0271    f1 = 42.184 - delta*(7.444 - 1.623*delta);
0272    f2 = 41.326 - delta*(5.848 - 0.902*delta);
0273  }
0274 }