54 :
G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),fIsInitialised(false)
86 G4cout <<
"Calling G4PenelopeAnnihilationModel::Initialise()" <<
G4endl;
93 G4cout <<
"Penelope Annihilation model is initialized " <<
G4endl
111 G4cout <<
"Calling G4PenelopeAnnihilationModel::InitialiseLocal()" <<
G4endl;
137 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" <<
143 G4cout <<
"Annihilation cross Section at " <<
energy/
keV <<
" keV for Z=" <<
Z <<
172 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" <<
G4endl;
180 if (kineticEnergy == 0.0)
184 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
186 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
192 fvect->push_back(firstGamma);
193 fvect->push_back(secondGamma);
201 G4double gamma21 = std::sqrt(gamma*gamma-1);
203 G4double chimin = 1.0/(ani+gamma21);
204 G4double rchi = (1.0-chimin)/chimin;
220 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
222 G4double dirx1 = sinTheta1 * std::cos(phi1);
223 G4double diry1 = sinTheta1 * std::sin(phi1);
226 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
228 G4double dirx2 = sinTheta2 * std::cos(phi2);
229 G4double diry2 = sinTheta2 * std::sin(phi2);
233 photon1Direction.
rotateUz(positronDirection);
238 fvect->push_back(aParticle1);
241 photon2Direction.
rotateUz(positronDirection);
246 fvect->push_back(aParticle2);
250 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
251 G4cout <<
"Energy balance from G4PenelopeAnnihilation" <<
G4endl;
252 G4cout <<
"Kinetic positron energy: " << kineticEnergy/
keV <<
" keV" <<
G4endl;
253 G4cout <<
"Total available energy: " << totalAvailableEnergy/
keV <<
" keV " <<
G4endl;
254 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
255 G4cout <<
"Photon energy 1: " << photon1Energy/
keV <<
" keV" <<
G4endl;
256 G4cout <<
"Photon energy 2: " << photon2Energy/
keV <<
" keV" <<
G4endl;
257 G4cout <<
"Total final state: " << (photon1Energy+photon2Energy)/
keV <<
259 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
263 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
264 if (energyDiff > 0.05*
keV)
265 G4cout <<
"Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
266 (photon1Energy+photon2Energy)/
keV <<
267 " keV (final) vs. " <<
268 totalAvailableEnergy/
keV <<
" keV (initial)" <<
G4endl;
289 - (gamma+3.0)/f1)/(gamma+1.0);
G4double epsilon(G4double density, G4double temperature)
G4double G4Log(G4double x)
static constexpr double twopi
static constexpr double barn
static constexpr double keV
static constexpr double eV
static constexpr double GeV
static constexpr double pi
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double ComputeCrossSectionPerElectron(G4double energy)
G4double fIntrinsicHighEnergyLimit
void SetParticle(const G4ParticleDefinition *)
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *) override
G4PenelopeAnnihilationModel(const G4ParticleDefinition *p=nullptr, const G4String &processName="PenAnnih")
const G4ParticleDefinition * fParticle
G4double fIntrinsicLowEnergyLimit
virtual ~G4PenelopeAnnihilationModel()
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
G4ParticleChangeForGamma * fParticleChange
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
int classic_electr_radius