ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhotoElectricAngularGeneratorPolarized.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PhotoElectricAngularGeneratorPolarized.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: 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 // 06 May 2011, Replace FILE with std::ifstream
41 //
42 // Class Description:
43 //
44 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation
45 //
46 // Class Description:
47 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution.
48 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961).
49 // For higher shells the L1 cross-section is used.
50 //
51 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using
52 // the inverse-transform method (James 1980). Instead a more general approach based on the one already
53 // used to sample bremsstrahlung 2BN cross section (G4Generator2BN, Peralta, 2005) was used.
54 //
55 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526 (1959)
56 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961)
57 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980)
58 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005)
59 //
60 //
61 // -------------------------------------------------------------------
62 //
63 //
64 
66 #include "G4PhysicalConstants.hh"
67 #include "G4RotationMatrix.hh"
68 #include "Randomize.hh"
69 #include "G4Exp.hh"
70 
72  :G4VEmAngularDistribution("AngularGenSauterGavrilaPolarized")
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 are required in order to generate Gavrila angular photoelectron distribution
84  for(G4int level = 0; level < 2; level++){
85 
86  char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
87  char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
88 
90  if(level == 0) filename = nameChar0;
91  if(level == 1) filename = nameChar1;
92 
93  char* path = std::getenv("G4LEDATA");
94  if (!path)
95  {
96  G4String excep = "G4EMDataSet - G4LEDATA environment variable not set";
97  G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
98  "em0006",FatalException,"G4LEDATA environment variable not set");
99  return;
100  }
101 
102  G4String pathString(path);
103  G4String dirFile = pathString + "/photoelectric_angular/" + filename;
104  std::ifstream infile(dirFile);
105  if (!infile.is_open())
106  {
107  G4String excep = "data file: " + dirFile + " not found";
108  G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
109  "em0003",FatalException,excep);
110  return;
111  }
112 
113  // Read parameters into tables. The parameters are function of incident electron energy and shell level
114  G4float aRead=0,cRead=0, beta=0;
115  for(G4int i=0 ; i<arrayDim ;i++){
116  //fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
117  infile >> beta >> aRead >> cRead;
118  aMajorantSurfaceParameterTable[i][level] = aRead;
119  cMajorantSurfaceParameterTable[i][level] = cRead;
120  }
121  infile.close();
122  }
123 }
124 
126 {}
127 
130  const G4DynamicParticle* dp,
131  G4double eKinEnergy,
132  G4int shellId,
133  const G4Material*)
134 {
135  // (shellId == 0) - K-shell - Polarized model for K-shell
136  // (shellId > 0) - L1-shell - Polarized model for L1 and higher shells
137 
138  // Calculate Lorentz term (gamma) and beta parameters
139  G4double gamma = 1 + eKinEnergy/electron_mass_c2;
140  G4double beta = std::sqrt((gamma - 1)*(gamma + 1))/gamma;
141 
142  const G4ThreeVector& direction = dp->GetMomentumDirection();
143  const G4ThreeVector& polarization = dp->GetPolarization();
144 
145  G4double theta, phi = 0;
146  // Majorant surface parameters
147  // function of the outgoing electron kinetic energy
148  G4double aBeta = 0;
149  G4double cBeta = 0;
150 
151  // For the outgoing kinetic energy
152  // find the current majorant surface parameters
153  PhotoElectronGetMajorantSurfaceAandCParameters(shellId, beta, &aBeta, &cBeta);
154 
155  // Generate pho and theta according to the shell level
156  // and beta parameter of the electron
157  PhotoElectronGeneratePhiAndTheta(shellId, beta, aBeta, cBeta, &phi, &theta);
158 
159  // Determine the rotation matrix
160  const G4RotationMatrix rotation =
161  PhotoElectronRotationMatrix(direction, polarization);
162 
163  // Compute final direction of the outgoing electron
164  fLocalDirection = PhotoElectronComputeFinalDirection(rotation, theta, phi);
165 
166  return fLocalDirection;
167 }
168 
169 void
171  G4int shellLevel, G4double beta, G4double aBeta, G4double cBeta,
172  G4double *pphi, G4double *ptheta) const
173 {
174  G4double rand1, rand2, rand3 = 0;
175  G4double phi = 0;
176  G4double theta = 0;
177  G4double crossSectionValue = 0;
178  G4double crossSectionMajorantFunctionValue = 0;
179  G4double maxBeta = 0;
180 
181  //G4cout << "shell= " << shellLevel << " beta= " << beta
182  // << " aBeta= " << aBeta << " cBeta= " << cBeta << G4endl;
183 
184  do {
185 
186  rand1 = G4UniformRand();
187  rand2 = G4UniformRand();
188  rand3 = G4UniformRand();
189 
190  phi=2*pi*rand1;
191 
192  if(shellLevel == 0){
193 
194  // Polarized Gavrila Cross-Section for K-shell (1959)
195  theta=std::sqrt(((G4Exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
196  crossSectionMajorantFunctionValue =
197  CrossSectionMajorantFunction(theta, cBeta);
198  crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
199 
200  } else {
201 
202  // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
203  theta = std::sqrt(((G4Exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
204  crossSectionMajorantFunctionValue =
205  CrossSectionMajorantFunction(theta, cBeta);
206  crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
207 
208  }
209 
210  maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
211  //G4cout << " crossSectionValue= " << crossSectionValue
212  // << " max= " << maxBeta << G4endl;
213  if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; }
214 
215  } while(maxBeta > crossSectionValue || theta > CLHEP::pi);
216 
217  *pphi = phi;
218  *ptheta = theta;
219 }
220 
221 G4double
223  G4double theta, G4double cBeta) const
224 {
225  // Compute Majorant Function
226  G4double crossSectionMajorantFunctionValue = 0;
227  crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
228  return crossSectionMajorantFunctionValue;
229 }
230 
231 G4double
233  G4double beta, G4double theta, G4double phi) const
234 {
235  //Double differential K shell cross-section (Gavrila 1959)
236 
237  G4double beta2 = beta*beta;
238  G4double oneBeta2 = 1 - beta2;
239  G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
240  G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
241  G4double cosTheta = std::cos(theta);
242  G4double sinTheta2 = std::sin(theta)*std::sin(theta);
243  G4double cosPhi2 = std::cos(phi)*std::cos(phi);
244  G4double oneBetaCosTheta = 1-beta*cosTheta;
245  G4double dsigma = 0;
246  G4double firstTerm = 0;
247  G4double secondTerm = 0;
248 
249  firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
250  (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
251  (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
252 
253  secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
254  (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
255  - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
256  + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
257  + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
258  (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
259 
260  dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) )*std::sin(theta);
261 
262  return dsigma;
263 }
264 
265 //
266 
267 G4double
269  G4double beta, G4double theta, G4double phi) const
270 {
271  //Double differential L1 shell cross-section (Gavrila 1961)
272 
273  G4double beta2 = beta*beta;
274  G4double oneBeta2 = 1-beta2;
275  G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
276  G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
277  G4double cosTheta = std::cos(theta);
278  G4double sinTheta2 =std::sin(theta)*std::sin(theta);
279  G4double cosPhi2 = std::cos(phi)*std::cos(phi);
280  G4double oneBetaCosTheta = 1-beta*cosTheta;
281 
282  G4double dsigma = 0;
283  G4double firstTerm = 0;
284  G4double secondTerm = 0;
285 
286  firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
287  * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
288  (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
289 
290  secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
291  (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
292  - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
293  + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
294  + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
295  (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
296 
297  dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) )*std::sin(theta);
298 
299  return dsigma;
300 }
301 
304  const G4ThreeVector& direction,
305  const G4ThreeVector& polarization)
306 {
307  G4double mK = direction.mag();
308  G4double mS = polarization.mag();
309  G4ThreeVector polarization2 = polarization;
310  const G4double kTolerance = 1e-6;
311 
312  if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){
313  G4ThreeVector d0 = direction.unit();
315  G4ThreeVector a0 = a1.unit();
318  G4ThreeVector b0 = d0.cross(a0);
320  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
321  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
322  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
323  polarization2 = c.unit();
324  mS = polarization2.mag();
325  }else
326  {
327  if ( polarization.howOrthogonal(direction) != 0)
328  {
329  polarization2 = polarization
330  - polarization.dot(direction)/direction.dot(direction) * direction;
331  }
332  }
333 
334  G4ThreeVector direction2 = direction/mK;
335  polarization2 = polarization2/mS;
336 
337  G4ThreeVector y = direction2.cross(polarization2);
338 
339  G4RotationMatrix R(polarization2,y,direction2);
340  return R;
341 }
342 
343 void
344 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
345 {
346  // This member function finds for a given shell and beta value
347  // of the outgoing electron the correct Majorant Surface parameters
348 
349  G4double aBeta,cBeta;
350  G4double bMin,bStep;
351  G4int indexMax;
352  G4int level = 0;
353  if(shellId > 0) { level = 1; }
354 
355  bMin = betaArray[0];
356  bStep = betaArray[1];
357  indexMax = (G4int)betaArray[2];
358  const G4double kBias = 1e-9;
359 
360  G4int k = (G4int)((beta-bMin+kBias)/bStep);
361 
362  if(k < 0)
363  k = 0;
364  if(k > indexMax)
365  k = indexMax;
366 
367  if(k == 0)
369  else if(k==indexMax)
371  else{
373  aBeta = std::max(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
374  }
375 
376  if(k == 0)
378  else if(k == indexMax)
380  else{
382  cBeta = std::max(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
383  }
384 
385  *majorantSurfaceParameterA = aBeta;
386  *majorantSurfaceParameterC = cBeta;
387 }
388 
390 {
391  //computes the photoelectron momentum unitary vector
392  G4double sint = std::sin(theta);
393  G4double px = std::cos(phi)*sint;
394  G4double py = std::sin(phi)*sint;
395  G4double pz = std::cos(theta);
396 
397  G4ThreeVector samplingDirection(px,py,pz);
398 
399  G4ThreeVector outgoingDirection = rotation*samplingDirection;
400  return outgoingDirection;
401 }
402 
404 {
405  G4cout << "\n" << G4endl;
406  G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
407  G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
408  G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
409  G4cout << "For higher shells the L1 cross-section is used." << G4endl;
410  G4cout << "(see Physics Reference Manual) \n" << G4endl;
411 }
412 
415  const G4ThreeVector& a) const
416 {
417  G4double dx = a.x();
418  G4double dy = a.y();
419  G4double dz = a.z();
420  G4double x = dx < 0.0 ? -dx : dx;
421  G4double y = dy < 0.0 ? -dy : dy;
422  G4double z = dz < 0.0 ? -dz : dz;
423  if (x < y) {
424  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
425  }else{
426  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
427  }
428 }