G4PhotoElectricAngularGeneratorPolarized.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // -------------------------------------------------------------------
00028 //
00029 // GEANT4 Class file
00030 //
00031 //
00032 // File name:     G4PhotoElectricAngularGeneratorPolarized
00033 //
00034 // Author: A. C. Farinha, L. Peralta, P. Rodrigues and A. Trindade
00035 // 
00036 // Creation date:
00037 //
00038 // Modifications: 
00039 // 10 January 2006      
00040 // 06 May 2011, Replace FILE with std::ifstream
00041 //
00042 // Class Description: 
00043 //
00044 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation 
00045 //
00046 // Class Description: 
00047 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution.
00048 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961).
00049 // For higher shells the L1 cross-section is used. 
00050 //
00051 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using
00052 // the inverse-transform method (James 1980). Instead a more general approach based on the one already 
00053 // used to sample bremsstrahlung 2BN cross section (G4Generator2BN, Peralta, 2005) was used.
00054 //
00055 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526   (1959)
00056 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961)
00057 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980)
00058 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005)
00059 //
00060 //
00061 // -------------------------------------------------------------------
00062 //
00063 //    
00064 
00065 #include "G4PhotoElectricAngularGeneratorPolarized.hh"
00066 #include "G4PhysicalConstants.hh"
00067 #include "G4RotationMatrix.hh"
00068 #include "Randomize.hh"
00069 
00070 G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized()
00071   :G4VEmAngularDistribution("AngularGenSauterGavrilaPolarized")
00072 {
00073   const G4int arrayDim = 980;
00074 
00075   //minimum electron beta parameter allowed
00076   betaArray[0] = 0.02;
00077   //beta step
00078   betaArray[1] = 0.001;
00079   //maximum index array for a and c tables
00080   betaArray[2] = arrayDim - 1;
00081 
00082   // read Majorant Surface Parameters. This are required in order to generate Gavrila angular photoelectron distribution
00083   for(G4int level = 0; level < 2; level++){
00084 
00085     char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
00086     char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
00087 
00088     G4String filename;
00089     if(level == 0) filename = nameChar0;
00090     if(level == 1) filename = nameChar1;
00091 
00092     char* path = getenv("G4LEDATA");
00093     if (!path)
00094       {
00095         G4String excep = "G4EMDataSet - G4LEDATA environment variable not set";
00096         G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
00097                     "em0006",FatalException,"G4LEDATA environment variable not set");
00098         return;
00099       }
00100 
00101     G4String pathString(path);
00102     G4String dirFile = pathString + "/photoelectric_angular/" + filename;
00103     std::ifstream infile(dirFile);
00104     if (!infile.is_open())
00105       {
00106         G4String excep = "data file: " + dirFile + " not found";
00107         G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
00108                     "em0003",FatalException,excep);
00109         return;
00110       }
00111 
00112     // Read parameters into tables. The parameters are function of incident electron energy and shell level
00113     G4float aRead=0,cRead=0, beta=0;
00114     for(G4int i=0 ; i<arrayDim ;i++){
00115       //fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
00116       infile >> beta >> aRead >> cRead;    
00117       aMajorantSurfaceParameterTable[i][level] = aRead;    
00118       cMajorantSurfaceParameterTable[i][level] = cRead;
00119     }
00120     infile.close();
00121   }
00122 }
00123 
00124 G4PhotoElectricAngularGeneratorPolarized::~G4PhotoElectricAngularGeneratorPolarized() 
00125 {}
00126 
00127 G4ThreeVector& 
00128 G4PhotoElectricAngularGeneratorPolarized::SampleDirection(
00129                                  const G4DynamicParticle* dp,
00130                                  G4double eKinEnergy, 
00131                                  G4int shellId, 
00132                                  const G4Material*)
00133 {
00134   // (shellId == 0) - K-shell  - Polarized model for K-shell
00135   // (shellId > 0)  - L1-shell - Polarized model for L1 and higher shells
00136 
00137   // Calculate Lorentz term (gamma) and beta parameters
00138   G4double gamma = 1 + eKinEnergy/electron_mass_c2;
00139   G4double beta  = std::sqrt((gamma - 1)*(gamma + 1))/gamma;
00140 
00141   const G4ThreeVector& direction = dp->GetMomentumDirection();
00142   const G4ThreeVector& polarization = dp->GetPolarization();
00143 
00144   G4double theta, phi = 0;
00145   // Majorant surface parameters 
00146   // function of the outgoing electron kinetic energy
00147   G4double aBeta = 0; 
00148   G4double cBeta = 0; 
00149 
00150   // For the outgoing kinetic energy 
00151   // find the current majorant surface parameters
00152   PhotoElectronGetMajorantSurfaceAandCParameters(shellId, beta, &aBeta, &cBeta);
00153 
00154   // Generate pho and theta according to the shell level 
00155   // and beta parameter of the electron
00156   PhotoElectronGeneratePhiAndTheta(shellId, beta, aBeta, cBeta, &phi, &theta);
00157 
00158   // Determine the rotation matrix
00159   const G4RotationMatrix rotation = 
00160     PhotoElectronRotationMatrix(direction, polarization);
00161 
00162   // Compute final direction of the outgoing electron
00163   fLocalDirection = PhotoElectronComputeFinalDirection(rotation, theta, phi);
00164 
00165   return fLocalDirection;
00166 }
00167 
00168 void 
00169 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(
00170       G4int shellLevel, G4double beta, G4double aBeta, G4double cBeta, 
00171       G4double *pphi, G4double *ptheta) const
00172 {
00173   G4double rand1, rand2, rand3 = 0;
00174   G4double phi = 0;
00175   G4double theta = 0;
00176   G4double crossSectionValue = 0;
00177   G4double crossSectionMajorantFunctionValue = 0;
00178   G4double maxBeta = 0;
00179 
00180   //G4cout << "shell= " << shellLevel << " beta= " << beta
00181   //     << " aBeta= " << aBeta << " cBeta= " << cBeta << G4endl;
00182 
00183   do {
00184 
00185     rand1 = G4UniformRand();
00186     rand2 = G4UniformRand();
00187     rand3 = G4UniformRand();
00188         
00189     phi=2*pi*rand1;
00190 
00191     if(shellLevel == 0){
00192 
00193       // Polarized Gavrila Cross-Section for K-shell (1959)
00194       theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
00195       crossSectionMajorantFunctionValue = 
00196         CrossSectionMajorantFunction(theta, cBeta);
00197       crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
00198 
00199     } else {
00200 
00201       //  Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
00202       theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
00203       crossSectionMajorantFunctionValue = 
00204         CrossSectionMajorantFunction(theta, cBeta);
00205       crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
00206 
00207     }
00208 
00209     maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
00210     //G4cout << " crossSectionValue= " << crossSectionValue 
00211     //     << " max= " << maxBeta << G4endl;
00212     if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; }
00213 
00214   } while(maxBeta > crossSectionValue || theta > CLHEP::pi);
00215 
00216   *pphi = phi;
00217   *ptheta = theta;
00218 }
00219 
00220 G4double 
00221 G4PhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(
00222          G4double theta, G4double cBeta) const
00223 {
00224   // Compute Majorant Function
00225   G4double crossSectionMajorantFunctionValue = 0; 
00226   crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
00227   return crossSectionMajorantFunctionValue;
00228 }
00229 
00230 G4double 
00231 G4PhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(
00232          G4double beta, G4double theta, G4double phi) const
00233 {
00234   //Double differential K shell cross-section (Gavrila 1959)
00235 
00236   G4double beta2 = beta*beta;
00237   G4double oneBeta2 = 1 - beta2;
00238   G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
00239   G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
00240   G4double cosTheta = std::cos(theta);
00241   G4double sinTheta2 = std::sin(theta)*std::sin(theta);
00242   G4double cosPhi2 = std::cos(phi)*std::cos(phi);
00243   G4double oneBetaCosTheta = 1-beta*cosTheta;
00244   G4double dsigma = 0;
00245   G4double firstTerm = 0;
00246   G4double secondTerm = 0;
00247 
00248   firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) * 
00249               (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
00250               (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
00251 
00252   secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
00253                (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
00254                - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
00255                + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
00256                + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 
00257                (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
00258 
00259   dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
00260 
00261   return dsigma;
00262 }
00263 
00264 //
00265 
00266 G4double 
00267 G4PhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(
00268         G4double beta, G4double theta, G4double phi) const
00269 {
00270   //Double differential L1 shell cross-section (Gavrila 1961)
00271 
00272   G4double beta2 = beta*beta;
00273   G4double oneBeta2 = 1-beta2;
00274   G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
00275   G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
00276   G4double cosTheta = std::cos(theta);
00277   G4double sinTheta2 =std::sin(theta)*std::sin(theta);
00278   G4double cosPhi2 = std::cos(phi)*std::cos(phi);
00279   G4double oneBetaCosTheta = 1-beta*cosTheta;
00280         
00281   G4double dsigma = 0;
00282   G4double firstTerm = 0;
00283   G4double secondTerm = 0;
00284 
00285   firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
00286               *  (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
00287               (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
00288 
00289   secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
00290                (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
00291                - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
00292                + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
00293                + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 
00294                (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
00295 
00296   dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
00297 
00298   return dsigma;
00299 }
00300 
00301 G4RotationMatrix 
00302 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(
00303            const G4ThreeVector& direction,
00304            const G4ThreeVector& polarization)
00305 {
00306   G4double mK = direction.mag();
00307   G4double mS = polarization.mag();
00308   G4ThreeVector polarization2 = polarization;
00309   const G4double kTolerance = 1e-6;
00310 
00311   if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){
00312     G4ThreeVector d0 = direction.unit();
00313     G4ThreeVector a1 = PerpendicularVector(d0); 
00314     G4ThreeVector a0 = a1.unit(); 
00315     G4double rand1 = G4UniformRand();
00316     G4double angle = twopi*rand1; 
00317     G4ThreeVector b0 = d0.cross(a0); 
00318     G4ThreeVector c;
00319     c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
00320     c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
00321     c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
00322     polarization2 = c.unit();
00323     mS = polarization2.mag();
00324   }else
00325     {
00326       if ( polarization.howOrthogonal(direction) != 0)
00327         {
00328           polarization2 = polarization 
00329             - polarization.dot(direction)/direction.dot(direction) * direction;
00330         }
00331     }
00332 
00333   G4ThreeVector direction2 = direction/mK;
00334   polarization2 = polarization2/mS;
00335 
00336   G4ThreeVector y = direction2.cross(polarization2);
00337     
00338   G4RotationMatrix R(polarization2,y,direction2);
00339   return R;
00340 }
00341 
00342 void 
00343 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
00344 {
00345   // This member function finds for a given shell and beta value 
00346   // of the outgoing electron the correct Majorant Surface parameters
00347 
00348   G4double aBeta,cBeta;
00349   G4double bMin,bStep;
00350   G4int indexMax;
00351   G4int level = 0;    
00352   if(shellId > 0) { level = 1; }
00353 
00354   bMin = betaArray[0];
00355   bStep = betaArray[1];
00356   indexMax = (G4int)betaArray[2];
00357   const G4double kBias = 1e-9;
00358 
00359   G4int k = (G4int)((beta-bMin+kBias)/bStep);    
00360     
00361   if(k < 0)
00362     k = 0;
00363   if(k > indexMax)
00364     k = indexMax; 
00365     
00366   if(k == 0) 
00367     aBeta = std::max(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
00368   else if(k==indexMax)
00369     aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
00370   else{
00371     aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
00372     aBeta = std::max(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
00373   }   
00374     
00375   if(k == 0)
00376     cBeta = std::max(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
00377   else if(k == indexMax)
00378     cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
00379   else{
00380     cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
00381     cBeta = std::max(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
00382   }
00383 
00384   *majorantSurfaceParameterA = aBeta;
00385   *majorantSurfaceParameterC = cBeta;
00386 }
00387 
00388 G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, G4double theta, G4double phi) const
00389 {
00390   //computes the photoelectron momentum unitary vector 
00391   G4double sint = std::sin(theta);
00392   G4double px = std::cos(phi)*sint;
00393   G4double py = std::sin(phi)*sint;
00394   G4double pz = std::cos(theta);
00395 
00396   G4ThreeVector samplingDirection(px,py,pz);
00397 
00398   G4ThreeVector outgoingDirection = rotation*samplingDirection;
00399   return outgoingDirection;
00400 }
00401 
00402 void G4PhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const
00403 {
00404   G4cout << "\n" << G4endl;
00405   G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
00406   G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
00407   G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
00408   G4cout << "For higher shells the L1 cross-section is used." << G4endl;
00409   G4cout << "(see Physics Reference Manual) \n" << G4endl;
00410 } 
00411 
00412 G4ThreeVector 
00413 G4PhotoElectricAngularGeneratorPolarized::PerpendicularVector(
00414          const G4ThreeVector& a) const
00415 {
00416   G4double dx = a.x();
00417   G4double dy = a.y();
00418   G4double dz = a.z();
00419   G4double x = dx < 0.0 ? -dx : dx;
00420   G4double y = dy < 0.0 ? -dy : dy;
00421   G4double z = dz < 0.0 ? -dz : dz;
00422   if (x < y) {
00423     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
00424   }else{
00425     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
00426   }
00427 }

Generated on Mon May 27 17:49:18 2013 for Geant4 by  doxygen 1.4.7