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 1.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.3.8.1 1999/12/07 20:51:02 gunter Exp $
                                                   >>   9 // GEANT4 tag $Name: geant4-01-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 // --------------------------------------------------------------
 36 //                                             <<  31  
 37 // Modified by Michel Maire and Vladimir Ivanc << 
 38 //                                             << 
 39 // ------------------------------------------- << 
 40 //                                             << 
 41 //....oooOO0OOooo........oooOO0OOooo........oo << 
 42 //....oooOO0OOooo........oooOO0OOooo........oo << 
 43                                                    32 
 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 // constructor and destructor
                                                   >>  39  
                                                   >>  40 G4eIonisation::G4eIonisation(const G4String& processName)
                                                   >>  41    : G4eEnergyLoss(processName),
                                                   >>  42      theMeanFreePathTable(NULL),
                                                   >>  43      LowestKineticEnergy(1.*keV),
                                                   >>  44      HighestKineticEnergy(100.*TeV),
                                                   >>  45      TotBin(100)
                                                   >>  46 { }
 49                                                    47 
 50 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51                                                    49 
 52 G4eIonisation::G4eIonisation(const G4String& n <<  50 G4eIonisation::~G4eIonisation() 
 53   : G4VEnergyLossProcess(name),                << 
 54     theElectron(G4Electron::Electron()),       << 
 55     isElectron(true),                          << 
 56     isInitialised(false)                       << 
 57 {                                                  51 {
 58   SetProcessSubType(fIonisation);              <<  52      if (theMeanFreePathTable) {
 59   SetSecondaryParticle(theElectron);           <<  53         theMeanFreePathTable->clearAndDestroy();
                                                   >>  54         delete theMeanFreePathTable;
                                                   >>  55      }
                                                   >>  56 
 60 }                                                  57 }
 61                                                    58 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63                                                    60 
 64 G4eIonisation::~G4eIonisation() = default;     <<  61 void G4eIonisation::SetPhysicsTableBining(G4double lowE, G4double highE, G4int nBins)
                                                   >>  62 {
                                                   >>  63   LowestKineticEnergy = lowE;  HighestKineticEnergy = highE; TotBin = nBins;
                                                   >>  64 }  
 65                                                    65 
 66 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 67                                                    67 
 68 G4double G4eIonisation::MinPrimaryEnergy(const <<  68 void G4eIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
 69            const G4Material*,                  <<  69 //  just call BuildLossTable+BuildLambdaTable
 70            G4double cut)                       << 
 71 {                                                  70 {
 72   G4double x = cut;                            <<  71  
 73   if(isElectron) x += cut;                     <<  72     BuildLossTable(aParticleType) ;
 74   return x;                                    <<  73 
                                                   >>  74   if(&aParticleType==G4Electron::Electron())
                                                   >>  75   {
                                                   >>  76     RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable ;
                                                   >>  77     CounterOfElectronProcess++;
                                                   >>  78   }
                                                   >>  79   else
                                                   >>  80   {
                                                   >>  81     RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable ;
                                                   >>  82     CounterOfPositronProcess++;
                                                   >>  83   }
                                                   >>  84  
                                                   >>  85     BuildLambdaTable(aParticleType) ;
                                                   >>  86  
                                                   >>  87     BuildDEDXTable(aParticleType) ;
                                                   >>  88                                              
                                                   >>  89   if(&aParticleType==G4Electron::Electron())
                                                   >>  90     PrintInfoDefinition();  
 75 }                                                  91 }
 76                                                    92 
 77 //....oooOO0OOooo........oooOO0OOooo........oo     93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 78                                                    94 
 79 G4bool G4eIonisation::IsApplicable(const G4Par <<  95 void G4eIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
 80 {                                                  96 {
 81   return (&p == theElectron || &p == G4Positro <<  97 // Build tables for the ionization energy loss
                                                   >>  98 //  the tables are built for *MATERIALS*
                                                   >>  99 
                                                   >> 100     const G4double twoln10 = 2.*log(10.);
                                                   >> 101     const G4double Factor = twopi_mc2_rcl2;
                                                   >> 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       // now comes the loop for the kinetic energy values
                                                   >> 146 
                                                   >> 147       for (G4int i = 0 ; i < TotBin ; i++)
                                                   >> 148          {
                                                   >> 149           LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;      
                                                   >> 150           tau = LowEdgeEnergy/ParticleMass ;
                                                   >> 151 
                                                   >> 152           // Seltzer-Berger formula 
                                                   >> 153           gamma = tau + 1.; gamma2 = gamma*gamma; 
                                                   >> 154           bg2 = tau*(tau+2.);
                                                   >> 155           beta2 = bg2/gamma2;
                                                   >> 156 
                                                   >> 157           // electron
                                                   >> 158           if (&aParticleType==G4Electron::Electron())
                                                   >> 159             {
                                                   >> 160               Tmax = LowEdgeEnergy/2.;  
                                                   >> 161               d = min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
                                                   >> 162               ionloss = log(2.*(tau+2.)/Eexcm2)-1.-beta2
                                                   >> 163                        + log((tau-d)*d)+tau/(tau-d)
                                                   >> 164                        + (0.5*d*d+(2.*tau+1.)*log(1.-d/tau))/gamma2;
                                                   >> 165             }
                                                   >> 166           else        //positron
                                                   >> 167             {
                                                   >> 168               Tmax = LowEdgeEnergy ;  
                                                   >> 169               d = min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
                                                   >> 170               d2=d*d/2.; d3=d*d*d/3.; d4=d*d*d*d/4.;
                                                   >> 171               y=1./(1.+gamma);
                                                   >> 172               ionloss = log(2.*(tau+2.)/Eexcm2)+log(tau*d)
                                                   >> 173                        - beta2*(tau+2.*d-y*(3.*d2+y*(d-d3+y*(d2-tau*d3+d4))))/tau;
                                                   >> 174             } 
                                                   >> 175 
                                                   >> 176           //density correction
                                                   >> 177           x = log(bg2)/twoln10;
                                                   >> 178           if (x < X0den) delta = 0.;
                                                   >> 179           else { delta = twoln10*x - Cden;
                                                   >> 180                  if (x < X1den) delta += Aden*pow((X1den-x),Mden);
                                                   >> 181                } 
                                                   >> 182 
                                                   >> 183           //now you can compute the total ionization loss
                                                   >> 184           ionloss -= delta ;
                                                   >> 185           ionloss *= Factor*ElectronDensity/beta2 ;
                                                   >> 186           if (ionloss <= 0.) ionloss = 0.;
                                                   >> 187    
                                                   >> 188           aVector->PutValue(i,ionloss) ;
                                                   >> 189          }          
                                                   >> 190       theLossTable->insert(aVector);
                                                   >> 191     }
 82 }                                                 192 }
 83                                                   193 
 84 //....oooOO0OOooo........oooOO0OOooo........oo    194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 85                                                   195 
 86 void G4eIonisation::InitialiseEnergyLossProces << 196 void G4eIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
 87         const G4ParticleDefinition* part,      << 
 88         const G4ParticleDefinition*)           << 
 89 {                                                 197 {
 90   if(!isInitialised) {                         << 198   // Build mean free path tables for the delta ray production process
 91     if(part != theElectron) { isElectron = fal << 199   //     tables are built for MATERIALS 
 92     if (nullptr == EmModel(0)) { SetEmModel(ne << 200 
 93     G4EmParameters* param = G4EmParameters::In << 201   G4double LowEdgeEnergy, Value, SIGMA;
 94     EmModel(0)->SetLowEnergyLimit(param->MinKi << 202 
 95     EmModel(0)->SetHighEnergyLimit(param->MaxK << 203   //create table
 96     if (nullptr == FluctModel()) {             << 204   const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
 97       SetFluctModel(G4EmStandUtil::ModelOfFluc << 205   G4int numOfMaterials = theMaterialTable->length();
                                                   >> 206 
                                                   >> 207   if (theMeanFreePathTable) { theMeanFreePathTable->clearAndDestroy();
                                                   >> 208                               delete theMeanFreePathTable;
                                                   >> 209                             }
                                                   >> 210 
                                                   >> 211   theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
                                                   >> 212 
                                                   >> 213   // get electron  cuts in kinetic energy
                                                   >> 214   G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
                                                   >> 215  
                                                   >> 216 
                                                   >> 217   // loop for materials 
                                                   >> 218 
                                                   >> 219  for (G4int J=0 ; J < numOfMaterials; J++)
                                                   >> 220     { 
                                                   >> 221      //create physics vector then fill it ....
                                                   >> 222 
                                                   >> 223      G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
                                                   >> 224                LowestKineticEnergy, HighestKineticEnergy, TotBin);
                                                   >> 225 
                                                   >> 226      // compute the (macroscopic) cross section first
                                                   >> 227  
                                                   >> 228      const G4Material* material= (*theMaterialTable)[J];        
                                                   >> 229      const 
                                                   >> 230      G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 231      const
                                                   >> 232      G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
                                                   >> 233      const
                                                   >> 234      G4int NumberOfElements = material->GetNumberOfElements() ;
                                                   >> 235  
                                                   >> 236      // get the electron kinetic energy cut for the actual material,
                                                   >> 237      // it will be used in ComputeMicroscopicCrossSection
                                                   >> 238      // (--> it will be the same for all the elements in this material )
                                                   >> 239      G4double DeltaThreshold = DeltaCutInKineticEnergy[J] ;
                                                   >> 240 
                                                   >> 241      for (G4int i = 0 ; i < TotBin ; i++)
                                                   >> 242         {
                                                   >> 243           LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
                                                   >> 244           SIGMA = 0.;           
                                                   >> 245           for (G4int iel=0; iel<NumberOfElements; iel++ )
                                                   >> 246              {
                                                   >> 247                 SIGMA += theAtomicNumDensityVector[iel]*
                                                   >> 248                          ComputeMicroscopicCrossSection( aParticleType,
                                                   >> 249                                                          LowEdgeEnergy,
                                                   >> 250                                       (*theElementVector)(iel)->GetZ(),
                                                   >> 251                                                        DeltaThreshold);
                                                   >> 252              }
                                                   >> 253 
                                                   >> 254           // mean free path = 1./macroscopic cross section
                                                   >> 255           Value = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;     
                                                   >> 256           aVector->PutValue(i, Value) ;
                                                   >> 257         }
                                                   >> 258      theMeanFreePathTable->insert(aVector);
 98     }                                             259     }
 99     AddEmModel(1, EmModel(), FluctModel());    << 260 
100     isInitialised = true;                      << 
101   }                                            << 
102 }                                                 261 }
103                                                   262 
104 //....oooOO0OOooo........oooOO0OOooo........oo    263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105                                                   264 
106 void G4eIonisation::ProcessDescription(std::os << 265 G4double G4eIonisation::ComputeMicroscopicCrossSection(
                                                   >> 266                                  const G4ParticleDefinition& aParticleType,
                                                   >> 267                                  G4double KineticEnergy,
                                                   >> 268                                  G4double AtomicNumber ,
                                                   >> 269                                  G4double DeltaThreshold)
                                                   >> 270 {
                                                   >> 271   // calculates the microscopic cross section
                                                   >> 272   //(it is called for elements , AtomicNumber = Z )
                                                   >> 273  
                                                   >> 274   G4double MaxKineticEnergyTransfer, TotalCrossSection(0.);
                                                   >> 275   
                                                   >> 276   ParticleMass = aParticleType.GetPDGMass();
                                                   >> 277   G4double TotalEnergy = KineticEnergy + ParticleMass;
                                                   >> 278 
                                                   >> 279   G4double betasquare = KineticEnergy*(TotalEnergy+ParticleMass)
                                                   >> 280                        /(TotalEnergy*TotalEnergy);
                                                   >> 281   G4double gamma = TotalEnergy/ParticleMass, gamma2 = gamma*gamma;
                                                   >> 282   G4double x=DeltaThreshold/KineticEnergy, x2 = x*x;
                                                   >> 283 
                                                   >> 284   if (&aParticleType==G4Electron::Electron())
                                                   >> 285                             MaxKineticEnergyTransfer = 0.5*KineticEnergy;
                                                   >> 286   else                      MaxKineticEnergyTransfer =     KineticEnergy;
                                                   >> 287 
                                                   >> 288   // now you can calculate the total cross section
                                                   >> 289 
                                                   >> 290  if (MaxKineticEnergyTransfer > DeltaThreshold)
                                                   >> 291    {
                                                   >> 292     if (&aParticleType==G4Electron::Electron())   //Moller (e-e-) scattering
                                                   >> 293       {
                                                   >> 294         TotalCrossSection  = (gamma-1.)*(gamma-1.)*(0.5-x)/gamma2 + 1./x
                                                   >> 295                             - 1./(1.-x)-(2.*gamma-1.)*log((1.-x)/x)/gamma2; 
                                                   >> 296         TotalCrossSection /= betasquare; 
                                                   >> 297       }
                                                   >> 298     else                                         //Bhabha (e+e-) scattering
                                                   >> 299       {
                                                   >> 300        G4double y=1./(1.+gamma), y2 =y*y, y12=1.-2.*y;
                                                   >> 301        G4double b1=2.-y2, b2=y12*(3.+y2), b4=y12*y12*y12, b3=b4+y12*y12;
                                                   >> 302        TotalCrossSection = (1./x-1.)/betasquare+b1*log(x)+b2*(1.-x)
                                                   >> 303                           - b3*(1.-x2)/2.+b4*(1.-x2*x)/3.;
                                                   >> 304       }   
                                                   >> 305     TotalCrossSection *= (twopi_mc2_rcl2*AtomicNumber/KineticEnergy);
                                                   >> 306    }
                                                   >> 307  return TotalCrossSection ;
                                                   >> 308 }
                                                   >> 309 
                                                   >> 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 
                                                   >> 311  
                                                   >> 312 G4VParticleChange* G4eIonisation::PostStepDoIt( const G4Track& trackData,   
                                                   >> 313                                                 const G4Step&  stepData)         
107 {                                                 314 {
108   out << "  Ionisation";                       << 315   aParticleChange.Initialize(trackData) ;
109   G4VEnergyLossProcess::ProcessDescription(out << 316   
                                                   >> 317   G4Material*               aMaterial = trackData.GetMaterial() ;
                                                   >> 318   const G4DynamicParticle*  aParticle = trackData.GetDynamicParticle() ;
                                                   >> 319 
                                                   >> 320   ParticleMass = aParticle->GetDefinition()->GetPDGMass();
                                                   >> 321   G4double KineticEnergy = aParticle->GetKineticEnergy();
                                                   >> 322   G4double TotalEnergy = KineticEnergy + ParticleMass;
                                                   >> 323   G4double Psquare = KineticEnergy*(TotalEnergy+ParticleMass);
                                                   >> 324   G4double TotalMomentum = sqrt(Psquare);
                                                   >> 325   G4double Esquare=TotalEnergy*TotalEnergy;
                                                   >> 326   G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
                                                   >> 327 
                                                   >> 328   //  get kinetic energy cut for the electron
                                                   >> 329   G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
                                                   >> 330   G4double DeltaThreshold = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
                                                   >> 331 
                                                   >> 332   // some kinematics
                                                   >> 333   G4double MaxKineticEnergyTransfer;
                                                   >> 334   if (Charge < 0.) MaxKineticEnergyTransfer = 0.5*KineticEnergy;
                                                   >> 335   else             MaxKineticEnergyTransfer =     KineticEnergy;
                                                   >> 336 
                                                   >> 337   // sampling kinetic energy of the delta ray 
                                                   >> 338 
                                                   >> 339   if (MaxKineticEnergyTransfer <= DeltaThreshold) // pathological case (should not happen,
                                                   >> 340                                                   // there is no change at all)
                                                   >> 341      return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
                                                   >> 342 
                                                   >> 343 
                                                   >> 344   // normal case
                                                   >> 345   G4double cc,y,y2,c2,b0,b1,b2,b3,b4,x,x1,grej,grejc;
                                                   >> 346  
                                                   >> 347   G4double tau = KineticEnergy/ParticleMass;
                                                   >> 348   G4double gamma = tau+1., gamma2=gamma*gamma;
                                                   >> 349   G4double xc = DeltaThreshold/KineticEnergy, xc1=1.-xc;
                                                   >> 350 
                                                   >> 351   if (Charge < 0.)  // Moller (e-e-) scattering
                                                   >> 352     { 
                                                   >> 353       b1=4./(9.*gamma2-10.*gamma+5.);
                                                   >> 354       b2=tau*tau*b1; b3=(2.*gamma2+2.*gamma-1.)*b1;
                                                   >> 355       cc=1.-2.*xc;       
                                                   >> 356       do { 
                                                   >> 357            x    = xc/(1.-cc*G4UniformRand()); x1 = 1.-x;
                                                   >> 358            grej = b2*x*x-b3*x/x1+b1*gamma2/(x1*x1);
                                                   >> 359          } while (G4UniformRand()>grej) ;
                                                   >> 360     }
                                                   >> 361   else             // Bhabha (e+e-) scattering
                                                   >> 362     {
                                                   >> 363       y=1./(gamma+1.); y2=y*y; cc=1.-2.*y;
                                                   >> 364       b1=2.-y2; b2=cc*(3.+y2);
                                                   >> 365       c2=cc*cc; b4=c2*cc; b3=c2+b4;
                                                   >> 366       b0=gamma2/(gamma2-1.);
                                                   >> 367       grejc=(((b4*xc-b3)*xc+b2)*xc-b1)*xc+b0;
                                                   >> 368       do {
                                                   >> 369            x    = xc/(1.-xc1*G4UniformRand());
                                                   >> 370            grej = ((((b4*x-b3)*x+b2)*x-b1)*x+b0)/grejc;
                                                   >> 371          } while (G4UniformRand()>grej);
                                                   >> 372     }
                                                   >> 373  
                                                   >> 374     G4double DeltaKineticEnergy = x * KineticEnergy;
                                                   >> 375 
                                                   >> 376   // protection :do not produce a secondary with 0. kinetic energy !
                                                   >> 377   if (DeltaKineticEnergy <= 0.)
                                                   >> 378       return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
                                                   >> 379 
                                                   >> 380   G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
                                                   >> 381                                                        2. * electron_mass_c2 ));
                                                   >> 382    
                                                   >> 383   G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
                                                   >> 384                       /(DeltaTotalMomentum * TotalMomentum);
                                                   >> 385 
                                                   >> 386   if (costheta < -1.) costheta = -1.;
                                                   >> 387   if (costheta > +1.) costheta = +1.;
                                                   >> 388 
                                                   >> 389   //  direction of the delta electron
                                                   >> 390 
                                                   >> 391   G4double phi = twopi * G4UniformRand(); 
                                                   >> 392   G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
                                                   >> 393   G4double dirx = sintheta * cos(phi), diry = sintheta * sin(phi), dirz = costheta;
                                                   >> 394 
                                                   >> 395   G4ThreeVector DeltaDirection(dirx,diry,dirz);
                                                   >> 396   DeltaDirection.rotateUz(ParticleDirection);
                                                   >> 397 
                                                   >> 398   // create G4DynamicParticle object for delta ray
                                                   >> 399 
                                                   >> 400   G4DynamicParticle* theDeltaRay = new G4DynamicParticle;
                                                   >> 401   theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
                                                   >> 402   theDeltaRay->SetMomentumDirection(
                                                   >> 403                    DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z()); 
                                                   >> 404   theDeltaRay->SetDefinition(G4Electron::Electron());
                                                   >> 405    
                                                   >> 406   // fill aParticleChange 
                                                   >> 407   // changed energy and momentum of the actual particle
                                                   >> 408   G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy;
                                                   >> 409   
                                                   >> 410   G4double Edep = 0. ;
                                                   >> 411 
                                                   >> 412   if (finalKineticEnergy > MinKineticEnergy)
                                                   >> 413     {
                                                   >> 414       G4double finalPx = TotalMomentum*ParticleDirection.x()
                                                   >> 415                         - DeltaTotalMomentum*DeltaDirection.x(); 
                                                   >> 416       G4double finalPy = TotalMomentum*ParticleDirection.y()
                                                   >> 417                         - DeltaTotalMomentum*DeltaDirection.y(); 
                                                   >> 418       G4double finalPz = TotalMomentum*ParticleDirection.z()
                                                   >> 419                         - DeltaTotalMomentum*DeltaDirection.z(); 
                                                   >> 420       G4double finalMomentum =
                                                   >> 421                 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
                                                   >> 422       finalPx /= finalMomentum ;
                                                   >> 423       finalPy /= finalMomentum ;
                                                   >> 424       finalPz /= finalMomentum ;
                                                   >> 425 
                                                   >> 426       aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
                                                   >> 427     }
                                                   >> 428   else
                                                   >> 429     {
                                                   >> 430       finalKineticEnergy = 0.;
                                                   >> 431       Edep = finalKineticEnergy ;
                                                   >> 432       if (Charge < 0.) aParticleChange.SetStatusChange(fStopAndKill);
                                                   >> 433       else             aParticleChange.SetStatusChange(fStopButAlive);
                                                   >> 434     }
                                                   >> 435       
                                                   >> 436   aParticleChange.SetEnergyChange( finalKineticEnergy );
                                                   >> 437   aParticleChange.SetNumberOfSecondaries(1);  
                                                   >> 438   aParticleChange.AddSecondary( theDeltaRay );
                                                   >> 439   aParticleChange.SetLocalEnergyDeposit (Edep);
                                                   >> 440       
                                                   >> 441   return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
110 }                                                 442 }
111                                                   443 
112 //....oooOO0OOooo........oooOO0OOooo........oo << 444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 445 
                                                   >> 446 void G4eIonisation::PrintInfoDefinition()
                                                   >> 447 {
                                                   >> 448   G4String comments = "delta cross sections from Moller+Bhabha. ";
                                                   >> 449            comments += "Good description from 1 KeV to 100 GeV.\n";
                                                   >> 450            comments += "        delta ray energy sampled from  differential Xsection.";
                                                   >> 451                      
                                                   >> 452   G4cout << endl << GetProcessName() << ":  " << comments
                                                   >> 453          << "\n        PhysicsTables from " << G4BestUnit(LowestKineticEnergy,"Energy")
                                                   >> 454          << " to " << G4BestUnit(HighestKineticEnergy,"Energy") 
                                                   >> 455          << " in " << TotBin << " bins. \n";
                                                   >> 456 }         
                                                   >> 457 
                                                   >> 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113                                                   459