include/G4HepEm/G4HepEmGammaInteractionCompton.icc

0001 #include "G4HepEmGammaInteractionCompton.hh"
0002
0003 #include "G4HepEmTLData.hh"
0004 #include "G4HepEmRandomEngine.hh"
0005 #include "G4HepEmData.hh"
0006 //#include "G4HepEmMatCutData.hh"
0007
0008 #include "G4HepEmElectronTrack.hh"
0009 #include "G4HepEmGammaTrack.hh"
0010 #include "G4HepEmConstants.hh"
0011 #include "G4HepEmRunUtils.hh"
0012
0013 #include "G4HepEmMath.hh"
0014
0015 #include <iostream>
0016
0017 void G4HepEmGammaInteractionCompton::Perform(G4HepEmTLData* tlData, struct G4HepEmData* /*hepEmData*/) {
0018   G4HepEmTrack* thePrimaryTrack = tlData->GetPrimaryGammaTrack()->GetTrack();
0019   const double       thePrimGmE = thePrimaryTrack->GetEKin();
0020   // low energy limit: both for the primary gamma and secondary e-
0021   const double theLowEnergyThreshold = 0.0001; // 100 eV
0022   if (thePrimGmE<theLowEnergyThreshold) {
0023     return;
0024   }
0025   // sample the post interaction photon energy and direction (already in lab frame)
0026   // note: we might need the original photon direction so we keep it here
0027   double*        thePrimGmDir = thePrimaryTrack->GetDirection();
0028   const double theOrgGmDir[3] = {thePrimGmDir[0], thePrimGmDir[1], thePrimGmDir[2]};
0029   // the 'thePrimGmDir' will be updated
0030   const double     thePostGmE = SamplePhotonEnergyAndDirection(thePrimGmE, thePrimGmDir, theOrgGmDir, tlData->GetRNGEngine());
0031   // compute the secondary e- energy and check aganints the threshold:
0032   //  - if below threshold: simple deposit the corresponding energy
0033   //  - compute the secondary e- direction otherwise and create the secondary track
0034   const double  theSecElE = thePrimGmE-thePostGmE;
0035   // keep track of energy deposits due to particles killed below threshold
0036   double theEnergyDeposit = 0.0;
0037   if (theSecElE > theLowEnergyThreshold) {
0038     // get a secondary e- track and sample/compute directions (all will be already in lab. frame)
0039     G4HepEmTrack* theSecTrack = tlData->AddSecondaryElectronTrack()->GetTrack();
0040     double*       theSecElDir = theSecTrack->GetDirection();
0041     theSecElDir[0] = thePrimGmE * theOrgGmDir[0] - thePostGmE * thePrimGmDir[0];
0042     theSecElDir[1] = thePrimGmE * theOrgGmDir[1] - thePostGmE * thePrimGmDir[1];
0043     theSecElDir[2] = thePrimGmE * theOrgGmDir[2] - thePostGmE * thePrimGmDir[2];
0044     // normalisation
0045     const double  norm = 1.0 / std::sqrt(theSecElDir[0] * theSecElDir[0] + theSecElDir[1] * theSecElDir[1] + theSecElDir[2] * theSecElDir[2]);
0046     theSecElDir[0] *= norm;
0047     theSecElDir[1] *= norm;
0048     theSecElDir[2] *= norm;
0049     // set other properties of the secondary track as well
0050     theSecTrack->SetEKin(theSecElE);
0051     theSecTrack->SetParentID(thePrimaryTrack->GetID());
0052   } else {
0053     theEnergyDeposit += theSecElE;
0054   }
0055   //
0056   // check the post interaction gamma energy aganints the threshold:
0057   //  - if below threshold: simple deposit the corresponding energy
0058   if (thePostGmE > theLowEnergyThreshold) {
0059     thePrimaryTrack->SetEKin(thePostGmE);
0060   } else {
0061     theEnergyDeposit += thePostGmE;
0062     thePrimaryTrack->SetEKin(0.0);
0063   }
0064   thePrimaryTrack->SetEnergyDeposit(theEnergyDeposit);
0065 }
0066
0067 double G4HepEmGammaInteractionCompton::SamplePhotonEnergyAndDirection(
0068     const double thePrimGmE, double* thePrimGmDir, const double* theOrgPrimGmDir, G4HepEmRandomEngine* rnge) {
0069   // sample the post interaction reduced photon energy according to the KN DCS
0070   const double kappa = thePrimGmE * kInvElectronMassC2;
0071   const double eps0  = 1. / (1. + 2. * kappa);
0072   const double eps02 = eps0 * eps0;
0073   const double al1   = -G4HepEmLog(eps0);
0074   const double al2   = al1 + 0.5 * (1. - eps02);
0075   double eps, eps2, gf;
0076   double oneMinusCost, sint2;
0077   double rndm[3];
0078   do {
0079     rnge->flatArray(3, rndm);
0080     if (al1 > al2*rndm[0]) {
0081       eps  = G4HepEmExp(-al1 * rndm[1]);
0082       eps2 = eps * eps;
0083     } else {
0084       eps2 = eps02 + (1. - eps02) * rndm[1];
0085       eps  = std::sqrt(eps2);
0086     }
0087     oneMinusCost = (1. - eps) / (eps * kappa);
0088     sint2    = oneMinusCost * (2. - oneMinusCost);
0089     gf       = 1. - eps * sint2 / (1. + eps2);
0090   } while (gf < rndm[2]);
0091   // compute the post interaction photon direction and transform to lab frame
0092   const double cost = 1.0 - oneMinusCost;
0093   const double sint = std::sqrt(G4HepEmMax(0., sint2));
0094   const double phi  = k2Pi * rnge->flat();
0095   // direction of the scattered gamma in the scattering frame
0096   thePrimGmDir[0]   = sint * std::cos(phi);
0097   thePrimGmDir[1]   = sint * std::sin(phi);
0098   thePrimGmDir[2]   = cost;
0099   // rotate to refernce frame (G4HepEmRunUtils function) to get it in lab. frame
0100   RotateToReferenceFrame(thePrimGmDir, theOrgPrimGmDir);
0101   // return with the post interaction gamma energy
0102   return thePrimGmE*eps;
0103 }