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