Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4eIonisation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4eIonisation.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4eIonisation.cc (Version 2.0)


                                                   >>   1 // This code implementation is the intellectual property of
                                                   >>   2 // the GEANT4 collaboration.
  1 //                                                  3 //
  2 // ******************************************* <<   4 // By copying, distributing or modifying the Program (or any work
  3 // * License and Disclaimer                    <<   5 // based on the Program) you indicate your acceptance of this statement,
  4 // *                                           <<   6 // and all its terms.
  5 // * The  Geant4 software  is  copyright of th << 
  6 // * the Geant4 Collaboration.  It is provided << 
  7 // * conditions of the Geant4 Software License << 
  8 // * LICENSE and available at  http://cern.ch/ << 
  9 // * include a list of copyright holders.      << 
 10 // *                                           << 
 11 // * Neither the authors of this software syst << 
 12 // * institutes,nor the agencies providing fin << 
 13 // * work  make  any representation or  warran << 
 14 // * regarding  this  software system or assum << 
 15 // * use.  Please see the license in the file  << 
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                           << 
 18 // * This  code  implementation is the result  << 
 19 // * technical work of the GEANT4 collaboratio << 
 20 // * By using,  copying,  modifying or  distri << 
 21 // * any work based  on the software)  you  ag << 
 22 // * use  in  resulting  scientific  publicati << 
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // ******************************************* << 
 25 //                                                  7 //
 26 // ------------------------------------------- <<   8 // $Id: G4eIonisation.cc,v 1.10 2000/05/23 14:42:21 urban Exp $
                                                   >>   9 // GEANT4 tag $Name: geant4-02-00 $
 27 //                                                 10 //
 28 // GEANT4 Class file                           <<  11 // 
                                                   >>  12 // -------------------------------------------------------------
                                                   >>  13 //      GEANT 4 class implementation file 
 29 //                                                 14 //
                                                   >>  15 //      For information related to this code contact:
                                                   >>  16 //      CERN, IT Division, ASD group
                                                   >>  17 //      History: based on object model of
                                                   >>  18 //      2nd December 1995, G.Cosmo
                                                   >>  19 //      ---------- G4eIonisation physics process -----------
                                                   >>  20 //                by Laszlo Urban, 20 March 1997 
                                                   >>  21 // **************************************************************
                                                   >>  22 // It is the first implementation of the NEW IONISATION PROCESS.
                                                   >>  23 // It calculates the ionisation of e+/e-.
                                                   >>  24 // **************************************************************
 30 //                                                 25 //
 31 // File name:     G4eIonisation                <<  26 // 07-04-98: remove 'tracking cut' of the ionizing particle, MMa 
 32 //                                             <<  27 // 04-09-98: new methods SetBining() PrintInfo()
 33 // Author:        Laszlo Urban                 <<  28 // 07-09-98: Cleanup
 34 //                                             <<  29 // 02/02/99: correction inDoIt , L.Urban
 35 // Creation date: 20.03.1997                   <<  30 // 10/02/00  modifications , new e.m. structure, L.Urban
 36 //                                             <<  31 // --------------------------------------------------------------
 37 // Modified by Michel Maire and Vladimir Ivanc <<  32  
 38 //                                             << 
 39 // ------------------------------------------- << 
 40 //                                             << 
 41 //....oooOO0OOooo........oooOO0OOooo........oo << 
 42 //....oooOO0OOooo........oooOO0OOooo........oo << 
 43                                                << 
 44 #include "G4eIonisation.hh"                        33 #include "G4eIonisation.hh"
 45 #include "G4Electron.hh"                       <<  34 #include "G4EnergyLossTables.hh"
 46 #include "G4MollerBhabhaModel.hh"              <<  35 #include "G4ios.hh"
 47 #include "G4EmParameters.hh"                   <<  36 #include "G4UnitsTable.hh"
 48 #include "G4EmStandUtil.hh"                    <<  37 
                                                   >>  38 G4double G4eIonisation::LowerBoundLambda = 1.*keV ;
                                                   >>  39 G4double G4eIonisation::UpperBoundLambda = 100.*TeV ;
                                                   >>  40 G4int    G4eIonisation::NbinLambda = 100 ;
                                                   >>  41 
                                                   >>  42 // constructor and destructor
                                                   >>  43  
                                                   >>  44 G4eIonisation::G4eIonisation(const G4String& processName)
                                                   >>  45    : G4VeEnergyLoss(processName),
                                                   >>  46      theMeanFreePathTable(NULL)
                                                   >>  47 { }
 49                                                    48 
 50 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51                                                    50 
 52 G4eIonisation::G4eIonisation(const G4String& n <<  51 G4eIonisation::~G4eIonisation() 
 53   : G4VEnergyLossProcess(name),                << 
 54     theElectron(G4Electron::Electron()),       << 
 55     isElectron(true),                          << 
 56     isInitialised(false)                       << 
 57 {                                                  52 {
 58   SetProcessSubType(fIonisation);              <<  53      if (theMeanFreePathTable) {
 59   SetSecondaryParticle(theElectron);           <<  54         theMeanFreePathTable->clearAndDestroy();
                                                   >>  55         delete theMeanFreePathTable;
                                                   >>  56      }
                                                   >>  57 
 60 }                                                  58 }
 61                                                    59 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63                                                    61 
 64 G4eIonisation::~G4eIonisation() = default;     <<  62 void G4eIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
                                                   >>  63 //  just call BuildLossTable+BuildLambdaTable
                                                   >>  64 {
                                                   >>  65     // get bining from EnergyLoss
                                                   >>  66     LowestKineticEnergy  = GetLowerBoundEloss() ;
                                                   >>  67     HighestKineticEnergy = GetUpperBoundEloss() ;
                                                   >>  68     TotBin               = GetNbinEloss() ;
                                                   >>  69  
                                                   >>  70     BuildLossTable(aParticleType) ;
                                                   >>  71 
                                                   >>  72   if(&aParticleType==G4Electron::Electron())
                                                   >>  73   {
                                                   >>  74     RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable ;
                                                   >>  75     CounterOfElectronProcess++;
                                                   >>  76   }
                                                   >>  77   else
                                                   >>  78   {
                                                   >>  79     RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable ;
                                                   >>  80     CounterOfPositronProcess++;
                                                   >>  81   }
                                                   >>  82  
                                                   >>  83     BuildLambdaTable(aParticleType) ;
                                                   >>  84  
                                                   >>  85     BuildDEDXTable(aParticleType) ;
                                                   >>  86                                              
                                                   >>  87   if(&aParticleType==G4Electron::Electron())
                                                   >>  88     PrintInfoDefinition();  
                                                   >>  89 }
 65                                                    90 
 66 //....oooOO0OOooo........oooOO0OOooo........oo     91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 67                                                    92 
 68 G4double G4eIonisation::MinPrimaryEnergy(const <<  93 void G4eIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
 69            const G4Material*,                  << 
 70            G4double cut)                       << 
 71 {                                                  94 {
 72   G4double x = cut;                            <<  95 // Build tables for the ionization energy loss
 73   if(isElectron) x += cut;                     <<  96 //  the tables are built for *MATERIALS*
 74   return x;                                    <<  97 
                                                   >>  98     const G4double twoln10 = 2.*log(10.);
                                                   >>  99     const G4double Factor = twopi_mc2_rcl2;
                                                   >> 100 
                                                   >> 101     static const G4double Tl = 0.2*keV ;
                                                   >> 102 
                                                   >> 103     G4double LowEdgeEnergy, ionloss;
                                                   >> 104     
                                                   >> 105     // material properties
                                                   >> 106     G4double ElectronDensity,Eexc,Eexcm2,Cden,Mden,Aden,X0den,X1den ;
                                                   >> 107     // some local variables
                                                   >> 108     G4double tau,Tmax,gamma,gamma2,bg2,beta2,d,d2,d3,d4,delta,x,y ;
                                                   >> 109 
                                                   >> 110     ParticleMass = aParticleType.GetPDGMass();
                                                   >> 111     G4double* ParticleCutInKineticEnergy = aParticleType.GetEnergyCuts() ;
                                                   >> 112 
                                                   >> 113     //  create table
                                                   >> 114     
                                                   >> 115     const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
                                                   >> 116     G4int numOfMaterials = theMaterialTable->length();
                                                   >> 117 
                                                   >> 118      if (theLossTable) { theLossTable->clearAndDestroy();
                                                   >> 119                          delete theLossTable;
                                                   >> 120                        }
                                                   >> 121                        
                                                   >> 122     theLossTable = new G4PhysicsTable(numOfMaterials);
                                                   >> 123 
                                                   >> 124 //  loop for materials
                                                   >> 125 
                                                   >> 126     for (G4int J=0; J<numOfMaterials; J++)
                                                   >> 127     {
                                                   >> 128       // create physics vector and fill it
                                                   >> 129 
                                                   >> 130       G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
                                                   >> 131                     LowestKineticEnergy, HighestKineticEnergy, TotBin);
                                                   >> 132   
                                                   >> 133       // get material parameters needed for the energy loss calculation   
                                                   >> 134       const G4Material* material= (*theMaterialTable)[J];
                                                   >> 135 
                                                   >> 136       ElectronDensity = material->GetElectronDensity();
                                                   >> 137       Eexc   = material->GetIonisation()->GetMeanExcitationEnergy();
                                                   >> 138       Eexc  /= ParticleMass; Eexcm2 = Eexc*Eexc;
                                                   >> 139       Cden   = material->GetIonisation()->GetCdensity();
                                                   >> 140       Mden   = material->GetIonisation()->GetMdensity();
                                                   >> 141       Aden   = material->GetIonisation()->GetAdensity();
                                                   >> 142       X0den  = material->GetIonisation()->GetX0density();
                                                   >> 143       X1den  = material->GetIonisation()->GetX1density();
                                                   >> 144 
                                                   >> 145       // for the lowenergy extrapolation
                                                   >> 146       G4double Zeff = material->GetTotNbOfElectPerVolume()/
                                                   >> 147                 material->GetTotNbOfAtomsPerVolume() ;
                                                   >> 148       G4double Th = 0.25*sqrt(Zeff)*keV ;
                                                   >> 149 
                                                   >> 150       G4double Tsav ;
                                                   >> 151 
                                                   >> 152       // now comes the loop for the kinetic energy values
                                                   >> 153 
                                                   >> 154       for (G4int i = 0 ; i < TotBin ; i++)
                                                   >> 155          {
                                                   >> 156           LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ; 
                                                   >> 157           //  low energy ?     
                                                   >> 158           if(LowEdgeEnergy < Th)
                                                   >> 159           {
                                                   >> 160             Tsav = LowEdgeEnergy ;
                                                   >> 161             LowEdgeEnergy = Th ;
                                                   >> 162           }
                                                   >> 163           else
                                                   >> 164             Tsav = 0. ;
                                                   >> 165 
                                                   >> 166           tau = LowEdgeEnergy/ParticleMass ;
                                                   >> 167 
                                                   >> 168           // Seltzer-Berger formula 
                                                   >> 169           gamma = tau + 1.; gamma2 = gamma*gamma; 
                                                   >> 170           bg2 = tau*(tau+2.);
                                                   >> 171           beta2 = bg2/gamma2;
                                                   >> 172 
                                                   >> 173           // electron
                                                   >> 174           if (&aParticleType==G4Electron::Electron())
                                                   >> 175             {
                                                   >> 176               Tmax = LowEdgeEnergy/2.;  
                                                   >> 177               d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
                                                   >> 178               ionloss = log(2.*(tau+2.)/Eexcm2)-1.-beta2
                                                   >> 179                        + log((tau-d)*d)+tau/(tau-d)
                                                   >> 180                        + (0.5*d*d+(2.*tau+1.)*log(1.-d/tau))/gamma2;
                                                   >> 181             }
                                                   >> 182           else        //positron
                                                   >> 183             {
                                                   >> 184               Tmax = LowEdgeEnergy ;  
                                                   >> 185               d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
                                                   >> 186               d2=d*d/2.; d3=d*d*d/3.; d4=d*d*d*d/4.;
                                                   >> 187               y=1./(1.+gamma);
                                                   >> 188               ionloss = log(2.*(tau+2.)/Eexcm2)+log(tau*d)
                                                   >> 189                        - beta2*(tau+2.*d-y*(3.*d2+y*(d-d3+y*(d2-tau*d3+d4))))/tau;
                                                   >> 190             } 
                                                   >> 191 
                                                   >> 192           //density correction
                                                   >> 193           x = log(bg2)/twoln10;
                                                   >> 194           if (x < X0den) delta = 0.;
                                                   >> 195           else { delta = twoln10*x - Cden;
                                                   >> 196                  if (x < X1den) delta += Aden*pow((X1den-x),Mden);
                                                   >> 197                } 
                                                   >> 198 
                                                   >> 199           //now you can compute the total ionization loss
                                                   >> 200           ionloss -= delta ;
                                                   >> 201           ionloss *= Factor*ElectronDensity/beta2 ;
                                                   >> 202           if (ionloss <= 0.) ionloss = 0.;
                                                   >> 203    
                                                   >> 204           // low energy ?
                                                   >> 205           if(Tsav > 0.)
                                                   >> 206           {
                                                   >> 207             if(Tsav >= Tl)
                                                   >> 208               ionloss *= sqrt(LowEdgeEnergy/Tsav) ;
                                                   >> 209             else
                                                   >> 210               ionloss *= sqrt(LowEdgeEnergy*Tsav)/Tl ;
                                                   >> 211           }
                                                   >> 212 
                                                   >> 213           aVector->PutValue(i,ionloss) ;
                                                   >> 214          }          
                                                   >> 215       theLossTable->insert(aVector);
                                                   >> 216     }
 75 }                                                 217 }
 76                                                   218 
 77 //....oooOO0OOooo........oooOO0OOooo........oo    219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 78                                                   220 
 79 G4bool G4eIonisation::IsApplicable(const G4Par << 221 void G4eIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
 80 {                                                 222 {
 81   return (&p == theElectron || &p == G4Positro << 223   // Build mean free path tables for the delta ray production process
                                                   >> 224   //     tables are built for MATERIALS 
                                                   >> 225 
                                                   >> 226   G4double LowEdgeEnergy, Value, SIGMA;
                                                   >> 227 
                                                   >> 228   //create table
                                                   >> 229   const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
                                                   >> 230   G4int numOfMaterials = theMaterialTable->length();
                                                   >> 231 
                                                   >> 232   if (theMeanFreePathTable) { theMeanFreePathTable->clearAndDestroy();
                                                   >> 233                               delete theMeanFreePathTable;
                                                   >> 234                             }
                                                   >> 235 
                                                   >> 236   theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
                                                   >> 237 
                                                   >> 238   // get electron  cuts in kinetic energy
                                                   >> 239   // The electron cuts needed in the case of the positron , too!
                                                   >> 240   // This is the reason why SetCut has to be called for electron first !!!!!!!
                                                   >> 241 
                                                   >> 242   if((G4Electron::Electron()->GetCutsInEnergy() == 0) &&
                                                   >> 243      ( &aParticleType == G4Positron::Positron()))
                                                   >> 244   {
                                                   >> 245      G4cout << " The ELECTRON energy cuts needed to compute energy loss/mean free path "
                                                   >> 246                " for POSITRON , too. " << G4endl;
                                                   >> 247      G4Exception(" Call SetCut for e- first !!!!!!") ;
                                                   >> 248   }
                                                   >> 249    
                                                   >> 250   G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
                                                   >> 251  
                                                   >> 252   // loop for materials 
                                                   >> 253 
                                                   >> 254  for (G4int J=0 ; J < numOfMaterials; J++)
                                                   >> 255     { 
                                                   >> 256      //create physics vector then fill it ....
                                                   >> 257 
                                                   >> 258      G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
                                                   >> 259                LowerBoundLambda, UpperBoundLambda, NbinLambda);
                                                   >> 260 
                                                   >> 261      // compute the (macroscopic) cross section first
                                                   >> 262  
                                                   >> 263      const G4Material* material= (*theMaterialTable)[J];        
                                                   >> 264      const 
                                                   >> 265      G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 266      const
                                                   >> 267      G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
                                                   >> 268      const
                                                   >> 269      G4int NumberOfElements = material->GetNumberOfElements() ;
                                                   >> 270  
                                                   >> 271      // get the electron kinetic energy cut for the actual material,
                                                   >> 272      // it will be used in ComputeMicroscopicCrossSection
                                                   >> 273      // (--> it will be the same for all the elements in this material )
                                                   >> 274      G4double DeltaThreshold = DeltaCutInKineticEnergy[J] ;
                                                   >> 275 
                                                   >> 276      for (G4int i = 0 ; i < NbinLambda ; i++)
                                                   >> 277         {
                                                   >> 278           LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
                                                   >> 279           SIGMA = 0.;           
                                                   >> 280           for (G4int iel=0; iel<NumberOfElements; iel++ )
                                                   >> 281              {
                                                   >> 282                 SIGMA += theAtomicNumDensityVector[iel]*
                                                   >> 283                          ComputeMicroscopicCrossSection( aParticleType,
                                                   >> 284                                                          LowEdgeEnergy,
                                                   >> 285                                       (*theElementVector)(iel)->GetZ(),
                                                   >> 286                                                        DeltaThreshold);
                                                   >> 287              }
                                                   >> 288 
                                                   >> 289           // mean free path = 1./macroscopic cross section
                                                   >> 290           Value = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;     
                                                   >> 291           aVector->PutValue(i, Value) ;
                                                   >> 292         }
                                                   >> 293      theMeanFreePathTable->insert(aVector);
                                                   >> 294     }
 82 }                                                 295 }
 83                                                   296 
 84 //....oooOO0OOooo........oooOO0OOooo........oo    297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 85                                                   298 
 86 void G4eIonisation::InitialiseEnergyLossProces << 299 G4double G4eIonisation::ComputeMicroscopicCrossSection(
 87         const G4ParticleDefinition* part,      << 300                                  const G4ParticleDefinition& aParticleType,
 88         const G4ParticleDefinition*)           << 301                                  G4double KineticEnergy,
                                                   >> 302                                  G4double AtomicNumber ,
                                                   >> 303                                  G4double DeltaThreshold)
                                                   >> 304 {
                                                   >> 305   // calculates the microscopic cross section
                                                   >> 306   //(it is called for elements , AtomicNumber = Z )
                                                   >> 307  
                                                   >> 308   G4double MaxKineticEnergyTransfer, TotalCrossSection(0.);
                                                   >> 309   
                                                   >> 310   ParticleMass = aParticleType.GetPDGMass();
                                                   >> 311   G4double TotalEnergy = KineticEnergy + ParticleMass;
                                                   >> 312 
                                                   >> 313   G4double betasquare = KineticEnergy*(TotalEnergy+ParticleMass)
                                                   >> 314                        /(TotalEnergy*TotalEnergy);
                                                   >> 315   G4double gamma = TotalEnergy/ParticleMass, gamma2 = gamma*gamma;
                                                   >> 316   G4double x=DeltaThreshold/KineticEnergy, x2 = x*x;
                                                   >> 317 
                                                   >> 318   if (&aParticleType==G4Electron::Electron())
                                                   >> 319                             MaxKineticEnergyTransfer = 0.5*KineticEnergy;
                                                   >> 320   else                      MaxKineticEnergyTransfer =     KineticEnergy;
                                                   >> 321 
                                                   >> 322   // now you can calculate the total cross section
                                                   >> 323 
                                                   >> 324  if (MaxKineticEnergyTransfer > DeltaThreshold)
                                                   >> 325    {
                                                   >> 326     if (&aParticleType==G4Electron::Electron())   //Moller (e-e-) scattering
                                                   >> 327       {
                                                   >> 328         TotalCrossSection  = (gamma-1.)*(gamma-1.)*(0.5-x)/gamma2 + 1./x
                                                   >> 329                             - 1./(1.-x)-(2.*gamma-1.)*log((1.-x)/x)/gamma2; 
                                                   >> 330         TotalCrossSection /= betasquare; 
                                                   >> 331       }
                                                   >> 332     else                                         //Bhabha (e+e-) scattering
                                                   >> 333       {
                                                   >> 334        G4double y=1./(1.+gamma), y2 =y*y, y12=1.-2.*y;
                                                   >> 335        G4double b1=2.-y2, b2=y12*(3.+y2), b4=y12*y12*y12, b3=b4+y12*y12;
                                                   >> 336        TotalCrossSection = (1./x-1.)/betasquare+b1*log(x)+b2*(1.-x)
                                                   >> 337                           - b3*(1.-x2)/2.+b4*(1.-x2*x)/3.;
                                                   >> 338       }   
                                                   >> 339     TotalCrossSection *= (twopi_mc2_rcl2*AtomicNumber/KineticEnergy);
                                                   >> 340    }
                                                   >> 341  return TotalCrossSection ;
                                                   >> 342 }
                                                   >> 343 
                                                   >> 344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 
                                                   >> 345  
                                                   >> 346 G4VParticleChange* G4eIonisation::PostStepDoIt( const G4Track& trackData,   
                                                   >> 347                                                 const G4Step&  stepData)         
 89 {                                                 348 {
 90   if(!isInitialised) {                         << 349   aParticleChange.Initialize(trackData) ;
 91     if(part != theElectron) { isElectron = fal << 350   
 92     if (nullptr == EmModel(0)) { SetEmModel(ne << 351   G4Material*               aMaterial = trackData.GetMaterial() ;
 93     G4EmParameters* param = G4EmParameters::In << 352   const G4DynamicParticle*  aParticle = trackData.GetDynamicParticle() ;
 94     EmModel(0)->SetLowEnergyLimit(param->MinKi << 353 
 95     EmModel(0)->SetHighEnergyLimit(param->MaxK << 354   ParticleMass = aParticle->GetDefinition()->GetPDGMass();
 96     if (nullptr == FluctModel()) {             << 355   G4double KineticEnergy = aParticle->GetKineticEnergy();
 97       SetFluctModel(G4EmStandUtil::ModelOfFluc << 356   G4double TotalEnergy = KineticEnergy + ParticleMass;
                                                   >> 357   G4double Psquare = KineticEnergy*(TotalEnergy+ParticleMass);
                                                   >> 358   G4double TotalMomentum = sqrt(Psquare);
                                                   >> 359   G4double Esquare=TotalEnergy*TotalEnergy;
                                                   >> 360   G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
                                                   >> 361 
                                                   >> 362   //  get kinetic energy cut for the electron
                                                   >> 363   G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
                                                   >> 364   G4double DeltaThreshold = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
                                                   >> 365 
                                                   >> 366   // some kinematics
                                                   >> 367   G4double MaxKineticEnergyTransfer;
                                                   >> 368   if (Charge < 0.) MaxKineticEnergyTransfer = 0.5*KineticEnergy;
                                                   >> 369   else             MaxKineticEnergyTransfer =     KineticEnergy;
                                                   >> 370 
                                                   >> 371   // sampling kinetic energy of the delta ray 
                                                   >> 372 
                                                   >> 373   if (MaxKineticEnergyTransfer <= DeltaThreshold) // pathological case (should not happen,
                                                   >> 374                                                   // there is no change at all)
                                                   >> 375      return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
                                                   >> 376 
                                                   >> 377 
                                                   >> 378   // normal case
                                                   >> 379   G4double cc,y,y2,c2,b0,b1,b2,b3,b4,x,x1,grej,grejc;
                                                   >> 380  
                                                   >> 381   G4double tau = KineticEnergy/ParticleMass;
                                                   >> 382   G4double gamma = tau+1., gamma2=gamma*gamma;
                                                   >> 383   G4double xc = DeltaThreshold/KineticEnergy, xc1=1.-xc;
                                                   >> 384 
                                                   >> 385   if (Charge < 0.)  // Moller (e-e-) scattering
                                                   >> 386     { 
                                                   >> 387       b1=4./(9.*gamma2-10.*gamma+5.);
                                                   >> 388       b2=tau*tau*b1; b3=(2.*gamma2+2.*gamma-1.)*b1;
                                                   >> 389       cc=1.-2.*xc;       
                                                   >> 390       do { 
                                                   >> 391            x    = xc/(1.-cc*G4UniformRand()); x1 = 1.-x;
                                                   >> 392            grej = b2*x*x-b3*x/x1+b1*gamma2/(x1*x1);
                                                   >> 393          } while (G4UniformRand()>grej) ;
 98     }                                             394     }
 99     AddEmModel(1, EmModel(), FluctModel());    << 395   else             // Bhabha (e+e-) scattering
100     isInitialised = true;                      << 396     {
101   }                                            << 397       y=1./(gamma+1.); y2=y*y; cc=1.-2.*y;
                                                   >> 398       b1=2.-y2; b2=cc*(3.+y2);
                                                   >> 399       c2=cc*cc; b4=c2*cc; b3=c2+b4;
                                                   >> 400       b0=gamma2/(gamma2-1.);
                                                   >> 401       grejc=(((b4*xc-b3)*xc+b2)*xc-b1)*xc+b0;
                                                   >> 402       do {
                                                   >> 403            x    = xc/(1.-xc1*G4UniformRand());
                                                   >> 404            grej = ((((b4*x-b3)*x+b2)*x-b1)*x+b0)/grejc;
                                                   >> 405          } while (G4UniformRand()>grej);
                                                   >> 406     }
                                                   >> 407  
                                                   >> 408     G4double DeltaKineticEnergy = x * KineticEnergy;
                                                   >> 409 
                                                   >> 410   // protection :do not produce a secondary with 0. kinetic energy !
                                                   >> 411   if (DeltaKineticEnergy <= 0.)
                                                   >> 412       return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
                                                   >> 413 
                                                   >> 414   G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
                                                   >> 415                                                        2. * electron_mass_c2 ));
                                                   >> 416    
                                                   >> 417   G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
                                                   >> 418                       /(DeltaTotalMomentum * TotalMomentum);
                                                   >> 419 
                                                   >> 420   if (costheta < -1.) costheta = -1.;
                                                   >> 421   if (costheta > +1.) costheta = +1.;
                                                   >> 422 
                                                   >> 423   //  direction of the delta electron
                                                   >> 424 
                                                   >> 425   G4double phi = twopi * G4UniformRand(); 
                                                   >> 426   G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
                                                   >> 427   G4double dirx = sintheta * cos(phi), diry = sintheta * sin(phi), dirz = costheta;
                                                   >> 428 
                                                   >> 429   G4ThreeVector DeltaDirection(dirx,diry,dirz);
                                                   >> 430   DeltaDirection.rotateUz(ParticleDirection);
                                                   >> 431 
                                                   >> 432   // create G4DynamicParticle object for delta ray
                                                   >> 433 
                                                   >> 434   G4DynamicParticle* theDeltaRay = new G4DynamicParticle;
                                                   >> 435   theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
                                                   >> 436   theDeltaRay->SetMomentumDirection(
                                                   >> 437                    DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z()); 
                                                   >> 438   theDeltaRay->SetDefinition(G4Electron::Electron());
                                                   >> 439    
                                                   >> 440   // fill aParticleChange 
                                                   >> 441   // changed energy and momentum of the actual particle
                                                   >> 442   G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy;
                                                   >> 443   
                                                   >> 444   G4double Edep = 0. ;
                                                   >> 445 
                                                   >> 446   if (finalKineticEnergy > MinKineticEnergy)
                                                   >> 447     {
                                                   >> 448       G4double finalPx = TotalMomentum*ParticleDirection.x()
                                                   >> 449                         - DeltaTotalMomentum*DeltaDirection.x(); 
                                                   >> 450       G4double finalPy = TotalMomentum*ParticleDirection.y()
                                                   >> 451                         - DeltaTotalMomentum*DeltaDirection.y(); 
                                                   >> 452       G4double finalPz = TotalMomentum*ParticleDirection.z()
                                                   >> 453                         - DeltaTotalMomentum*DeltaDirection.z(); 
                                                   >> 454       G4double finalMomentum =
                                                   >> 455                 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
                                                   >> 456       finalPx /= finalMomentum ;
                                                   >> 457       finalPy /= finalMomentum ;
                                                   >> 458       finalPz /= finalMomentum ;
                                                   >> 459 
                                                   >> 460       aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
                                                   >> 461     }
                                                   >> 462   else
                                                   >> 463     {
                                                   >> 464       finalKineticEnergy = 0.;
                                                   >> 465       Edep = finalKineticEnergy ;
                                                   >> 466       if (Charge < 0.) aParticleChange.SetStatusChange(fStopAndKill);
                                                   >> 467       else             aParticleChange.SetStatusChange(fStopButAlive);
                                                   >> 468     }
                                                   >> 469       
                                                   >> 470   aParticleChange.SetEnergyChange( finalKineticEnergy );
                                                   >> 471   aParticleChange.SetNumberOfSecondaries(1);  
                                                   >> 472   aParticleChange.AddSecondary( theDeltaRay );
                                                   >> 473   aParticleChange.SetLocalEnergyDeposit (Edep);
                                                   >> 474       
                                                   >> 475   return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
102 }                                                 476 }
103                                                   477 
104 //....oooOO0OOooo........oooOO0OOooo........oo    478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105                                                   479 
106 void G4eIonisation::ProcessDescription(std::os << 480 void G4eIonisation::PrintInfoDefinition()
107 {                                                 481 {
108   out << "  Ionisation";                       << 482   G4String comments = "delta cross sections from Moller+Bhabha. ";
109   G4VEnergyLossProcess::ProcessDescription(out << 483            comments += "Good description from 1 KeV to 100 GeV.\n";
110 }                                              << 484            comments += "        delta ray energy sampled from  differential Xsection.";
                                                   >> 485                      
                                                   >> 486   G4cout << G4endl << GetProcessName() << ":  " << comments
                                                   >> 487          << "\n        PhysicsTables from " << G4BestUnit(LowerBoundLambda,"Energy")
                                                   >> 488          << " to " << G4BestUnit(UpperBoundLambda,"Energy") 
                                                   >> 489          << " in " << NbinLambda << " bins. \n";
                                                   >> 490 }         
111                                                   491 
112 //....oooOO0OOooo........oooOO0OOooo........oo << 492 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113                                                   493