| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 // $Id: G4eeToTwoGammaModel.cc,v 1.8 2005/04/29 18:02:35 vnivanch Exp $
>> 24 // GEANT4 tag $Name: geant4-08-00-patch-01 $
26 // 25 //
27 // ------------------------------------------- 26 // -------------------------------------------------------------------
28 // 27 //
29 // GEANT4 Class file 28 // GEANT4 Class file
30 // 29 //
31 // 30 //
32 // File name: G4eeToTwoGammaModel 31 // File name: G4eeToTwoGammaModel
33 // 32 //
34 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 // 34 //
36 // Creation date: 02.08.2004 35 // Creation date: 02.08.2004
37 // 36 //
38 // Modifications: 37 // Modifications:
39 // 08-04-05 Major optimisation of internal int << 38 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
40 // 18-04-05 Compute CrossSectionPerVolume (V.I << 39 // 18-04-05 Compute CrossSectionPerVolume (V.Ivantchenko)
41 // 06-02-06 ComputeCrossSectionPerElectron, Co <<
42 // 29-06-06 Fix problem for zero energy incide <<
43 // 20-10-06 Add theGamma as a member (V.Ivanch <<
44 // 18-01-20 Introduce thermal model of annihil <<
45 // 40 //
46 // 41 //
47 // Class Description: 42 // Class Description:
48 // 43 //
49 // Implementation of e+ annihilation into 2 ga 44 // Implementation of e+ annihilation into 2 gamma
50 // 45 //
51 // The secondaries Gamma energies are sampled 46 // The secondaries Gamma energies are sampled using the Heitler cross section.
52 // 47 //
53 // A modified version of the random number tec 48 // A modified version of the random number techniques of Butcher & Messel
54 // is used (Nuc Phys 20(1960),15). 49 // is used (Nuc Phys 20(1960),15).
55 // 50 //
56 // GEANT4 internal units. 51 // GEANT4 internal units.
57 // 52 //
58 // Note 1: The initial electron is assumed fre << 53 // Note 1: The initial electron is assumed free and at rest.
59 // is not defined <<
60 // 54 //
61 // Note 2: The annihilation processes producin 55 // Note 2: The annihilation processes producing one or more than two photons are
62 // ignored, as negligible compared to 56 // ignored, as negligible compared to the two photons process.
63 57
>> 58
>> 59
64 // 60 //
65 // ------------------------------------------- 61 // -------------------------------------------------------------------
66 // 62 //
67 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 65
70 #include "G4eeToTwoGammaModel.hh" 66 #include "G4eeToTwoGammaModel.hh"
71 #include "G4PhysicalConstants.hh" <<
72 #include "G4SystemOfUnits.hh" <<
73 #include "G4TrackStatus.hh" <<
74 #include "G4Electron.hh" 67 #include "G4Electron.hh"
75 #include "G4Positron.hh" 68 #include "G4Positron.hh"
76 #include "G4Gamma.hh" 69 #include "G4Gamma.hh"
77 #include "Randomize.hh" 70 #include "Randomize.hh"
78 #include "G4RandomDirection.hh" <<
79 #include "G4ParticleChangeForGamma.hh" <<
80 #include "G4EmParameters.hh" <<
81 #include "G4Log.hh" <<
82 #include "G4Exp.hh" <<
83 71
84 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 73
>> 74 using namespace std;
>> 75
86 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const 76 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
87 const 77 const G4String& nam)
88 : G4VEmModel(nam), 78 : G4VEmModel(nam),
89 pi_rcl2(CLHEP::pi*CLHEP::classic_electr_ra << 79 pi_rcl2(pi*classic_electr_radius*classic_electr_radius)
90 { 80 {
91 theGamma = G4Gamma::Gamma(); <<
92 fParticleChange = nullptr; <<
93 } 81 }
94 82
95 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 84
97 G4eeToTwoGammaModel::~G4eeToTwoGammaModel() = << 85 G4eeToTwoGammaModel::~G4eeToTwoGammaModel()
>> 86 {}
98 87
99 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 89
101 void G4eeToTwoGammaModel::Initialise(const G4P 90 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
102 const G4D 91 const G4DataVector&)
103 { << 92 {}
104 if (nullptr != fParticleChange) { return; } <<
105 fParticleChange = GetParticleChangeForGamma( <<
106 } <<
107 93
108 //....oooOO0OOooo........oooOO0OOooo........oo << 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 95
110 G4double << 96 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(const G4Material* material,
111 G4eeToTwoGammaModel::ComputeCrossSectionPerEle << 97 const G4ParticleDefinition*,
>> 98 G4double kineticEnergy,
>> 99 G4double,
>> 100 G4double)
112 { 101 {
113 // Calculates the cross section per electron << 102 // Calculates the cross section per atom of annihilation into two photons
114 // from the Heilter formula. 103 // from the Heilter formula.
>> 104 G4double eDensity = material->GetElectronDensity();
115 105
116 G4double ekin = std::max(CLHEP::eV, kinetic << 106 G4double tau = kineticEnergy/electron_mass_c2;
117 <<
118 G4double tau = ekin/CLHEP::electron_mass_c <<
119 G4double gam = tau + 1.0; 107 G4double gam = tau + 1.0;
120 G4double gamma2= gam*gam; 108 G4double gamma2= gam*gam;
121 G4double bg2 = tau * (tau+2.0); 109 G4double bg2 = tau * (tau+2.0);
122 G4double bg = std::sqrt(bg2); << 110 G4double bg = sqrt(bg2);
123 111
124 G4double cross = pi_rcl2*((gamma2+4*gam+1.)* << 112 G4double cross = pi_rcl2*eDensity*((gamma2+4*gam+1.)*log(gam+bg) - (gam+3.)*bg)
125 / (bg2*(gam+1.)); 113 / (bg2*(gam+1.));
126 return cross; << 114 return cross;
127 } <<
128 <<
129 //....oooOO0OOooo........oooOO0OOooo........oo <<
130 <<
131 G4double G4eeToTwoGammaModel::ComputeCrossSect <<
132 const G4Pa <<
133 G4double k <<
134 G4double, G4double, G4double) <<
135 { <<
136 // Calculates the cross section per atom of <<
137 return Z*ComputeCrossSectionPerElectron(kine <<
138 } <<
139 <<
140 //....oooOO0OOooo........oooOO0OOooo........oo <<
141 <<
142 G4double G4eeToTwoGammaModel::CrossSectionPerV <<
143 const G4Material* material, <<
144 const G4ParticleDefinition*, <<
145 G4double kineticEnergy, <<
146 G4double, G4double) <<
147 { <<
148 // Calculates the cross section per volume o <<
149 return material->GetElectronDensity()*Comput <<
150 } 115 }
151 116
152 //....oooOO0OOooo........oooOO0OOooo........oo << 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 <<
154 // Polarisation of gamma according to M.H.L.Pr <<
155 // Nature 4065 (1947) 435. <<
156 118
157 void G4eeToTwoGammaModel::SampleSecondaries(st << 119 vector<G4DynamicParticle*>* G4eeToTwoGammaModel::SampleSecondaries(
158 const G4MaterialCutsCouple*, << 120 const G4MaterialCutsCouple*,
159 const G4DynamicParticle* dp, << 121 const G4DynamicParticle* dp,
160 G4double, << 122 G4double,
161 G4double) << 123 G4double)
162 { << 124 {
163 // kill primary positron << 125 G4double PositKinEnergy = dp->GetKineticEnergy();
164 fParticleChange->SetProposedKineticEnergy(0. << 126 G4ThreeVector PositDirection = dp->GetMomentumDirection();
165 fParticleChange->ProposeTrackStatus(fStopAnd << 127
166 << 128 G4double tau = PositKinEnergy/electron_mass_c2;
167 // Case at rest not considered anymore insid << 129 G4double gam = tau + 1.0;
168 G4LorentzVector lv(dp->GetMomentum(), << 130 G4double tau2 = tau + 2.0;
169 dp->GetKineticEnergy() + 2*CLHEP::ele << 131 G4double sqgrate = sqrt(tau/tau2)*0.5;
170 G4double eGammaCMS = 0.5 * lv.mag(); << 132 G4double sqg2m1 = sqrt(tau*tau2);
171 << 133
172 G4ThreeVector dir1 = G4RandomDirection(); << 134 // limits of the energy sampling
173 G4double phi = CLHEP::twopi * G4UniformRand( << 135 G4double epsilmin = 0.5 - sqgrate;
174 G4double cosphi = std::cos(phi); << 136 G4double epsilmax = 0.5 + sqgrate;
175 G4double sinphi = std::sin(phi); << 137 G4double epsilqot = epsilmax/epsilmin;
176 G4ThreeVector pol1(cosphi, sinphi, 0.0); << 138
177 pol1.rotateUz(dir1); << 139 //
178 G4LorentzVector lv1(eGammaCMS*dir1, eGammaCM << 140 // sample the energy rate of the created gammas
179 << 141 //
180 G4ThreeVector pol2(-sinphi, cosphi, 0.0); << 142 G4double epsil, greject;
181 pol2.rotateUz(dir1); << 143
182 << 144 do {
183 // transformation to lab system << 145 epsil = epsilmin*pow(epsilqot,G4UniformRand());
184 lv1.boost(lv.boostVector()); << 146 greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
185 lv -= lv1; << 147 } while( greject < G4UniformRand() );
186 << 148
187 //!!! boost of polarisation vector is not ye << 149 //
188 << 150 // scattered Gamma angles. ( Z - axis along the parent positron)
189 // use constructors optimal for massless par << 151 //
190 auto aGamma1 = new G4DynamicParticle(G4Gamma << 152
191 aGamma1->SetPolarization(pol1); << 153 G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
192 auto aGamma2 = new G4DynamicParticle(G4Gamma << 154 G4double sint = sqrt((1.+cost)*(1.-cost));
193 aGamma2->SetPolarization(pol2); << 155 G4double phi = twopi * G4UniformRand();
194 << 156
195 vdp->push_back(aGamma1); << 157 G4double dirx = sint*cos(phi) , diry = sint*sin(phi) , dirz = cost;
196 vdp->push_back(aGamma2); << 158
>> 159 //
>> 160 // kinematic of the created pair
>> 161 //
>> 162
>> 163 G4double TotalAvailableEnergy = PositKinEnergy + 2.0*electron_mass_c2;
>> 164 G4double Phot1Energy = epsil*TotalAvailableEnergy;
>> 165
>> 166 vector<G4DynamicParticle*>* vdp = new vector<G4DynamicParticle*>;
>> 167
>> 168 G4ThreeVector Phot1Direction (dirx, diry, dirz);
>> 169 Phot1Direction.rotateUz(PositDirection);
>> 170 G4DynamicParticle* aParticle1 = new G4DynamicParticle (G4Gamma::Gamma(),
>> 171 Phot1Direction, Phot1Energy);
>> 172 vdp->push_back(aParticle1);
>> 173
>> 174 G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
>> 175 G4double Eratio= Phot1Energy/Phot2Energy;
>> 176 G4double PositP= sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
>> 177 G4ThreeVector Phot2Direction (-dirx*Eratio, -diry*Eratio,
>> 178 (PositP-dirz*Phot1Energy)/Phot2Energy);
>> 179 Phot2Direction.unit();
>> 180 Phot2Direction.rotateUz(PositDirection);
>> 181 // create G4DynamicParticle object for the particle2
>> 182 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
>> 183 Phot2Direction, Phot2Energy);
>> 184 vdp->push_back(aParticle2);
>> 185 return vdp;
197 } 186 }
198 187
199 //....oooOO0OOooo........oooOO0OOooo........oo 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200 189