| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eplusTo2or3GammaModel
33 //
34 // Author: Vladimir Ivanchenko and Omrame
35 //
36 // Creation date: 29.03.2018
37 //
38 //
39 // -------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oo
42 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
64
65 G4PhysicsVector* G4eplusTo2or3GammaModel::fCro
66 G4PhysicsVector* G4eplusTo2or3GammaModel::f3GP
67
68 G4eplusTo2or3GammaModel::G4eplusTo2or3GammaMod
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/e
84 fCrossSection = new G4PhysicsLogVector(emi
85 f3GProbability = new G4PhysicsLogVector(em
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 G4eplusTo2or3GammaModel::~G4eplusTo2or3GammaMo
92 {
93 if (IsMaster()) {
94 delete fCrossSection;
95 delete f3GProbability;
96 fCrossSection = nullptr;
97 f3GProbability = nullptr;
98 }
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void G4eplusTo2or3GammaModel::Initialise(const
104 const G4DataVector& cuts)
105 {
106 // here particle change is set for the tripl
107 if (nullptr == fParticleChange) {
108 fParticleChange = GetParticleChangeForGamm
109 }
110 // initialialise 3-gamma model before new ru
111 f3GModel->Initialise(p, cuts);
112 fGammaTh = G4EmParameters::Instance()->Lowes
113
114 // initialise vectors
115 if (IsMaster()) {
116 std::size_t num = fCrossSection->GetVector
117 for (std::size_t i=0; i<num; ++i) {
118 G4double e = fCrossSection->Energy(i);
119 G4double cs2 = ComputeCrossSectionPerEle
120 G4double cs3 = f3GModel->ComputeCrossSec
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........oo
132
133 G4double
134 G4eplusTo2or3GammaModel::ComputeCrossSectionPe
135 {
136 // Calculates the cross section per electron
137 // photons from the Heilter formula with the
138 // annihilation channel. (A.A.) rho is chang
139
140 G4double ekin = std::max(CLHEP::eV, kinEne
141 G4double tau = ekin/CLHEP::electron_mass_
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+
147 - (gam+3.)/(std::sqrt(gam*gam - 1.));
148 G4double eGammaCMS = CLHEP::electron_mass_c2
149 fDelta = std::max(fDeltaMin, fGammaTh/eGamma
150 f3GModel->SetDelta(fDelta);
151
152 static const G4double pir2 =
153 CLHEP::pi*CLHEP::classic_electr_radius*CLH
154 G4double cross = (pir2*rho + alpha_rcl2*2.*G
155
156 return cross;
157 }
158
159 //....oooOO0OOooo........oooOO0OOooo........oo
160
161 G4double G4eplusTo2or3GammaModel::ComputeCross
162 const G4Pa
163 G4double k
164 G4double, G4double, G4double)
165 {
166 // Calculates the cross section per atom of
167 G4double cross = Z*fCrossSection->Value(kine
168 return cross;
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oo
172
173 G4double G4eplusTo2or3GammaModel::CrossSection
174 const G4Material* material,
175 const G4ParticleDefinition*,
176 G4double kineticEnergy,
177 G4double, G4double)
178 {
179 // Calculates the cross section per volume o
180 G4double eDensity = material->GetElectronDen
181 G4double cross = eDensity*fCrossSection->Val
182 return cross;
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oo
186
187 // Polarisation of gamma according to M.H.L.Pr
188 // Nature 4065 (1947) 435.
189
190 void G4eplusTo2or3GammaModel::SampleSecondarie
191 std::vector<G4DynamicParticle*>*
192 const G4M
193 const G4D
194 G4double,
195 {
196 // kill primary positron
197 fParticleChange->SetProposedKineticEnergy(0.
198 fParticleChange->ProposeTrackStatus(fStopAnd
199
200 // Case at rest not considered anymore
201 G4double posiKinEnergy = dp->GetKineticEnerg
202 G4LorentzVector lv(dp->GetMomentum(),
203 posiKinEnergy + 2*CLHEP::electron_mas
204 G4double eGammaCMS = 0.5 * lv.mag();
205
206 if (G4UniformRand() < f3GProbability->Value(
207 fDelta = std::max(fDeltaMin, fGammaTh/eGam
208 f3GModel->SetDelta(fDelta);
209 f3GModel->SampleSecondaries(vdp, couple, d
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, eGammaCM
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 ye
229
230 // use constructors optimal for massless par
231 auto aGamma1 = new G4DynamicParticle(G4Gamma
232 aGamma1->SetPolarization(pol1);
233 auto aGamma2 = new G4DynamicParticle(G4Gamma
234 aGamma2->SetPolarization(pol2);
235
236 vdp->push_back(aGamma1);
237 vdp->push_back(aGamma2);
238 }
239
240 //....oooOO0OOooo........oooOO0OOooo........oo
241