Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PhotoElectricAngularGeneratorPolarized.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 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4PhotoElectricAngularGeneratorPolarized.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PhotoElectricAngularGeneratorPolarized.cc (Version 10.0.p1)


  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