| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4PhotoElectricAngularGenera
33 //
34 // Author: A. C. Farinha, L. Peralta, P. Rodri
35 //
36 // Creation date:
37 //
38 // Modifications:
39 // 10 January 2006
40 // 06 May 2011, Replace FILE with std::ifstrea
41 //
42 // Class Description:
43 //
44 // Concrete class for PhotoElectric Electron A
45 //
46 // Class Description:
47 // PhotoElectric Electron Angular Generator ba
48 // Includes polarization effects for K and L1
49 // For higher shells the L1 cross-section is u
50 //
51 // The Gavrila photoelectron angular distribut
52 // the inverse-transform method (James 1980).
53 // used to sample bremsstrahlung 2BN cross sec
54 //
55 // M. Gavrila, "Relativistic K-Shell Photoeffe
56 // M. Gavrila, "Relativistic L-Shell Photoeffe
57 // F. James, Rept. on Prog. in Phys. 43, 1145
58 // L. Peralta et al., "A new low-energy bremss
59 //
60 //
61 // -------------------------------------------
62 //
63 //
64
65 #include "G4PhotoElectricAngularGeneratorPolar
66 #include "G4PhysicalConstants.hh"
67 #include "G4RotationMatrix.hh"
68 #include "Randomize.hh"
69 #include "G4Exp.hh"
70
71 G4PhotoElectricAngularGeneratorPolarized::G4Ph
72 :G4VEmAngularDistribution("AngularGenSauterG
73 {
74 const G4int arrayDim = 980;
75
76 //minimum electron beta parameter allowed
77 betaArray[0] = 0.02;
78 //beta step
79 betaArray[1] = 0.001;
80 //maximum index array for a and c tables
81 betaArray[2] = arrayDim - 1;
82
83 // read Majorant Surface Parameters. This ar
84 //Gavrila angular photoelectron distribution
85 for(G4int level = 0; level < 2; level++){
86 char nameChar0[100] = "ftab0.dat"; // K-sh
87 char nameChar1[100] = "ftab1.dat"; // L-sh
88
89 G4String filename;
90 if(level == 0) filename = nameChar0;
91 if(level == 1) filename = nameChar1;
92
93 const char* path = G4FindDataDir("G4LEDATA
94 if (!path)
95 {
96 G4String excep = "G4EMDataSet - G4LEDA
97 G4Exception("G4PhotoElectricAngularGen
98 "em0006",FatalException,"G4LEDATA envi
99 return;
100 }
101
102 G4String pathString(path);
103 G4String dirFile = pathString + "/photoele
104 std::ifstream infile(dirFile);
105 if (!infile.is_open())
106 {
107 G4String excep = "data file: " + dirFile + "
108 G4Exception("G4PhotoElectricAngularGen
109 "em0003",FatalException,excep);
110 return;
111 }
112
113 // Read parameters into tables. The parame
114 // energy and shell level
115 G4float aRead=0,cRead=0, beta=0;
116 for(G4int i=0 ; i<arrayDim ;++i){
117 infile >> beta >> aRead >> cRead;
118 aMajorantSurfaceParameterTable[i][level]
119 cMajorantSurfaceParameterTable[i][level]
120 }
121 infile.close();
122 }
123 }
124
125 //....oooOO0OOooo........oooOO0OOooo........oo
126
127 G4PhotoElectricAngularGeneratorPolarized::~G4P
128 {}
129
130 //....oooOO0OOooo........oooOO0OOooo........oo
131
132 G4ThreeVector&
133 G4PhotoElectricAngularGeneratorPolarized::Samp
134 const G4Dynam
135 G4double eKinEnergy,
136 G4int shellId,
137 const G4Material*)
138 {
139 // (shellId == 0) - K-shell - Polarized mod
140 // (shellId > 0) - L1-shell - Polarized mod
141
142 // Calculate Lorentz term (gamma) and beta p
143 G4double gamma = 1 + eKinEnergy/electron_mas
144 G4double beta = std::sqrt((gamma - 1)*(gamm
145
146 const G4ThreeVector& direction = dp->GetMome
147 const G4ThreeVector& polarization = dp->GetP
148
149 G4double theta, phi = 0;
150 // Majorant surface parameters
151 // function of the outgoing electron kinetic
152 G4double aBeta = 0;
153 G4double cBeta = 0;
154
155 // For the outgoing kinetic energy
156 // find the current majorant surface paramet
157 PhotoElectronGetMajorantSurfaceAandCParamete
158
159 // Generate pho and theta according to the s
160 // and beta parameter of the electron
161 PhotoElectronGeneratePhiAndTheta(shellId, be
162
163 // Determine the rotation matrix
164 const G4RotationMatrix rotation =
165 PhotoElectronRotationMatrix(direction, pol
166
167 // Compute final direction of the outgoing e
168 fLocalDirection = PhotoElectronComputeFinalD
169
170 return fLocalDirection;
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oo
174
175 void
176 G4PhotoElectricAngularGeneratorPolarized::Phot
177 G4int shellLevel, G4double beta, G4doubl
178 G4double *pphi, G4double *ptheta) const
179 {
180 G4double rand1, rand2, rand3 = 0;
181 G4double phi = 0;
182 G4double theta = 0;
183 G4double crossSectionValue = 0;
184 G4double crossSectionMajorantFunctionValue =
185 G4double maxBeta = 0;
186
187 do {
188
189 rand1 = G4UniformRand();
190 rand2 = G4UniformRand();
191 rand3 = G4UniformRand();
192
193 phi=2*pi*rand1;
194
195 if(shellLevel == 0){
196 // Polarized Gavrila Cross-Section for K
197 theta=std::sqrt(((G4Exp(rand2*std::log(1
198 crossSectionMajorantFunctionValue =
199 CrossSectionMajorantFunction(theta, cBeta);
200 crossSectionValue = DSigmaKshellGavrila1
201 } else {
202 // Polarized Gavrila Cross-Section for
203 theta = std::sqrt(((G4Exp(rand2*std::log
204 crossSectionMajorantFunctionValue =
205 CrossSectionMajorantFunction(theta, cBeta);
206 crossSectionValue = DSigmaL1shellGavrila
207 }
208
209 maxBeta=rand3*aBeta*crossSectionMajorantFu
210 if(crossSectionValue < 0.0) { crossSection
211
212 } while(maxBeta > crossSectionValue || theta
213
214 *pphi = phi;
215 *ptheta = theta;
216 }
217
218 //....oooOO0OOooo........oooOO0OOooo........oo
219
220 G4double
221 G4PhotoElectricAngularGeneratorPolarized::Cros
222 G4double theta, G4double cBeta) const
223 {
224 // Compute Majorant Function
225 G4double crossSectionMajorantFunctionValue =
226 crossSectionMajorantFunctionValue = theta/(1
227 return crossSectionMajorantFunctionValue;
228 }
229
230 //....oooOO0OOooo........oooOO0OOooo........oo
231
232 G4double
233 G4PhotoElectricAngularGeneratorPolarized::DSig
234 G4double beta, G4double theta, G4doub
235 {
236 //Double differential K shell cross-section
237
238 G4double beta2 = beta*beta;
239 G4double oneBeta2 = 1 - beta2;
240 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
241 G4double oneBeta2_to_3_2 = std::pow(oneBeta2
242 G4double cosTheta = std::cos(theta);
243 G4double sinTheta2 = std::sin(theta)*std::si
244 G4double cosPhi2 = std::cos(phi)*std::cos(ph
245 G4double oneBetaCosTheta = 1-beta*cosTheta;
246 G4double dsigma = 0;
247 G4double firstTerm = 0;
248 G4double secondTerm = 0;
249
250 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe
251 (sinTheta2 * cosPhi2)/std::pow(o
252 (1-sqrtOneBeta2)/(4*oneBeta2_to_
253
254 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st
255 (4*beta2/sqrtOneBeta2 * sinThet
256 - 4*(1-sqrtOneBeta2)/oneBeta2 *
257 + 4*beta2*(1-sqrtOneBeta2)/oneB
258 + (1-sqrtOneBeta2)/(4*beta2*one
259 (1-sqrtOneBeta2)/oneBeta2_to_3_
260
261 dsigma = ( firstTerm*(1-pi*fine_structure_co
262
263 return dsigma;
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oo
267
268 G4double
269 G4PhotoElectricAngularGeneratorPolarized::DSig
270 G4double beta, G4double theta, G4doubl
271 {
272 //Double differential L1 shell cross-section
273 // 26Oct2022: included factor (1/8) in dsigm
274 G4double beta2 = beta*beta;
275 G4double oneBeta2 = 1-beta2;
276 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
277 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1
278 G4double cosTheta = std::cos(theta);
279 G4double sinTheta2 =std::sin(theta)*std::sin
280 G4double cosPhi2 = std::cos(phi)*std::cos(ph
281 G4double oneBetaCosTheta = 1-beta*cosTheta;
282
283 G4double dsigma = 0;
284 G4double firstTerm = 0;
285 G4double secondTerm = 0;
286
287 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe
288 * (sinTheta2 * cosPhi2)/std::po
289 (1-sqrtOneBeta2)/(4*oneBeta2_to_
290
291 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st
292 (4*beta2/sqrtOneBeta2 * sinThet
293 - 4*(1-sqrtOneBeta2)/oneBeta2 *
294 + 4*beta2*(1-sqrtOneBeta2)/oneB
295 + (1-sqrtOneBeta2)/(4*beta2*one
296 (1-sqrtOneBeta2)/oneBeta2_to_3_
297
298 dsigma = ( firstTerm*(1-pi*fine_structure_co
299
300 return dsigma;
301 }
302
303 //....oooOO0OOooo........oooOO0OOooo........oo
304
305 G4RotationMatrix
306 G4PhotoElectricAngularGeneratorPolarized::Phot
307 const G4ThreeVector& direction,
308 const G4ThreeVector& polarization)
309 {
310 G4double mK = direction.mag();
311 G4double mS = polarization.mag();
312 G4ThreeVector polarization2 = polarization;
313 const G4double kTolerance = 1e-6;
314
315 if(!(polarization.isOrthogonal(direction,kTo
316 G4ThreeVector d0 = direction.unit();
317 G4ThreeVector a1 = PerpendicularVector(d0)
318 G4ThreeVector a0 = a1.unit();
319 G4double rand1 = G4UniformRand();
320 G4double angle = twopi*rand1;
321 G4ThreeVector b0 = d0.cross(a0);
322 G4ThreeVector c;
323 c.setX(std::cos(angle)*(a0.x())+std::sin(a
324 c.setY(std::cos(angle)*(a0.y())+std::sin(a
325 c.setZ(std::cos(angle)*(a0.z())+std::sin(a
326 polarization2 = c.unit();
327 mS = polarization2.mag();
328 }else
329 {
330 if ( polarization.howOrthogonal(directio
331 {
332 polarization2 = polarization
333 - polarization.dot(direction)/direction.
334 }
335 }
336
337 G4ThreeVector direction2 = direction/mK;
338 polarization2 = polarization2/mS;
339
340 G4ThreeVector y = direction2.cross(polarizat
341
342 G4RotationMatrix R(polarization2,y,direction
343 return R;
344 }
345
346 //....oooOO0OOooo........oooOO0OOooo........oo
347
348 void
349 G4PhotoElectricAngularGeneratorPolarized::Phot
350 {
351 // This member function finds for a given sh
352 // of the outgoing electron the correct Majo
353 G4double aBeta,cBeta;
354 G4double bMin,bStep;
355 G4int indexMax;
356 G4int level = 0;
357 if(shellId > 0) { level = 1; }
358
359 bMin = betaArray[0];
360 bStep = betaArray[1];
361 indexMax = (G4int)betaArray[2];
362 const G4double kBias = 1e-9;
363
364 G4int k = (G4int)((beta-bMin+kBias)/bStep);
365
366 if(k < 0)
367 k = 0;
368 if(k > indexMax)
369 k = indexMax;
370
371 if(k == 0)
372 aBeta = std::max(aMajorantSurfaceParameter
373 else if(k==indexMax)
374 aBeta = std::max(aMajorantSurfaceParameter
375 else{
376 aBeta = std::max(aMajorantSurfaceParameter
377 aBeta = std::max(aBeta,aMajorantSurfacePar
378 }
379
380 if(k == 0)
381 cBeta = std::max(cMajorantSurfaceParameter
382 else if(k == indexMax)
383 cBeta = std::max(cMajorantSurfaceParameter
384 else{
385 cBeta = std::max(cMajorantSurfaceParameter
386 cBeta = std::max(cBeta,cMajorantSurfacePar
387 }
388
389 *majorantSurfaceParameterA = aBeta;
390 *majorantSurfaceParameterC = cBeta;
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oo
394
395 G4ThreeVector G4PhotoElectricAngularGeneratorP
396 {
397 //computes the photoelectron momentum unitar
398 G4double sint = std::sin(theta);
399 G4double px = std::cos(phi)*sint;
400 G4double py = std::sin(phi)*sint;
401 G4double pz = std::cos(theta);
402
403 G4ThreeVector samplingDirection(px,py,pz);
404
405 G4ThreeVector outgoingDirection = rotation*s
406 return outgoingDirection;
407 }
408
409 //....oooOO0OOooo........oooOO0OOooo........oo
410
411 void G4PhotoElectricAngularGeneratorPolarized:
412 {
413 G4cout << "\n" << G4endl;
414 G4cout << "Polarized Photoelectric Angular G
415 G4cout << "PhotoElectric Electron Angular Ge
416 G4cout << "Includes polarization effects for
417 G4cout << "For higher shells the L1 cross-se
418 G4cout << "(see Physics Reference Manual) \n
419 }
420
421 //....oooOO0OOooo........oooOO0OOooo........oo
422
423 G4ThreeVector
424 G4PhotoElectricAngularGeneratorPolarized::Perp
425 const G4ThreeVector& a) const
426 {
427 G4double dx = a.x();
428 G4double dy = a.y();
429 G4double dz = a.z();
430 G4double x = dx < 0.0 ? -dx : dx;
431 G4double y = dy < 0.0 ? -dy : dy;
432 G4double z = dz < 0.0 ? -dz : dz;
433 if (x < y) {
434 return x < z ? G4ThreeVector(-dy,dx,0) : G
435 }else{
436 return y < z ? G4ThreeVector(dz,0,-dx) : G
437 }
438 }
439