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 6.1)


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