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 // $Id: G4eeToTwoGammaModel.cc,v 1.14 2007/05/23 08:47:35 vnivanch Exp $ >> 27 // GEANT4 tag $Name: geant4-09-00-patch-01 $ 26 // 28 // 27 // ------------------------------------------- 29 // ------------------------------------------------------------------- 28 // 30 // 29 // GEANT4 Class file 31 // GEANT4 Class file 30 // 32 // 31 // 33 // 32 // File name: G4eeToTwoGammaModel 34 // File name: G4eeToTwoGammaModel 33 // 35 // 34 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of Michel Maire code 35 // 37 // 36 // Creation date: 02.08.2004 38 // Creation date: 02.08.2004 37 // 39 // 38 // Modifications: 40 // Modifications: 39 // 08-04-05 Major optimisation of internal int 41 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko) 40 // 18-04-05 Compute CrossSectionPerVolume (V.I 42 // 18-04-05 Compute CrossSectionPerVolume (V.Ivanchenko) 41 // 06-02-06 ComputeCrossSectionPerElectron, Co 43 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma) 42 // 29-06-06 Fix problem for zero energy incide 44 // 29-06-06 Fix problem for zero energy incident positron (V.Ivanchenko) 43 // 20-10-06 Add theGamma as a member (V.Ivanch 45 // 20-10-06 Add theGamma as a member (V.Ivanchenko) 44 // 18-01-20 Introduce thermal model of annihil << 45 // 46 // 46 // 47 // 47 // Class Description: 48 // Class Description: 48 // 49 // 49 // Implementation of e+ annihilation into 2 ga 50 // Implementation of e+ annihilation into 2 gamma 50 // 51 // 51 // The secondaries Gamma energies are sampled 52 // The secondaries Gamma energies are sampled using the Heitler cross section. 52 // 53 // 53 // A modified version of the random number tec 54 // A modified version of the random number techniques of Butcher & Messel 54 // is used (Nuc Phys 20(1960),15). 55 // is used (Nuc Phys 20(1960),15). 55 // 56 // 56 // GEANT4 internal units. 57 // GEANT4 internal units. 57 // 58 // 58 // Note 1: The initial electron is assumed fre << 59 // Note 1: The initial electron is assumed free and at rest. 59 // is not defined << 60 // 60 // 61 // Note 2: The annihilation processes producin 61 // Note 2: The annihilation processes producing one or more than two photons are 62 // ignored, as negligible compared to 62 // ignored, as negligible compared to the two photons process. 63 63 >> 64 >> 65 64 // 66 // 65 // ------------------------------------------- 67 // ------------------------------------------------------------------- 66 // 68 // 67 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 68 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 71 70 #include "G4eeToTwoGammaModel.hh" 72 #include "G4eeToTwoGammaModel.hh" 71 #include "G4PhysicalConstants.hh" << 72 #include "G4SystemOfUnits.hh" << 73 #include "G4TrackStatus.hh" 73 #include "G4TrackStatus.hh" 74 #include "G4Electron.hh" 74 #include "G4Electron.hh" 75 #include "G4Positron.hh" 75 #include "G4Positron.hh" 76 #include "G4Gamma.hh" 76 #include "G4Gamma.hh" 77 #include "Randomize.hh" 77 #include "Randomize.hh" 78 #include "G4RandomDirection.hh" << 79 #include "G4ParticleChangeForGamma.hh" 78 #include "G4ParticleChangeForGamma.hh" 80 #include "G4EmParameters.hh" << 81 #include "G4Log.hh" << 82 #include "G4Exp.hh" << 83 79 84 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 81 >> 82 using namespace std; >> 83 86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const 84 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*, 87 const 85 const G4String& nam) 88 : G4VEmModel(nam), 86 : G4VEmModel(nam), 89 pi_rcl2(CLHEP::pi*CLHEP::classic_electr_ra << 87 pi_rcl2(pi*classic_electr_radius*classic_electr_radius), >> 88 isInitialised(false) 90 { 89 { 91 theGamma = G4Gamma::Gamma(); 90 theGamma = G4Gamma::Gamma(); 92 fParticleChange = nullptr; << 93 } 91 } 94 92 95 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 96 94 97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = << 95 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() >> 96 {} 98 97 99 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 99 101 void G4eeToTwoGammaModel::Initialise(const G4P 100 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*, 102 const G4D 101 const G4DataVector&) 103 { 102 { 104 if (nullptr != fParticleChange) { return; } << 103 if(isInitialised) return; 105 fParticleChange = GetParticleChangeForGamma( << 104 >> 105 if(pParticleChange) >> 106 fParticleChange = >> 107 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); >> 108 else >> 109 fParticleChange = new G4ParticleChangeForGamma(); >> 110 >> 111 isInitialised = true; 106 } 112 } 107 113 108 //....oooOO0OOooo........oooOO0OOooo........oo 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 109 115 110 G4double << 116 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerElectron( 111 G4eeToTwoGammaModel::ComputeCrossSectionPerEle << 117 const G4ParticleDefinition*, >> 118 G4double kineticEnergy, >> 119 G4double, G4double) 112 { 120 { 113 // Calculates the cross section per electron 121 // Calculates the cross section per electron of annihilation into two photons 114 // from the Heilter formula. 122 // from the Heilter formula. 115 << 123 116 G4double ekin = std::max(CLHEP::eV, kinetic << 124 G4double tau = kineticEnergy/electron_mass_c2; 117 << 118 G4double tau = ekin/CLHEP::electron_mass_c << 119 G4double gam = tau + 1.0; 125 G4double gam = tau + 1.0; 120 G4double gamma2= gam*gam; 126 G4double gamma2= gam*gam; 121 G4double bg2 = tau * (tau+2.0); 127 G4double bg2 = tau * (tau+2.0); 122 G4double bg = std::sqrt(bg2); << 128 G4double bg = sqrt(bg2); 123 129 124 G4double cross = pi_rcl2*((gamma2+4*gam+1.)* << 130 G4double cross = pi_rcl2*((gamma2+4*gam+1.)*log(gam+bg) - (gam+3.)*bg) 125 / (bg2*(gam+1.)); 131 / (bg2*(gam+1.)); 126 return cross; 132 return cross; 127 } 133 } 128 134 129 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 130 136 131 G4double G4eeToTwoGammaModel::ComputeCrossSect 137 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom( 132 const G4Pa << 138 const G4ParticleDefinition* p, 133 G4double k 139 G4double kineticEnergy, G4double Z, 134 G4double, G4double, G4double) 140 G4double, G4double, G4double) 135 { 141 { 136 // Calculates the cross section per atom of 142 // Calculates the cross section per atom of annihilation into two photons 137 return Z*ComputeCrossSectionPerElectron(kine << 143 >> 144 G4double cross = Z*ComputeCrossSectionPerElectron(p,kineticEnergy); >> 145 return cross; 138 } 146 } 139 147 140 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 141 149 142 G4double G4eeToTwoGammaModel::CrossSectionPerV 150 G4double G4eeToTwoGammaModel::CrossSectionPerVolume( 143 const G4Material* material, 151 const G4Material* material, 144 const G4ParticleDefinition*, << 152 const G4ParticleDefinition* p, 145 G4double kineticEnergy, 153 G4double kineticEnergy, 146 G4double, G4double) 154 G4double, G4double) 147 { 155 { 148 // Calculates the cross section per volume o 156 // Calculates the cross section per volume of annihilation into two photons 149 return material->GetElectronDensity()*Comput << 157 >> 158 G4double eDensity = material->GetElectronDensity(); >> 159 G4double cross = eDensity*ComputeCrossSectionPerElectron(p,kineticEnergy); >> 160 return cross; 150 } 161 } 151 162 152 //....oooOO0OOooo........oooOO0OOooo........oo 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 153 164 154 // Polarisation of gamma according to M.H.L.Pr << 165 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp, 155 // Nature 4065 (1947) 435. << 156 << 157 void G4eeToTwoGammaModel::SampleSecondaries(st << 158 const G4MaterialCutsCouple*, 166 const G4MaterialCutsCouple*, 159 const G4DynamicParticle* dp, 167 const G4DynamicParticle* dp, 160 G4double, 168 G4double, 161 G4double) 169 G4double) 162 { 170 { 163 // kill primary positron << 171 G4double PositKinEnergy = dp->GetKineticEnergy(); 164 fParticleChange->SetProposedKineticEnergy(0. << 165 fParticleChange->ProposeTrackStatus(fStopAnd << 166 172 167 // Case at rest not considered anymore insid << 173 // Case at rest 168 G4LorentzVector lv(dp->GetMomentum(), << 174 if(PositKinEnergy == 0.0) { 169 dp->GetKineticEnergy() + 2*CLHEP::ele << 175 G4double cost = 2.*G4UniformRand()-1.; 170 G4double eGammaCMS = 0.5 * lv.mag(); << 176 G4double sint = sqrt((1. - cost)*(1. + cost)); 171 << 177 G4double phi = twopi * G4UniformRand(); 172 G4ThreeVector dir1 = G4RandomDirection(); << 178 G4ThreeVector dir (sint*cos(phi), sint*sin(phi), cost); 173 G4double phi = CLHEP::twopi * G4UniformRand( << 179 G4DynamicParticle* aGamma1 = new G4DynamicParticle(theGamma, 174 G4double cosphi = std::cos(phi); << 180 dir, electron_mass_c2); 175 G4double sinphi = std::sin(phi); << 181 G4DynamicParticle* aGamma2 = new G4DynamicParticle(theGamma, 176 G4ThreeVector pol1(cosphi, sinphi, 0.0); << 182 -dir, electron_mass_c2); 177 pol1.rotateUz(dir1); << 183 vdp->push_back(aGamma1); 178 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCM << 184 vdp->push_back(aGamma2); 179 << 185 180 G4ThreeVector pol2(-sinphi, cosphi, 0.0); << 186 } else { 181 pol2.rotateUz(dir1); << 187 182 << 188 G4ThreeVector PositDirection = dp->GetMomentumDirection(); 183 // transformation to lab system << 189 184 lv1.boost(lv.boostVector()); << 190 G4double tau = PositKinEnergy/electron_mass_c2; 185 lv -= lv1; << 191 G4double gam = tau + 1.0; 186 << 192 G4double tau2 = tau + 2.0; 187 //!!! boost of polarisation vector is not ye << 193 G4double sqgrate = sqrt(tau/tau2)*0.5; 188 << 194 G4double sqg2m1 = sqrt(tau*tau2); 189 // use constructors optimal for massless par << 195 190 auto aGamma1 = new G4DynamicParticle(G4Gamma << 196 // limits of the energy sampling 191 aGamma1->SetPolarization(pol1); << 197 G4double epsilmin = 0.5 - sqgrate; 192 auto aGamma2 = new G4DynamicParticle(G4Gamma << 198 G4double epsilmax = 0.5 + sqgrate; 193 aGamma2->SetPolarization(pol2); << 199 G4double epsilqot = epsilmax/epsilmin; 194 << 200 195 vdp->push_back(aGamma1); << 201 // 196 vdp->push_back(aGamma2); << 202 // sample the energy rate of the created gammas >> 203 // >> 204 G4double epsil, greject; >> 205 >> 206 do { >> 207 epsil = epsilmin*pow(epsilqot,G4UniformRand()); >> 208 greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2); >> 209 } while( greject < G4UniformRand() ); >> 210 >> 211 // >> 212 // scattered Gamma angles. ( Z - axis along the parent positron) >> 213 // >> 214 >> 215 G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1); >> 216 if(std::abs(cost) > 1.0) { >> 217 G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost >> 218 << " positron Ekin(MeV)= " << PositKinEnergy >> 219 << " gamma epsil= " << epsil >> 220 << G4endl; >> 221 if(cost > 1.0) cost = 1.0; >> 222 else cost = -1.0; >> 223 } >> 224 G4double sint = sqrt((1.+cost)*(1.-cost)); >> 225 G4double phi = twopi * G4UniformRand(); >> 226 >> 227 G4double dirx = sint*cos(phi) , diry = sint*sin(phi) , dirz = cost; >> 228 >> 229 // >> 230 // kinematic of the created pair >> 231 // >> 232 >> 233 G4double TotalAvailableEnergy = PositKinEnergy + 2.0*electron_mass_c2; >> 234 G4double Phot1Energy = epsil*TotalAvailableEnergy; >> 235 >> 236 G4ThreeVector Phot1Direction (dirx, diry, dirz); >> 237 Phot1Direction.rotateUz(PositDirection); >> 238 G4DynamicParticle* aGamma1 = >> 239 new G4DynamicParticle (theGamma,Phot1Direction, Phot1Energy); >> 240 vdp->push_back(aGamma1); >> 241 >> 242 G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy; >> 243 G4double PositP= sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2)); >> 244 G4ThreeVector dir = PositDirection*PositP - Phot1Direction*Phot1Energy; >> 245 G4ThreeVector Phot2Direction = dir.unit(); >> 246 >> 247 // create G4DynamicParticle object for the particle2 >> 248 G4DynamicParticle* aGamma2= >> 249 new G4DynamicParticle (theGamma,Phot2Direction, Phot2Energy); >> 250 vdp->push_back(aGamma2); >> 251 /* >> 252 G4cout << "Annihilation in fly: e0= " << PositKinEnergy >> 253 << " m= " << electron_mass_c2 >> 254 << " e1= " << Phot1Energy >> 255 << " e2= " << Phot2Energy << " dir= " << dir >> 256 << " -> " << Phot1Direction << " " >> 257 << Phot2Direction << G4endl; >> 258 */ >> 259 } >> 260 fParticleChange->SetProposedKineticEnergy(0.); >> 261 fParticleChange->ProposeTrackStatus(fStopAndKill); 197 } 262 } 198 263 199 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 200 265