| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eeToTwoGammaModel
33 //
34 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 //
36 // Creation date: 02.08.2004
37 //
38 // Modifications:
39 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
40 // 18-04-05 Compute CrossSectionPerVolume (V.Ivanchenko)
41 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
42 // 29-06-06 Fix problem for zero energy incident positron (V.Ivanchenko)
43 // 20-10-06 Add theGamma as a member (V.Ivanchenko)
44 // 18-01-20 Introduce thermal model of annihilation at rest (J.Allison)
45 //
46 //
47 // Class Description:
48 //
49 // Implementation of e+ annihilation into 2 gamma
50 //
51 // The secondaries Gamma energies are sampled using the Heitler cross section.
52 //
53 // A modified version of the random number techniques of Butcher & Messel
54 // is used (Nuc Phys 20(1960),15).
55 //
56 // GEANT4 internal units.
57 //
58 // Note 1: The initial electron is assumed free and at rest if atomic PDF
59 // is not defined
60 //
61 // Note 2: The annihilation processes producing one or more than two photons are
62 // ignored, as negligible compared to the two photons process.
63
64 //
65 // -------------------------------------------------------------------
66 //
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70 #include "G4eeToTwoGammaModel.hh"
71 #include "G4PhysicalConstants.hh"
72 #include "G4SystemOfUnits.hh"
73 #include "G4TrackStatus.hh"
74 #include "G4Electron.hh"
75 #include "G4Positron.hh"
76 #include "G4Gamma.hh"
77 #include "Randomize.hh"
78 #include "G4RandomDirection.hh"
79 #include "G4ParticleChangeForGamma.hh"
80 #include "G4EmParameters.hh"
81 #include "G4Log.hh"
82 #include "G4Exp.hh"
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
87 const G4String& nam)
88 : G4VEmModel(nam),
89 pi_rcl2(CLHEP::pi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius)
90 {
91 theGamma = G4Gamma::Gamma();
92 fParticleChange = nullptr;
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = default;
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
102 const G4DataVector&)
103 {
104 if (nullptr != fParticleChange) { return; }
105 fParticleChange = GetParticleChangeForGamma();
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109
110 G4double
111 G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(G4double kineticEnergy)
112 {
113 // Calculates the cross section per electron of annihilation into two photons
114 // from the Heilter formula.
115
116 G4double ekin = std::max(CLHEP::eV, kineticEnergy);
117
118 G4double tau = ekin/CLHEP::electron_mass_c2;
119 G4double gam = tau + 1.0;
120 G4double gamma2= gam*gam;
121 G4double bg2 = tau * (tau+2.0);
122 G4double bg = std::sqrt(bg2);
123
124 G4double cross = pi_rcl2*((gamma2+4*gam+1.)*G4Log(gam+bg) - (gam+3.)*bg)
125 / (bg2*(gam+1.));
126 return cross;
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130
131 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
132 const G4ParticleDefinition*,
133 G4double kineticEnergy, G4double Z,
134 G4double, G4double, G4double)
135 {
136 // Calculates the cross section per atom of annihilation into two photons
137 return Z*ComputeCrossSectionPerElectron(kineticEnergy);
138 }
139
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141
142 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
143 const G4Material* material,
144 const G4ParticleDefinition*,
145 G4double kineticEnergy,
146 G4double, G4double)
147 {
148 // Calculates the cross section per volume of annihilation into two photons
149 return material->GetElectronDensity()*ComputeCrossSectionPerElectron(kineticEnergy);
150 }
151
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153
154 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
155 // Nature 4065 (1947) 435.
156
157 void G4eeToTwoGammaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
158 const G4MaterialCutsCouple*,
159 const G4DynamicParticle* dp,
160 G4double,
161 G4double)
162 {
163 // kill primary positron
164 fParticleChange->SetProposedKineticEnergy(0.0);
165 fParticleChange->ProposeTrackStatus(fStopAndKill);
166
167 // Case at rest not considered anymore inside this model
168 G4LorentzVector lv(dp->GetMomentum(),
169 dp->GetKineticEnergy() + 2*CLHEP::electron_mass_c2);
170 G4double eGammaCMS = 0.5 * lv.mag();
171
172 G4ThreeVector dir1 = G4RandomDirection();
173 G4double phi = CLHEP::twopi * G4UniformRand();
174 G4double cosphi = std::cos(phi);
175 G4double sinphi = std::sin(phi);
176 G4ThreeVector pol1(cosphi, sinphi, 0.0);
177 pol1.rotateUz(dir1);
178 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCMS);
179
180 G4ThreeVector pol2(-sinphi, cosphi, 0.0);
181 pol2.rotateUz(dir1);
182
183 // transformation to lab system
184 lv1.boost(lv.boostVector());
185 lv -= lv1;
186
187 //!!! boost of polarisation vector is not yet implemented
188
189 // use constructors optimal for massless particle
190 auto aGamma1 = new G4DynamicParticle(G4Gamma::Gamma(), lv1.vect());
191 aGamma1->SetPolarization(pol1);
192 auto aGamma2 = new G4DynamicParticle(G4Gamma::Gamma(), lv.vect());
193 aGamma2->SetPolarization(pol2);
194
195 vdp->push_back(aGamma1);
196 vdp->push_back(aGamma2);
197 }
198
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200