Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermoreGammaConversionModel.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/lowenergy/src/G4LivermoreGammaConversionModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermoreGammaConversionModel.cc (Version 9.6.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Author: Sebastien Incerti                       26 // Author: Sebastien Incerti
 27 //         22 January 2012                         27 //         22 January 2012
 28 //         on base of G4LivermoreGammaConversi <<  28 //         on base of G4LivermoreGammaConversionModel
 29 //         and G4LivermoreRayleighModel (MT ve << 
 30 //                                             << 
 31 // Modifications: Zhuxin Li@CENBG              << 
 32 //                11 March 2020                << 
 33 //                derives from G4PairProductio << 
 34 // ------------------------------------------- << 
 35                                                    29 
 36 #include "G4LivermoreGammaConversionModel.hh"      30 #include "G4LivermoreGammaConversionModel.hh"
 37                                                << 
 38 #include "G4AutoLock.hh"                       << 
 39 #include "G4Electron.hh"                       << 
 40 #include "G4EmParameters.hh"                   << 
 41 #include "G4Exp.hh"                            << 
 42 #include "G4ParticleChangeForGamma.hh"         << 
 43 #include "G4PhysicalConstants.hh"                  31 #include "G4PhysicalConstants.hh"
 44 #include "G4PhysicsFreeVector.hh"              << 
 45 #include "G4SystemOfUnits.hh"                      32 #include "G4SystemOfUnits.hh"
 46                                                    33 
 47 namespace                                      << 
 48 {                                              << 
 49 G4Mutex LivermoreGammaConversionModelMutex = G << 
 50 }                                              << 
 51                                                << 
 52 //....oooOO0OOooo........oooOO0OOooo........oo << 
 53 //....oooOO0OOooo........oooOO0OOooo........oo     34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    35 
 55 G4PhysicsFreeVector* G4LivermoreGammaConversio <<  36 using namespace std;
 56 G4String G4LivermoreGammaConversionModel::gDat <<  37 
                                                   >>  38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 57                                                    39 
 58 G4LivermoreGammaConversionModel::G4LivermoreGa <<  40 G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*,
 59                                                <<  41                  const G4String& nam)
 60   : G4PairProductionRelModel(p, nam)           <<  42 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),maxZ(99)
 61 {                                                  43 {
 62   fParticleChange = nullptr;                   <<  44   fParticleChange = 0;
 63   lowEnergyLimit = 2. * CLHEP::electron_mass_c <<  45 
 64   verboseLevel = 0;                            <<  46   lowEnergyLimit = 2.0*electron_mass_c2;
                                                   >>  47   data.resize(maxZ+1,0); 
                                                   >>  48      
                                                   >>  49   verboseLevel= 0;
 65   // Verbosity scale for debugging purposes:       50   // Verbosity scale for debugging purposes:
 66   // 0 = nothing                               <<  51   // 0 = nothing 
 67   // 1 = calculation of cross sections, file o     52   // 1 = calculation of cross sections, file openings...
 68   // 2 = entering in methods                       53   // 2 = entering in methods
 69   if (verboseLevel > 0) {                      <<  54 
                                                   >>  55   if(verboseLevel > 0) 
                                                   >>  56   {
 70     G4cout << "G4LivermoreGammaConversionModel     57     G4cout << "G4LivermoreGammaConversionModel is constructed " << G4endl;
 71   }                                                58   }
 72 }                                                  59 }
 73                                                    60 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    62 
 76 G4LivermoreGammaConversionModel::~G4LivermoreG     63 G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel()
 77 {                                              <<  64 {  
 78   if (IsMaster()) {                            <<  65   for(G4int i=0; i<=maxZ; ++i) { delete data[i]; }
 79     for (G4int i = 0; i <= maxZ; ++i) {        << 
 80       if (data[i]) {                           << 
 81         delete data[i];                        << 
 82         data[i] = nullptr;                     << 
 83       }                                        << 
 84     }                                          << 
 85   }                                            << 
 86 }                                                  66 }
 87                                                    67 
 88 //....oooOO0OOooo........oooOO0OOooo........oo     68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 89                                                    69 
 90 void G4LivermoreGammaConversionModel::Initiali <<  70 void 
 91                                                <<  71 G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
                                                   >>  72               const G4DataVector& cuts)
 92 {                                                  73 {
 93   G4PairProductionRelModel::Initialise(particl <<  74   if (verboseLevel > 1) 
 94   if (verboseLevel > 1) {                      <<  75   {
 95     G4cout << "Calling Initialise() of G4Liver     76     G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel." << G4endl
 96            << "Energy range: " << LowEnergyLim <<  77      << "Energy range: "
 97            << " GeV isMater: " << IsMaster() < <<  78      << LowEnergyLimit() / MeV << " MeV - "
 98   }                                            <<  79      << HighEnergyLimit() / GeV << " GeV"
 99                                                <<  80      << G4endl;
100   if (IsMaster()) {                            <<  81   }
101     // Initialise element selector             <<  82 
102     InitialiseElementSelectors(particle, cuts) <<  83   // Initialise element selector
103                                                <<  84   
104     // Access to elements                      <<  85   InitialiseElementSelectors(particle, cuts);
105     const G4ElementTable* elemTable = G4Elemen <<  86 
106     std::size_t numElems = (*elemTable).size() <<  87   // Access to elements
107     for (std::size_t ie = 0; ie < numElems; ++ <<  88   
108       const G4Element* elem = (*elemTable)[ie] <<  89   char* path = getenv("G4LEDATA");
109       const G4int Z = std::min(maxZ, elem->Get <<  90 
110       if (data[Z] == nullptr) {                <<  91   G4ProductionCutsTable* theCoupleTable =
111         ReadData(Z);                           <<  92     G4ProductionCutsTable::GetProductionCutsTable();
112       }                                        <<  93   G4int numOfCouples = theCoupleTable->GetTableSize();
                                                   >>  94   
                                                   >>  95   for(G4int i=0; i<numOfCouples; ++i) 
                                                   >>  96   {
                                                   >>  97     const G4Material* material = 
                                                   >>  98       theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
                                                   >>  99     const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 100     G4int nelm = material->GetNumberOfElements();
                                                   >> 101     
                                                   >> 102     for (G4int j=0; j<nelm; ++j) 
                                                   >> 103     {
                                                   >> 104         
                                                   >> 105       G4int Z = (G4int)(*theElementVector)[j]->GetZ();
                                                   >> 106       if(Z < 1)          { Z = 1; }
                                                   >> 107       else if(Z > maxZ)  { Z = maxZ; }
                                                   >> 108       if(!data[Z]) { ReadData(Z, path); }
113     }                                             109     }
114   }                                               110   }
115   if (isInitialised) {                         << 111   //
116     return;                                    << 112   
117   }                                            << 113   if(isInitialised) { return; }
118   fParticleChange = GetParticleChangeForGamma(    114   fParticleChange = GetParticleChangeForGamma();
119   isInitialised = true;                           115   isInitialised = true;
120 }                                                 116 }
121                                                   117 
122 //....oooOO0OOooo........oooOO0OOooo........oo    118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123                                                   119 
124 const G4String& G4LivermoreGammaConversionMode << 120 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path)
125 {                                                 121 {
126   // no check in this method - environment var << 122   if (verboseLevel > 1) 
127   if (gDataDirectory.empty()) {                << 123   {
128     auto param = G4EmParameters::Instance();   << 124     G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel" 
129     std::ostringstream ost;                    << 125      << G4endl;
130     if (param->LivermoreDataDir() == "livermor << 126   }
131       ost << param->GetDirLEDATA() << "/liverm << 127 
132       useSpline = true;                        << 128   if(data[Z]) { return; }
                                                   >> 129   
                                                   >> 130   const char* datadir = path;
                                                   >> 131 
                                                   >> 132   if(!datadir) 
                                                   >> 133   {
                                                   >> 134     datadir = getenv("G4LEDATA");
                                                   >> 135     if(!datadir) 
                                                   >> 136     {
                                                   >> 137       G4Exception("G4LivermoreGammaConversionModel::ReadData()",
                                                   >> 138       "em0006",FatalException,
                                                   >> 139       "Environment variable G4LEDATA not defined");
                                                   >> 140       return;
133     }                                             141     }
134     else {                                     << 
135       ost << param->GetDirLEDATA() << "/epics2 << 
136     }                                          << 
137     gDataDirectory = ost.str();                << 
138   }                                            << 
139   return gDataDirectory;                       << 
140 }                                              << 
141                                                << 
142 //....oooOO0OOooo........oooOO0OOooo........oo << 
143                                                << 
144 void G4LivermoreGammaConversionModel::ReadData << 
145 {                                              << 
146   if (verboseLevel > 1) {                      << 
147     G4cout << "Calling ReadData() of G4Livermo << 
148   }                                            << 
149                                                << 
150   if (data[Z] != nullptr) {                    << 
151     return;                                    << 
152   }                                               142   }
153                                                   143 
                                                   >> 144   //
                                                   >> 145   
                                                   >> 146   data[Z] = new G4LPhysicsFreeVector();
                                                   >> 147   
                                                   >> 148   // Activation of spline interpolation
                                                   >> 149   data[Z] ->SetSpline(true);
                                                   >> 150   //
                                                   >> 151   
154   std::ostringstream ost;                         152   std::ostringstream ost;
155   ost << FindDirectoryPath() << "pp-cs-" << Z  << 153   ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
156                                                << 
157   data[Z] = new G4PhysicsFreeVector(useSpline) << 
158                                                << 
159   std::ifstream fin(ost.str().c_str());           154   std::ifstream fin(ost.str().c_str());
160                                                << 155   
161   if (!fin.is_open()) {                        << 156   if( !fin.is_open()) 
                                                   >> 157   {
162     G4ExceptionDescription ed;                    158     G4ExceptionDescription ed;
163     ed << "G4LivermoreGammaConversionModel dat << 159     ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
164        << G4endl;                              << 160        << "> is not opened!" << G4endl;
165     G4Exception("G4LivermoreGammaConversionMod << 161     G4Exception("G4LivermoreGammaConversionModel::ReadData()",
166                 "G4LEDATA version should be G4 << 162     "em0003",FatalException,
                                                   >> 163     ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
167     return;                                       164     return;
168   }                                            << 165   } 
169   else {                                       << 166   
170     if (verboseLevel > 1) {                    << 167   else 
171       G4cout << "File " << ost.str() << " is o << 168   {
172     }                                          << 169     
173                                                << 170     if(verboseLevel > 3) { G4cout << "File " << ost.str() 
                                                   >> 171        << " is opened by G4LivermoreGammaConversionModel" << G4endl;}
                                                   >> 172     
174     data[Z]->Retrieve(fin, true);                 173     data[Z]->Retrieve(fin, true);
175   }                                            << 174   } 
176   // Activation of spline interpolation        << 175   
177   if (useSpline) data[Z]->FillSecondDerivative << 176   
178 }                                                 177 }
179                                                   178 
180 //....oooOO0OOooo........oooOO0OOooo........oo    179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
181                                                   180 
182 G4double                                       << 181 G4double 
183 G4LivermoreGammaConversionModel::ComputeCrossS << 182 G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
184                                                << 183                   G4double GammaEnergy,
185                                                << 184                   G4double Z, G4double,
                                                   >> 185                   G4double, G4double)
186 {                                                 186 {
187   if (verboseLevel > 1) {                      << 187   if (verboseLevel > 1) 
188     G4cout << "G4LivermoreGammaConversionModel << 188   {
                                                   >> 189     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel" 
                                                   >> 190      << G4endl;
189   }                                               191   }
190                                                   192 
191   if (GammaEnergy < lowEnergyLimit) {          << 193   if (GammaEnergy < lowEnergyLimit) { return 0.0; } 
192     return 0.0;                                << 
193   }                                            << 
194                                                   194 
195   G4double xs = 0.0;                              195   G4double xs = 0.0;
                                                   >> 196   
                                                   >> 197   G4int intZ=G4int(Z);
196                                                   198 
197   G4int intZ = std::max(1, std::min(G4lrint(Z) << 199   if(intZ < 1 || intZ > maxZ) { return xs; }
198                                                   200 
199   G4PhysicsFreeVector* pv = data[intZ];        << 201   G4LPhysicsFreeVector* pv = data[intZ];
200                                                   202 
201   // if element was not initialised            << 203   // element was not initialised
202   // do initialisation safely for MT mode      << 204   if(!pv) 
203   if (pv == nullptr) {                         << 205   {
204     InitialiseForElement(particle, intZ);      << 206     char* path = getenv("G4LEDATA");
                                                   >> 207     ReadData(intZ, path);
205     pv = data[intZ];                              208     pv = data[intZ];
206     if (pv == nullptr) {                       << 209     if(!pv) { return xs; }
207       return xs;                               << 
208     }                                          << 
209   }                                               210   }
210   // x-section is taken from the table            211   // x-section is taken from the table
211   xs = pv->Value(GammaEnergy);                 << 212   xs = pv->Value(GammaEnergy); 
212                                                   213 
213   if (verboseLevel > 0) {                      << 214   if(verboseLevel > 0)
214     G4cout << "*** Gamma conversion xs for Z=" << 215   {
215            << "  cs=" << xs / millibarn << " m << 216     G4int n = pv->GetVectorLength() - 1;
                                                   >> 217     G4cout  <<  "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << GammaEnergy/MeV << G4endl;
                                                   >> 218     G4cout  <<  "  cs (Geant4 internal unit)=" << xs << G4endl;
                                                   >> 219     G4cout  <<  "    -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
                                                   >> 220     G4cout  <<  "    -> last  cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
                                                   >> 221     G4cout  <<  "*********************************************************" << G4endl;
216   }                                               222   }
                                                   >> 223 
217   return xs;                                      224   return xs;
                                                   >> 225 
218 }                                                 226 }
219                                                   227 
220 //....oooOO0OOooo........oooOO0OOooo........oo << 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221                                                   229 
222 void G4LivermoreGammaConversionModel::Initiali << 230 void G4LivermoreGammaConversionModel::SampleSecondaries(
                                                   >> 231                                  std::vector<G4DynamicParticle*>* fvect,
                                                   >> 232          const G4MaterialCutsCouple* couple,
                                                   >> 233          const G4DynamicParticle* aDynamicGamma,
                                                   >> 234          G4double, G4double)
223 {                                                 235 {
224   G4AutoLock l(&LivermoreGammaConversionModelM << 236 
225   if (data[Z] == nullptr) {                    << 237 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
226     ReadData(Z);                               << 238 // cross sections with Coulomb correction. A modified version of the random
227   }                                            << 239 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
228   l.unlock();                                  << 240 
                                                   >> 241 // Note 1 : Effects due to the breakdown of the Born approximation at low
                                                   >> 242 // energy are ignored.
                                                   >> 243 // Note 2 : The differential cross section implicitly takes account of
                                                   >> 244 // pair creation in both nuclear and atomic electron fields. However triplet
                                                   >> 245 // prodution is not generated.
                                                   >> 246 
                                                   >> 247   if (verboseLevel > 1)
                                                   >> 248     G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl;
                                                   >> 249 
                                                   >> 250   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
                                                   >> 251   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
                                                   >> 252 
                                                   >> 253   G4double epsilon ;
                                                   >> 254   G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
                                                   >> 255 
                                                   >> 256   // Do it fast if photon energy < 2. MeV
                                                   >> 257   if (photonEnergy < smallEnergy )
                                                   >> 258   {
                                                   >> 259     epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
                                                   >> 260   }
                                                   >> 261   else
                                                   >> 262   {
                                                   >> 263     // Select randomly one element in the current material
                                                   >> 264 
                                                   >> 265     const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
                                                   >> 266     const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
                                                   >> 267 
                                                   >> 268     if (element == 0)
                                                   >> 269       {
                                                   >> 270   G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0" 
                                                   >> 271          << G4endl;
                                                   >> 272   return;
                                                   >> 273       }
                                                   >> 274     G4IonisParamElm* ionisation = element->GetIonisation();
                                                   >> 275     if (ionisation == 0)
                                                   >> 276       {
                                                   >> 277   G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0" 
                                                   >> 278          << G4endl;
                                                   >> 279   return;
                                                   >> 280       }
                                                   >> 281 
                                                   >> 282     // Extract Coulomb factor for this Element
                                                   >> 283     G4double fZ = 8. * (ionisation->GetlogZ3());
                                                   >> 284     if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
                                                   >> 285 
                                                   >> 286     // Limits of the screening variable
                                                   >> 287     G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
                                                   >> 288     G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
                                                   >> 289     G4double screenMin = std::min(4.*screenFactor,screenMax) ;
                                                   >> 290 
                                                   >> 291     // Limits of the energy sampling
                                                   >> 292     G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
                                                   >> 293     G4double epsilonMin = std::max(epsilon0Local,epsilon1);
                                                   >> 294     G4double epsilonRange = 0.5 - epsilonMin ;
                                                   >> 295 
                                                   >> 296     // Sample the energy rate of the created electron (or positron)
                                                   >> 297     G4double screen;
                                                   >> 298     G4double gReject ;
                                                   >> 299 
                                                   >> 300     G4double f10 = ScreenFunction1(screenMin) - fZ;
                                                   >> 301     G4double f20 = ScreenFunction2(screenMin) - fZ;
                                                   >> 302     G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
                                                   >> 303     G4double normF2 = std::max(1.5 * f20,0.);
                                                   >> 304 
                                                   >> 305     do 
                                                   >> 306       {
                                                   >> 307   if (normF1 / (normF1 + normF2) > G4UniformRand() )
                                                   >> 308     {
                                                   >> 309       epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
                                                   >> 310       screen = screenFactor / (epsilon * (1. - epsilon));
                                                   >> 311       gReject = (ScreenFunction1(screen) - fZ) / f10 ;
                                                   >> 312     }
                                                   >> 313   else
                                                   >> 314     {
                                                   >> 315       epsilon = epsilonMin + epsilonRange * G4UniformRand();
                                                   >> 316       screen = screenFactor / (epsilon * (1 - epsilon));
                                                   >> 317       gReject = (ScreenFunction2(screen) - fZ) / f20 ;
                                                   >> 318     }
                                                   >> 319       } while ( gReject < G4UniformRand() );
                                                   >> 320     
                                                   >> 321   }   //  End of epsilon sampling
                                                   >> 322 
                                                   >> 323   // Fix charges randomly
                                                   >> 324 
                                                   >> 325   G4double electronTotEnergy;
                                                   >> 326   G4double positronTotEnergy;
                                                   >> 327 
                                                   >> 328   if (G4UniformRand() > 0.5)
                                                   >> 329     {
                                                   >> 330       electronTotEnergy = (1. - epsilon) * photonEnergy;
                                                   >> 331       positronTotEnergy = epsilon * photonEnergy;
                                                   >> 332     }
                                                   >> 333   else
                                                   >> 334     {
                                                   >> 335       positronTotEnergy = (1. - epsilon) * photonEnergy;
                                                   >> 336       electronTotEnergy = epsilon * photonEnergy;
                                                   >> 337     }
                                                   >> 338 
                                                   >> 339   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
                                                   >> 340   // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
                                                   >> 341   // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
                                                   >> 342 
                                                   >> 343   G4double u;
                                                   >> 344   const G4double a1 = 0.625;
                                                   >> 345   G4double a2 = 3. * a1;
                                                   >> 346   //  G4double d = 27. ;
                                                   >> 347 
                                                   >> 348   //  if (9. / (9. + d) > G4UniformRand())
                                                   >> 349   if (0.25 > G4UniformRand())
                                                   >> 350     {
                                                   >> 351       u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
                                                   >> 352     }
                                                   >> 353   else
                                                   >> 354     {
                                                   >> 355       u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
                                                   >> 356     }
                                                   >> 357 
                                                   >> 358   G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
                                                   >> 359   G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
                                                   >> 360   G4double phi  = twopi * G4UniformRand();
                                                   >> 361 
                                                   >> 362   G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
                                                   >> 363   G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
                                                   >> 364   
                                                   >> 365   
                                                   >> 366   // Kinematics of the created pair:
                                                   >> 367   // the electron and positron are assumed to have a symetric angular 
                                                   >> 368   // distribution with respect to the Z axis along the parent photon
                                                   >> 369   
                                                   >> 370   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
                                                   >> 371   
                                                   >> 372   G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
                                                   >> 373   electronDirection.rotateUz(photonDirection);
                                                   >> 374       
                                                   >> 375   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
                                                   >> 376               electronDirection,
                                                   >> 377               electronKineEnergy);
                                                   >> 378 
                                                   >> 379   // The e+ is always created 
                                                   >> 380   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
                                                   >> 381 
                                                   >> 382   G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
                                                   >> 383   positronDirection.rotateUz(photonDirection);   
                                                   >> 384   
                                                   >> 385   // Create G4DynamicParticle object for the particle2 
                                                   >> 386   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
                                                   >> 387                    positronDirection, 
                                                   >> 388                    positronKineEnergy);
                                                   >> 389   // Fill output vector
                                                   >> 390   fvect->push_back(particle1);
                                                   >> 391   fvect->push_back(particle2);
                                                   >> 392 
                                                   >> 393   // kill incident photon
                                                   >> 394   fParticleChange->SetProposedKineticEnergy(0.);
                                                   >> 395   fParticleChange->ProposeTrackStatus(fStopAndKill);   
                                                   >> 396 
229 }                                                 397 }
230                                                   398 
231 //....oooOO0OOooo........oooOO0OOooo........oo << 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 400 
                                                   >> 401 G4double 
                                                   >> 402 G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable)
                                                   >> 403 {
                                                   >> 404   // Compute the value of the screening function 3*phi1 - phi2
                                                   >> 405 
                                                   >> 406   G4double value;
                                                   >> 407   
                                                   >> 408   if (screenVariable > 1.)
                                                   >> 409     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
                                                   >> 410   else
                                                   >> 411     value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
                                                   >> 412   
                                                   >> 413   return value;
                                                   >> 414 } 
                                                   >> 415 
                                                   >> 416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 417 
                                                   >> 418 G4double 
                                                   >> 419 G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable)
                                                   >> 420 {
                                                   >> 421   // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
                                                   >> 422   
                                                   >> 423   G4double value;
                                                   >> 424   
                                                   >> 425   if (screenVariable > 1.)
                                                   >> 426     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
                                                   >> 427   else
                                                   >> 428     value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
                                                   >> 429   
                                                   >> 430   return value;
                                                   >> 431 } 
                                                   >> 432 
232                                                   433