#include <G4PenelopeAnnihilationModel.hh>
Inheritance diagram for G4PenelopeAnnihilationModel:
Public Member Functions | |
| G4PenelopeAnnihilationModel (const G4ParticleDefinition *p=0, const G4String &processName="PenAnnih") | |
| virtual | ~G4PenelopeAnnihilationModel () |
| virtual void | Initialise (const G4ParticleDefinition *, const G4DataVector &) |
| virtual G4double | ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) |
| virtual void | SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) |
| void | SetVerbosityLevel (G4int lev) |
| G4int | GetVerbosityLevel () |
Protected Attributes | |
| G4ParticleChangeForGamma * | fParticleChange |
Definition at line 55 of file G4PenelopeAnnihilationModel.hh.
| G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel | ( | const G4ParticleDefinition * | p = 0, |
|
| const G4String & | processName = "PenAnnih" | |||
| ) |
Definition at line 50 of file G4PenelopeAnnihilationModel.cc.
References G4INCL::Math::pi, and G4VEmModel::SetHighEnergyLimit().
00052 :G4VEmModel(nam),fParticleChange(0),isInitialised(false) 00053 { 00054 fIntrinsicLowEnergyLimit = 0.0; 00055 fIntrinsicHighEnergyLimit = 100.0*GeV; 00056 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 00057 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 00058 00059 //Calculate variable that will be used later on 00060 fPielr2 = pi*classic_electr_radius*classic_electr_radius; 00061 00062 verboseLevel= 0; 00063 // Verbosity scale: 00064 // 0 = nothing 00065 // 1 = warning for energy non-conservation 00066 // 2 = details of energy budget 00067 // 3 = calculation of cross sections, file openings, sampling of atoms 00068 // 4 = entering in methods 00069 00070 }
| G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel | ( | ) | [virtual] |
| G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
| G4double | kinEnergy, | |||
| G4double | Z, | |||
| G4double | A = 0, |
|||
| G4double | cut = 0, |
|||
| G4double | emax = DBL_MAX | |||
| ) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 100 of file G4PenelopeAnnihilationModel.cc.
References G4cout, and G4endl.
00105 { 00106 if (verboseLevel > 3) 00107 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" << 00108 G4endl; 00109 00110 G4double cs = Z*ComputeCrossSectionPerElectron(energy); 00111 00112 if (verboseLevel > 2) 00113 G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z << 00114 " = " << cs/barn << " barn" << G4endl; 00115 return cs; 00116 }
| G4int G4PenelopeAnnihilationModel::GetVerbosityLevel | ( | ) | [inline] |
| void G4PenelopeAnnihilationModel::Initialise | ( | const G4ParticleDefinition * | , | |
| const G4DataVector & | ||||
| ) | [virtual] |
Implements G4VEmModel.
Definition at line 79 of file G4PenelopeAnnihilationModel.cc.
References fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), and G4VEmModel::LowEnergyLimit().
00081 { 00082 if (verboseLevel > 3) 00083 G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl; 00084 00085 if(verboseLevel > 0) { 00086 G4cout << "Penelope Annihilation model is initialized " << G4endl 00087 << "Energy range: " 00088 << LowEnergyLimit() / keV << " keV - " 00089 << HighEnergyLimit() / GeV << " GeV" 00090 << G4endl; 00091 } 00092 00093 if(isInitialised) return; 00094 fParticleChange = GetParticleChangeForGamma(); 00095 isInitialised = true; 00096 }
| void G4PenelopeAnnihilationModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
| const G4MaterialCutsCouple * | , | |||
| const G4DynamicParticle * | , | |||
| G4double | tmin, | |||
| G4double | maxEnergy | |||
| ) | [virtual] |
Implements G4VEmModel.
Definition at line 120 of file G4PenelopeAnnihilationModel.cc.
References fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4INCL::Math::pi, G4VParticleChange::ProposeTrackStatus(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
00125 { 00126 // 00127 // Penelope model to sample final state for positron annihilation. 00128 // Target eletrons are assumed to be free and at rest. Binding effects enabling 00129 // one-photon annihilation are neglected. 00130 // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV 00131 // and isotropic angular distribution. 00132 // For annihilation in flight, it is used the theory from 00133 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954) 00134 // The two photons can have different energy. The efficiency of the sampling algorithm 00135 // of the photon energy from the dSigma/dE distribution is practically 100% for 00136 // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy 00137 // of about 10 MeV. 00138 // The angle theta is kinematically linked to the photon energy, to ensure momentum 00139 // conservation. The angle phi is sampled isotropically for the first gamma. 00140 // 00141 if (verboseLevel > 3) 00142 G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl; 00143 00144 G4double kineticEnergy = aDynamicPositron->GetKineticEnergy(); 00145 00146 // kill primary 00147 fParticleChange->SetProposedKineticEnergy(0.); 00148 fParticleChange->ProposeTrackStatus(fStopAndKill); 00149 00150 if (kineticEnergy == 0.0) 00151 { 00152 //Old AtRestDoIt 00153 G4double cosTheta = -1.0+2.0*G4UniformRand(); 00154 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta); 00155 G4double phi = twopi*G4UniformRand(); 00156 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta); 00157 G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(), 00158 direction, electron_mass_c2); 00159 G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(), 00160 -direction, electron_mass_c2); 00161 00162 fvect->push_back(firstGamma); 00163 fvect->push_back(secondGamma); 00164 return; 00165 } 00166 00167 //This is the "PostStep" case (annihilation in flight) 00168 G4ParticleMomentum positronDirection = 00169 aDynamicPositron->GetMomentumDirection(); 00170 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2; 00171 G4double gamma21 = std::sqrt(gamma*gamma-1); 00172 G4double ani = 1.0+gamma; 00173 G4double chimin = 1.0/(ani+gamma21); 00174 G4double rchi = (1.0-chimin)/chimin; 00175 G4double gt0 = ani*ani-2.0; 00176 G4double test=0.0; 00177 G4double epsilon = 0; 00178 do{ 00179 epsilon = chimin*std::pow(rchi,G4UniformRand()); 00180 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon); 00181 test = G4UniformRand()*gt0-reject; 00182 }while(test>0); 00183 00184 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2; 00185 G4double photon1Energy = epsilon*totalAvailableEnergy; 00186 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy; 00187 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21; 00188 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21; 00189 00190 //G4double localEnergyDeposit = 0.; 00191 00192 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1); 00193 G4double phi1 = twopi * G4UniformRand(); 00194 G4double dirx1 = sinTheta1 * std::cos(phi1); 00195 G4double diry1 = sinTheta1 * std::sin(phi1); 00196 G4double dirz1 = cosTheta1; 00197 00198 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2); 00199 G4double phi2 = phi1+pi; 00200 G4double dirx2 = sinTheta2 * std::cos(phi2); 00201 G4double diry2 = sinTheta2 * std::sin(phi2); 00202 G4double dirz2 = cosTheta2; 00203 00204 G4ThreeVector photon1Direction (dirx1,diry1,dirz1); 00205 photon1Direction.rotateUz(positronDirection); 00206 // create G4DynamicParticle object for the particle1 00207 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(), 00208 photon1Direction, 00209 photon1Energy); 00210 fvect->push_back(aParticle1); 00211 00212 G4ThreeVector photon2Direction(dirx2,diry2,dirz2); 00213 photon2Direction.rotateUz(positronDirection); 00214 // create G4DynamicParticle object for the particle2 00215 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(), 00216 photon2Direction, 00217 photon2Energy); 00218 fvect->push_back(aParticle2); 00219 00220 if (verboseLevel > 1) 00221 { 00222 G4cout << "-----------------------------------------------------------" << G4endl; 00223 G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl; 00224 G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl; 00225 G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl; 00226 G4cout << "-----------------------------------------------------------" << G4endl; 00227 G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl; 00228 G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl; 00229 G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV << 00230 " keV" << G4endl; 00231 G4cout << "-----------------------------------------------------------" << G4endl; 00232 } 00233 if (verboseLevel > 0) 00234 { 00235 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy); 00236 if (energyDiff > 0.05*keV) 00237 G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " << 00238 (photon1Energy+photon2Energy)/keV << 00239 " keV (final) vs. " << 00240 totalAvailableEnergy/keV << " keV (initial)" << G4endl; 00241 } 00242 return; 00243 }
| void G4PenelopeAnnihilationModel::SetVerbosityLevel | ( | G4int | lev | ) | [inline] |
Definition at line 80 of file G4PenelopeAnnihilationModel.hh.
Referenced by Initialise(), and SampleSecondaries().
1.4.7