| 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 // $Id: G4eeToTwoGammaModel.cc 101198 2016-11-09 09:34:52Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class file 30 // GEANT4 Class file
30 // 31 //
31 // 32 //
32 // File name: G4eeToTwoGammaModel 33 // File name: G4eeToTwoGammaModel
33 // 34 //
34 // Author: Vladimir Ivanchenko on base 35 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 // 36 //
36 // Creation date: 02.08.2004 37 // Creation date: 02.08.2004
37 // 38 //
38 // Modifications: 39 // Modifications:
39 // 08-04-05 Major optimisation of internal int 40 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
40 // 18-04-05 Compute CrossSectionPerVolume (V.I 41 // 18-04-05 Compute CrossSectionPerVolume (V.Ivanchenko)
41 // 06-02-06 ComputeCrossSectionPerElectron, Co 42 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
42 // 29-06-06 Fix problem for zero energy incide 43 // 29-06-06 Fix problem for zero energy incident positron (V.Ivanchenko)
43 // 20-10-06 Add theGamma as a member (V.Ivanch 44 // 20-10-06 Add theGamma as a member (V.Ivanchenko)
44 // 18-01-20 Introduce thermal model of annihil <<
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.
59 // is not defined <<
60 // 59 //
61 // Note 2: The annihilation processes producin 60 // Note 2: The annihilation processes producing one or more than two photons are
62 // ignored, as negligible compared to 61 // ignored, as negligible compared to the two photons process.
63 62
>> 63
>> 64
64 // 65 //
65 // ------------------------------------------- 66 // -------------------------------------------------------------------
66 // 67 //
67 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 70
70 #include "G4eeToTwoGammaModel.hh" 71 #include "G4eeToTwoGammaModel.hh"
71 #include "G4PhysicalConstants.hh" 72 #include "G4PhysicalConstants.hh"
72 #include "G4SystemOfUnits.hh" 73 #include "G4SystemOfUnits.hh"
73 #include "G4TrackStatus.hh" 74 #include "G4TrackStatus.hh"
74 #include "G4Electron.hh" 75 #include "G4Electron.hh"
75 #include "G4Positron.hh" 76 #include "G4Positron.hh"
76 #include "G4Gamma.hh" 77 #include "G4Gamma.hh"
77 #include "Randomize.hh" 78 #include "Randomize.hh"
78 #include "G4RandomDirection.hh" <<
79 #include "G4ParticleChangeForGamma.hh" 79 #include "G4ParticleChangeForGamma.hh"
80 #include "G4EmParameters.hh" <<
81 #include "G4Log.hh" 80 #include "G4Log.hh"
82 #include "G4Exp.hh" 81 #include "G4Exp.hh"
83 82
84 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 84
>> 85 using namespace std;
>> 86
86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const 87 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
87 const 88 const G4String& nam)
88 : G4VEmModel(nam), 89 : G4VEmModel(nam),
89 pi_rcl2(CLHEP::pi*CLHEP::classic_electr_ra << 90 pi_rcl2(pi*classic_electr_radius*classic_electr_radius),
>> 91 isInitialised(false)
90 { 92 {
91 theGamma = G4Gamma::Gamma(); 93 theGamma = G4Gamma::Gamma();
92 fParticleChange = nullptr; 94 fParticleChange = nullptr;
93 } 95 }
94 96
95 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 98
97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = << 99 G4eeToTwoGammaModel::~G4eeToTwoGammaModel()
>> 100 {}
98 101
99 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 103
101 void G4eeToTwoGammaModel::Initialise(const G4P 104 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
102 const G4D 105 const G4DataVector&)
103 { 106 {
104 if (nullptr != fParticleChange) { return; } << 107 if(isInitialised) { return; }
105 fParticleChange = GetParticleChangeForGamma( 108 fParticleChange = GetParticleChangeForGamma();
>> 109 isInitialised = true;
106 } 110 }
107 111
108 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 113
110 G4double << 114 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(
111 G4eeToTwoGammaModel::ComputeCrossSectionPerEle << 115 const G4ParticleDefinition*,
>> 116 G4double kineticEnergy,
>> 117 G4double, G4double)
112 { 118 {
113 // Calculates the cross section per electron 119 // Calculates the cross section per electron of annihilation into two photons
114 // from the Heilter formula. 120 // from the Heilter formula.
115 121
116 G4double ekin = std::max(CLHEP::eV, kinetic << 122 G4double ekin = std::max(eV,kineticEnergy);
117 123
118 G4double tau = ekin/CLHEP::electron_mass_c << 124 G4double tau = ekin/electron_mass_c2;
119 G4double gam = tau + 1.0; 125 G4double gam = tau + 1.0;
120 G4double gamma2= gam*gam; 126 G4double gamma2= gam*gam;
121 G4double bg2 = tau * (tau+2.0); 127 G4double bg2 = tau * (tau+2.0);
122 G4double bg = std::sqrt(bg2); << 128 G4double bg = sqrt(bg2);
123 129
124 G4double cross = pi_rcl2*((gamma2+4*gam+1.)* 130 G4double cross = pi_rcl2*((gamma2+4*gam+1.)*G4Log(gam+bg) - (gam+3.)*bg)
125 / (bg2*(gam+1.)); 131 / (bg2*(gam+1.));
126 return cross; 132 return cross;
127 } 133 }
128 134
129 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 136
131 G4double G4eeToTwoGammaModel::ComputeCrossSect 137 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
132 const G4Pa << 138 const G4ParticleDefinition* p,
133 G4double k 139 G4double kineticEnergy, G4double Z,
134 G4double, G4double, G4double) 140 G4double, G4double, G4double)
135 { 141 {
136 // Calculates the cross section per atom of 142 // Calculates the cross section per atom of annihilation into two photons
137 return Z*ComputeCrossSectionPerElectron(kine << 143
>> 144 G4double cross = Z*ComputeCrossSectionPerElectron(p,kineticEnergy);
>> 145 return cross;
138 } 146 }
139 147
140 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 149
142 G4double G4eeToTwoGammaModel::CrossSectionPerV 150 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
143 const G4Material* material, 151 const G4Material* material,
144 const G4ParticleDefinition*, << 152 const G4ParticleDefinition* p,
145 G4double kineticEnergy, 153 G4double kineticEnergy,
146 G4double, G4double) 154 G4double, G4double)
147 { 155 {
148 // Calculates the cross section per volume o 156 // Calculates the cross section per volume of annihilation into two photons
149 return material->GetElectronDensity()*Comput << 157
>> 158 G4double eDensity = material->GetElectronDensity();
>> 159 G4double cross = eDensity*ComputeCrossSectionPerElectron(p,kineticEnergy);
>> 160 return cross;
150 } 161 }
151 162
152 //....oooOO0OOooo........oooOO0OOooo........oo 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153 164
154 // Polarisation of gamma according to M.H.L.Pr << 165 // Polarisation of gamma according to M.H.L.Pryce and J.C.Ward,
155 // Nature 4065 (1947) 435. 166 // Nature 4065 (1947) 435.
156 167
157 void G4eeToTwoGammaModel::SampleSecondaries(st << 168 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
158 const G4MaterialCutsCouple*, 169 const G4MaterialCutsCouple*,
159 const G4DynamicParticle* dp, 170 const G4DynamicParticle* dp,
160 G4double, 171 G4double,
161 G4double) 172 G4double)
162 { 173 {
163 // kill primary positron << 174 G4double PositKinEnergy = dp->GetKineticEnergy();
164 fParticleChange->SetProposedKineticEnergy(0. << 175 G4DynamicParticle *aGamma1, *aGamma2;
165 fParticleChange->ProposeTrackStatus(fStopAnd <<
166 <<
167 // Case at rest not considered anymore insid <<
168 G4LorentzVector lv(dp->GetMomentum(), <<
169 dp->GetKineticEnergy() + 2*CLHEP::ele <<
170 G4double eGammaCMS = 0.5 * lv.mag(); <<
171 <<
172 G4ThreeVector dir1 = G4RandomDirection(); <<
173 G4double phi = CLHEP::twopi * G4UniformRand( <<
174 G4double cosphi = std::cos(phi); <<
175 G4double sinphi = std::sin(phi); <<
176 G4ThreeVector pol1(cosphi, sinphi, 0.0); <<
177 pol1.rotateUz(dir1); <<
178 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCM <<
179 <<
180 G4ThreeVector pol2(-sinphi, cosphi, 0.0); <<
181 pol2.rotateUz(dir1); <<
182 176
183 // transformation to lab system << 177 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
184 lv1.boost(lv.boostVector()); << 178
185 lv -= lv1; << 179 // Case at rest
186 << 180 if(PositKinEnergy == 0.0) {
187 //!!! boost of polarisation vector is not ye << 181 G4double cost = 2.*rndmEngine->flat()-1.;
188 << 182 G4double sint = sqrt((1. - cost)*(1. + cost));
189 // use constructors optimal for massless par << 183 G4double phi = twopi * rndmEngine->flat();
190 auto aGamma1 = new G4DynamicParticle(G4Gamma << 184 G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
191 aGamma1->SetPolarization(pol1); << 185 phi = twopi * rndmEngine->flat();
192 auto aGamma2 = new G4DynamicParticle(G4Gamma << 186 G4double cosphi = cos(phi);
193 aGamma2->SetPolarization(pol2); << 187 G4double sinphi = sin(phi);
>> 188 G4ThreeVector pol(cosphi, sinphi, 0.0);
>> 189 pol.rotateUz(dir);
>> 190 aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
>> 191 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
>> 192 aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
>> 193 pol.set(-sinphi, cosphi, 0.0);
>> 194 pol.rotateUz(dir);
>> 195 aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
>> 196
>> 197 } else {
>> 198
>> 199 G4ThreeVector PositDirection = dp->GetMomentumDirection();
>> 200
>> 201 G4double tau = PositKinEnergy/electron_mass_c2;
>> 202 G4double gam = tau + 1.0;
>> 203 G4double tau2 = tau + 2.0;
>> 204 G4double sqgrate = sqrt(tau/tau2)*0.5;
>> 205 G4double sqg2m1 = sqrt(tau*tau2);
>> 206
>> 207 // limits of the energy sampling
>> 208 G4double epsilmin = 0.5 - sqgrate;
>> 209 G4double epsilmax = 0.5 + sqgrate;
>> 210 G4double epsilqot = epsilmax/epsilmin;
>> 211
>> 212 //
>> 213 // sample the energy rate of the created gammas
>> 214 //
>> 215 G4double epsil, greject;
>> 216
>> 217 do {
>> 218 epsil = epsilmin*G4Exp(G4Log(epsilqot)*rndmEngine->flat());
>> 219 greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
>> 220 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
>> 221 } while( greject < rndmEngine->flat());
>> 222
>> 223 //
>> 224 // scattered Gamma angles. ( Z - axis along the parent positron)
>> 225 //
>> 226
>> 227 G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
>> 228 if(std::abs(cost) > 1.0) {
>> 229 G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost
>> 230 << " positron Ekin(MeV)= " << PositKinEnergy
>> 231 << " gamma epsil= " << epsil
>> 232 << G4endl;
>> 233 if(cost > 1.0) cost = 1.0;
>> 234 else cost = -1.0;
>> 235 }
>> 236 G4double sint = sqrt((1.+cost)*(1.-cost));
>> 237 G4double phi = twopi * rndmEngine->flat();
>> 238
>> 239 //
>> 240 // kinematic of the created pair
>> 241 //
>> 242
>> 243 G4double TotalAvailableEnergy = PositKinEnergy + 2.0*electron_mass_c2;
>> 244 G4double Phot1Energy = epsil*TotalAvailableEnergy;
>> 245
>> 246 G4ThreeVector Phot1Direction(sint*cos(phi), sint*sin(phi), cost);
>> 247 Phot1Direction.rotateUz(PositDirection);
>> 248 aGamma1 = new G4DynamicParticle (theGamma,Phot1Direction, Phot1Energy);
>> 249 phi = twopi * rndmEngine->flat();
>> 250 G4double cosphi = cos(phi);
>> 251 G4double sinphi = sin(phi);
>> 252 G4ThreeVector pol(cosphi, sinphi, 0.0);
>> 253 pol.rotateUz(Phot1Direction);
>> 254 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
>> 255
>> 256 G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
>> 257 G4double PositP= sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
>> 258 G4ThreeVector dir = PositDirection*PositP - Phot1Direction*Phot1Energy;
>> 259 G4ThreeVector Phot2Direction = dir.unit();
>> 260
>> 261 // create G4DynamicParticle object for the particle2
>> 262 aGamma2 = new G4DynamicParticle (theGamma,Phot2Direction, Phot2Energy);
>> 263
>> 264 //!!! likely problematic direction to be checked
>> 265 pol.set(-sinphi, cosphi, 0.0);
>> 266 pol.rotateUz(Phot1Direction);
>> 267 cost = pol*Phot2Direction;
>> 268 pol -= cost*Phot2Direction;
>> 269 pol = pol.unit();
>> 270 aGamma2->SetPolarization(pol.x(),pol.y(),pol.z());
>> 271 }
>> 272 /*
>> 273 G4cout << "Annihilation in fly: e0= " << PositKinEnergy
>> 274 << " m= " << electron_mass_c2
>> 275 << " e1= " << Phot1Energy
>> 276 << " e2= " << Phot2Energy << " dir= " << dir
>> 277 << " -> " << Phot1Direction << " "
>> 278 << Phot2Direction << G4endl;
>> 279 */
194 280
195 vdp->push_back(aGamma1); 281 vdp->push_back(aGamma1);
196 vdp->push_back(aGamma2); 282 vdp->push_back(aGamma2);
>> 283 fParticleChange->SetProposedKineticEnergy(0.);
>> 284 fParticleChange->ProposeTrackStatus(fStopAndKill);
197 } 285 }
198 286
199 //....oooOO0OOooo........oooOO0OOooo........oo 287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200 288