Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // ------------------------------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // 32 // File name: G4eplusTo3GammaOKVIModel 33 // 34 // Authors: Andrei Alkin, Vladimir Ivanchenko, Omrame Kadri 35 // 36 // Creation date: 29.03.2018 37 // 38 // ------------------------------------------------------------------- 39 // 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 42 43 #include "G4eplusTo3GammaOKVIModel.hh" 44 #include "G4PhysicalConstants.hh" 45 #include "G4SystemOfUnits.hh" 46 #include "G4EmParameters.hh" 47 #include "G4TrackStatus.hh" 48 #include "G4Electron.hh" 49 #include "G4Positron.hh" 50 #include "G4Gamma.hh" 51 #include "Randomize.hh" 52 #include "G4ParticleChangeForGamma.hh" 53 #include "G4Log.hh" 54 #include "G4Exp.hh" 55 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 57 58 using namespace std; 59 60 G4eplusTo3GammaOKVIModel::G4eplusTo3GammaOKVIModel(const G4ParticleDefinition*, 61 const G4String& nam) 62 : G4VEmModel(nam), fDelta(0.001) 63 { 64 theGamma = G4Gamma::Gamma(); 65 } 66 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 68 69 G4eplusTo3GammaOKVIModel::~G4eplusTo3GammaOKVIModel() = default; 70 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 72 73 void G4eplusTo3GammaOKVIModel::Initialise(const G4ParticleDefinition*, 74 const G4DataVector&) 75 {} 76 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 78 79 // (A.A.) F_{ijk} calculation method 80 G4double G4eplusTo3GammaOKVIModel::ComputeF(G4double fr1, G4double fr2, 81 G4double fr3, G4double kinEnergy) 82 { 83 G4double ekin = std::max(eV,kinEnergy); 84 G4double tau = ekin/electron_mass_c2; 85 G4double gam = tau + 1.0; 86 G4double gamma2 = gam*gam; 87 G4double bg2 = tau * (tau+2.0); 88 G4double bg = sqrt(bg2); 89 90 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.) 91 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.; 92 G4double border; 93 94 if(ekin < 500*MeV) { 95 border = 1. - (electron_mass_c2)/(2*(ekin + electron_mass_c2)); 96 } else { 97 border = 1. - (100*electron_mass_c2)/(2*(ekin + electron_mass_c2)); 98 } 99 100 border = std::min(border, 0.9999); 101 102 if (fr1>border) { fr1 = border; } 103 if (fr2>border) { fr2 = border; } 104 if (fr3>border) { fr3 = border; } 105 106 G4double fr1s = fr1*fr1; // "s" for "squared" 107 G4double fr2s = fr2*fr2; 108 G4double fr3s = fr3*fr3; 109 110 G4double aa = (1.-fr1)*(1.-fr2); 111 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2); 112 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa); 113 114 G4double fres = -rho*(1./fr1s + 1./fr2s) 115 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2))) 116 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add; 117 118 return fres; 119 } 120 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 122 123 // (A.A.) F_{ijk} calculation method 124 G4double G4eplusTo3GammaOKVIModel::ComputeF0(G4double fr1, G4double fr2, 125 G4double fr3) 126 { 127 G4double tau = 0.0; 128 G4double gam = tau + 1.0; 129 G4double gamma2 = gam*gam; 130 G4double bg2 = tau * (tau+2.0); 131 G4double bg = sqrt(bg2); 132 133 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.) 134 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.; 135 G4double border = 0.5; 136 137 if (fr1>border) { fr1 = border; } 138 if (fr2>border) { fr2 = border; } 139 if (fr3>border) { fr3 = border; } 140 141 G4double fr1s = fr1*fr1; // "s" for "squared" 142 G4double fr2s = fr2*fr2; 143 G4double fr3s = fr3*fr3; 144 145 G4double aa = (1.-fr1)*(1.-fr2); 146 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2); 147 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa); 148 149 G4double fres = -rho*(1./fr1s + 1./fr2s) 150 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2))) 151 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add; 152 153 return fres; 154 } 155 156 //(A.A.) diff x-sections for maximum search and rejection 157 G4double G4eplusTo3GammaOKVIModel::ComputeFS(G4double fr1, 158 G4double fr2, G4double fr3, G4double kinEnergy) 159 { 160 G4double ekin = std::max(eV,kinEnergy); 161 G4double tau = ekin/electron_mass_c2; 162 G4double gam = tau + 1.0; 163 164 G4double fsum = fr1*fr1*(ComputeF(fr1,fr2,fr3,ekin) + 165 ComputeF(fr3,fr1,fr2,ekin) + 166 ComputeF(fr2,fr3,fr1,ekin)); 167 168 G4double dcross = fsum/((3*fr1*fr1*(gam+1.))); 169 170 return dcross; 171 } 172 173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 174 175 G4double 176 G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerElectron(G4double kinEnergy) 177 { 178 // Calculates the cross section per electron of annihilation into 3 photons 179 // from the Heilter formula. 180 181 G4double ekin = std::max(CLHEP::eV, kinEnergy); 182 G4double tau = ekin/CLHEP::electron_mass_c2; 183 G4double gam = tau + 1.0; 184 G4double gamma2 = gam*gam; 185 G4double bg2 = tau * (tau+2.0); 186 G4double bg = sqrt(bg2); 187 188 G4double rho = (gamma2+4*gam+1.)*G4Log(gam+bg)/(gamma2-1.) 189 - (gam+3.)/(sqrt(gam*gam - 1.)); 190 191 G4double cross = alpha_rcl2*(4.2 - (2.*G4Log(fDelta)+1.)*rho*rho)/(gam+1.); 192 return cross; 193 } 194 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 196 197 G4double G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerAtom( 198 const G4ParticleDefinition*, 199 G4double kineticEnergy, G4double Z, 200 G4double, G4double, G4double) 201 { 202 G4double cross = Z*ComputeCrossSectionPerElectron(kineticEnergy); 203 return cross; 204 } 205 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 207 208 G4double G4eplusTo3GammaOKVIModel::CrossSectionPerVolume( 209 const G4Material* material, 210 const G4ParticleDefinition*, 211 G4double kineticEnergy, 212 G4double, G4double) 213 { 214 // Calculates the cross section per volume of annihilation into two photons 215 216 G4double eDensity = material->GetElectronDensity(); 217 G4double cross = eDensity*ComputeCrossSectionPerElectron(kineticEnergy); 218 return cross; 219 } 220 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 222 223 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward, 224 // Nature 4065 (1947) 435. 225 226 void 227 G4eplusTo3GammaOKVIModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp, 228 const G4MaterialCutsCouple*, 229 const G4DynamicParticle* dp, 230 G4double, G4double) 231 { 232 // let us perform sampling in C.M.S. reference frame of e- at rest and e+ on fly 233 G4double posiKinEnergy = dp->GetKineticEnergy(); 234 G4LorentzVector lv(dp->GetMomentum(), 235 posiKinEnergy + 2*CLHEP::electron_mass_c2); 236 G4double eGammaCMS = 0.5 * lv.mag(); 237 238 // the limit value fDelta is defined by a class, which call this method 239 // thickness of border defined by C.M.S. energy 240 G4double border = 241 1.0 - std::min(std::max(CLHEP::electron_mass_c2/eGammaCMS, fDelta), 0.1); 242 243 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 244 245 G4ThreeVector posiDirection = dp->GetMomentumDirection(); 246 247 // (A.A.) LIMITS FOR 1st GAMMA 248 G4double xmin = 0.01; 249 G4double xmax = 0.667; // CHANGE to 3/2 250 251 G4double d1, d0, x1, x2, dmax, x2min; 252 253 // (A.A.) sampling of x1 x2 x3 (whole cycle of rejection) 254 do { 255 x1 = 1./((1./xmin) - ((1./xmin)-(1./xmax))*rndmEngine->flat()); 256 dmax = ComputeFS(eGammaCMS, x1, 1.-x1, border); 257 x2min = 1. - x1; 258 x2 = 1 - rndmEngine->flat()*(1. - x2min); 259 d1 = dmax*rndmEngine->flat(); 260 d0 = ComputeFS(eGammaCMS, x1, x2, 2.-x1-x2); 261 } 262 while(d0 < d1); 263 264 G4double x3 = 2 - x1 - x2; 265 // 266 // angles between Gammas 267 // 268 G4double psi13 = 2*std::asin(std::sqrt(std::abs((x1+x3-1.)/(x1*x3)))); 269 G4double psi12 = 2*std::asin(std::sqrt(std::abs((x1+x2-1.)/(x1*x2)))); 270 // 271 // kinematic of the created pair 272 // 273 G4double phot1Energy = x1*eGammaCMS; 274 G4double phot2Energy = x2*eGammaCMS; 275 G4double phot3Energy = x3*eGammaCMS; 276 277 // DIRECTIONS 278 279 // The azimuthal angles of q1 and q3 with respect to some plane 280 // through the beam axis are generated at random. 281 282 G4ThreeVector phot1Direction(0, 0, 1); 283 G4ThreeVector phot2Direction(0, std::sin(psi12), std::cos(psi12)); 284 G4ThreeVector phot3Direction(0, std::sin(psi13), std::cos(psi13)); 285 286 G4LorentzVector lv1(phot1Energy*phot1Direction, phot1Energy); 287 G4LorentzVector lv2(phot2Energy*phot2Direction, phot2Energy); 288 G4LorentzVector lv3(phot3Energy*phot3Direction, phot3Energy); 289 290 auto boostV = lv.boostVector(); 291 lv1.boost(boostV); 292 lv2.boost(boostV); 293 lv3.boost(boostV); 294 295 auto aGamma1 = new G4DynamicParticle (theGamma, lv1.vect()); 296 auto aGamma2 = new G4DynamicParticle (theGamma, lv2.vect()); 297 auto aGamma3 = new G4DynamicParticle (theGamma, lv3.vect()); 298 299 //!!! POLARIZATION - not yet implemented 300 301 vdp->push_back(aGamma1); 302 vdp->push_back(aGamma2); 303 vdp->push_back(aGamma3); 304 } 305 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 307