| 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.31 2003/04/29 04:59:48 kurasige Exp $
>> 25 // GEANT4 tag $Name: geant4-05-01-patch-01 $
28 // 26 //
29 // GEANT4 Class file <<
30 // 27 //
>> 28 // ------------ G4eBremsstrahlung physics process --------
>> 29 // by Michel Maire, 24 July 1996
31 // 30 //
32 // File name: G4eBremsstrahlung << 31 // 26-09-96 extension of the total crosssection above 100 GeV, M.Maire
33 // << 32 // 1-10-96 new type G4OrderedTable; ComputePartialSumSigma(), M.Maire
34 // Author: Michel Maire << 33 // 16-10-96 DoIt() call to the non static GetEnergyCuts(), L.Urban
35 // << 34 // 13-12-96 Sign corrected in grejmax and greject
36 // Creation date: 26.06.1996 << 35 // error definition of screenvar, L.Urban
37 // << 36 // 20-03-97 new energy loss+ionisation+brems scheme, L.Urban
38 // Modified by Michel Maire, Vladimir Ivanchen << 37 // 07-04-98 remove 'tracking cut' of the diffracted particle, MMa
39 // << 38 // 13-08-98 new methods SetBining() PrintInfo()
40 // ------------------------------------------- << 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 // 11-11-02 fix of division by 0 (VI)
>> 50 // 16-01-03 Migrade to cut per region (V.Ivanchenko)
>> 51 // 26-04-03 fix problems of retrieve tables (V.Ivanchenko)
41 // 52 //
>> 53 // --------------------------------------------------------------
>> 54
42 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 57
45 #include "G4eBremsstrahlung.hh" 58 #include "G4eBremsstrahlung.hh"
46 #include "G4SystemOfUnits.hh" << 59 #include "G4EnergyLossTables.hh"
47 #include "G4Gamma.hh" << 60 #include "G4ios.hh"
48 #include "G4SeltzerBergerModel.hh" <<
49 #include "G4eBremsstrahlungRelModel.hh" <<
50 #include "G4UnitsTable.hh" 61 #include "G4UnitsTable.hh"
51 #include "G4LossTableManager.hh" <<
52 <<
53 #include "G4ProductionCutsTable.hh" 62 #include "G4ProductionCutsTable.hh"
54 #include "G4MaterialCutsCouple.hh" <<
55 #include "G4EmParameters.hh" <<
56 63
57 //....oooOO0OOooo........oooOO0OOooo........oo << 64 G4double G4eBremsstrahlung::LowerBoundLambda = 1.*keV;
>> 65 G4double G4eBremsstrahlung::UpperBoundLambda = 100.*TeV;
>> 66 G4int G4eBremsstrahlung::NbinLambda = 100;
>> 67 G4double G4eBremsstrahlung::probsup = 1.00;
>> 68 G4bool G4eBremsstrahlung::LPMflag = true;
58 69
59 using namespace std; << 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 71
>> 72 // constructor
60 73
61 G4eBremsstrahlung::G4eBremsstrahlung(const G4S << 74 G4eBremsstrahlung::G4eBremsstrahlung(const G4String& processName)
62 G4VEnergyLossProcess(name) << 75 : G4VeEnergyLoss(processName), // initialization
>> 76 theMeanFreePathTable(NULL)
>> 77 { // MinThreshold = 10*keV;
>> 78 MinThreshold = 1*keV; }
>> 79
>> 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 81
>> 82 // destructor
>> 83
>> 84 G4eBremsstrahlung::~G4eBremsstrahlung()
63 { 85 {
64 SetProcessSubType(fBremsstrahlung); << 86 if (theMeanFreePathTable) {
65 SetSecondaryParticle(G4Gamma::Gamma()); << 87 theMeanFreePathTable->clearAndDestroy();
66 SetIonisation(false); << 88 delete theMeanFreePathTable;
67 SetCrossSectionType(fEmTwoPeaks); << 89 }
>> 90
>> 91 PartialSumSigma.clearAndDestroy();
>> 92 }
>> 93
>> 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 95
>> 96 void G4eBremsstrahlung::SetLowerBoundLambda(G4double val)
>> 97 {LowerBoundLambda = val;}
>> 98
>> 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 100
>> 101 void G4eBremsstrahlung::SetUpperBoundLambda(G4double val)
>> 102 {UpperBoundLambda = val;}
>> 103
>> 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 105
>> 106 void G4eBremsstrahlung::SetNbinLambda(G4int n)
>> 107 {NbinLambda = n;}
>> 108
>> 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 110
>> 111 G4double G4eBremsstrahlung::GetLowerBoundLambda()
>> 112 {return LowerBoundLambda;}
>> 113
>> 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 115
>> 116 G4double G4eBremsstrahlung::GetUpperBoundLambda()
>> 117 {return UpperBoundLambda;}
>> 118
>> 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 120
>> 121 G4int G4eBremsstrahlung::GetNbinLambda()
>> 122 {return NbinLambda;}
>> 123
>> 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 125
>> 126 void G4eBremsstrahlung::SetLPMflag(G4bool val)
>> 127 {LPMflag = val;}
>> 128
>> 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 130
>> 131 G4bool G4eBremsstrahlung::GetLPMflag()
>> 132 {return LPMflag;}
>> 133
>> 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 135
>> 136 void G4eBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
>> 137 {
>> 138 if( !CutsWhereModified() && theLossTable) return;
>> 139
>> 140 LowestKineticEnergy = GetLowerBoundEloss();
>> 141 HighestKineticEnergy = GetUpperBoundEloss();
>> 142 TotBin = GetNbinEloss();
>> 143
>> 144 BuildLossTable(aParticleType);
>> 145
>> 146 if (&aParticleType==G4Electron::Electron())
>> 147 {
>> 148 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 149 CounterOfElectronProcess++;
>> 150 }
>> 151 else
>> 152 {
>> 153 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 154 CounterOfPositronProcess++;
>> 155 }
>> 156
>> 157 BuildLambdaTable(aParticleType);
>> 158
>> 159 BuildDEDXTable (aParticleType);
>> 160
>> 161 if (&aParticleType==G4Electron::Electron()) PrintInfoDefinition();
68 } 162 }
69 163
70 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 165
72 G4eBremsstrahlung::~G4eBremsstrahlung() = defa << 166 void G4eBremsstrahlung::BuildLossTable(const G4ParticleDefinition& aParticleType)
>> 167 {
>> 168 G4double KineticEnergy,TotalEnergy,bremloss,Z,x,
>> 169 losslim,loss,rate,natom,Cut;
>> 170
>> 171 const G4double MinKinEnergy = 1.*keV;
>> 172 const G4double Thigh = 100.*GeV;
>> 173 const G4double Cuthigh = 50.*GeV;
>> 174 const G4double Factorhigh = 36./(1450.*GeV);
>> 175 const G4double coef1 = -0.5, coef2 = 2./9.;
>> 176
>> 177 G4double particleMass = aParticleType.GetPDGMass() ;
>> 178
>> 179 const G4ProductionCutsTable* theCoupleTable=
>> 180 G4ProductionCutsTable::GetProductionCutsTable();
>> 181 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 182
>> 183 if (theLossTable) {theLossTable->clearAndDestroy(); delete theLossTable;}
>> 184 theLossTable = new G4PhysicsTable(numOfCouples);
>> 185
>> 186 secondaryEnergyCuts = theCoupleTable->GetEnergyCutsVector(0);
>> 187
>> 188 // loop for materials
>> 189 //
>> 190 for (size_t J=0; J<numOfCouples; J++)
>> 191 {
>> 192
>> 193 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 194 LowestKineticEnergy,HighestKineticEnergy,TotBin);
>> 195
>> 196 // get elements in the material
>> 197 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
>> 198 const G4Material* material= couple->GetMaterial();
>> 199
>> 200 const G4ElementVector* theElementVector = material->GetElementVector();
>> 201 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
>> 202 const G4int NumberOfElements = material->GetNumberOfElements();
>> 203
>> 204 // loop for the kinetic energy values
>> 205 for (G4int i=0; i<TotBin; i++)
>> 206 {
>> 207 KineticEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 208 TotalEnergy = KineticEnergy+particleMass ;
>> 209 Cut = SecondaryEnergyThreshold(J);
>> 210 if (Cut < MinThreshold) Cut = MinThreshold;
>> 211 if (Cut > KineticEnergy) Cut = KineticEnergy;
>> 212
>> 213 bremloss = 0.;
>> 214
>> 215 if (KineticEnergy>MinKinEnergy)
>> 216 {
>> 217 // loop for elements in the material
>> 218 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 219 {
>> 220 Z=(*theElementVector)[iel]->GetZ();
>> 221 natom = theAtomicNumDensityVector[iel] ;
>> 222 if (KineticEnergy <= Thigh)
>> 223 {
>> 224 //loss for MinKinEnergy<KineticEnergy<=100 GeV
>> 225 x=log(TotalEnergy/particleMass);
>> 226 loss = ComputeBremLoss(Z,natom,KineticEnergy,Cut,x) ;
>> 227 if (&aParticleType==G4Positron::Positron())
>> 228 loss *= ComputePositronCorrFactorLoss(Z,KineticEnergy,Cut) ;
>> 229 }
>> 230 else
>> 231 {
>> 232 // extrapolation for KineticEnergy>100 GeV
>> 233 x=log(Thigh/particleMass) ;
>> 234 if (Cut<Thigh)
>> 235 {
>> 236 losslim = ComputeBremLoss(Z,natom,Thigh,Cut,x) ;
>> 237 if (&aParticleType==G4Positron::Positron())
>> 238 loss *= ComputePositronCorrFactorLoss(Z,Thigh,Cut) ;
>> 239 rate = Cut/TotalEnergy ;
>> 240 loss = losslim*(1.+coef1*rate+coef2*rate*rate) ;
>> 241 rate = Cut/Thigh ;
>> 242 loss /= (1.+coef1*rate+coef2*rate*rate) ;
>> 243 }
>> 244 else
>> 245 {
>> 246 losslim = ComputeBremLoss(Z,natom,Thigh,Cuthigh,x) ;
>> 247 if (&aParticleType==G4Positron::Positron())
>> 248 loss *= ComputePositronCorrFactorLoss(Z,Thigh,Cuthigh) ;
>> 249 rate = Cut/TotalEnergy ;
>> 250 loss = losslim*(1.+coef1*rate+coef2*rate*rate) ;
>> 251 loss *= Factorhigh*Cut ;
>> 252 }
>> 253
>> 254 }
>> 255 bremloss += natom*loss;
>> 256 }
>> 257
>> 258 }
>> 259
>> 260 static const G4double MigdalConstant = classic_electr_radius
>> 261 *electron_Compton_length
>> 262 *electron_Compton_length/pi;
>> 263 G4double TotalEnergy = KineticEnergy+electron_mass_c2 ;
>> 264 G4double kp2 = MigdalConstant*TotalEnergy*TotalEnergy*
>> 265 (material->GetElectronDensity()) ;
>> 266
>> 267 // now compute the correction due to the supression(s)
>> 268 G4double kmin = 1.*eV ;
>> 269 G4double kmax = Cut ;
>> 270
>> 271 if(kmax > kmin)
>> 272 {
>> 273
>> 274 G4double floss = 0. ;
>> 275 G4int nmax = 100 ;
>> 276
>> 277 G4double vmin=log(kmin);
>> 278 G4double vmax=log(kmax) ;
>> 279 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(HighestKineticEnergy)-vmin)) ;
>> 280 G4double u,fac,c,v,dv ;
>> 281 if(nn > 0)
>> 282 {
>> 283 dv = (vmax-vmin)/nn ;
>> 284 v = vmin-dv ;
>> 285 for(G4int n=0; n<=nn; n++)
>> 286 {
>> 287 v += dv; u = exp(v);
>> 288 fac = u*SupressionFunction(material,KineticEnergy,u);
>> 289 probsup = 1.;
>> 290 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
>> 291 if ((n==0)||(n==nn)) c=0.5;
>> 292 else c=1. ;
>> 293 fac *= c ;
>> 294 floss += fac ;
>> 295 }
>> 296 floss *=dv/(kmax-kmin);
>> 297 }
>> 298 else floss = 1.;
>> 299 if(floss > 1.) floss = 1.;
>> 300 // correct the loss
>> 301 bremloss *= floss;
>> 302 }
>> 303
>> 304 if(bremloss < 0.) bremloss = 0.;
>> 305 aVector->PutValue(i,bremloss);
>> 306 }
>> 307
>> 308 theLossTable->insert(aVector);
>> 309 }
>> 310
>> 311 }
73 312
74 //....oooOO0OOooo........oooOO0OOooo........oo 313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 314
76 G4bool G4eBremsstrahlung::IsApplicable(const G << 315 G4double G4eBremsstrahlung::ComputeBremLoss(G4double Z,G4double natom,
>> 316 G4double T,G4double Cut,G4double x)
>> 317
>> 318 // compute loss due to soft brems
77 { 319 {
78 return (&p == G4Electron::Electron() || &p = << 320 static const G4double beta=1.00,ksi=2.00 ;
>> 321 static const G4double clossh = 0.254 , closslow = 1./3. , alosslow = 1. ;
>> 322 static const G4double Tlim= 10.*MeV ;
>> 323
>> 324 static const G4double xlim = 1.2 ;
>> 325 static const G4int NZ = 8 ;
>> 326 static const G4int Nloss = 11 ;
>> 327 static const G4double ZZ[NZ] =
>> 328 {2.,4.,6.,14.,26.,50.,82.,92.};
>> 329 static const G4double coefloss[NZ][Nloss] = {
>> 330 // Z=2
>> 331 { 0.98916, 0.47564, -0.2505, -0.45186, 0.14462,
>> 332 0.21307, -0.013738, -0.045689, -0.0042914, 0.0034429,
>> 333 0.00064189},
>> 334
>> 335 // Z=4
>> 336 { 1.0626, 0.37662, -0.23646, -0.45188, 0.14295,
>> 337 0.22906, -0.011041, -0.051398, -0.0055123, 0.0039919,
>> 338 0.00078003},
>> 339 // Z=6
>> 340 { 1.0954, 0.315, -0.24011, -0.43849, 0.15017,
>> 341 0.23001, -0.012846, -0.052555, -0.0055114, 0.0041283,
>> 342 0.00080318},
>> 343
>> 344 // Z=14
>> 345 { 1.1649, 0.18976, -0.24972, -0.30124, 0.1555,
>> 346 0.13565, -0.024765, -0.027047, -0.00059821, 0.0019373,
>> 347 0.00027647},
>> 348
>> 349 // Z=26
>> 350 { 1.2261, 0.14272, -0.25672, -0.28407, 0.13874,
>> 351 0.13586, -0.020562, -0.026722, -0.00089557, 0.0018665,
>> 352 0.00026981},
>> 353
>> 354 // Z=50
>> 355 { 1.3147, 0.020049, -0.35543, -0.13927, 0.17666,
>> 356 0.073746, -0.036076, -0.013407, 0.0025727, 0.00084005,
>> 357 -1.4082e-05},
>> 358
>> 359 // Z=82
>> 360 { 1.3986, -0.10586, -0.49187, -0.0048846, 0.23621,
>> 361 0.031652, -0.052938, -0.0076639, 0.0048181, 0.00056486,
>> 362 -0.00011995},
>> 363
>> 364 // Z=92
>> 365 { 1.4217, -0.116, -0.55497, -0.044075, 0.27506,
>> 366 0.081364, -0.058143, -0.023402, 0.0031322, 0.0020201,
>> 367 0.00017519}
>> 368
>> 369 } ;
>> 370 static G4double aaa=0.414 ;
>> 371 static G4double bbb=0.345 ;
>> 372 static G4double ccc=0.460 ;
>> 373
>> 374 G4int iz = 0;
>> 375 G4double delz = 1.e6;
>> 376 for (G4int ii=0; ii<NZ; ii++)
>> 377 {
>> 378 if(abs(Z-ZZ[ii]) < delz) { iz = ii; delz = abs(Z-ZZ[ii]);}
>> 379 }
>> 380
>> 381 G4double xx = log10(T);
>> 382 G4double fl = 1.;
>> 383
>> 384 if (xx <= xlim)
>> 385 {
>> 386 fl = coefloss[iz][Nloss-1];
>> 387 for (G4int j=Nloss-2; j>=0; j--) fl = fl*xx+coefloss[iz][j];
>> 388 if (fl < 0.) fl = 0.;
>> 389 }
>> 390
>> 391 G4double loss;
>> 392 G4double E = T+electron_mass_c2 ;
>> 393
>> 394 loss = Z*(Z+ksi)*E*E/(T+E)*exp(beta*log(Cut/T))*(2.-clossh*exp(log(Z)/4.));
>> 395 if (T <= Tlim) loss /= exp(closslow*log(Tlim/T));
>> 396 if( T <= Cut) loss *= exp(alosslow*log(T/Cut));
>> 397
>> 398 // correction ................................
>> 399 loss *= (aaa+bbb*T/Tlim)/(1.+ccc*T/Tlim);
>> 400 loss *= fl;
>> 401 loss /= Avogadro;
>> 402
>> 403 return loss;
>> 404 }
>> 405
>> 406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 407
>> 408 G4double G4eBremsstrahlung::ComputePositronCorrFactorLoss(
>> 409 G4double Z,G4double KineticEnergy,G4double GammaCut)
>> 410
>> 411 //calculates the correction factor for the energy loss due to bremsstrahlung for positrons
>> 412 //the same correction is in the (discrete) bremsstrahlung
>> 413
>> 414 {
>> 415 static const G4double K = 132.9416*eV ;
>> 416 static const G4double a1=4.15e-1, a3=2.10e-3, a5=54.0e-5 ;
>> 417
>> 418 G4double x = log(KineticEnergy/(K*Z*Z)), x2 = x*x, x3 = x2*x;
>> 419 G4double eta = 0.5+atan(a1*x+a3*x3+a5*x3*x2)/pi;
>> 420 G4double e0 = GammaCut/KineticEnergy;
>> 421
>> 422 G4double factor(0.);
>> 423 if (e0!=1.0) { factor=log(1.-e0)/eta; factor=exp(factor);}
>> 424 factor = eta*(1.-factor)/e0;
>> 425
>> 426 return factor;
>> 427 }
>> 428
>> 429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 430
>> 431 void G4eBremsstrahlung::BuildLambdaTable(
>> 432 const G4ParticleDefinition& ParticleType)
>> 433 {
>> 434 G4double LowEdgeEnergy , Value;
>> 435 G4double FixedEnergy = (LowerBoundLambda + UpperBoundLambda)/2.;
>> 436
>> 437 //create table
>> 438 //
>> 439 const G4ProductionCutsTable* theCoupleTable=
>> 440 G4ProductionCutsTable::GetProductionCutsTable();
>> 441 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 442
>> 443 //create table
>> 444 if (theMeanFreePathTable) {theMeanFreePathTable->clearAndDestroy();
>> 445 delete theMeanFreePathTable;
>> 446 }
>> 447 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
>> 448
>> 449 PartialSumSigma.clearAndDestroy();
>> 450 PartialSumSigma.resize(numOfCouples);
>> 451
>> 452 G4PhysicsLogVector* ptrVector;
>> 453 for ( size_t J=0; J<numOfCouples; J++ )
>> 454 {
>> 455 //create physics vector then fill it ....
>> 456 ptrVector = new G4PhysicsLogVector(LowerBoundLambda, UpperBoundLambda,
>> 457 NbinLambda ) ;
>> 458
>> 459 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
>> 460
>> 461 for ( G4int i = 0 ; i < NbinLambda ; i++ )
>> 462 {
>> 463 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 464 Value = ComputeMeanFreePath( &ParticleType, LowEdgeEnergy,couple);
>> 465 ptrVector->PutValue( i , Value ) ;
>> 466 }
>> 467
>> 468 theMeanFreePathTable->insertAt( J , ptrVector );
>> 469
>> 470 // Compute the PartialSumSigma table at a given fixed energy
>> 471 ComputePartialSumSigma( &ParticleType, FixedEnergy, couple) ;
>> 472 }
>> 473
>> 474 }
>> 475
>> 476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 477
>> 478 G4double G4eBremsstrahlung::ComputeMeanFreePath(
>> 479 const G4ParticleDefinition* ParticleType,
>> 480 G4double KineticEnergy,
>> 481 const G4MaterialCutsCouple* couple)
>> 482 {
>> 483 const G4Material* aMaterial = couple->GetMaterial();
>> 484 const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
>> 485 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 486 G4double GammaEnergyCut = SecondaryEnergyThreshold(couple->GetIndex());
>> 487 if (GammaEnergyCut < MinThreshold) GammaEnergyCut = MinThreshold;
>> 488
>> 489 G4double SIGMA = 0;
>> 490
>> 491 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
>> 492 {
>> 493 SIGMA += theAtomNumDensityVector[i] *
>> 494 ComputeCrossSectionPerAtom( ParticleType, KineticEnergy,
>> 495 (*theElementVector)[i]->GetZ(),
>> 496 GammaEnergyCut );
>> 497 }
>> 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 G4double vmin=log(kmin);
>> 516 G4double vmax=log(kmax) ;
>> 517 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(HighestKineticEnergy)-vmin));
>> 518 G4double u,fac,c,v,dv,y ;
>> 519 if(nn > 0)
>> 520 {
>> 521 dv = (vmax-vmin)/nn ;
>> 522 v = vmin-dv ;
>> 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},
89 582
90 G4double emax = param->MaxKinEnergy(); << 583 // Z=6
91 G4VEmFluctuationModel* fm = nullptr; << 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},
92 587
93 if (nullptr == EmModel(0)) { SetEmModel(ne << 588 // Z=14
94 G4double energyLimit = std::min(EmModel(0) << 589 { 0.55058, 0.25629, 0.35854, -0.080656, -0.054308,
95 EmModel(0)->SetHighEnergyLimit(energyLimit << 590 -0.049933, -0.00064246, 0.016597, 0.0021789, -0.001327,
96 EmModel(0)->SetSecondaryThreshold(param->B << 591 -0.00025983},
97 AddEmModel(1, EmModel(0), fm); <<
98 592
99 if(emax > energyLimit) { << 593 // Z=26
100 if (nullptr == EmModel(1)) { << 594 { 0.5791, 0.26152, 0.38953, -0.17104, -0.099172,
101 SetEmModel(new G4eBremsstrahlungRelModel()); << 595 0.024596, 0.023718, -0.0039205, -0.0036658, 0.00041749,
102 } << 596 0.00023408},
103 EmModel(1)->SetLowEnergyLimit(energyLimi << 597
104 EmModel(1)->SetHighEnergyLimit(emax); << 598 // Z=50
105 EmModel(1)->SetSecondaryThreshold(param- << 599 { 0.62085, 0.27045, 0.39073, -0.37916, -0.18878,
106 AddEmModel(1, EmModel(1), fm); << 600 0.23905, 0.095028, -0.068744, -0.023809, 0.0062408,
>> 601 0.0020407},
>> 602
>> 603 // Z=82
>> 604 { 0.66053, 0.24513, 0.35404, -0.47275, -0.22837,
>> 605 0.35647, 0.13203, -0.1049, -0.034851, 0.0095046,
>> 606 0.0030535},
>> 607
>> 608 // Z=92
>> 609 { 0.67143, 0.23079, 0.32256, -0.46248, -0.20013,
>> 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 G4MaterialCutsCouple* couple)
>> 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 const G4Material* aMaterial= couple->GetMaterial();
>> 688 G4int NbOfElements = aMaterial->GetNumberOfElements();
>> 689 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 690 const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 691 size_t index = couple->GetIndex();
>> 692 G4double GammaEnergyCut = SecondaryEnergyThreshold(index);
>> 693
>> 694 PartialSumSigma[index] = 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[index]->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
>> 758 //G4double LPMEnergy = LPMconstant*(aMaterial->GetRadlen()) ;
>> 759
>> 760 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
>> 761 const G4Material* aMaterial = couple->GetMaterial();
>> 762 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
>> 763 G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();
>> 764
>> 765 G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
>> 766 G4ParticleMomentum ParticleDirection = aDynamicParticle->GetMomentumDirection();
>> 767
>> 768 // Gamma production cut in this material
>> 769
>> 770 G4double GammaEnergyCut = SecondaryEnergyThreshold(couple->GetIndex());
>> 771 if (GammaEnergyCut < MinThreshold) GammaEnergyCut = MinThreshold;
>> 772
>> 773 // check against insufficient energy
>> 774 if (KineticEnergy < GammaEnergyCut)
>> 775 {
>> 776 aParticleChange.SetMomentumChange( ParticleDirection );
>> 777 aParticleChange.SetEnergyChange( KineticEnergy );
>> 778 aParticleChange.SetLocalEnergyDeposit (0.);
>> 779 aParticleChange.SetNumberOfSecondaries(0);
>> 780 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 781 }
>> 782
>> 783 // select randomly one element constituing the material
>> 784 G4Element* anElement = SelectRandomAtom(couple);
>> 785
>> 786 // Extract Z factors for this Element
>> 787 G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
>> 788 G4double FZ = lnZ* (4.- 0.55*lnZ);
>> 789 G4double ZZ = anElement->GetIonisation()->GetZZ3();
>> 790
>> 791 // limits of the energy sampling
>> 792 G4double TotalEnergy = KineticEnergy + electron_mass_c2;
>> 793 G4double xmin = GammaEnergyCut/KineticEnergy;
>> 794 G4double epsilmin = GammaEnergyCut/TotalEnergy;
>> 795 G4double epsilmax = KineticEnergy/TotalEnergy;
>> 796
>> 797 // Migdal factor
>> 798 G4double
>> 799 MigdalFactor = (aMaterial->GetElectronDensity())*MigdalConstant
>> 800 /(epsilmax*epsilmax);
>> 801
>> 802 //
>> 803 G4double x, epsil, greject, migdal, grejmax;
>> 804 G4double U = log(KineticEnergy/electron_mass_c2), U2 = U*U;
>> 805
>> 806 //
>> 807 // sample the energy rate of the emitted gamma for electron kinetic energy > 1 MeV
>> 808 //
>> 809
>> 810 do {
>> 811 if (KineticEnergy > 1.*MeV)
>> 812 {
>> 813 // parameters
>> 814 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12),
>> 815 ah2 = ah20 + ZZ* (ah21 + ZZ* ah22),
>> 816 ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
>> 817
>> 818 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12),
>> 819 bh2 = bh20 + ZZ* (bh21 + ZZ* bh22),
>> 820 bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
>> 821
>> 822 G4double ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
>> 823 G4double bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
>> 824
>> 825 // limit of the screening variable
>> 826 G4double screenfac =
>> 827 136.*electron_mass_c2/((anElement->GetIonisation()->GetZ3())*TotalEnergy);
>> 828 G4double screenmin = screenfac*epsilmin/(1.-epsilmin);
>> 829
>> 830 // Compute the maximum of the rejection function
>> 831 G4double F1 = G4std::max(ScreenFunction1(screenmin) - FZ ,0.);
>> 832 G4double F2 = G4std::max(ScreenFunction2(screenmin) - FZ ,0.);
>> 833 grejmax = (F1 - epsilmin* (F1*ah - bh*epsilmin*F2))/(42.392 - FZ);
>> 834
>> 835 // sample the energy rate of the emitted Gamma
>> 836 G4double screenvar;
>> 837
>> 838
>> 839 do {
>> 840
>> 841 x = pow(xmin, G4UniformRand());
>> 842 epsil = x*KineticEnergy/TotalEnergy;
>> 843 screenvar = screenfac*epsil/(1-epsil);
>> 844 F1 = G4std::max(ScreenFunction1(screenvar) - FZ ,0.);
>> 845 F2 = G4std::max(ScreenFunction2(screenvar) - FZ ,0.);
>> 846 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
>> 847 greject = migdal*(F1 - epsil* (ah*F1 - bh*epsil*F2))/(42.392 - FZ);
>> 848 } while( greject < G4UniformRand()*grejmax );
>> 849
>> 850 }
>> 851
>> 852 else
>> 853 {
>> 854 // sample the energy rate of the emitted gamma for electron kinetic energy < 1 MeV
>> 855 //
>> 856 // parameters
>> 857 G4double al0 = al00 + ZZ* (al01 + ZZ* al02),
>> 858 al1 = al10 + ZZ* (al11 + ZZ* al12),
>> 859 al2 = al20 + ZZ* (al21 + ZZ* al22);
>> 860
>> 861 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02),
>> 862 bl1 = bl10 + ZZ* (bl11 + ZZ* bl12),
>> 863 bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
>> 864
>> 865 G4double al = al0 + al1*U + al2*U2;
>> 866 G4double bl = bl0 + bl1*U + bl2*U2;
>> 867
>> 868 // Compute the maximum of the rejection function
>> 869 grejmax = G4std::max(1. + xmin* (al + bl*xmin), 1.+al+bl);
>> 870 G4double xm = -al/(2.*bl);
>> 871 if ((xmin < xm)&&(xm < 1.)) grejmax = G4std::max(grejmax, 1.+ xm* (al + bl*xm));
>> 872
>> 873 // sample the energy rate of the emitted Gamma
>> 874
>> 875 do { x = pow(xmin, G4UniformRand());
>> 876 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
>> 877 greject = migdal*(1. + x* (al + bl*x));
>> 878 } while( greject < G4UniformRand()*grejmax );
>> 879 }
>> 880
>> 881 GammaEnergy = x*KineticEnergy;
>> 882
>> 883 if(LPMflag)
>> 884 {
>> 885 // take into account the supression due to the LPM effect
>> 886 if (G4UniformRand() <= SupressionFunction(aMaterial,KineticEnergy,GammaEnergy))
>> 887 LPMOK = true ;
>> 888 }
>> 889 else
>> 890 LPMOK = true ;
>> 891
>> 892 } while (!LPMOK) ;
>> 893
>> 894
>> 895 //protection: DO NOT PRODUCE a gamma with energy 0. !
>> 896 if (GammaEnergy <= 0.)
>> 897 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 898
>> 899 //
>> 900 // angles of the emitted gamma. ( Z - axis along the parent particle)
>> 901 //
>> 902 // universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
>> 903 // derived from Tsai distribution (Rev Mod Phys 49,421(1977))
>> 904
>> 905 G4double u;
>> 906 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
>> 907
>> 908 if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1 ;
>> 909 else u = - log(G4UniformRand()*G4UniformRand())/a2 ;
>> 910
>> 911 G4double Teta = u*electron_mass_c2/TotalEnergy ;
>> 912 G4double Phi = twopi * G4UniformRand() ;
>> 913 G4double dirx = sin(Teta)*cos(Phi) , diry = sin(Teta)*sin(Phi) , dirz = cos(Teta) ;
>> 914
>> 915 G4ThreeVector GammaDirection ( dirx, diry, dirz);
>> 916 GammaDirection.rotateUz(ParticleDirection);
>> 917
>> 918 // create G4DynamicParticle object for the Gamma
>> 919 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
>> 920 GammaDirection, GammaEnergy);
>> 921
>> 922 aParticleChange.SetNumberOfSecondaries(1);
>> 923 aParticleChange.AddSecondary(aGamma);
>> 924
>> 925 //
>> 926 // Update the incident particle
>> 927 //
>> 928
>> 929 G4double NewKinEnergy = KineticEnergy - GammaEnergy;
>> 930 if (NewKinEnergy > 0.)
>> 931 {
>> 932 aParticleChange.SetMomentumChange( ParticleDirection );
>> 933 aParticleChange.SetEnergyChange( NewKinEnergy );
>> 934 aParticleChange.SetLocalEnergyDeposit (0.);
>> 935 }
>> 936 else
>> 937 {
>> 938 aParticleChange.SetEnergyChange( 0. );
>> 939 aParticleChange.SetLocalEnergyDeposit (0.);
>> 940 if (charge<0.) aParticleChange.SetStatusChange(fStopAndKill);
>> 941 else aParticleChange.SetStatusChange(fStopButAlive);
>> 942 }
>> 943
>> 944 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
110 } 945 }
111 946
112 //....oooOO0OOooo........oooOO0OOooo........oo 947 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 948
114 void G4eBremsstrahlung::StreamProcessInfo(std: << 949 G4Element* G4eBremsstrahlung::SelectRandomAtom(const G4MaterialCutsCouple* couple) const
115 { 950 {
116 if(nullptr != EmModel(0)) { << 951 // select randomly 1 element within the material
117 G4EmParameters* param = G4EmParameters::In << 952
118 G4double eth = param->BremsstrahlungTh(); << 953 size_t index = couple->GetIndex();
119 out << " LPM flag: " << param->LPM() << 954 const G4Material* aMaterial = couple->GetMaterial();
120 << EmModel(0)->HighEnergyLimit()/GeV << " Ge << 955 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
121 if(eth < DBL_MAX) { << 956 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
122 out << ", VertexHighEnergyTh(GeV)= " << << 957
>> 958 G4double rval = G4UniformRand()*((*PartialSumSigma[index])[NumberOfElements-1]);
>> 959 for ( G4int i=0; i < NumberOfElements; i++ )
>> 960 {
>> 961 if (rval <= (*PartialSumSigma[index])[i]) return ((*theElementVector)[i]);
>> 962 }
>> 963 G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
>> 964 << "' has no elements, NULL pointer returned." << G4endl;
>> 965 return 0;
>> 966 }
>> 967
>> 968 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 969
>> 970
>> 971 G4double G4eBremsstrahlung::SupressionFunction(const G4Material* aMaterial,
>> 972 G4double KineticEnergy,G4double GammaEnergy)
>> 973 {
>> 974 // supression due to the LPM effect+polarisation of the medium/
>> 975 // supression due to the polarisation alone
>> 976 const G4double MigdalConstant = classic_electr_radius*
>> 977 electron_Compton_length*
>> 978 electron_Compton_length/pi ;
>> 979
>> 980 const G4double LPMconstant = fine_structure_const*electron_mass_c2*
>> 981 electron_mass_c2/(8.*pi*hbarc) ;
>> 982 G4double TotalEnergy,TotalEnergySquare,LPMEnergy,LPMGammaEnergyLimit,
>> 983 LPMGammaEnergyLimit2,GammaEnergySquare,sp,s2lpm,supr,w,splim,Cnorm ;
>> 984
>> 985 TotalEnergy = KineticEnergy+electron_mass_c2 ;
>> 986 TotalEnergySquare = TotalEnergy*TotalEnergy ;
>> 987
>> 988 LPMEnergy = LPMconstant*(aMaterial->GetRadlen()) ;
>> 989 LPMGammaEnergyLimit = TotalEnergySquare/LPMEnergy ;
>> 990 GammaEnergySquare = GammaEnergy*GammaEnergy ;
>> 991
>> 992 LPMGammaEnergyLimit2 = LPMGammaEnergyLimit*LPMGammaEnergyLimit ;
>> 993 splim = LPMGammaEnergyLimit2/(LPMGammaEnergyLimit2+MigdalConstant*TotalEnergySquare*
>> 994 (aMaterial->GetElectronDensity())) ;
>> 995 w = 1.+1./splim ;
>> 996
>> 997 sp = GammaEnergySquare/(GammaEnergySquare+MigdalConstant*TotalEnergySquare*
>> 998 (aMaterial->GetElectronDensity())) ;
>> 999 if (LPMflag)
>> 1000 {
>> 1001 s2lpm = LPMEnergy*GammaEnergy/TotalEnergySquare;
>> 1002 if (s2lpm < 1.)
>> 1003 {
>> 1004 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
>> 1005 else w = s2lpm*(1.+1./sp);
>> 1006
>> 1007 Cnorm = 2./(sqrt(w*w+4.)-w) ;
>> 1008 supr = Cnorm*(sqrt(w*w+4.*s2lpm)-w)/2. ;
>> 1009 }
>> 1010 else supr = sp;
123 } 1011 }
124 out << G4endl; << 1012 else supr = sp;
>> 1013
>> 1014 supr /= sp;
>> 1015 return supr;
>> 1016 }
>> 1017
>> 1018 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1019
>> 1020 G4bool G4eBremsstrahlung::StorePhysicsTable(G4ParticleDefinition* particle,
>> 1021 const G4String& directory,
>> 1022 G4bool ascii)
>> 1023 {
>> 1024 G4String filename;
>> 1025
>> 1026 // store stopping power table
>> 1027 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 1028 if ( !theLossTable->StorePhysicsTable(filename, ascii) ){
>> 1029 G4cout << " FAIL theLossTable->StorePhysicsTable in " << filename
>> 1030 << G4endl;
>> 1031 return false;
>> 1032 }
>> 1033
>> 1034 // store mean free path table
>> 1035 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 1036 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 1037 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 1038 << G4endl;
>> 1039 return false;
>> 1040 }
>> 1041
>> 1042 // store PartialSumSigma table (G4OrderedTable)
>> 1043 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 1044 if ( !PartialSumSigma.Store(filename, ascii) ){
>> 1045 G4cout << " FAIL PartialSumSigma.store in " << filename
>> 1046 << G4endl;
>> 1047 return false;
>> 1048 }
>> 1049
>> 1050 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 1051 << ": Success to store the PhysicsTables in "
>> 1052 << directory << G4endl;
>> 1053 return true;
>> 1054 }
>> 1055
>> 1056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1057
>> 1058 G4bool G4eBremsstrahlung::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 1059 const G4String& directory,
>> 1060 G4bool ascii)
>> 1061 {
>> 1062 // delete theLossTable and theMeanFreePathTable
>> 1063 if (theLossTable != 0) {
>> 1064 theLossTable->clearAndDestroy();
>> 1065 delete theLossTable;
>> 1066 }
>> 1067 if (theMeanFreePathTable != 0) {
>> 1068 theMeanFreePathTable->clearAndDestroy();
>> 1069 delete theMeanFreePathTable;
125 } 1070 }
>> 1071
>> 1072
>> 1073 // get bining from EnergyLoss
>> 1074 LowestKineticEnergy = GetLowerBoundEloss();
>> 1075 HighestKineticEnergy = GetUpperBoundEloss();
>> 1076 TotBin = GetNbinEloss();
>> 1077
>> 1078 G4String filename;
>> 1079 const G4ProductionCutsTable* theCoupleTable=
>> 1080 G4ProductionCutsTable::GetProductionCutsTable();
>> 1081 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 1082
>> 1083 secondaryEnergyCuts = theCoupleTable->GetEnergyCutsVector(0);
>> 1084
>> 1085 // retreive stopping power table
>> 1086 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 1087 theLossTable = new G4PhysicsTable(numOfCouples);
>> 1088 if ( !theLossTable->RetrievePhysicsTable(filename, ascii) ){
>> 1089 G4cout << " FAIL theLossTable0->RetrievePhysicsTable in " << filename
>> 1090 << G4endl;
>> 1091 return false;
>> 1092 }
>> 1093
>> 1094 // retreive mean free path table
>> 1095 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 1096 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
>> 1097 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 1098 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 1099 << G4endl;
>> 1100 return false;
>> 1101 }
>> 1102
>> 1103 // retrieve PartialSumSigma table (G4OrderedTable)
>> 1104 PartialSumSigma.clearAndDestroy();
>> 1105 PartialSumSigma.reserve(numOfCouples);
>> 1106 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 1107 if ( !PartialSumSigma.Retrieve(filename, ascii) ){
>> 1108 G4cout << " FAIL PartialSumSigma.retrieve in " << filename
>> 1109 << G4endl;
>> 1110 return false;
>> 1111 }
>> 1112
>> 1113 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 1114 << ": Success to retrieve the PhysicsTables from "
>> 1115 << directory << G4endl;
>> 1116
>> 1117 if (particle==G4Electron::Electron())
>> 1118 {
>> 1119 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 1120 CounterOfElectronProcess++;
>> 1121 }
>> 1122 else
>> 1123 {
>> 1124 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 1125 CounterOfPositronProcess++;
>> 1126 }
>> 1127
>> 1128 BuildDEDXTable (*particle);
>> 1129
>> 1130 if (particle==G4Electron::Electron()) PrintInfoDefinition();
>> 1131 return true;
126 } 1132 }
127 1133
128 //....oooOO0OOooo........oooOO0OOooo........oo << 1134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 1135
130 void G4eBremsstrahlung::ProcessDescription(std << 1136 void G4eBremsstrahlung::PrintInfoDefinition()
131 { 1137 {
132 out << " Bremsstrahlung"; << 1138 G4String comments = "Total cross sections from a NEW parametrisation"
133 G4VEnergyLossProcess::ProcessDescription(out << 1139 " based on the EEDL data library. "
>> 1140 "\n Good description from 1 KeV to 100 GeV.\n"
>> 1141 " log scale extrapolation above 100 GeV \n"
>> 1142 " Gamma energy sampled from a parametrised formula.";
>> 1143
>> 1144 G4cout << G4endl << GetProcessName() << ": " << comments
>> 1145 << "\n PhysicsTables from "
>> 1146 << G4BestUnit(LowerBoundLambda,"Energy")
>> 1147 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 1148 << " in " << NbinLambda << " bins. \n";
134 } 1149 }
135 1150
136 //....oooOO0OOooo........oooOO0OOooo........oo << 1151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 1152