Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4GammaConversion.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/G4GammaConversion.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4GammaConversion.cc (Version 4.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 //                                                 23 //
 27 //                                             <<  24 // $Id: G4GammaConversion.cc,v 1.14 2001/10/01 15:00:29 maire Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-04-00 $
                                                   >>  26 //
 28 //------------------ G4GammaConversion physics     27 //------------------ G4GammaConversion physics process -------------------------
 29 //                   by Michel Maire, 24 May 1     28 //                   by Michel Maire, 24 May 1996
 30 //                                             <<  29 // 
 31 // Modified by Michel Maire and Vladimir Ivanc <<  30 // 11-06-96, Added SelectRandomAtom() method, M.Maire
 32 //                                             <<  31 // 21-06-96, SetCuts implementation, M.Maire
                                                   >>  32 // 24-06-96, simplification in ComputeCrossSectionPerAtom, M.Maire
                                                   >>  33 // 24-06-96, in DoIt : change the particleType stuff, M.Maire
                                                   >>  34 // 25-06-96, modification in the generation of the teta angle, M.Maire
                                                   >>  35 // 16-09-96, minors optimisations in DoIt. Thanks to P.Urban
                                                   >>  36 //           dynamical array PartialSumSigma
                                                   >>  37 // 13-12-96, fast sampling of epsil below 2 MeV, L.Urban
                                                   >>  38 // 14-01-97, crossection table + meanfreepath table.
                                                   >>  39 //           PartialSumSigma removed, M.Maire
                                                   >>  40 // 14-01-97, in DoIt the positron is always created, even with Ekine=0,
                                                   >>  41 //           for further annihilation, M.Maire
                                                   >>  42 // 14-03-97, new Physics scheme for geant4alpha, M.Maire
                                                   >>  43 // 28-03-97, protection in BuildPhysicsTable, M.Maire
                                                   >>  44 // 19-06-97, correction in ComputeCrossSectionPerAtom, L.Urban
                                                   >>  45 // 04-06-98, in DoIt, secondary production condition:
                                                   >>  46 //             range>G4std::min(threshold,safety)
                                                   >>  47 // 13-08-98, new methods SetBining() PrintInfo()
                                                   >>  48 // 28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation
                                                   >>  49 // 11-07-01, PostStepDoIt - sampling epsil: power(rndm,0.333333)
                                                   >>  50 // 13-07-01, DoIt: suppression of production cut for the (e-,e+) (mma)
                                                   >>  51 // 06-08-01, new methods Store/Retrieve PhysicsTable (mma)
                                                   >>  52 // 06-08-01, BuildThePhysicsTable() called from constructor (mma)
                                                   >>  53 // 17-09-01, migration of Materials to pure STL (mma)
                                                   >>  54 // 20-09-01, DoIt: fminimalEnergy = 1*eV (mma)
                                                   >>  55 // 01-10-01, come back to BuildPhysicsTable(const G4ParticleDefinition&)         
 33 // -------------------------------------------     56 // -----------------------------------------------------------------------------
 34                                                    57 
 35 #include "G4GammaConversion.hh"                    58 #include "G4GammaConversion.hh"
 36 #include "G4PhysicalConstants.hh"              <<  59 #include "G4UnitsTable.hh"
 37 #include "G4SystemOfUnits.hh"                  << 
 38 #include "G4PairProductionRelModel.hh"         << 
 39 #include "G4Electron.hh"                       << 
 40 #include "G4EmParameters.hh"                   << 
 41                                                    60 
 42 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  62  
                                                   >>  63 // constructor
                                                   >>  64  
                                                   >>  65 G4GammaConversion::G4GammaConversion(const G4String& processName)
                                                   >>  66   : G4VDiscreteProcess (processName),            // initialization
                                                   >>  67     theCrossSectionTable(NULL),
                                                   >>  68     theMeanFreePathTable(NULL),  
                                                   >>  69     LowestEnergyLimit (2*electron_mass_c2),
                                                   >>  70     HighestEnergyLimit(100*GeV),
                                                   >>  71     NumbBinTable(100),
                                                   >>  72     fminimalEnergy(1*eV)
                                                   >>  73 {}
                                                   >>  74 
                                                   >>  75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  76  
                                                   >>  77 // destructor
                                                   >>  78  
                                                   >>  79 G4GammaConversion::~G4GammaConversion()
                                                   >>  80 {
                                                   >>  81    if (theCrossSectionTable) {
                                                   >>  82       theCrossSectionTable->clearAndDestroy();
                                                   >>  83       delete theCrossSectionTable;
                                                   >>  84    }
                                                   >>  85 
                                                   >>  86    if (theMeanFreePathTable) {
                                                   >>  87       theMeanFreePathTable->clearAndDestroy();
                                                   >>  88       delete theMeanFreePathTable;
                                                   >>  89    }
                                                   >>  90 }
                                                   >>  91  
                                                   >>  92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 43                                                    93 
 44 G4GammaConversion::G4GammaConversion(const G4S <<  94 
 45   G4ProcessType type):G4VEmProcess (processNam <<  95 void G4GammaConversion::SetPhysicsTableBining(
                                                   >>  96                                   G4double lowE, G4double highE, G4int nBins)
                                                   >>  97 {
                                                   >>  98   LowestEnergyLimit = lowE; HighestEnergyLimit = highE; NumbBinTable = nBins;
                                                   >>  99 }  
                                                   >> 100 
                                                   >> 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 102  
                                                   >> 103 void G4GammaConversion::BuildPhysicsTable(const G4ParticleDefinition&)
                                                   >> 104 // Build cross section and mean free path tables
 46 {                                                 105 {
 47   SetMinKinEnergy(2.0*CLHEP::electron_mass_c2) << 106    G4double LowEdgeEnergy, Value;
 48   SetProcessSubType(fGammaConversion);         << 107    G4PhysicsLogVector* ptrVector;
 49   SetStartFromNullFlag(true);                  << 108 
 50   SetBuildTableFlag(true);                     << 109 // Build cross section per atom tables for the e+e- pair creation
 51   SetSecondaryParticle(G4Electron::Electron()) << 110 
 52   SetLambdaBinning(220);                       << 111    if (theCrossSectionTable) {
                                                   >> 112          theCrossSectionTable->clearAndDestroy(); delete theCrossSectionTable;}
                                                   >> 113 
                                                   >> 114    theCrossSectionTable = new G4PhysicsTable(G4Element::GetNumberOfElements());
                                                   >> 115    const G4ElementTable* theElementTable = G4Element::GetElementTable();
                                                   >> 116    G4double AtomicNumber;
                                                   >> 117    size_t J;
                                                   >> 118 
                                                   >> 119    for ( J=0 ; J < G4Element::GetNumberOfElements(); J++ )  
                                                   >> 120        { 
                                                   >> 121         //create physics vector then fill it ....
                                                   >> 122         ptrVector = new G4PhysicsLogVector(LowestEnergyLimit,HighestEnergyLimit,
                                                   >> 123                                            NumbBinTable );
                                                   >> 124         AtomicNumber = (*theElementTable)[J]->GetZ();
                                                   >> 125  
                                                   >> 126         for ( G4int i = 0 ; i < NumbBinTable ; i++ )      
                                                   >> 127            {
                                                   >> 128              LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
                                                   >> 129              Value = ComputeCrossSectionPerAtom( LowEdgeEnergy, AtomicNumber);  
                                                   >> 130              ptrVector->PutValue( i , Value ) ;
                                                   >> 131            }
                                                   >> 132 
                                                   >> 133         theCrossSectionTable->insertAt( J , ptrVector ) ;
                                                   >> 134 
                                                   >> 135       }
                                                   >> 136 
                                                   >> 137 // Build mean free path table for the e+e- pair creation
                                                   >> 138 
                                                   >> 139    if (theMeanFreePathTable) 
                                                   >> 140      { theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
                                                   >> 141 
                                                   >> 142    theMeanFreePathTable= new G4PhysicsTable(G4Material::GetNumberOfMaterials());
                                                   >> 143    const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
                                                   >> 144    G4Material* material;
                                                   >> 145 
                                                   >> 146    for ( J=0 ; J < G4Material::GetNumberOfMaterials(); J++ )
                                                   >> 147        { 
                                                   >> 148         //create physics vector then fill it ....
                                                   >> 149         ptrVector = new G4PhysicsLogVector(LowestEnergyLimit,HighestEnergyLimit,
                                                   >> 150                                            NumbBinTable);
                                                   >> 151         material = (*theMaterialTable)[J];
                                                   >> 152  
                                                   >> 153         for ( G4int i = 0 ; i < NumbBinTable ; i++ )      
                                                   >> 154            {
                                                   >> 155              LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
                                                   >> 156              Value = ComputeMeanFreePath( LowEdgeEnergy, material);  
                                                   >> 157              ptrVector->PutValue( i , Value ) ;
                                                   >> 158            }
                                                   >> 159 
                                                   >> 160         theMeanFreePathTable->insertAt( J , ptrVector ) ;
                                                   >> 161 
                                                   >> 162       }
                                                   >> 163    
                                                   >> 164    PrintInfoDefinition();        
 53 }                                                 165 }
 54                                                   166 
 55 //....oooOO0OOooo........oooOO0OOooo........oo    167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 168 
                                                   >> 169 G4double G4GammaConversion::ComputeCrossSectionPerAtom
                                                   >> 170                               (G4double GammaEnergy, G4double AtomicNumber)
 56                                                   171  
 57 G4GammaConversion::~G4GammaConversion() = defa << 172 // Calculates the microscopic cross section in GEANT4 internal units.
                                                   >> 173 // A parametrized formula from L. Urban is used to estimate
                                                   >> 174 // the total cross section.
                                                   >> 175 // It gives a good description of the data from 1.5 MeV to 100 GeV.
                                                   >> 176 // below 1.5 MeV: sigma=sigma(1.5MeV)*(GammaEnergy-2electronmass)
                                                   >> 177 //                                   *(GammaEnergy-2electronmass) 
                                                   >> 178 
                                                   >> 179 {
                                                   >> 180  G4double GammaEnergyLimit = 1.5*MeV;
                                                   >> 181  G4double CrossSection = 0.0 ;
                                                   >> 182  if ( AtomicNumber < 1. ) return CrossSection;
                                                   >> 183  if ( GammaEnergy < 2*electron_mass_c2 ) return CrossSection;
                                                   >> 184  
                                                   >> 185  static const G4double
                                                   >> 186     a0= 8.7842e+2*microbarn, a1=-1.9625e+3*microbarn, a2= 1.2949e+3*microbarn,
                                                   >> 187     a3=-2.0028e+2*microbarn, a4= 1.2575e+1*microbarn, a5=-2.8333e-1*microbarn;
                                                   >> 188 
                                                   >> 189  static const G4double
                                                   >> 190     b0=-1.0342e+1*microbarn, b1= 1.7692e+1*microbarn, b2=-8.2381   *microbarn,
                                                   >> 191     b3= 1.3063   *microbarn, b4=-9.0815e-2*microbarn, b5= 2.3586e-3*microbarn;
 58                                                   192 
 59 //....oooOO0OOooo........oooOO0OOooo........oo << 193  static const G4double
                                                   >> 194     c0=-4.5263e+2*microbarn, c1= 1.1161e+3*microbarn, c2=-8.6749e+2*microbarn,
                                                   >> 195     c3= 2.1773e+2*microbarn, c4=-2.0467e+1*microbarn, c5= 6.5372e-1*microbarn;
 60                                                   196 
 61 G4bool G4GammaConversion::IsApplicable(const G << 197  G4double GammaEnergySave = GammaEnergy ;
                                                   >> 198  if (GammaEnergy < GammaEnergyLimit) GammaEnergy = GammaEnergyLimit ;
                                                   >> 199 
                                                   >> 200  G4double X=log(GammaEnergy/electron_mass_c2),X2=X*X, X3=X2*X, X4=X3*X, X5=X4*X;
                                                   >> 201 
                                                   >> 202  G4double F1 = a0 + a1*X + a2*X2 + a3*X3 + a4*X4 + a5*X5,
                                                   >> 203           F2 = b0 + b1*X + b2*X2 + b3*X3 + b4*X4 + b5*X5,
                                                   >> 204           F3 = c0 + c1*X + c2*X2 + c3*X3 + c4*X4 + c5*X5;     
                                                   >> 205 
                                                   >> 206  CrossSection = (AtomicNumber+1.)*
                                                   >> 207                 (F1*AtomicNumber + F2*AtomicNumber*AtomicNumber + F3);
                                                   >> 208 
                                                   >> 209  if (GammaEnergySave < GammaEnergyLimit)
                                                   >> 210    {
                                                   >> 211      X=GammaEnergySave-2.*electron_mass_c2;
                                                   >> 212      CrossSection *= X*X;
                                                   >> 213    }
                                                   >> 214 
                                                   >> 215  if (CrossSection < 0.) CrossSection = 0.;
                                                   >> 216 
                                                   >> 217  return CrossSection;
                                                   >> 218 }
                                                   >> 219 
                                                   >> 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 221  
                                                   >> 222 G4VParticleChange* G4GammaConversion::PostStepDoIt(const G4Track& aTrack,
                                                   >> 223                                                   const G4Step&  aStep)
                                                   >> 224 //
                                                   >> 225 // The secondaries e+e- energies are sampled using the Bethe - Heitler
                                                   >> 226 // cross sections with Coulomb correction.
                                                   >> 227 // A modified version of the random number techniques of Butcher & Messel
                                                   >> 228 // is used (Nuc Phys 20(1960),15).
                                                   >> 229 //
                                                   >> 230 // GEANT4 internal units.
                                                   >> 231 //
                                                   >> 232 // Note 1 : Effects due to the breakdown of the Born approximation at
                                                   >> 233 //          low energy are ignored.
                                                   >> 234 // Note 2 : The differential cross section implicitly takes account of 
                                                   >> 235 //          pair creation in both nuclear and atomic electron fields.
                                                   >> 236 //          However triplet prodution is not generated.
                                                   >> 237  
 62 {                                                 238 {
 63   return (&p == G4Gamma::Gamma());             << 239    aParticleChange.Initialize(aTrack);
                                                   >> 240    G4Material* aMaterial = aTrack.GetMaterial();
                                                   >> 241 
                                                   >> 242    const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
                                                   >> 243    G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
                                                   >> 244    G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
                                                   >> 245 
                                                   >> 246 
                                                   >> 247    G4double epsil ;
                                                   >> 248    G4double epsil0 = electron_mass_c2/GammaEnergy ;
                                                   >> 249 
                                                   >> 250   // do it fast if GammaEnergy < 2. MeV
                                                   >> 251    const G4double Egsmall=2.*MeV; 
                                                   >> 252    if (GammaEnergy<Egsmall) { epsil = epsil0 + (0.5-epsil0)*G4UniformRand(); }
                                                   >> 253 
                                                   >> 254    else
                                                   >> 255    {  // now comes the case with GammaEnergy >= 2. MeV
                                                   >> 256 
                                                   >> 257    // select randomly one element constituing the material  
                                                   >> 258    G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
                                                   >> 259 
                                                   >> 260    // Extract Coulomb factor for this Element
                                                   >> 261    G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
                                                   >> 262    if (GammaEnergy > 50.*MeV) FZ += 8.*(anElement->GetfCoulomb());
                                                   >> 263 
                                                   >> 264    // limits of the screening variable
                                                   >> 265    G4double screenfac = 136.*epsil0/(anElement->GetIonisation()->GetZ3());
                                                   >> 266    G4double screenmax = exp ((42.24 - FZ)/8.368) - 0.952 ;
                                                   >> 267    G4double screenmin = G4std::min(4.*screenfac,screenmax);
                                                   >> 268 
                                                   >> 269    // limits of the energy sampling
                                                   >> 270    G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ;
                                                   >> 271    G4double epsilmin = G4std::max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin;
                                                   >> 272 
                                                   >> 273    //
                                                   >> 274    // sample the energy rate of the created electron (or positron) 
                                                   >> 275    //
                                                   >> 276    //G4double epsil, screenvar, greject ;
                                                   >> 277    G4double  screenvar, greject ;
                                                   >> 278 
                                                   >> 279    G4double F10 = ScreenFunction1(screenmin) - FZ;
                                                   >> 280    G4double F20 = ScreenFunction2(screenmin) - FZ;
                                                   >> 281    G4double NormF1 = G4std::max(F10*epsilrange*epsilrange,0.); 
                                                   >> 282    G4double NormF2 = G4std::max(1.5*F20,0.);
                                                   >> 283 
                                                   >> 284    do {
                                                   >> 285         if ( NormF1/(NormF1+NormF2) > G4UniformRand() )
                                                   >> 286              { epsil = 0.5 - epsilrange*pow(G4UniformRand(), 0.333333);
                                                   >> 287                screenvar = screenfac/(epsil*(1-epsil));
                                                   >> 288                greject = (ScreenFunction1(screenvar) - FZ)/F10;
                                                   >> 289              } 
                                                   >> 290         else { epsil = epsilmin + epsilrange*G4UniformRand();
                                                   >> 291                screenvar = screenfac/(epsil*(1-epsil));
                                                   >> 292                greject = (ScreenFunction2(screenvar) - FZ)/F20;
                                                   >> 293              }
                                                   >> 294 
                                                   >> 295    } while( greject < G4UniformRand() );
                                                   >> 296 
                                                   >> 297    }   //  end of epsil sampling
                                                   >> 298    
                                                   >> 299    //
                                                   >> 300    // fixe charges randomly
                                                   >> 301    //
                                                   >> 302 
                                                   >> 303    G4double ElectTotEnergy, PositTotEnergy;
                                                   >> 304    if (RandBit::shootBit())
                                                   >> 305      {
                                                   >> 306        ElectTotEnergy = (1.-epsil)*GammaEnergy;
                                                   >> 307        PositTotEnergy = epsil*GammaEnergy;
                                                   >> 308      }
                                                   >> 309    else
                                                   >> 310      {
                                                   >> 311        PositTotEnergy = (1.-epsil)*GammaEnergy;
                                                   >> 312        ElectTotEnergy = epsil*GammaEnergy;
                                                   >> 313      }
                                                   >> 314 
                                                   >> 315    //
                                                   >> 316    // scattered electron (positron) angles. ( Z - axis along the parent photon)
                                                   >> 317    //
                                                   >> 318    //  universal distribution suggested by L. Urban 
                                                   >> 319    // (Geant3 manual (1993) Phys211),
                                                   >> 320    //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
                                                   >> 321 
                                                   >> 322    G4double u;
                                                   >> 323    const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
                                                   >> 324 
                                                   >> 325    if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
                                                   >> 326    else                            u= - log(G4UniformRand()*G4UniformRand())/a2;
                                                   >> 327 
                                                   >> 328    G4double Teta = u*electron_mass_c2/GammaEnergy;
                                                   >> 329    G4double Phi  = twopi * G4UniformRand();
                                                   >> 330    G4double dirx=sin(Teta)*cos(Phi), diry=sin(Teta)*sin(Phi), dirz=cos(Teta);
                                                   >> 331  
                                                   >> 332    //
                                                   >> 333    // kinematic of the created pair
                                                   >> 334    //
                                                   >> 335    // the electron and positron are assumed to have a symetric
                                                   >> 336    // angular distribution with respect to the Z axis along the parent photon.
                                                   >> 337 
                                                   >> 338    aParticleChange.SetNumberOfSecondaries(2) ; 
                                                   >> 339 
                                                   >> 340    G4double ElectKineEnergy = G4std::max(0.,ElectTotEnergy - electron_mass_c2);
                                                   >> 341    G4double localEnergyDeposit = 0.;
                                                   >> 342 
                                                   >> 343    if (ElectKineEnergy > fminimalEnergy)
                                                   >> 344      {
                                                   >> 345        G4ThreeVector ElectDirection (dirx, diry, dirz);
                                                   >> 346        ElectDirection.rotateUz(GammaDirection);   
                                                   >> 347  
                                                   >> 348        // create G4DynamicParticle object for the particle1  
                                                   >> 349        G4DynamicParticle* aParticle1= new G4DynamicParticle(
                                                   >> 350                         G4Electron::Electron(),ElectDirection,ElectKineEnergy);
                                                   >> 351        aParticleChange.AddSecondary(aParticle1); 
                                                   >> 352      }
                                                   >> 353    else
                                                   >> 354      { localEnergyDeposit += ElectKineEnergy;}   
                                                   >> 355 
                                                   >> 356    // the e+ is always created (even with Ekine=0) for further annihilation.
                                                   >> 357 
                                                   >> 358    G4double PositKineEnergy = G4std::max(0.,PositTotEnergy - electron_mass_c2);
                                                   >> 359    if (PositKineEnergy < fminimalEnergy)
                                                   >> 360      { localEnergyDeposit += PositKineEnergy; PositKineEnergy = 0.;}
                                                   >> 361 
                                                   >> 362    G4ThreeVector PositDirection (-dirx, -diry, dirz);
                                                   >> 363    PositDirection.rotateUz(GammaDirection);   
                                                   >> 364  
                                                   >> 365    // create G4DynamicParticle object for the particle2 
                                                   >> 366    G4DynamicParticle* aParticle2= new G4DynamicParticle(
                                                   >> 367                       G4Positron::Positron(),PositDirection,PositKineEnergy);
                                                   >> 368    aParticleChange.AddSecondary(aParticle2); 
                                                   >> 369 
                                                   >> 370    aParticleChange.SetLocalEnergyDeposit(localEnergyDeposit);
                                                   >> 371 
                                                   >> 372    //
                                                   >> 373    // Kill the incident photon 
                                                   >> 374    //
                                                   >> 375 
                                                   >> 376    aParticleChange.SetMomentumChange( 0., 0., 0. );
                                                   >> 377    aParticleChange.SetEnergyChange( 0. ); 
                                                   >> 378    aParticleChange.SetStatusChange( fStopAndKill );
                                                   >> 379 
                                                   >> 380    //  Reset NbOfInteractionLengthLeft and return aParticleChange
                                                   >> 381    return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
 64 }                                                 382 }
 65                                                   383 
 66 //....oooOO0OOooo........oooOO0OOooo........oo << 384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67                                                   385 
 68 void G4GammaConversion::InitialiseProcess(cons << 386 G4Element* G4GammaConversion::SelectRandomAtom(
                                                   >> 387                                          const G4DynamicParticle* aDynamicGamma,
                                                   >> 388                                                G4Material* aMaterial)
 69 {                                                 389 {
 70   if(!isInitialised) {                         << 390   // select randomly 1 element within the material 
 71     isInitialised = true;                      << 
 72     G4EmParameters* param = G4EmParameters::In << 
 73     G4double emin = std::max(param->MinKinEner << 
 74     G4double emax = param->MaxKinEnergy();     << 
 75                                                   391 
 76     SetMinKinEnergy(emin);                     << 392   const G4int NumberOfElements            = aMaterial->GetNumberOfElements();
                                                   >> 393   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
                                                   >> 394   if (NumberOfElements == 1) return (*theElementVector)[0];
 77                                                   395 
 78     if(nullptr == EmModel(0)) { SetEmModel(new << 396   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
 79     EmModel(0)->SetLowEnergyLimit(emin);       << 397 
 80     EmModel(0)->SetHighEnergyLimit(emax);      << 398   G4double PartialSumSigma = 0. ;
 81     AddEmModel(1, EmModel(0));                 << 399   G4double rval = G4UniformRand()/MeanFreePath;
 82   }                                            << 400  
                                                   >> 401   for ( G4int i=0 ; i < NumberOfElements ; i++ )
                                                   >> 402       { PartialSumSigma += NbOfAtomsPerVolume[i] *
                                                   >> 403                   GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]);
                                                   >> 404         if (rval <= PartialSumSigma) return ((*theElementVector)[i]);
                                                   >> 405       }
                                                   >> 406   G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
                                                   >> 407        << "' has no elements, NULL pointer returned." << G4endl;
                                                   >> 408   return NULL;
 83 }                                                 409 }
 84                                                   410 
 85 //....oooOO0OOooo........oooOO0OOooo........oo << 411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 86                                                   412 
 87 G4double G4GammaConversion::MinPrimaryEnergy(c << 413 G4bool G4GammaConversion::StorePhysicsTable(G4ParticleDefinition* particle,
 88                const G4Material*)              << 414                       const G4String& directory, 
                                                   >> 415                       G4bool          ascii)
 89 {                                                 416 {
 90   return 2*CLHEP::electron_mass_c2;            << 417   G4String filename;
                                                   >> 418 
                                                   >> 419   // store cross section table
                                                   >> 420   filename = GetPhysicsTableFileName(particle,directory,"CrossSection",ascii);
                                                   >> 421   if ( !theCrossSectionTable->StorePhysicsTable(filename, ascii) ){
                                                   >> 422     G4cout << " FAIL theCrossSectionTable->StorePhysicsTable in " << filename
                                                   >> 423            << G4endl;
                                                   >> 424     return false;
                                                   >> 425   }
                                                   >> 426 
                                                   >> 427   // store mean free path table
                                                   >> 428   filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
                                                   >> 429   if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
                                                   >> 430     G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
                                                   >> 431            << G4endl;
                                                   >> 432     return false;
                                                   >> 433   }
                                                   >> 434   
                                                   >> 435   G4cout << GetProcessName() << " for " << particle->GetParticleName()
                                                   >> 436          << ": Success to store the PhysicsTables in "  
                                                   >> 437          << directory << G4endl;
                                                   >> 438   return true;
 91 }                                                 439 }
 92                                                   440 
 93 //....oooOO0OOooo........oooOO0OOooo........oo    441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 94                                                   442 
 95 void G4GammaConversion::ProcessDescription(std << 443 G4bool G4GammaConversion::RetrievePhysicsTable(G4ParticleDefinition* particle,
                                                   >> 444                    const G4String& directory, 
                                                   >> 445                          G4bool          ascii)
 96 {                                                 446 {
 97   out << "  Gamma conversion";                 << 447   // delete theCrossSectionTable and theMeanFreePathTable
 98   G4VEmProcess::ProcessDescription(out);       << 448   if (theCrossSectionTable != 0) {
                                                   >> 449     theCrossSectionTable->clearAndDestroy();
                                                   >> 450     delete theCrossSectionTable;
                                                   >> 451   }
                                                   >> 452   if (theMeanFreePathTable != 0) {
                                                   >> 453     theMeanFreePathTable->clearAndDestroy();
                                                   >> 454     delete theMeanFreePathTable;
                                                   >> 455   }
                                                   >> 456 
                                                   >> 457   G4String filename;
                                                   >> 458 
                                                   >> 459   // retreive cross section table
                                                   >> 460   filename = GetPhysicsTableFileName(particle,directory,"CrossSection",ascii);
                                                   >> 461   theCrossSectionTable = new G4PhysicsTable(G4Element::GetNumberOfElements());
                                                   >> 462   if ( !theCrossSectionTable->RetrievePhysicsTable(filename, ascii) ){
                                                   >> 463     G4cout << " FAIL theCrossSectionTable->RetrievePhysicsTable in " << filename
                                                   >> 464            << G4endl;  
                                                   >> 465     return false;
                                                   >> 466   }
                                                   >> 467 
                                                   >> 468   // retreive mean free path table
                                                   >> 469   filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
                                                   >> 470   theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
                                                   >> 471   if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
                                                   >> 472     G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
                                                   >> 473            << G4endl;  
                                                   >> 474     return false;
                                                   >> 475   }
                                                   >> 476   
                                                   >> 477   G4cout << GetProcessName() << " for " << particle->GetParticleName()
                                                   >> 478          << ": Success to retrieve the PhysicsTables from "
                                                   >> 479          << directory << G4endl;
                                                   >> 480   return true;
 99 }                                                 481 }
                                                   >> 482  
                                                   >> 483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100                                                   484 
101 //....oooOO0OOooo........oooOO0OOooo........oo << 485 void G4GammaConversion::PrintInfoDefinition()
                                                   >> 486 {
                                                   >> 487   G4String comments = "Total cross sections from a parametrisation. ";
                                                   >> 488            comments += "Good description from 1.5 MeV to 100 GeV for all Z. \n";
                                                   >> 489            comments += "        e+e- energies according Bethe-Heitler";
                                                   >> 490                      
                                                   >> 491   G4cout << G4endl << GetProcessName() << ":  " << comments
                                                   >> 492          << "\n        PhysicsTables from " 
                                                   >> 493              << G4BestUnit(LowestEnergyLimit, "Energy")
                                                   >> 494          << " to " << G4BestUnit(HighestEnergyLimit,"Energy") 
                                                   >> 495          << " in " << NumbBinTable << " bins. \n";
                                                   >> 496 }         
                                                   >> 497 
                                                   >> 498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102                                                   499