Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhotoElectricAngularGeneratorPolarized.cc
Go to the documentation of this file.
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 
89  G4String filename;
90  if(level == 0) filename = nameChar0;
91  if(level == 1) filename = nameChar1;
92 
93  char* path = 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
170 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(
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
222 G4PhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(
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
232 G4PhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(
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) );
261 
262  return dsigma;
263 }
264 
265 //
266 
267 G4double
268 G4PhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(
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) );
298 
299  return dsigma;
300 }
301 
303 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(
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();
316  G4double rand1 = G4UniformRand();
317  G4double angle = twopi*rand1;
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)
368  aBeta = std::max(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
369  else if(k==indexMax)
370  aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
371  else{
372  aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
373  aBeta = std::max(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
374  }
375 
376  if(k == 0)
377  cBeta = std::max(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
378  else if(k == indexMax)
379  cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
380  else{
381  cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
382  cBeta = std::max(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
383  }
384 
385  *majorantSurfaceParameterA = aBeta;
386  *majorantSurfaceParameterC = cBeta;
387 }
388 
389 G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, G4double theta, G4double phi) const
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 }
const G4double a0
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
float G4float
Definition: G4Types.hh:77
static G4double angle[DIM]
tuple x
Definition: test.py:50
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double eKinEnergy, G4int shellId, const G4Material *mat=0)
int G4int
Definition: G4Types.hh:78
void setY(double)
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
double z() const
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
float electron_mass_c2
Definition: hepunit.py:274
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ThreeVector PerpendicularVector(const G4ThreeVector &a) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Hep3Vector unit() const
double y() const
const G4ThreeVector & GetPolarization() const
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
double mag() const
static constexpr double pi
Definition: SystemOfUnits.h:54