include/G4HepEm/G4HepEmPositronInteractionAnnihilation.icc

0001
0002 #include "G4HepEmPositronInteractionAnnihilation.hh"
0003
0004 #include "G4HepEmTLData.hh"
0005
0006 #include "G4HepEmElectronTrack.hh"
0007 #include "G4HepEmGammaTrack.hh"
0008 #include "G4HepEmConstants.hh"
0009 #include "G4HepEmRunUtils.hh"
0010 #include "G4HepEmMath.hh"
0011
0012 //#include <iostream>
0013
0014
0015 void G4HepEmPositronInteractionAnnihilation::Perform(G4HepEmTLData* tlData, bool isatrest) {
0016    if (isatrest) {
0017      AnnihilateAtRest(tlData);
0018    } else {
0019      AnnihilateInFlight(tlData);
0020    }
0021 }
0022
0023 void G4HepEmPositronInteractionAnnihilation::AnnihilateAtRest(G4HepEmTLData* tlData) {
0024   // compute kinematics of the first gamma (isotropic direction)
0025   const double cost = 2. * tlData->GetRNGEngine()->flat() - 1.;
0026   const double sint = std::sqrt((1. - cost)*(1. + cost));
0027   const double  phi = k2Pi * tlData->GetRNGEngine()->flat();
0028   // get 2 secondary gamma track
0029   G4HepEmTrack*    secGamma1 = tlData->AddSecondaryGammaTrack()->GetTrack();
0030   double*       secGamma1Dir = secGamma1->GetDirection();
0031   G4HepEmTrack*    secGamma2 = tlData->AddSecondaryGammaTrack()->GetTrack();
0032   double*       secGamma2Dir = secGamma2->GetDirection();
0033   secGamma1Dir[0]  = sint*std::cos(phi);
0034   secGamma1Dir[1]  = sint*std::sin(phi);
0035   secGamma1Dir[2]  = cost;
0036   // compute the kinematics of the second gamma (conservation==> -first)
0037   secGamma2Dir[0]  = -secGamma1Dir[0];
0038   secGamma2Dir[1]  = -secGamma1Dir[1];
0039   secGamma2Dir[2]  = -secGamma1Dir[2];
0040   //G4HepEmTrack* thePrimaryTrack = tlData->GetPrimaryElectronTrack()->GetTrack();
0041   //const int theParentID = thePrimaryTrack->GetID();
0042   // ekin should have been set by the caller
0043   //thePrimaryTrack->SetEKin(0.0);
0044   const int theParentID = tlData->GetPrimaryElectronTrack()->GetTrack()->GetID();
0045   secGamma1->SetEKin(kElectronMassC2);
0046   secGamma1->SetParentID(theParentID);
0047   secGamma2->SetEKin(kElectronMassC2);
0048   secGamma2->SetParentID(theParentID);
0049 }
0050
0051 void G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight(
0052     const double thePrimEkin, const double *thePrimDir, double *theGamma1Ekin, double *theGamma1Dir,
0053     double *theGamma2Ekin, double *theGamma2Dir, G4HepEmRandomEngine* rnge) {
0054   // compute kinetic limits
0055   const double tau     = thePrimEkin*kInvElectronMassC2;
0056   const double gam     = tau + 1.0;
0057   const double tau2    = tau + 2.0;
0058   const double sqgrate = std::sqrt(tau/tau2)*0.5;
0059   //
0060   const double epsmin  = 0.5 - sqgrate;
0061   const double epsmax  = 0.5 + sqgrate;
0062   const double epsqot  = epsmax/epsmin;
0063   // sampling of the energy rate of the gammas
0064   const double tau4    = tau2*tau2;
0065   double eps   = 0.0;
0066   double rfunc = 0.0;
0067   double rndArray[2];
0068   do {
0069     rnge->flatArray(2, rndArray);
0070     eps   = epsmin*G4HepEmExp(G4HepEmLog(epsqot)*rndArray[0]);
0071     rfunc = 1. - eps + (2.*gam*eps-1.)/(eps*tau4);
0072   } while( rfunc < rndArray[1]);
0073   // compute direction of the gammas
0074   const double sqg2m1 = std::sqrt(tau*tau2);
0075   const double   cost = G4HepEmMin(1., G4HepEmMax(-1., (eps*tau2-1.)/(eps*sqg2m1)));
0076   const double   sint = std::sqrt((1.+cost)*(1.-cost));
0077   const double    phi = k2Pi * rnge->flat();
0078   // kinematics of the first gamma
0079   const double initEt = thePrimEkin + 2.*kElectronMassC2;
0080   const double ekinG1 = eps*initEt;
0081   *theGamma1Ekin = ekinG1;
0082   theGamma1Dir[0] = sint*std::cos(phi);
0083   theGamma1Dir[1] = sint*std::sin(phi);
0084   theGamma1Dir[2] = cost;
0085   // use the G4HepEmRunUtils function
0086   RotateToReferenceFrame(theGamma1Dir, thePrimDir);
0087   // kinematics of the second gamma (direction <== conservation)
0088   *theGamma2Ekin = initEt-ekinG1;
0089   const double initPt = std::sqrt(thePrimEkin*(thePrimEkin+2*kElectronMassC2));
0090   const double     px = initPt*thePrimDir[0] - theGamma1Dir[0]*ekinG1;
0091   const double     py = initPt*thePrimDir[1] - theGamma1Dir[1]*ekinG1;
0092   const double     pz = initPt*thePrimDir[2] - theGamma1Dir[2]*ekinG1;
0093   const double   norm = 1.0 / std::sqrt(px*px + py*py + pz*pz);
0094   theGamma2Dir[0] = px*norm;
0095   theGamma2Dir[1] = py*norm;
0096   theGamma2Dir[2] = pz*norm;
0097 }
0098
0099 void G4HepEmPositronInteractionAnnihilation::AnnihilateInFlight(G4HepEmTLData* tlData) {
0100   // get the primary e+ track
0101   G4HepEmElectronTrack* thePrimaryElTrack = tlData->GetPrimaryElectronTrack();
0102   G4HepEmTrack* thePrimaryTrack = thePrimaryElTrack->GetTrack();
0103   double            thePrimEkin = thePrimaryTrack->GetEKin();
0104   const double*      thePrimDir = thePrimaryTrack->GetDirection();
0105
0106   G4HepEmTrack*  gTr1 = tlData->AddSecondaryGammaTrack()->GetTrack();
0107   G4HepEmTrack*  gTr2 = tlData->AddSecondaryGammaTrack()->GetTrack();
0108   double gamE1, gamE2;
0109   SampleEnergyAndDirectionsInFlight(thePrimEkin, thePrimDir, &gamE1, gTr1->GetDirection(), &gamE2, gTr2->GetDirection(), tlData->GetRNGEngine());
0110
0111   gTr1->SetEKin(gamE1);
0112   gTr2->SetEKin(gamE2);
0113   const int theParentID = tlData->GetPrimaryElectronTrack()->GetTrack()->GetID();
0114   gTr1->SetParentID(theParentID);
0115   gTr2->SetParentID(theParentID);
0116   //
0117   // set primary e+ track kinetic energy to zero ==> killed
0118   thePrimaryTrack->SetEKin(0.0);
0119 }