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: G4BetheHeitler5DModel.cc 32 // File name: G4BetheHeitler5DModel.cc 33 // 33 // 34 // Authors: 34 // Authors: 35 // Igor Semeniouk and Denis Bernard, 35 // Igor Semeniouk and Denis Bernard, 36 // LLR, Ecole polytechnique & CNRS/IN2P3, 9112 36 // LLR, Ecole polytechnique & CNRS/IN2P3, 91128 Palaiseau, France 37 // 37 // 38 // Acknowledgement of the support of the Frenc 38 // Acknowledgement of the support of the French National Research Agency 39 // (ANR-13-BS05-0002). 39 // (ANR-13-BS05-0002). 40 // 40 // 41 // Reference: Nucl. Instrum. Meth. A 899 (2018 41 // Reference: Nucl. Instrum. Meth. A 899 (2018) 85 (arXiv:1802.08253 [hep-ph]) 42 // Nucl. Instrum. Meth., A 936 (201 << 43 // 42 // 44 // Class Description: 43 // Class Description: 45 // 44 // 46 // Generates the conversion of a high-energy p << 45 // Generates the conversion of a high-energy photon to an e+e- pair, either in the field of an 47 // atomic electron (triplet) or nucleus (nucle 46 // atomic electron (triplet) or nucleus (nuclear). 48 // Samples the five-dimensional (5D) different 47 // Samples the five-dimensional (5D) differential cross-section analytical expression: 49 // . Non polarized conversion: 48 // . Non polarized conversion: 50 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lon 49 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lond. Ser. A 146 (1934) 83. 51 // . Polarized conversion: 50 // . Polarized conversion: 52 // T. H. Berlin and L. Madansky, Phys. Rev. 51 // T. H. Berlin and L. Madansky, Phys. Rev. 78 (1950) 623, 53 // M. M. May, Phys. Rev. 84 (1951) 265, 52 // M. M. May, Phys. Rev. 84 (1951) 265, 54 // J. M. Jauch and F. Rohrlich, The theory o 53 // J. M. Jauch and F. Rohrlich, The theory of photons and electrons, 1976. 55 // 54 // 56 // All the above expressions are named "Bethe- 55 // All the above expressions are named "Bethe-Heitler" here. 57 // 56 // 58 // Bethe & Heitler, put in Feynman diagram par 57 // Bethe & Heitler, put in Feynman diagram parlance, compute only the two dominant diagrams of 59 // the first order Born development, which is 58 // the first order Born development, which is an excellent approximation for nuclear conversion 60 // and for high-energy triplet conversion. 59 // and for high-energy triplet conversion. 61 // 60 // 62 // Only the linear polarisation of the incomin 61 // Only the linear polarisation of the incoming photon takes part in these expressions. 63 // The circular polarisation of the incoming p << 62 // The circular polarisation of the incoming photon does not (take part) and no polarisation 64 // is transfered to the final leptons. 63 // is transfered to the final leptons. 65 // 64 // 66 // In case conversion takes place in the field << 65 // In case conversion takes place in the field of an isolated nucleus or electron, the bare 67 // Bethe-Heitler expression is used. 66 // Bethe-Heitler expression is used. 68 // 67 // 69 // In case the nucleus or the electron are par << 68 // In case the nucleus or the electron are part of an atom, the screening of the target field 70 // by the other electrons of the atom is descr 69 // by the other electrons of the atom is described by a simple form factor, function of q2: 71 // . nuclear: N.F. Mott, H.S.W. Massey, The Th 70 // . nuclear: N.F. Mott, H.S.W. Massey, The Theory of Atomic Collisions, 1934. 72 // . triplet: J.A. Wheeler and W.E. Lamb, Phys 71 // . triplet: J.A. Wheeler and W.E. Lamb, Phys. Rev. 55 (1939) 858. 73 // 72 // 74 // The nuclear form factor that affects the pr 73 // The nuclear form factor that affects the probability of very large-q2 events, is not considered. 75 // 74 // 76 // In principle the code is valid from thresho << 75 // In principle the code is valid from threshold, that is from 2 * m_e c^2 for nuclear and from 77 // 4 * m_e c^2 for triplet, up to infinity, wh << 76 // 4 * m_e c^2 for triplet, up to infinity, while in pratice the divergence of the differential 78 // cross section at small q2 and, at high-ener << 77 // cross section at small q2 and, at high-energy, at small polar angle, make it break down at 79 // some point that depends on machine precisio 78 // some point that depends on machine precision. 80 // 79 // 81 // Very-high-energy (above a few tens of TeV) 80 // Very-high-energy (above a few tens of TeV) LPM suppression effects in the normalized differential 82 // cross-section are not considered. 81 // cross-section are not considered. 83 // 82 // 84 // The 5D differential cross section is sample << 83 // The 5D differential cross section is sampled without any high-energy nor small 85 // angle approximation(s). 84 // angle approximation(s). 86 // The generation is strictly energy-momentum << 85 // The generation is strictly energy-momentum conserving when all particles in the final state 87 // are taken into account, that is, including 86 // are taken into account, that is, including the recoiling target. 88 // (In contrast with the BH expressions taken 87 // (In contrast with the BH expressions taken at face values, for which the electron energy is 89 // taken to be EMinus = GammaEnergy - EPlus) 88 // taken to be EMinus = GammaEnergy - EPlus) 90 // 89 // 91 // Tests include the examination of 1D distrib 90 // Tests include the examination of 1D distributions: see TestEm15 92 // 91 // 93 // Total cross sections are not computed (we i 92 // Total cross sections are not computed (we inherit from other classes). 94 // We just convert a photon on a target when a 93 // We just convert a photon on a target when asked to do so. 95 // 94 // 96 // Pure nuclear, pure triplet and 1/Z triplet/ << 95 // Pure nuclear, pure triplet and 1/Z triplet/nuclear mixture can be generated. 97 // 96 // 98 // ------------------------------------------- 97 // ------------------------------------------------------------------- 99 98 100 #include "G4BetheHeitler5DModel.hh" 99 #include "G4BetheHeitler5DModel.hh" 101 #include "G4EmParameters.hh" 100 #include "G4EmParameters.hh" 102 101 103 #include "G4PhysicalConstants.hh" 102 #include "G4PhysicalConstants.hh" 104 #include "G4SystemOfUnits.hh" 103 #include "G4SystemOfUnits.hh" 105 #include "G4Electron.hh" 104 #include "G4Electron.hh" 106 #include "G4Positron.hh" 105 #include "G4Positron.hh" 107 #include "G4Gamma.hh" 106 #include "G4Gamma.hh" 108 #include "G4IonTable.hh" 107 #include "G4IonTable.hh" 109 #include "G4NucleiProperties.hh" 108 #include "G4NucleiProperties.hh" 110 109 111 #include "Randomize.hh" 110 #include "Randomize.hh" 112 #include "G4ParticleChangeForGamma.hh" 111 #include "G4ParticleChangeForGamma.hh" 113 #include "G4Pow.hh" 112 #include "G4Pow.hh" 114 #include "G4Log.hh" 113 #include "G4Log.hh" 115 #include "G4Exp.hh" 114 #include "G4Exp.hh" 116 115 117 #include "G4LorentzVector.hh" 116 #include "G4LorentzVector.hh" 118 #include "G4ThreeVector.hh" 117 #include "G4ThreeVector.hh" 119 #include "G4RotationMatrix.hh" << 120 << 121 #include <cassert> << 122 << 123 const G4int kEPair = 0; << 124 const G4int kMuPair = 1; << 125 << 126 118 127 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 128 120 129 G4BetheHeitler5DModel::G4BetheHeitler5DModel(c 121 G4BetheHeitler5DModel::G4BetheHeitler5DModel(const G4ParticleDefinition* pd, 130 c 122 const G4String& nam) 131 : G4PairProductionRelModel(pd, nam), << 123 : G4BetheHeitlerModel(pd, nam), fVerbose(1), fConversionType(0), iraw(false) 132 fLepton1(G4Electron::Definition()),fLepton << 133 fTheMuPlus(nullptr),fTheMuMinus(nullptr), << 134 fVerbose(1), << 135 fConversionType(0), << 136 fConvMode(kEPair), << 137 iraw(false) << 138 { 124 { >> 125 SetLowEnergyLimit(2*CLHEP::electron_mass_c2); 139 theIonTable = G4IonTable::GetIonTable(); 126 theIonTable = G4IonTable::GetIonTable(); 140 //Q: Do we need this on Model << 127 // Verbosity levels: ( Can redefine as needed, but some consideration ) 141 SetLowEnergyLimit(2*fTheElectron->GetPDGMass << 128 // 0 = nothing >> 129 // > 2 print results >> 130 // > 3 print rejection warning from transformation (fix bug from gammaray .. ) >> 131 // > 4 print photon direction & polarisation 142 } 132 } 143 133 144 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 145 135 146 G4BetheHeitler5DModel::~G4BetheHeitler5DModel( << 136 G4BetheHeitler5DModel::~G4BetheHeitler5DModel() >> 137 {} 147 138 148 //....oooOO0OOooo........oooOO0OOooo........oo 139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 140 150 void G4BetheHeitler5DModel::Initialise(const G 141 void G4BetheHeitler5DModel::Initialise(const G4ParticleDefinition* part, 151 const G4DataVector& vec) 142 const G4DataVector& vec) 152 { 143 { 153 G4PairProductionRelModel::Initialise(part, v << 144 G4BetheHeitlerModel::Initialise(part, vec); 154 145 155 G4EmParameters* theManager = G4EmParameters: 146 G4EmParameters* theManager = G4EmParameters::Instance(); 156 // place to initialise model parameters 147 // place to initialise model parameters 157 // Verbosity levels: ( Can redefine as neede << 158 // 0 = nothing << 159 // > 2 print results << 160 // > 3 print rejection warning from transfor << 161 // > 4 print photon direction & polarisation << 162 fVerbose = theManager->Verbose(); 148 fVerbose = theManager->Verbose(); 163 fConversionType = theManager->GetConversionT << 149 fConversionType = theManager->GetConversionType(); 164 //////////////////////////////////////////// 150 ////////////////////////////////////////////////////////////// 165 // iraw : 151 // iraw : 166 // true : isolated electron or nucleus 152 // true : isolated electron or nucleus. 167 // false : inside atom -> screening for 153 // false : inside atom -> screening form factor 168 iraw = theManager->OnIsolated(); 154 iraw = theManager->OnIsolated(); 169 // G4cout << "BH5DModel::Initialise verbose 155 // G4cout << "BH5DModel::Initialise verbose " << fVerbose 170 // << " isolated " << iraw << " ctype "<< 156 // << " isolated " << iraw << " ctype "<< fConversionType << G4endl; 171 << 172 //Q: Do we need this on Model << 173 // The Leptons defined via SetLeptonPair(..) << 174 SetLowEnergyLimit(2*CLHEP::electron_mass_c2) << 175 } 157 } 176 158 177 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 160 179 void G4BetheHeitler5DModel::SetLeptonPair(cons << 161 180 const G4ParticleDefinition* p2) << 162 // >> 163 // Converting from pair coordinate >> 164 // >> 165 void >> 166 G4BetheHeitler5DModel::BoostG4LorentzVector(const G4LorentzVector& p, >> 167 const G4LorentzVector& q, >> 168 G4LorentzVector& res) const 181 { 169 { 182 G4int pdg1 = p1->GetPDGEncoding(); << 170 // p : 4-vector which will be boosted 183 G4int pdg2 = p2->GetPDGEncoding(); << 171 // q : 4-vector of new origin in the old coordinates 184 G4int pdg = std::abs(pdg1); << 172 const G4double pq = p.x()*q.x() + p.y()*q.y() + p.z()*q.z(); 185 if ( pdg1 != -pdg2 || (pdg != 11 && pdg != 1 << 173 const G4double qq = q.x()*q.x() + q.y()*q.y() + q.z()*q.z(); 186 G4ExceptionDescription ed; << 174 const G4double mass2 = q.t()*q.t()-qq; 187 ed << " Wrong pair of leptons: " << p1->Ge << 175 if ( mass2 > 0.0 ) { 188 << " and " << p1->GetParticleName(); << 176 const G4double mass = std::sqrt(q.t()*q.t()-qq); 189 G4Exception("G4BetheHeitler5DModel::SetLep << 177 const G4double lf = ((q.t()-mass)*pq/qq+p.t())/mass; 190 FatalErrorInArgument, ed, ""); << 178 res.set( (p.x()+q.x()*lf), (p.y()+q.y()*lf), (p.z()+q.z()*lf), >> 179 ((p.t()*q.t()+pq)/mass) ); 191 } else { 180 } else { 192 if ( pdg == 11 ) { << 181 res = p; 193 SetConversionMode(kEPair); << 182 if ( fVerbose > 3 ) { 194 if( pdg1 == 11 ) { << 183 G4cout << "G4BetheHeitler5DModel::BoostG4LorentzVector Warning point not converted" 195 fLepton1 = p1; << 184 << G4endl << "secondary in " << p 196 fLepton2 = p2; << 185 << G4endl << "Pair1 " << q << G4endl; 197 } else { << 198 fLepton1 = p2; << 199 fLepton2 = p1; << 200 } << 201 if (fVerbose > 0) << 202 G4cout << "G4BetheHeitler5DModel::SetLeptonP << 203 << G4endl; << 204 } else { << 205 SetConversionMode(kMuPair); << 206 if( pdg1 == 13 ) { << 207 fLepton1 = p1; << 208 fLepton2 = p2; << 209 } else { << 210 fLepton1 = p2; << 211 fLepton2 = p1; << 212 } << 213 fTheMuPlus = fLepton2; << 214 fTheMuMinus= fLepton1; << 215 if (fVerbose > 0) << 216 G4cout << "G4BetheHeitler5DModel::SetLeptonP << 217 << G4endl; << 218 } 186 } 219 } 187 } 220 } 188 } 221 189 >> 190 // assuming that q.x=q.y=0.0 >> 191 void >> 192 G4BetheHeitler5DModel::BoostG4LorentzVector(const G4LorentzVector& p, >> 193 const G4double qz, >> 194 const G4double qt, >> 195 const G4double lffac, >> 196 const G4double imass, >> 197 G4LorentzVector& res) const >> 198 { >> 199 // p : 4-vector which will be boosted >> 200 // q : 4-vector of new origin in the old coordinates >> 201 const G4double pq = p.z()*qz; >> 202 const G4double lf = (lffac*pq+p.t())*imass; >> 203 res.setZ(p.z()+qz*lf); >> 204 res.setT((p.t()*qt+pq)*imass); >> 205 } >> 206 222 //....oooOO0OOooo........oooOO0OOooo........oo 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 223 208 224 G4double G4BetheHeitler5DModel::MaxDiffCrossSe << 209 G4double G4BetheHeitler5DModel::MaxDiffCrossSection(const G4double* par, 225 << 210 G4double Z, 226 << 211 G4double e, 227 212 G4double loge) const 228 { 213 { 229 const G4double Q = e/par[9]; 214 const G4double Q = e/par[9]; 230 return par[0] * G4Exp((par[2]+loge*par[4])*l 215 return par[0] * G4Exp((par[2]+loge*par[4])*loge) 231 / (par[1]+ G4Exp(par[3]*loge)+G4Exp(p 216 / (par[1]+ G4Exp(par[3]*loge)+G4Exp(par[5]*loge)) 232 * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/ 217 * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/(1+Q)); 233 } 218 } 234 219 235 //....oooOO0OOooo........oooOO0OOooo........oo 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 236 221 237 void << 222 void 238 G4BetheHeitler5DModel::SampleSecondaries(std:: 223 G4BetheHeitler5DModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 239 const 224 const G4MaterialCutsCouple* couple, 240 const 225 const G4DynamicParticle* aDynamicGamma, 241 G4dou 226 G4double, G4double) 242 { 227 { 243 // MeV 228 // MeV 244 static const G4double ElectronMass = CLHEP 229 static const G4double ElectronMass = CLHEP::electron_mass_c2; 245 << 230 static const G4double ElectronMass2 = ElectronMass*ElectronMass; 246 const G4double LeptonMass = fLepton1->GetPDG << 247 const G4double LeptonMass2 = LeptonMass*Lep << 248 << 249 static const G4double alpha0 = CLHEP 231 static const G4double alpha0 = CLHEP::fine_structure_const; 250 // mm << 232 // mm 251 static const G4double r0 = CLHEP 233 static const G4double r0 = CLHEP::classic_electr_radius; 252 // mbarn 234 // mbarn 253 static const G4double r02 = r0*r0 235 static const G4double r02 = r0*r0*1.e+25; 254 static const G4double twoPi = CLHEP 236 static const G4double twoPi = CLHEP::twopi; 255 static const G4double factor = alpha 237 static const G4double factor = alpha0 * r02 / (twoPi*twoPi); 256 // static const G4double factor1 = p 238 // static const G4double factor1 = pow((6.0 * pi),(1.0/3.0))/(8.*alpha0*ElectronMass); 257 static const G4double factor1 = 2.661 239 static const G4double factor1 = 2.66134007899/(8.*alpha0*ElectronMass); 258 // 240 // 259 G4double PairInvMassMin = 2.*LeptonMass; << 241 static const G4double PairInvMassMin = 2.*ElectronMass; 260 G4double TrThreshold = 2.0 * ( (LeptonMass2 << 261 << 262 // 242 // 263 static const G4double nu[2][10] = { << 243 static const G4double nu[10] = { 0.0227436, 0.0582046, 3.0322675, 2.8275065, 264 //electron << 244 -0.0034004, 1.1212766, 1.8989468, 68.3492750, 265 { 0.0227436, 0.0582046, 3.0322675, 2.82750 << 245 0.0211186, 14.4 }; 266 1.1212766, 1.8989468, 68.3492750, 0.0211 << 246 static const G4double tr[10] = { 0.0332350, 4.3942537, 2.8515925, 2.6351695, 267 //muon << 247 -0.0031510, 1.5737305, 1.8104647, 20.6434021, 268 {0.67810E-06, 0.86037E+05, 2.0008395, 1.67 << 248 -0.0272586, 28.9}; 269 1.4222, 0.0, 263230.0, 0.0521, 51.1338} << 249 // 270 }; << 250 static const G4double para[3][2] = { {11., -16.},{-1.17, -2.95},{-2., -0.5} }; 271 static const G4double tr[2][10] = { << 272 //electron << 273 { 0.0332350, 4.3942537, 2.8515925, 2.6351 << 274 1.5737305, 1.8104647, 20.6434021, -0.027 << 275 //muon << 276 {0.10382E-03, 0.14408E+17, 4.1368679, 3.26 << 277 0.0000, 0.0, 0.0, 0.0000, 1.0000} << 278 }; << 279 // << 280 static const G4double para[2][3][2] = { << 281 //electron << 282 { {11., -16.},{-1.17, -2.95},{-2., -0.5} } << 283 //muon << 284 { {17.5, 1.},{-1.17, -2.95},{2., 6.} } << 285 }; << 286 // 251 // 287 static const G4double correctionIndex = 1.4; 252 static const G4double correctionIndex = 1.4; 288 // 253 // 289 const G4double GammaEnergy = aDynamicGamma- 254 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 290 // Protection, Will not be true tot cross se << 291 if ( GammaEnergy <= PairInvMassMin) { return << 292 << 293 const G4double GammaEnergy2 = GammaEnergy*Ga 255 const G4double GammaEnergy2 = GammaEnergy*GammaEnergy; 294 << 256 // Will not be true tot cross section = 0 295 //////////////////////////////////////////// << 257 if ( GammaEnergy <= 2.0*ElectronMass) { return; } 296 const G4ParticleMomentum GammaDirection = << 258 // 297 aDynamicGamma->GetMomentumDirection(); << 259 const G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection(); 298 G4ThreeVector GammaPolarization = aDynamicGa 260 G4ThreeVector GammaPolarization = aDynamicGamma->GetPolarization(); 299 261 300 // The protection polarization perpendicular 262 // The protection polarization perpendicular to the direction vector, 301 // as it done in G4LivermorePolarizedGammaCo 263 // as it done in G4LivermorePolarizedGammaConversionModel, 302 // assuming Direction is unitary vector 264 // assuming Direction is unitary vector 303 // (projection to plane) p_proj = p - (p o 265 // (projection to plane) p_proj = p - (p o d)/(d o d) x d 304 if ( GammaPolarization.howOrthogonal(GammaDi 266 if ( GammaPolarization.howOrthogonal(GammaDirection) != 0) { 305 GammaPolarization -= GammaPolarization.dot 267 GammaPolarization -= GammaPolarization.dot(GammaDirection) * GammaDirection; 306 } 268 } 307 // End of Protection 269 // End of Protection 308 // 270 // 309 const G4double GammaPolarizationMag = GammaP 271 const G4double GammaPolarizationMag = GammaPolarization.mag(); 310 << 311 //////////////////////////////////////////// 272 ////////////////////////////////////////////////////////////// 312 // target element 273 // target element 313 // select randomly one element constituting 274 // select randomly one element constituting the material 314 const G4Element* anElement = SelectTargetAt << 275 const G4Element* anElement = SelectRandomAtom(couple, fTheGamma, GammaEnergy); 315 aDyna << 316 // Atomic number 276 // Atomic number 317 const G4int Z = anElement->GetZasInt() 277 const G4int Z = anElement->GetZasInt(); 318 const G4int A = SelectIsotopeNumber(an 278 const G4int A = SelectIsotopeNumber(anElement); 319 const G4double iZ13 = 1./anElement->GetIonis 279 const G4double iZ13 = 1./anElement->GetIonisation()->GetZ3(); 320 const G4double targetMass = G4NucleiProperti 280 const G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z); 321 281 322 const G4double NuThreshold = 2.0 * ( (Lept << 323 // No conversion possible below nuclear thre << 324 if ( GammaEnergy <= NuThreshold) { return; } << 325 << 326 CLHEP::HepRandomEngine* rndmEngine = G4Rando 282 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 327 283 328 // itriplet : true -- triplet, false -- nucl 284 // itriplet : true -- triplet, false -- nuclear. 329 G4bool itriplet = false; 285 G4bool itriplet = false; 330 if (fConversionType == 1) { 286 if (fConversionType == 1) { 331 itriplet = false; 287 itriplet = false; 332 } else if (fConversionType == 2) { 288 } else if (fConversionType == 2) { 333 itriplet = true; 289 itriplet = true; 334 if ( GammaEnergy <= TrThreshold ) return; << 290 if ( GammaEnergy <= 4.0*ElectronMass ) return; 335 } else if ( GammaEnergy > TrThreshold ) { << 291 } else if ( GammaEnergy > 4.0*ElectronMass ) { 336 // choose triplet or nuclear from a triple 292 // choose triplet or nuclear from a triplet/nuclear=1/Z 337 // total cross section ratio. 293 // total cross section ratio. 338 // approximate at low energies ! 294 // approximate at low energies ! 339 if(rndmEngine->flat()*(Z+1) < 1.) { 295 if(rndmEngine->flat()*(Z+1) < 1.) { 340 itriplet = true; 296 itriplet = true; 341 } 297 } 342 } 298 } 343 << 344 // 299 // 345 const G4double RecoilMass = itriplet ? Elec 300 const G4double RecoilMass = itriplet ? ElectronMass : targetMass; 346 const G4double RecoilMass2 = RecoilMass*Reco 301 const G4double RecoilMass2 = RecoilMass*RecoilMass; 347 const G4double sCMS = 2.*RecoilMass*G 302 const G4double sCMS = 2.*RecoilMass*GammaEnergy + RecoilMass2; 348 const G4double sCMSPlusRM2 = sCMS + RecoilMa 303 const G4double sCMSPlusRM2 = sCMS + RecoilMass2; 349 const G4double sqrts = std::sqrt(sCMS) 304 const G4double sqrts = std::sqrt(sCMS); 350 const G4double isqrts2 = 1./(2.*sqrts); 305 const G4double isqrts2 = 1./(2.*sqrts); 351 // 306 // 352 const G4double PairInvMassMax = sqrts-Reco 307 const G4double PairInvMassMax = sqrts-RecoilMass; 353 const G4double PairInvMassRange = PairInvMas 308 const G4double PairInvMassRange = PairInvMassMax/PairInvMassMin; 354 const G4double lnPairInvMassRange = G4Log(Pa 309 const G4double lnPairInvMassRange = G4Log(PairInvMassRange); 355 310 356 // initial state. Defines z axis of "0" fram 311 // initial state. Defines z axis of "0" frame as along photon propagation. 357 // Since CMS(0., 0., GammaEnergy, GammaEnerg << 312 // create 4-vectors: gamma0 + target0 and CMS=gamma0+target0 358 const G4double betaCMS = G4LorentzVector(0.0 << 313 // Since CMS(0., 0., GammaEnergy, GammaEnergy+RecoilMass) set some constants 359 << 314 // for the special boost that makes use of the form of CMS 4-vector >> 315 const G4double CMSqz = GammaEnergy; >> 316 const G4double CMSt = GammaEnergy+RecoilMass; >> 317 const G4double iCMSmass = 1./std::sqrt(RecoilMass*(RecoilMass+2.*GammaEnergy)); >> 318 const G4double CMSfact = (CMSt-1./iCMSmass)/(CMSqz*CMSqz); 360 // maximum value of pdf 319 // maximum value of pdf 361 const G4double EffectiveZ = iraw ? 0.5 : Z; 320 const G4double EffectiveZ = iraw ? 0.5 : Z; 362 const G4double Threshold = itriplet ? TrThr << 321 const G4double Threshold = itriplet ? 4.*ElectronMass : 2.*ElectronMass; 363 const G4double AvailableEnergy = GammaEne 322 const G4double AvailableEnergy = GammaEnergy - Threshold; 364 const G4double LogAvailableEnergy = G4Log(Av 323 const G4double LogAvailableEnergy = G4Log(AvailableEnergy); 365 // 324 // 366 const G4double MaxDiffCross = itriplet 325 const G4double MaxDiffCross = itriplet 367 ? MaxDiffCrossSection(tr[fConvMode], << 326 ? MaxDiffCrossSection(tr, EffectiveZ, AvailableEnergy, LogAvailableEnergy) 368 EffectiveZ, AvailableEnergy, LogAvaila << 327 : MaxDiffCrossSection(nu, EffectiveZ, AvailableEnergy, LogAvailableEnergy); 369 : MaxDiffCrossSection(nu[fConvMode], << 370 EffectiveZ, AvailableEnergy, LogAva << 371 // 328 // 372 // 50% safety marging factor 329 // 50% safety marging factor 373 const G4double ymax = 1.5 * MaxDiffCross; 330 const G4double ymax = 1.5 * MaxDiffCross; 374 // x1 bounds 331 // x1 bounds 375 const G4double xu1 = (LogAvailableEnergy > << 332 const G4double xu1 = (LogAvailableEnergy > para[2][0]) 376 ? para[fConvMode][0][0] + << 333 ? para[0][0] + para[1][0]*LogAvailableEnergy 377 para[fConvMode][1][0]*LogAvailableEner << 334 : para[0][0] + para[2][0]*para[1][0]; 378 : para[fConvMode][0][0] << 335 const G4double xl1 = (LogAvailableEnergy > para[2][1]) 379 para[fConvMode][2][0]*para[fConvMode][ << 336 ? para[0][1] + para[1][1]*LogAvailableEnergy 380 const G4double xl1 = (LogAvailableEnergy > << 337 : para[0][1] + para[2][1]*para[1][1]; 381 ? para[fConvMode][0][1] << 338 // 382 para[fConvMode][1][1]*LogAvailableEner << 339 G4LorentzVector Recoil0; 383 : para[fConvMode][0][1] << 340 G4LorentzVector Positron0; 384 para[fConvMode][2][1]*para[fConvMode][ << 341 G4LorentzVector Electron0; 385 // << 386 G4LorentzVector Recoil; << 387 G4LorentzVector LeptonPlus; << 388 G4LorentzVector LeptonMinus; << 389 G4double pdf = 0.; 342 G4double pdf = 0.; 390 << 391 G4double rndmv6[6] = {0.0}; << 392 const G4double corrFac = 1.0/(correctionInde << 393 const G4double expLowLim = -20.; << 394 const G4double logLowLim = G4Exp(expLowLim/c << 395 G4double z0, z1, z2, x0, x1; << 396 G4double betheheitler, sinTheta, cosTheta, d << 397 // START Sampling 343 // START Sampling 398 do { 344 do { 399 << 345 ////////////////////////////////////////////////// 400 rndmEngine->flatArray(6, rndmv6); << 401 << 402 ////////////////////////////////////////// << 403 // pdf pow(x,c) with c = 1.4 346 // pdf pow(x,c) with c = 1.4 404 // integral y = pow(x,(c+1))/(c+1) @ x = 1 347 // integral y = pow(x,(c+1))/(c+1) @ x = 1 => y = 1 /(1+c) 405 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 348 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 1.0 ) */ ) /( c + 1.0) ) 406 ////////////////////////////////////////// << 349 ////////////////////////////////////////////////// 407 << 350 const G4double X1 = 408 z0 = (rndmv6[0] > logLowLim) ? G4Log(rndmv << 351 G4Exp(G4Log(rndmEngine->flat())/(correctionIndex + 1.0)); 409 G4double X1 = (z0 > expLowLim) ? G4Exp(z0) << 352 410 z1 = xl1 + (xu1 - xl1)*rndmv6[1]; << 353 const G4double x0 = G4Exp(xl1 + (xu1 - xl1)*rndmEngine->flat()); 411 if (z1 > expLowLim) { << 354 const G4double dum0 = 1./(1.+x0); 412 x0 = G4Exp(z1); << 355 const G4double cosTheta = (x0-1.)*dum0; 413 dum0 = 1.0/(1.0 + x0); << 356 const G4double sinTheta = std::sqrt(4.*x0)*dum0; 414 x1 = dum0*x0; << 357 415 cosTheta = -1.0 + 2.0*x1; << 358 const G4double PairInvMass = PairInvMassMin*G4Exp(X1*X1*lnPairInvMassRange); 416 sinTheta = 2*std::sqrt(x1*(1.0 - x1)); << 359 417 } else { << 360 G4double rndmv3[3]; 418 x0 = 0.0; << 361 rndmEngine->flatArray(3, rndmv3); 419 dum0 = 1.0; << 362 //-------------------------------------------------------------------------- 420 cosTheta = -1.0; << 363 // const G4double ThetaLept = pi*rndmv3[0]; 421 sinTheta = 0.0; << 364 // const G4double cosThetaLept = std::cos(ThetaLept); 422 } << 365 // const G4double sinThetaLept = std::sin(ThetaLept); 423 << 366 // 424 z2 = X1*X1*lnPairInvMassRange; << 367 // const G4double PhiLept = twoPi*rndmv3[1]-pi; 425 const G4double PairInvMass = PairInvMassMi << 368 // const G4double cosPhiLept = std::cos(PhiLept); 426 << 369 // const G4double sinPhiLept = std::sin(PhiLept); >> 370 // >> 371 // const G4double Phi = twoPi*rndmv3[2]-pi; >> 372 // const G4double cosPhi = std::cos(Phi); >> 373 // const G4double sinPhi = std::sin(Phi); >> 374 //--------------------------------------------------------------------------- 427 // cos and sin theta-lepton 375 // cos and sin theta-lepton 428 const G4double cosThetaLept = std::cos(pi* << 376 const G4double cosThetaLept = std::cos(pi*rndmv3[0]); 429 // sin(ThetaLept) is always in [0,+1] if T 377 // sin(ThetaLept) is always in [0,+1] if ThetaLept is in [0,pi] 430 const G4double sinThetaLept = std::sqrt((1 << 378 const G4double sinThetaLept = std::sqrt((1.-cosThetaLept)*(1.+cosThetaLept)); 431 // cos and sin phi-lepton 379 // cos and sin phi-lepton 432 const G4double cosPhiLept = std::cos(two << 380 const G4double cosPhiLept = std::cos(twoPi*rndmv3[1]-pi); 433 // sin(PhiLept) is in [-1,0] if PhiLept in 381 // sin(PhiLept) is in [-1,0] if PhiLept in [-pi,0) and 434 // is in [0,+1] if PhiLept in 382 // is in [0,+1] if PhiLept in [0,+pi] 435 const G4double sinPhiLept = std::copysig << 383 const G4double sinPhiLept = std::copysign(std::sqrt((1.-cosPhiLept)*(1.+cosPhiLept)),rndmv3[1]-0.5); 436 // cos and sin phi 384 // cos and sin phi 437 const G4double cosPhi = std::cos(two << 385 const G4double cosPhi = std::cos(twoPi*rndmv3[2]-pi); 438 const G4double sinPhi = std::copysig << 386 const G4double sinPhi = std::copysign(std::sqrt((1.-cosPhi)*(1.+cosPhi)),rndmv3[2]-0.5); 439 387 440 ////////////////////////////////////////// << 441 // frames: 388 // frames: 442 // 3 : the laboratory Lorentz frame, Geant << 443 // 0 : the laboratory Lorentz frame, axes 389 // 0 : the laboratory Lorentz frame, axes along photon direction and polarisation 444 // 1 : the center-of-mass Lorentz frame 390 // 1 : the center-of-mass Lorentz frame 445 // 2 : the pair Lorentz frame 391 // 2 : the pair Lorentz frame 446 ////////////////////////////////////////// << 392 // 3 : the laboratory Lorentz frame, Geant4 axes definition 447 393 448 // in the center-of-mass frame 394 // in the center-of-mass frame 449 << 450 const G4double RecEnergyCMS = (sCMSPlusRM 395 const G4double RecEnergyCMS = (sCMSPlusRM2-PairInvMass*PairInvMass)*isqrts2; 451 const G4double LeptonEnergy2 = PairInvMass 396 const G4double LeptonEnergy2 = PairInvMass*0.5; 452 << 397 // Denis ** correction 453 // New way of calucaltion thePRecoil to av << 398 // const G4double thePRecoil = std::sqrt( (RecEnergyCMS-RecoilMass) 454 G4double abp = std::max((2.0*GammaEnergy*R << 399 // *(RecEnergyCMS+RecoilMass)); 455 PairInvMass*PairInvMass + 2.0*PairI << 400 const G4double ap1 = 2.0*GammaEnergy*RecoilMass - 456 (2.0*GammaEnergy*R << 401 PairInvMass*PairInvMass + 2.0*PairInvMass*RecoilMass; 457 PairInvMass*PairInvMass - 2.0*PairI << 402 const G4double bp1 = 2.0*GammaEnergy*RecoilMass - 458 << 403 PairInvMass*PairInvMass - 2.0*PairInvMass*RecoilMass; 459 G4double thePRecoil = std::sqrt(abp) * isq << 404 >> 405 if (bp1 <= 0.0 ) { >> 406 if ( fVerbose > 3 ) { >> 407 G4cout >> 408 << "G4BetheHeitler5DModel::SampleSecondaries Warning bp1 " >> 409 << bp1 << "point rejected" << G4endl >> 410 << "GammaEnergy " << GammaEnergy << G4endl >> 411 << "PairInvMass " << PairInvMass << G4endl; >> 412 } >> 413 pdf = -1.0; // force next iteration >> 414 continue; >> 415 } >> 416 const G4double thePRecoil = std::sqrt(ap1 * bp1) * isqrts2; >> 417 // Denis ** correction 460 418 461 // back to the center-of-mass frame 419 // back to the center-of-mass frame 462 Recoil.set( thePRecoil*sinTheta*cosPhi, << 420 const G4LorentzVector Recoil1( thePRecoil*sinTheta*cosPhi, 463 thePRecoil*sinTheta*sinPhi, 421 thePRecoil*sinTheta*sinPhi, 464 thePRecoil*cosTheta, 422 thePRecoil*cosTheta, 465 RecEnergyCMS); 423 RecEnergyCMS); 466 424 467 // in the pair frame << 425 const G4LorentzVector Pair1(-Recoil1.x(), 468 const G4double thePLepton = std::sqrt( << 426 -Recoil1.y(), 469 * << 427 -Recoil1.z(), >> 428 sqrts-RecEnergyCMS); 470 429 471 LeptonPlus.set(thePLepton*sinThetaLept*cos << 430 // in the pair frame 472 thePLepton*sinThetaLept*sinPhiLept, << 431 const G4double thePLepton = std::sqrt( (LeptonEnergy2-ElectronMass) 473 thePLepton*cosThetaLept, << 432 *(LeptonEnergy2+ElectronMass)); 474 LeptonEnergy2); << 433 475 << 434 const G4LorentzVector Positron2( thePLepton*sinThetaLept*cosPhiLept, 476 LeptonMinus.set(-LeptonPlus.x(), << 435 thePLepton*sinThetaLept*sinPhiLept, 477 -LeptonPlus.y(), << 436 thePLepton*cosThetaLept, 478 -LeptonPlus.z(), << 437 LeptonEnergy2); 479 LeptonEnergy2); << 438 >> 439 const G4LorentzVector Electron2(-Positron2.x(), >> 440 -Positron2.y(), >> 441 -Positron2.z(), >> 442 LeptonEnergy2); 480 443 481 444 482 // Normalisation of final state phase spac 445 // Normalisation of final state phase space: 483 // Section 47 of Particle Data Group, Chin 446 // Section 47 of Particle Data Group, Chin. Phys. C, 40, 100001 (2016) 484 // const G4double Norme = Recoil1.vect( << 447 const G4double Norme = Recoil1.vect().mag() * Positron2.vect().mag(); 485 const G4double Norme = Recoil.vect().mag() << 448 // 486 << 449 G4LorentzVector Positron1; 487 // e+, e- to CMS frame from pair frame << 450 G4LorentzVector Electron1; 488 << 451 BoostG4LorentzVector(Positron2, Pair1, Positron1); 489 // boost vector from Pair to CMS << 452 BoostG4LorentzVector(Electron2, Pair1, Electron1); 490 const G4ThreeVector pair2cms = << 453 491 G4LorentzVector( -Recoil.x(), -Recoil.y(), << 492 sqrts-RecEnergyCMS).boostVector(); << 493 << 494 LeptonPlus.boost(pair2cms); << 495 LeptonMinus.boost(pair2cms); << 496 << 497 // back to the laboratory frame (make use 454 // back to the laboratory frame (make use of the CMS(0,0,Eg,Eg+RM)) form 498 << 455 Recoil0.setX(Recoil1.x()); 499 Recoil.boostZ(betaCMS); << 456 Recoil0.setY(Recoil1.y()); 500 LeptonPlus.boostZ(betaCMS); << 457 BoostG4LorentzVector(Recoil1 , CMSqz, CMSt, CMSfact, iCMSmass, Recoil0); 501 LeptonMinus.boostZ(betaCMS); << 458 >> 459 Positron0.setX(Positron1.x()); >> 460 Positron0.setY(Positron1.y()); >> 461 BoostG4LorentzVector(Positron1, CMSqz, CMSt, CMSfact, iCMSmass, Positron0); >> 462 >> 463 Electron0.setX(Electron1.x()); >> 464 Electron0.setY(Electron1.y()); >> 465 BoostG4LorentzVector(Electron1, CMSqz, CMSt, CMSfact, iCMSmass, Electron0); 502 466 503 // Jacobian factors 467 // Jacobian factors 504 const G4double Jacob0 = x0*dum0*dum0; 468 const G4double Jacob0 = x0*dum0*dum0; 505 const G4double Jacob1 = 2.*X1*lnPairInvMas 469 const G4double Jacob1 = 2.*X1*lnPairInvMassRange*PairInvMass; 506 const G4double Jacob2 = std::abs(sinThetaL 470 const G4double Jacob2 = std::abs(sinThetaLept); 507 471 508 // there is no probability to have a lepto << 472 const G4double EPlus = Positron0.t(); 509 // X and Y components of momentum may be z << 473 const G4double PPlus = Positron0.vect().mag(); 510 const G4double EPlus = LeptonPlus.t(); << 474 const G4double sinThetaPlus = Positron0.vect().perp()/PPlus; 511 const G4double PPlus = LeptonPlus.vect().m << 475 const G4double cosThetaPlus = Positron0.vect().cosTheta(); 512 const G4double pPX = LeptonPlus.x(); << 476 513 const G4double pPY = LeptonPlus.y(); << 477 const G4double pPX = Positron0.x(); 514 const G4double pPZ = LeptonPlus.z(); << 478 const G4double pPY = Positron0.y(); 515 G4double sinPhiPlus = 1.0; << 479 const G4double dum1 = 1./std::sqrt( pPX*pPX + pPY*pPY ); 516 G4double cosPhiPlus = 0.0; << 480 const G4double cosPhiPlus = pPX*dum1; 517 G4double sinThetaPlus = 0.0; << 481 const G4double sinPhiPlus = pPY*dum1; 518 G4double cosThetaPlus = pPZ/PPlus; << 519 if (cosThetaPlus < 1.0 && cosThetaPlus > - << 520 sinThetaPlus = std::sqrt((1.0 - cosTheta << 521 sinPhiPlus = pPY/(PPlus*sinThetaPlus); << 522 cosPhiPlus = pPX/(PPlus*sinThetaPlus); << 523 } << 524 482 525 // denominators: 483 // denominators: 526 // the two cancelling leading terms for fo 484 // the two cancelling leading terms for forward emission at high energy, removed 527 const G4double elMassCTP = LeptonMass*cosT << 485 const G4double elMassCTP = ElectronMass*cosThetaPlus; 528 const G4double ePlusSTP = EPlus*sinThetaP << 486 const G4double ePlusSTP = EPlus*sinThetaPlus; 529 const G4double DPlus = (elMassCTP*elMa 487 const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP) 530 /(EPlus + PPlus* 488 /(EPlus + PPlus*cosThetaPlus); 531 489 532 // there is no probability to have a lepto << 490 const G4double EMinus = Electron0.t(); 533 // X and Y components of momentum may be z << 491 const G4double PMinus = Electron0.vect().mag(); 534 const G4double EMinus = LeptonMinus.t(); << 492 const G4double sinThetaMinus = Electron0.vect().perp()/PMinus; 535 const G4double PMinus = LeptonMinus.vect() << 493 const G4double cosThetaMinus = Electron0.vect().cosTheta(); 536 const G4double ePX = LeptonMinus.x(); << 494 537 const G4double ePY = LeptonMinus.y(); << 495 const G4double ePX = Electron0.x(); 538 const G4double ePZ = LeptonMinus.z(); << 496 const G4double ePY = Electron0.y(); 539 G4double sinPhiMinus = 0.0; << 497 const G4double dum2 = 1./std::sqrt( ePX*ePX + ePY*ePY ); 540 G4double cosPhiMinus = 1.0; << 498 const G4double cosPhiMinus = ePX*dum2; 541 G4double sinThetaMinus = 0.0; << 499 const G4double sinPhiMinus = ePY*dum2; 542 G4double cosThetaMinus = ePZ/PMinus; << 500 543 if (cosThetaMinus < 1.0 && cosThetaMinus > << 501 const G4double elMassCTM = ElectronMass*cosThetaMinus; 544 sinThetaMinus = std::sqrt((1.0 - cosThet << 502 const G4double eMinSTM = EMinus*sinThetaMinus; 545 sinPhiMinus = ePY/(PMinus*sinThetaMinus) << 503 const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM) 546 cosPhiMinus = ePX/(PMinus*sinThetaMinus) << 547 } << 548 << 549 const G4double elMassCTM = LeptonMass*cosT << 550 const G4double eMinSTM = EMinus*sinTheta << 551 const G4double DMinus = (elMassCTM*elMa << 552 /(EMinus + PMinu 504 /(EMinus + PMinus*cosThetaMinus); 553 505 554 // cos(phiMinus-PhiPlus) 506 // cos(phiMinus-PhiPlus) 555 const G4double cosdPhi = cosPhiPlus*cosPhi 507 const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus; 556 const G4double PRec = Recoil.vect().mag << 508 const G4double PRec = Recoil0.vect().mag(); 557 const G4double q2 = PRec*PRec; 509 const G4double q2 = PRec*PRec; 558 << 510 const G4double BigPhi = -ElectronMass2 / (GammaEnergy*GammaEnergy2 * q2*q2); 559 const G4double BigPhi = -LeptonMass2 / (G << 560 511 561 G4double FormFactor = 1.; 512 G4double FormFactor = 1.; 562 if (!iraw) { 513 if (!iraw) { 563 if (itriplet) { 514 if (itriplet) { 564 const G4double qun = factor1*iZ13*iZ13; 515 const G4double qun = factor1*iZ13*iZ13; 565 const G4double nun = qun * PRec; 516 const G4double nun = qun * PRec; 566 if (nun < 1.) { 517 if (nun < 1.) { 567 FormFactor = (nun < 0.01) ? (13.8-5 << 518 FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun 568 : std::sq << 519 : std::sqrt(1-(nun-1)*(nun-1)); 569 } // else FormFactor = 1 by default 520 } // else FormFactor = 1 by default 570 } else { 521 } else { 571 const G4double dum3 = 217.*PRec*iZ13; 522 const G4double dum3 = 217.*PRec*iZ13; 572 const G4double AFF = 1./(1. + dum3*dum3); 523 const G4double AFF = 1./(1. + dum3*dum3); 573 FormFactor = (1.-AFF)*(1-AFF); 524 FormFactor = (1.-AFF)*(1-AFF); 574 } 525 } 575 } // else FormFactor = 1 by default 526 } // else FormFactor = 1 by default 576 527 >> 528 G4double betheheitler; 577 if (GammaPolarizationMag==0.) { 529 if (GammaPolarizationMag==0.) { 578 const G4double pPlusSTP = PPlus*sinThe 530 const G4double pPlusSTP = PPlus*sinThetaPlus; 579 const G4double pMinusSTM = PMinus*sinTh 531 const G4double pMinusSTM = PMinus*sinThetaMinus; 580 const G4double pPlusSTPperDP = pPlusSTP 532 const G4double pPlusSTPperDP = pPlusSTP/DPlus; 581 const G4double pMinusSTMperDM = pMinusST 533 const G4double pMinusSTMperDM = pMinusSTM/DMinus; 582 const G4double dunpol = BigPhi*( << 534 const G4double dunpol = BigPhi*( 583 pPlusSTPperDP *pPlusSTPperDP 535 pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2) 584 + pMinusSTMperDM*pMinusSTMperD 536 + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2) 585 + 2.*pPlusSTPperDP*pMinusSTMpe 537 + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi 586 *(4.*EPlus*EMinus + q2 - 2 538 *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2) 587 - 2.*GammaEnergy2*(pPlusSTP*pP 539 - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus)); 588 betheheitler = dunpol * factor; 540 betheheitler = dunpol * factor; 589 } else { 541 } else { 590 const G4double pPlusSTP = PPlus*sinThet 542 const G4double pPlusSTP = PPlus*sinThetaPlus; 591 const G4double pMinusSTM = PMinus*sinThe 543 const G4double pMinusSTM = PMinus*sinThetaMinus; 592 const G4double pPlusSTPCPPperDP = pPlus 544 const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus; 593 const G4double pMinusSTMCPMperDM = pMinu 545 const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus; 594 const G4double caa = 2.*(EPlus*pMinusSTM 546 const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP); 595 const G4double cbb = pMinusSTMCPMperDM-p 547 const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP; 596 const G4double ccc = (pPlusSTP*pPlusSTP << 548 const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM 597 +2.*pPlusSTP*pMinusS 549 +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus); 598 const G4double dtot= 2.*BigPhi*( caa*caa 550 const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc); 599 betheheitler = dtot * factor; 551 betheheitler = dtot * factor; 600 } 552 } 601 // 553 // 602 const G4double cross = Norme * Jacob0 * J << 554 const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler 603 * FormFactor * Recoi 555 * FormFactor * RecoilMass / sqrts; 604 pdf = cross * (xu1 - xl1) / G4Exp(correcti 556 pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1; 605 } while ( pdf < ymax * rndmv6[5] ); << 557 } while ( pdf < ymax * rndmEngine->flat() ); 606 // END of Sampling 558 // END of Sampling 607 << 559 608 if ( fVerbose > 2 ) { 560 if ( fVerbose > 2 ) { 609 G4double recul = std::sqrt(Recoil.x()*Reco << 561 G4double recul = std::sqrt(Recoil0.x()*Recoil0.x()+Recoil0.y()*Recoil0.y() 610 +Recoil.z()*Reco << 562 +Recoil0.z()*Recoil0.z()); 611 G4cout << "BetheHeitler5DModel GammaEnergy 563 G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy 612 << " PDF= " << pdf << " ymax= " << ymax << 564 << " PDF= " << pdf << " ymax= " << ymax 613 << " recul= " << recul << G4endl; 565 << " recul= " << recul << G4endl; 614 } 566 } 615 << 567 616 // back to Geant4 system 568 // back to Geant4 system 617 569 618 if ( fVerbose > 4 ) { 570 if ( fVerbose > 4 ) { 619 G4cout << "BetheHeitler5DModel GammaDirect 571 G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl; 620 G4cout << "BetheHeitler5DModel GammaPolari 572 G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl; 621 G4cout << "BetheHeitler5DModel GammaEnergy 573 G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl; 622 G4cout << "BetheHeitler5DModel Conv " 574 G4cout << "BetheHeitler5DModel Conv " 623 << (itriplet ? "triplet" : "nucl") << G4e 575 << (itriplet ? "triplet" : "nucl") << G4endl; 624 } 576 } 625 577 626 if (GammaPolarizationMag == 0.0) { 578 if (GammaPolarizationMag == 0.0) { 627 // set polarization axis orthohonal to dir 579 // set polarization axis orthohonal to direction 628 GammaPolarization = GammaDirection.orthogo 580 GammaPolarization = GammaDirection.orthogonal().unit(); 629 } else { 581 } else { 630 // GammaPolarization not a unit vector 582 // GammaPolarization not a unit vector 631 GammaPolarization /= GammaPolarizationMag; 583 GammaPolarization /= GammaPolarizationMag; 632 } 584 } 633 585 634 // The unit norm vector that is orthogonal t 586 // The unit norm vector that is orthogonal to the two others 635 G4ThreeVector yGrec = GammaDirection.cross(G 587 G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization); 636 << 588 // rotation 637 // rotation from gamma ref. sys. to World << 589 G4ThreeVector Rot = Recoil0.x()*GammaPolarization + Recoil0.y()*yGrec 638 G4RotationMatrix GtoW(GammaPolarization,yGre << 590 + Recoil0.z()*GammaDirection; 639 << 591 Recoil0.setVect(Rot); 640 Recoil.transform(GtoW); << 592 Rot = Positron0.x()*GammaPolarization + Positron0.y()*yGrec 641 LeptonPlus.transform(GtoW); << 593 + Positron0.z()*GammaDirection; 642 LeptonMinus.transform(GtoW); << 594 Positron0.setVect(Rot); 643 << 595 Rot = Electron0.x()*GammaPolarization + Electron0.y()*yGrec >> 596 + Electron0.z()*GammaDirection; >> 597 Electron0.setVect(Rot); >> 598 // 644 if ( fVerbose > 2 ) { 599 if ( fVerbose > 2 ) { 645 G4cout << "BetheHeitler5DModel Recoil " << << 600 G4cout << "BetheHeitler5DModel Recoil0 " << Recoil0.x() << " " << Recoil0.y() << " " << Recoil0.z() 646 << " " << Recoil.t() << " " << G4endl; << 601 << " " << Recoil0.t() << " " << G4endl; 647 G4cout << "BetheHeitler5DModel LeptonPlus << 602 G4cout << "BetheHeitler5DModel Positron0 " << Positron0.x() << " " << Positron0.y() << " " 648 << LeptonPlus.z() << " " << LeptonPlus.t( << 603 << Positron0.z() << " " << Positron0.t() << " " << G4endl; 649 G4cout << "BetheHeitler5DModel LeptonMinus << 604 G4cout << "BetheHeitler5DModel Electron0 " << Electron0.x() << " " << Electron0.y() << " " 650 << LeptonMinus.z() << " " << LeptonMinus. << 605 << Electron0.z() << " " << Electron0.t() << " " << G4endl; 651 } 606 } 652 << 607 653 // Create secondaries 608 // Create secondaries 654 auto aParticle1 = new G4DynamicParticle(fLep << 609 655 auto aParticle2 = new G4DynamicParticle(fLep << 610 // electron 656 << 611 G4DynamicParticle* aParticle1 = new G4DynamicParticle(fTheElectron,Electron0); >> 612 // positron >> 613 G4DynamicParticle* aParticle2 = new G4DynamicParticle(fThePositron,Positron0); 657 // create G4DynamicParticle object for the p 614 // create G4DynamicParticle object for the particle3 ( recoil ) 658 G4ParticleDefinition* RecoilPart; 615 G4ParticleDefinition* RecoilPart; 659 if (itriplet) { 616 if (itriplet) { 660 // triplet 617 // triplet 661 RecoilPart = fTheElectron; 618 RecoilPart = fTheElectron; 662 } else{ 619 } else{ 663 RecoilPart = theIonTable->GetIon(Z, A, 0); 620 RecoilPart = theIonTable->GetIon(Z, A, 0); 664 } 621 } 665 auto aParticle3 = new G4DynamicParticle(Reco << 622 G4DynamicParticle* aParticle3 = new G4DynamicParticle(RecoilPart,Recoil0); 666 << 623 667 // Fill output vector 624 // Fill output vector 668 fvect->push_back(aParticle1); 625 fvect->push_back(aParticle1); 669 fvect->push_back(aParticle2); 626 fvect->push_back(aParticle2); 670 fvect->push_back(aParticle3); 627 fvect->push_back(aParticle3); 671 628 672 // kill incident photon 629 // kill incident photon 673 fParticleChange->SetProposedKineticEnergy(0. 630 fParticleChange->SetProposedKineticEnergy(0.); 674 fParticleChange->ProposeTrackStatus(fStopAnd 631 fParticleChange->ProposeTrackStatus(fStopAndKill); 675 } 632 } 676 633 677 //....oooOO0OOooo........oooOO0OOooo........oo 634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 678 635