G4PenelopeIonisationXSHandler.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id$
//
// Author: Luciano Pandola
//
// History:
// --------
// 08 Mar 2012   L Pandola    First complete implementation
//

#include "G4PenelopeIonisationXSHandler.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4ParticleDefinition.hh"
#include "G4Electron.hh"
#include "G4Positron.hh"
#include "G4PenelopeOscillatorManager.hh"
#include "G4PenelopeOscillator.hh"
#include "G4PenelopeCrossSection.hh"
#include "G4PhysicsFreeVector.hh"
#include "G4PhysicsLogVector.hh" 

G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler(size_t nb)
  :XSTableElectron(0),XSTablePositron(0),
   theDeltaTable(0),energyGrid(0)
{
  nBins = nb;
  G4double LowEnergyLimit = 100.0*eV;
  G4double HighEnergyLimit = 100.0*GeV;
  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
  XSTableElectron = new 
    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
  XSTablePositron = new 
    std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;

  theDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  energyGrid = new G4PhysicsLogVector(LowEnergyLimit,
                                      HighEnergyLimit, 
                                      nBins-1); //one hidden bin is added

  verboseLevel = 0;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4PenelopeIonisationXSHandler::~G4PenelopeIonisationXSHandler()
{
  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i;
  if (XSTableElectron)
    {
      for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++)
        {
          G4PenelopeCrossSection* tab = i->second;
          delete tab;
        }
      delete XSTableElectron;
      XSTableElectron = 0;
    }

  if (XSTablePositron)
    {
      for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++)
        {
          G4PenelopeCrossSection* tab = i->second;
          delete tab;
        }
      delete XSTablePositron;
      XSTablePositron = 0;
    }

  std::map<const G4Material*,G4PhysicsFreeVector*>::iterator k;
  if (theDeltaTable)
    {
      for (k=theDeltaTable->begin();k!=theDeltaTable->end();k++)        
        delete k->second;
      delete theDeltaTable;
      theDeltaTable = 0;
    }
    
  if (energyGrid)
    delete energyGrid;

  if (verboseLevel > 2)
    G4cout << "G4PenelopeIonisationXSHandler. Tables have been cleared" 
           << G4endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4PenelopeCrossSection* 
G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
                                                                const G4Material* mat,
                                                                G4double cut)
{
  if (part != G4Electron::Electron() && part != G4Positron::Positron())
    {
      G4ExceptionDescription ed;
      ed << "Invalid particle: " << part->GetParticleName() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
                  "em0001",FatalException,ed);
      return NULL;
    }

  if (part == G4Electron::Electron())
    {
      if (!XSTableElectron)
        {
          G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
                      "em0028",FatalException,  
                      "The Cross Section Table for e- was not initialized correctly!");
          return NULL;
        }
      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
      if (XSTableElectron->count(theKey)) //table already built 
        return XSTableElectron->find(theKey)->second;
      else
        {
          BuildXSTable(mat,cut,part);
          if (XSTableElectron->count(theKey)) //now it should be ok!
            return XSTableElectron->find(theKey)->second;
          else
            {
              G4ExceptionDescription ed;
              ed << "Unable to build e- table for " << mat->GetName() << G4endl;
              G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
                          "em0029",FatalException,ed);     
            }
        }
    }

  if (part == G4Positron::Positron())
    {
      if (!XSTablePositron)
        {
          G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
                      "em0028",FatalException,  
                      "The Cross Section Table for e+ was not initialized correctly!");
          return NULL;
        }
      std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
      if (XSTablePositron->count(theKey)) //table already built 
        return XSTablePositron->find(theKey)->second;
      else
        {
          BuildXSTable(mat,cut,part);
          if (XSTablePositron->count(theKey)) //now it should be ok!
            return XSTablePositron->find(theKey)->second;
          else
            {
              G4ExceptionDescription ed;
              ed << "Unable to build e+ table for " << mat->GetName() << G4endl;
              G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
                          "em0029",FatalException,ed);             
            }
        }
    }
  return NULL;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void G4PenelopeIonisationXSHandler::BuildXSTable(const G4Material* mat,G4double cut,
                                             const G4ParticleDefinition* part)
{
  //
  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
  //and for the given material/cut couple. The calculation is done as sum over the 
  //individual shells.
  //Equivalent of subroutines EINaT and PINaT of Penelope
  //
  if (verboseLevel > 2)
    {
      G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl;
      G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
    }

  //Tables have been already created (checked by GetCrossSectionTableForCouple)
  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
  size_t numberOfOscillators = theTable->size();

  if (energyGrid->GetVectorLength() != nBins) 
    {
      G4ExceptionDescription ed;
      ed << "Energy Grid looks not initialized" << G4endl;
      ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
                  "em2030",FatalException,ed);
    }

  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators);
 
  //loop on the energy grid
  for (size_t bin=0;bin<nBins;bin++)
    {
       G4double energy = energyGrid->GetLowEdgeEnergy(bin);
       G4double XH0=0, XH1=0, XH2=0;
       G4double XS0=0, XS1=0, XS2=0;
   
       //oscillator loop
       for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
         {
           G4DataVector* tempStorage = 0;

           G4PenelopeOscillator* theOsc = (*theTable)[iosc];
           G4double delta = GetDensityCorrection(mat,energy);
           if (part == G4Electron::Electron())       
             tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
           else if (part == G4Positron::Positron())
             tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
           //check results are all right
           if (!tempStorage)
             {
               G4ExceptionDescription ed;
               ed << "Problem in calculating the shell XS for shell # " 
                  << iosc << G4endl;
               G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
                           "em2031",FatalException,ed);                    
               delete XSEntry;
               return;
             }
           if (tempStorage->size() != 6)
             {
               G4ExceptionDescription ed;
               ed << "Problem in calculating the shell XS " << G4endl;
               ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
               G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
                           "em2031",FatalException,ed); 
             }
           G4double stre = theOsc->GetOscillatorStrength();

           XH0 += stre*(*tempStorage)[0];
           XH1 += stre*(*tempStorage)[1];
           XH2 += stre*(*tempStorage)[2];
           XS0 += stre*(*tempStorage)[3];
           XS1 += stre*(*tempStorage)[4];
           XS2 += stre*(*tempStorage)[5];
           XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]);
           if (tempStorage)
             {
               delete tempStorage;
               tempStorage = 0;
             }
         }       
       XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
    }

  //Insert in the appropriate table
  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
  if (part == G4Electron::Electron())      
    XSTableElectron->insert(std::make_pair(theKey,XSEntry));
  else if (part == G4Positron::Positron())
    XSTablePositron->insert(std::make_pair(theKey,XSEntry));
  else
    delete XSEntry;
  
  return;
}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4PenelopeIonisationXSHandler::GetDensityCorrection(const G4Material* mat,
                                                         G4double energy)
{
  G4double result = 0;
  if (!theDeltaTable)
    {
      G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
                  "em2032",FatalException,
                  "Delta Table not initialized. Was Initialise() run?");
      return 0;
    }
  if (energy <= 0*eV)
    {
      G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl;
      G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
      return 0;
    }
  G4double logene = std::log(energy);

  //check if the material has been built
  if (!(theDeltaTable->count(mat)))
    BuildDeltaTable(mat);

  if (theDeltaTable->count(mat))
    {
      G4PhysicsFreeVector* vec = theDeltaTable->find(mat)->second;
      result = vec->Value(logene); //the table has delta vs. ln(E)      
    }
  else
    {
      G4ExceptionDescription ed;      
      ed << "Unable to build table for " << mat->GetName() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
                  "em2033",FatalException,ed);
    }

  return result;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void G4PenelopeIonisationXSHandler::BuildDeltaTable(const G4Material* mat)
{
  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
  G4double plasmaSq = oscManager->GetPlasmaEnergySquared(mat);
  G4double totalZ = oscManager->GetTotalZ(mat);
  size_t numberOfOscillators = theTable->size();

  if (energyGrid->GetVectorLength() != nBins) 
    {
      G4ExceptionDescription ed;
      ed << "Energy Grid for Delta table looks not initialized" << G4endl;
      ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
      G4Exception("G4PenelopeIonisationXSHandler::BuildDeltaTable()",
                  "em2030",FatalException,ed);
    }

  G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(nBins);

  //loop on the energy grid
  for (size_t bin=0;bin<nBins;bin++)
    {
      G4double delta = 0.;
      G4double energy = energyGrid->GetLowEdgeEnergy(bin);

      //Here calculate delta
      G4double gam = 1.0+(energy/electron_mass_c2);
      G4double gamSq = gam*gam;

      G4double TST = totalZ/(gamSq*plasmaSq);
      G4double wl2 = 0;
      G4double fdel = 0;

      //loop on oscillators
      for (size_t i=0;i<numberOfOscillators;i++)
        {
          G4PenelopeOscillator* theOsc = (*theTable)[i];
          G4double wri = theOsc->GetResonanceEnergy();
          fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
        }      
      if (fdel >= TST) //if fdel < TST, delta = 0
        {
          //get last oscillator
          G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1]; 
          wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy();
          
          //First iteration
          G4bool loopAgain = false;
          do
            {
              loopAgain = false;
              wl2 += wl2;
              fdel = 0.;
              for (size_t i=0;i<numberOfOscillators;i++)
                {
                  G4PenelopeOscillator* theOscLocal1 = (*theTable)[i];
                  G4double wri = theOscLocal1->GetResonanceEnergy();
                  fdel += theOscLocal1->GetOscillatorStrength()/(wri*wri+wl2);
                }
              if (fdel > TST)
                loopAgain = true;
            }while(loopAgain);

          G4double wl2l = 0;
          G4double wl2u = wl2;
          //second iteration
          do
            {        
              loopAgain = false;
              wl2 = 0.5*(wl2l+wl2u);
              fdel = 0;
              for (size_t i=0;i<numberOfOscillators;i++)
                {
                  G4PenelopeOscillator* theOscLocal2 = (*theTable)[i];
                  G4double wri = theOscLocal2->GetResonanceEnergy();
                  fdel += theOscLocal2->GetOscillatorStrength()/(wri*wri+wl2);
                }
              if (fdel > TST)
                wl2l = wl2;
              else
                wl2u = wl2;
              if ((wl2u-wl2l)>1e-12*wl2)
                loopAgain = true;
            }while(loopAgain);
          
          //Eventually get density correction
          delta = 0.;
          for (size_t i=0;i<numberOfOscillators;i++)
            {
              G4PenelopeOscillator* theOscLocal3 = (*theTable)[i];
              G4double wri = theOscLocal3->GetResonanceEnergy();
              delta += theOscLocal3->GetOscillatorStrength()*
                std::log(1.0+(wl2/(wri*wri)));               
            }
          delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
        }
      energy = std::max(1e-9*eV,energy); //prevents log(0)
      theVector->PutValue(bin,std::log(energy),delta);
    }
  theDeltaTable->insert(std::make_pair(mat,theVector));
  return;
}
                                                        
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc,
                                                                                  G4double energy,
                                                                                  G4double cut,
                                                                                  G4double delta)
{
  //
  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 
  //the given oscillator/cut and at the given energy.
  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
  //Equivalent of subroutines EINaT1 of Penelope
  //
  // Results are _per target electron_
  //
  G4DataVector* result = new G4DataVector();
  for (size_t i=0;i<6;i++)
    result->push_back(0.);
  G4double ionEnergy = theOsc->GetIonisationEnergy();
  
  //return a set of zero's if the energy it too low to excite the current oscillator
  if (energy < ionEnergy)
    return result;

  G4double H0=0.,H1=0.,H2=0.;
  G4double S0=0.,S1=0.,S2=0.;

  //Define useful constants to be used in the calculation
  G4double gamma = 1.0+energy/electron_mass_c2;
  G4double gammaSq = gamma*gamma;
  G4double beta = (gammaSq-1.0)/gammaSq;
  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
  G4double constant = pielr2*2.0*electron_mass_c2/beta;
  G4double XHDT0 = std::log(gammaSq)-beta;

  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
  G4double cp = std::sqrt(cpSq);
  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));

  //
  // Distant interactions
  //
  G4double resEne = theOsc->GetResonanceEnergy();
  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
  if (energy > resEne)
    {
      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
      G4double cp1 = std::sqrt(cp1Sq);
      
      //Distant longitudinal interactions
      G4double QM = 0;
      if (resEne > 1e-6*energy)
        QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
      else
        {
          QM = resEne*resEne/(beta*2.0*electron_mass_c2);
          QM = QM*(1.0-0.5*QM/electron_mass_c2);
        }
      G4double SDL1 = 0;
      if (QM < cutoffEne)
        SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
      
      //Distant transverse interactions
      if (SDL1)
        {
          G4double SDT1 = std::max(XHDT0-delta,0.0);
          G4double SD1 = SDL1+SDT1;
          if (cut > resEne)
            {
              S1 = SD1; //XS1
              S0 = SD1/resEne; //XS0
              S2 = SD1*resEne; //XS2
            }
          else
            {
              H1 = SD1; //XH1
              H0 = SD1/resEne; //XH0
              H2 = SD1*resEne; //XH2
            }
        }
    }
  //
  // Close collisions (Moller's cross section)
  //
  G4double wl = std::max(cut,cutoffEne);
  G4double ee = energy + ionEnergy;
  G4double wu = 0.5*ee;
  if (wl < wu-(1e-5*eV))
    {
      H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
        (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + 
        amol*(wu-wl)/(ee*ee);
      H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
        (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 
        amol*(wu*wu-wl*wl)/(2.0*ee*ee);
      H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
        (wl*(2.0*ee-wl)/(ee-wl)) + 
        (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +       
        amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
      wu = wl;
    }
  wl = cutoffEne;
  
  if (wl > wu-(1e-5*eV))
    {
      (*result)[0] = constant*H0;
      (*result)[1] = constant*H1;
      (*result)[2] = constant*H2;
      (*result)[3] = constant*S0;
      (*result)[4] = constant*S1;
      (*result)[5] = constant*S2;
      return result;
    }

  S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 
    (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
    amol*(wu-wl)/(ee*ee);
  S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 
    (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 
    amol*(wu*wu-wl*wl)/(2.0*ee*ee);
  S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 
    (wl*(2.0*ee-wl)/(ee-wl)) + 
    (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + 
    amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);

  (*result)[0] = constant*H0;
  (*result)[1] = constant*H1;
  (*result)[2] = constant*H2;
  (*result)[3] = constant*S0;
  (*result)[4] = constant*S1;
  (*result)[5] = constant*S2;
  return result;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc,
                                                                                  G4double energy,
                                                                                  G4double cut,
                                                                                  G4double delta)
{
  //
  //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 
  //the given oscillator/cut and at the given energy.
  //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
  //Equivalent of subroutines PINaT1 of Penelope
  //
  // Results are _per target electron_
  //
  G4DataVector* result = new G4DataVector();
  for (size_t i=0;i<6;i++)
    result->push_back(0.);
  G4double ionEnergy = theOsc->GetIonisationEnergy();
  
  //return a set of zero's if the energy it too low to excite the current oscillator
  if (energy < ionEnergy)
    return result;

  G4double H0=0.,H1=0.,H2=0.;
  G4double S0=0.,S1=0.,S2=0.;

  //Define useful constants to be used in the calculation
  G4double gamma = 1.0+energy/electron_mass_c2;
  G4double gammaSq = gamma*gamma;
  G4double beta = (gammaSq-1.0)/gammaSq;
  G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
  G4double constant = pielr2*2.0*electron_mass_c2/beta;
  G4double XHDT0 = std::log(gammaSq)-beta;

  G4double cpSq = energy*(energy+2.0*electron_mass_c2);
  G4double cp = std::sqrt(cpSq);
  G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
  G4double g12 = (gamma+1.0)*(gamma+1.0);
  //Bhabha coefficients
  G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
  G4double bha2 = amol*(3.0+1.0/g12);
  G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
  G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;

  //
  // Distant interactions
  //
  G4double resEne = theOsc->GetResonanceEnergy();
  G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
  if (energy > resEne)
    {
      G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
      G4double cp1 = std::sqrt(cp1Sq);
      
      //Distant longitudinal interactions
      G4double QM = 0;
      if (resEne > 1e-6*energy)
        QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
      else
        {
          QM = resEne*resEne/(beta*2.0*electron_mass_c2);
          QM = QM*(1.0-0.5*QM/electron_mass_c2);
        }
      G4double SDL1 = 0;
      if (QM < cutoffEne)
        SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
      
      //Distant transverse interactions
      if (SDL1)
        {
          G4double SDT1 = std::max(XHDT0-delta,0.0);
          G4double SD1 = SDL1+SDT1;
          if (cut > resEne)
            {
              S1 = SD1; //XS1
              S0 = SD1/resEne; //XS0
              S2 = SD1*resEne; //XS2
            }
          else
            {
              H1 = SD1; //XH1
              H0 = SD1/resEne; //XH0
              H2 = SD1*resEne; //XH2
            }
        }
    }

  //
  // Close collisions (Bhabha's cross section)
  //
  G4double wl = std::max(cut,cutoffEne);
  G4double wu = energy; 
  G4double energySq = energy*energy;
  if (wl < wu-(1e-5*eV))
    {
      G4double wlSq = wl*wl;
      G4double wuSq = wu*wu;
      H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy  
        + bha2*(wu-wl)/energySq  
        - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
        + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
      H1 += std::log(wu/wl) - bha1*(wu-wl)/energy
        + bha2*(wuSq-wlSq)/(2.0*energySq)
        - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
        + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
      H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
        + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
        - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
        + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
      wu = wl;
    }
  wl = cutoffEne;
  
  if (wl > wu-(1e-5*eV))
    {
      (*result)[0] = constant*H0;
      (*result)[1] = constant*H1;
      (*result)[2] = constant*H2;
      (*result)[3] = constant*S0;
      (*result)[4] = constant*S1;
      (*result)[5] = constant*S2;
      return result;
    }

  G4double wlSq = wl*wl;
  G4double wuSq = wu*wu;

  S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy 
    + bha2*(wu-wl)/energySq  
    - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
    + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);

  S1 += std::log(wu/wl) - bha1*(wu-wl)/energy
    + bha2*(wuSq-wlSq)/(2.0*energySq)
    - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
    + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);

  S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
    + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
    - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
    + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);

 (*result)[0] = constant*H0;
 (*result)[1] = constant*H1;
 (*result)[2] = constant*H2;
 (*result)[3] = constant*S0;
 (*result)[4] = constant*S1;
 (*result)[5] = constant*S2;

 return result;
}

---