EIC code displayed by LXR master/​include/​G4HepEm/​G4HepEmElectronManager.icc	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-07-13 08:21:04

 
 #include "G4HepEmElectronManager.hh"
 
 #include "G4HepEmData.hh"
 #include "G4HepEmParameters.hh"
 #include "G4HepEmTLData.hh"
 #include "G4HepEmRandomEngine.hh"
 
 #include "G4HepEmConstants.hh"
 #include "G4HepEmMatCutData.hh"
 #include "G4HepEmMaterialData.hh"
 #include "G4HepEmElectronData.hh"
 
 #include "G4HepEmMath.hh"
 
 #include "G4HepEmRunUtils.hh"
 #include "G4HepEmTrack.hh"
 #include "G4HepEmElectronTrack.hh"
 #include "G4HepEmMSCTrackData.hh"
 #include "G4HepEmGammaTrack.hh"
 #include "G4HepEmElectronInteractionIoni.hh"
 #include "G4HepEmElectronInteractionBrem.hh"
 #include "G4HepEmElectronEnergyLossFluctuation.hh"
 #include "G4HepEmElectronInteractionUMSC.hh"
 #include "G4HepEmPositronInteractionAnnihilation.hh"
 
 // tlData GetPrimaryElectronTrack needs to be set needs to be set based on the G4Track;
 
 //#define NOMSC
 //#define NOFLUCTUATION
 
 
 // Note: pStepLength will be set here i.e. this is the first access to it that
 //       will clear the previous step value.
 void G4HepEmElectronManager::HowFar(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmTLData* tlData) {
   G4HepEmElectronTrack* theElTrack = tlData->GetPrimaryElectronTrack();
   G4HepEmTrack* theTrack = theElTrack->GetTrack();
   // Sample the `number-of-interaction-left`
   for (int ip=0; ip<4; ++ip) {
     if (theTrack->GetNumIALeft(ip)<=0.) {
       theTrack->SetNumIALeft(-G4HepEmLog(tlData->GetRNGEngine()->flat()), ip);
     }
   }
   HowFar(hepEmData, hepEmPars, theElTrack, tlData->GetRNGEngine());
 }
 
 
 void G4HepEmElectronManager::HowFarToDiscreteInteraction(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack) {
   int indxWinnerProcess = -1;  // init to continous
   // === 1. Continuous energy loss limit
   double pStepLength     = kALargeValue;
   G4HepEmTrack* theTrack = theElTrack->GetTrack();
   const double   theEkin = theTrack->GetEKin();
   const double  theLEkin = theTrack->GetLogEKin();
   const int       theIMC = theTrack->GetMCIndex();
   const bool  isElectron = (theTrack->GetCharge() < 0.0);
 
   const G4HepEmElectronData* theElectronData = isElectron
                                                ? hepEmData->fTheElectronData
                                                : hepEmData->fThePositronData;
   //
   const double range  = GetRestRange(theElectronData, theIMC, theEkin, theLEkin);
   theElTrack->SetRange(range);
   const int indxRegion = hepEmData->fTheMatCutData->fMatCutData[theIMC].fG4RegionIndex;
   const double frange  = hepEmPars->fParametersPerRegion[indxRegion].fFinalRange;
   const double drange  = hepEmPars->fParametersPerRegion[indxRegion].fDRoverRange;
   pStepLength = (range > frange)
                 ? range*drange + frange*(1.0-drange)*(2.0-frange/range)
                 : range;
 //  std::cout << " pStepLength = " << pStepLength << " range = " << range << " frange = " << frange << std::endl;
   // === 2. Discrete limits due to eestricted Ioni and Brem (accounting e-loss)
   const int theImat = (hepEmData->fTheMatCutData->fMatCutData[theIMC]).fHepEmMatIndex;
   double mxSecs[4];
   // ioni, brem and annihilation to 2 gammas (only for e+), electron/positron nuclear
   mxSecs[0] = GetRestMacXSecForStepping(theElectronData, theIMC, theEkin, theLEkin, true);
   mxSecs[1] = GetRestMacXSecForStepping(theElectronData, theIMC, theEkin, theLEkin, false);
   mxSecs[2] = (isElectron)
               ? 0.0
               : ComputeMacXsecAnnihilationForStepping(theEkin, hepEmData->fTheMaterialData->fMaterialData[theImat].fElectronDensity);
   mxSecs[3] = GetMacXSecNuclearForStepping(theElectronData, theImat, theEkin, theLEkin);
   // compute mfp and see if we need to sample the `number-of-interaction-left`
   // before we use it to get the current discrete proposed step length
   for (int ip=0; ip<4; ++ip) {
     const double mxsec = mxSecs[ip];
     const double   mfp = (mxsec>0.) ? 1./mxsec : kALargeValue;
     // save the mac-xsec for the update of the `number-of-interaction-left`:
     // the `number-of-intercation-left` should be updated in the along-step-action
     // after the MSC has changed the step.
     theTrack->SetMFP(mfp, ip);
     // sample the proposed step length
     const double dStepLimit = mfp*theTrack->GetNumIALeft(ip);
     if (dStepLimit<pStepLength) {
       pStepLength = dStepLimit;
       indxWinnerProcess = ip;
     }
   }
 
   theElTrack->SetPStepLength(pStepLength);
   theTrack->SetWinnerProcessIndex(indxWinnerProcess);
   theTrack->SetGStepLength(pStepLength);
 }
 
 void G4HepEmElectronManager::HowFarToMSC(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack, G4HepEmRandomEngine* rnge) {
   //
   // Now MSC is called to see:
   // - if it limits the (true, i.e. physical) step length further
   // Then we perform the physical --> geometric step length conversion here:
   //  - it provides the projection of the transport vector along the original
   //    direction.
   // Note, this later also part of the MSC model and might also limit the true
   // step length further (see my note inside ConvertTrueToGeometricLength though).
   // Therefore, the check if MSC limited the step must be done after
   // the physical -->  geometric (i.e. true to geom.) conversion.
   //
 #ifndef NOMSC
   double pStepLength     = theElTrack->GetPStepLength();
   const double range     = theElTrack->GetRange();
   G4HepEmTrack* theTrack = theElTrack->GetTrack();
   const double   theEkin = theTrack->GetEKin();
   const double  theLEkin = theTrack->GetLogEKin();
   const int       theIMC = theTrack->GetMCIndex();
   const bool  isElectron = (theTrack->GetCharge() < 0.0);
 
   const G4HepEmElectronData* theElectronData = isElectron
                                                ? hepEmData->fTheElectronData
                                                : hepEmData->fThePositronData;
 
   const G4HepEmMCCData& theMatCutData = hepEmData->fTheMatCutData->fMatCutData[theIMC];
   const int theImat = theMatCutData.fHepEmMatIndex;
   const int theIreg = theMatCutData.fG4RegionIndex;
   G4HepEmMSCTrackData* mscData = theElTrack->GetMSCTrackData();
   // init some mscData for the case if we skipp calling msc due to very small step
   mscData->fTrueStepLength      = pStepLength;
   mscData->fZPathLength         = pStepLength;
   mscData->fIsActive            = false;
   mscData->SetDisplacement(0., 0., 0.);
   mscData->SetNewDirection(0., 0., 1.);
   // no msc in case of very small steps
   const double kGeomMinLength = 5.E-8; // 0.05 [nm]
   if (pStepLength > kGeomMinLength && theEkin > 1.0E-3) {
     mscData->fIsActive = true;
     // compute the fist transport mean free path
     mscData->fLambtr1  = GetTransportMFP(theElectronData, theImat, theEkin, theLEkin);
     G4HepEmElectronInteractionUMSC::StepLimit(hepEmData, hepEmPars, mscData, theEkin, theImat, theIreg, range,
                                               theTrack->GetSafety(), theTrack->GetOnBoundary(), isElectron, rnge);
     // If msc limited the true step length, then the G4HepEmMSCTrackData::fTrueStepLength member of
     // the input electron track is < pStepLengt. Otherwise its = pStepLengt.
     // Call the True --> Geometric conversion since that might limits further the true step Length:
     //   - convert the physical step length to geometrical one. The result will be
     //     written into mscData::fZPathLength.
     ConvertTrueToGeometricLength(hepEmData, mscData, theEkin, range, theIMC, isElectron);
     // check now if msc limited the step:
     const double mscTruStepLength = mscData->fTrueStepLength;
     if (mscTruStepLength < pStepLength) {
       // indicate continuous step limit as msc limited the step and set the new pStepLength
       theTrack->SetWinnerProcessIndex(-2);
       pStepLength = mscTruStepLength;
       theElTrack->SetPStepLength(pStepLength);
     }
     // set geometrical step length (protect agains wrong conversion, i.e. if gL > pL)
     theTrack->SetGStepLength(G4HepEmMin(mscData->fZPathLength, pStepLength));
   }
 #endif
 }
 
 void G4HepEmElectronManager::HowFar(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack, G4HepEmRandomEngine* rnge) {
   HowFarToDiscreteInteraction(hepEmData, hepEmPars, theElTrack);
   HowFarToMSC(hepEmData, hepEmPars, theElTrack, rnge);
 }
 
 void G4HepEmElectronManager::UpdatePStepLength(G4HepEmElectronTrack* theElTrack) {
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   const double gStepLength = theTrack->GetGStepLength();
   double pStepLength       = gStepLength;
   // call MSC::ConvertGeometricToTrueLength that will provide the true (i.e. physical)
   // step length in the G4HepEmMSCTrackData::fTrueStepLength member.
   // NOTE: in case the step was NOT limited by boundary, we know the true step length since
   //       the particle went as far as we expected.
 #ifndef NOMSC
   const double        theRange = theElTrack->GetRange();
   G4HepEmMSCTrackData* mscData = theElTrack->GetMSCTrackData();
   if (mscData->fIsActive) {
     pStepLength = mscData->fTrueStepLength;
     // if we hit boundary or stopped before we wanted for any reasons: convert geom. -> true
     if (gStepLength < mscData->fZPathLength) {
       // the converted geom --> true step Length will be written into mscData::fTrueStepLength
       ConvertGeometricToTrueLength(mscData, theRange, gStepLength);
       // protect against wrong true --> geom --> true conevrsion: physical step
       // cannot be longer than before converted to geometrical
       pStepLength = G4HepEmMin(pStepLength, mscData->fTrueStepLength);
       // store the final true step length value
       mscData->fTrueStepLength = pStepLength;
     }
     // optimisation: do not sample msc and dispalcement in case of last (rangeing out) or short steps
     const double kGeomMinLength = 5.E-8; // 0.05 [nm]
     if (pStepLength <= kGeomMinLength || theRange <= pStepLength) {
       mscData->fIsActive = false;
     }
   }
 #endif
   // set the results of the geom ---> true in the primary e- etrack
   theElTrack->SetPStepLength(pStepLength);
 }
 
 void G4HepEmElectronManager::UpdateNumIALeft(G4HepEmElectronTrack* theElTrack) {
   const double pStepLength = theElTrack->GetPStepLength();
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   double*    numInterALeft = theTrack->GetNumIALeft();
   double*       preStepMFP = theTrack->GetMFP();
   numInterALeft[0] -= pStepLength/preStepMFP[0];
   numInterALeft[1] -= pStepLength/preStepMFP[1];
   numInterALeft[2] -= pStepLength/preStepMFP[2];
   numInterALeft[3] -= pStepLength/preStepMFP[3];
 }
 
 bool G4HepEmElectronManager::ApplyMeanEnergyLoss(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack) {
   const double pStepLength = theElTrack->GetPStepLength();
 
   G4HepEmTrack* theTrack = theElTrack->GetTrack();
   const bool  isElectron = (theTrack->GetCharge() < 0.0);
   const double   theEkin = theTrack->GetEKin();
   const double  theRange = theElTrack->GetRange();
   // 0. stop tracking when reached the end (i.e. it has been ranged out by the limit)
   // @TODO: actually the tracking cut is around 1 keV and the min-table energy is 100 eV so the second should never
   //        under standard EM constructor configurations
   if (pStepLength >= theRange || theEkin <= hepEmPars->fMinLossTableEnergy) {
     // stop and deposit the remaining energy
     theTrack->SetEnergyDeposit(theEkin);
     theTrack->SetEKin(0.0);
     return true;
   }
   // 1. try linear energy loss approximation:
   const G4HepEmElectronData* elData = isElectron
                                       ? hepEmData->fTheElectronData
                                       : hepEmData->fThePositronData;
    // NOTE: this is the pre-step IMC !!!
   const int      theIMC = theTrack->GetMCIndex();
   const double theLEkin = theTrack->GetLogEKin();
   double eloss = pStepLength*GetRestDEDX(elData, theIMC, theEkin, theLEkin);
   // 2. use integral if linear energy loss is over the limit fraction
   const int indxRegion = hepEmData->fTheMatCutData->fMatCutData[theIMC].fG4RegionIndex;
   const double parLinELossLimit = hepEmPars->fParametersPerRegion[indxRegion].fLinELossLimit;
   if (eloss > theEkin*parLinELossLimit) {
     const double postStepRange = theRange - pStepLength;
     eloss = theEkin - GetInvRange(elData, theIMC, postStepRange);
   }
   eloss = G4HepEmMax(eloss, 0.0);
   if (eloss >= theEkin) {
     eloss = theEkin;
     theTrack->SetEKin(0);
     theTrack->SetEnergyDeposit(eloss);
     return true;
   }
   // 3. Compute the energy after the mean energy loss.
   double theEkinAfterMeanEloss = theEkin - eloss;
   theTrack->SetEKin(theEkinAfterMeanEloss);
   theTrack->SetEnergyDeposit(eloss);
   return false;
 }
 
 void G4HepEmElectronManager::SampleMSC(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack, G4HepEmRandomEngine* rnge) {
 #ifndef NOMSC
   const double pStepLength = theElTrack->GetPStepLength();
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   const bool    isElectron = (theTrack->GetCharge() < 0.0);
   const int         theIMC = theTrack->GetMCIndex();
   const double preStepEkin = theElTrack->GetPreStepEKin();
   const double    theRange = theElTrack->GetRange();
 
   const G4HepEmElectronData* elData = isElectron
                                       ? hepEmData->fTheElectronData
                                       : hepEmData->fThePositronData;
 
   //
   //
   // Complete here the MSC part by computing the net angular deflection and dispalcement
   //
   // Smaple scattering in MSC and compute the new direction and displacement vectors (if any)
   // The new direction and dispalcement vectors, proposed by MSC, are given in mscData::fDirection and
   // mscData::fDisplacement.
   const double kTLimitMinfix = 1.0E-8; // 0.01 [nm] 1.0E-8 [mm]
   const double kTauSmall     = 1.0e-16;
   G4HepEmMSCTrackData* mscData = theElTrack->GetMSCTrackData();
   if (mscData->fIsActive && (pStepLength > G4HepEmMax(kTLimitMinfix, kTauSmall*mscData->fLambtr1))) {
     // only to make sure that we also use E2 = E1 under the same condition as in G4
     double postStepEkin  = preStepEkin;
     double postStepLEkin = theElTrack->GetPreStepLogEKin();
 #if G4VERSION_NUM >= 1100
     if (pStepLength > theRange*0.01)
 #endif
     {
       // At this point, we have only computed the mean energy loss. Query the
       // track and if energy loss fluctuations are *not* active, we safe a log
       // call in the next step.
       postStepEkin  = theTrack->GetEKin();
       postStepLEkin = theTrack->GetLogEKin();
     }
     // sample msc scattering:
     // - compute the fist transport mean free path at the post-step energy point
     const int           theImat = (hepEmData->fTheMatCutData->fMatCutData[theIMC]).fHepEmMatIndex;
     const double postStepTr1mfp = GetTransportMFP(elData, theImat, postStepEkin, postStepLEkin);
     // - sample scattering: including net angular deflection and lateral dispacement that will be
     //                      written into mscData::fDirection and mscData::fDisplacement
     const bool isPosCor   = hepEmPars->fIsMSCPositronCor;
     const bool isDisplace = hepEmPars->fIsMSCDisplacement;
     G4HepEmElectronInteractionUMSC::SampleScattering(hepEmData, mscData, pStepLength, preStepEkin, mscData->fLambtr1, postStepEkin, postStepTr1mfp,
                                     theImat, isElectron, isPosCor, isDisplace, rnge);
     // NOTE: displacement will be applied in the caller where we have access to the required Geant4 functionality
     //       (and if its length is longer than a small minimal length and we are not ended up on boundary)
     //
     // rotate direction and displacement vectors (if any) and update new direction of the primary
     if (!(mscData->fIsNoScatteringInMSC)) {
       RotateToReferenceFrame(mscData->fDirection, theTrack->GetDirection());
       if (mscData->fIsDisplace) {
         RotateToReferenceFrame(mscData->fDisplacement, theTrack->GetDirection());
       }
       // upadte new direction
       theTrack->SetDirection(mscData->fDirection);
     }
   }
 #endif
 }
 
 bool G4HepEmElectronManager::SampleLossFluctuations(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack, G4HepEmRandomEngine* rnge) {
   const double pStepLength = theElTrack->GetPStepLength();
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   const bool    isElectron = (theTrack->GetCharge() < 0.0);
   const int         theIMC = theTrack->GetMCIndex();
 
   const double thePreStepEkin = theElTrack->GetPreStepEKin();
 
   // Previously ApplyMeanEnergyLoss computed the mean energy loss and stored the
   // result into the track.
   double finalEkin = theTrack->GetEKin();
   double eloss     = theTrack->GetEnergyDeposit();
   // sample energy loss fluctuations
 #ifndef NOFLUCTUATION
   const int iregion = hepEmData->fTheMatCutData->fMatCutData[theIMC].fG4RegionIndex;;
   const int isFluctuation = hepEmPars->fParametersPerRegion[iregion].fIsELossFluctuation;
   const double kFluctParMinEnergy  = 1.E-5; // 10 eV
   if (isFluctuation && eloss > kFluctParMinEnergy) {
     const G4HepEmMCCData& theMatCutData = hepEmData->fTheMatCutData->fMatCutData[theIMC];
     const double elCut   = theMatCutData.fSecElProdCutE;
     const int    theImat = theMatCutData.fHepEmMatIndex;
     const double meanExE = hepEmData->fTheMaterialData->fMaterialData[theImat].fMeanExEnergy;
     //
     const double tmax = isElectron ? 0.5*thePreStepEkin : thePreStepEkin;
     const double tcut = G4HepEmMin(elCut, tmax);
     eloss = G4HepEmElectronEnergyLossFluctuation::SampleEnergyLossFLuctuation(thePreStepEkin, tcut, tmax,
                                                   meanExE, pStepLength, eloss, rnge);
     eloss = G4HepEmMax(eloss, 0.0);
     // Update the final kinetic energy after loss fluctuations.
     finalEkin = thePreStepEkin - eloss;
   }
 #endif
   //
   // Check if the final kinetic energy drops below the tracking cut and stop.
   if (finalEkin <= hepEmPars->fElectronTrackingCut) {
     eloss     = thePreStepEkin;
     finalEkin = 0.0;
     theTrack->SetEKin(finalEkin);
     theTrack->SetEnergyDeposit(eloss);
     return true;
   }
   theTrack->SetEKin(finalEkin);
   theTrack->SetEnergyDeposit(eloss);
   return false;
 }
 
 // Here I can have my own transportation to be called BUT at the moment I cannot
 // skip the G4Transportation if I do it by myself !!!
 
 // Note: energy deposit will be set here i.e. this is the first access to it that
 //       will clear the previous step value.
 bool G4HepEmElectronManager::PerformContinuous(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmElectronTrack* theElTrack, G4HepEmRandomEngine* rnge) {
   theElTrack->SavePreStepEKin();
   //
   // === 1. MSC should be invoked to obtain the physics step Length
   UpdatePStepLength(theElTrack);
   const double pStepLength = theElTrack->GetPStepLength();
 
   if (pStepLength<=0.0) {
     return false;
   }
   // compute the energy loss first based on the new step length: it will be needed in the
   // MSC scatteirng and displacement computation here as well (that is done only if not
   // the last step with the particle).
   // But update the number of interaction length left before.
   //
   // === 2. The `number-of-interaction-left` needs to be updated based on the actual
   //        physical step Length
   UpdateNumIALeft(theElTrack);
   //
   // === 3. Continuous energy loss needs to be computed
   bool stopped = ApplyMeanEnergyLoss(hepEmData, hepEmPars, theElTrack);
   if (stopped) {
     return true;
   }
 
   // === 4. Sample MSC direction change and displacement.
   SampleMSC(hepEmData, hepEmPars, theElTrack, rnge);
 
   // === 5. Sample loss fluctuations.
   return SampleLossFluctuations(hepEmData, hepEmPars, theElTrack, rnge);
 }
 
 
 bool G4HepEmElectronManager::CheckDelta(struct G4HepEmData* hepEmData, G4HepEmTrack* theTrack, double rand) {
   const bool isElectron = (theTrack->GetCharge() < 0.0);
   const G4HepEmElectronData* elData = isElectron
                                       ? hepEmData->fTheElectronData
                                       : hepEmData->fThePositronData;
   const int iDProc      = theTrack->GetWinnerProcessIndex();
   const int theIMC      = theTrack->GetMCIndex();
   const int theMatIndex = hepEmData->fTheMatCutData->fMatCutData[theIMC].fHepEmMatIndex;
   const double theEkin  = theTrack->GetEKin();
   const double theLEkin = theTrack->GetLogEKin();
   const double mxsec = (iDProc<2 ? GetRestMacXSec(elData, theIMC, theEkin, theLEkin, iDProc==0) :
                        (iDProc<3 ? ComputeMacXsecAnnihilation(theEkin, hepEmData->fTheMaterialData->fMaterialData[theMatIndex].fElectronDensity)
                                  : GetMacXSecNuclear(elData, theMatIndex, theEkin, theLEkin))
                        );
   return mxsec <= 0.0 || rand > mxsec*theTrack->GetMFP(iDProc);
 }
 
 void G4HepEmElectronManager::PerformDiscrete(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmTLData* tlData) {
   G4HepEmElectronTrack* theElTrack = tlData->GetPrimaryElectronTrack();
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   const bool isElectron = (theTrack->GetCharge() < 0.0);
 
   // 1. check if discrete process limited the step return otherwise (i.e. if
   //      continous or boundary process limited the step)
   const int iDProc = theTrack->GetWinnerProcessIndex();
   if (iDProc < 0 || theTrack->GetOnBoundary()) {
     return;
   }
   // reset number of interaction left for the winner discrete process
   theTrack->SetNumIALeft(-1.0, iDProc);
 
   // 2. check if delta interaction happens instead of the real discrete process
   if (CheckDelta(hepEmData, theTrack, tlData->GetRNGEngine()->flat())) {
     return;
   }
 
   // 3. perform the discrete part of the winner interaction
   const double theEkin = theTrack->GetEKin();
   switch (iDProc) {
     case 0: // invoke ioni (for e-/e+):
             G4HepEmElectronInteractionIoni::Perform(tlData, hepEmData, isElectron);
             break;
     case 1: // invoke brem (for e-/e+): either SB- or Rel-Brem
             G4HepEmElectronInteractionBrem::Perform(tlData, hepEmData, isElectron, theEkin < hepEmPars->fElectronBremModelLim);
             break;
     case 2: // invoke annihilation (in-flight) for e+
             G4HepEmPositronInteractionAnnihilation::Perform(tlData, false);
             break;
     case 3: // electron/positorn - nuclear interaction is not handled by HepEm: do nothing
             break;
   }
 }
 
 void G4HepEmElectronManager::Perform(struct G4HepEmData* hepEmData, struct G4HepEmParameters* hepEmPars, G4HepEmTLData* tlData) {
   G4HepEmElectronTrack* theElTrack = tlData->GetPrimaryElectronTrack();
   G4HepEmTrack*   theTrack = theElTrack->GetTrack();
   // Set default values to cover all early returns due to protection against
   // zero step lengths, conversion errors, etc.
   theTrack->SetEnergyDeposit(0);
   theElTrack->SetPStepLength(theTrack->GetGStepLength());
   const bool isElectron = (theTrack->GetCharge() < 0.0);
 
   if (theTrack->GetGStepLength()<=0.) return;
 
   bool stopped = PerformContinuous(hepEmData, hepEmPars, theElTrack, tlData->GetRNGEngine());
   if (stopped) {
     // call annihilation for e+ !!!
     if (!isElectron) {
       G4HepEmPositronInteractionAnnihilation::Perform(tlData, true);
     }
     return;
   }
 
   // === 4. Discrete part of the interaction (if any)
   PerformDiscrete(hepEmData, hepEmPars, tlData);
 }
 
 
 double  G4HepEmElectronManager::GetRestRange(const struct G4HepEmElectronData* elData, const int imc, const double ekin, const double lekin) {
   const int numELossData = elData->fELossEnergyGridSize;
   const int iRangeStarts = 5*numELossData*imc;
   // use the G4HepEmRunUtils function for interpolation
   const double     range = GetSplineLog(numELossData, elData->fELossEnergyGrid, &(elData->fELossData[iRangeStarts]), ekin, lekin, elData->fELossLogMinEkin, elData->fELossEILDelta);
   return G4HepEmMax(0.0, range);
 }
 
 
 double  G4HepEmElectronManager::GetRestDEDX(const struct G4HepEmElectronData* elData, const int imc, const double ekin, const double lekin) {
   const int numELossData = elData->fELossEnergyGridSize;
   const int  iDEDXStarts = numELossData*(5*imc + 2); // 5*imc*numELossData is where range-start + 2*numELossData
   // use the G4HepEmRunUtils function for interpolation
   const double      dedx = GetSplineLog(numELossData, elData->fELossEnergyGrid, &(elData->fELossData[iDEDXStarts]), ekin, lekin, elData->fELossLogMinEkin, elData->fELossEILDelta);
   return G4HepEmMax(0.0, dedx);
 }
 
 
 double  G4HepEmElectronManager::GetInvRange(const struct G4HepEmElectronData* elData, int imc, double range) {
   const int numELossData = elData->fELossEnergyGridSize;
   const int iRangeStarts = 5*numELossData*imc;
   // low-energy approximation
   const double minRange = elData->fELossData[iRangeStarts];
   if (range<minRange) {
     const double dum = range/minRange;
     return G4HepEmMax(0.0, elData->fELossEnergyGrid[0]*dum*dum);
   }
   // use the G4HepEmRunUtils function for finding the range bin index and for interpolation
   // find `i`, lower index of the range such that R_{i} <= r < R_{i+1}
   const int     iRlow = FindLowerBinIndex(&(elData->fELossData[iRangeStarts]), numELossData, range, 2);
   // interpolate: x,y and sd
   const double energy = GetSpline(&(elData->fELossData[iRangeStarts]), elData->fELossEnergyGrid, &(elData->fELossData[iRangeStarts+4*numELossData]), range, iRlow, 2);
   return G4HepEmMax(0.0, energy);
 }
 
 
 double  G4HepEmElectronManager::GetRestMacXSec(const struct G4HepEmElectronData* elData, const int imc, const double ekin, const double lekin, bool isioni) {
   const int iIoniStarts = elData->fResMacXSecStartIndexPerMatCut[imc];
   const int numIoniData = elData->fResMacXSecData[iIoniStarts]; // x3 for the 3 values and +5 at the beginning
   const int      iStart = (isioni) ? iIoniStarts : iIoniStarts + 3*numIoniData + 5;
   const int     numData = elData->fResMacXSecData[iStart];
   const double  minEKin = elData->fResMacXSecData[iStart+5];
   if (ekin<minEKin) {return 0.0; }
   // use the G4HepEmRunUtils function for interpolation
   const double    mxsec = GetSplineLog(numData, &(elData->fResMacXSecData[iStart+5]), ekin, lekin, elData->fResMacXSecData[iStart+3],elData->fResMacXSecData[iStart+4]);
   return G4HepEmMax(0.0, mxsec);
 }
 
 double  G4HepEmElectronManager::GetMacXSecNuclear(const struct G4HepEmElectronData* elData, const int imat, const double ekin, const double lekin) {
   if (ekin < elData->fENucEnergyGrid[0]) { return 0.0; }
   const int numEKin   = elData->fENucEnergyGridSize; // #energy point per material
   const int iStartMat = imat*2*numEKin; // value and second derivative at each energy point
   // use the G4HepEmRunUtils function for interpolation
   const double mxsec = GetSplineLog(numEKin, elData->fENucEnergyGrid, &(elData->fENucMacXsecData[iStartMat]), ekin, lekin, elData->fENucLogMinEkin, elData->fENucEILDelta);
   return G4HepEmMax(0.0, mxsec);
 }
 
 
 double  G4HepEmElectronManager::GetRestMacXSecForStepping(const struct G4HepEmElectronData* elData, const int imc, double ekin, double lekin, bool isioni) {
   constexpr double log08 = -0.22314355131420971;
   const int  iIoniStarts = elData->fResMacXSecStartIndexPerMatCut[imc];
   const int  numIoniData = elData->fResMacXSecData[iIoniStarts]; // x3 for the 3 values and +5 at the beginning
   const int       iStart = (isioni) ? iIoniStarts : iIoniStarts + 3*numIoniData + 5;
   const int      numData = elData->fResMacXSecData[iStart];
   const double mxsecMinE = elData->fResMacXSecData[iStart+5];
   const double mxsecMaxE = elData->fResMacXSecData[iStart+1];
   const double mxsecMaxV = elData->fResMacXSecData[iStart+2];
   if (ekin > mxsecMaxE) {
     // compute reduced energy: we assume that 1/lambda is higher at lower energy so we provide an overestimate
     const double ekinReduced = 0.8 * ekin;
     if (ekinReduced < mxsecMaxE) {
       return G4HepEmMax(0.0, mxsecMaxV);
     } else {
       // otherwise we are still on the right side of the maximum so provide 1/lambda at this reduced energy
       ekin   = ekinReduced;
       lekin += log08;
     }
   }
   if (ekin<mxsecMinE) {return 0.0; }
   // use the G4HepEmRunUtils function for interpolation
   const double mxsec = GetSplineLog(numData, &(elData->fResMacXSecData[iStart+5]), ekin, lekin, elData->fResMacXSecData[iStart+3], elData->fResMacXSecData[iStart+4]);
   return G4HepEmMax(0.0, mxsec);
 }
 
 double  G4HepEmElectronManager::GetMacXSecNuclearForStepping(const struct G4HepEmElectronData* elData, const int imat, const double ekin, const double lekin) {
   // assuming increasing macroscopic cross section with increasing energy
   return GetMacXSecNuclear(elData, imat, ekin, lekin);
 }
 
 
 double  G4HepEmElectronManager::GetTransportMFP(const struct G4HepEmElectronData* elData, const int im, const double ekin, const double lekin) {
   const int numEkin = elData->fELossEnergyGridSize;
   const int iStarts = 2*numEkin*im;
   // use the G4HepEmRunUtils function for interpolation
   const double tr1mxsec = G4HepEmMax(0.0, GetSplineLog(numEkin, elData->fELossEnergyGrid, &(elData->fTr1MacXSecData[iStarts]), ekin, lekin, elData->fELossLogMinEkin, elData->fELossEILDelta));
   return tr1mxsec > 0. ? 1./tr1mxsec : kALargeValue;
 }
 
 
 double G4HepEmElectronManager::ComputeMacXsecAnnihilation(const double ekin, const double electronDensity) {
   // Heitler model for e+e- -> 2g annihilation
   const double tau   = ekin*kInvElectronMassC2;
   const double gam   = tau + 1.0;
   const double gam2  = gam*gam;
   const double bg2   = tau * (tau+2.0);
   const double bg    = std::sqrt(bg2);
   return electronDensity*kPir02*((gam2+4.*gam+1.)*G4HepEmLog(gam+bg) - (gam+3.)*bg) / (bg2*(gam+1.));
 }
 
 
 double G4HepEmElectronManager::ComputeMacXsecAnnihilationForStepping(const double ekin, const double electronDensity) {
   // compute mxsec for the reduced energy (assuming that the mac-xsec decreasing with ekin)
   return ComputeMacXsecAnnihilation(0.8*ekin, electronDensity);
 }
 
 
 void G4HepEmElectronManager::ConvertTrueToGeometricLength(const G4HepEmData* hepEmData, G4HepEmMSCTrackData* mscData,
                                                          double ekin, double range, int imc, bool iselectron) {
   mscData->fPar1 = -1.;
   mscData->fPar2 =  0.;
   mscData->fPar3 =  0.;
   // NOTE: in theory, the continuous e-loss limits the step length such its at most the range.
   //       So this line below should never have any effects! Since this is the only part that
   //       might limits the true step length in the true --> geom conversion, we might drop this
   //       that makes possible the decision if MSC limited the step even before this conversion.
   mscData->fTrueStepLength = G4HepEmMin(mscData->fTrueStepLength, range);
   //  do the true -> geom transformation
   mscData->fZPathLength = mscData->fTrueStepLength;
   // z = t for very small true-path-length
   const double kTlimitMinfix2 = 1.0E-6; // 1 [nm]
   if (mscData->fTrueStepLength < kTlimitMinfix2) {
     return;
   }
   //
   const double kTauSmall  = 1.0e-16;
   const double kDtrl      = 0.05;
   const double tau        = mscData->fTrueStepLength / mscData->fLambtr1;
   if (tau < kTauSmall) {
     mscData->fZPathLength = G4HepEmMin(mscData->fTrueStepLength, mscData->fLambtr1);
   } else if (mscData->fTrueStepLength < range * kDtrl) {
     const double kTauLim  = 1.0e-6;
     mscData->fZPathLength = (tau < kTauLim)
                             ? mscData->fTrueStepLength * (1. - 0.5 * tau)
                             : mscData->fLambtr1 * (1. - G4HepEmExp(-tau));
   } else if (ekin < kElectronMassC2 || mscData->fTrueStepLength == range) {
     mscData->fPar1        = 1. / range;
     mscData->fPar2        = 1. / (mscData->fPar1 * mscData->fLambtr1);
     mscData->fPar3        = 1. + mscData->fPar2;
     mscData->fZPathLength = 1. / (mscData->fPar1 * mscData->fPar3);
     if (mscData->fTrueStepLength < range) {
       mscData->fZPathLength *= (1. - G4HepEmPow(1. - mscData->fTrueStepLength/range, mscData->fPar3));
     }
   } else {
     const double rfin     = G4HepEmMax(range - mscData->fTrueStepLength, 0.01 * range);
     const G4HepEmElectronData* elData = iselectron ? hepEmData->fTheElectronData : hepEmData->fThePositronData;
     const double t1       = GetInvRange(elData, imc, rfin);
     const int    imat     = (hepEmData->fTheMatCutData->fMatCutData[imc]).fHepEmMatIndex;
     const double lambda1  = GetTransportMFP(elData, imat, t1, G4HepEmLog(t1));
     mscData->fPar1        = (mscData->fLambtr1 - lambda1) / (mscData->fLambtr1 * mscData->fTrueStepLength); // alpha
     mscData->fPar2        = 1. / (mscData->fPar1 * mscData->fLambtr1);
     mscData->fPar3        = 1. + mscData->fPar2;
     mscData->fZPathLength = (1. - G4HepEmPow(lambda1/mscData->fLambtr1, mscData->fPar3))/(mscData->fPar1 * mscData->fPar3);
   }
   mscData->fZPathLength = G4HepEmMin(mscData->fZPathLength, mscData->fLambtr1);
 }
 
 
 void G4HepEmElectronManager::ConvertGeometricToTrueLength(G4HepEmMSCTrackData* mscData, double range, double gStepToConvert) {
   // step was not defined by transportation: i.e. physics so we know everything,
   // i.e. fTrueStepLength is known because the particle went as far as we expected
   // NOTE: this is done directly now in the caller
 //  if (!onBoundary) {
 //    return;
 //  }
 //   return;
   // else: convert geom -> true by using the mean value
   // get the geometrical step length
   mscData->fZPathLength = gStepToConvert;
   // t = z for very small step
   const double kTLimitMinfix2 = 1.0E-6; // 1 [nm]
   if (gStepToConvert < kTLimitMinfix2) {
     mscData->fTrueStepLength = gStepToConvert;
   } else {
     // recalculation
     const double kTauSmall = 1.0e-16;
     double tlength = gStepToConvert;
     if (gStepToConvert > mscData->fLambtr1 * kTauSmall) {
       if (mscData->fPar1 < 0.) {
         tlength = -mscData->fLambtr1 * G4HepEmLog(1. - gStepToConvert / mscData->fLambtr1);
       } else {
         const double dum = mscData->fPar1 * mscData->fPar3 * gStepToConvert;
         if (dum < 1.) {
           tlength = (1. - G4HepEmPow(1.-dum, 1./mscData->fPar3)) / mscData->fPar1;
         } else {
           tlength = range;
         }
       }
       if (tlength < gStepToConvert) {
         tlength = gStepToConvert;
       }
     }
     // we make sure in the caller that the results of g->t conversion is not longer
     // than the original true step length
     mscData->fTrueStepLength = tlength;
   }
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team