Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDPhotoElectricAngularGeneratorPolarized.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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     const char* path = G4FindDataDir("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 }
421