| 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: G4RDPhotoElectricAngularGene
33 //
34 // Author: A. C. Farinha, L. Peralta, P. Rodri
35 //
36 // Creation date:
37 //
38 // Modifications:
39 // 10 January 2006
40 //
41 // Class Description:
42 //
43 // Concrete class for PhotoElectric Electron A
44 //
45 // Class Description:
46 // PhotoElectric Electron Angular Generator ba
47 // Includes polarization effects for K and L1
48 // For higher shells the L1 cross-section is u
49 //
50 // The Gavrila photoelectron angular distribut
51 // the inverse-transform method (James 1980).
52 // used to sample bremsstrahlung 2BN cross sec
53 //
54 // M. Gavrila, "Relativistic K-Shell Photoeffe
55 // M. Gavrila, "Relativistic L-Shell Photoeffe
56 // F. James, Rept. on Prog. in Phys. 43, 1145
57 // L. Peralta et al., "A new low-energy bremss
58 //
59 //
60 // -------------------------------------------
61 //
62 //
63
64 #include "G4RDPhotoElectricAngularGeneratorPol
65 #include "G4RotationMatrix.hh"
66 #include "Randomize.hh"
67 #include "G4PhysicalConstants.hh"
68
69 //
70
71 G4RDPhotoElectricAngularGeneratorPolarized::G4
72 {
73 const G4int arrayDim = 980;
74
75 //minimum electron beta parameter allowed
76 betaArray[0] = 0.02;
77 //beta step
78 betaArray[1] = 0.001;
79 //maximum index array for a and c tables
80 betaArray[2] = arrayDim - 1;
81
82 // read Majorant Surface Parameters. This ar
83 for(G4int level = 0; level < 2; level++){
84
85 char nameChar0[100] = "ftab0.dat"; // K-sh
86 char nameChar1[100] = "ftab1.dat"; // L-sh
87
88 G4String filename;
89 if(level == 0) filename = nameChar0;
90 if(level == 1) filename = nameChar1;
91
92 const char* path = G4FindDataDir("G4LEDATA
93 if (!path)
94 {
95 G4String excep = "G4LEDATA environment
96 G4Exception("G4RDPhotoElectricAngularG
97 "InvalidSetup", FatalExcep
98 }
99
100 G4String pathString(path);
101 G4String dirFile = pathString + "/photoele
102 FILE *infile;
103 infile = fopen(dirFile,"r");
104 if (infile == 0)
105 {
106 G4String excep = "Data file: " + dirFile + "
107 G4Exception("G4RDPhotoElectricAngularGenerat
108 "DataNotFound", FatalExcep
109 }
110
111 // Read parameters into tables. The parame
112 G4float aRead,cRead, beta;
113 for(G4int i=0 ; i<arrayDim ;i++){
114 fscanf(infile,"%f\t %e\t %e",&beta,&aRea
115 aMajorantSurfaceParameterTable[i][level]
116 cMajorantSurfaceParameterTable[i][level]
117 }
118 fclose(infile);
119
120 }
121 }
122
123 //
124
125 G4RDPhotoElectricAngularGeneratorPolarized::~G
126 {;}
127
128 //
129
130 G4ThreeVector G4RDPhotoElectricAngularGenerato
131 const G4ThreeVector& pol
132 {
133 // Calculate Lorentz term (gamma) and beta p
134 G4double gamma = 1. + eKineticEnergy/elect
135 G4double beta = std::sqrt(gamma*gamma-1.)/g
136
137 G4double theta, phi = 0;
138 G4double aBeta = 0; // Majorant surface para
139 G4double cBeta = 0; // Majorant surface para
140
141 G4int shellLevel = 0;
142 if(shellId < 2) shellLevel = 0; // K-shell
143 if(shellId >= 2) shellLevel = 1; // L1-shell
144
145 // For the outgoing kinetic energy find the
146 PhotoElectronGetMajorantSurfaceAandCParamete
147
148 // Generate pho and theta according to the s
149 PhotoElectronGeneratePhiAndTheta(shellLevel,
150
151 // Determine the rotation matrix
152 G4RotationMatrix rotation = PhotoElectronRot
153
154 // Compute final direction of the outgoing e
155 G4ThreeVector final_direction = PhotoElectro
156
157 return final_direction;
158 }
159
160 //
161
162 void G4RDPhotoElectricAngularGeneratorPolarize
163 const G4double aBeta
164 G4double *pphi, G4do
165 {
166 G4double rand1, rand2, rand3 = 0;
167 G4double phi = 0;
168 G4double theta = 0;
169 G4double crossSectionValue = 0;
170 G4double crossSectionMajorantFunctionValue =
171 G4double maxBeta = 0;
172
173 do {
174
175 rand1 = G4UniformRand();
176 rand2 = G4UniformRand();
177 rand3 = G4UniformRand();
178
179 phi=2*pi*rand1;
180
181 if(shellLevel == 0){
182
183 // Polarized Gavrila Cross-Section for K
184 theta=std::sqrt(((std::exp(rand2*std::lo
185 crossSectionMajorantFunctionValue = Cros
186 crossSectionValue = DSigmaKshellGavrila1
187
188 } else {
189
190 // Polarized Gavrila Cross-Section for
191 theta = std::sqrt(((std::exp(rand2*std::
192 crossSectionMajorantFunctionValue = Cros
193 crossSectionValue = DSigmaL1shellGavrila
194
195 }
196
197 maxBeta=rand3*aBeta*crossSectionMajorantFu
198
199 }while(maxBeta > crossSectionValue);
200
201 *pphi = phi;
202 *ptheta = theta;
203 }
204
205 //
206
207 G4double G4RDPhotoElectricAngularGeneratorPola
208 {
209 // Compute Majorant Function
210 G4double crossSectionMajorantFunctionValue =
211 crossSectionMajorantFunctionValue = theta/(1
212 return crossSectionMajorantFunctionValue;
213 }
214
215 //
216
217 G4double G4RDPhotoElectricAngularGeneratorPola
218 {
219
220 //Double differential K shell cross-section
221
222 G4double beta2 = beta*beta;
223 G4double oneBeta2 = 1 - beta2;
224 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
225 G4double oneBeta2_to_3_2 = std::pow(oneBeta2
226 G4double cosTheta = std::cos(theta);
227 G4double sinTheta2 = std::sin(theta)*std::si
228 G4double cosPhi2 = std::cos(phi)*std::cos(ph
229 G4double oneBetaCosTheta = 1-beta*cosTheta;
230 G4double dsigma = 0;
231 G4double firstTerm = 0;
232 G4double secondTerm = 0;
233
234 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe
235 (sinTheta2 * cosPhi2)/std::pow(o
236 (1-sqrtOneBeta2)/(4*oneBeta2_to_
237
238 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st
239 (4*beta2/sqrtOneBeta2 * sinThet
240 - 4*(1-sqrtOneBeta2)/oneBeta2 *
241 + 4*beta2*(1-sqrtOneBeta2)/oneB
242 + (1-sqrtOneBeta2)/(4*beta2*one
243 (1-sqrtOneBeta2)/oneBeta2_to_3_
244
245 dsigma = ( firstTerm*(1-pi*fine_structure_co
246
247 return dsigma;
248 }
249
250 //
251
252 G4double G4RDPhotoElectricAngularGeneratorPola
253 {
254
255 //Double differential L1 shell cross-section
256
257 G4double beta2 = beta*beta;
258 G4double oneBeta2 = 1-beta2;
259 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
260 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1
261 G4double cosTheta = std::cos(theta);
262 G4double sinTheta2 =std::sin(theta)*std::sin
263 G4double cosPhi2 = std::cos(phi)*std::cos(ph
264 G4double oneBetaCosTheta = 1-beta*cosTheta;
265
266 G4double dsigma = 0;
267 G4double firstTerm = 0;
268 G4double secondTerm = 0;
269
270 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe
271 * (sinTheta2 * cosPhi2)/std::po
272 (1-sqrtOneBeta2)/(4*oneBeta2_to_
273
274 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st
275 (4*beta2/sqrtOneBeta2 * sinThet
276 - 4*(1-sqrtOneBeta2)/oneBeta2 *
277 + 4*beta2*(1-sqrtOneBeta2)/oneB
278 + (1-sqrtOneBeta2)/(4*beta2*one
279 (1-sqrtOneBeta2)/oneBeta2_to_3_
280
281 dsigma = ( firstTerm*(1-pi*fine_structure_co
282
283 return dsigma;
284 }
285
286 G4double G4RDPhotoElectricAngularGeneratorPola
287 {
288 if (arg1 > arg2)
289 return arg1;
290 else
291 return arg2;
292 }
293
294 //
295
296 G4RotationMatrix G4RDPhotoElectricAngularGener
297 const G4ThreeVector
298 {
299 G4double mK = direction.mag();
300 G4double mS = polarization.mag();
301 G4ThreeVector polarization2 = polarization;
302 const G4double kTolerance = 1e-6;
303
304 if(!(polarization.isOrthogonal(direction,kTo
305 G4ThreeVector d0 = direction.unit();
306 G4ThreeVector a1 = SetPerpendicularVector(
307 G4ThreeVector a0 = a1.unit();
308 G4double rand1 = G4UniformRand();
309 G4double angle = twopi*rand1;
310 G4ThreeVector b0 = d0.cross(a0);
311 G4ThreeVector c;
312 c.setX(std::cos(angle)*(a0.x())+std::sin(a
313 c.setY(std::cos(angle)*(a0.y())+std::sin(a
314 c.setZ(std::cos(angle)*(a0.z())+std::sin(a
315 polarization2 = c.unit();
316 mS = polarization2.mag();
317 }else
318 {
319 if ( polarization.howOrthogonal(directio
320 {
321 polarization2 = polarization - polarizatio
322 }
323 }
324
325 G4ThreeVector direction2 = direction/mK;
326 polarization2 = polarization2/mS;
327
328 G4ThreeVector y = direction2.cross(polarizat
329
330 G4RotationMatrix R(polarization2,y,direction
331 return R;
332 }
333
334 void G4RDPhotoElectricAngularGeneratorPolarize
335 {
336 // This member function finds for a given sh
337
338 G4double aBeta,cBeta;
339 G4double bMin,bStep;
340 G4int indexMax;
341 G4int level = shellLevel;
342 if(shellLevel > 1) level = 1; // protection
343
344 bMin = betaArray[0];
345 bStep = betaArray[1];
346 indexMax = (G4int)betaArray[2];
347 const G4double kBias = 1e-9;
348
349 G4int k = (G4int)((beta-bMin+kBias)/bStep);
350
351 if(k < 0)
352 k = 0;
353 if(k > indexMax)
354 k = indexMax;
355
356 if(k == 0)
357 aBeta = GetMax(aMajorantSurfaceParameterTa
358 else if(k==indexMax)
359 aBeta = GetMax(aMajorantSurfaceParameterTa
360 else{
361 aBeta = GetMax(aMajorantSurfaceParameterTa
362 aBeta = GetMax(aBeta,aMajorantSurfaceParam
363 }
364
365 if(k == 0)
366 cBeta = GetMax(cMajorantSurfaceParameterTa
367 else if(k == indexMax)
368 cBeta = GetMax(cMajorantSurfaceParameterTa
369 else{
370 cBeta = GetMax(cMajorantSurfaceParameterTa
371 cBeta = GetMax(cBeta,cMajorantSurfaceParam
372 }
373
374 *majorantSurfaceParameterA = aBeta;
375 *majorantSurfaceParameterC = cBeta;
376
377 }
378
379
380 //
381 G4ThreeVector G4RDPhotoElectricAngularGenerato
382 {
383
384 //computes the photoelectron momentum unitar
385 G4double px = std::cos(phi)*std::sin(theta);
386 G4double py = std::sin(phi)*std::sin(theta);
387 G4double pz = std::cos(theta);
388
389 G4ThreeVector samplingDirection(px,py,pz);
390
391 G4ThreeVector outgoingDirection = rotation*s
392 return outgoingDirection;
393 }
394
395 //
396
397 void G4RDPhotoElectricAngularGeneratorPolarize
398 {
399 G4cout << "\n" << G4endl;
400 G4cout << "Polarized Photoelectric Angular G
401 G4cout << "PhotoElectric Electron Angular Ge
402 G4cout << "Includes polarization effects for
403 G4cout << "For higher shells the L1 cross-se
404 G4cout << "(see Physics Reference Manual) \n
405 }
406
407 G4ThreeVector G4RDPhotoElectricAngularGenerato
408 {
409 G4double dx = a.x();
410 G4double dy = a.y();
411 G4double dz = a.z();
412 G4double x = dx < 0.0 ? -dx : dx;
413 G4double y = dy < 0.0 ? -dy : dy;
414 G4double z = dz < 0.0 ? -dz : dz;
415 if (x < y) {
416 return x < z ? G4ThreeVector(-dy,dx,0) : G
417 }else{
418 return y < z ? G4ThreeVector(dz,0,-dx) : G
419 }
420 }
421