| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4eBremParametrizedModel.cc 91726 2015-08-03 15:41:36Z gcosmo $
>> 27 // GEANT4 tag $Name: geant4-09-04 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4eBremParametrizedModel 34 // File name: G4eBremParametrizedModel
33 // 35 //
34 // Author: Andreas Schaelicke 36 // Author: Andreas Schaelicke
35 // 37 //
36 // Creation date: 06.04.2011 38 // Creation date: 06.04.2011
37 // 39 //
38 // Modifications: 40 // Modifications:
39 // 41 //
40 // Main References: 42 // Main References:
41 // - based on G4eBremsstrahlungModel and G4eB 43 // - based on G4eBremsstrahlungModel and G4eBremsstrahlungRelModel
42 // ------------------------------------------- 44 // -------------------------------------------------------------------
43 // 45 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 48
47 #include "G4eBremParametrizedModel.hh" 49 #include "G4eBremParametrizedModel.hh"
48 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
50 #include "G4Electron.hh" 52 #include "G4Electron.hh"
51 #include "G4Positron.hh" 53 #include "G4Positron.hh"
52 #include "G4Gamma.hh" 54 #include "G4Gamma.hh"
53 #include "Randomize.hh" 55 #include "Randomize.hh"
54 #include "G4Material.hh" 56 #include "G4Material.hh"
55 #include "G4Element.hh" 57 #include "G4Element.hh"
56 #include "G4ElementVector.hh" 58 #include "G4ElementVector.hh"
57 #include "G4ProductionCutsTable.hh" 59 #include "G4ProductionCutsTable.hh"
58 #include "G4ParticleChangeForLoss.hh" 60 #include "G4ParticleChangeForLoss.hh"
59 #include "G4LossTableManager.hh" 61 #include "G4LossTableManager.hh"
60 #include "G4ModifiedTsai.hh" 62 #include "G4ModifiedTsai.hh"
61 #include "G4Exp.hh" 63 #include "G4Exp.hh"
62 #include "G4Log.hh" 64 #include "G4Log.hh"
63 65
64 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 67
66 const G4double G4eBremParametrizedModel::xgi[] 68 const G4double G4eBremParametrizedModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
67 0.5917, 0.7628, 0.8983, 0.9801 } 69 0.5917, 0.7628, 0.8983, 0.9801 };
68 const G4double G4eBremParametrizedModel::wgi[] 70 const G4double G4eBremParametrizedModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
69 0.1813, 0.1569, 0.1112, 0.0506 } 71 0.1813, 0.1569, 0.1112, 0.0506 };
70 72
71 static const G4double tlow = 1.*CLHEP::MeV; 73 static const G4double tlow = 1.*CLHEP::MeV;
72 74
73 // 75 //
74 // GEANT4 internal units. 76 // GEANT4 internal units.
75 // 77 //
76 static const G4double 78 static const G4double
77 ah10 = 4.67733E+00, ah11 =-6.19012E-01, a 79 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
78 ah20 =-7.34101E+00, ah21 = 1.00462E+00, a 80 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
79 ah30 = 2.93119E+00, ah31 =-4.03761E-01, a 81 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
80 82
81 static const G4double 83 static const G4double
82 bh10 = 4.23071E+00, bh11 =-6.10995E-01, b 84 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
83 bh20 =-7.12527E+00, bh21 = 9.69160E-01, b 85 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
84 bh30 = 2.69925E+00, bh31 =-3.63283E-01, b 86 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
85 87
86 static const G4double 88 static const G4double
87 al00 =-2.05398E+00, al01 = 2.38815E-02, a 89 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
88 al10 =-7.69748E-02, al11 =-6.91499E-02, a 90 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
89 al20 = 4.06463E-02, al21 =-1.01281E-02, a 91 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
90 92
91 static const G4double 93 static const G4double
92 bl00 = 1.04133E+00, bl01 =-9.43291E-03, b 94 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
93 bl10 = 1.19253E-01, bl11 = 4.07467E-02, b 95 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
94 bl20 =-1.59391E-02, bl21 = 7.27752E-03, b 96 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
95 97
96 using namespace std; 98 using namespace std;
97 99
98 G4eBremParametrizedModel::G4eBremParametrizedM 100 G4eBremParametrizedModel::G4eBremParametrizedModel(const G4ParticleDefinition* p,
99 const G4String& nam) 101 const G4String& nam)
100 : G4VEmModel(nam), 102 : G4VEmModel(nam),
101 particle(nullptr), << 103 particle(0),
>> 104 isElectron(true),
102 fMigdalConstant(classic_electr_radius*elec 105 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
103 bremFactor(fine_structure_const*classic_el 106 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
104 isInitialised(false), << 107 isInitialised(false)
105 isElectron(true) <<
106 { 108 {
107 theGamma = G4Gamma::Gamma(); 109 theGamma = G4Gamma::Gamma();
108 110
109 minThreshold = 0.1*keV; 111 minThreshold = 0.1*keV;
110 lowKinEnergy = 10.*MeV; 112 lowKinEnergy = 10.*MeV;
111 SetLowEnergyLimit(lowKinEnergy); 113 SetLowEnergyLimit(lowKinEnergy);
112 114
113 nist = G4NistManager::Instance(); 115 nist = G4NistManager::Instance();
114 116
115 SetAngularDistribution(new G4ModifiedTsai()) 117 SetAngularDistribution(new G4ModifiedTsai());
116 118
117 particleMass = kinEnergy = totalEnergy = cur 119 particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel = Finel
118 = densityFactor = densityCorr = fMax = fCo 120 = densityFactor = densityCorr = fMax = fCoulomb = 0.;
119 121
120 InitialiseConstants(); 122 InitialiseConstants();
121 if(nullptr != p) { SetParticle(p); } << 123 if(p) { SetParticle(p); }
122 } 124 }
123 125
124 //....oooOO0OOooo........oooOO0OOooo........oo 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 127
126 void G4eBremParametrizedModel::InitialiseConst 128 void G4eBremParametrizedModel::InitialiseConstants()
127 { 129 {
128 facFel = G4Log(184.15); 130 facFel = G4Log(184.15);
129 facFinel = G4Log(1194.); 131 facFinel = G4Log(1194.);
130 } 132 }
131 133
132 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 135
134 G4eBremParametrizedModel::~G4eBremParametrized << 136 G4eBremParametrizedModel::~G4eBremParametrizedModel()
>> 137 {
>> 138 }
135 139
136 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 141
138 void G4eBremParametrizedModel::SetParticle(con 142 void G4eBremParametrizedModel::SetParticle(const G4ParticleDefinition* p)
139 { 143 {
140 particle = p; 144 particle = p;
141 particleMass = p->GetPDGMass(); 145 particleMass = p->GetPDGMass();
142 if(p == G4Electron::Electron()) { isElectron 146 if(p == G4Electron::Electron()) { isElectron = true; }
143 else { isElectron 147 else { isElectron = false;}
144 } 148 }
145 149
146 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147 151
148 G4double G4eBremParametrizedModel::MinEnergyCu 152 G4double G4eBremParametrizedModel::MinEnergyCut(const G4ParticleDefinition*,
149 const G4MaterialCutsCouple*) 153 const G4MaterialCutsCouple*)
150 { 154 {
151 return minThreshold; 155 return minThreshold;
152 } 156 }
153 157
154 //....oooOO0OOooo........oooOO0OOooo........oo 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 159
156 void G4eBremParametrizedModel::SetupForMateria 160 void G4eBremParametrizedModel::SetupForMaterial(const G4ParticleDefinition*,
157 const G4Material* mat, 161 const G4Material* mat,
158 G4double kineticEnergy) 162 G4double kineticEnergy)
159 { 163 {
160 densityFactor = mat->GetElectronDensity()*fM 164 densityFactor = mat->GetElectronDensity()*fMigdalConstant;
161 165
162 // calculate threshold for density effect 166 // calculate threshold for density effect
163 kinEnergy = kineticEnergy; 167 kinEnergy = kineticEnergy;
164 totalEnergy = kineticEnergy + particleMass; 168 totalEnergy = kineticEnergy + particleMass;
165 densityCorr = densityFactor*totalEnergy*tota 169 densityCorr = densityFactor*totalEnergy*totalEnergy;
166 } 170 }
167 171
168 172
169 //....oooOO0OOooo........oooOO0OOooo........oo 173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170 174
171 void G4eBremParametrizedModel::Initialise(cons 175 void G4eBremParametrizedModel::Initialise(const G4ParticleDefinition* p,
172 const G4DataVector& cuts) 176 const G4DataVector& cuts)
173 { 177 {
174 if(p) { SetParticle(p); } 178 if(p) { SetParticle(p); }
175 179
176 lowKinEnergy = LowEnergyLimit(); 180 lowKinEnergy = LowEnergyLimit();
177 181
178 currentZ = 0.; 182 currentZ = 0.;
179 183
180 if(IsMaster()) { InitialiseElementSelectors( 184 if(IsMaster()) { InitialiseElementSelectors(p, cuts); }
181 185
182 if(isInitialised) { return; } 186 if(isInitialised) { return; }
183 fParticleChange = GetParticleChangeForLoss() 187 fParticleChange = GetParticleChangeForLoss();
184 isInitialised = true; 188 isInitialised = true;
185 } 189 }
186 190
187 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188 192
189 void G4eBremParametrizedModel::InitialiseLocal 193 void G4eBremParametrizedModel::InitialiseLocal(const G4ParticleDefinition*,
190 G4VEmModel* masterModel) 194 G4VEmModel* masterModel)
191 { 195 {
192 SetElementSelectors(masterModel->GetElementS 196 SetElementSelectors(masterModel->GetElementSelectors());
193 } 197 }
194 198
195 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 200
197 G4double G4eBremParametrizedModel::ComputeDEDX 201 G4double G4eBremParametrizedModel::ComputeDEDXPerVolume(
198 const G4Material* material, 202 const G4Material* material,
199 c 203 const G4ParticleDefinition* p,
200 G4double kineticEnergy, 204 G4double kineticEnergy,
201 G4double cutEnergy) 205 G4double cutEnergy)
202 { 206 {
203 if(!particle) { SetParticle(p); } 207 if(!particle) { SetParticle(p); }
204 if(kineticEnergy < lowKinEnergy) { return 0. 208 if(kineticEnergy < lowKinEnergy) { return 0.0; }
205 G4double cut = std::min(cutEnergy, kineticEn 209 G4double cut = std::min(cutEnergy, kineticEnergy);
206 if(cut == 0.0) { return 0.0; } 210 if(cut == 0.0) { return 0.0; }
207 211
208 SetupForMaterial(particle, material,kineticE 212 SetupForMaterial(particle, material,kineticEnergy);
209 213
210 const G4ElementVector* theElementVector = ma 214 const G4ElementVector* theElementVector = material->GetElementVector();
211 const G4double* theAtomicNumDensityVector = 215 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
212 216
213 G4double dedx = 0.0; 217 G4double dedx = 0.0;
214 218
215 // loop for elements in the material 219 // loop for elements in the material
216 for (size_t i=0; i<material->GetNumberOfElem 220 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
217 221
218 G4VEmModel::SetCurrentElement((*theElement 222 G4VEmModel::SetCurrentElement((*theElementVector)[i]);
219 SetCurrentElement((*theElementVector)[i]-> 223 SetCurrentElement((*theElementVector)[i]->GetZ());
220 224
221 dedx += theAtomicNumDensityVector[i]*curre 225 dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
222 } 226 }
223 dedx *= bremFactor; 227 dedx *= bremFactor;
224 228
225 return dedx; 229 return dedx;
226 } 230 }
227 231
228 //....oooOO0OOooo........oooOO0OOooo........oo 232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 233
230 G4double G4eBremParametrizedModel::ComputeBrem 234 G4double G4eBremParametrizedModel::ComputeBremLoss(G4double cut)
231 { 235 {
232 G4double loss = 0.0; 236 G4double loss = 0.0;
233 237
234 // number of intervals and integration step 238 // number of intervals and integration step
235 G4double vcut = cut/totalEnergy; 239 G4double vcut = cut/totalEnergy;
236 G4int n = (G4int)(20*vcut) + 3; 240 G4int n = (G4int)(20*vcut) + 3;
237 G4double delta = vcut/G4double(n); 241 G4double delta = vcut/G4double(n);
238 242
239 G4double e0 = 0.0; 243 G4double e0 = 0.0;
240 G4double xs; 244 G4double xs;
241 245
242 // integration 246 // integration
243 for(G4int l=0; l<n; l++) { 247 for(G4int l=0; l<n; l++) {
244 248
245 for(G4int i=0; i<8; i++) { 249 for(G4int i=0; i<8; i++) {
246 250
247 G4double eg = (e0 + xgi[i]*delta)*totalE 251 G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
248 252
249 xs = ComputeDXSectionPerAtom(eg); 253 xs = ComputeDXSectionPerAtom(eg);
250 254
251 loss += wgi[i]*xs/(1.0 + densityCorr/(eg 255 loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
252 } 256 }
253 e0 += delta; 257 e0 += delta;
254 } 258 }
255 259
256 loss *= delta*totalEnergy; 260 loss *= delta*totalEnergy;
257 261
258 return loss; 262 return loss;
259 } 263 }
260 264
261 //....oooOO0OOooo........oooOO0OOooo........oo 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 266
263 G4double G4eBremParametrizedModel::ComputeCros 267 G4double G4eBremParametrizedModel::ComputeCrossSectionPerAtom(
264 268 const G4ParticleDefinition* p,
265 G4double kineticEnergy, 269 G4double kineticEnergy,
266 G4double Z, G4double, 270 G4double Z, G4double,
267 G4double cutEnergy, 271 G4double cutEnergy,
268 G4double maxEnergy) 272 G4double maxEnergy)
269 { 273 {
270 if(!particle) { SetParticle(p); } 274 if(!particle) { SetParticle(p); }
271 if(kineticEnergy < lowKinEnergy) { return 0. 275 if(kineticEnergy < lowKinEnergy) { return 0.0; }
272 G4double cut = std::min(cutEnergy, kineticE 276 G4double cut = std::min(cutEnergy, kineticEnergy);
273 G4double tmax = std::min(maxEnergy, kineticE 277 G4double tmax = std::min(maxEnergy, kineticEnergy);
274 278
275 if(cut >= tmax) { return 0.0; } 279 if(cut >= tmax) { return 0.0; }
276 280
277 SetCurrentElement(Z); 281 SetCurrentElement(Z);
278 282
279 G4double cross = ComputeXSectionPerAtom(cut) 283 G4double cross = ComputeXSectionPerAtom(cut);
280 284
281 // allow partial integration 285 // allow partial integration
282 if(tmax < kinEnergy) { cross -= ComputeXSect 286 if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
283 287
284 cross *= Z*Z*bremFactor; 288 cross *= Z*Z*bremFactor;
285 289
286 return cross; 290 return cross;
287 } 291 }
288 292
289 //....oooOO0OOooo........oooOO0OOooo........oo 293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 294
291 G4double G4eBremParametrizedModel::ComputeXSec 295 G4double G4eBremParametrizedModel::ComputeXSectionPerAtom(G4double cut)
292 { 296 {
293 G4double cross = 0.0; 297 G4double cross = 0.0;
294 298
295 // number of intervals and integration step 299 // number of intervals and integration step
296 G4double vcut = G4Log(cut/totalEnergy); 300 G4double vcut = G4Log(cut/totalEnergy);
297 G4double vmax = G4Log(kinEnergy/totalEnergy) 301 G4double vmax = G4Log(kinEnergy/totalEnergy);
298 G4int n = (G4int)(0.45*(vmax - vcut)) + 4; 302 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
299 // n=1; // integration test 303 // n=1; // integration test
300 G4double delta = (vmax - vcut)/G4double(n); 304 G4double delta = (vmax - vcut)/G4double(n);
301 305
302 G4double e0 = vcut; 306 G4double e0 = vcut;
303 G4double xs; 307 G4double xs;
304 308
305 // integration 309 // integration
306 for(G4int l=0; l<n; l++) { 310 for(G4int l=0; l<n; l++) {
307 311
308 for(G4int i=0; i<8; i++) { 312 for(G4int i=0; i<8; i++) {
309 313
310 G4double eg = G4Exp(e0 + xgi[i]*delta)*t 314 G4double eg = G4Exp(e0 + xgi[i]*delta)*totalEnergy;
311 315
312 xs = ComputeDXSectionPerAtom(eg); 316 xs = ComputeDXSectionPerAtom(eg);
313 317
314 cross += wgi[i]*xs/(1.0 + densityCorr/(e 318 cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
315 } 319 }
316 e0 += delta; 320 e0 += delta;
317 } 321 }
318 322
319 cross *= delta; 323 cross *= delta;
320 324
321 return cross; 325 return cross;
322 } 326 }
323 327
324 //....oooOO0OOooo........oooOO0OOooo........oo 328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325 329
326 // compute the value of the screening function 330 // compute the value of the screening function 3*PHI1 - PHI2
327 331
328 G4double G4eBremParametrizedModel::ScreenFunct 332 G4double G4eBremParametrizedModel::ScreenFunction1(G4double ScreenVariable)
329 { 333 {
330 G4double screenVal; 334 G4double screenVal;
331 335
332 if (ScreenVariable > 1.) 336 if (ScreenVariable > 1.)
333 screenVal = 42.24 - 8.368*G4Log(ScreenVari 337 screenVal = 42.24 - 8.368*G4Log(ScreenVariable+0.952);
334 else 338 else
335 screenVal = 42.392 - ScreenVariable* (7.79 339 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
336 340
337 return screenVal; 341 return screenVal;
338 } 342 }
339 343
340 //....oooOO0OOooo........oooOO0OOooo........oo 344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
341 345
342 // compute the value of the screening function 346 // compute the value of the screening function 1.5*PHI1 - 0.5*PHI2
343 347
344 G4double G4eBremParametrizedModel::ScreenFunct 348 G4double G4eBremParametrizedModel::ScreenFunction2(G4double ScreenVariable)
345 { 349 {
346 G4double screenVal; 350 G4double screenVal;
347 351
348 if (ScreenVariable > 1.) 352 if (ScreenVariable > 1.)
349 screenVal = 42.24 - 8.368*G4Log(ScreenVari 353 screenVal = 42.24 - 8.368*G4Log(ScreenVariable+0.952);
350 else 354 else
351 screenVal = 41.734 - ScreenVariable* (6.48 355 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
352 356
353 return screenVal; 357 return screenVal;
354 } 358 }
355 359
356 //....oooOO0OOooo........oooOO0OOooo........oo 360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357 361
358 // Parametrized cross section 362 // Parametrized cross section
359 G4double G4eBremParametrizedModel::ComputePara 363 G4double G4eBremParametrizedModel::ComputeParametrizedDXSectionPerAtom(
360 G4double ki 364 G4double kineticEnergy,
361 G4double gammaEnergy, G4double Z) 365 G4double gammaEnergy, G4double Z)
362 { 366 {
363 SetCurrentElement(Z); 367 SetCurrentElement(Z);
364 G4double FZ = lnZ* (4.- 0.55*lnZ); 368 G4double FZ = lnZ* (4.- 0.55*lnZ);
365 G4double Z3 = z13; 369 G4double Z3 = z13;
366 G4double ZZ = z13*nist->GetZ13(G4lrint(Z)+1) 370 G4double ZZ = z13*nist->GetZ13(G4lrint(Z)+1);
367 371
368 totalEnergy = kineticEnergy + electron_mass_ 372 totalEnergy = kineticEnergy + electron_mass_c2;
369 373
370 // G4double x, epsil, greject, migdal, grej 374 // G4double x, epsil, greject, migdal, grejmax, q;
371 G4double epsil, greject; 375 G4double epsil, greject;
372 G4double U = G4Log(kineticEnergy/electron_m 376 G4double U = G4Log(kineticEnergy/electron_mass_c2);
373 G4double U2 = U*U; 377 G4double U2 = U*U;
374 378
375 // precalculated parameters 379 // precalculated parameters
376 G4double ah, bh; 380 G4double ah, bh;
377 381
378 if (kineticEnergy > tlow) { 382 if (kineticEnergy > tlow) {
379 383
380 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12 384 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
381 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22 385 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
382 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32 386 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
383 387
384 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12 388 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
385 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22 389 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
386 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32 390 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
387 391
388 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U 392 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
389 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U 393 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
390 394
391 // limit of the screening variable 395 // limit of the screening variable
392 G4double screenfac = 396 G4double screenfac =
393 136.*electron_mass_c2/(Z3*totalEnergy); 397 136.*electron_mass_c2/(Z3*totalEnergy);
394 398
395 epsil = gammaEnergy/totalEnergy; // 399 epsil = gammaEnergy/totalEnergy; // epsil = x*kineticEnergy/totalEnergy;
396 G4double screenvar = screenfac*epsil/( 400 G4double screenvar = screenfac*epsil/(1.0-epsil);
397 G4double F1 = max(ScreenFunction1(scre 401 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
398 G4double F2 = max(ScreenFunction2(scre 402 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
399 403
400 404
401 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2) 405 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.; // 1./(42.392 - FZ);
402 406
403 std::cout << " yy = "<<epsil<<std::endl; 407 std::cout << " yy = "<<epsil<<std::endl;
404 std::cout << " F1/(...) "<<F1/(42.392 - FZ 408 std::cout << " F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
405 std::cout << " F2/(...) "<<F2/(42.392 - FZ 409 std::cout << " F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
406 std::cout << " (42.392 - FZ) " << (42.392 410 std::cout << " (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
407 411
408 } else { 412 } else {
409 413
410 G4double al0 = al00 + ZZ* (al01 + ZZ* al02 414 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
411 G4double al1 = al10 + ZZ* (al11 + ZZ* al12 415 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
412 G4double al2 = al20 + ZZ* (al21 + ZZ* al22 416 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
413 417
414 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02 418 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
415 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12 419 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
416 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22 420 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
417 421
418 ah = al0 + al1*U + al2*U2; 422 ah = al0 + al1*U + al2*U2;
419 bh = bl0 + bl1*U + bl2*U2; 423 bh = bl0 + bl1*U + bl2*U2;
420 424
421 G4double x=gammaEnergy/kineticEnergy; 425 G4double x=gammaEnergy/kineticEnergy;
422 greject=(1. + x* (ah + bh*x)); 426 greject=(1. + x* (ah + bh*x));
423 427
424 /* 428 /*
425 // Compute the maximum of the rejection fu 429 // Compute the maximum of the rejection function
426 grejmax = max(1. + xmin* (ah + bh*xmin), 1 430 grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
427 G4double xm = -ah/(2.*bh); 431 G4double xm = -ah/(2.*bh);
428 if ( xmin < xm && xm < xmax) grejmax = max 432 if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
429 */ 433 */
430 } 434 }
431 435
432 return greject; 436 return greject;
433 } 437 }
434 438
435 //....oooOO0OOooo........oooOO0OOooo........oo 439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436 440
437 G4double G4eBremParametrizedModel::ComputeDXSe 441 G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
438 { 442 {
439 443
440 if(gammaEnergy < 0.0) { return 0.0; } 444 if(gammaEnergy < 0.0) { return 0.0; }
441 445
442 G4double y = gammaEnergy/totalEnergy; 446 G4double y = gammaEnergy/totalEnergy;
443 447
444 G4double main=0.; 448 G4double main=0.;
445 //secondTerm=0.; 449 //secondTerm=0.;
446 450
447 // ** form factors complete screening case * 451 // ** form factors complete screening case **
448 // only valid for high energies (and if LPM 452 // only valid for high energies (and if LPM suppression does not play a role)
449 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoul 453 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
450 // secondTerm = (1.-y)/12.*(1.+1./currentZ) 454 // secondTerm = (1.-y)/12.*(1.+1./currentZ);
451 455
452 std::cout<<" F1(0) "<<ScreenFunction1(0.) << 456 std::cout<<" F1(0) "<<ScreenFunction1(0.) <<std::endl;
453 std::cout<<" F1(0) "<<ScreenFunction2(0.) << 457 std::cout<<" F1(0) "<<ScreenFunction2(0.) <<std::endl;
454 std::cout<<"Ekin = "<<kinEnergy<<std::endl; 458 std::cout<<"Ekin = "<<kinEnergy<<std::endl;
455 std::cout<<"Z = "<<currentZ<<std::endl; 459 std::cout<<"Z = "<<currentZ<<std::endl;
456 std::cout<<"main = "<<main<<std::endl; 460 std::cout<<"main = "<<main<<std::endl;
457 std::cout<<" y = "<<y<<std::endl; 461 std::cout<<" y = "<<y<<std::endl;
458 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) 462 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) <<std::endl;
459 463
460 G4double main2 = ComputeParametrizedDXSectio 464 G4double main2 = ComputeParametrizedDXSectionPerAtom(kinEnergy,gammaEnergy,currentZ);
461 std::cout<<"main2 = "<<main2<<std::endl; 465 std::cout<<"main2 = "<<main2<<std::endl;
462 std::cout<<"main2tot = "<<main2 * ( (Fel-fCo 466 std::cout<<"main2tot = "<<main2 * ( (Fel-fCoulomb) + Finel/currentZ )/(Fel-fCoulomb);
463 467
464 G4double cross = main2; //main+secondTerm; 468 G4double cross = main2; //main+secondTerm;
465 return cross; 469 return cross;
466 } 470 }
467 471
468 //....oooOO0OOooo........oooOO0OOooo........oo 472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469 473
470 void G4eBremParametrizedModel::SampleSecondari 474 void G4eBremParametrizedModel::SampleSecondaries(
471 std::vector<G4DynamicParticle*>* 475 std::vector<G4DynamicParticle*>* vdp,
472 const G4MaterialCutsCouple* coup 476 const G4MaterialCutsCouple* couple,
473 const G4DynamicParticle* dp, 477 const G4DynamicParticle* dp,
474 G4double cutEnergy, 478 G4double cutEnergy,
475 G4double maxEnergy) 479 G4double maxEnergy)
476 { 480 {
477 G4double kineticEnergy = dp->GetKineticEnerg 481 G4double kineticEnergy = dp->GetKineticEnergy();
478 if(kineticEnergy < lowKinEnergy) { return; } 482 if(kineticEnergy < lowKinEnergy) { return; }
479 G4double cut = std::min(cutEnergy, kineticE 483 G4double cut = std::min(cutEnergy, kineticEnergy);
480 G4double emax = std::min(maxEnergy, kineticE 484 G4double emax = std::min(maxEnergy, kineticEnergy);
481 if(cut >= emax) { return; } 485 if(cut >= emax) { return; }
482 486
483 SetupForMaterial(particle, couple->GetMateri 487 SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
484 488
485 const G4Element* elm = SelectTargetAtom(coup << 489 const G4Element* elm =
486 dp-> << 490 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
487 SetCurrentElement(elm->GetZ()); 491 SetCurrentElement(elm->GetZ());
488 492
489 kinEnergy = kineticEnergy; 493 kinEnergy = kineticEnergy;
490 totalEnergy = kineticEnergy + particleMass; 494 totalEnergy = kineticEnergy + particleMass;
491 densityCorr = densityFactor*totalEnergy*tota 495 densityCorr = densityFactor*totalEnergy*totalEnergy;
492 496
493 G4double xmin = G4Log(cut*cut + densityCorr) 497 G4double xmin = G4Log(cut*cut + densityCorr);
494 G4double xmax = G4Log(emax*emax + densityCo 498 G4double xmax = G4Log(emax*emax + densityCorr);
495 G4double gammaEnergy, f, x; 499 G4double gammaEnergy, f, x;
496 500
497 CLHEP::HepRandomEngine* rndmEngine = G4Rando 501 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
498 502
499 do { 503 do {
500 x = G4Exp(xmin + rndmEngine->flat()*(xmax 504 x = G4Exp(xmin + rndmEngine->flat()*(xmax - xmin)) - densityCorr;
501 if(x < 0.0) x = 0.0; 505 if(x < 0.0) x = 0.0;
502 gammaEnergy = sqrt(x); 506 gammaEnergy = sqrt(x);
503 f = ComputeDXSectionPerAtom(gammaEnergy); 507 f = ComputeDXSectionPerAtom(gammaEnergy);
504 508
505 if ( f > fMax ) { 509 if ( f > fMax ) {
506 G4cout << "### G4eBremParametrizedModel 510 G4cout << "### G4eBremParametrizedModel Warning: Majoranta exceeded! "
507 << f << " > " << fMax 511 << f << " > " << fMax
508 << " Egamma(MeV)= " << gammaEnergy 512 << " Egamma(MeV)= " << gammaEnergy
509 << " E(mEV)= " << kineticEnergy 513 << " E(mEV)= " << kineticEnergy
510 << G4endl; 514 << G4endl;
511 } 515 }
512 516
513 // Loop checking, 03-Aug-2015, Vladimir Iv 517 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
514 } while (f < fMax*rndmEngine->flat()); 518 } while (f < fMax*rndmEngine->flat());
515 519
516 // 520 //
517 // angles of the emitted gamma. ( Z - axis a 521 // angles of the emitted gamma. ( Z - axis along the parent particle)
518 // use general interface 522 // use general interface
519 // 523 //
520 G4ThreeVector gammaDirection = 524 G4ThreeVector gammaDirection =
521 GetAngularDistribution()->SampleDirection( 525 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
522 G4lrint(currentZ), 526 G4lrint(currentZ),
523 couple->GetMaterial()); 527 couple->GetMaterial());
524 528
525 // create G4DynamicParticle object for the G 529 // create G4DynamicParticle object for the Gamma
526 auto gamma = new G4DynamicParticle(theGamma, << 530 G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
>> 531 gammaEnergy);
527 vdp->push_back(gamma); 532 vdp->push_back(gamma);
528 533
529 G4double totMomentum = sqrt(kineticEnergy*(t 534 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
530 G4ThreeVector direction = (totMomentum*dp->G 535 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
531 - gammaEnergy*gammaDirection).unit( 536 - gammaEnergy*gammaDirection).unit();
532 537
533 // energy of primary 538 // energy of primary
534 G4double finalE = kineticEnergy - gammaEnerg 539 G4double finalE = kineticEnergy - gammaEnergy;
535 540
536 // stop tracking and create new secondary in 541 // stop tracking and create new secondary instead of primary
537 if(gammaEnergy > SecondaryThreshold()) { 542 if(gammaEnergy > SecondaryThreshold()) {
538 fParticleChange->ProposeTrackStatus(fStopA 543 fParticleChange->ProposeTrackStatus(fStopAndKill);
539 fParticleChange->SetProposedKineticEnergy( 544 fParticleChange->SetProposedKineticEnergy(0.0);
540 auto el = << 545 G4DynamicParticle* el =
541 new G4DynamicParticle(const_cast<G4Parti 546 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
542 direction, finalE); 547 direction, finalE);
543 vdp->push_back(el); 548 vdp->push_back(el);
544 549
545 // continue tracking 550 // continue tracking
546 } else { 551 } else {
547 fParticleChange->SetProposedMomentumDirect 552 fParticleChange->SetProposedMomentumDirection(direction);
548 fParticleChange->SetProposedKineticEnergy( 553 fParticleChange->SetProposedKineticEnergy(finalE);
549 } 554 }
550 } 555 }
551 556
552 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 558
554 559
555 560