include/G4HepEm/G4HepEmElectronInteractionIoni.icc

0001
0002 #include "G4HepEmElectronInteractionIoni.hh"
0003
0004 #include "G4HepEmTLData.hh"
0005 #include "G4HepEmRandomEngine.hh"
0006 #include "G4HepEmData.hh"
0007 #include "G4HepEmMatCutData.hh"
0008
0009 #include "G4HepEmElectronTrack.hh"
0010 #include "G4HepEmConstants.hh"
0011 #include "G4HepEmRunUtils.hh"
0012 #include "G4HepEmMath.hh"
0013
0014
0015
0016 #include <iostream>
0017
0018
0019 void G4HepEmElectronInteractionIoni::Perform(G4HepEmTLData* tlData, struct G4HepEmData* hepEmData, bool iselectron) {
0020   G4HepEmElectronTrack* thePrimaryElTrack = tlData->GetPrimaryElectronTrack();
0021   G4HepEmTrack* thePrimaryTrack = thePrimaryElTrack->GetTrack();
0022   double    thePrimEkin = thePrimaryTrack->GetEKin();
0023   const int   theMCIndx = thePrimaryTrack->GetMCIndex();
0024   const double theElCut = hepEmData->fTheMatCutData->fMatCutData[theMCIndx].fSecElProdCutE;
0025
0026   const double maxETransfer = (iselectron) ? 0.5*thePrimEkin : thePrimEkin;
0027   if (maxETransfer <= theElCut) return;
0028
0029   //
0030   // sample energy transfer and compute direction
0031   const double  deltaEkin = (iselectron)
0032                             ? SampleETransferMoller(theElCut, thePrimEkin, tlData->GetRNGEngine())
0033                             : SampleETransferBhabha(theElCut, thePrimEkin, tlData->GetRNGEngine());
0034   // get a secondary e- track and sample/compute directions (all will be already in lab. frame)
0035   G4HepEmTrack* theSecTrack = tlData->AddSecondaryElectronTrack()->GetTrack();
0036   double*     theSecElecDir = theSecTrack->GetDirection();
0037   double*    thePrimElecDir = thePrimaryTrack->GetDirection();
0038   //
0039   // == Sample/copute the emitted secondary e- and post interaction primary e-/e+ directions
0040   SampleDirections(thePrimEkin, deltaEkin, theSecElecDir, thePrimElecDir, tlData->GetRNGEngine());
0041   //
0042   // == Update primary e-/e+ and the secondary e- track properties
0043   //    note: the directions are already set in SampleDirections
0044   thePrimaryTrack->SetEKin(thePrimEkin - deltaEkin);
0045   theSecTrack->SetEKin(deltaEkin);
0046   theSecTrack->SetParentID(thePrimaryTrack->GetID());
0047 }
0048
0049
0050 double G4HepEmElectronInteractionIoni::SampleETransferMoller(const double elCut, const double primEkin,
0051                                                              G4HepEmRandomEngine* rnge) {
0052   const double tmin    = elCut;
0053   const double tmax    = 0.5*primEkin;
0054   const double xmin    = tmin / primEkin;
0055   const double xmax    = tmax / primEkin;
0056   const double gamma   = primEkin * kInvElectronMassC2 + 1.0;
0057   const double gamma2  = gamma * gamma;
0058   const double xminmax = xmin * xmax;
0059   // Moller (e-e-) scattering
0060   const double gg      = (2.0 * gamma - 1.0) / gamma2;
0061   const double y       = 1. - xmax;
0062   const double gf      = 1.0 - gg * xmax + xmax * xmax * (1.0 - gg + (1.0 - gg * y) / (y * y));
0063   //
0064   double dum;
0065   double rndArray[2];
0066   double deltaEkin  = 0.;
0067   do {
0068     rnge->flatArray(2, rndArray);
0069     deltaEkin       = xminmax / (xmin * (1.0 - rndArray[0]) + xmax * rndArray[0]);
0070     const double xx = 1.0 - deltaEkin;
0071     dum             = 1.0 - gg * deltaEkin + deltaEkin * deltaEkin * (1.0 - gg + (1.0 - gg * xx) / (xx * xx));
0072   } while (gf * rndArray[1] > dum);
0073   return deltaEkin * primEkin;
0074 }
0075
0076 double G4HepEmElectronInteractionIoni::SampleETransferBhabha(const double elCut, const double primEkin,
0077                                                              G4HepEmRandomEngine* rnge) {
0078   const double tmin    = elCut;
0079   const double tmax    = primEkin;
0080   const double xmin    = tmin / primEkin;
0081   const double xmax    = tmax / primEkin;
0082   const double gamma   = primEkin * kInvElectronMassC2 + 1.0;
0083   const double gamma2  = gamma * gamma;
0084   const double beta2   = 1. - 1. / gamma2;
0085   const double xminmax = xmin * xmax;
0086   // Bhabha (e+e-) scattering
0087   const double y       = 1.0 / (1.0 + gamma);
0088   const double y2      = y * y;
0089   const double y12     = 1.0 - 2.0 * y;
0090   const double b1      = 2.0 - y2;
0091   const double b2      = y12 * (3.0 + y2);
0092   const double y122    = y12 * y12;
0093   const double b4      = y122 * y12;
0094   const double b3      = b4 + y122;
0095   const double xmax2   = xmax * xmax;
0096   const double gf      = 1.0 + (xmax2 * b4 - xmin * xmin * xmin * b3 + xmax2 * b2 - xmin * b1) * beta2;
0097   //
0098   double dum;
0099   double rndArray[2];
0100   double deltaEkin  = 0.;
0101   do {
0102     rnge->flatArray(2, rndArray);
0103     deltaEkin       = xminmax / (xmin * (1.0 - rndArray[0]) + xmax * rndArray[0]);
0104     const double xx = deltaEkin * deltaEkin;
0105     dum             = 1.0 + (xx * xx * b4 - deltaEkin * xx * b3 + xx * b2 - deltaEkin * b1) * beta2;
0106   } while (gf * rndArray[1] > dum);
0107   return deltaEkin * primEkin;
0108 }
0109
0110
0111 void G4HepEmElectronInteractionIoni::SampleDirections(const double thePrimEkin, const double deltaEkin,
0112                                                       double* theSecElecDir, double* thePrimElecDir,
0113                                                       G4HepEmRandomEngine* rnge) {
0114     const double elInitETot = thePrimEkin + kElectronMassC2;
0115     const double elInitPTot = std::sqrt(thePrimEkin * (elInitETot + kElectronMassC2));
0116     const double  deltaPTot = std::sqrt(deltaEkin * (deltaEkin + 2.0 * kElectronMassC2));
0117     const double       cost = deltaEkin * (elInitETot + kElectronMassC2) / (deltaPTot * elInitPTot);
0118     // check cosTheta limit
0119     const double   cosTheta = G4HepEmMax(-1.0, G4HepEmMin(cost, 1.0));
0120     const double   sinTheta = std::sqrt((1.0 - cosTheta) * (1.0 + cosTheta));
0121     const double        phi = k2Pi * rnge->flat();     // spherical symmetry
0122     //
0123     theSecElecDir[0]  = sinTheta * std::cos(phi);
0124     theSecElecDir[1]  = sinTheta * std::sin(phi);
0125     theSecElecDir[2]  = cosTheta;
0126     // rotate to refernce frame (G4HepEmRunUtils function) to get it in lab. frame
0127     RotateToReferenceFrame(theSecElecDir, thePrimElecDir);
0128     // go for the post-interaction primary electron/positiorn direction in lab. farme
0129     // (compute from momentum vector conservation)
0130     thePrimElecDir[0] = elInitPTot * thePrimElecDir[0] - deltaPTot * theSecElecDir[0];
0131     thePrimElecDir[1] = elInitPTot * thePrimElecDir[1] - deltaPTot * theSecElecDir[1];
0132     thePrimElecDir[2] = elInitPTot * thePrimElecDir[2] - deltaPTot * theSecElecDir[2];
0133     // normalisation
0134     const double  norm = 1.0 / std::sqrt(thePrimElecDir[0] * thePrimElecDir[0] + thePrimElecDir[1] * thePrimElecDir[1] + thePrimElecDir[2] * thePrimElecDir[2]);
0135     thePrimElecDir[0] *= norm;
0136     thePrimElecDir[1] *= norm;
0137     thePrimElecDir[2] *= norm;
0138 }