| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4eplusTo3GammaOKVIModel 32 // File name: G4eplusTo3GammaOKVIModel
33 // 33 //
34 // Authors: Andrei Alkin, Vladimir Ivanchenko 34 // Authors: Andrei Alkin, Vladimir Ivanchenko, Omrame Kadri
35 // 35 //
36 // Creation date: 29.03.2018 36 // Creation date: 29.03.2018
37 // 37 //
38 // ------------------------------------------- 38 // -------------------------------------------------------------------
39 // 39 //
40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 42
43 #include "G4eplusTo3GammaOKVIModel.hh" 43 #include "G4eplusTo3GammaOKVIModel.hh"
44 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
46 #include "G4EmParameters.hh" 46 #include "G4EmParameters.hh"
47 #include "G4TrackStatus.hh" 47 #include "G4TrackStatus.hh"
48 #include "G4Electron.hh" 48 #include "G4Electron.hh"
49 #include "G4Positron.hh" 49 #include "G4Positron.hh"
50 #include "G4Gamma.hh" 50 #include "G4Gamma.hh"
51 #include "Randomize.hh" 51 #include "Randomize.hh"
52 #include "G4ParticleChangeForGamma.hh" 52 #include "G4ParticleChangeForGamma.hh"
53 #include "G4Log.hh" 53 #include "G4Log.hh"
54 #include "G4Exp.hh" 54 #include "G4Exp.hh"
55 55
56 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 57
58 using namespace std; 58 using namespace std;
59 59
60 G4eplusTo3GammaOKVIModel::G4eplusTo3GammaOKVIM 60 G4eplusTo3GammaOKVIModel::G4eplusTo3GammaOKVIModel(const G4ParticleDefinition*,
61 61 const G4String& nam)
62 : G4VEmModel(nam), fDelta(0.001) 62 : G4VEmModel(nam), fDelta(0.001)
63 { 63 {
64 theGamma = G4Gamma::Gamma(); 64 theGamma = G4Gamma::Gamma();
>> 65 fParticleChange = nullptr;
65 } 66 }
66 67
67 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 69
69 G4eplusTo3GammaOKVIModel::~G4eplusTo3GammaOKVI << 70 G4eplusTo3GammaOKVIModel::~G4eplusTo3GammaOKVIModel()
>> 71 {}
70 72
71 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 74
73 void G4eplusTo3GammaOKVIModel::Initialise(cons 75 void G4eplusTo3GammaOKVIModel::Initialise(const G4ParticleDefinition*,
74 const G4DataVector&) 76 const G4DataVector&)
75 {} << 77 {
>> 78 // here particle change is set for the triplet model
>> 79 if(fParticleChange) { return; }
>> 80 fParticleChange = GetParticleChangeForGamma();
>> 81 }
76 82
77 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 84
79 // (A.A.) F_{ijk} calculation method 85 // (A.A.) F_{ijk} calculation method
80 G4double G4eplusTo3GammaOKVIModel::ComputeF(G4 86 G4double G4eplusTo3GammaOKVIModel::ComputeF(G4double fr1, G4double fr2,
81 G4 87 G4double fr3, G4double kinEnergy)
82 { 88 {
83 G4double ekin = std::max(eV,kinEnergy); 89 G4double ekin = std::max(eV,kinEnergy);
84 G4double tau = ekin/electron_mass_c2; 90 G4double tau = ekin/electron_mass_c2;
85 G4double gam = tau + 1.0; 91 G4double gam = tau + 1.0;
86 G4double gamma2 = gam*gam; 92 G4double gamma2 = gam*gam;
87 G4double bg2 = tau * (tau+2.0); 93 G4double bg2 = tau * (tau+2.0);
88 G4double bg = sqrt(bg2); 94 G4double bg = sqrt(bg2);
89 95
90 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+ 96 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
91 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.; 97 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.;
92 G4double border; 98 G4double border;
93 99
94 if(ekin < 500*MeV) { 100 if(ekin < 500*MeV) {
95 border = 1. - (electron_mass_c2)/(2*(ekin 101 border = 1. - (electron_mass_c2)/(2*(ekin + electron_mass_c2));
96 } else { 102 } else {
97 border = 1. - (100*electron_mass_c2)/(2*(e 103 border = 1. - (100*electron_mass_c2)/(2*(ekin + electron_mass_c2));
98 } 104 }
99 105
100 border = std::min(border, 0.9999); 106 border = std::min(border, 0.9999);
101 107
102 if (fr1>border) { fr1 = border; } 108 if (fr1>border) { fr1 = border; }
103 if (fr2>border) { fr2 = border; } 109 if (fr2>border) { fr2 = border; }
104 if (fr3>border) { fr3 = border; } 110 if (fr3>border) { fr3 = border; }
105 111
106 G4double fr1s = fr1*fr1; // "s" for "square 112 G4double fr1s = fr1*fr1; // "s" for "squared"
107 G4double fr2s = fr2*fr2; 113 G4double fr2s = fr2*fr2;
108 G4double fr3s = fr3*fr3; 114 G4double fr3s = fr3*fr3;
109 115
110 G4double aa = (1.-fr1)*(1.-fr2); 116 G4double aa = (1.-fr1)*(1.-fr2);
111 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2); 117 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2);
112 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2) 118 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa);
113 119
114 G4double fres = -rho*(1./fr1s + 1./fr2s) 120 G4double fres = -rho*(1./fr1s + 1./fr2s)
115 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/ 121 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2)))
116 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*( 122 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add;
117 123
118 return fres; 124 return fres;
119 } 125 }
120 126
121 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 128
123 // (A.A.) F_{ijk} calculation method 129 // (A.A.) F_{ijk} calculation method
124 G4double G4eplusTo3GammaOKVIModel::ComputeF0(G 130 G4double G4eplusTo3GammaOKVIModel::ComputeF0(G4double fr1, G4double fr2,
125 G 131 G4double fr3)
126 { 132 {
127 G4double tau = 0.0; 133 G4double tau = 0.0;
128 G4double gam = tau + 1.0; 134 G4double gam = tau + 1.0;
129 G4double gamma2 = gam*gam; 135 G4double gamma2 = gam*gam;
130 G4double bg2 = tau * (tau+2.0); 136 G4double bg2 = tau * (tau+2.0);
131 G4double bg = sqrt(bg2); 137 G4double bg = sqrt(bg2);
132 138
133 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+ 139 G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
134 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.; 140 - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.;
135 G4double border = 0.5; 141 G4double border = 0.5;
136 142
137 if (fr1>border) { fr1 = border; } 143 if (fr1>border) { fr1 = border; }
138 if (fr2>border) { fr2 = border; } 144 if (fr2>border) { fr2 = border; }
139 if (fr3>border) { fr3 = border; } 145 if (fr3>border) { fr3 = border; }
140 146
141 G4double fr1s = fr1*fr1; // "s" for "square 147 G4double fr1s = fr1*fr1; // "s" for "squared"
142 G4double fr2s = fr2*fr2; 148 G4double fr2s = fr2*fr2;
143 G4double fr3s = fr3*fr3; 149 G4double fr3s = fr3*fr3;
144 150
145 G4double aa = (1.-fr1)*(1.-fr2); 151 G4double aa = (1.-fr1)*(1.-fr2);
146 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2); 152 G4double ab = fr3s + (fr1-fr2)*(fr1-fr2);
147 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2) 153 G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa);
148 154
149 G4double fres = -rho*(1./fr1s + 1./fr2s) 155 G4double fres = -rho*(1./fr1s + 1./fr2s)
150 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/ 156 + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2)))
151 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*( 157 + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add;
152 158
153 return fres; 159 return fres;
154 } 160 }
155 161
156 //(A.A.) diff x-sections for maximum search an 162 //(A.A.) diff x-sections for maximum search and rejection
157 G4double G4eplusTo3GammaOKVIModel::ComputeFS(G 163 G4double G4eplusTo3GammaOKVIModel::ComputeFS(G4double fr1,
158 G4double fr2, G4double fr3, G4double 164 G4double fr2, G4double fr3, G4double kinEnergy)
159 { 165 {
160 G4double ekin = std::max(eV,kinEnergy); 166 G4double ekin = std::max(eV,kinEnergy);
161 G4double tau = ekin/electron_mass_c2; 167 G4double tau = ekin/electron_mass_c2;
162 G4double gam = tau + 1.0; 168 G4double gam = tau + 1.0;
163 169
164 G4double fsum = fr1*fr1*(ComputeF(fr1,fr2,fr 170 G4double fsum = fr1*fr1*(ComputeF(fr1,fr2,fr3,ekin) +
165 ComputeF(fr3,fr1,fr2,ekin) + 171 ComputeF(fr3,fr1,fr2,ekin) +
166 ComputeF(fr2,fr3,fr1,ekin)); 172 ComputeF(fr2,fr3,fr1,ekin));
167 173
168 G4double dcross = fsum/((3*fr1*fr1*(gam+1.)) 174 G4double dcross = fsum/((3*fr1*fr1*(gam+1.)));
169 175
170 return dcross; 176 return dcross;
171 } 177 }
172 178
173 //....oooOO0OOooo........oooOO0OOooo........oo 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174 180
175 G4double 181 G4double
176 G4eplusTo3GammaOKVIModel::ComputeCrossSectionP 182 G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerElectron(G4double kinEnergy)
177 { 183 {
178 // Calculates the cross section per electron 184 // Calculates the cross section per electron of annihilation into 3 photons
179 // from the Heilter formula. 185 // from the Heilter formula.
180 186
181 G4double ekin = std::max(CLHEP::eV, kinEne << 187 G4double ekin = std::max(eV,kinEnergy);
182 G4double tau = ekin/CLHEP::electron_mass_ << 188 G4double tau = ekin/electron_mass_c2;
183 G4double gam = tau + 1.0; 189 G4double gam = tau + 1.0;
184 G4double gamma2 = gam*gam; 190 G4double gamma2 = gam*gam;
185 G4double bg2 = tau * (tau+2.0); 191 G4double bg2 = tau * (tau+2.0);
186 G4double bg = sqrt(bg2); 192 G4double bg = sqrt(bg2);
187 193
188 G4double rho = (gamma2+4*gam+1.)*G4Log(gam+b 194 G4double rho = (gamma2+4*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
189 - (gam+3.)/(sqrt(gam*gam - 1.)); 195 - (gam+3.)/(sqrt(gam*gam - 1.));
190 196
191 G4double cross = alpha_rcl2*(4.2 - (2.*G4Log 197 G4double cross = alpha_rcl2*(4.2 - (2.*G4Log(fDelta)+1.)*rho*rho)/(gam+1.);
192 return cross; 198 return cross;
193 } 199 }
194 200
195 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 202
197 G4double G4eplusTo3GammaOKVIModel::ComputeCros 203 G4double G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerAtom(
198 const G4Pa 204 const G4ParticleDefinition*,
199 G4double k 205 G4double kineticEnergy, G4double Z,
200 G4double, G4double, G4double) 206 G4double, G4double, G4double)
201 { 207 {
>> 208 // Calculates the cross section per atom of annihilation into two photons
>> 209
>> 210
>> 211
202 G4double cross = Z*ComputeCrossSectionPerEle 212 G4double cross = Z*ComputeCrossSectionPerElectron(kineticEnergy);
203 return cross; 213 return cross;
204 } 214 }
205 215
206 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
207 217
208 G4double G4eplusTo3GammaOKVIModel::CrossSectio 218 G4double G4eplusTo3GammaOKVIModel::CrossSectionPerVolume(
209 const G4Material* material, 219 const G4Material* material,
210 const G4ParticleDefinition*, 220 const G4ParticleDefinition*,
211 G4double kineticEnergy, 221 G4double kineticEnergy,
212 G4double, G4double) 222 G4double, G4double)
213 { 223 {
214 // Calculates the cross section per volume o 224 // Calculates the cross section per volume of annihilation into two photons
215 225
216 G4double eDensity = material->GetElectronDen 226 G4double eDensity = material->GetElectronDensity();
217 G4double cross = eDensity*ComputeCrossSectio 227 G4double cross = eDensity*ComputeCrossSectionPerElectron(kineticEnergy);
218 return cross; 228 return cross;
219 } 229 }
220 230
221 //....oooOO0OOooo........oooOO0OOooo........oo 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 232
223 // Polarisation of gamma according to M.H.L.Pr 233 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
224 // Nature 4065 (1947) 435. 234 // Nature 4065 (1947) 435.
225 235
226 void 236 void
227 G4eplusTo3GammaOKVIModel::SampleSecondaries(ve 237 G4eplusTo3GammaOKVIModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
228 const G4MaterialCutsCouple*, 238 const G4MaterialCutsCouple*,
229 const G4DynamicParticle* dp, 239 const G4DynamicParticle* dp,
230 G4double, G4double) 240 G4double, G4double)
231 { 241 {
232 // let us perform sampling in C.M.S. referen << 242
233 G4double posiKinEnergy = dp->GetKineticEnerg 243 G4double posiKinEnergy = dp->GetKineticEnergy();
234 G4LorentzVector lv(dp->GetMomentum(), << 244 G4DynamicParticle *aGamma1, *aGamma2;
235 posiKinEnergy + 2*CLHEP::electron_mas << 245 G4DynamicParticle* aGamma3 = nullptr;
236 G4double eGammaCMS = 0.5 * lv.mag(); << 246 G4double border;
237 << 247
238 // the limit value fDelta is defined by a cl << 248 if(posiKinEnergy < 500*MeV) {
239 // thickness of border defined by C.M.S. ene << 249 border = 1. - (electron_mass_c2)/(2*(posiKinEnergy + electron_mass_c2));
240 G4double border = << 250 } else {
241 1.0 - std::min(std::max(CLHEP::electron_ma << 251 border = 1. - (100*electron_mass_c2)/(2*(posiKinEnergy + electron_mass_c2));
>> 252 }
>> 253 border = std::min(border, 0.9999);
242 254
243 CLHEP::HepRandomEngine* rndmEngine = G4Rando 255 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
244 256
245 G4ThreeVector posiDirection = dp->GetMomentu << 257 // Case at rest
246 << 258 if(posiKinEnergy == 0.0) {
247 // (A.A.) LIMITS FOR 1st GAMMA << 259 G4double cost = 2.*rndmEngine->flat()-1.;
248 G4double xmin = 0.01; << 260 G4double sint = sqrt((1. - cost)*(1. + cost));
249 G4double xmax = 0.667; // CHANGE to 3/2 << 261 G4double phi = twopi * rndmEngine->flat();
>> 262 G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
>> 263 phi = twopi * rndmEngine->flat();
>> 264 G4double cosphi = cos(phi);
>> 265 G4double sinphi = sin(phi);
>> 266 G4ThreeVector pol(cosphi, sinphi, 0.0);
>> 267 pol.rotateUz(dir);
>> 268 aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
>> 269 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
>> 270 aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
>> 271 pol.set(-sinphi, cosphi, 0.0);
>> 272 pol.rotateUz(dir);
>> 273 aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
>> 274
>> 275 } else {
>> 276
>> 277 G4ThreeVector posiDirection = dp->GetMomentumDirection();
>> 278
>> 279 // (A.A.) LIMITS FOR 1st GAMMA
>> 280 G4double xmin = 0.01;
>> 281 G4double xmax = 0.667; // CHANGE to 3/2
250 282
251 G4double d1, d0, x1, x2, dmax, x2min; << 283 G4double d1, d0, x1, x2, dmax, x2min;
>> 284
>> 285 // (A.A.) sampling of x1 x2 x3 (whole cycle of rejection)
>> 286 do {
>> 287 x1 = 1/((1/xmin) - ((1/xmin)-(1/xmax))*rndmEngine->flat());
>> 288 dmax = ComputeFS(posiKinEnergy, x1,1.-x1,border);
>> 289 x2min = 1.-x1;
>> 290 x2 = 1 - rndmEngine->flat()*(1-x2min);
>> 291 d1 = dmax*rndmEngine->flat();
>> 292 d0 = ComputeFS(posiKinEnergy,x1,x2,2-x1-x2);
>> 293 }
>> 294 while(d0 < d1);
>> 295
>> 296 G4double x3 = 2 - x1 - x2;
>> 297
>> 298 //
>> 299 // angles between Gammas
>> 300 //
>> 301
>> 302 G4double psi13 = 2*asin(sqrt(std::abs((x1+x3-1)/(x1*x3))));
>> 303 G4double psi12 = 2*asin(sqrt(std::abs((x1+x2-1)/(x1*x2))));
>> 304
>> 305 // sin^t
>> 306
>> 307 //G4double phi = twopi * rndmEngine->flat();
>> 308 //G4double psi = acos(x3); // Angle of the plane
>> 309
>> 310 //
>> 311 // kinematic of the created pair
>> 312 //
>> 313
>> 314 G4double TotalAvailableEnergy = posiKinEnergy + 2.0*electron_mass_c2;
>> 315
>> 316 G4double phot1Energy = 0.5*x1*TotalAvailableEnergy;
>> 317 G4double phot2Energy = 0.5*x2*TotalAvailableEnergy;
>> 318 G4double phot3Energy = 0.5*x3*TotalAvailableEnergy;
252 319
253 // (A.A.) sampling of x1 x2 x3 (whole cycle <<
254 do { <<
255 x1 = 1./((1./xmin) - ((1./xmin)-(1./xmax)) <<
256 dmax = ComputeFS(eGammaCMS, x1, 1.-x1, bor <<
257 x2min = 1. - x1; <<
258 x2 = 1 - rndmEngine->flat()*(1. - x2min); <<
259 d1 = dmax*rndmEngine->flat(); <<
260 d0 = ComputeFS(eGammaCMS, x1, x2, 2.-x1-x2 <<
261 } <<
262 while(d0 < d1); <<
263 <<
264 G4double x3 = 2 - x1 - x2; <<
265 // <<
266 // angles between Gammas <<
267 // <<
268 G4double psi13 = 2*std::asin(std::sqrt(std:: <<
269 G4double psi12 = 2*std::asin(std::sqrt(std:: <<
270 // <<
271 // kinematic of the created pair <<
272 // <<
273 G4double phot1Energy = x1*eGammaCMS; <<
274 G4double phot2Energy = x2*eGammaCMS; <<
275 G4double phot3Energy = x3*eGammaCMS; <<
276 320
277 // DIRECTIONS << 321 // DIRECTIONS
278 322
279 // The azimuthal angles of q1 and q3 with re << 323 // The azimuthal angles of ql and q3 with respect to some plane
280 // through the beam axis are generated at ra << 324 // through the beam axis are generated at random.
>> 325
>> 326 G4ThreeVector phot1Direction(0, 0, 1);
>> 327 G4ThreeVector phot2Direction(0, sin(psi12), cos(psi12));
>> 328 G4ThreeVector phot3Direction(0, sin(psi13), cos(psi13));
>> 329
>> 330 phot1Direction.rotateUz(posiDirection);
>> 331 phot2Direction.rotateUz(posiDirection);
>> 332 phot3Direction.rotateUz(posiDirection);
>> 333
>> 334 aGamma1 = new G4DynamicParticle (theGamma,phot1Direction, phot1Energy);
>> 335 aGamma2 = new G4DynamicParticle (theGamma,phot2Direction, phot2Energy);
>> 336 aGamma3 = new G4DynamicParticle (theGamma,phot3Direction, phot3Energy);
281 337
282 G4ThreeVector phot1Direction(0, 0, 1); <<
283 G4ThreeVector phot2Direction(0, std::sin(psi <<
284 G4ThreeVector phot3Direction(0, std::sin(psi <<
285 <<
286 G4LorentzVector lv1(phot1Energy*phot1Directi <<
287 G4LorentzVector lv2(phot2Energy*phot2Directi <<
288 G4LorentzVector lv3(phot3Energy*phot3Directi <<
289 <<
290 auto boostV = lv.boostVector(); <<
291 lv1.boost(boostV); <<
292 lv2.boost(boostV); <<
293 lv3.boost(boostV); <<
294 <<
295 auto aGamma1 = new G4DynamicParticle (theGam <<
296 auto aGamma2 = new G4DynamicParticle (theGam <<
297 auto aGamma3 = new G4DynamicParticle (theGam <<
298 338
299 //!!! POLARIZATION - not yet implemented << 339 //POLARIZATION - ???
>> 340 /*
>> 341
>> 342
>> 343 phi = twopi * rndmEngine->flat();
>> 344 G4double cosphi = cos(phi);
>> 345 G4double sinphi = sin(phi);
>> 346 G4ThreeVector pol(cosphi, sinphi, 0.0);
>> 347 pol.rotateUz(phot1Direction);
>> 348 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
>> 349
>> 350 G4double phot2Energy =(1.-epsil)*TotalAvailableEnergy;
>> 351 G4double posiP= sqrt(posiKinEnergy*(posiKinEnergy+2.*electron_mass_c2));
>> 352 G4ThreeVector dir = posiDirection*posiP - phot1Direction*phot1Energy;
>> 353 G4ThreeVector phot2Direction = dir.unit();
>> 354
>> 355 // create G4DynamicParticle object for the particle2
>> 356 aGamma2 = new G4DynamicParticle (theGamma,phot2Direction, phot2Energy);
>> 357
>> 358 //!!! likely problematic direction to be checked
>> 359 pol.set(-sinphi, cosphi, 0.0);
>> 360 pol.rotateUz(phot1Direction);
>> 361 cost = pol*phot2Direction;
>> 362 pol -= cost*phot2Direction;
>> 363 pol = pol.unit();
>> 364 aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
>> 365
>> 366 */
>> 367
>> 368 }
>> 369 /*
>> 370 G4cout << "Annihilation in fly: e0= " << posiKinEnergy
>> 371 << " m= " << electron_mass_c2
>> 372 << " e1= " << phot1Energy
>> 373 << " e2= " << phot2Energy << " dir= " << dir
>> 374 << " -> " << phot1Direction << " "
>> 375 << phot2Direction << G4endl;
>> 376 */
300 377
301 vdp->push_back(aGamma1); 378 vdp->push_back(aGamma1);
302 vdp->push_back(aGamma2); 379 vdp->push_back(aGamma2);
303 vdp->push_back(aGamma3); << 380 if(aGamma3 != nullptr) { vdp->push_back(aGamma3); }
>> 381
>> 382 // kill primary positron
>> 383 fParticleChange->SetProposedKineticEnergy(0.0);
>> 384 fParticleChange->ProposeTrackStatus(fStopAndKill);
304 } 385 }
305 386
306 //....oooOO0OOooo........oooOO0OOooo........oo 387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 388