Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 
65 #include "G4RotationMatrix.hh"
66 #include "Randomize.hh"
67 #include "G4PhysicalConstants.hh"
68 
69 //
70 
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 = 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 
126 {;}
127 
128 //
129 
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 
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 
208 {
209  // Compute Majorant Function
210  G4double crossSectionMajorantFunctionValue = 0;
211  crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
212  return crossSectionMajorantFunctionValue;
213 }
214 
215 //
216 
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 
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 
287 {
288  if (arg1 > arg2)
289  return arg1;
290  else
291  return arg2;
292 }
293 
294 //
295 
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();
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)
358  else if(k==indexMax)
360  else{
362  aBeta = GetMax(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
363  }
364 
365  if(k == 0)
367  else if(k == indexMax)
369  else{
371  cBeta = GetMax(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
372  }
373 
374  *majorantSurfaceParameterA = aBeta;
375  *majorantSurfaceParameterC = cBeta;
376 
377 }
378 
379 
380 //
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 
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 
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 }
Float_t x
Definition: compare.C:6
const XML_Char * name
Definition: expat.h:151
void PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta, const G4double aBeta, const G4double cBeta, G4double *pphi, G4double *ptheta) const
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
float G4float
Definition: G4Types.hh:77
Double_t z
double z() const
double dot(const Hep3Vector &) const
Double_t beta
G4double DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const
G4double CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
void setX(double)
fclose(fg1)
double G4double
Definition: G4Types.hh:76
void setZ(double)
static constexpr double electron_mass_c2
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
Double_t R
G4ThreeVector GetPhotoElectronDirection(const G4ThreeVector &direction, const G4double kineticEnergy, const G4ThreeVector &polarization, const G4int shellId) const
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
G4double GetMax(const G4double arg1, const G4double arg2) const
G4double DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector unit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
const G4double a0
double mag() const
void PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
int G4int
Definition: G4Types.hh:78
G4RotationMatrix PhotoElectronRotationMatrix(const G4ThreeVector &direction, const G4ThreeVector &polarization) const
G4ThreeVector SetPerpendicularVector(const G4ThreeVector &a) const
G4GLOB_DLL std::ostream G4cout
double x() const
static constexpr double pi
Definition: G4SIunits.hh:75
double y() const
static constexpr double fine_structure_const
G4ThreeVector PhotoElectronComputeFinalDirection(const G4RotationMatrix &rotation, const G4double theta, const G4double phi) const
void setY(double)