// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: G4PenelopeAnnihilationModel.cc,v 1.3 2009/04/17 10:29:20 vnivanch Exp $ // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ // // Author: Luciano Pandola // // History: // -------- // 29 Oct 2008 L Pandola Migration from process to model // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries: // - apply internal high-energy limit only in constructor // - do not apply low-energy limit (default is 0) // - do not use G4ElementSelector #include "G4PenelopeAnnihilationModel.hh" #include "G4ParticleDefinition.hh" #include "G4MaterialCutsCouple.hh" #include "G4ProductionCutsTable.hh" #include "G4DynamicParticle.hh" #include "G4Gamma.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*, const G4String& nam) :G4VEmModel(nam),isInitialised(false) { fIntrinsicLowEnergyLimit = 0.0; fIntrinsicHighEnergyLimit = 100.0*GeV; // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); SetHighEnergyLimit(fIntrinsicHighEnergyLimit); //Calculate variable that will be used later on fPielr2 = pi*classic_electr_radius*classic_electr_radius; verboseLevel= 0; // Verbosity scale: // 0 = nothing // 1 = warning for energy non-conservation // 2 = details of energy budget // 3 = calculation of cross sections, file openings, sampling of atoms // 4 = entering in methods } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel() {;} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*, const G4DataVector&) { if (verboseLevel > 3) G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl; if(verboseLevel > 0) { G4cout << "Penelope Annihilation model is initialized " << G4endl << "Energy range: " << LowEnergyLimit() / keV << " keV - " << HighEnergyLimit() / GeV << " GeV" << G4endl; } if(isInitialised) return; fParticleChange = GetParticleChangeForGamma(); isInitialised = true; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom( const G4ParticleDefinition*, G4double energy, G4double Z, G4double, G4double, G4double) { if (verboseLevel > 3) G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << G4endl; G4double cs = Z*ComputeCrossSectionPerElectron(energy); if (verboseLevel > 2) G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << " = " << cs/barn << " barn" << G4endl; return cs; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector* fvect, const G4MaterialCutsCouple*, const G4DynamicParticle* aDynamicPositron, G4double, G4double) { // // Penelope model to sample final state for positron annihilation. // Target eletrons are assumed to be free and at rest. Binding effects enabling // one-photon annihilation are neglected. // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV // and isotropic angular distribution. // For annihilation in flight, it is used the theory from // W. Heitler, The quantum theory of radiation, Oxford University Press (1954) // The two photons can have different energy. The efficiency of the sampling algorithm // of the photon energy from the dSigma/dE distribution is practically 100% for // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy // of about 10 MeV. // The angle theta is kinematically linked to the photon energy, to ensure momentum // conservation. The angle phi is sampled isotropically for the first gamma. // if (verboseLevel > 3) G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl; G4double kineticEnergy = aDynamicPositron->GetKineticEnergy(); // kill primary fParticleChange->SetProposedKineticEnergy(0.); fParticleChange->ProposeTrackStatus(fStopAndKill); if (kineticEnergy == 0.0) { //Old AtRestDoIt G4double cosTheta = -1.0+2.0*G4UniformRand(); G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta); G4double phi = twopi*G4UniformRand(); G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta); G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(), direction, electron_mass_c2); G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(), -direction, electron_mass_c2); fvect->push_back(firstGamma); fvect->push_back(secondGamma); return; } //This is the "PostStep" case (annihilation in flight) G4ParticleMomentum positronDirection = aDynamicPositron->GetMomentumDirection(); G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2; G4double gamma21 = std::sqrt(gamma*gamma-1); G4double ani = 1.0+gamma; G4double chimin = 1.0/(ani+gamma21); G4double rchi = (1.0-chimin)/chimin; G4double gt0 = ani*ani-2.0; G4double test=0.0; G4double epsilon = 0; do{ epsilon = chimin*std::pow(rchi,G4UniformRand()); G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon); test = G4UniformRand()*gt0-reject; }while(test>0); G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2; G4double photon1Energy = epsilon*totalAvailableEnergy; G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy; G4double cosTheta1 = (ani-1.0/epsilon)/gamma21; G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21; //G4double localEnergyDeposit = 0.; G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1); G4double phi1 = twopi * G4UniformRand(); G4double dirx1 = sinTheta1 * std::cos(phi1); G4double diry1 = sinTheta1 * std::sin(phi1); G4double dirz1 = cosTheta1; G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2); G4double phi2 = phi1+pi; G4double dirx2 = sinTheta2 * std::cos(phi2); G4double diry2 = sinTheta2 * std::sin(phi2); G4double dirz2 = cosTheta2; G4ThreeVector photon1Direction (dirx1,diry1,dirz1); photon1Direction.rotateUz(positronDirection); // create G4DynamicParticle object for the particle1 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(), photon1Direction, photon1Energy); fvect->push_back(aParticle1); G4ThreeVector photon2Direction(dirx2,diry2,dirz2); photon2Direction.rotateUz(positronDirection); // create G4DynamicParticle object for the particle2 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(), photon2Direction, photon2Energy); fvect->push_back(aParticle2); if (verboseLevel > 1) { G4cout << "-----------------------------------------------------------" << G4endl; G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl; G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl; G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl; G4cout << "-----------------------------------------------------------" << G4endl; G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl; G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl; G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << " keV" << G4endl; G4cout << "-----------------------------------------------------------" << G4endl; } if (verboseLevel > 0) { G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy); if (energyDiff > 0.05*keV) G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << (photon1Energy+photon2Energy)/keV << " keV (final) vs. " << totalAvailableEnergy/keV << " keV (initial)" << G4endl; } return; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy) { // // Penelope model to calculate cross section for positron annihilation. // The annihilation cross section per electron is calculated according // to the Heitler formula // W. Heitler, The quantum theory of radiation, Oxford University Press (1954) // in the assumptions of electrons free and at rest. // G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2; G4double gamma2 = gamma*gamma; G4double f2 = gamma2-1.0; G4double f1 = std::sqrt(f2); G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2 - (gamma+3.0)/f1)/(gamma+1.0); return crossSection; }