| 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: G4eBremsstrahlung.cc,v 1.24 2001/11/09 13:59:46 maire Exp $
>> 25 // GEANT4 tag $Name: geant4-04-00 $
28 // 26 //
29 // GEANT4 Class file <<
30 // <<
31 // <<
32 // File name: G4eBremsstrahlung <<
33 // <<
34 // Author: Michel Maire <<
35 // <<
36 // Creation date: 26.06.1996 <<
37 // <<
38 // Modified by Michel Maire, Vladimir Ivanchen <<
39 // <<
40 // ------------------------------------------- <<
41 // 27 //
>> 28 // ------------ G4eBremsstrahlung physics process --------
>> 29 // by Michel Maire, 24 July 1996
>> 30 //
>> 31 // 26-09-96 extension of the total crosssection above 100 GeV, M.Maire
>> 32 // 1-10-96 new type G4OrderedTable; ComputePartialSumSigma(), M.Maire
>> 33 // 16-10-96 DoIt() call to the non static GetEnergyCuts(), L.Urban
>> 34 // 13-12-96 Sign corrected in grejmax and greject
>> 35 // error definition of screenvar, L.Urban
>> 36 // 20-03-97 new energy loss+ionisation+brems scheme, L.Urban
>> 37 // 07-04-98 remove 'tracking cut' of the diffracted particle, MMa
>> 38 // 13-08-98 new methods SetBining() PrintInfo()
>> 39 // 03-03-99 Bug fixed in LPM effect, L.Urban
>> 40 // 10-02-00 modifications , new e.m. structure, L.Urban
>> 41 // 07-08-00 new cross section/en.loss parametrisation, LPM flag , L.Urban
>> 42 // 21-09-00 corrections in the LPM implementation, L.Urban
>> 43 // 28-05-01 V.Ivanchenko minor changes to provide ANSI -wall compilation
>> 44 // 09-08-01 new methods Store/Retrieve PhysicsTable (mma)
>> 45 // 17-09-01 migration of Materials to pure STL (mma)
>> 46 // 21-09-01 completion of RetrievePhysicsTable() (mma)
>> 47 // 29-10-01 all static functions no more inlined (mma)
>> 48 // 08-11-01 particleMass becomes a local variable
>> 49 // --------------------------------------------------------------
>> 50
42 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 << 53
45 #include "G4eBremsstrahlung.hh" 54 #include "G4eBremsstrahlung.hh"
46 #include "G4SystemOfUnits.hh" << 55 #include "G4EnergyLossTables.hh"
47 #include "G4Gamma.hh" << 56 #include "G4ios.hh"
48 #include "G4SeltzerBergerModel.hh" <<
49 #include "G4eBremsstrahlungRelModel.hh" <<
50 #include "G4UnitsTable.hh" 57 #include "G4UnitsTable.hh"
51 #include "G4LossTableManager.hh" <<
52 58
53 #include "G4ProductionCutsTable.hh" << 59 G4double G4eBremsstrahlung::LowerBoundLambda = 1.*keV;
54 #include "G4MaterialCutsCouple.hh" << 60 G4double G4eBremsstrahlung::UpperBoundLambda = 100.*TeV;
55 #include "G4EmParameters.hh" << 61 G4int G4eBremsstrahlung::NbinLambda = 100;
>> 62 G4double G4eBremsstrahlung::probsup = 1.00;
>> 63 G4bool G4eBremsstrahlung::LPMflag = true;
56 64
57 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 66
>> 67 // constructor
>> 68
>> 69 G4eBremsstrahlung::G4eBremsstrahlung(const G4String& processName)
>> 70 : G4VeEnergyLoss(processName), // initialization
>> 71 theMeanFreePathTable(NULL)
>> 72 { // MinThreshold = 10*keV;
>> 73 MinThreshold = 1*keV; }
58 74
59 using namespace std; << 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 76
>> 77 // destructor
60 78
61 G4eBremsstrahlung::G4eBremsstrahlung(const G4S << 79 G4eBremsstrahlung::~G4eBremsstrahlung()
62 G4VEnergyLossProcess(name) << 80 {
>> 81 if (theMeanFreePathTable) {
>> 82 theMeanFreePathTable->clearAndDestroy();
>> 83 delete theMeanFreePathTable;
>> 84 }
>> 85
>> 86 if (&PartialSumSigma) {
>> 87 PartialSumSigma.clearAndDestroy();
>> 88 }
>> 89 }
>> 90
>> 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 92
>> 93 void G4eBremsstrahlung::SetLowerBoundLambda(G4double val)
>> 94 {LowerBoundLambda = val;}
>> 95
>> 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 97
>> 98 void G4eBremsstrahlung::SetUpperBoundLambda(G4double val)
>> 99 {UpperBoundLambda = val;}
>> 100
>> 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 102
>> 103 void G4eBremsstrahlung::SetNbinLambda(G4int n)
>> 104 {NbinLambda = n;}
>> 105
>> 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 107
>> 108 G4double G4eBremsstrahlung::GetLowerBoundLambda()
>> 109 {return LowerBoundLambda;}
>> 110
>> 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 112
>> 113 G4double G4eBremsstrahlung::GetUpperBoundLambda()
>> 114 {return UpperBoundLambda;}
>> 115
>> 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 117
>> 118 G4int G4eBremsstrahlung::GetNbinLambda()
>> 119 {return NbinLambda;}
>> 120
>> 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 122
>> 123 void G4eBremsstrahlung::SetLPMflag(G4bool val)
>> 124 {LPMflag = val;}
>> 125
>> 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 127
>> 128 G4bool G4eBremsstrahlung::GetLPMflag()
>> 129 {return LPMflag;}
>> 130
>> 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 132
>> 133 void G4eBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
>> 134 // just call BuildLossTable+BuildLambdaTable
>> 135 {
>> 136 // get bining from EnergyLoss
>> 137 LowestKineticEnergy = GetLowerBoundEloss();
>> 138 HighestKineticEnergy = GetUpperBoundEloss();
>> 139 TotBin = GetNbinEloss();
>> 140
>> 141 BuildLossTable(aParticleType);
>> 142
>> 143 if (&aParticleType==G4Electron::Electron())
>> 144 {
>> 145 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 146 CounterOfElectronProcess++;
>> 147 }
>> 148 else
>> 149 {
>> 150 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 151 CounterOfPositronProcess++;
>> 152 }
>> 153
>> 154 BuildLambdaTable(aParticleType);
>> 155
>> 156 BuildDEDXTable (aParticleType);
>> 157
>> 158 if (&aParticleType==G4Electron::Electron()) PrintInfoDefinition();
>> 159 }
>> 160
>> 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 162
>> 163 void G4eBremsstrahlung::BuildLossTable(const G4ParticleDefinition& aParticleType)
>> 164 // Build table for energy loss due to soft brems
>> 165 // tables are built for *MATERIALS*
>> 166 {
>> 167 G4double KineticEnergy,TotalEnergy,bremloss,Z,x,
>> 168 losslim,loss,rate,natom,Cut;
>> 169
>> 170 const G4double MinKinEnergy = 1.*keV;
>> 171 const G4double Thigh = 100.*GeV;
>> 172 const G4double Cuthigh = 50.*GeV;
>> 173 const G4double Factorhigh = 36./(1450.*GeV);
>> 174 const G4double coef1 = -0.5, coef2 = 2./9.;
>> 175
>> 176 G4double particleMass = aParticleType.GetPDGMass() ;
>> 177 G4double* GammaCutInKineticEnergy = G4Gamma::Gamma()->GetEnergyCuts();
>> 178
>> 179 // create table
>> 180
>> 181 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 182 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
>> 183
>> 184 if (theLossTable) { theLossTable->clearAndDestroy();
>> 185 delete theLossTable;
>> 186 }
>> 187
>> 188 theLossTable = new G4PhysicsTable(numOfMaterials);
>> 189
>> 190 // loop for materials
>> 191
>> 192 for (G4int J=0; J<numOfMaterials; J++)
>> 193 {
>> 194 // create physics vector and fill it
>> 195
>> 196 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 197 LowestKineticEnergy,HighestKineticEnergy,TotBin);
>> 198
>> 199 // get elements in the material
>> 200 const G4Material* material = (*theMaterialTable)[J];
>> 201
>> 202 const G4ElementVector* theElementVector = material->GetElementVector();
>> 203 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
>> 204 const G4int NumberOfElements = material->GetNumberOfElements();
>> 205
>> 206 // loop for the kinetic energy values
>> 207 for (G4int i=0; i<TotBin; i++)
>> 208 {
>> 209 KineticEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 210 TotalEnergy = KineticEnergy+particleMass ;
>> 211 Cut = GammaCutInKineticEnergy[J] ;
>> 212 if (Cut < MinThreshold) Cut = MinThreshold;
>> 213 if (Cut > KineticEnergy) Cut = KineticEnergy;
>> 214
>> 215 bremloss = 0.;
>> 216
>> 217 if (KineticEnergy>MinKinEnergy)
>> 218 {
>> 219 // loop for elements in the material
>> 220 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 221 {
>> 222 Z=(*theElementVector)[iel]->GetZ();
>> 223 natom = theAtomicNumDensityVector[iel] ;
>> 224 if (KineticEnergy <= Thigh)
>> 225 {
>> 226 //loss for MinKinEnergy<KineticEnergy<=100 GeV
>> 227 x=log(TotalEnergy/particleMass);
>> 228 loss = ComputeBremLoss(Z,natom,KineticEnergy,Cut,x) ;
>> 229 if (&aParticleType==G4Positron::Positron())
>> 230 loss *= ComputePositronCorrFactorLoss(Z,KineticEnergy,Cut) ;
>> 231 }
>> 232 else
>> 233 {
>> 234 // extrapolation for KineticEnergy>100 GeV
>> 235 x=log(Thigh/particleMass) ;
>> 236 if (Cut<Thigh)
>> 237 {
>> 238 losslim = ComputeBremLoss(Z,natom,Thigh,Cut,x) ;
>> 239 if (&aParticleType==G4Positron::Positron())
>> 240 loss *= ComputePositronCorrFactorLoss(Z,Thigh,Cut) ;
>> 241 rate = Cut/TotalEnergy ;
>> 242 loss = losslim*(1.+coef1*rate+coef2*rate*rate) ;
>> 243 rate = Cut/Thigh ;
>> 244 loss /= (1.+coef1*rate+coef2*rate*rate) ;
>> 245 }
>> 246 else
>> 247 {
>> 248 losslim = ComputeBremLoss(Z,natom,Thigh,Cuthigh,x) ;
>> 249 if (&aParticleType==G4Positron::Positron())
>> 250 loss *= ComputePositronCorrFactorLoss(Z,Thigh,Cuthigh) ;
>> 251 rate = Cut/TotalEnergy ;
>> 252 loss = losslim*(1.+coef1*rate+coef2*rate*rate) ;
>> 253 loss *= Factorhigh*Cut ;
>> 254 }
>> 255
>> 256 }
>> 257 bremloss += natom*loss;
>> 258 }
>> 259
>> 260 }
>> 261
>> 262 static const G4double MigdalConstant = classic_electr_radius
>> 263 *electron_Compton_length
>> 264 *electron_Compton_length/pi;
>> 265 G4double TotalEnergy = KineticEnergy+electron_mass_c2 ;
>> 266 G4double kp2 = MigdalConstant*TotalEnergy*TotalEnergy*
>> 267 (material->GetElectronDensity()) ;
>> 268
>> 269 // now compute the correction due to the supression(s)
>> 270 G4double kmin = 1.*eV ;
>> 271 G4double kmax = Cut ;
>> 272
>> 273 if(kmax > kmin)
>> 274 {
>> 275
>> 276 G4double floss = 0. ;
>> 277 G4int nmax = 100 ;
>> 278 G4int nn ;
>> 279 G4double vmin=log(kmin);
>> 280 G4double vmax=log(kmax) ;
>> 281 nn = int(nmax*(vmax-vmin)/(log(HighestKineticEnergy)-vmin)) ;
>> 282 G4double u,fac,c,v,dv ;
>> 283 dv = (vmax-vmin)/nn ;
>> 284 v = vmin-dv ;
>> 285 if(nn > 0)
>> 286 {
>> 287 for(G4int n=0; n<=nn; n++)
>> 288 {
>> 289 v += dv; u = exp(v);
>> 290 fac = u*SupressionFunction(material,KineticEnergy,u);
>> 291 probsup = 1.;
>> 292 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 293 if ((n==0)||(n==nn)) c=0.5;
>> 294 else c=1. ;
>> 295 fac *= c ;
>> 296 floss += fac ;
>> 297 }
>> 298 floss *=dv/(kmax-kmin);
>> 299 }
>> 300 else floss = 1.;
>> 301 if(floss > 1.) floss = 1.;
>> 302 // correct the loss
>> 303 bremloss *= floss;
>> 304 }
>> 305
>> 306 if(bremloss < 0.) bremloss = 0.;
>> 307 aVector->PutValue(i,bremloss);
>> 308 }
>> 309
>> 310 theLossTable->insert(aVector);
>> 311 }
>> 312
>> 313 }
>> 314
>> 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 316
>> 317 G4double G4eBremsstrahlung::ComputeBremLoss(G4double Z,G4double natom,
>> 318 G4double T,G4double Cut,G4double x)
>> 319
>> 320 // compute loss due to soft brems
63 { 321 {
64 SetProcessSubType(fBremsstrahlung); << 322 static const G4double beta=1.00,ksi=2.00 ;
65 SetSecondaryParticle(G4Gamma::Gamma()); << 323 static const G4double clossh = 0.254 , closslow = 1./3. , alosslow = 1. ;
66 SetIonisation(false); << 324 static const G4double Tlim= 10.*MeV ;
67 SetCrossSectionType(fEmTwoPeaks); << 325
>> 326 static const G4double xlim = 1.2 ;
>> 327 static const G4int NZ = 8 ;
>> 328 static const G4int Nloss = 11 ;
>> 329 static const G4double ZZ[NZ] =
>> 330 {2.,4.,6.,14.,26.,50.,82.,92.};
>> 331 static const G4double coefloss[NZ][Nloss] = {
>> 332 // Z=2
>> 333 { 0.98916, 0.47564, -0.2505, -0.45186, 0.14462,
>> 334 0.21307, -0.013738, -0.045689, -0.0042914, 0.0034429,
>> 335 0.00064189},
>> 336
>> 337 // Z=4
>> 338 { 1.0626, 0.37662, -0.23646, -0.45188, 0.14295,
>> 339 0.22906, -0.011041, -0.051398, -0.0055123, 0.0039919,
>> 340 0.00078003},
>> 341 // Z=6
>> 342 { 1.0954, 0.315, -0.24011, -0.43849, 0.15017,
>> 343 0.23001, -0.012846, -0.052555, -0.0055114, 0.0041283,
>> 344 0.00080318},
>> 345
>> 346 // Z=14
>> 347 { 1.1649, 0.18976, -0.24972, -0.30124, 0.1555,
>> 348 0.13565, -0.024765, -0.027047, -0.00059821, 0.0019373,
>> 349 0.00027647},
>> 350
>> 351 // Z=26
>> 352 { 1.2261, 0.14272, -0.25672, -0.28407, 0.13874,
>> 353 0.13586, -0.020562, -0.026722, -0.00089557, 0.0018665,
>> 354 0.00026981},
>> 355
>> 356 // Z=50
>> 357 { 1.3147, 0.020049, -0.35543, -0.13927, 0.17666,
>> 358 0.073746, -0.036076, -0.013407, 0.0025727, 0.00084005,
>> 359 -1.4082e-05},
>> 360
>> 361 // Z=82
>> 362 { 1.3986, -0.10586, -0.49187, -0.0048846, 0.23621,
>> 363 0.031652, -0.052938, -0.0076639, 0.0048181, 0.00056486,
>> 364 -0.00011995},
>> 365
>> 366 // Z=92
>> 367 { 1.4217, -0.116, -0.55497, -0.044075, 0.27506,
>> 368 0.081364, -0.058143, -0.023402, 0.0031322, 0.0020201,
>> 369 0.00017519}
>> 370
>> 371 } ;
>> 372 static G4double aaa=0.414 ;
>> 373 static G4double bbb=0.345 ;
>> 374 static G4double ccc=0.460 ;
>> 375
>> 376 G4int iz = 0;
>> 377 G4double delz = 1.e6;
>> 378 for (G4int ii=0; ii<NZ; ii++)
>> 379 {
>> 380 if(abs(Z-ZZ[ii]) < delz) { iz = ii; delz = abs(Z-ZZ[ii]);}
>> 381 }
>> 382
>> 383 G4double xx = log10(T);
>> 384 G4double fl = 1.;
>> 385
>> 386 if (xx <= xlim)
>> 387 {
>> 388 fl = coefloss[iz][Nloss-1];
>> 389 for (G4int j=Nloss-2; j>=0; j--) fl = fl*xx+coefloss[iz][j];
>> 390 if (fl < 0.) fl = 0.;
>> 391 }
>> 392
>> 393 G4double loss;
>> 394 G4double E = T+electron_mass_c2 ;
>> 395
>> 396 loss = Z*(Z+ksi)*E*E/(T+E)*exp(beta*log(Cut/T))*(2.-clossh*exp(log(Z)/4.));
>> 397 if (T <= Tlim) loss /= exp(closslow*log(Tlim/T));
>> 398 if( T <= Cut) loss *= exp(alosslow*log(T/Cut));
>> 399
>> 400 // correction ................................
>> 401 loss *= (aaa+bbb*T/Tlim)/(1.+ccc*T/Tlim);
>> 402 loss *= fl;
>> 403 loss /= Avogadro;
>> 404
>> 405 return loss;
68 } 406 }
69 407
70 //....oooOO0OOooo........oooOO0OOooo........oo 408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 409
72 G4eBremsstrahlung::~G4eBremsstrahlung() = defa << 410 G4double G4eBremsstrahlung::ComputePositronCorrFactorLoss(
>> 411 G4double Z,G4double KineticEnergy,G4double GammaCut)
>> 412
>> 413 //calculates the correction factor for the energy loss due to bremsstrahlung for positrons
>> 414 //the same correction is in the (discrete) bremsstrahlung
>> 415
>> 416 {
>> 417 static const G4double K = 132.9416*eV ;
>> 418 static const G4double a1=4.15e-1, a3=2.10e-3, a5=54.0e-5 ;
>> 419
>> 420 G4double x = log(KineticEnergy/(K*Z*Z)), x2 = x*x, x3 = x2*x;
>> 421 G4double eta = 0.5+atan(a1*x+a3*x3+a5*x3*x2)/pi;
>> 422 G4double e0 = GammaCut/KineticEnergy;
>> 423
>> 424 G4double factor(0.);
>> 425 if (e0!=1.0) { factor=log(1.-e0)/eta; factor=exp(factor);}
>> 426 factor = eta*(1.-factor)/e0;
>> 427
>> 428 return factor;
>> 429 }
>> 430
>> 431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 432
>> 433 void G4eBremsstrahlung::BuildLambdaTable(
>> 434 const G4ParticleDefinition& ParticleType)
>> 435
>> 436 // Build mean free path tables for the gamma emission by e- or e+.
>> 437 // tables are Build for MATERIALS.
>> 438 {
>> 439 G4double LowEdgeEnergy , Value;
>> 440 G4double FixedEnergy = (LowerBoundLambda + UpperBoundLambda)/2.;
>> 441
>> 442 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 443
>> 444 //create table
>> 445 if (theMeanFreePathTable) {theMeanFreePathTable->clearAndDestroy();
>> 446 delete theMeanFreePathTable;
>> 447 }
>> 448 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 449
>> 450 PartialSumSigma.clearAndDestroy();
>> 451 PartialSumSigma.resize(G4Material::GetNumberOfMaterials());
>> 452
>> 453 G4PhysicsLogVector* ptrVector;
>> 454 for ( size_t J=0 ; J < G4Material::GetNumberOfMaterials(); J++ )
>> 455 {
>> 456 //create physics vector then fill it ....
>> 457 ptrVector = new G4PhysicsLogVector(LowerBoundLambda, UpperBoundLambda,
>> 458 NbinLambda ) ;
>> 459
>> 460 const G4Material* material= (*theMaterialTable)[J];
>> 461
>> 462 for ( G4int i = 0 ; i < NbinLambda ; i++ )
>> 463 {
>> 464 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 465 Value = ComputeMeanFreePath( &ParticleType, LowEdgeEnergy,
>> 466 material );
>> 467 ptrVector->PutValue( i , Value ) ;
>> 468 }
>> 469
>> 470 theMeanFreePathTable->insertAt( J , ptrVector );
>> 471
>> 472 // Compute the PartialSumSigma table at a given fixed energy
>> 473 ComputePartialSumSigma( &ParticleType, FixedEnergy, material) ;
>> 474 }
>> 475
>> 476 }
73 477
74 //....oooOO0OOooo........oooOO0OOooo........oo 478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 479
76 G4bool G4eBremsstrahlung::IsApplicable(const G << 480 G4double G4eBremsstrahlung::ComputeMeanFreePath(
>> 481 const G4ParticleDefinition* ParticleType,
>> 482 G4double KineticEnergy,
>> 483 const G4Material* aMaterial)
77 { 484 {
78 return (&p == G4Electron::Electron() || &p = << 485 const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
>> 486 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 487 G4double GammaEnergyCut = G4Gamma::Gamma()->GetEnergyThreshold(aMaterial);
>> 488 if (GammaEnergyCut < MinThreshold) GammaEnergyCut = MinThreshold;
>> 489
>> 490 G4double SIGMA = 0;
>> 491
>> 492 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
>> 493 {
>> 494 SIGMA += theAtomNumDensityVector[i] *
>> 495 ComputeCrossSectionPerAtom( ParticleType, KineticEnergy,
>> 496 (*theElementVector)[i]->GetZ(),
>> 497 GammaEnergyCut );
>> 498 }
>> 499 // now compute the correction due to the supression(s)
>> 500
>> 501 G4double kmax = KineticEnergy ;
>> 502 G4double kmin = GammaEnergyCut ;
>> 503
>> 504 static const G4double MigdalConstant = classic_electr_radius
>> 505 *electron_Compton_length
>> 506 *electron_Compton_length/pi;
>> 507 G4double TotalEnergy = KineticEnergy+electron_mass_c2 ;
>> 508 G4double kp2 = MigdalConstant*TotalEnergy*TotalEnergy*
>> 509 (aMaterial->GetElectronDensity()) ;
>> 510
>> 511 if(kmax > kmin)
>> 512 {
>> 513 G4double fsig = 0.;
>> 514 G4int nmax = 100 ;
>> 515 G4int nn ;
>> 516 G4double vmin=log(kmin);
>> 517 G4double vmax=log(kmax) ;
>> 518 nn = int(nmax*(vmax-vmin)/(log(HighestKineticEnergy)-vmin));
>> 519 G4double u,fac,c,v,dv,y ;
>> 520 dv = (vmax-vmin)/nn ;
>> 521 v = vmin-dv ;
>> 522 if(nn > 0)
>> 523 {
>> 524 for(G4int n=0; n<=nn; n++)
>> 525 {
>> 526 v += dv; u = exp(v);
>> 527 fac = SupressionFunction(aMaterial,KineticEnergy,u);
>> 528 y = u/kmax;
>> 529 fac *= (4.-4.*y+3.*y*y)/3.;
>> 530 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 531 if ((n==0)||(n==nn)) c=0.5;
>> 532 else c=1. ;
>> 533 fac *= c;
>> 534 fsig += fac;
>> 535 }
>> 536 y = kmin/kmax ;
>> 537 fsig *=dv/(-4.*log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
>> 538 }
>> 539 else fsig = 1.;
>> 540 if (fsig > 1.) fsig = 1.;
>> 541 // correct the cross section
>> 542 SIGMA *= fsig;
>> 543 }
>> 544
>> 545 return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
79 } 546 }
80 547
81 //....oooOO0OOooo........oooOO0OOooo........oo 548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 549
83 void << 550 G4double G4eBremsstrahlung::ComputeCrossSectionPerAtom(
84 G4eBremsstrahlung::InitialiseEnergyLossProcess << 551 const G4ParticleDefinition* ParticleType,
85 const G4ParticleDefinition*) << 552 G4double KineticEnergy, G4double AtomicNumber,
>> 553 G4double GammaEnergyCut)
>> 554
>> 555 // Calculates the cross section per atom in GEANT4 internal units.
>> 556 //
>> 557
86 { 558 {
87 if(!isInitialised) { << 559 G4double CrossSection = 0.0 ;
88 G4EmParameters* param = G4EmParameters::In << 560 if ( KineticEnergy < 1*keV ) return CrossSection;
>> 561 if ( KineticEnergy <= GammaEnergyCut ) return CrossSection;
>> 562
>> 563 static const G4double ksi=2.0, alfa=1.00;
>> 564 static const G4double csigh = 0.127, csiglow = 0.25, asiglow = 0.020*MeV ;
>> 565 static const G4double Tlim = 10.*MeV ;
>> 566
>> 567 static const G4double xlim = 1.2 ;
>> 568 static const G4int NZ = 8 ;
>> 569 static const G4int Nsig = 11 ;
>> 570 static const G4double ZZ[NZ] =
>> 571 {2.,4.,6.,14.,26.,50.,82.,92.} ;
>> 572 static const G4double coefsig[NZ][Nsig] = {
>> 573 // Z=2
>> 574 { 0.4638, 0.37748, 0.32249, -0.060362, -0.065004,
>> 575 -0.033457, -0.004583, 0.011954, 0.0030404, -0.0010077,
>> 576 -0.00028131},
>> 577
>> 578 // Z=4
>> 579 { 0.50008, 0.33483, 0.34364, -0.086262, -0.055361,
>> 580 -0.028168, -0.0056172, 0.011129, 0.0027528, -0.00092265,
>> 581 -0.00024348},
>> 582
>> 583 // Z=6
>> 584 { 0.51587, 0.31095, 0.34996, -0.11623, -0.056167,
>> 585 -0.0087154, 0.00053943, 0.0054092, 0.00077685, -0.00039635,
>> 586 -6.7818e-05},
>> 587
>> 588 // Z=14
>> 589 { 0.55058, 0.25629, 0.35854, -0.080656, -0.054308,
>> 590 -0.049933, -0.00064246, 0.016597, 0.0021789, -0.001327,
>> 591 -0.00025983},
89 592
90 G4double emax = param->MaxKinEnergy(); << 593 // Z=26
91 G4VEmFluctuationModel* fm = nullptr; << 594 { 0.5791, 0.26152, 0.38953, -0.17104, -0.099172,
>> 595 0.024596, 0.023718, -0.0039205, -0.0036658, 0.00041749,
>> 596 0.00023408},
92 597
93 if (nullptr == EmModel(0)) { SetEmModel(ne << 598 // Z=50
94 G4double energyLimit = std::min(EmModel(0) << 599 { 0.62085, 0.27045, 0.39073, -0.37916, -0.18878,
95 EmModel(0)->SetHighEnergyLimit(energyLimit << 600 0.23905, 0.095028, -0.068744, -0.023809, 0.0062408,
96 EmModel(0)->SetSecondaryThreshold(param->B << 601 0.0020407},
97 AddEmModel(1, EmModel(0), fm); <<
98 602
99 if(emax > energyLimit) { << 603 // Z=82
100 if (nullptr == EmModel(1)) { << 604 { 0.66053, 0.24513, 0.35404, -0.47275, -0.22837,
101 SetEmModel(new G4eBremsstrahlungRelModel()); << 605 0.35647, 0.13203, -0.1049, -0.034851, 0.0095046,
102 } << 606 0.0030535},
103 EmModel(1)->SetLowEnergyLimit(energyLimi << 607
104 EmModel(1)->SetHighEnergyLimit(emax); << 608 // Z=92
105 EmModel(1)->SetSecondaryThreshold(param- << 609 { 0.67143, 0.23079, 0.32256, -0.46248, -0.20013,
106 AddEmModel(1, EmModel(1), fm); << 610 0.3506, 0.11779, -0.1024, -0.032013, 0.0092279,
>> 611 0.0028592}
>> 612
>> 613 } ;
>> 614
>> 615 G4int iz = 0 ;
>> 616 G4double delz = 1.e6 ;
>> 617 for (G4int ii=0; ii<NZ; ii++)
>> 618 {
>> 619 if(abs(AtomicNumber-ZZ[ii]) < delz)
>> 620 {
>> 621 iz = ii ;
>> 622 delz = abs(AtomicNumber-ZZ[ii]) ;
107 } 623 }
108 isInitialised = true; <<
109 } 624 }
>> 625
>> 626 G4double xx = log10(KineticEnergy) ;
>> 627 G4double fs = 1. ;
>> 628
>> 629 if(xx <= xlim)
>> 630 {
>> 631 fs = coefsig[iz][Nsig-1] ;
>> 632 for (G4int j=Nsig-2; j>=0; j--)
>> 633 {
>> 634 fs = fs*xx+coefsig[iz][j] ;
>> 635 }
>> 636 if(fs < 0.) fs = 0. ;
>> 637 }
>> 638
>> 639
>> 640 CrossSection = AtomicNumber*(AtomicNumber+ksi)*
>> 641 (1.-csigh*exp(log(AtomicNumber)/4.))*
>> 642 pow(log(KineticEnergy/GammaEnergyCut),alfa) ;
>> 643
>> 644 if(KineticEnergy <= Tlim)
>> 645 CrossSection *= exp(csiglow*log(Tlim/KineticEnergy))*
>> 646 (1.+asiglow/(sqrt(AtomicNumber)*KineticEnergy)) ;
>> 647
>> 648 if (ParticleType == G4Positron::Positron())
>> 649 CrossSection *= ComputePositronCorrFactorSigma(AtomicNumber, KineticEnergy,
>> 650 GammaEnergyCut);
>> 651 CrossSection *= fs ;
>> 652 CrossSection /= Avogadro ;
>> 653
>> 654 if (CrossSection < 0.) CrossSection = 0.;
>> 655 return CrossSection;
>> 656 }
>> 657
>> 658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 659
>> 660 G4double G4eBremsstrahlung::ComputePositronCorrFactorSigma( G4double AtomicNumber,
>> 661 G4double KineticEnergy, G4double GammaEnergyCut)
>> 662
>> 663 // Calculates the correction factor for the total cross section of the positron bremsstrahl.
>> 664 // Eta is the ratio of positron to electron energy loss by bremstrahlung.
>> 665 // A parametrized formula from L. Urban is used to estimate eta. It is a fit to the results
>> 666 // of L. Kim & al: Phys Rev. A33,3002 (1986)
>> 667
>> 668 {
>> 669 static const G4double K = 132.9416*eV;
>> 670 static const G4double a1 = 4.15e-1, a3 = 2.10e-3, a5 = 54.0e-5;
>> 671
>> 672 G4double x = log(KineticEnergy/(K*AtomicNumber*AtomicNumber));
>> 673 G4double eta = 0.5 + atan(a1*x + a3*x*x*x + a5*x*x*x*x*x)/pi ;
>> 674 G4double alfa = (1. - eta)/eta;
>> 675 return eta*pow((1. - GammaEnergyCut/KineticEnergy) , alfa);
>> 676 }
>> 677
>> 678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 679
>> 680 void G4eBremsstrahlung::ComputePartialSumSigma(const G4ParticleDefinition* ParticleType,
>> 681 G4double KineticEnergy,
>> 682 const G4Material* aMaterial)
>> 683
>> 684 // Build the table of cross section per element. The table is built for MATERIALS.
>> 685 // This table is used by DoIt to select randomly an element in the material.
>> 686 {
>> 687 G4int Imate = aMaterial->GetIndex();
>> 688 G4int NbOfElements = aMaterial->GetNumberOfElements();
>> 689 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 690 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 691 G4double GammaEnergyCut = G4Gamma::Gamma()->GetEnergyThreshold(aMaterial);
>> 692
>> 693
>> 694 PartialSumSigma[Imate] = new G4DataVector();
>> 695
>> 696 G4double SIGMA = 0. ;
>> 697
>> 698 for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ )
>> 699 {
>> 700 SIGMA += theAtomNumDensityVector[Ielem] *
>> 701 ComputeCrossSectionPerAtom( ParticleType, KineticEnergy,
>> 702 (*theElementVector)[Ielem]->GetZ(),
>> 703 GammaEnergyCut );
>> 704 PartialSumSigma[Imate]->push_back(SIGMA);
>> 705 }
>> 706 }
>> 707
>> 708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 709
>> 710 G4VParticleChange* G4eBremsstrahlung::PostStepDoIt(const G4Track& trackData,
>> 711 const G4Step& stepData)
>> 712 //
>> 713 // The emitted gamma energy is sampled using a parametrized formula from L. Urban.
>> 714 // This parametrization is derived from :
>> 715 // cross-section values of Seltzer and Berger for electron energies 1 keV - 10 GeV,
>> 716 // screened Bethe Heilter differential cross section above 10 GeV,
>> 717 // Migdal corrections in both case.
>> 718 // Seltzer & Berger: Nim B 12:95 (1985)
>> 719 // Nelson, Hirayama & Rogers: Technical report 265 SLAC (1985)
>> 720 // Migdal: Phys Rev 103:1811 (1956); Messel & Crawford: Pergamon Press (1970)
>> 721 //
>> 722 // A modified version of the random number techniques of Butcher & Messel is used
>> 723 // (Nuc Phys 20(1960),15).
>> 724 //
>> 725 // GEANT4 internal units.
>> 726 //
>> 727 {
>> 728 static const G4double
>> 729 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
>> 730 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
>> 731 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
>> 732
>> 733 static const G4double
>> 734 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
>> 735 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
>> 736 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
>> 737
>> 738 static const G4double
>> 739 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
>> 740 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
>> 741 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
>> 742
>> 743 static const G4double
>> 744 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
>> 745 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
>> 746 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
>> 747
>> 748 static const G4double MigdalConstant = classic_electr_radius
>> 749 *electron_Compton_length
>> 750 *electron_Compton_length/pi;
>> 751 // const G4double LPMconstant = fine_structure_const*electron_mass_c2*
>> 752 // electron_mass_c2/(8.*pi*hbarc) ;
>> 753 G4double GammaEnergy ;
>> 754 G4bool LPMOK = false ;
>> 755
>> 756 aParticleChange.Initialize(trackData);
>> 757 G4Material* aMaterial=trackData.GetMaterial() ;
>> 758
>> 759 //G4double LPMEnergy = LPMconstant*(aMaterial->GetRadlen()) ;
>> 760
>> 761 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
>> 762 G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();
>> 763
>> 764 G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
>> 765 G4ParticleMomentum ParticleDirection = aDynamicParticle->GetMomentumDirection();
>> 766
>> 767 // Gamma production cut in this material
>> 768 G4double GammaEnergyCut = G4Gamma::Gamma()->GetEnergyThreshold(aMaterial);
>> 769 if (GammaEnergyCut < MinThreshold) GammaEnergyCut = MinThreshold;
>> 770
>> 771 // check against insufficient energy
>> 772 if (KineticEnergy < GammaEnergyCut)
>> 773 {
>> 774 aParticleChange.SetMomentumChange( ParticleDirection );
>> 775 aParticleChange.SetEnergyChange( KineticEnergy );
>> 776 aParticleChange.SetLocalEnergyDeposit (0.);
>> 777 aParticleChange.SetNumberOfSecondaries(0);
>> 778 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 779 }
>> 780
>> 781 // select randomly one element constituing the material
>> 782 G4Element* anElement = SelectRandomAtom(aMaterial);
>> 783
>> 784 // Extract Z factors for this Element
>> 785 G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
>> 786 G4double FZ = lnZ* (4.- 0.55*lnZ);
>> 787 G4double ZZ = anElement->GetIonisation()->GetZZ3();
>> 788
>> 789 // limits of the energy sampling
>> 790 G4double TotalEnergy = KineticEnergy + electron_mass_c2;
>> 791 //G4double TotalEnergysquare = TotalEnergy*TotalEnergy ;
>> 792 //G4double LPMGammaEnergyLimit = TotalEnergysquare/LPMEnergy ;
>> 793 G4double xmin = GammaEnergyCut/KineticEnergy, epsilmin = GammaEnergyCut/TotalEnergy;
>> 794 G4double epsilmax = KineticEnergy/TotalEnergy;
>> 795
>> 796 // Migdal factor
>> 797 G4double
>> 798 MigdalFactor = (aMaterial->GetElectronDensity())*MigdalConstant
>> 799 /(epsilmax*epsilmax);
>> 800
>> 801 //
>> 802 G4double x, epsil, greject, migdal, grejmax;
>> 803 G4double U = log(KineticEnergy/electron_mass_c2), U2 = U*U;
>> 804
>> 805 //
>> 806 // sample the energy rate of the emitted gamma for electron kinetic energy > 1 MeV
>> 807 //
>> 808
>> 809 do {
>> 810 if (KineticEnergy > 1.*MeV)
>> 811 {
>> 812 // parameters
>> 813 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12),
>> 814 ah2 = ah20 + ZZ* (ah21 + ZZ* ah22),
>> 815 ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
>> 816
>> 817 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12),
>> 818 bh2 = bh20 + ZZ* (bh21 + ZZ* bh22),
>> 819 bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
>> 820
>> 821 G4double ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
>> 822 G4double bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
>> 823
>> 824 // limit of the screening variable
>> 825 G4double screenfac =
>> 826 136.*electron_mass_c2/((anElement->GetIonisation()->GetZ3())*TotalEnergy);
>> 827 G4double screenmin = screenfac*epsilmin/(1.-epsilmin);
>> 828
>> 829 // Compute the maximum of the rejection function
>> 830 G4double F1 = G4std::max(ScreenFunction1(screenmin) - FZ ,0.);
>> 831 G4double F2 = G4std::max(ScreenFunction2(screenmin) - FZ ,0.);
>> 832 grejmax = (F1 - epsilmin* (F1*ah - bh*epsilmin*F2))/(42.392 - FZ);
>> 833
>> 834 // sample the energy rate of the emitted Gamma
>> 835 G4double screenvar;
>> 836
>> 837
>> 838 do {
>> 839
>> 840 x = pow(xmin, G4UniformRand());
>> 841 epsil = x*KineticEnergy/TotalEnergy;
>> 842 screenvar = screenfac*epsil/(1-epsil);
>> 843 F1 = G4std::max(ScreenFunction1(screenvar) - FZ ,0.);
>> 844 F2 = G4std::max(ScreenFunction2(screenvar) - FZ ,0.);
>> 845 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
>> 846 greject = migdal*(F1 - epsil* (ah*F1 - bh*epsil*F2))/(42.392 - FZ);
>> 847 } while( greject < G4UniformRand()*grejmax );
>> 848
>> 849 }
>> 850
>> 851 else
>> 852 {
>> 853 // sample the energy rate of the emitted gamma for electron kinetic energy < 1 MeV
>> 854 //
>> 855 // parameters
>> 856 G4double al0 = al00 + ZZ* (al01 + ZZ* al02),
>> 857 al1 = al10 + ZZ* (al11 + ZZ* al12),
>> 858 al2 = al20 + ZZ* (al21 + ZZ* al22);
>> 859
>> 860 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02),
>> 861 bl1 = bl10 + ZZ* (bl11 + ZZ* bl12),
>> 862 bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
>> 863
>> 864 G4double al = al0 + al1*U + al2*U2;
>> 865 G4double bl = bl0 + bl1*U + bl2*U2;
>> 866
>> 867 // Compute the maximum of the rejection function
>> 868 grejmax = G4std::max(1. + xmin* (al + bl*xmin), 1.+al+bl);
>> 869 G4double xm = -al/(2.*bl);
>> 870 if ((xmin < xm)&&(xm < 1.)) grejmax = G4std::max(grejmax, 1.+ xm* (al + bl*xm));
>> 871
>> 872 // sample the energy rate of the emitted Gamma
>> 873
>> 874 do { x = pow(xmin, G4UniformRand());
>> 875 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
>> 876 greject = migdal*(1. + x* (al + bl*x));
>> 877 } while( greject < G4UniformRand()*grejmax );
>> 878 }
>> 879
>> 880 GammaEnergy = x*KineticEnergy;
>> 881
>> 882 if(LPMflag)
>> 883 {
>> 884 // take into account the supression due to the LPM effect
>> 885 if (G4UniformRand() <= SupressionFunction(aMaterial,KineticEnergy,GammaEnergy))
>> 886 LPMOK = true ;
>> 887 }
>> 888 else
>> 889 LPMOK = true ;
>> 890
>> 891 } while (!LPMOK) ;
>> 892
>> 893
>> 894 //protection: DO NOT PRODUCE a gamma with energy 0. !
>> 895 if (GammaEnergy <= 0.)
>> 896 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 897
>> 898 //
>> 899 // angles of the emitted gamma. ( Z - axis along the parent particle)
>> 900 //
>> 901 // universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
>> 902 // derived from Tsai distribution (Rev Mod Phys 49,421(1977))
>> 903
>> 904 G4double u;
>> 905 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
>> 906
>> 907 if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1 ;
>> 908 else u = - log(G4UniformRand()*G4UniformRand())/a2 ;
>> 909
>> 910 G4double Teta = u*electron_mass_c2/TotalEnergy ;
>> 911 G4double Phi = twopi * G4UniformRand() ;
>> 912 G4double dirx = sin(Teta)*cos(Phi) , diry = sin(Teta)*sin(Phi) , dirz = cos(Teta) ;
>> 913
>> 914 G4ThreeVector GammaDirection ( dirx, diry, dirz);
>> 915 GammaDirection.rotateUz(ParticleDirection);
>> 916
>> 917 // create G4DynamicParticle object for the Gamma
>> 918 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
>> 919 GammaDirection, GammaEnergy);
>> 920
>> 921 aParticleChange.SetNumberOfSecondaries(1);
>> 922 aParticleChange.AddSecondary(aGamma);
>> 923
>> 924 //
>> 925 // Update the incident particle
>> 926 //
>> 927
>> 928 G4double NewKinEnergy = KineticEnergy - GammaEnergy;
>> 929 if (NewKinEnergy > 0.)
>> 930 {
>> 931 aParticleChange.SetMomentumChange( ParticleDirection );
>> 932 aParticleChange.SetEnergyChange( NewKinEnergy );
>> 933 aParticleChange.SetLocalEnergyDeposit (0.);
>> 934 }
>> 935 else
>> 936 {
>> 937 aParticleChange.SetEnergyChange( 0. );
>> 938 aParticleChange.SetLocalEnergyDeposit (0.);
>> 939 if (charge<0.) aParticleChange.SetStatusChange(fStopAndKill);
>> 940 else aParticleChange.SetStatusChange(fStopButAlive);
>> 941 }
>> 942
>> 943 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
110 } 944 }
111 945
112 //....oooOO0OOooo........oooOO0OOooo........oo 946 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 947
114 void G4eBremsstrahlung::StreamProcessInfo(std: << 948 G4Element* G4eBremsstrahlung::SelectRandomAtom(G4Material* aMaterial) const
115 { 949 {
116 if(nullptr != EmModel(0)) { << 950 // select randomly 1 element within the material
117 G4EmParameters* param = G4EmParameters::In << 951
118 G4double eth = param->BremsstrahlungTh(); << 952 const G4int Index = aMaterial->GetIndex();
119 out << " LPM flag: " << param->LPM() << 953 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
120 << EmModel(0)->HighEnergyLimit()/GeV << " Ge << 954 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
121 if(eth < DBL_MAX) { << 955
122 out << ", VertexHighEnergyTh(GeV)= " << << 956 G4double rval = G4UniformRand()*((*PartialSumSigma[Index])[NumberOfElements-1]);
>> 957 for ( G4int i=0; i < NumberOfElements; i++ )
>> 958 if (rval <= (*PartialSumSigma[Index])[i]) return ((*theElementVector)[i]);
>> 959 G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
>> 960 << "' has no elements, NULL pointer returned." << G4endl;
>> 961 return NULL;
>> 962 }
>> 963
>> 964 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 965
>> 966
>> 967 G4double G4eBremsstrahlung::SupressionFunction(const G4Material* aMaterial,
>> 968 G4double KineticEnergy,G4double GammaEnergy)
>> 969 {
>> 970 // supression due to the LPM effect+polarisation of the medium/
>> 971 // supression due to the polarisation alone
>> 972 const G4double MigdalConstant = classic_electr_radius*
>> 973 electron_Compton_length*
>> 974 electron_Compton_length/pi ;
>> 975
>> 976 const G4double LPMconstant = fine_structure_const*electron_mass_c2*
>> 977 electron_mass_c2/(8.*pi*hbarc) ;
>> 978 G4double TotalEnergy,TotalEnergySquare,LPMEnergy,LPMGammaEnergyLimit,
>> 979 LPMGammaEnergyLimit2,GammaEnergySquare,sp,s2lpm,supr,w,splim,Cnorm ;
>> 980
>> 981 TotalEnergy = KineticEnergy+electron_mass_c2 ;
>> 982 TotalEnergySquare = TotalEnergy*TotalEnergy ;
>> 983
>> 984 LPMEnergy = LPMconstant*(aMaterial->GetRadlen()) ;
>> 985 LPMGammaEnergyLimit = TotalEnergySquare/LPMEnergy ;
>> 986 GammaEnergySquare = GammaEnergy*GammaEnergy ;
>> 987
>> 988 LPMGammaEnergyLimit2 = LPMGammaEnergyLimit*LPMGammaEnergyLimit ;
>> 989 splim = LPMGammaEnergyLimit2/(LPMGammaEnergyLimit2+MigdalConstant*TotalEnergySquare*
>> 990 (aMaterial->GetElectronDensity())) ;
>> 991 w = 1.+1./splim ;
>> 992 Cnorm = 2./(sqrt(w*w+4.)-w) ;
>> 993
>> 994 sp = GammaEnergySquare/(GammaEnergySquare+MigdalConstant*TotalEnergySquare*
>> 995 (aMaterial->GetElectronDensity())) ;
>> 996 if (LPMflag)
>> 997 {
>> 998 s2lpm = LPMEnergy*GammaEnergy/TotalEnergySquare;
>> 999 if (s2lpm < 1.)
>> 1000 {
>> 1001 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
>> 1002 else w = s2lpm*(1.+1./sp);
>> 1003 supr = Cnorm*(sqrt(w*w+4.*s2lpm)-w)/2. ;
>> 1004 }
>> 1005 else supr = sp;
123 } 1006 }
124 out << G4endl; << 1007 else supr = sp;
>> 1008
>> 1009 supr /= sp;
>> 1010 return supr;
>> 1011 }
>> 1012
>> 1013 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1014
>> 1015 G4bool G4eBremsstrahlung::StorePhysicsTable(G4ParticleDefinition* particle,
>> 1016 const G4String& directory,
>> 1017 G4bool ascii)
>> 1018 {
>> 1019 G4String filename;
>> 1020
>> 1021 // store stopping power table
>> 1022 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 1023 if ( !theLossTable->StorePhysicsTable(filename, ascii) ){
>> 1024 G4cout << " FAIL theLossTable->StorePhysicsTable in " << filename
>> 1025 << G4endl;
>> 1026 return false;
125 } 1027 }
>> 1028
>> 1029 // store mean free path table
>> 1030 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 1031 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 1032 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 1033 << G4endl;
>> 1034 return false;
>> 1035 }
>> 1036
>> 1037 // store PartialSumSigma table (G4OrderedTable)
>> 1038 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 1039 if ( !PartialSumSigma.Store(filename, ascii) ){
>> 1040 G4cout << " FAIL PartialSumSigma.store in " << filename
>> 1041 << G4endl;
>> 1042 return false;
>> 1043 }
>> 1044
>> 1045 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 1046 << ": Success to store the PhysicsTables in "
>> 1047 << directory << G4endl;
>> 1048 return true;
126 } 1049 }
127 1050
128 //....oooOO0OOooo........oooOO0OOooo........oo << 1051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 1052
130 void G4eBremsstrahlung::ProcessDescription(std << 1053 G4bool G4eBremsstrahlung::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 1054 const G4String& directory,
>> 1055 G4bool ascii)
131 { 1056 {
132 out << " Bremsstrahlung"; << 1057 // delete theLossTable and theMeanFreePathTable
133 G4VEnergyLossProcess::ProcessDescription(out << 1058 if (theLossTable != 0) {
>> 1059 theLossTable->clearAndDestroy();
>> 1060 delete theLossTable;
>> 1061 }
>> 1062 if (theMeanFreePathTable != 0) {
>> 1063 theMeanFreePathTable->clearAndDestroy();
>> 1064 delete theMeanFreePathTable;
>> 1065 }
>> 1066
>> 1067 if (&PartialSumSigma != 0) PartialSumSigma.clear();
>> 1068
>> 1069 // get bining from EnergyLoss
>> 1070 LowestKineticEnergy = GetLowerBoundEloss();
>> 1071 HighestKineticEnergy = GetUpperBoundEloss();
>> 1072 TotBin = GetNbinEloss();
>> 1073
>> 1074 G4String filename;
>> 1075
>> 1076 // retreive stopping power table
>> 1077 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 1078 theLossTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 1079 if ( !theLossTable->RetrievePhysicsTable(filename, ascii) ){
>> 1080 G4cout << " FAIL theLossTable0->RetrievePhysicsTable in " << filename
>> 1081 << G4endl;
>> 1082 return false;
>> 1083 }
>> 1084
>> 1085 // retreive mean free path table
>> 1086 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 1087 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 1088 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 1089 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 1090 << G4endl;
>> 1091 return false;
>> 1092 }
>> 1093
>> 1094 // retrieve PartialSumSigma table (G4OrderedTable)
>> 1095 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 1096 if ( !PartialSumSigma.Retrieve(filename, ascii) ){
>> 1097 G4cout << " FAIL PartialSumSigma.retrieve in " << filename
>> 1098 << G4endl;
>> 1099 return false;
>> 1100 }
>> 1101
>> 1102 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 1103 << ": Success to retrieve the PhysicsTables from "
>> 1104 << directory << G4endl;
>> 1105
>> 1106 if (particle==G4Electron::Electron())
>> 1107 {
>> 1108 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 1109 CounterOfElectronProcess++;
>> 1110 }
>> 1111 else
>> 1112 {
>> 1113 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 1114 CounterOfPositronProcess++;
>> 1115 }
>> 1116
>> 1117 BuildDEDXTable (*particle);
>> 1118
>> 1119 return true;
134 } 1120 }
>> 1121
>> 1122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 1123
136 //....oooOO0OOooo........oooOO0OOooo........oo << 1124 void G4eBremsstrahlung::PrintInfoDefinition()
>> 1125 {
>> 1126 G4String comments = "Total cross sections from a NEW parametrisation"
>> 1127 " based on the EEDL data library. "
>> 1128 "\n Good description from 1 KeV to 100 GeV.\n"
>> 1129 " log scale extrapolation above 100 GeV \n"
>> 1130 " Gamma energy sampled from a parametrised formula.";
>> 1131
>> 1132 G4cout << G4endl << GetProcessName() << ": " << comments
>> 1133 << "\n PhysicsTables from "
>> 1134 << G4BestUnit(LowerBoundLambda,"Energy")
>> 1135 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 1136 << " in " << NbinLambda << " bins. \n";
>> 1137 }
>> 1138
>> 1139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 1140