include/G4HepEm/G4HepEmGammaInteractionPhotoelectric.icc

0001
0002 #include "G4HepEmGammaInteractionPhotoelectric.hh"
0003
0004
0005 #include  "G4HepEmTLData.hh"
0006 #include  "G4HepEmData.hh"
0007 #include  "G4HepEmElementData.hh"
0008 #include  "G4HepEmMaterialData.hh"
0009 #include  "G4HepEmMatCutData.hh"
0010 #include  "G4HepEmRunUtils.hh"
0011 #include  "G4HepEmConstants.hh"
0012
0013 void G4HepEmGammaInteractionPhotoelectric::Perform(G4HepEmTLData* tlData, struct G4HepEmData* hepEmData) {
0014   G4HepEmGammaTrack* theGammaTrack = tlData->GetPrimaryGammaTrack();
0015   G4HepEmTrack*    thePrimaryTrack = theGammaTrack->GetTrack();
0016   const double*        theGammaDir = thePrimaryTrack->GetDirection();
0017   const double           theGammaE = thePrimaryTrack->GetEKin();
0018
0019   const int theMCIndx = thePrimaryTrack->GetMCIndex();
0020   const double  mxsec = theGammaTrack->GetPEmxSec();
0021
0022   const double bindingEnergy = SelectElementBindingEnergy(hepEmData, theMCIndx, mxsec, theGammaE, tlData->GetRNGEngine());
0023
0024   const double theLowEnergyThreshold = 0.000001; // 1 eV
0025   const double photoElecE = theGammaE - bindingEnergy;
0026   if (photoElecE > theLowEnergyThreshold) {
0027     G4HepEmTrack* theSecTrack = tlData->AddSecondaryElectronTrack()->GetTrack();
0028     double*       theSecElDir = theSecTrack->GetDirection();
0029     SamplePhotoElectronDirection(photoElecE, theGammaDir, theSecElDir, tlData->GetRNGEngine());
0030     theSecTrack->SetEKin(photoElecE);
0031     theSecTrack->SetParentID(thePrimaryTrack->GetID());
0032     thePrimaryTrack->SetEnergyDeposit(bindingEnergy);
0033   } else {
0034     thePrimaryTrack->SetEnergyDeposit(theGammaE);
0035   }
0036   thePrimaryTrack->SetEKin(0.0);
0037 }
0038
0039 double G4HepEmGammaInteractionPhotoelectric::SelectElementBindingEnergy(const struct G4HepEmData* hepEmData, const int imc, const double mxsec, const double ekin, G4HepEmRandomEngine *rnge) {
0040   const int theMatIndx = hepEmData->fTheMatCutData->fMatCutData[imc].fHepEmMatIndex;
0041   const G4HepEmMatData& theMData = hepEmData->fTheMaterialData->fMaterialData[theMatIndx];
0042
0043   // Possible optimization: if ekin minus the minimum binding energy of all elements in the material
0044   // is already smaller than the electron cut, we could skip selecting an element.
0045
0046   int ielem = 0;
0047   if (theMData.fNumOfElement > 1) {
0048     const double x = rnge->flat() * mxsec;
0049     double sum = 0;
0050     double invE = 1 / ekin;
0051     for (int i = 0; i < theMData.fNumOfElement; i++) {
0052       const G4HepEmElemData& theEData = hepEmData->fTheElementData->fElementData[theMData.fElementVect[i]];
0053       int interval = 0;
0054       if (ekin >= theEData.fSandiaEnergies[0]) {
0055         // Optimization: linear search starting with intervals for higher energies.
0056         for (int i = theEData.fNumOfSandiaIntervals - 1; i >= 0; i--) {
0057           if (ekin >= theEData.fSandiaEnergies[i]) {
0058             interval = i;
0059             break;
0060           }
0061         }
0062       }
0063       const double* sandiaCof = &theEData.fSandiaCoefficients[4 * interval];
0064       sum += theMData.fNumOfAtomsPerVolumeVect[i] *
0065         invE * (sandiaCof[0] + invE * (sandiaCof[1] + invE * (sandiaCof[2] + invE * sandiaCof[3])));
0066       if (x <= sum) {
0067         ielem = i;
0068         break;
0069       }
0070     }
0071   }
0072
0073   // Only check the k-shell binding energy of the sampled element. If the
0074   // gamma energy is below that, the energy of the secondary photo electron
0075   // will be so low that it won't travel far and we can skip generating it.
0076
0077   return hepEmData->fTheElementData->fElementData[theMData.fElementVect[ielem]].fKShellBindingEnergy;
0078 }
0079
0080 void G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(const double kinE, const double* theGammaDir, double* theDir, G4HepEmRandomEngine* rnge) {
0081   // -- Sample from SauterGavrila angular distribution, code taken from Geant4:
0082   // Initial algorithm according Penelope 2008 manual and
0083   // F.Sauter Ann. Physik 9, 217(1931); 11, 454(1931).
0084   // Modified according Penelope 2014 manual
0085   const double tau = kinE * kInvElectronMassC2;
0086   const double gamma = 1.0 + tau;
0087   const double beta = std::sqrt(tau*(tau + 2.0))/gamma;
0088
0089   // ac corresponds to "A" of Eq. (2.31)
0090   //
0091   const double ac = (1.0 - beta)/beta;
0092   const double a1 = 0.5*beta*gamma*tau*(gamma-2.0);
0093   const double a2 = ac + 2.0;
0094   // gtmax = maximum of the rejection function according to Eq. (2.28),
0095   // obtained for tsam=0
0096   const double gtmax = 2.0*(a1 + 1.0/ac);
0097
0098   double tsam = 0.0;
0099   double gtr  = 0.0;
0100
0101   //2) sampling. Eq. (2.31) of Penelope Manual
0102   // tsam = 1-std::cos(theta)
0103   // gtr = rejection function according to Eq. (2.28)
0104   double rndm[2];
0105   do {
0106     rnge->flatArray(2, rndm);
0107     tsam = 2.0*ac * (2.0*rndm[0] + a2*std::sqrt(rndm[0])) / (a2*a2 - 4.0*rndm[0]);
0108     gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
0109     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
0110   } while(rndm[1]*gtmax > gtr);
0111
0112   const double costheta = 1.0 - tsam;
0113
0114   const double sint = std::sqrt(tsam*(2.0 - tsam));
0115   const double phi  = k2Pi*rnge->flat();
0116
0117   theDir[0] = sint * std::cos(phi);
0118   theDir[1] = sint * std::sin(phi);
0119   theDir[2] = costheta;
0120   RotateToReferenceFrame(theDir, theGammaDir);
0121 }