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: G4PenelopeAnnihilationModel.cc 74452 2013-10-07 15:08:00Z gcosmo $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 29 Oct 2008 L Pandola Migration from p 32 // 29 Oct 2008 L Pandola Migration from process to model 32 // 15 Apr 2009 V Ivanchenko Cleanup initiali 33 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of 33 // secondaries: 34 // secondaries: 34 // - apply internal high-ener 35 // - apply internal high-energy limit only in constructor 35 // - do not apply low-energy 36 // - do not apply low-energy limit (default is 0) 36 // - do not use G4ElementSele 37 // - do not use G4ElementSelector 37 // 02 Oct 2013 L.Pandola Migration to MT 38 // 02 Oct 2013 L.Pandola Migration to MT 38 39 39 #include "G4PenelopeAnnihilationModel.hh" 40 #include "G4PenelopeAnnihilationModel.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 42 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh" 43 #include "G4MaterialCutsCouple.hh" 44 #include "G4MaterialCutsCouple.hh" 44 #include "G4ProductionCutsTable.hh" 45 #include "G4ProductionCutsTable.hh" 45 #include "G4DynamicParticle.hh" 46 #include "G4DynamicParticle.hh" 46 #include "G4Gamma.hh" 47 #include "G4Gamma.hh" 47 48 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 49 50 50 G4double G4PenelopeAnnihilationModel::fPielr2 51 G4double G4PenelopeAnnihilationModel::fPielr2 = 0; 51 52 52 G4PenelopeAnnihilationModel::G4PenelopeAnnihil 53 G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition* part, 53 c 54 const G4String& nam) 54 :G4VEmModel(nam),fParticleChange(nullptr),fP << 55 :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false) 55 { 56 { 56 fIntrinsicLowEnergyLimit = 0.0; 57 fIntrinsicLowEnergyLimit = 0.0; 57 fIntrinsicHighEnergyLimit = 100.0*GeV; 58 fIntrinsicHighEnergyLimit = 100.0*GeV; >> 59 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 58 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 60 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 59 61 60 if (part) 62 if (part) 61 SetParticle(part); 63 SetParticle(part); 62 64 63 //Calculate variable that will be used later 65 //Calculate variable that will be used later on 64 fPielr2 = pi*classic_electr_radius*classic_e 66 fPielr2 = pi*classic_electr_radius*classic_electr_radius; 65 67 66 fVerboseLevel= 0; << 68 verboseLevel= 0; 67 // Verbosity scale: 69 // Verbosity scale: 68 // 0 = nothing 70 // 0 = nothing 69 // 1 = warning for energy non-conservation 71 // 1 = warning for energy non-conservation 70 // 2 = details of energy budget 72 // 2 = details of energy budget 71 // 3 = calculation of cross sections, file o 73 // 3 = calculation of cross sections, file openings, sampling of atoms 72 // 4 = entering in methods 74 // 4 = entering in methods >> 75 73 } 76 } 74 77 75 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 76 79 77 G4PenelopeAnnihilationModel::~G4PenelopeAnnihi 80 G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel() 78 {;} 81 {;} 79 82 80 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 84 82 void G4PenelopeAnnihilationModel::Initialise(c 85 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition* part, 83 const G4DataVector&) 86 const G4DataVector&) 84 { 87 { 85 if (fVerboseLevel > 3) << 88 if (verboseLevel > 3) 86 G4cout << "Calling G4PenelopeAnnihilationM 89 G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl; 87 SetParticle(part); 90 SetParticle(part); 88 91 89 if (IsMaster() && part == fParticle) 92 if (IsMaster() && part == fParticle) 90 { 93 { 91 94 92 if(fVerboseLevel > 0) { << 95 if(verboseLevel > 0) { 93 G4cout << "Penelope Annihilation model is in 96 G4cout << "Penelope Annihilation model is initialized " << G4endl 94 << "Energy range: " 97 << "Energy range: " 95 << LowEnergyLimit() / keV << " keV - 98 << LowEnergyLimit() / keV << " keV - " 96 << HighEnergyLimit() / GeV << " GeV" 99 << HighEnergyLimit() / GeV << " GeV" 97 << G4endl; 100 << G4endl; 98 } 101 } 99 } 102 } 100 103 101 if(fIsInitialised) return; << 104 if(isInitialised) return; 102 fParticleChange = GetParticleChangeForGamma( 105 fParticleChange = GetParticleChangeForGamma(); 103 fIsInitialised = true; << 106 isInitialised = true; 104 } 107 } 105 108 106 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 void G4PenelopeAnnihilationModel::InitialiseLo 110 void G4PenelopeAnnihilationModel::InitialiseLocal(const G4ParticleDefinition* part, 108 G4VEmModel* masterModel) 111 G4VEmModel* masterModel) 109 { 112 { 110 if (fVerboseLevel > 3) << 113 if (verboseLevel > 3) 111 G4cout << "Calling G4PenelopeAnnihilationM 114 G4cout << "Calling G4PenelopeAnnihilationModel::InitialiseLocal()" << G4endl; 112 115 113 // 116 // 114 //Check that particle matches: one might hav 117 //Check that particle matches: one might have multiple master models (e.g. 115 //for e+ and e-). 118 //for e+ and e-). 116 // 119 // 117 if (part == fParticle) 120 if (part == fParticle) 118 { 121 { 119 //Get the const table pointers from the 122 //Get the const table pointers from the master to the workers 120 const G4PenelopeAnnihilationModel* theMo 123 const G4PenelopeAnnihilationModel* theModel = 121 static_cast<G4PenelopeAnnihilationMode 124 static_cast<G4PenelopeAnnihilationModel*> (masterModel); 122 125 123 //Same verbosity for all workers, as the 126 //Same verbosity for all workers, as the master 124 fVerboseLevel = theModel->fVerboseLevel; << 127 verboseLevel = theModel->verboseLevel; 125 } 128 } 126 } 129 } 127 130 128 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 129 132 130 G4double G4PenelopeAnnihilationModel::ComputeC 133 G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom( 131 const G 134 const G4ParticleDefinition*, 132 G 135 G4double energy, 133 G 136 G4double Z, G4double, 134 G 137 G4double, G4double) 135 { 138 { 136 if (fVerboseLevel > 3) << 139 if (verboseLevel > 3) 137 G4cout << "Calling ComputeCrossSectionPerA 140 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << 138 G4endl; 141 G4endl; 139 142 140 G4double cs = Z*ComputeCrossSectionPerElectr 143 G4double cs = Z*ComputeCrossSectionPerElectron(energy); 141 144 142 if (fVerboseLevel > 2) << 145 if (verboseLevel > 2) 143 G4cout << "Annihilation cross Section at " 146 G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << 144 " = " << cs/barn << " barn" << G4endl; 147 " = " << cs/barn << " barn" << G4endl; 145 return cs; 148 return cs; 146 } 149 } 147 150 148 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 152 150 void G4PenelopeAnnihilationModel::SampleSecond 153 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 151 const G4MaterialCutsCouple*, 154 const G4MaterialCutsCouple*, 152 const G4DynamicParticle* aDyna 155 const G4DynamicParticle* aDynamicPositron, 153 G4double, 156 G4double, 154 G4double) 157 G4double) 155 { 158 { 156 // 159 // 157 // Penelope model to sample final state for 160 // Penelope model to sample final state for positron annihilation. 158 // Target eletrons are assumed to be free an 161 // Target eletrons are assumed to be free and at rest. Binding effects enabling 159 // one-photon annihilation are neglected. 162 // one-photon annihilation are neglected. 160 // For annihilation at rest, two back-to-bac 163 // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV 161 // and isotropic angular distribution. 164 // and isotropic angular distribution. 162 // For annihilation in flight, it is used th 165 // For annihilation in flight, it is used the theory from 163 // W. Heitler, The quantum theory of radiat 166 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954) 164 // The two photons can have different energy 167 // The two photons can have different energy. The efficiency of the sampling algorithm 165 // of the photon energy from the dSigma/dE d 168 // of the photon energy from the dSigma/dE distribution is practically 100% for 166 // positrons of kinetic energy < 10 keV. It 169 // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy 167 // of about 10 MeV. 170 // of about 10 MeV. 168 // The angle theta is kinematically linked t 171 // The angle theta is kinematically linked to the photon energy, to ensure momentum 169 // conservation. The angle phi is sampled is 172 // conservation. The angle phi is sampled isotropically for the first gamma. 170 // 173 // 171 if (fVerboseLevel > 3) << 174 if (verboseLevel > 3) 172 G4cout << "Calling SamplingSecondaries() o 175 G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl; 173 176 174 G4double kineticEnergy = aDynamicPositron->G 177 G4double kineticEnergy = aDynamicPositron->GetKineticEnergy(); 175 178 176 // kill primary 179 // kill primary 177 fParticleChange->SetProposedKineticEnergy(0. 180 fParticleChange->SetProposedKineticEnergy(0.); 178 fParticleChange->ProposeTrackStatus(fStopAnd 181 fParticleChange->ProposeTrackStatus(fStopAndKill); 179 182 180 if (kineticEnergy == 0.0) 183 if (kineticEnergy == 0.0) 181 { 184 { 182 //Old AtRestDoIt 185 //Old AtRestDoIt 183 G4double cosTheta = -1.0+2.0*G4UniformRa 186 G4double cosTheta = -1.0+2.0*G4UniformRand(); 184 G4double sinTheta = std::sqrt(1.0-cosThe 187 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta); 185 G4double phi = twopi*G4UniformRand(); 188 G4double phi = twopi*G4UniformRand(); 186 G4ThreeVector direction (sinTheta*std::c 189 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta); 187 G4DynamicParticle* firstGamma = new G4Dy 190 G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(), 188 direction, electron_mass_c2 191 direction, electron_mass_c2); 189 G4DynamicParticle* secondGamma = new G4D 192 G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(), 190 -direction, electron_mass_ 193 -direction, electron_mass_c2); 191 194 192 fvect->push_back(firstGamma); 195 fvect->push_back(firstGamma); 193 fvect->push_back(secondGamma); 196 fvect->push_back(secondGamma); 194 return; 197 return; 195 } 198 } 196 199 197 //This is the "PostStep" case (annihilation 200 //This is the "PostStep" case (annihilation in flight) 198 G4ParticleMomentum positronDirection = 201 G4ParticleMomentum positronDirection = 199 aDynamicPositron->GetMomentumDirection(); 202 aDynamicPositron->GetMomentumDirection(); 200 G4double gamma = 1.0 + std::max(kineticEnerg 203 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2; 201 G4double gamma21 = std::sqrt(gamma*gamma-1); 204 G4double gamma21 = std::sqrt(gamma*gamma-1); 202 G4double ani = 1.0+gamma; 205 G4double ani = 1.0+gamma; 203 G4double chimin = 1.0/(ani+gamma21); 206 G4double chimin = 1.0/(ani+gamma21); 204 G4double rchi = (1.0-chimin)/chimin; 207 G4double rchi = (1.0-chimin)/chimin; 205 G4double gt0 = ani*ani-2.0; 208 G4double gt0 = ani*ani-2.0; 206 G4double test=0.0; 209 G4double test=0.0; 207 G4double epsilon = 0; 210 G4double epsilon = 0; 208 do{ 211 do{ 209 epsilon = chimin*std::pow(rchi,G4UniformRa 212 epsilon = chimin*std::pow(rchi,G4UniformRand()); 210 G4double reject = ani*ani*(1.0-epsilon)+2. 213 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon); 211 test = G4UniformRand()*gt0-reject; 214 test = G4UniformRand()*gt0-reject; 212 }while(test>0); 215 }while(test>0); 213 216 214 G4double totalAvailableEnergy = kineticEnerg 217 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2; 215 G4double photon1Energy = epsilon*totalAvaila 218 G4double photon1Energy = epsilon*totalAvailableEnergy; 216 G4double photon2Energy = (1.0-epsilon)*total 219 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy; 217 G4double cosTheta1 = (ani-1.0/epsilon)/gamma 220 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21; 218 G4double cosTheta2 = (ani-1.0/(1.0-epsilon)) 221 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21; 219 222 >> 223 //G4double localEnergyDeposit = 0.; >> 224 220 G4double sinTheta1 = std::sqrt(1.-cosTheta1* 225 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1); 221 G4double phi1 = twopi * G4UniformRand(); 226 G4double phi1 = twopi * G4UniformRand(); 222 G4double dirx1 = sinTheta1 * std::cos(phi1); 227 G4double dirx1 = sinTheta1 * std::cos(phi1); 223 G4double diry1 = sinTheta1 * std::sin(phi1); 228 G4double diry1 = sinTheta1 * std::sin(phi1); 224 G4double dirz1 = cosTheta1; 229 G4double dirz1 = cosTheta1; 225 230 226 G4double sinTheta2 = std::sqrt(1.-cosTheta2* 231 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2); 227 G4double phi2 = phi1+pi; 232 G4double phi2 = phi1+pi; 228 G4double dirx2 = sinTheta2 * std::cos(phi2); 233 G4double dirx2 = sinTheta2 * std::cos(phi2); 229 G4double diry2 = sinTheta2 * std::sin(phi2); 234 G4double diry2 = sinTheta2 * std::sin(phi2); 230 G4double dirz2 = cosTheta2; 235 G4double dirz2 = cosTheta2; 231 236 232 G4ThreeVector photon1Direction (dirx1,diry1, 237 G4ThreeVector photon1Direction (dirx1,diry1,dirz1); 233 photon1Direction.rotateUz(positronDirection) 238 photon1Direction.rotateUz(positronDirection); 234 // create G4DynamicParticle object for the p 239 // create G4DynamicParticle object for the particle1 235 G4DynamicParticle* aParticle1= new G4Dynamic 240 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(), 236 photon1Direction, 241 photon1Direction, 237 photon1Energy); 242 photon1Energy); 238 fvect->push_back(aParticle1); 243 fvect->push_back(aParticle1); 239 244 240 G4ThreeVector photon2Direction(dirx2,diry2,d 245 G4ThreeVector photon2Direction(dirx2,diry2,dirz2); 241 photon2Direction.rotateUz(positronDirection) 246 photon2Direction.rotateUz(positronDirection); 242 // create G4DynamicParticle object for the p << 247 // create G4DynamicParticle object for the particle2 243 G4DynamicParticle* aParticle2= new G4Dynamic 248 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(), 244 photon2Direction, 249 photon2Direction, 245 photon2Energy); 250 photon2Energy); 246 fvect->push_back(aParticle2); 251 fvect->push_back(aParticle2); 247 252 248 if (fVerboseLevel > 1) << 253 if (verboseLevel > 1) 249 { 254 { 250 G4cout << "----------------------------- 255 G4cout << "-----------------------------------------------------------" << G4endl; 251 G4cout << "Energy balance from G4Penelop 256 G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl; 252 G4cout << "Kinetic positron energy: " << 257 G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl; 253 G4cout << "Total available energy: " << 258 G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl; 254 G4cout << "----------------------------- 259 G4cout << "-----------------------------------------------------------" << G4endl; 255 G4cout << "Photon energy 1: " << photon1 260 G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl; 256 G4cout << "Photon energy 2: " << photon2 261 G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl; 257 G4cout << "Total final state: " << (phot 262 G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << 258 " keV" << G4endl; 263 " keV" << G4endl; 259 G4cout << "----------------------------- 264 G4cout << "-----------------------------------------------------------" << G4endl; 260 } 265 } 261 if (fVerboseLevel > 0) << 266 if (verboseLevel > 0) 262 { 267 { 263 G4double energyDiff = std::fabs(totalAva 268 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy); 264 if (energyDiff > 0.05*keV) 269 if (energyDiff > 0.05*keV) 265 G4cout << "Warning from G4PenelopeAnnihilati 270 G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << 266 (photon1Energy+photon2Energy)/keV << 271 (photon1Energy+photon2Energy)/keV << 267 " keV (final) vs. " << 272 " keV (final) vs. " << 268 totalAvailableEnergy/keV << " keV (initial 273 totalAvailableEnergy/keV << " keV (initial)" << G4endl; 269 } 274 } 270 return; 275 return; 271 } 276 } 272 277 273 //....oooOO0OOooo........oooOO0OOooo........oo 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 274 279 275 G4double G4PenelopeAnnihilationModel:: Compute 280 G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy) 276 { 281 { 277 // 282 // 278 // Penelope model to calculate cross section 283 // Penelope model to calculate cross section for positron annihilation. 279 // The annihilation cross section per electr 284 // The annihilation cross section per electron is calculated according 280 // to the Heitler formula 285 // to the Heitler formula 281 // W. Heitler, The quantum theory of radiat 286 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954) 282 // in the assumptions of electrons free and 287 // in the assumptions of electrons free and at rest. 283 // 288 // 284 G4double gamma = 1.0+std::max(energy,1.0*eV) 289 G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2; 285 G4double gamma2 = gamma*gamma; 290 G4double gamma2 = gamma*gamma; 286 G4double f2 = gamma2-1.0; 291 G4double f2 = gamma2-1.0; 287 G4double f1 = std::sqrt(f2); 292 G4double f1 = std::sqrt(f2); 288 G4double crossSection = fPielr2*((gamma2+4.0 << 293 G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2 289 - (gamma+3.0)/f1)/(gamma+1.0); 294 - (gamma+3.0)/f1)/(gamma+1.0); 290 return crossSection; 295 return crossSection; 291 } 296 } 292 297 293 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 294 299 295 void G4PenelopeAnnihilationModel::SetParticle( 300 void G4PenelopeAnnihilationModel::SetParticle(const G4ParticleDefinition* p) 296 { 301 { 297 if(!fParticle) { 302 if(!fParticle) { 298 fParticle = p; 303 fParticle = p; 299 } 304 } 300 } 305 } 301 306