Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // 32 // File name: G4eplusTo3GammaOKVIModel 33 // 34 // Authors: Andrei Alkin, Vladimir Ivanchenko 35 // 36 // Creation date: 29.03.2018 37 // 38 // ------------------------------------------- 39 // 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oo 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........oo 57 58 using namespace std; 59 60 G4eplusTo3GammaOKVIModel::G4eplusTo3GammaOKVIM 61 62 : G4VEmModel(nam), fDelta(0.001) 63 { 64 theGamma = G4Gamma::Gamma(); 65 } 66 67 //....oooOO0OOooo........oooOO0OOooo........oo 68 69 G4eplusTo3GammaOKVIModel::~G4eplusTo3GammaOKVI 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 72 73 void G4eplusTo3GammaOKVIModel::Initialise(cons 74 const G4DataVector&) 75 {} 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 78 79 // (A.A.) F_{ijk} calculation method 80 G4double G4eplusTo3GammaOKVIModel::ComputeF(G4 81 G4 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+ 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 96 } else { 97 border = 1. - (100*electron_mass_c2)/(2*(e 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 "square 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) 113 114 G4double fres = -rho*(1./fr1s + 1./fr2s) 115 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/ 116 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*( 117 118 return fres; 119 } 120 121 //....oooOO0OOooo........oooOO0OOooo........oo 122 123 // (A.A.) F_{ijk} calculation method 124 G4double G4eplusTo3GammaOKVIModel::ComputeF0(G 125 G 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+ 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 "square 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) 148 149 G4double fres = -rho*(1./fr1s + 1./fr2s) 150 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/ 151 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*( 152 153 return fres; 154 } 155 156 //(A.A.) diff x-sections for maximum search an 157 G4double G4eplusTo3GammaOKVIModel::ComputeFS(G 158 G4double fr2, G4double fr3, G4double 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,fr 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........oo 174 175 G4double 176 G4eplusTo3GammaOKVIModel::ComputeCrossSectionP 177 { 178 // Calculates the cross section per electron 179 // from the Heilter formula. 180 181 G4double ekin = std::max(CLHEP::eV, kinEne 182 G4double tau = ekin/CLHEP::electron_mass_ 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+b 189 - (gam+3.)/(sqrt(gam*gam - 1.)); 190 191 G4double cross = alpha_rcl2*(4.2 - (2.*G4Log 192 return cross; 193 } 194 195 //....oooOO0OOooo........oooOO0OOooo........oo 196 197 G4double G4eplusTo3GammaOKVIModel::ComputeCros 198 const G4Pa 199 G4double k 200 G4double, G4double, G4double) 201 { 202 G4double cross = Z*ComputeCrossSectionPerEle 203 return cross; 204 } 205 206 //....oooOO0OOooo........oooOO0OOooo........oo 207 208 G4double G4eplusTo3GammaOKVIModel::CrossSectio 209 const G4Material* material, 210 const G4ParticleDefinition*, 211 G4double kineticEnergy, 212 G4double, G4double) 213 { 214 // Calculates the cross section per volume o 215 216 G4double eDensity = material->GetElectronDen 217 G4double cross = eDensity*ComputeCrossSectio 218 return cross; 219 } 220 221 //....oooOO0OOooo........oooOO0OOooo........oo 222 223 // Polarisation of gamma according to M.H.L.Pr 224 // Nature 4065 (1947) 435. 225 226 void 227 G4eplusTo3GammaOKVIModel::SampleSecondaries(ve 228 const G4MaterialCutsCouple*, 229 const G4DynamicParticle* dp, 230 G4double, G4double) 231 { 232 // let us perform sampling in C.M.S. referen 233 G4double posiKinEnergy = dp->GetKineticEnerg 234 G4LorentzVector lv(dp->GetMomentum(), 235 posiKinEnergy + 2*CLHEP::electron_mas 236 G4double eGammaCMS = 0.5 * lv.mag(); 237 238 // the limit value fDelta is defined by a cl 239 // thickness of border defined by C.M.S. ene 240 G4double border = 241 1.0 - std::min(std::max(CLHEP::electron_ma 242 243 CLHEP::HepRandomEngine* rndmEngine = G4Rando 244 245 G4ThreeVector posiDirection = dp->GetMomentu 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 254 do { 255 x1 = 1./((1./xmin) - ((1./xmin)-(1./xmax)) 256 dmax = ComputeFS(eGammaCMS, x1, 1.-x1, bor 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:: 269 G4double psi12 = 2*std::asin(std::sqrt(std:: 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 re 280 // through the beam axis are generated at ra 281 282 G4ThreeVector phot1Direction(0, 0, 1); 283 G4ThreeVector phot2Direction(0, std::sin(psi 284 G4ThreeVector phot3Direction(0, std::sin(psi 285 286 G4LorentzVector lv1(phot1Energy*phot1Directi 287 G4LorentzVector lv2(phot2Energy*phot2Directi 288 G4LorentzVector lv3(phot3Energy*phot3Directi 289 290 auto boostV = lv.boostVector(); 291 lv1.boost(boostV); 292 lv2.boost(boostV); 293 lv3.boost(boostV); 294 295 auto aGamma1 = new G4DynamicParticle (theGam 296 auto aGamma2 = new G4DynamicParticle (theGam 297 auto aGamma3 = new G4DynamicParticle (theGam 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........oo 307