ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eeToTwoGammaModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eeToTwoGammaModel.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: 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 //
45 //
46 // Class Description:
47 //
48 // Implementation of e+ annihilation into 2 gamma
49 //
50 // The secondaries Gamma energies are sampled using the Heitler cross section.
51 //
52 // A modified version of the random number techniques of Butcher & Messel
53 // is used (Nuc Phys 20(1960),15).
54 //
55 // GEANT4 internal units.
56 //
57 // Note 1: The initial electron is assumed free and at rest.
58 //
59 // Note 2: The annihilation processes producing one or more than two photons are
60 // ignored, as negligible compared to the two photons process.
61 
62 
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 "G4ParticleChangeForGamma.hh"
79 #include "G4Log.hh"
80 #include "G4Exp.hh"
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 using namespace std;
85 
86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
87  const G4String& nam)
88  : G4VEmModel(nam),
89  pi_rcl2(pi*classic_electr_radius*classic_electr_radius)
90 {
91  theGamma = G4Gamma::Gamma();
92  fParticleChange = nullptr;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 
97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel()
98 {}
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
103  const G4DataVector&)
104 {
105  if(fParticleChange) { return; }
106  fParticleChange = GetParticleChangeForGamma();
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 G4double
112 G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(G4double kineticEnergy)
113 {
114  // Calculates the cross section per electron of annihilation into two photons
115  // from the Heilter formula.
116 
117  G4double ekin = std::max(eV,kineticEnergy);
118 
119  G4double tau = ekin/electron_mass_c2;
120  G4double gam = tau + 1.0;
121  G4double gamma2= gam*gam;
122  G4double bg2 = tau * (tau+2.0);
123  G4double bg = sqrt(bg2);
124 
125  G4double cross = pi_rcl2*((gamma2+4*gam+1.)*G4Log(gam+bg) - (gam+3.)*bg)
126  / (bg2*(gam+1.));
127  return cross;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 
132 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
133  const G4ParticleDefinition*,
134  G4double kineticEnergy, G4double Z,
135  G4double, G4double, G4double)
136 {
137  // Calculates the cross section per atom of annihilation into two photons
138 
139  G4double cross = Z*ComputeCrossSectionPerElectron(kineticEnergy);
140  return cross;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
146  const G4Material* material,
147  const G4ParticleDefinition*,
148  G4double kineticEnergy,
149  G4double, G4double)
150 {
151  // Calculates the cross section per volume of annihilation into two photons
152 
153  G4double eDensity = material->GetElectronDensity();
154  G4double cross = eDensity*ComputeCrossSectionPerElectron(kineticEnergy);
155  return cross;
156 }
157 
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159 
160 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
161 // Nature 4065 (1947) 435.
162 
163 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
164  const G4MaterialCutsCouple*,
165  const G4DynamicParticle* dp,
166  G4double,
167  G4double)
168 {
169  G4double posiKinEnergy = dp->GetKineticEnergy();
170  G4DynamicParticle *aGamma1, *aGamma2;
171 
172  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
173 
174  // Case at rest
175  if(posiKinEnergy == 0.0) {
176  G4double cost = 2.*rndmEngine->flat()-1.;
177  G4double sint = sqrt((1. - cost)*(1. + cost));
178  G4double phi = twopi * rndmEngine->flat();
179  G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
180  phi = twopi * rndmEngine->flat();
181  G4double cosphi = cos(phi);
182  G4double sinphi = sin(phi);
183  G4ThreeVector pol(cosphi, sinphi, 0.0);
184  pol.rotateUz(dir);
185  aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
186  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
187  aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
188  pol.set(-sinphi, cosphi, 0.0);
189  pol.rotateUz(dir);
190  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
191  /*
192  G4cout << "Annihilation at rest fly: e0= " << " dir= " << dir
193  << G4endl;
194  */
195  } else {
196 
197  G4ThreeVector posiDirection = dp->GetMomentumDirection();
198 
199  G4double tau = posiKinEnergy/electron_mass_c2;
200  G4double gam = tau + 1.0;
201  G4double tau2 = tau + 2.0;
202  G4double sqgrate = sqrt(tau/tau2)*0.5;
203  G4double sqg2m1 = sqrt(tau*tau2);
204 
205  // limits of the energy sampling
206  G4double epsilmin = 0.5 - sqgrate;
207  G4double epsilmax = 0.5 + sqgrate;
208  G4double epsilqot = epsilmax/epsilmin;
209 
210  //
211  // sample the energy rate of the created gammas
212  //
213  G4double epsil, greject;
214 
215  do {
216  epsil = epsilmin*G4Exp(G4Log(epsilqot)*rndmEngine->flat());
217  greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
218  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
219  } while( greject < rndmEngine->flat());
220 
221  //
222  // scattered Gamma angles. ( Z - axis along the parent positron)
223  //
224 
225  G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
226  if(std::abs(cost) > 1.0) {
227  G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost
228  << " positron Ekin(MeV)= " << posiKinEnergy
229  << " gamma epsil= " << epsil
230  << G4endl;
231  if(cost > 1.0) cost = 1.0;
232  else cost = -1.0;
233  }
234  G4double sint = sqrt((1.+cost)*(1.-cost));
235  G4double phi = twopi * rndmEngine->flat();
236 
237  //
238  // kinematic of the created pair
239  //
240 
241  G4double totalEnergy = posiKinEnergy + 2.0*electron_mass_c2;
242  G4double phot1Energy = epsil*totalEnergy;
243 
244  G4ThreeVector phot1Direction(sint*cos(phi), sint*sin(phi), cost);
245  phot1Direction.rotateUz(posiDirection);
246  aGamma1 = new G4DynamicParticle (theGamma,phot1Direction, phot1Energy);
247  phi = twopi * rndmEngine->flat();
248  G4double cosphi = cos(phi);
249  G4double sinphi = sin(phi);
250  G4ThreeVector pol(cosphi, sinphi, 0.0);
251  pol.rotateUz(phot1Direction);
252  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
253 
254  G4double phot2Energy =(1.-epsil)*totalEnergy;
255  G4double posiP= sqrt(posiKinEnergy*(posiKinEnergy+2.*electron_mass_c2));
256  G4ThreeVector dir = posiDirection*posiP - phot1Direction*phot1Energy;
257  G4ThreeVector phot2Direction = dir.unit();
258 
259  // create G4DynamicParticle object for the particle2
260  aGamma2 = new G4DynamicParticle (theGamma, phot2Direction, phot2Energy);
261 
263  pol.set(-sinphi, cosphi, 0.0);
264  pol.rotateUz(phot1Direction);
265  cost = pol*phot2Direction;
266  pol -= cost*phot2Direction;
267  pol = pol.unit();
268  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
269  /*
270  G4cout << "Annihilation on fly: e0= " << posiKinEnergy
271  << " m= " << electron_mass_c2
272  << " e1= " << phot1Energy
273  << " e2= " << phot2Energy << " dir= " << dir
274  << " -> " << phot1Direction << " "
275  << phot2Direction << G4endl;
276  */
277  }
278 
279  vdp->push_back(aGamma1);
280  vdp->push_back(aGamma2);
281 
282  // kill primary positron
283  fParticleChange->SetProposedKineticEnergy(0.0);
284  fParticleChange->ProposeTrackStatus(fStopAndKill);
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....