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