ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDPhotoElectricAngularGeneratorPolarized.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RDPhotoElectricAngularGeneratorPolarized.cc
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: G4RDPhotoElectricAngularGeneratorPolarized
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 (G4RDGenerator2BN, 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 "G4RDPhotoElectricAngularGeneratorPolarized.hh"
65 #include "G4RotationMatrix.hh"
66 #include "Randomize.hh"
67 #include "G4PhysicalConstants.hh"
68 
69 //
70 
71 G4RDPhotoElectricAngularGeneratorPolarized::G4RDPhotoElectricAngularGeneratorPolarized(const G4String& name):G4RDVPhotoElectricAngularDistribution(name)
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 are required in order to generate Gavrila angular photoelectron distribution
83  for(G4int level = 0; level < 2; level++){
84 
85  char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
86  char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
87 
88  G4String filename;
89  if(level == 0) filename = nameChar0;
90  if(level == 1) filename = nameChar1;
91 
92  char* path = std::getenv("G4LEDATA");
93  if (!path)
94  {
95  G4String excep = "G4LEDATA environment variable not set!";
96  G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()",
97  "InvalidSetup", FatalException, excep);
98  }
99 
100  G4String pathString(path);
101  G4String dirFile = pathString + "/photoelectric_angular/" + filename;
102  FILE *infile;
103  infile = fopen(dirFile,"r");
104  if (infile == 0)
105  {
106  G4String excep = "Data file: " + dirFile + " not found";
107  G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()",
108  "DataNotFound", FatalException, excep);
109  }
110 
111  // Read parameters into tables. The parameters are function of incident electron energy and shell level
112  G4float aRead,cRead, beta;
113  for(G4int i=0 ; i<arrayDim ;i++){
114  fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
115  aMajorantSurfaceParameterTable[i][level] = aRead;
116  cMajorantSurfaceParameterTable[i][level] = cRead;
117  }
118  fclose(infile);
119 
120  }
121 }
122 
123 //
124 
125 G4RDPhotoElectricAngularGeneratorPolarized::~G4RDPhotoElectricAngularGeneratorPolarized()
126 {;}
127 
128 //
129 
130 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::GetPhotoElectronDirection(const G4ThreeVector& direction, const G4double eKineticEnergy,
131  const G4ThreeVector& polarization, const G4int shellId) const
132 {
133  // Calculate Lorentz term (gamma) and beta parameters
134  G4double gamma = 1. + eKineticEnergy/electron_mass_c2;
135  G4double beta = std::sqrt(gamma*gamma-1.)/gamma;
136 
137  G4double theta, phi = 0;
138  G4double aBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
139  G4double cBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
140 
141  G4int shellLevel = 0;
142  if(shellId < 2) shellLevel = 0; // K-shell // Polarized model for K-shell
143  if(shellId >= 2) shellLevel = 1; // L1-shell // Polarized model for L1 and higher shells
144 
145  // For the outgoing kinetic energy find the current majorant surface parameters
146  PhotoElectronGetMajorantSurfaceAandCParameters( shellLevel, beta, &aBeta, &cBeta);
147 
148  // Generate pho and theta according to the shell level and beta parameter of the electron
149  PhotoElectronGeneratePhiAndTheta(shellLevel, beta, aBeta, cBeta, &phi, &theta);
150 
151  // Determine the rotation matrix
152  G4RotationMatrix rotation = PhotoElectronRotationMatrix(direction, polarization);
153 
154  // Compute final direction of the outgoing electron
155  G4ThreeVector final_direction = PhotoElectronComputeFinalDirection(rotation, theta, phi);
156 
157  return final_direction;
158 }
159 
160 //
161 
162 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta,
163  const G4double aBeta, const G4double cBeta,
164  G4double *pphi, G4double *ptheta) const
165 {
166  G4double rand1, rand2, rand3 = 0;
167  G4double phi = 0;
168  G4double theta = 0;
169  G4double crossSectionValue = 0;
170  G4double crossSectionMajorantFunctionValue = 0;
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-shell (1959)
184  theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
185  crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
186  crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
187 
188  } else {
189 
190  // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
191  theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
192  crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
193  crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
194 
195  }
196 
197  maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
198 
199  }while(maxBeta > crossSectionValue);
200 
201  *pphi = phi;
202  *ptheta = theta;
203 }
204 
205 //
206 
207 G4double G4RDPhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const
208 {
209  // Compute Majorant Function
210  G4double crossSectionMajorantFunctionValue = 0;
211  crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
212  return crossSectionMajorantFunctionValue;
213 }
214 
215 //
216 
217 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const
218 {
219 
220  //Double differential K shell cross-section (Gavrila 1959)
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,1.5);
226  G4double cosTheta = std::cos(theta);
227  G4double sinTheta2 = std::sin(theta)*std::sin(theta);
228  G4double cosPhi2 = std::cos(phi)*std::cos(phi);
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(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
235  (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
236  (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
237 
238  secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
239  (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
240  - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
241  + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
242  + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
243  (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
244 
245  dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
246 
247  return dsigma;
248 }
249 
250 //
251 
252 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const
253 {
254 
255  //Double differential L1 shell cross-section (Gavrila 1961)
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.5);
261  G4double cosTheta = std::cos(theta);
262  G4double sinTheta2 =std::sin(theta)*std::sin(theta);
263  G4double cosPhi2 = std::cos(phi)*std::cos(phi);
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(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
271  * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
272  (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
273 
274  secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
275  (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
276  - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
277  + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
278  + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
279  (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
280 
281  dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
282 
283  return dsigma;
284 }
285 
286 G4double G4RDPhotoElectricAngularGeneratorPolarized::GetMax(const G4double arg1, const G4double arg2) const
287 {
288  if (arg1 > arg2)
289  return arg1;
290  else
291  return arg2;
292 }
293 
294 //
295 
296 G4RotationMatrix G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(const G4ThreeVector& direction,
297  const G4ThreeVector& polarization) const
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,kTolerance)) || mS == 0){
305  G4ThreeVector d0 = direction.unit();
306  G4ThreeVector a1 = SetPerpendicularVector(d0);
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(angle)*b0.x());
313  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
314  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
315  polarization2 = c.unit();
316  mS = polarization2.mag();
317  }else
318  {
319  if ( polarization.howOrthogonal(direction) != 0)
320  {
321  polarization2 = polarization - polarization.dot(direction)/direction.dot(direction) * direction;
322  }
323  }
324 
325  G4ThreeVector direction2 = direction/mK;
326  polarization2 = polarization2/mS;
327 
328  G4ThreeVector y = direction2.cross(polarization2);
329 
330  G4RotationMatrix R(polarization2,y,direction2);
331  return R;
332 }
333 
334 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta,G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
335 {
336  // This member function finds for a given shell and beta value of the outgoing electron the correct Majorant Surface parameters
337 
338  G4double aBeta,cBeta;
339  G4double bMin,bStep;
340  G4int indexMax;
341  G4int level = shellLevel;
342  if(shellLevel > 1) level = 1; // protection since only K and L1 polarized double differential cross-sections were implemented
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(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
358  else if(k==indexMax)
359  aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
360  else{
361  aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
362  aBeta = GetMax(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
363  }
364 
365  if(k == 0)
366  cBeta = GetMax(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
367  else if(k == indexMax)
368  cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
369  else{
370  cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
371  cBeta = GetMax(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
372  }
373 
374  *majorantSurfaceParameterA = aBeta;
375  *majorantSurfaceParameterC = cBeta;
376 
377 }
378 
379 
380 //
381 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, const G4double theta, const G4double phi) const
382 {
383 
384  //computes the photoelectron momentum unitary vector
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*samplingDirection;
392  return outgoingDirection;
393 }
394 
395 //
396 
397 void G4RDPhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const
398 {
399  G4cout << "\n" << G4endl;
400  G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
401  G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
402  G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
403  G4cout << "For higher shells the L1 cross-section is used." << G4endl;
404  G4cout << "(see Physics Reference Manual) \n" << G4endl;
405 }
406 
407 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::SetPerpendicularVector(const G4ThreeVector& a) const
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) : G4ThreeVector(0,-dz,dy);
417  }else{
418  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
419  }
420 }