| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4PhotoElectricAngularGenera 32 // File name: G4PhotoElectricAngularGeneratorPolarized
33 // 33 //
34 // Author: A. C. Farinha, L. Peralta, P. Rodri 34 // Author: A. C. Farinha, L. Peralta, P. Rodrigues and A. Trindade
35 // << 35 //
36 // Creation date: 36 // Creation date:
37 // 37 //
38 // Modifications: << 38 // Modifications:
39 // 10 January 2006 << 39 // 10 January 2006
40 // 06 May 2011, Replace FILE with std::ifstrea 40 // 06 May 2011, Replace FILE with std::ifstream
41 // 41 //
42 // Class Description: << 42 // Class Description:
43 // 43 //
44 // Concrete class for PhotoElectric Electron A << 44 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation
45 // 45 //
46 // Class Description: << 46 // Class Description:
47 // PhotoElectric Electron Angular Generator ba 47 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution.
48 // Includes polarization effects for K and L1 48 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961).
49 // For higher shells the L1 cross-section is u << 49 // For higher shells the L1 cross-section is used.
50 // 50 //
51 // The Gavrila photoelectron angular distribut 51 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using
52 // the inverse-transform method (James 1980). << 52 // the inverse-transform method (James 1980). Instead a more general approach based on the one already
53 // used to sample bremsstrahlung 2BN cross sec 53 // used to sample bremsstrahlung 2BN cross section (G4Generator2BN, Peralta, 2005) was used.
54 // 54 //
55 // M. Gavrila, "Relativistic K-Shell Photoeffe 55 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526 (1959)
56 // M. Gavrila, "Relativistic L-Shell Photoeffe 56 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961)
57 // F. James, Rept. on Prog. in Phys. 43, 1145 57 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980)
58 // L. Peralta et al., "A new low-energy bremss 58 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005)
59 // 59 //
60 // 60 //
61 // ------------------------------------------- 61 // -------------------------------------------------------------------
62 // 62 //
63 // << 63 //
64 64
65 #include "G4PhotoElectricAngularGeneratorPolar 65 #include "G4PhotoElectricAngularGeneratorPolarized.hh"
66 #include "G4PhysicalConstants.hh" 66 #include "G4PhysicalConstants.hh"
67 #include "G4RotationMatrix.hh" 67 #include "G4RotationMatrix.hh"
68 #include "Randomize.hh" 68 #include "Randomize.hh"
69 #include "G4Exp.hh" <<
70 69
71 G4PhotoElectricAngularGeneratorPolarized::G4Ph 70 G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized()
72 :G4VEmAngularDistribution("AngularGenSauterG 71 :G4VEmAngularDistribution("AngularGenSauterGavrilaPolarized")
73 { 72 {
74 const G4int arrayDim = 980; 73 const G4int arrayDim = 980;
75 74
76 //minimum electron beta parameter allowed 75 //minimum electron beta parameter allowed
77 betaArray[0] = 0.02; 76 betaArray[0] = 0.02;
78 //beta step 77 //beta step
79 betaArray[1] = 0.001; 78 betaArray[1] = 0.001;
80 //maximum index array for a and c tables 79 //maximum index array for a and c tables
81 betaArray[2] = arrayDim - 1; 80 betaArray[2] = arrayDim - 1;
82 81
83 // read Majorant Surface Parameters. This ar << 82 // read Majorant Surface Parameters. This are required in order to generate Gavrila angular photoelectron distribution
84 //Gavrila angular photoelectron distribution <<
85 for(G4int level = 0; level < 2; level++){ 83 for(G4int level = 0; level < 2; level++){
>> 84
86 char nameChar0[100] = "ftab0.dat"; // K-sh 85 char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
87 char nameChar1[100] = "ftab1.dat"; // L-sh 86 char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
88 87
89 G4String filename; 88 G4String filename;
90 if(level == 0) filename = nameChar0; 89 if(level == 0) filename = nameChar0;
91 if(level == 1) filename = nameChar1; 90 if(level == 1) filename = nameChar1;
92 91
93 const char* path = G4FindDataDir("G4LEDATA << 92 char* path = getenv("G4LEDATA");
94 if (!path) 93 if (!path)
95 { 94 {
96 G4String excep = "G4EMDataSet - G4LEDA 95 G4String excep = "G4EMDataSet - G4LEDATA environment variable not set";
97 G4Exception("G4PhotoElectricAngularGen 96 G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
98 "em0006",FatalException,"G4LEDATA envi 97 "em0006",FatalException,"G4LEDATA environment variable not set");
99 return; 98 return;
100 } 99 }
101 100
102 G4String pathString(path); 101 G4String pathString(path);
103 G4String dirFile = pathString + "/photoele 102 G4String dirFile = pathString + "/photoelectric_angular/" + filename;
104 std::ifstream infile(dirFile); 103 std::ifstream infile(dirFile);
105 if (!infile.is_open()) 104 if (!infile.is_open())
106 { 105 {
107 G4String excep = "data file: " + dirFile + " 106 G4String excep = "data file: " + dirFile + " not found";
108 G4Exception("G4PhotoElectricAngularGen 107 G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
109 "em0003",FatalException,excep); 108 "em0003",FatalException,excep);
110 return; 109 return;
111 } 110 }
112 111
113 // Read parameters into tables. The parame << 112 // Read parameters into tables. The parameters are function of incident electron energy and shell level
114 // energy and shell level <<
115 G4float aRead=0,cRead=0, beta=0; 113 G4float aRead=0,cRead=0, beta=0;
116 for(G4int i=0 ; i<arrayDim ;++i){ << 114 for(G4int i=0 ; i<arrayDim ;i++){
117 infile >> beta >> aRead >> cRead; << 115 //fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
118 aMajorantSurfaceParameterTable[i][level] << 116 infile >> beta >> aRead >> cRead;
>> 117 aMajorantSurfaceParameterTable[i][level] = aRead;
119 cMajorantSurfaceParameterTable[i][level] 118 cMajorantSurfaceParameterTable[i][level] = cRead;
120 } 119 }
121 infile.close(); 120 infile.close();
122 } 121 }
123 } 122 }
124 123
125 //....oooOO0OOooo........oooOO0OOooo........oo << 124 G4PhotoElectricAngularGeneratorPolarized::~G4PhotoElectricAngularGeneratorPolarized()
126 <<
127 G4PhotoElectricAngularGeneratorPolarized::~G4P <<
128 {} 125 {}
129 126
130 //....oooOO0OOooo........oooOO0OOooo........oo << 127 G4ThreeVector&
131 <<
132 G4ThreeVector& <<
133 G4PhotoElectricAngularGeneratorPolarized::Samp 128 G4PhotoElectricAngularGeneratorPolarized::SampleDirection(
134 const G4Dynam 129 const G4DynamicParticle* dp,
135 G4double eKinEnergy, << 130 G4double eKinEnergy,
136 G4int shellId, << 131 G4int shellId,
137 const G4Material*) 132 const G4Material*)
138 { 133 {
139 // (shellId == 0) - K-shell - Polarized mod 134 // (shellId == 0) - K-shell - Polarized model for K-shell
140 // (shellId > 0) - L1-shell - Polarized mod 135 // (shellId > 0) - L1-shell - Polarized model for L1 and higher shells
141 136
142 // Calculate Lorentz term (gamma) and beta p 137 // Calculate Lorentz term (gamma) and beta parameters
143 G4double gamma = 1 + eKinEnergy/electron_mas 138 G4double gamma = 1 + eKinEnergy/electron_mass_c2;
144 G4double beta = std::sqrt((gamma - 1)*(gamm 139 G4double beta = std::sqrt((gamma - 1)*(gamma + 1))/gamma;
145 140
146 const G4ThreeVector& direction = dp->GetMome 141 const G4ThreeVector& direction = dp->GetMomentumDirection();
147 const G4ThreeVector& polarization = dp->GetP 142 const G4ThreeVector& polarization = dp->GetPolarization();
148 143
149 G4double theta, phi = 0; 144 G4double theta, phi = 0;
150 // Majorant surface parameters << 145 // Majorant surface parameters
151 // function of the outgoing electron kinetic 146 // function of the outgoing electron kinetic energy
152 G4double aBeta = 0; << 147 G4double aBeta = 0;
153 G4double cBeta = 0; << 148 G4double cBeta = 0;
154 149
155 // For the outgoing kinetic energy << 150 // For the outgoing kinetic energy
156 // find the current majorant surface paramet 151 // find the current majorant surface parameters
157 PhotoElectronGetMajorantSurfaceAandCParamete 152 PhotoElectronGetMajorantSurfaceAandCParameters(shellId, beta, &aBeta, &cBeta);
158 153
159 // Generate pho and theta according to the s << 154 // Generate pho and theta according to the shell level
160 // and beta parameter of the electron 155 // and beta parameter of the electron
161 PhotoElectronGeneratePhiAndTheta(shellId, be 156 PhotoElectronGeneratePhiAndTheta(shellId, beta, aBeta, cBeta, &phi, &theta);
162 157
163 // Determine the rotation matrix 158 // Determine the rotation matrix
164 const G4RotationMatrix rotation = << 159 const G4RotationMatrix rotation =
165 PhotoElectronRotationMatrix(direction, pol 160 PhotoElectronRotationMatrix(direction, polarization);
166 161
167 // Compute final direction of the outgoing e 162 // Compute final direction of the outgoing electron
168 fLocalDirection = PhotoElectronComputeFinalD 163 fLocalDirection = PhotoElectronComputeFinalDirection(rotation, theta, phi);
169 164
170 return fLocalDirection; 165 return fLocalDirection;
171 } 166 }
172 167
173 //....oooOO0OOooo........oooOO0OOooo........oo << 168 void
174 <<
175 void <<
176 G4PhotoElectricAngularGeneratorPolarized::Phot 169 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(
177 G4int shellLevel, G4double beta, G4doubl << 170 G4int shellLevel, G4double beta, G4double aBeta, G4double cBeta,
178 G4double *pphi, G4double *ptheta) const 171 G4double *pphi, G4double *ptheta) const
179 { 172 {
180 G4double rand1, rand2, rand3 = 0; 173 G4double rand1, rand2, rand3 = 0;
181 G4double phi = 0; 174 G4double phi = 0;
182 G4double theta = 0; 175 G4double theta = 0;
183 G4double crossSectionValue = 0; 176 G4double crossSectionValue = 0;
184 G4double crossSectionMajorantFunctionValue = 177 G4double crossSectionMajorantFunctionValue = 0;
185 G4double maxBeta = 0; 178 G4double maxBeta = 0;
186 179
>> 180 //G4cout << "shell= " << shellLevel << " beta= " << beta
>> 181 // << " aBeta= " << aBeta << " cBeta= " << cBeta << G4endl;
>> 182
187 do { 183 do {
188 184
189 rand1 = G4UniformRand(); 185 rand1 = G4UniformRand();
190 rand2 = G4UniformRand(); 186 rand2 = G4UniformRand();
191 rand3 = G4UniformRand(); 187 rand3 = G4UniformRand();
192 << 188
193 phi=2*pi*rand1; 189 phi=2*pi*rand1;
194 190
195 if(shellLevel == 0){ 191 if(shellLevel == 0){
>> 192
196 // Polarized Gavrila Cross-Section for K 193 // Polarized Gavrila Cross-Section for K-shell (1959)
197 theta=std::sqrt(((G4Exp(rand2*std::log(1 << 194 theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
198 crossSectionMajorantFunctionValue = << 195 crossSectionMajorantFunctionValue =
199 CrossSectionMajorantFunction(theta, cBeta); 196 CrossSectionMajorantFunction(theta, cBeta);
200 crossSectionValue = DSigmaKshellGavrila1 197 crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
>> 198
201 } else { 199 } else {
>> 200
202 // Polarized Gavrila Cross-Section for 201 // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
203 theta = std::sqrt(((G4Exp(rand2*std::log << 202 theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
204 crossSectionMajorantFunctionValue = << 203 crossSectionMajorantFunctionValue =
205 CrossSectionMajorantFunction(theta, cBeta); 204 CrossSectionMajorantFunction(theta, cBeta);
206 crossSectionValue = DSigmaL1shellGavrila 205 crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
>> 206
207 } 207 }
208 208
209 maxBeta=rand3*aBeta*crossSectionMajorantFu 209 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
>> 210 //G4cout << " crossSectionValue= " << crossSectionValue
>> 211 // << " max= " << maxBeta << G4endl;
210 if(crossSectionValue < 0.0) { crossSection 212 if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; }
211 213
212 } while(maxBeta > crossSectionValue || theta 214 } while(maxBeta > crossSectionValue || theta > CLHEP::pi);
213 215
214 *pphi = phi; 216 *pphi = phi;
215 *ptheta = theta; 217 *ptheta = theta;
216 } 218 }
217 219
218 //....oooOO0OOooo........oooOO0OOooo........oo << 220 G4double
219 <<
220 G4double <<
221 G4PhotoElectricAngularGeneratorPolarized::Cros 221 G4PhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(
222 G4double theta, G4double cBeta) const 222 G4double theta, G4double cBeta) const
223 { 223 {
224 // Compute Majorant Function 224 // Compute Majorant Function
225 G4double crossSectionMajorantFunctionValue = << 225 G4double crossSectionMajorantFunctionValue = 0;
226 crossSectionMajorantFunctionValue = theta/(1 226 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
227 return crossSectionMajorantFunctionValue; 227 return crossSectionMajorantFunctionValue;
228 } 228 }
229 229
230 //....oooOO0OOooo........oooOO0OOooo........oo << 230 G4double
231 <<
232 G4double <<
233 G4PhotoElectricAngularGeneratorPolarized::DSig 231 G4PhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(
234 G4double beta, G4double theta, G4doub 232 G4double beta, G4double theta, G4double phi) const
235 { 233 {
236 //Double differential K shell cross-section 234 //Double differential K shell cross-section (Gavrila 1959)
237 235
238 G4double beta2 = beta*beta; 236 G4double beta2 = beta*beta;
239 G4double oneBeta2 = 1 - beta2; 237 G4double oneBeta2 = 1 - beta2;
240 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 238 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
241 G4double oneBeta2_to_3_2 = std::pow(oneBeta2 239 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
242 G4double cosTheta = std::cos(theta); 240 G4double cosTheta = std::cos(theta);
243 G4double sinTheta2 = std::sin(theta)*std::si 241 G4double sinTheta2 = std::sin(theta)*std::sin(theta);
244 G4double cosPhi2 = std::cos(phi)*std::cos(ph 242 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
245 G4double oneBetaCosTheta = 1-beta*cosTheta; 243 G4double oneBetaCosTheta = 1-beta*cosTheta;
246 G4double dsigma = 0; 244 G4double dsigma = 0;
247 G4double firstTerm = 0; 245 G4double firstTerm = 0;
248 G4double secondTerm = 0; 246 G4double secondTerm = 0;
249 247
250 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe << 248 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
251 (sinTheta2 * cosPhi2)/std::pow(o 249 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
252 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 250 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
253 251
254 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 252 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
255 (4*beta2/sqrtOneBeta2 * sinThet 253 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
256 - 4*(1-sqrtOneBeta2)/oneBeta2 * 254 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
257 + 4*beta2*(1-sqrtOneBeta2)/oneB 255 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
258 + (1-sqrtOneBeta2)/(4*beta2*one << 256 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
259 (1-sqrtOneBeta2)/oneBeta2_to_3_ 257 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
260 258
261 dsigma = ( firstTerm*(1-pi*fine_structure_co << 259 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
262 260
263 return dsigma; 261 return dsigma;
264 } 262 }
265 263
266 //....oooOO0OOooo........oooOO0OOooo........oo << 264 //
267 265
268 G4double << 266 G4double
269 G4PhotoElectricAngularGeneratorPolarized::DSig 267 G4PhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(
270 G4double beta, G4double theta, G4doubl 268 G4double beta, G4double theta, G4double phi) const
271 { 269 {
272 //Double differential L1 shell cross-section 270 //Double differential L1 shell cross-section (Gavrila 1961)
273 // 26Oct2022: included factor (1/8) in dsigm << 271
274 G4double beta2 = beta*beta; 272 G4double beta2 = beta*beta;
275 G4double oneBeta2 = 1-beta2; 273 G4double oneBeta2 = 1-beta2;
276 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 274 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
277 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1 275 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
278 G4double cosTheta = std::cos(theta); 276 G4double cosTheta = std::cos(theta);
279 G4double sinTheta2 =std::sin(theta)*std::sin 277 G4double sinTheta2 =std::sin(theta)*std::sin(theta);
280 G4double cosPhi2 = std::cos(phi)*std::cos(ph 278 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
281 G4double oneBetaCosTheta = 1-beta*cosTheta; 279 G4double oneBetaCosTheta = 1-beta*cosTheta;
282 << 280
283 G4double dsigma = 0; 281 G4double dsigma = 0;
284 G4double firstTerm = 0; 282 G4double firstTerm = 0;
285 G4double secondTerm = 0; 283 G4double secondTerm = 0;
286 284
287 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe 285 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
288 * (sinTheta2 * cosPhi2)/std::po 286 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
289 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 287 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
290 288
291 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 289 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
292 (4*beta2/sqrtOneBeta2 * sinThet 290 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
293 - 4*(1-sqrtOneBeta2)/oneBeta2 * 291 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
294 + 4*beta2*(1-sqrtOneBeta2)/oneB 292 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
295 + (1-sqrtOneBeta2)/(4*beta2*one << 293 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
296 (1-sqrtOneBeta2)/oneBeta2_to_3_ 294 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
297 295
298 dsigma = ( firstTerm*(1-pi*fine_structure_co << 296 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
299 297
300 return dsigma; 298 return dsigma;
301 } 299 }
302 300
303 //....oooOO0OOooo........oooOO0OOooo........oo << 301 G4RotationMatrix
304 <<
305 G4RotationMatrix <<
306 G4PhotoElectricAngularGeneratorPolarized::Phot 302 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(
307 const G4ThreeVector& direction, 303 const G4ThreeVector& direction,
308 const G4ThreeVector& polarization) 304 const G4ThreeVector& polarization)
309 { 305 {
310 G4double mK = direction.mag(); 306 G4double mK = direction.mag();
311 G4double mS = polarization.mag(); 307 G4double mS = polarization.mag();
312 G4ThreeVector polarization2 = polarization; 308 G4ThreeVector polarization2 = polarization;
313 const G4double kTolerance = 1e-6; 309 const G4double kTolerance = 1e-6;
314 310
315 if(!(polarization.isOrthogonal(direction,kTo 311 if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){
316 G4ThreeVector d0 = direction.unit(); 312 G4ThreeVector d0 = direction.unit();
317 G4ThreeVector a1 = PerpendicularVector(d0) << 313 G4ThreeVector a1 = PerpendicularVector(d0);
318 G4ThreeVector a0 = a1.unit(); << 314 G4ThreeVector a0 = a1.unit();
319 G4double rand1 = G4UniformRand(); 315 G4double rand1 = G4UniformRand();
320 G4double angle = twopi*rand1; << 316 G4double angle = twopi*rand1;
321 G4ThreeVector b0 = d0.cross(a0); << 317 G4ThreeVector b0 = d0.cross(a0);
322 G4ThreeVector c; 318 G4ThreeVector c;
323 c.setX(std::cos(angle)*(a0.x())+std::sin(a 319 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
324 c.setY(std::cos(angle)*(a0.y())+std::sin(a 320 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
325 c.setZ(std::cos(angle)*(a0.z())+std::sin(a 321 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
326 polarization2 = c.unit(); 322 polarization2 = c.unit();
327 mS = polarization2.mag(); 323 mS = polarization2.mag();
328 }else 324 }else
329 { 325 {
330 if ( polarization.howOrthogonal(directio 326 if ( polarization.howOrthogonal(direction) != 0)
331 { 327 {
332 polarization2 = polarization << 328 polarization2 = polarization
333 - polarization.dot(direction)/direction. 329 - polarization.dot(direction)/direction.dot(direction) * direction;
334 } 330 }
335 } 331 }
336 332
337 G4ThreeVector direction2 = direction/mK; 333 G4ThreeVector direction2 = direction/mK;
338 polarization2 = polarization2/mS; 334 polarization2 = polarization2/mS;
339 335
340 G4ThreeVector y = direction2.cross(polarizat 336 G4ThreeVector y = direction2.cross(polarization2);
341 << 337
342 G4RotationMatrix R(polarization2,y,direction 338 G4RotationMatrix R(polarization2,y,direction2);
343 return R; 339 return R;
344 } 340 }
345 341
346 //....oooOO0OOooo........oooOO0OOooo........oo << 342 void
347 <<
348 void <<
349 G4PhotoElectricAngularGeneratorPolarized::Phot 343 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
350 { 344 {
351 // This member function finds for a given sh << 345 // This member function finds for a given shell and beta value
352 // of the outgoing electron the correct Majo 346 // of the outgoing electron the correct Majorant Surface parameters
>> 347
353 G4double aBeta,cBeta; 348 G4double aBeta,cBeta;
354 G4double bMin,bStep; 349 G4double bMin,bStep;
355 G4int indexMax; 350 G4int indexMax;
356 G4int level = 0; << 351 G4int level = 0;
357 if(shellId > 0) { level = 1; } 352 if(shellId > 0) { level = 1; }
358 353
359 bMin = betaArray[0]; 354 bMin = betaArray[0];
360 bStep = betaArray[1]; 355 bStep = betaArray[1];
361 indexMax = (G4int)betaArray[2]; 356 indexMax = (G4int)betaArray[2];
362 const G4double kBias = 1e-9; 357 const G4double kBias = 1e-9;
363 358
364 G4int k = (G4int)((beta-bMin+kBias)/bStep); << 359 G4int k = (G4int)((beta-bMin+kBias)/bStep);
365 << 360
366 if(k < 0) 361 if(k < 0)
367 k = 0; 362 k = 0;
368 if(k > indexMax) 363 if(k > indexMax)
369 k = indexMax; << 364 k = indexMax;
370 << 365
371 if(k == 0) << 366 if(k == 0)
372 aBeta = std::max(aMajorantSurfaceParameter 367 aBeta = std::max(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
373 else if(k==indexMax) 368 else if(k==indexMax)
374 aBeta = std::max(aMajorantSurfaceParameter 369 aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
375 else{ 370 else{
376 aBeta = std::max(aMajorantSurfaceParameter 371 aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
377 aBeta = std::max(aBeta,aMajorantSurfacePar 372 aBeta = std::max(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
378 } << 373 }
379 << 374
380 if(k == 0) 375 if(k == 0)
381 cBeta = std::max(cMajorantSurfaceParameter 376 cBeta = std::max(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
382 else if(k == indexMax) 377 else if(k == indexMax)
383 cBeta = std::max(cMajorantSurfaceParameter 378 cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
384 else{ 379 else{
385 cBeta = std::max(cMajorantSurfaceParameter 380 cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
386 cBeta = std::max(cBeta,cMajorantSurfacePar 381 cBeta = std::max(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
387 } 382 }
388 383
389 *majorantSurfaceParameterA = aBeta; 384 *majorantSurfaceParameterA = aBeta;
390 *majorantSurfaceParameterC = cBeta; 385 *majorantSurfaceParameterC = cBeta;
391 } 386 }
392 387
393 //....oooOO0OOooo........oooOO0OOooo........oo <<
394 <<
395 G4ThreeVector G4PhotoElectricAngularGeneratorP 388 G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, G4double theta, G4double phi) const
396 { 389 {
397 //computes the photoelectron momentum unitar << 390 //computes the photoelectron momentum unitary vector
398 G4double sint = std::sin(theta); 391 G4double sint = std::sin(theta);
399 G4double px = std::cos(phi)*sint; 392 G4double px = std::cos(phi)*sint;
400 G4double py = std::sin(phi)*sint; 393 G4double py = std::sin(phi)*sint;
401 G4double pz = std::cos(theta); 394 G4double pz = std::cos(theta);
402 395
403 G4ThreeVector samplingDirection(px,py,pz); 396 G4ThreeVector samplingDirection(px,py,pz);
404 397
405 G4ThreeVector outgoingDirection = rotation*s 398 G4ThreeVector outgoingDirection = rotation*samplingDirection;
406 return outgoingDirection; 399 return outgoingDirection;
407 } 400 }
408 401
409 //....oooOO0OOooo........oooOO0OOooo........oo <<
410 <<
411 void G4PhotoElectricAngularGeneratorPolarized: 402 void G4PhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const
412 { 403 {
413 G4cout << "\n" << G4endl; 404 G4cout << "\n" << G4endl;
414 G4cout << "Polarized Photoelectric Angular G 405 G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
415 G4cout << "PhotoElectric Electron Angular Ge 406 G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
416 G4cout << "Includes polarization effects for 407 G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
417 G4cout << "For higher shells the L1 cross-se 408 G4cout << "For higher shells the L1 cross-section is used." << G4endl;
418 G4cout << "(see Physics Reference Manual) \n 409 G4cout << "(see Physics Reference Manual) \n" << G4endl;
419 } << 410 }
420 <<
421 //....oooOO0OOooo........oooOO0OOooo........oo <<
422 411
423 G4ThreeVector << 412 G4ThreeVector
424 G4PhotoElectricAngularGeneratorPolarized::Perp 413 G4PhotoElectricAngularGeneratorPolarized::PerpendicularVector(
425 const G4ThreeVector& a) const 414 const G4ThreeVector& a) const
426 { 415 {
427 G4double dx = a.x(); 416 G4double dx = a.x();
428 G4double dy = a.y(); 417 G4double dy = a.y();
429 G4double dz = a.z(); 418 G4double dz = a.z();
430 G4double x = dx < 0.0 ? -dx : dx; 419 G4double x = dx < 0.0 ? -dx : dx;
431 G4double y = dy < 0.0 ? -dy : dy; 420 G4double y = dy < 0.0 ? -dy : dy;
432 G4double z = dz < 0.0 ? -dz : dz; 421 G4double z = dz < 0.0 ? -dz : dz;
433 if (x < y) { 422 if (x < y) {
434 return x < z ? G4ThreeVector(-dy,dx,0) : G 423 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
435 }else{ 424 }else{
436 return y < z ? G4ThreeVector(dz,0,-dx) : G 425 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
437 } 426 }
438 } 427 }
439 428