00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
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
00076 betaArray[0] = 0.02;
00077
00078 betaArray[1] = 0.001;
00079
00080 betaArray[2] = arrayDim - 1;
00081
00082
00083 for(G4int level = 0; level < 2; level++){
00084
00085 char nameChar0[100] = "ftab0.dat";
00086 char nameChar1[100] = "ftab1.dat";
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
00113 G4float aRead=0,cRead=0, beta=0;
00114 for(G4int i=0 ; i<arrayDim ;i++){
00115
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
00135
00136
00137
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
00146
00147 G4double aBeta = 0;
00148 G4double cBeta = 0;
00149
00150
00151
00152 PhotoElectronGetMajorantSurfaceAandCParameters(shellId, beta, &aBeta, &cBeta);
00153
00154
00155
00156 PhotoElectronGeneratePhiAndTheta(shellId, beta, aBeta, cBeta, &phi, &theta);
00157
00158
00159 const G4RotationMatrix rotation =
00160 PhotoElectronRotationMatrix(direction, polarization);
00161
00162
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
00181
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
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
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
00211
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
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
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
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
00346
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
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 }