source: trunk/source/processes/electromagnetic/lowenergy/src/G4PhotoElectricAngularGeneratorPolarized.cc@ 1058

Last change on this file since 1058 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 16.0 KB
Line 
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//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4PhotoElectricAngularGeneratorPolarized
33//
34// Author: A. C. Farinha, L. Peralta, P. Rodrigues and A. Trindade
35//
36// Creation date:
37//
38// Modifications:
39// 10 January 2006
40//
41// Class Description:
42//
43// Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation
44//
45// Class Description:
46// PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution.
47// Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961).
48// For higher shells the L1 cross-section is used.
49//
50// The Gavrila photoelectron angular distribution is a complex function which can not be sampled using
51// the inverse-transform method (James 1980). Instead a more general approach based on the one already
52// used to sample bremsstrahlung 2BN cross section (G4Generator2BN, Peralta, 2005) was used.
53//
54// M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526 (1959)
55// M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961)
56// F. James, Rept. on Prog. in Phys. 43, 1145 (1980)
57// L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005)
58//
59//
60// -------------------------------------------------------------------
61//
62//
63
64#include "G4PhotoElectricAngularGeneratorPolarized.hh"
65#include "G4RotationMatrix.hh"
66#include "Randomize.hh"
67
68//
69
70G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized(const G4String& name):G4VPhotoElectricAngularDistribution(name)
71{
72 const G4int arrayDim = 980;
73
74 //minimum electron beta parameter allowed
75 betaArray[0] = 0.02;
76 //beta step
77 betaArray[1] = 0.001;
78 //maximum index array for a and c tables
79 betaArray[2] = arrayDim - 1;
80
81 // read Majorant Surface Parameters. This are required in order to generate Gavrila angular photoelectron distribution
82 for(G4int level = 0; level < 2; level++){
83
84 char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
85 char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
86
87 G4String filename;
88 if(level == 0) filename = nameChar0;
89 if(level == 1) filename = nameChar1;
90
91 char* path = getenv("G4LEDATA");
92 if (!path)
93 {
94 G4String excep = "G4EMDataSet - G4LEDATA environment variable not set";
95 G4Exception(excep);
96 }
97
98 G4String pathString(path);
99 G4String dirFile = pathString + "/photoelectric_angular/" + filename;
100 FILE *infile;
101 infile = fopen(dirFile,"r");
102 if (infile == 0)
103 {
104 G4String excep = "G4PhotoElectricAngularGeneratorPolarized - data file: " + dirFile + " not found";
105 G4Exception(excep);
106 }
107
108 // Read parameters into tables. The parameters are function of incident electron energy and shell level
109 G4float aRead,cRead, beta;
110 for(G4int i=0 ; i<arrayDim ;i++){
111 fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
112 aMajorantSurfaceParameterTable[i][level] = aRead;
113 cMajorantSurfaceParameterTable[i][level] = cRead;
114 }
115 fclose(infile);
116
117 }
118}
119
120//
121
122G4PhotoElectricAngularGeneratorPolarized::~G4PhotoElectricAngularGeneratorPolarized()
123{;}
124
125//
126
127G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::GetPhotoElectronDirection(const G4ThreeVector& direction, const G4double eKineticEnergy,
128 const G4ThreeVector& polarization, const G4int shellId) const
129{
130 // Calculate Lorentz term (gamma) and beta parameters
131 G4double gamma = 1. + eKineticEnergy/electron_mass_c2;
132 G4double beta = std::sqrt(gamma*gamma-1.)/gamma;
133
134 G4double theta, phi = 0;
135 G4double aBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
136 G4double cBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
137
138 G4int shellLevel = 0;
139 if(shellId < 2) shellLevel = 0; // K-shell // Polarized model for K-shell
140 if(shellId >= 2) shellLevel = 1; // L1-shell // Polarized model for L1 and higher shells
141
142 // For the outgoing kinetic energy find the current majorant surface parameters
143 PhotoElectronGetMajorantSurfaceAandCParameters( shellLevel, beta, &aBeta, &cBeta);
144
145 // Generate pho and theta according to the shell level and beta parameter of the electron
146 PhotoElectronGeneratePhiAndTheta(shellLevel, beta, aBeta, cBeta, &phi, &theta);
147
148 // Determine the rotation matrix
149 G4RotationMatrix rotation = PhotoElectronRotationMatrix(direction, polarization);
150
151 // Compute final direction of the outgoing electron
152 G4ThreeVector final_direction = PhotoElectronComputeFinalDirection(rotation, theta, phi);
153
154 return final_direction;
155}
156
157//
158
159void G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta,
160 const G4double aBeta, const G4double cBeta,
161 G4double *pphi, G4double *ptheta) const
162{
163 G4double rand1, rand2, rand3 = 0;
164 G4double phi = 0;
165 G4double theta = 0;
166 G4double crossSectionValue = 0;
167 G4double crossSectionMajorantFunctionValue = 0;
168 G4double maxBeta = 0;
169
170 do {
171
172 rand1 = G4UniformRand();
173 rand2 = G4UniformRand();
174 rand3 = G4UniformRand();
175
176 phi=2*pi*rand1;
177
178 if(shellLevel == 0){
179
180 // Polarized Gavrila Cross-Section for K-shell (1959)
181 theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
182 crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
183 crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
184
185 } else {
186
187 // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
188 theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
189 crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
190 crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
191
192 }
193
194 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
195
196 }while(maxBeta > crossSectionValue);
197
198 *pphi = phi;
199 *ptheta = theta;
200}
201
202//
203
204G4double G4PhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const
205{
206 // Compute Majorant Function
207 G4double crossSectionMajorantFunctionValue = 0;
208 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
209 return crossSectionMajorantFunctionValue;
210}
211
212//
213
214G4double G4PhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const
215{
216
217 //Double differential K shell cross-section (Gavrila 1959)
218
219 G4double beta2 = beta*beta;
220 G4double oneBeta2 = 1 - beta2;
221 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
222 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
223 G4double cosTheta = std::cos(theta);
224 G4double sinTheta2 = std::sin(theta)*std::sin(theta);
225 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
226 G4double oneBetaCosTheta = 1-beta*cosTheta;
227 G4double dsigma = 0;
228 G4double firstTerm = 0;
229 G4double secondTerm = 0;
230
231 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
232 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
233 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
234
235 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
236 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
237 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
238 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
239 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
240 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
241
242 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
243
244 return dsigma;
245}
246
247//
248
249G4double G4PhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const
250{
251
252 //Double differential L1 shell cross-section (Gavrila 1961)
253
254 G4double beta2 = beta*beta;
255 G4double oneBeta2 = 1-beta2;
256 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
257 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
258 G4double cosTheta = std::cos(theta);
259 G4double sinTheta2 =std::sin(theta)*std::sin(theta);
260 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
261 G4double oneBetaCosTheta = 1-beta*cosTheta;
262
263 G4double dsigma = 0;
264 G4double firstTerm = 0;
265 G4double secondTerm = 0;
266
267 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
268 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
269 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
270
271 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
272 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
273 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
274 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
275 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
276 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
277
278 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
279
280 return dsigma;
281}
282
283G4double G4PhotoElectricAngularGeneratorPolarized::GetMax(const G4double arg1, const G4double arg2) const
284{
285 if (arg1 > arg2)
286 return arg1;
287 else
288 return arg2;
289}
290
291//
292
293G4RotationMatrix G4PhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(const G4ThreeVector& direction,
294 const G4ThreeVector& polarization) const
295{
296 G4double mK = direction.mag();
297 G4double mS = polarization.mag();
298 G4ThreeVector polarization2 = polarization;
299 const G4double kTolerance = 1e-6;
300
301 if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){
302 G4ThreeVector d0 = direction.unit();
303 G4ThreeVector a1 = SetPerpendicularVector(d0);
304 G4ThreeVector a0 = a1.unit();
305 G4double rand1 = G4UniformRand();
306 G4double angle = twopi*rand1;
307 G4ThreeVector b0 = d0.cross(a0);
308 G4ThreeVector c;
309 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
310 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
311 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
312 polarization2 = c.unit();
313 mS = polarization2.mag();
314 }else
315 {
316 if ( polarization.howOrthogonal(direction) != 0)
317 {
318 polarization2 = polarization - polarization.dot(direction)/direction.dot(direction) * direction;
319 }
320 }
321
322 G4ThreeVector direction2 = direction/mK;
323 polarization2 = polarization2/mS;
324
325 G4ThreeVector y = direction2.cross(polarization2);
326
327 G4RotationMatrix R(polarization2,y,direction2);
328 return R;
329}
330
331void G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta,G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
332{
333 // This member function finds for a given shell and beta value of the outgoing electron the correct Majorant Surface parameters
334
335 G4double aBeta,cBeta;
336 G4double bMin,bStep;
337 G4int indexMax;
338 G4int level = shellLevel;
339 if(shellLevel > 1) level = 1; // protection since only K and L1 polarized double differential cross-sections were implemented
340
341 bMin = betaArray[0];
342 bStep = betaArray[1];
343 indexMax = (G4int)betaArray[2];
344 const G4double kBias = 1e-9;
345
346 G4int k = (G4int)((beta-bMin+kBias)/bStep);
347
348 if(k < 0)
349 k = 0;
350 if(k > indexMax)
351 k = indexMax;
352
353 if(k == 0)
354 aBeta = GetMax(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
355 else if(k==indexMax)
356 aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
357 else{
358 aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
359 aBeta = GetMax(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
360 }
361
362 if(k == 0)
363 cBeta = GetMax(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
364 else if(k == indexMax)
365 cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
366 else{
367 cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
368 cBeta = GetMax(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
369 }
370
371 *majorantSurfaceParameterA = aBeta;
372 *majorantSurfaceParameterC = cBeta;
373
374}
375
376
377//
378G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, const G4double theta, const G4double phi) const
379{
380
381 //computes the photoelectron momentum unitary vector
382 G4double px = std::cos(phi)*std::sin(theta);
383 G4double py = std::sin(phi)*std::sin(theta);
384 G4double pz = std::cos(theta);
385
386 G4ThreeVector samplingDirection(px,py,pz);
387
388 G4ThreeVector outgoingDirection = rotation*samplingDirection;
389 return outgoingDirection;
390}
391
392//
393
394void G4PhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const
395{
396 G4cout << "\n" << G4endl;
397 G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
398 G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
399 G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
400 G4cout << "For higher shells the L1 cross-section is used." << G4endl;
401 G4cout << "(see Physics Reference Manual) \n" << G4endl;
402}
403
404G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::SetPerpendicularVector(const G4ThreeVector& a) const
405{
406 G4double dx = a.x();
407 G4double dy = a.y();
408 G4double dz = a.z();
409 G4double x = dx < 0.0 ? -dx : dx;
410 G4double y = dy < 0.0 ? -dy : dy;
411 G4double z = dz < 0.0 ? -dz : dz;
412 if (x < y) {
413 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
414 }else{
415 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
416 }
417}
Note: See TracBrowser for help on using the repository browser.