1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | #include "G4AdjointBremsstrahlungModel.hh" |
---|
27 | #include "G4AdjointCSManager.hh" |
---|
28 | #include "G4Integrator.hh" |
---|
29 | #include "G4TrackStatus.hh" |
---|
30 | #include "G4ParticleChange.hh" |
---|
31 | #include "G4AdjointElectron.hh" |
---|
32 | #include "G4Timer.hh" |
---|
33 | |
---|
34 | //////////////////////////////////////////////////////////////////////////////// |
---|
35 | // |
---|
36 | G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(): |
---|
37 | G4VEmAdjointModel("AdjointBremModel"), |
---|
38 | probsup(1.0), |
---|
39 | MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi), |
---|
40 | LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)), |
---|
41 | theLPMflag(true) |
---|
42 | |
---|
43 | { isElectron= true; |
---|
44 | SetUseMatrix(true); |
---|
45 | SetUseMatrixPerElement(false); |
---|
46 | SetApplyCutInRange(true); |
---|
47 | SetIsIonisation(false); |
---|
48 | highKinEnergy= 100.*TeV; |
---|
49 | lowKinEnergy = 1.0*keV; |
---|
50 | theTimer =new G4Timer(); |
---|
51 | |
---|
52 | theTimer->Start(); |
---|
53 | InitialiseParameters(); |
---|
54 | theTimer->Stop(); |
---|
55 | G4cout<<"Time elapsed in second for the initialidation of AdjointBrem "<<theTimer->GetRealElapsed()<<std::endl; |
---|
56 | |
---|
57 | ModeldCS="MODEL1"; |
---|
58 | |
---|
59 | } |
---|
60 | //////////////////////////////////////////////////////////////////////////////// |
---|
61 | // |
---|
62 | G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel() |
---|
63 | {;} |
---|
64 | //////////////////////////////////////////////////////////////////////////////// |
---|
65 | // |
---|
66 | /*G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond( |
---|
67 | const G4Material* aMaterial, |
---|
68 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
69 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
70 | ) |
---|
71 | { |
---|
72 | |
---|
73 | static const G4double |
---|
74 | ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02, |
---|
75 | ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02, |
---|
76 | ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02; |
---|
77 | |
---|
78 | static const G4double |
---|
79 | bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02, |
---|
80 | bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02, |
---|
81 | bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03; |
---|
82 | |
---|
83 | static const G4double |
---|
84 | al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04, |
---|
85 | al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03, |
---|
86 | al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04; |
---|
87 | |
---|
88 | static const G4double |
---|
89 | bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04, |
---|
90 | bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03, |
---|
91 | bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04; |
---|
92 | |
---|
93 | static const G4double tlow = 1.*MeV; |
---|
94 | |
---|
95 | G4double dCrossEprod=0.; |
---|
96 | G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); |
---|
97 | G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); |
---|
98 | |
---|
99 | |
---|
100 | if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ |
---|
101 | |
---|
102 | G4double cross = 0.0; |
---|
103 | |
---|
104 | |
---|
105 | |
---|
106 | G4double E1=kinEnergyProd; |
---|
107 | G4double E2=kinEnergyProd*1.000000001; |
---|
108 | G4double dE=(E2-E1); |
---|
109 | |
---|
110 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
111 | const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector(); |
---|
112 | G4double dum=0.; |
---|
113 | |
---|
114 | for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) { |
---|
115 | |
---|
116 | G4double fac= |
---|
117 | |
---|
118 | cross += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(), |
---|
119 | kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1); |
---|
120 | |
---|
121 | |
---|
122 | |
---|
123 | } |
---|
124 | dCrossEprod=(cross1-cross2)/dE; //first term |
---|
125 | |
---|
126 | //Now come the correction |
---|
127 | //----------------------- |
---|
128 | |
---|
129 | //First compute fsig for E1 |
---|
130 | //------------------------- |
---|
131 | |
---|
132 | |
---|
133 | G4double totalEnergy = kinEnergyProj+electron_mass_c2 ; |
---|
134 | G4double kp2 = MigdalConstant*totalEnergy*totalEnergy |
---|
135 | *(aMaterial->GetElectronDensity()); |
---|
136 | |
---|
137 | G4double fsig = 0.; |
---|
138 | G4int nmax = 100; |
---|
139 | G4double vmin=std::log(E1); |
---|
140 | G4double vmax=std::log(kinEnergyProj) ; |
---|
141 | G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin)); |
---|
142 | G4double u,fac,c,v,dv,y ; |
---|
143 | if(nn > 0) { |
---|
144 | |
---|
145 | dv = (vmax-vmin)/nn ; |
---|
146 | v = vmin-dv ; |
---|
147 | for(G4int n=0; n<=nn; n++) { |
---|
148 | |
---|
149 | v += dv; |
---|
150 | u = std::exp(v); |
---|
151 | fac = SupressionFunction(aMaterial, kinEnergyProj, u); |
---|
152 | y = u/kinEnergyProj; |
---|
153 | fac *= (4.-4.*y+3.*y*y)/3.; |
---|
154 | fac *= probsup*(u*u/(u*u+kp2))+1.-probsup; |
---|
155 | |
---|
156 | if ((n==0)||(n==nn)) c=0.5; |
---|
157 | else c=1. ; |
---|
158 | |
---|
159 | fac *= c; |
---|
160 | fsig += fac; |
---|
161 | } |
---|
162 | y = E1/kinEnergyProj ; |
---|
163 | fsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y)); |
---|
164 | |
---|
165 | } |
---|
166 | else { |
---|
167 | fsig = 1.; |
---|
168 | } |
---|
169 | if (fsig > 1.) fsig = 1.; |
---|
170 | |
---|
171 | dCrossEprod*=fsig; |
---|
172 | //return dCrossEprod; |
---|
173 | //Now we compute dfsig |
---|
174 | //------------------------- |
---|
175 | G4double dfsig = 0.; |
---|
176 | nn=20; |
---|
177 | vmax=std::log(E2) ; |
---|
178 | dv = (vmax-vmin)/nn ; |
---|
179 | v = vmin-dv ; |
---|
180 | for(G4int n=0; n<=nn; n++) { |
---|
181 | v += dv; |
---|
182 | u = std::exp(v); |
---|
183 | fac = SupressionFunction(aMaterial, kinEnergyProj, u); |
---|
184 | y = u/kinEnergyProj; |
---|
185 | fac *= (4.-4.*y+3.*y*y)/3.; |
---|
186 | fac *= probsup*(u*u/(u*u+kp2))+1.-probsup; |
---|
187 | |
---|
188 | if ((n==0)||(n==nn)) c=0.5; |
---|
189 | else c=1. ; |
---|
190 | |
---|
191 | fac *= c; |
---|
192 | dfsig += fac; |
---|
193 | } |
---|
194 | y = E1/kinEnergyProj; |
---|
195 | dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y)); |
---|
196 | |
---|
197 | dCrossEprod+=dfsig*cross1/dE; |
---|
198 | |
---|
199 | |
---|
200 | |
---|
201 | |
---|
202 | |
---|
203 | } |
---|
204 | return dCrossEprod; |
---|
205 | |
---|
206 | } |
---|
207 | */ |
---|
208 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial, |
---|
209 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
210 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
211 | ) |
---|
212 | {if (ModeldCS=="MODEL2") return DiffCrossSectionPerVolumePrimToSecond2(aMaterial, |
---|
213 | kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
214 | kinEnergyProd); |
---|
215 | if (ModeldCS=="MODEL3") return DiffCrossSectionPerVolumePrimToSecond3(aMaterial, |
---|
216 | kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
217 | kinEnergyProd); |
---|
218 | return DiffCrossSectionPerVolumePrimToSecond1(aMaterial, |
---|
219 | kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
220 | kinEnergyProd); |
---|
221 | } |
---|
222 | //////////////////////////////////////////////////////////////////////////////// |
---|
223 | // the one used till now |
---|
224 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond1( |
---|
225 | const G4Material* aMaterial, |
---|
226 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
227 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
228 | ) |
---|
229 | { |
---|
230 | G4double dCrossEprod=0.; |
---|
231 | G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); |
---|
232 | G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); |
---|
233 | |
---|
234 | |
---|
235 | if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ |
---|
236 | |
---|
237 | G4double cross1 = 0.0; |
---|
238 | G4double cross2 = 0.0; |
---|
239 | |
---|
240 | |
---|
241 | G4double E1=kinEnergyProd; |
---|
242 | G4double E2=kinEnergyProd*1.01; |
---|
243 | G4double dE=(E2-E1); |
---|
244 | |
---|
245 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
246 | const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector(); |
---|
247 | G4double dum=0.; |
---|
248 | |
---|
249 | for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) { |
---|
250 | |
---|
251 | cross1 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(), |
---|
252 | kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1); |
---|
253 | |
---|
254 | cross2 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(), |
---|
255 | kinEnergyProj, (*theElementVector)[i]->GetZ(), dum, E2); |
---|
256 | |
---|
257 | } |
---|
258 | dCrossEprod=(cross1-cross2)/dE; //first term |
---|
259 | |
---|
260 | //Now come the correction |
---|
261 | //----------------------- |
---|
262 | |
---|
263 | //First compute fsig for E1 |
---|
264 | //------------------------- |
---|
265 | |
---|
266 | |
---|
267 | G4double totalEnergy = kinEnergyProj+electron_mass_c2 ; |
---|
268 | G4double kp2 = MigdalConstant*totalEnergy*totalEnergy |
---|
269 | *(aMaterial->GetElectronDensity()); |
---|
270 | |
---|
271 | G4double fsig1 = 0.; |
---|
272 | G4int nmax = 100; |
---|
273 | G4double vmin=std::log(E1); |
---|
274 | G4double vmax=std::log(kinEnergyProj) ; |
---|
275 | G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin)); |
---|
276 | G4double u,fac,c,v,dv,y ; |
---|
277 | if(nn > 0) { |
---|
278 | |
---|
279 | dv = (vmax-vmin)/nn ; |
---|
280 | v = vmin-dv ; |
---|
281 | for(G4int n=0; n<=nn; n++) { |
---|
282 | |
---|
283 | v += dv; |
---|
284 | u = std::exp(v); |
---|
285 | fac = SupressionFunction(aMaterial, kinEnergyProj, u); |
---|
286 | y = u/kinEnergyProj; |
---|
287 | fac *= (4.-4.*y+3.*y*y)/3.; |
---|
288 | fac *= probsup*(u*u/(u*u+kp2))+1.-probsup; |
---|
289 | |
---|
290 | if ((n==0)||(n==nn)) c=0.5; |
---|
291 | else c=1. ; |
---|
292 | |
---|
293 | fac *= c; |
---|
294 | fsig1 += fac; |
---|
295 | } |
---|
296 | y = E1/kinEnergyProj ; |
---|
297 | fsig1 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y)); |
---|
298 | |
---|
299 | } |
---|
300 | else { |
---|
301 | fsig1 = 1.; |
---|
302 | } |
---|
303 | if (fsig1 > 1.) fsig1 = 1.; |
---|
304 | |
---|
305 | dCrossEprod*=fsig1; |
---|
306 | |
---|
307 | |
---|
308 | G4double fsig2 = 0.; |
---|
309 | vmin=std::log(E2); |
---|
310 | nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin)); |
---|
311 | if(nn > 0) { |
---|
312 | |
---|
313 | dv = (vmax-vmin)/nn ; |
---|
314 | v = vmin-dv ; |
---|
315 | for(G4int n=0; n<=nn; n++) { |
---|
316 | |
---|
317 | v += dv; |
---|
318 | u = std::exp(v); |
---|
319 | fac = SupressionFunction(aMaterial, kinEnergyProj, u); |
---|
320 | y = u/kinEnergyProj; |
---|
321 | fac *= (4.-4.*y+3.*y*y)/3.; |
---|
322 | fac *= probsup*(u*u/(u*u+kp2))+1.-probsup; |
---|
323 | |
---|
324 | if ((n==0)||(n==nn)) c=0.5; |
---|
325 | else c=1. ; |
---|
326 | |
---|
327 | fac *= c; |
---|
328 | fsig2 += fac; |
---|
329 | } |
---|
330 | y = E2/kinEnergyProj ; |
---|
331 | fsig2 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y)); |
---|
332 | |
---|
333 | } |
---|
334 | else { |
---|
335 | fsig2 = 1.; |
---|
336 | } |
---|
337 | if (fsig2 > 1.) fsig2 = 1.; |
---|
338 | |
---|
339 | |
---|
340 | G4double dfsig=(fsig2-fsig1); |
---|
341 | dCrossEprod+=dfsig*cross1/dE; |
---|
342 | |
---|
343 | dCrossEprod=(fsig1*cross1-fsig2*cross2)/dE; |
---|
344 | |
---|
345 | |
---|
346 | |
---|
347 | |
---|
348 | |
---|
349 | /*if (fsig < 1.){ |
---|
350 | //Now we compute dfsig |
---|
351 | //------------------------- |
---|
352 | G4double dfsig = 0.; |
---|
353 | nn=20; |
---|
354 | vmax=std::log(E2) ; |
---|
355 | dv = (vmax-vmin)/nn ; |
---|
356 | v = vmin-dv ; |
---|
357 | for(G4int n=0; n<=nn; n++) { |
---|
358 | v += dv; |
---|
359 | u = std::exp(v); |
---|
360 | fac = SupressionFunction(aMaterial, kinEnergyProj, u); |
---|
361 | y = u/kinEnergyProj; |
---|
362 | fac *= (4.-4.*y+3.*y*y)/3.; |
---|
363 | fac *= probsup*(u*u/(u*u+kp2))+1.-probsup; |
---|
364 | |
---|
365 | if ((n==0)||(n==nn)) c=0.5; |
---|
366 | else c=1. ; |
---|
367 | |
---|
368 | fac *= c; |
---|
369 | dfsig += fac; |
---|
370 | } |
---|
371 | y = E1/kinEnergyProj; |
---|
372 | dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y)); |
---|
373 | dCrossEprod+=dfsig*cross1/dE; |
---|
374 | |
---|
375 | } |
---|
376 | */ |
---|
377 | |
---|
378 | |
---|
379 | |
---|
380 | |
---|
381 | |
---|
382 | |
---|
383 | |
---|
384 | } |
---|
385 | return dCrossEprod; |
---|
386 | |
---|
387 | } |
---|
388 | |
---|
389 | |
---|
390 | //////////////////////////////////////////////////////////////////////////////// |
---|
391 | // |
---|
392 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2( |
---|
393 | const G4Material* aMaterial, |
---|
394 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
395 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
396 | ) |
---|
397 | { |
---|
398 | G4double dCrossEprod=0.; |
---|
399 | G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); |
---|
400 | G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); |
---|
401 | |
---|
402 | |
---|
403 | if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ |
---|
404 | |
---|
405 | G4double dEdX1 = 0.0; |
---|
406 | G4double dEdX2 = 0.0; |
---|
407 | |
---|
408 | |
---|
409 | G4double E1=kinEnergyProd; |
---|
410 | G4double E2=kinEnergyProd*1.001; |
---|
411 | G4double dE=(E2-E1); |
---|
412 | //G4double dum=0.; |
---|
413 | |
---|
414 | dEdX1 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E1); |
---|
415 | dEdX2 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E2); |
---|
416 | dCrossEprod=(dEdX2-dEdX1)/dE/E1; |
---|
417 | |
---|
418 | |
---|
419 | |
---|
420 | |
---|
421 | |
---|
422 | |
---|
423 | |
---|
424 | } |
---|
425 | return dCrossEprod; |
---|
426 | |
---|
427 | } |
---|
428 | //////////////////////////////////////////////////////////////////////////////// |
---|
429 | // |
---|
430 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond3( |
---|
431 | const G4Material* aMaterial, |
---|
432 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
433 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
434 | ) |
---|
435 | { |
---|
436 | |
---|
437 | return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial, |
---|
438 | kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
439 | kinEnergyProd); |
---|
440 | |
---|
441 | } |
---|
442 | |
---|
443 | //////////////////////////////////////////////////////////////////////////////// |
---|
444 | // |
---|
445 | G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material, |
---|
446 | G4double kineticEnergy, G4double gammaEnergy) |
---|
447 | { |
---|
448 | // supression due to the LPM effect+polarisation of the medium/ |
---|
449 | // supression due to the polarisation alone |
---|
450 | |
---|
451 | |
---|
452 | G4double totEnergy = kineticEnergy+electron_mass_c2 ; |
---|
453 | G4double totEnergySquare = totEnergy*totEnergy ; |
---|
454 | |
---|
455 | G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ; |
---|
456 | |
---|
457 | G4double gammaEnergySquare = gammaEnergy*gammaEnergy ; |
---|
458 | |
---|
459 | G4double electronDensity = material->GetElectronDensity(); |
---|
460 | |
---|
461 | G4double sp = gammaEnergySquare/ |
---|
462 | (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity); |
---|
463 | |
---|
464 | G4double supr = 1.0; |
---|
465 | |
---|
466 | if (theLPMflag) { |
---|
467 | |
---|
468 | G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare; |
---|
469 | |
---|
470 | if (s2lpm < 1.) { |
---|
471 | |
---|
472 | G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ; |
---|
473 | G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit; |
---|
474 | G4double splim = LPMgEnergyLimit2/ |
---|
475 | (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity); |
---|
476 | G4double w = 1.+1./splim ; |
---|
477 | |
---|
478 | if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp); |
---|
479 | else w = s2lpm*(1.+1./sp); |
---|
480 | |
---|
481 | supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ; |
---|
482 | supr /= sp; |
---|
483 | } |
---|
484 | |
---|
485 | } |
---|
486 | return supr; |
---|
487 | } |
---|
488 | |
---|
489 | //////////////////////////////////////////////////////////////////////////////// |
---|
490 | // |
---|
491 | void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack, |
---|
492 | G4bool IsScatProjToProjCase, |
---|
493 | G4ParticleChange* fParticleChange) |
---|
494 | { |
---|
495 | |
---|
496 | //G4cout<<"Adjoint Brem"<<std::endl; |
---|
497 | const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); |
---|
498 | |
---|
499 | size_t ind=0; |
---|
500 | |
---|
501 | if (UseMatrixPerElement ) { //Select Material |
---|
502 | std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase; |
---|
503 | if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase; |
---|
504 | G4double rand_var= G4UniformRand(); |
---|
505 | G4double SumCS=0.; |
---|
506 | for (size_t i=0;i<CS_Vs_Element->size();i++){ |
---|
507 | SumCS+=(*CS_Vs_Element)[i]; |
---|
508 | if (rand_var<=SumCS/lastCS){ |
---|
509 | ind=i; |
---|
510 | break; |
---|
511 | } |
---|
512 | } |
---|
513 | } |
---|
514 | else { |
---|
515 | ind = currentMaterialIndex; |
---|
516 | } |
---|
517 | |
---|
518 | |
---|
519 | //Elastic inverse scattering modified compared to general G4VEmAdjointModel |
---|
520 | //--------------------------- |
---|
521 | G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); |
---|
522 | G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); |
---|
523 | //G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum(); |
---|
524 | if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ |
---|
525 | return; |
---|
526 | } |
---|
527 | |
---|
528 | //Sample secondary energy |
---|
529 | //----------------------- |
---|
530 | |
---|
531 | G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind, |
---|
532 | adjointPrimKinEnergy, |
---|
533 | IsScatProjToProjCase); |
---|
534 | |
---|
535 | |
---|
536 | |
---|
537 | |
---|
538 | //Weight correction |
---|
539 | //----------------------- |
---|
540 | CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy); |
---|
541 | |
---|
542 | |
---|
543 | //Kinematic |
---|
544 | //--------- |
---|
545 | |
---|
546 | G4double projectileM0 = electron_mass_c2; |
---|
547 | G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; |
---|
548 | G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; |
---|
549 | G4double projectileP = std::sqrt(projectileP2); |
---|
550 | |
---|
551 | |
---|
552 | //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel |
---|
553 | //------------------------------------------------ |
---|
554 | G4double u; |
---|
555 | const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; |
---|
556 | |
---|
557 | if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; |
---|
558 | else u = - std::log(G4UniformRand()*G4UniformRand())/a2; |
---|
559 | |
---|
560 | G4double theta = u*electron_mass_c2/projectileTotalEnergy; |
---|
561 | |
---|
562 | G4double sint = std::sin(theta); |
---|
563 | G4double cost = std::cos(theta); |
---|
564 | |
---|
565 | G4double phi = twopi * G4UniformRand() ; |
---|
566 | |
---|
567 | G4ThreeVector projectileMomentum; |
---|
568 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame |
---|
569 | if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- |
---|
570 | G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); |
---|
571 | G4ThreeVector dirProd=projectileMomentum-gammaMomentum; |
---|
572 | G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); |
---|
573 | G4double sint1 = std::sqrt(1.-cost1*cost1); |
---|
574 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; |
---|
575 | |
---|
576 | } |
---|
577 | |
---|
578 | projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); |
---|
579 | |
---|
580 | |
---|
581 | |
---|
582 | if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary |
---|
583 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
584 | fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); |
---|
585 | //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl; |
---|
586 | } |
---|
587 | else { |
---|
588 | fParticleChange->ProposeEnergy(projectileKinEnergy); |
---|
589 | fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); |
---|
590 | //G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl; |
---|
591 | } |
---|
592 | } |
---|
593 | //////////////////////////////////////////////////////////////////////////////// |
---|
594 | // |
---|
595 | void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel) |
---|
596 | {theDirectBremModel=aModel; |
---|
597 | DefineDirectEMModel(aModel); |
---|
598 | } |
---|
599 | //////////////////////////////////////////////////////////////////////////////// |
---|
600 | // |
---|
601 | void G4AdjointBremsstrahlungModel::InitialiseParameters() |
---|
602 | { |
---|
603 | static const G4double |
---|
604 | ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02, |
---|
605 | ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02, |
---|
606 | ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02; |
---|
607 | |
---|
608 | static const G4double |
---|
609 | bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02, |
---|
610 | bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02, |
---|
611 | bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03; |
---|
612 | |
---|
613 | /* static const G4double |
---|
614 | al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04, |
---|
615 | al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03, |
---|
616 | al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04; |
---|
617 | |
---|
618 | static const G4double |
---|
619 | bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04, |
---|
620 | bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03, |
---|
621 | bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/ |
---|
622 | |
---|
623 | |
---|
624 | const G4ElementTable* theElementTable = G4Element::GetElementTable(); |
---|
625 | FZ.clear(); |
---|
626 | ah1.clear(); |
---|
627 | ah2.clear(); |
---|
628 | ah3.clear(); |
---|
629 | |
---|
630 | bh1.clear(); |
---|
631 | bh2.clear(); |
---|
632 | bh3.clear(); |
---|
633 | |
---|
634 | al0.clear(); |
---|
635 | al1.clear(); |
---|
636 | al2.clear(); |
---|
637 | |
---|
638 | bl0.clear(); |
---|
639 | bl1.clear(); |
---|
640 | bl2.clear(); |
---|
641 | SigmaPerAtom.clear(); |
---|
642 | |
---|
643 | for (size_t j=0; j<theElementTable->size();j++){ |
---|
644 | |
---|
645 | G4Element* anElement=(*theElementTable)[j]; |
---|
646 | G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3()); |
---|
647 | FZ.push_back(lnZ* (4.- 0.55*lnZ)); |
---|
648 | G4double ZZ = anElement->GetIonisation()->GetZZ3(); |
---|
649 | |
---|
650 | ah1.push_back(ah10 + ZZ* (ah11 + ZZ* ah12)); |
---|
651 | ah2.push_back(ah20 + ZZ* (ah21 + ZZ* ah22)); |
---|
652 | ah3.push_back(ah30 + ZZ* (ah31 + ZZ* ah32)); |
---|
653 | |
---|
654 | bh1.push_back(bh10 + ZZ* (bh11 + ZZ* bh12)); |
---|
655 | bh2.push_back(bh20 + ZZ* (bh21 + ZZ* bh22)); |
---|
656 | bh3.push_back(bh30 + ZZ* (bh31 + ZZ* bh32)); |
---|
657 | /*SigmaPerAtom.push_back(theDirectEMModel->ComputeCrossSectionPerAtom( |
---|
658 | theDirectPrimaryPartDef,GetHighEnergyLimit()/2., |
---|
659 | anElement->GetZ(),1.,GetLowEnergyLimit(),1.e20));*/ |
---|
660 | |
---|
661 | |
---|
662 | |
---|
663 | } |
---|
664 | } |
---|