| 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: G4eplusTo2or3GammaModel
33 //
34 // Author: Vladimir Ivanchenko and Omrame Kadri
35 //
36 // Creation date: 29.03.2018
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43
44
45 #include "G4eplusTo2or3GammaModel.hh"
46 #include "G4eplusTo3GammaOKVIModel.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4EmParameters.hh"
50 #include "G4TrackStatus.hh"
51 #include "G4Electron.hh"
52 #include "G4Positron.hh"
53 #include "G4Gamma.hh"
54 #include "G4DataVector.hh"
55 #include "G4PhysicsVector.hh"
56 #include "G4PhysicsLogVector.hh"
57 #include "G4RandomDirection.hh"
58 #include "Randomize.hh"
59 #include "G4ParticleChangeForGamma.hh"
60 #include "G4Log.hh"
61 #include "G4Exp.hh"
62
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65 G4PhysicsVector* G4eplusTo2or3GammaModel::fCrossSection = nullptr;
66 G4PhysicsVector* G4eplusTo2or3GammaModel::f3GProbability = nullptr;
67
68 G4eplusTo2or3GammaModel::G4eplusTo2or3GammaModel()
69 : G4VEmModel("eplusTo2or3gamma"),
70 fDeltaMin(0.001),
71 fDelta(fDeltaMin),
72 fGammaTh(CLHEP::MeV)
73 {
74 theGamma = G4Gamma::Gamma();
75 fParticleChange = nullptr;
76 f3GModel = new G4eplusTo3GammaOKVIModel();
77 SetTripletModel(f3GModel);
78
79 // instantiate vectors once
80 if (nullptr == fCrossSection) {
81 G4double emin = 10*CLHEP::eV;
82 G4double emax = 100*CLHEP::TeV;
83 G4int nbins = 20*G4lrint(std::log10(emax/emin));
84 fCrossSection = new G4PhysicsLogVector(emin, emax, nbins, true);
85 f3GProbability = new G4PhysicsLogVector(emin, emax, nbins, true);
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91 G4eplusTo2or3GammaModel::~G4eplusTo2or3GammaModel()
92 {
93 if (IsMaster()) {
94 delete fCrossSection;
95 delete f3GProbability;
96 fCrossSection = nullptr;
97 f3GProbability = nullptr;
98 }
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103 void G4eplusTo2or3GammaModel::Initialise(const G4ParticleDefinition* p,
104 const G4DataVector& cuts)
105 {
106 // here particle change is set for the triplet model
107 if (nullptr == fParticleChange) {
108 fParticleChange = GetParticleChangeForGamma();
109 }
110 // initialialise 3-gamma model before new run
111 f3GModel->Initialise(p, cuts);
112 fGammaTh = G4EmParameters::Instance()->LowestTripletEnergy();
113
114 // initialise vectors
115 if (IsMaster()) {
116 std::size_t num = fCrossSection->GetVectorLength();
117 for (std::size_t i=0; i<num; ++i) {
118 G4double e = fCrossSection->Energy(i);
119 G4double cs2 = ComputeCrossSectionPerElectron(e);
120 G4double cs3 = f3GModel->ComputeCrossSectionPerElectron(e);
121 cs2 += cs3;
122 fCrossSection->PutValue(i, cs2);
123 G4double y = (cs2 > 0.0) ? cs3/cs2 : 0.0;
124 f3GProbability->PutValue(i, y);
125 }
126 fCrossSection->FillSecondDerivatives();
127 f3GProbability->FillSecondDerivatives();
128 }
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133 G4double
134 G4eplusTo2or3GammaModel::ComputeCrossSectionPerElectron(G4double kinEnergy)
135 {
136 // Calculates the cross section per electron of annihilation into two
137 // photons from the Heilter formula with the radiation correction to 3 gamma
138 // annihilation channel. (A.A.) rho is changed
139
140 G4double ekin = std::max(CLHEP::eV, kinEnergy);
141 G4double tau = ekin/CLHEP::electron_mass_c2;
142 G4double gam = tau + 1.0;
143 G4double gamma2 = gam*gam;
144 G4double bg2 = tau * (tau+2.0);
145 G4double bg = std::sqrt(bg2);
146 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
147 - (gam+3.)/(std::sqrt(gam*gam - 1.));
148 G4double eGammaCMS = CLHEP::electron_mass_c2 * std::sqrt(0.5*(tau + 2.0));
149 fDelta = std::max(fDeltaMin, fGammaTh/eGammaCMS);
150 f3GModel->SetDelta(fDelta);
151
152 static const G4double pir2 =
153 CLHEP::pi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
154 G4double cross = (pir2*rho + alpha_rcl2*2.*G4Log(fDelta)*rho*rho)/(gam+1.);
155
156 return cross;
157 }
158
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
160
161 G4double G4eplusTo2or3GammaModel::ComputeCrossSectionPerAtom(
162 const G4ParticleDefinition*,
163 G4double kineticEnergy, G4double Z,
164 G4double, G4double, G4double)
165 {
166 // Calculates the cross section per atom of annihilation into two photons
167 G4double cross = Z*fCrossSection->Value(kineticEnergy);
168 return cross;
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172
173 G4double G4eplusTo2or3GammaModel::CrossSectionPerVolume(
174 const G4Material* material,
175 const G4ParticleDefinition*,
176 G4double kineticEnergy,
177 G4double, G4double)
178 {
179 // Calculates the cross section per volume of annihilation into two photons
180 G4double eDensity = material->GetElectronDensity();
181 G4double cross = eDensity*fCrossSection->Value(kineticEnergy);
182 return cross;
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
187 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
188 // Nature 4065 (1947) 435.
189
190 void G4eplusTo2or3GammaModel::SampleSecondaries(
191 std::vector<G4DynamicParticle*>* vdp,
192 const G4MaterialCutsCouple* couple,
193 const G4DynamicParticle* dp,
194 G4double, G4double)
195 {
196 // kill primary positron
197 fParticleChange->SetProposedKineticEnergy(0.0);
198 fParticleChange->ProposeTrackStatus(fStopAndKill);
199
200 // Case at rest not considered anymore
201 G4double posiKinEnergy = dp->GetKineticEnergy();
202 G4LorentzVector lv(dp->GetMomentum(),
203 posiKinEnergy + 2*CLHEP::electron_mass_c2);
204 G4double eGammaCMS = 0.5 * lv.mag();
205
206 if (G4UniformRand() < f3GProbability->Value(posiKinEnergy)) {
207 fDelta = std::max(fDeltaMin, fGammaTh/eGammaCMS);
208 f3GModel->SetDelta(fDelta);
209 f3GModel->SampleSecondaries(vdp, couple, dp);
210 return;
211 }
212
213 G4ThreeVector dir1 = G4RandomDirection();
214 G4double phi = CLHEP::twopi * G4UniformRand();
215 G4double cosphi = std::cos(phi);
216 G4double sinphi = std::sin(phi);
217 G4ThreeVector pol1(cosphi, sinphi, 0.0);
218 pol1.rotateUz(dir1);
219 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCMS);
220
221 G4ThreeVector pol2(-sinphi, cosphi, 0.0);
222 pol2.rotateUz(dir1);
223
224 // transformation to lab system
225 lv1.boost(lv.boostVector());
226 lv -= lv1;
227
228 //!!! boost of polarisation vector is not yet implemented
229
230 // use constructors optimal for massless particle
231 auto aGamma1 = new G4DynamicParticle(G4Gamma::Gamma(), lv1.vect());
232 aGamma1->SetPolarization(pol1);
233 auto aGamma2 = new G4DynamicParticle(G4Gamma::Gamma(), lv.vect());
234 aGamma2->SetPolarization(pol2);
235
236 vdp->push_back(aGamma1);
237 vdp->push_back(aGamma2);
238 }
239
240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
241