| 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: G4eeToTwoGammaModel 32 // File name: G4eeToTwoGammaModel
33 // 33 //
34 // Author: Vladimir Ivanchenko on base 34 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 // 35 //
36 // Creation date: 02.08.2004 36 // Creation date: 02.08.2004
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 08-04-05 Major optimisation of internal int 39 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
40 // 18-04-05 Compute CrossSectionPerVolume (V.I 40 // 18-04-05 Compute CrossSectionPerVolume (V.Ivanchenko)
41 // 06-02-06 ComputeCrossSectionPerElectron, Co 41 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
42 // 29-06-06 Fix problem for zero energy incide 42 // 29-06-06 Fix problem for zero energy incident positron (V.Ivanchenko)
43 // 20-10-06 Add theGamma as a member (V.Ivanch 43 // 20-10-06 Add theGamma as a member (V.Ivanchenko)
44 // 18-01-20 Introduce thermal model of annihil 44 // 18-01-20 Introduce thermal model of annihilation at rest (J.Allison)
45 // 45 //
46 // 46 //
47 // Class Description: 47 // Class Description:
48 // 48 //
49 // Implementation of e+ annihilation into 2 ga 49 // Implementation of e+ annihilation into 2 gamma
50 // 50 //
51 // The secondaries Gamma energies are sampled 51 // The secondaries Gamma energies are sampled using the Heitler cross section.
52 // 52 //
53 // A modified version of the random number tec 53 // A modified version of the random number techniques of Butcher & Messel
54 // is used (Nuc Phys 20(1960),15). 54 // is used (Nuc Phys 20(1960),15).
55 // 55 //
56 // GEANT4 internal units. 56 // GEANT4 internal units.
57 // 57 //
58 // Note 1: The initial electron is assumed fre 58 // Note 1: The initial electron is assumed free and at rest if atomic PDF
59 // is not defined 59 // is not defined
60 // 60 //
61 // Note 2: The annihilation processes producin 61 // Note 2: The annihilation processes producing one or more than two photons are
62 // ignored, as negligible compared to 62 // ignored, as negligible compared to the two photons process.
63 63
64 // 64 //
65 // ------------------------------------------- 65 // -------------------------------------------------------------------
66 // 66 //
67 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 69
70 #include "G4eeToTwoGammaModel.hh" 70 #include "G4eeToTwoGammaModel.hh"
71 #include "G4PhysicalConstants.hh" 71 #include "G4PhysicalConstants.hh"
72 #include "G4SystemOfUnits.hh" 72 #include "G4SystemOfUnits.hh"
73 #include "G4TrackStatus.hh" 73 #include "G4TrackStatus.hh"
74 #include "G4Electron.hh" 74 #include "G4Electron.hh"
75 #include "G4Positron.hh" 75 #include "G4Positron.hh"
76 #include "G4Gamma.hh" 76 #include "G4Gamma.hh"
77 #include "Randomize.hh" 77 #include "Randomize.hh"
78 #include "G4RandomDirection.hh" 78 #include "G4RandomDirection.hh"
79 #include "G4ParticleChangeForGamma.hh" 79 #include "G4ParticleChangeForGamma.hh"
80 #include "G4EmParameters.hh" 80 #include "G4EmParameters.hh"
81 #include "G4Log.hh" 81 #include "G4Log.hh"
82 #include "G4Exp.hh" 82 #include "G4Exp.hh"
83 83
84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 85
>> 86 using namespace std;
>> 87
>> 88 G4bool G4eeToTwoGammaModel::fSampleAtomicPDF = false;
>> 89
86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const 90 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
87 const 91 const G4String& nam)
88 : G4VEmModel(nam), 92 : G4VEmModel(nam),
89 pi_rcl2(CLHEP::pi*CLHEP::classic_electr_ra << 93 pi_rcl2(pi*classic_electr_radius*classic_electr_radius)
90 { 94 {
91 theGamma = G4Gamma::Gamma(); 95 theGamma = G4Gamma::Gamma();
92 fParticleChange = nullptr; 96 fParticleChange = nullptr;
93 } 97 }
94 98
95 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 100
97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = 101 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = default;
98 102
99 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 104
101 void G4eeToTwoGammaModel::Initialise(const G4P 105 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
102 const G4D 106 const G4DataVector&)
103 { 107 {
104 if (nullptr != fParticleChange) { return; } << 108 if(IsMaster()) {
>> 109 G4int verbose = G4EmParameters::Instance()->Verbose();
>> 110 // redo initialisation for each new run
>> 111 fSampleAtomicPDF = false;
>> 112 const auto& materialTable = G4Material::GetMaterialTable();
>> 113 for (const auto& material: *materialTable) {
>> 114 const G4double meanEnergyPerIonPair = material->GetIonisation()->GetMeanEnergyPerIonPair();
>> 115 if (meanEnergyPerIonPair > 0.) {
>> 116 fSampleAtomicPDF = true;
>> 117 if(verbose > 0) {
>> 118 G4cout << "### G4eeToTwoGammaModel: for " << material->GetName() << " mean energy per ion pair is "
>> 119 << meanEnergyPerIonPair/CLHEP::eV << " eV" << G4endl;
>> 120 }
>> 121 }
>> 122 }
>> 123 }
>> 124 // If no materials have meanEnergyPerIonPair set. This is probably the usual
>> 125 // case, since most applications are not senstive to the slight
>> 126 // non-collinearity of gammas in eeToTwoGamma. Do not issue any warning.
>> 127
>> 128 if(fParticleChange) { return; }
105 fParticleChange = GetParticleChangeForGamma( 129 fParticleChange = GetParticleChangeForGamma();
106 } 130 }
107 131
108 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 133
110 G4double 134 G4double
111 G4eeToTwoGammaModel::ComputeCrossSectionPerEle 135 G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(G4double kineticEnergy)
112 { 136 {
113 // Calculates the cross section per electron 137 // Calculates the cross section per electron of annihilation into two photons
114 // from the Heilter formula. 138 // from the Heilter formula.
115 139
116 G4double ekin = std::max(CLHEP::eV, kinetic << 140 G4double ekin = std::max(eV,kineticEnergy);
117 141
118 G4double tau = ekin/CLHEP::electron_mass_c << 142 G4double tau = ekin/electron_mass_c2;
119 G4double gam = tau + 1.0; 143 G4double gam = tau + 1.0;
120 G4double gamma2= gam*gam; 144 G4double gamma2= gam*gam;
121 G4double bg2 = tau * (tau+2.0); 145 G4double bg2 = tau * (tau+2.0);
122 G4double bg = std::sqrt(bg2); << 146 G4double bg = sqrt(bg2);
123 147
124 G4double cross = pi_rcl2*((gamma2+4*gam+1.)* 148 G4double cross = pi_rcl2*((gamma2+4*gam+1.)*G4Log(gam+bg) - (gam+3.)*bg)
125 / (bg2*(gam+1.)); 149 / (bg2*(gam+1.));
126 return cross; 150 return cross;
127 } 151 }
128 152
129 //....oooOO0OOooo........oooOO0OOooo........oo 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 154
131 G4double G4eeToTwoGammaModel::ComputeCrossSect 155 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
132 const G4Pa 156 const G4ParticleDefinition*,
133 G4double k 157 G4double kineticEnergy, G4double Z,
134 G4double, G4double, G4double) 158 G4double, G4double, G4double)
135 { 159 {
136 // Calculates the cross section per atom of 160 // Calculates the cross section per atom of annihilation into two photons
137 return Z*ComputeCrossSectionPerElectron(kine 161 return Z*ComputeCrossSectionPerElectron(kineticEnergy);
138 } 162 }
139 163
140 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 165
142 G4double G4eeToTwoGammaModel::CrossSectionPerV 166 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
143 const G4Material* material, 167 const G4Material* material,
144 const G4ParticleDefinition*, 168 const G4ParticleDefinition*,
145 G4double kineticEnergy, 169 G4double kineticEnergy,
146 G4double, G4double) 170 G4double, G4double)
147 { 171 {
148 // Calculates the cross section per volume o 172 // Calculates the cross section per volume of annihilation into two photons
149 return material->GetElectronDensity()*Comput 173 return material->GetElectronDensity()*ComputeCrossSectionPerElectron(kineticEnergy);
150 } 174 }
151 175
152 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153 177
154 // Polarisation of gamma according to M.H.L.Pr 178 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
155 // Nature 4065 (1947) 435. 179 // Nature 4065 (1947) 435.
156 180
157 void G4eeToTwoGammaModel::SampleSecondaries(st << 181 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
158 const G4MaterialCutsCouple*, << 182 const G4MaterialCutsCouple* pCutsCouple,
159 const G4DynamicParticle* dp, 183 const G4DynamicParticle* dp,
160 G4double, 184 G4double,
161 G4double) 185 G4double)
162 { 186 {
163 // kill primary positron << 187 G4double posiKinEnergy = dp->GetKineticEnergy();
164 fParticleChange->SetProposedKineticEnergy(0. << 188 G4DynamicParticle *aGamma1, *aGamma2;
165 fParticleChange->ProposeTrackStatus(fStopAnd <<
166 189
167 // Case at rest not considered anymore insid << 190 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
168 G4LorentzVector lv(dp->GetMomentum(), << 191
169 dp->GetKineticEnergy() + 2*CLHEP::ele << 192 // Case at rest
170 G4double eGammaCMS = 0.5 * lv.mag(); << 193 if(posiKinEnergy == 0.0) {
171 << 194
172 G4ThreeVector dir1 = G4RandomDirection(); << 195 const G4double eGamma = electron_mass_c2;
173 G4double phi = CLHEP::twopi * G4UniformRand( << 196
174 G4double cosphi = std::cos(phi); << 197 // In rest frame of positronium gammas are back to back
175 G4double sinphi = std::sin(phi); << 198 const G4ThreeVector& dir1 = G4RandomDirection();
176 G4ThreeVector pol1(cosphi, sinphi, 0.0); << 199 const G4ThreeVector& dir2 = -dir1;
177 pol1.rotateUz(dir1); << 200 aGamma1 = new G4DynamicParticle(G4Gamma::Gamma(),dir1,eGamma);
178 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCM << 201 aGamma2 = new G4DynamicParticle(G4Gamma::Gamma(),dir2,eGamma);
179 << 202
180 G4ThreeVector pol2(-sinphi, cosphi, 0.0); << 203 // In rest frame the gammas are polarised perpendicular to each other - see
181 pol2.rotateUz(dir1); << 204 // Pryce and Ward, Nature No 4065 (1947) p.435.
182 << 205 // Snyder et al, Physical Review 73 (1948) p.440.
183 // transformation to lab system << 206 G4ThreeVector pol1 = (G4RandomDirection().cross(dir1)).unit();
184 lv1.boost(lv.boostVector()); << 207 G4ThreeVector pol2 = (pol1.cross(dir2)).unit();
185 lv -= lv1; << 208
186 << 209 // But the positronium is moving...
187 //!!! boost of polarisation vector is not ye << 210 // A positron in matter slows down and combines with an atomic electron to
188 << 211 // make a neutral “atom” called positronium, about half the size of a normal
189 // use constructors optimal for massless par << 212 // atom. I expect that when the energy of the positron is small enough,
190 auto aGamma1 = new G4DynamicParticle(G4Gamma << 213 // less than the binding energy of positronium (6.8 eV), it is
191 aGamma1->SetPolarization(pol1); << 214 // energetically favourable for an electron from the outer orbitals of a
192 auto aGamma2 = new G4DynamicParticle(G4Gamma << 215 // nearby atom or molecule to transfer and bind to the positron, as in an
193 aGamma2->SetPolarization(pol2); << 216 // ionic bond, leaving behind a mildly ionised nearby atom/molecule. I
>> 217 // would expect the positronium to come away with a kinetic energy of a
>> 218 // few eV on average. In its para (spin 0) state it annihilates into two
>> 219 // photons, which in the rest frame of the positronium are collinear
>> 220 // (back-to-back) due to momentum conservation. Because of the motion of the
>> 221 // positronium, photons will be not quite back-to-back in the laboratory.
>> 222
>> 223 // The positroniuim acquires an energy of order its binding energy and
>> 224 // doesn't have time to thermalise. Nevertheless, here we approximate its
>> 225 // energy distribution by a Maxwell-Boltzman with mean energy <KE>. In terms
>> 226 // of a more familiar concept of temperature, and the law of equipartition
>> 227 // of energy of translational motion, <KE>=3kT/2. Each component of velocity
>> 228 // has a distribution exp(-mv^2/2kT), which is a Gaussian of mean zero
>> 229 // and variance kT/m=2<KE>/3m, where m is the positronium mass.
>> 230
>> 231 // We take <KE> = material->GetIonisation()->GetMeanEnergyPerIonPair().
>> 232
>> 233 if(fSampleAtomicPDF) {
>> 234 const G4Material* material = pCutsCouple->GetMaterial();
>> 235 const G4double meanEnergyPerIonPair = material->GetIonisation()->GetMeanEnergyPerIonPair();
>> 236 const G4double& meanKE = meanEnergyPerIonPair; // Just an alias
>> 237 if (meanKE > 0.) { // Positronium haas motion
>> 238 // Mass of positronium
>> 239 const G4double mass = 2.*electron_mass_c2;
>> 240 // Mean <KE>=3kT/2, as described above
>> 241 // const G4double T = 2.*meanKE/(3.*k_Boltzmann);
>> 242 // Component velocities: Gaussian, variance kT/m=2<KE>/3m.
>> 243 const G4double sigmav = std::sqrt(2.*meanKE/(3.*mass));
>> 244 // This is in units where c=1
>> 245 const G4double vx = G4RandGauss::shoot(0.,sigmav);
>> 246 const G4double vy = G4RandGauss::shoot(0.,sigmav);
>> 247 const G4double vz = G4RandGauss::shoot(0.,sigmav);
>> 248 const G4ThreeVector v(vx,vy,vz); // In unit where c=1
>> 249 const G4ThreeVector& beta = v; // so beta=v/c=v
>> 250
>> 251 aGamma1->Set4Momentum(aGamma1->Get4Momentum().boost(beta));
>> 252 aGamma2->Set4Momentum(aGamma2->Get4Momentum().boost(beta));
>> 253
>> 254 // Rotate polarisation vectors
>> 255 const G4ThreeVector& newDir1 = aGamma1->GetMomentumDirection();
>> 256 const G4ThreeVector& newDir2 = aGamma2->GetMomentumDirection();
>> 257 const G4ThreeVector& axis1 = dir1.cross(newDir1); // No need to be unit
>> 258 const G4ThreeVector& axis2 = dir2.cross(newDir2); // No need to be unit
>> 259 const G4double& angle1 = std::acos(dir1*newDir1);
>> 260 const G4double& angle2 = std::acos(dir2*newDir2);
>> 261 if (axis1 != G4ThreeVector()) pol1.rotate(axis1,angle1);
>> 262 if (axis2 != G4ThreeVector()) pol2.rotate(axis2,angle2);
>> 263 }
>> 264 }
>> 265 aGamma1->SetPolarization(pol1.x(),pol1.y(),pol1.z());
>> 266 aGamma2->SetPolarization(pol2.x(),pol2.y(),pol2.z());
>> 267
>> 268 } else { // Positron interacts in flight
>> 269
>> 270 G4ThreeVector posiDirection = dp->GetMomentumDirection();
>> 271
>> 272 G4double tau = posiKinEnergy/electron_mass_c2;
>> 273 G4double gam = tau + 1.0;
>> 274 G4double tau2 = tau + 2.0;
>> 275 G4double sqgrate = sqrt(tau/tau2)*0.5;
>> 276 G4double sqg2m1 = sqrt(tau*tau2);
>> 277
>> 278 // limits of the energy sampling
>> 279 G4double epsilmin = 0.5 - sqgrate;
>> 280 G4double epsilmax = 0.5 + sqgrate;
>> 281 G4double epsilqot = epsilmax/epsilmin;
>> 282
>> 283 //
>> 284 // sample the energy rate of the created gammas
>> 285 //
>> 286 G4double epsil, greject;
>> 287
>> 288 do {
>> 289 epsil = epsilmin*G4Exp(G4Log(epsilqot)*rndmEngine->flat());
>> 290 greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
>> 291 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
>> 292 } while( greject < rndmEngine->flat());
>> 293
>> 294 //
>> 295 // scattered Gamma angles. ( Z - axis along the parent positron)
>> 296 //
>> 297
>> 298 G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
>> 299 if(std::abs(cost) > 1.0) {
>> 300 G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost
>> 301 << " positron Ekin(MeV)= " << posiKinEnergy
>> 302 << " gamma epsil= " << epsil
>> 303 << G4endl;
>> 304 if(cost > 1.0) cost = 1.0;
>> 305 else cost = -1.0;
>> 306 }
>> 307 G4double sint = sqrt((1.+cost)*(1.-cost));
>> 308 G4double phi = twopi * rndmEngine->flat();
>> 309
>> 310 //
>> 311 // kinematic of the created pair
>> 312 //
>> 313
>> 314 G4double totalEnergy = posiKinEnergy + 2.0*electron_mass_c2;
>> 315 G4double phot1Energy = epsil*totalEnergy;
>> 316
>> 317 G4ThreeVector phot1Direction(sint*cos(phi), sint*sin(phi), cost);
>> 318 phot1Direction.rotateUz(posiDirection);
>> 319 aGamma1 = new G4DynamicParticle (theGamma,phot1Direction, phot1Energy);
>> 320 phi = twopi * rndmEngine->flat();
>> 321 G4double cosphi = cos(phi);
>> 322 G4double sinphi = sin(phi);
>> 323 G4ThreeVector pol(cosphi, sinphi, 0.0);
>> 324 pol.rotateUz(phot1Direction);
>> 325 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
>> 326
>> 327 G4double phot2Energy =(1.-epsil)*totalEnergy;
>> 328 G4double posiP= sqrt(posiKinEnergy*(posiKinEnergy+2.*electron_mass_c2));
>> 329 G4ThreeVector dir = posiDirection*posiP - phot1Direction*phot1Energy;
>> 330 G4ThreeVector phot2Direction = dir.unit();
>> 331
>> 332 // create G4DynamicParticle object for the particle2
>> 333 aGamma2 = new G4DynamicParticle (theGamma, phot2Direction, phot2Energy);
>> 334
>> 335 //!!! likely problematic direction to be checked
>> 336 pol.set(-sinphi, cosphi, 0.0);
>> 337 pol.rotateUz(phot1Direction);
>> 338 cost = pol*phot2Direction;
>> 339 pol -= cost*phot2Direction;
>> 340 pol = pol.unit();
>> 341 aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
>> 342 /*
>> 343 G4cout << "Annihilation on fly: e0= " << posiKinEnergy
>> 344 << " m= " << electron_mass_c2
>> 345 << " e1= " << phot1Energy
>> 346 << " e2= " << phot2Energy << " dir= " << dir
>> 347 << " -> " << phot1Direction << " "
>> 348 << phot2Direction << G4endl;
>> 349 */
>> 350 }
194 351
195 vdp->push_back(aGamma1); 352 vdp->push_back(aGamma1);
196 vdp->push_back(aGamma2); 353 vdp->push_back(aGamma2);
>> 354
>> 355 // kill primary positron
>> 356 fParticleChange->SetProposedKineticEnergy(0.0);
>> 357 fParticleChange->ProposeTrackStatus(fStopAndKill);
197 } 358 }
198 359
199 //....oooOO0OOooo........oooOO0OOooo........oo 360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200 361