ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eplusTo2GammaOKVIModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eplusTo2GammaOKVIModel.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: G4eplusTo2GammaOKVIModel
33 //
34 // Author: Vladimir Ivanchenko and Omrame Kadri
35 //
36 // Creation date: 29.03.2018
37 //
38 // -------------------------------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 
44 #include "G4eplusTo2GammaOKVIModel.hh"
45 #include "G4eplusTo3GammaOKVIModel.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4EmParameters.hh"
49 #include "G4TrackStatus.hh"
50 #include "G4Electron.hh"
51 #include "G4Positron.hh"
52 #include "G4Gamma.hh"
53 #include "G4DataVector.hh"
54 #include "G4PhysicsVector.hh"
55 #include "G4PhysicsLogVector.hh"
56 #include "Randomize.hh"
57 #include "G4ParticleChangeForGamma.hh"
58 #include "G4Log.hh"
59 #include "G4Exp.hh"
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 
63 using namespace std;
64 
65 G4PhysicsVector* G4eplusTo2GammaOKVIModel::fCrossSection = nullptr;
66 G4PhysicsVector* G4eplusTo2GammaOKVIModel::fCrossSection3G = nullptr;
67 G4PhysicsVector* G4eplusTo2GammaOKVIModel::f3GProbability = nullptr;
68 
69 G4eplusTo2GammaOKVIModel::G4eplusTo2GammaOKVIModel(const G4ParticleDefinition*,
70  const G4String& nam)
71  : G4VEmModel(nam),
72  fDelta(0.001),
73  fGammaTh(MeV)
74 {
75  theGamma = G4Gamma::Gamma();
76  fParticleChange = nullptr;
77  fCuts = nullptr;
78  f3GModel = new G4eplusTo3GammaOKVIModel();
79  SetTripletModel(f3GModel);
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 G4eplusTo2GammaOKVIModel::~G4eplusTo2GammaOKVIModel()
85 {}
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
89 void G4eplusTo2GammaOKVIModel::Initialise(const G4ParticleDefinition* p,
90  const G4DataVector& cuts)
91 {
92  f3GModel->Initialise(p, cuts);
93  fCuts = &cuts;
94  fGammaTh = G4EmParameters::Instance()->LowestTripletEnergy();
95  f3GModel->SetDelta(fDelta);
96 
97  if(IsMaster()) {
98  if(!fCrossSection) {
99  G4double emin = 10*eV;
100  G4double emax = 100*TeV;
101  G4int nbins = 20*G4lrint(std::log10(emax/emin));
102  fCrossSection = new G4PhysicsLogVector(emin, emax, nbins);
103  fCrossSection3G = new G4PhysicsLogVector(emin, emax, nbins);
104  f3GProbability = new G4PhysicsLogVector(emin, emax, nbins);
105  fCrossSection->SetSpline(true);
106  fCrossSection3G->SetSpline(true);
107  f3GProbability->SetSpline(true);
108  for(G4int i=0; i<= nbins; ++i) {
109  G4double e = fCrossSection->Energy(i);
110  G4double cs2 = ComputeCrossSectionPerElectron(e);
111  G4double cs3 = f3GModel->ComputeCrossSectionPerElectron(e);
112  cs2 += cs3;
113  fCrossSection->PutValue(i, cs2);
114  fCrossSection3G->PutValue(i, cs3);
115  f3GProbability->PutValue(i, cs3/cs2);
116  }
117  }
118  }
119  // here particle change is set for the triplet model
120  if(fParticleChange) { return; }
121  fParticleChange = GetParticleChangeForGamma();
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 
126 G4double
127 G4eplusTo2GammaOKVIModel::ComputeCrossSectionPerElectron(G4double kinEnergy)
128 {
129  // Calculates the cross section per electron of annihilation into two
130  // photons from the Heilter formula with the radiation correction to 3 gamma
131  // annihilation channel. (A.A.) rho is changed
132 
133  G4double ekin = std::max(eV,kinEnergy);
134  G4double tau = ekin/electron_mass_c2;
135  G4double gam = tau + 1.0;
136  G4double gamma2 = gam*gam;
137  G4double bg2 = tau * (tau+2.0);
138  G4double bg = sqrt(bg2);
139  G4double rho = (gamma2+4.*gam+1.)*G4Log(gam+bg)/(gamma2-1.)
140  - (gam+3.)/(sqrt(gam*gam - 1.));
141 
142  static const G4double pir2 = pi*classic_electr_radius*classic_electr_radius;
143  G4double cross = (pir2*rho + alpha_rcl2*2.*G4Log(fDelta)*rho*rho)/(gam+1.);
144 
145  return cross;
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 G4double G4eplusTo2GammaOKVIModel::ComputeCrossSectionPerAtom(
151  const G4ParticleDefinition*,
152  G4double kineticEnergy, G4double Z,
153  G4double, G4double, G4double)
154 {
155  // Calculates the cross section per atom of annihilation into two photons
156  G4double cross = Z*fCrossSection->Value(kineticEnergy);
157  return cross;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 G4double G4eplusTo2GammaOKVIModel::CrossSectionPerVolume(
163  const G4Material* material,
164  const G4ParticleDefinition*,
165  G4double kineticEnergy,
166  G4double, G4double)
167 {
168  // Calculates the cross section per volume of annihilation into two photons
169  G4double eDensity = material->GetElectronDensity();
170  G4double cross = eDensity*fCrossSection->Value(kineticEnergy);
171  return cross;
172 }
173 
174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175 
176 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
177 // Nature 4065 (1947) 435.
178 
179 void
180 G4eplusTo2GammaOKVIModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
181  const G4MaterialCutsCouple* mcc,
182  const G4DynamicParticle* dp,
183  G4double, G4double)
184 {
185  G4double posiKinEnergy = dp->GetKineticEnergy();
186  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
187 
188  if(rndmEngine->flat() < f3GProbability->Value(posiKinEnergy)) {
189  G4double cutd = std::max(fGammaTh,(*fCuts)[mcc->GetIndex()])
190  /(posiKinEnergy + electron_mass_c2);
191  // check cut to avoid production of 3d gamma below
192  if(cutd > fDelta) {
193  G4double cs30 = fCrossSection3G->Value(posiKinEnergy);
194  f3GModel->SetDelta(cutd);
195  G4double cs3 = f3GModel->ComputeCrossSectionPerElectron(posiKinEnergy);
196  if(rndmEngine->flat()*cs30 < cs3) {
197  f3GModel->SampleSecondaries(vdp, mcc, dp);
198  return;
199  }
200  } else {
201  f3GModel->SampleSecondaries(vdp, mcc, dp);
202  return;
203  }
204  }
205 
206  G4DynamicParticle *aGamma1, *aGamma2;
207 
208  // Case at rest
209  if(posiKinEnergy == 0.0) {
210  G4double cost = 2.*rndmEngine->flat()-1.;
211  G4double sint = sqrt((1. - cost)*(1. + cost));
212  G4double phi = twopi * rndmEngine->flat();
213  G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
214  phi = twopi * rndmEngine->flat();
215  G4double cosphi = cos(phi);
216  G4double sinphi = sin(phi);
217  G4ThreeVector pol(cosphi, sinphi, 0.0);
218  pol.rotateUz(dir);
219  aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
220  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
221  aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
222  pol.set(-sinphi, cosphi, 0.0);
223  pol.rotateUz(dir);
224  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
225 
226  } else {
227 
228  G4ThreeVector posiDirection = dp->GetMomentumDirection();
229 
230  G4double tau = posiKinEnergy/electron_mass_c2;
231  G4double gam = tau + 1.0;
232  G4double tau2 = tau + 2.0;
233  G4double sqgrate = sqrt(tau/tau2)*0.5;
234  G4double sqg2m1 = sqrt(tau*tau2);
235 
236  // limits of the energy sampling
237  G4double epsilmin = 0.5 - sqgrate;
238  G4double epsilmax = 0.5 + sqgrate;
239  G4double epsilqot = epsilmax/epsilmin;
240 
241  //
242  // sample the energy rate of the created gammas
243  //
244  G4double epsil, greject;
245 
246  do {
247  epsil = epsilmin*G4Exp(G4Log(epsilqot)*rndmEngine->flat());
248  greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
249  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
250  } while( greject < rndmEngine->flat());
251 
252  //
253  // scattered Gamma angles. ( Z - axis along the parent positron)
254  //
255 
256  G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
257  if(std::abs(cost) > 1.0) {
258  G4cout << "### G4eplusTo2GammaOKVIModel WARNING cost= " << cost
259  << " positron Ekin(MeV)= " << posiKinEnergy
260  << " gamma epsil= " << epsil
261  << G4endl;
262  if(cost > 1.0) cost = 1.0;
263  else cost = -1.0;
264  }
265  G4double sint = sqrt((1.+cost)*(1.-cost));
266  G4double phi = twopi * rndmEngine->flat();
267 
268  //
269  // kinematic of the created pair
270  //
271 
272  G4double TotalAvailableEnergy = posiKinEnergy + 2.0*electron_mass_c2;
273  G4double phot1Energy = epsil*TotalAvailableEnergy;
274 
275  G4ThreeVector phot1Direction(sint*cos(phi), sint*sin(phi), cost);
276  phot1Direction.rotateUz(posiDirection);
277  aGamma1 = new G4DynamicParticle (theGamma,phot1Direction, phot1Energy);
278  phi = twopi * rndmEngine->flat();
279  G4double cosphi = cos(phi);
280  G4double sinphi = sin(phi);
281  G4ThreeVector pol(cosphi, sinphi, 0.0);
282  pol.rotateUz(phot1Direction);
283  aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
284 
285  G4double phot2Energy =(1.-epsil)*TotalAvailableEnergy;
286  G4double posiP= sqrt(posiKinEnergy*(posiKinEnergy+2.*electron_mass_c2));
287  G4ThreeVector dir = posiDirection*posiP - phot1Direction*phot1Energy;
288  G4ThreeVector phot2Direction = dir.unit();
289 
290  // create G4DynamicParticle object for the particle2
291  aGamma2 = new G4DynamicParticle (theGamma,phot2Direction, phot2Energy);
292 
294  pol.set(-sinphi, cosphi, 0.0);
295  pol.rotateUz(phot1Direction);
296  cost = pol*phot2Direction;
297  pol -= cost*phot2Direction;
298  pol = pol.unit();
299  aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
300  }
301  /*
302  G4cout << "Annihilation in fly: e0= " << posiKinEnergy
303  << " m= " << electron_mass_c2
304  << " e1= " << phot1Energy
305  << " e2= " << phot2Energy << " dir= " << dir
306  << " -> " << phot1Direction << " "
307  << phot2Direction << G4endl;
308  */
309 
310  vdp->push_back(aGamma1);
311  vdp->push_back(aGamma2);
312 
313  // kill primary positron
314  fParticleChange->SetProposedKineticEnergy(0.0);
315  fParticleChange->ProposeTrackStatus(fStopAndKill);
316 }
317 
318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....