ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eplusTo3GammaOKVIModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eplusTo3GammaOKVIModel.cc
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: G4eplusTo3GammaOKVIModel
33 //
34 // Authors: Andrei Alkin, Vladimir Ivanchenko, Omrame Kadri
35 //
36 // Creation date: 29.03.2018
37 //
38 // -------------------------------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 #include "G4eplusTo3GammaOKVIModel.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4EmParameters.hh"
47 #include "G4TrackStatus.hh"
48 #include "G4Electron.hh"
49 #include "G4Positron.hh"
50 #include "G4Gamma.hh"
51 #include "Randomize.hh"
52 #include "G4ParticleChangeForGamma.hh"
53 #include "G4Log.hh"
54 #include "G4Exp.hh"
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 
58 using namespace std;
59 
60 G4eplusTo3GammaOKVIModel::G4eplusTo3GammaOKVIModel(const G4ParticleDefinition*,
61  const G4String& nam)
62  : G4VEmModel(nam), fDelta(0.001)
63 {
64  theGamma = G4Gamma::Gamma();
65  fParticleChange = nullptr;
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
70 G4eplusTo3GammaOKVIModel::~G4eplusTo3GammaOKVIModel()
71 {}
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 void G4eplusTo3GammaOKVIModel::Initialise(const G4ParticleDefinition*,
76  const G4DataVector&)
77 {
78  // here particle change is set for the triplet model
79  if(fParticleChange) { return; }
80  fParticleChange = GetParticleChangeForGamma();
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 // (A.A.) F_{ijk} calculation method
86 G4double G4eplusTo3GammaOKVIModel::ComputeF(G4double fr1, G4double fr2,
87  G4double fr3, G4double kinEnergy)
88 {
89  G4double ekin = std::max(eV,kinEnergy);
90  G4double tau = ekin/electron_mass_c2;
91  G4double gam = tau + 1.0;
92  G4double gamma2 = gam*gam;
93  G4double bg2 = tau * (tau+2.0);
94  G4double bg = sqrt(bg2);
95 
96  G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
97  - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.;
98  G4double border;
99 
100  if(ekin < 500*MeV) {
101  border = 1. - (electron_mass_c2)/(2*(ekin + electron_mass_c2));
102  } else {
103  border = 1. - (100*electron_mass_c2)/(2*(ekin + electron_mass_c2));
104  }
105 
106  border = std::min(border, 0.9999);
107 
108  if (fr1>border) { fr1 = border; }
109  if (fr2>border) { fr2 = border; }
110  if (fr3>border) { fr3 = border; }
111 
112  G4double fr1s = fr1*fr1; // "s" for "squared"
113  G4double fr2s = fr2*fr2;
114  G4double fr3s = fr3*fr3;
115 
116  G4double aa = (1.-fr1)*(1.-fr2);
117  G4double ab = fr3s + (fr1-fr2)*(fr1-fr2);
118  G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa);
119 
120  G4double fres = -rho*(1./fr1s + 1./fr2s)
121  + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2)))
122  + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add;
123 
124  return fres;
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 
129 // (A.A.) F_{ijk} calculation method
130 G4double G4eplusTo3GammaOKVIModel::ComputeF0(G4double fr1, G4double fr2,
131  G4double fr3)
132 {
133  G4double tau = 0.0;
134  G4double gam = tau + 1.0;
135  G4double gamma2 = gam*gam;
136  G4double bg2 = tau * (tau+2.0);
137  G4double bg = sqrt(bg2);
138 
139  G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
140  - (gam+3.)/(sqrt(gam*gam - 1.)) + 1.;
141  G4double border = 0.5;
142 
143  if (fr1>border) { fr1 = border; }
144  if (fr2>border) { fr2 = border; }
145  if (fr3>border) { fr3 = border; }
146 
147  G4double fr1s = fr1*fr1; // "s" for "squared"
148  G4double fr2s = fr2*fr2;
149  G4double fr3s = fr3*fr3;
150 
151  G4double aa = (1.-fr1)*(1.-fr2);
152  G4double ab = fr3s + (fr1-fr2)*(fr1-fr2);
153  G4double add= ((1.-fr1)*(1.-fr1) + (1.-fr2)*(1.-fr2))/(fr3s*aa);
154 
155  G4double fres = -rho*(1./fr1s + 1./fr2s)
156  + (ab/(2.*(fr1*fr2*aa)))*(G4Log(2.*gam*aa/(fr1*fr2)))
157  + (ab/(2.*fr1*fr2*(1-fr3)))*G4Log(2.*gam*(1.-fr3)/(fr1*fr2)) - add;
158 
159  return fres;
160 }
161 
162 //(A.A.) diff x-sections for maximum search and rejection
163 G4double G4eplusTo3GammaOKVIModel::ComputeFS(G4double fr1,
164  G4double fr2, G4double fr3, G4double kinEnergy)
165 {
166  G4double ekin = std::max(eV,kinEnergy);
167  G4double tau = ekin/electron_mass_c2;
168  G4double gam = tau + 1.0;
169 
170  G4double fsum = fr1*fr1*(ComputeF(fr1,fr2,fr3,ekin) +
171  ComputeF(fr3,fr1,fr2,ekin) +
172  ComputeF(fr2,fr3,fr1,ekin));
173 
174  G4double dcross = fsum/((3*fr1*fr1*(gam+1.)));
175 
176  return dcross;
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 
181 G4double
182 G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerElectron(G4double kinEnergy)
183 {
184  // Calculates the cross section per electron of annihilation into 3 photons
185  // from the Heilter formula.
186 
187  G4double ekin = std::max(eV,kinEnergy);
188  G4double tau = ekin/electron_mass_c2;
189  G4double gam = tau + 1.0;
190  G4double gamma2 = gam*gam;
191  G4double bg2 = tau * (tau+2.0);
192  G4double bg = sqrt(bg2);
193 
194  G4double rho = (gamma2+4*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
195  - (gam+3.)/(sqrt(gam*gam - 1.));
196 
197  G4double cross = alpha_rcl2*(4.2 - (2.*G4Log(fDelta)+1.)*rho*rho)/(gam+1.);
198  return cross;
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202 
203 G4double G4eplusTo3GammaOKVIModel::ComputeCrossSectionPerAtom(
204  const G4ParticleDefinition*,
205  G4double kineticEnergy, G4double Z,
206  G4double, G4double, G4double)
207 {
208  // Calculates the cross section per atom of annihilation into two photons
209 
210 
211 
212  G4double cross = Z*ComputeCrossSectionPerElectron(kineticEnergy);
213  return cross;
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 
218 G4double G4eplusTo3GammaOKVIModel::CrossSectionPerVolume(
219  const G4Material* material,
220  const G4ParticleDefinition*,
221  G4double kineticEnergy,
222  G4double, G4double)
223 {
224  // Calculates the cross section per volume of annihilation into two photons
225 
226  G4double eDensity = material->GetElectronDensity();
227  G4double cross = eDensity*ComputeCrossSectionPerElectron(kineticEnergy);
228  return cross;
229 }
230 
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232 
233 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
234 // Nature 4065 (1947) 435.
235 
236 void
237 G4eplusTo3GammaOKVIModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
238  const G4MaterialCutsCouple*,
239  const G4DynamicParticle* dp,
240  G4double, G4double)
241 {
242 
243  G4double posiKinEnergy = dp->GetKineticEnergy();
244  G4DynamicParticle *aGamma1, *aGamma2, *aGamma3;
245  G4double border;
246 
247  if(posiKinEnergy < 500*MeV) {
248  border = 1. - (electron_mass_c2)/(2*(posiKinEnergy + electron_mass_c2));
249  } else {
250  border = 1. - (100*electron_mass_c2)/(2*(posiKinEnergy + electron_mass_c2));
251  }
252  border = std::min(border, 0.9999);
253 
254  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
255 
256  // Case at rest
257  if(posiKinEnergy == 0.0) {
258  G4double cost = 2.*rndmEngine->flat()-1.;
259  G4double sint = sqrt((1. - cost)*(1. + cost));
260  G4double phi = twopi * rndmEngine->flat();
261  G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
262  phi = twopi * rndmEngine->flat();
263  G4double cosphi = cos(phi);
264  G4double sinphi = sin(phi);
265  G4ThreeVector pol(cosphi, sinphi, 0.0);
266  pol.rotateUz(dir);
267  aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
268  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
269  aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
270  pol.set(-sinphi, cosphi, 0.0);
271  pol.rotateUz(dir);
272  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
273 
274  } else {
275 
276  G4ThreeVector posiDirection = dp->GetMomentumDirection();
277 
278  // (A.A.) LIMITS FOR 1st GAMMA
279  G4double xmin = 0.01;
280  G4double xmax = 0.667; // CHANGE to 3/2
281 
282  G4double d1, d0, x1, x2, dmax, x2min;
283 
284  // (A.A.) sampling of x1 x2 x3 (whole cycle of rejection)
285  do {
286  x1 = 1/((1/xmin) - ((1/xmin)-(1/xmax))*rndmEngine->flat());
287  dmax = ComputeFS(posiKinEnergy, x1,1.-x1,border);
288  x2min = 1.-x1;
289  x2 = 1 - rndmEngine->flat()*(1-x2min);
290  d1 = dmax*rndmEngine->flat();
291  d0 = ComputeFS(posiKinEnergy,x1,x2,2-x1-x2);
292  }
293  while(d0 < d1);
294 
295  G4double x3 = 2 - x1 - x2;
296 
297  //
298  // angles between Gammas
299  //
300 
301  G4double psi13 = 2*asin(sqrt(std::abs((x1+x3-1)/(x1*x3))));
302  G4double psi12 = 2*asin(sqrt(std::abs((x1+x2-1)/(x1*x2))));
303 
304  // sin^t
305 
306  //G4double phi = twopi * rndmEngine->flat();
307  //G4double psi = acos(x3); // Angle of the plane
308 
309  //
310  // kinematic of the created pair
311  //
312 
313  G4double TotalAvailableEnergy = posiKinEnergy + 2.0*electron_mass_c2;
314 
315  G4double phot1Energy = 0.5*x1*TotalAvailableEnergy;
316  G4double phot2Energy = 0.5*x2*TotalAvailableEnergy;
317  G4double phot3Energy = 0.5*x3*TotalAvailableEnergy;
318 
319 
320  // DIRECTIONS
321 
322  // The azimuthal angles of ql and q3 with respect to some plane
323  // through the beam axis are generated at random.
324 
325  G4ThreeVector phot1Direction(0, 0, 1);
326  G4ThreeVector phot2Direction(0, sin(psi12), cos(psi12));
327  G4ThreeVector phot3Direction(0, sin(psi13), cos(psi13));
328 
329  phot1Direction.rotateUz(posiDirection);
330  phot2Direction.rotateUz(posiDirection);
331  phot3Direction.rotateUz(posiDirection);
332 
333  aGamma1 = new G4DynamicParticle (theGamma,phot1Direction, phot1Energy);
334  aGamma2 = new G4DynamicParticle (theGamma,phot2Direction, phot2Energy);
335  aGamma3 = new G4DynamicParticle (theGamma,phot3Direction, phot3Energy);
336 
337 
338  //POLARIZATION - ???
339  /*
340 
341 
342  phi = twopi * rndmEngine->flat();
343  G4double cosphi = cos(phi);
344  G4double sinphi = sin(phi);
345  G4ThreeVector pol(cosphi, sinphi, 0.0);
346  pol.rotateUz(phot1Direction);
347  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
348 
349  G4double phot2Energy =(1.-epsil)*TotalAvailableEnergy;
350  G4double posiP= sqrt(posiKinEnergy*(posiKinEnergy+2.*electron_mass_c2));
351  G4ThreeVector dir = posiDirection*posiP - phot1Direction*phot1Energy;
352  G4ThreeVector phot2Direction = dir.unit();
353 
354  // create G4DynamicParticle object for the particle2
355  aGamma2 = new G4DynamicParticle (theGamma,phot2Direction, phot2Energy);
356 
358  pol.set(-sinphi, cosphi, 0.0);
359  pol.rotateUz(phot1Direction);
360  cost = pol*phot2Direction;
361  pol -= cost*phot2Direction;
362  pol = pol.unit();
363  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
364 
365  */
366 
367  }
368  /*
369  G4cout << "Annihilation in fly: e0= " << posiKinEnergy
370  << " m= " << electron_mass_c2
371  << " e1= " << phot1Energy
372  << " e2= " << phot2Energy << " dir= " << dir
373  << " -> " << phot1Direction << " "
374  << phot2Direction << G4endl;
375  */
376 
377  vdp->push_back(aGamma1);
378  vdp->push_back(aGamma2);
379  vdp->push_back(aGamma3);
380 
381  // kill primary positron
382  fParticleChange->SetProposedKineticEnergy(0.0);
383  fParticleChange->ProposeTrackStatus(fStopAndKill);
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....