| 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 // ------------ G4AnnihiToMuPair physi 28 // ------------ G4AnnihiToMuPair physics process ------
29 // by H.Burkhardt, S. Kelner and R. Ko 29 // by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002
30 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
31 // 31 //
32 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......//
33 // 33 //
34 // 27.01.03 : first implementation (hbu) 34 // 27.01.03 : first implementation (hbu)
35 // 04.02.03 : cosmetic simplifications (mma) 35 // 04.02.03 : cosmetic simplifications (mma)
36 // 25.10.04 : migrade to new interfaces of Par 36 // 25.10.04 : migrade to new interfaces of ParticleChange (vi)
37 // 28.02.18 : cross section now including SSS 37 // 28.02.18 : cross section now including SSS threshold factor
38 // 38 //
39 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 40
41 #include "G4AnnihiToMuPair.hh" 41 #include "G4AnnihiToMuPair.hh"
42 42
43 #include "G4Exp.hh" << 43 #include "G4ios.hh"
44 #include "G4LossTableManager.hh" << 44 #include "Randomize.hh"
45 #include "G4Material.hh" <<
46 #include "G4MuonMinus.hh" <<
47 #include "G4MuonPlus.hh" <<
48 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh"
49 #include "G4Positron.hh" <<
50 #include "G4Step.hh" <<
51 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh"
52 #include "G4TauMinus.hh" << 47
>> 48 #include "G4Positron.hh"
>> 49 #include "G4MuonPlus.hh"
>> 50 #include "G4MuonMinus.hh"
53 #include "G4TauPlus.hh" 51 #include "G4TauPlus.hh"
54 #include "G4ios.hh" << 52 #include "G4TauMinus.hh"
55 #include "Randomize.hh" << 53 #include "G4Material.hh"
>> 54 #include "G4Step.hh"
>> 55 #include "G4LossTableManager.hh"
>> 56 #include "G4Exp.hh"
56 57
57 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 59
>> 60 using namespace std;
>> 61
59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4Str 62 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName,
60 G4ProcessType type):G4VDiscreteProcess (pr 63 G4ProcessType type):G4VDiscreteProcess (processName, type)
61 { 64 {
62 //e+ Energy threshold 65 //e+ Energy threshold
63 if(processName == "AnnihiToTauPair") { 66 if(processName == "AnnihiToTauPair") {
64 SetProcessSubType(fAnnihilationToTauTau); 67 SetProcessSubType(fAnnihilationToTauTau);
65 part1 = G4TauPlus::TauPlus(); 68 part1 = G4TauPlus::TauPlus();
66 part2 = G4TauMinus::TauMinus(); 69 part2 = G4TauMinus::TauMinus();
67 fInfo = "e+e->tau+tau-"; 70 fInfo = "e+e->tau+tau-";
68 } else { 71 } else {
69 SetProcessSubType(fAnnihilationToMuMu); 72 SetProcessSubType(fAnnihilationToMuMu);
70 part1 = G4MuonPlus::MuonPlus(); 73 part1 = G4MuonPlus::MuonPlus();
71 part2 = G4MuonMinus::MuonMinus(); 74 part2 = G4MuonMinus::MuonMinus();
72 } 75 }
73 fMass = part1->GetPDGMass(); 76 fMass = part1->GetPDGMass();
74 fLowEnergyLimit = 2. * fMass * fMass / CLHEP << 77 fLowEnergyLimit =
75 << 78 2.*fMass*fMass/CLHEP::electron_mass_c2 - CLHEP::electron_mass_c2;
76 // model is ok up to 1000 TeV due to neglect << 79
77 fHighEnergyLimit = 1000. * TeV; << 80 //model is ok up to 1000 TeV due to neglected Z-interference
78 << 81 fHighEnergyLimit = 1000.*TeV;
>> 82
79 fCurrentSigma = 0.0; 83 fCurrentSigma = 0.0;
80 fCrossSecFactor = 1.; 84 fCrossSecFactor = 1.;
81 fManager = G4LossTableManager::Instance(); 85 fManager = G4LossTableManager::Instance();
82 fManager->Register(this); 86 fManager->Register(this);
83 } 87 }
84 88
85 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 90
87 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empt 91 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor
88 { 92 {
89 fManager->DeRegister(this); 93 fManager->DeRegister(this);
90 } 94 }
91 95
92 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 97
94 G4bool G4AnnihiToMuPair::IsApplicable(const G4 98 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle)
95 { 99 {
96 return ( &particle == G4Positron::Positron() 100 return ( &particle == G4Positron::Positron() );
97 } 101 }
98 102
99 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 104
101 void G4AnnihiToMuPair::BuildPhysicsTable(const 105 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&)
102 { 106 {
103 PrintInfoDefinition(); 107 PrintInfoDefinition();
104 } 108 }
105 109
106 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 111
108 void G4AnnihiToMuPair::SetCrossSecFactor(G4dou 112 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac)
109 // Set the factor to artificially increase the 113 // Set the factor to artificially increase the cross section
110 { 114 {
111 fCrossSecFactor = fac; 115 fCrossSecFactor = fac;
112 //G4cout << "The cross section for AnnihiToM 116 //G4cout << "The cross section for AnnihiToMuPair is artificially "
113 // << "increased by the CrossSecFactor 117 // << "increased by the CrossSecFactor=" << fCrossSecFactor << G4endl;
114 } 118 }
115 119
116 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 121
118 G4double G4AnnihiToMuPair::ComputeCrossSection 122 G4double G4AnnihiToMuPair::ComputeCrossSectionPerElectron(const G4double e)
119 // Calculates the microscopic cross section in 123 // Calculates the microscopic cross section in GEANT4 internal units.
120 // It gives a good description from threshold 124 // It gives a good description from threshold to 1000 GeV
121 { 125 {
122 G4double rmuon = CLHEP::elm_coupling/fMass; 126 G4double rmuon = CLHEP::elm_coupling/fMass; //classical particle radius
123 G4double sig0 = CLHEP::pi*rmuon*rmuon/3.; 127 G4double sig0 = CLHEP::pi*rmuon*rmuon/3.; //constant in crossSection
124 const G4double pial = CLHEP::pi*CLHEP::fine_ 128 const G4double pial = CLHEP::pi*CLHEP::fine_structure_const; // pi * alphaQED
125 129
126 if (e <= fLowEnergyLimit) return 0.0; 130 if (e <= fLowEnergyLimit) return 0.0;
127 131
128 const G4double xi = fLowEnergyLimit/e; 132 const G4double xi = fLowEnergyLimit/e;
129 const G4double piaxi = pial * std::sqrt(xi); 133 const G4double piaxi = pial * std::sqrt(xi);
130 G4double sigma = sig0 * xi * (1. + xi*0.5); 134 G4double sigma = sig0 * xi * (1. + xi*0.5);
131 //G4cout << "### xi= " << xi << " piaxi=" << 135 //G4cout << "### xi= " << xi << " piaxi=" << piaxi << G4endl;
132 136
133 // argument of the exponent below 0.1 or abo 137 // argument of the exponent below 0.1 or above 10
134 // Sigma per electron * number of electrons 138 // Sigma per electron * number of electrons per atom
135 if(xi <= 1.0 - 100*piaxi*piaxi) { 139 if(xi <= 1.0 - 100*piaxi*piaxi) {
136 sigma *= std::sqrt(1.0 - xi); 140 sigma *= std::sqrt(1.0 - xi);
137 } << 141 } else if( xi >= 1.0 - 0.01*piaxi*piaxi) {
138 else if (xi >= 1.0 - 0.01 * piaxi * piaxi) { <<
139 sigma *= piaxi; 142 sigma *= piaxi;
>> 143 } else {
>> 144 sigma *= piaxi/(1. - G4Exp( -piaxi/std::sqrt(1-xi) ));
140 } 145 }
141 else { << 146 //G4cout << "### sigma= " << sigma << G4endl;
142 sigma *= piaxi / (1. - G4Exp(-piaxi / std: <<
143 } <<
144 // G4cout << "### sigma= " << sigma << G4end <<
145 return sigma; 147 return sigma;
146 } 148 }
147 149
148 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 151
150 G4double G4AnnihiToMuPair::ComputeCrossSection 152 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(const G4double energy,
151 153 const G4double Z)
152 { 154 {
153 return ComputeCrossSectionPerElectron(energy 155 return ComputeCrossSectionPerElectron(energy)*Z;
154 } 156 }
155 157
156 //....oooOO0OOooo........oooOO0OOooo........oo 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 159
158 G4double G4AnnihiToMuPair::CrossSectionPerVolu 160 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double energy,
159 const G4Material* aMaterial) 161 const G4Material* aMaterial)
160 { 162 {
161 return ComputeCrossSectionPerElectron(energy 163 return ComputeCrossSectionPerElectron(energy)*aMaterial->GetTotNbOfElectPerVolume();
162 } 164 }
163 165
164 //....oooOO0OOooo........oooOO0OOooo........oo 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165 167
166 G4double G4AnnihiToMuPair::GetMeanFreePath(con 168 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
167 G4d 169 G4double, G4ForceCondition*)
168 // returns the positron mean free path in GEAN 170 // returns the positron mean free path in GEANT4 internal units
169 { 171 {
170 const G4DynamicParticle* aDynamicPositron = 172 const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
171 G4double energy = aDynamicPositron->GetTotal 173 G4double energy = aDynamicPositron->GetTotalEnergy();
172 const G4Material* aMaterial = aTrack.GetMate 174 const G4Material* aMaterial = aTrack.GetMaterial();
173 175
174 // cross section before step 176 // cross section before step
175 fCurrentSigma = CrossSectionPerVolume(energy 177 fCurrentSigma = CrossSectionPerVolume(energy, aMaterial);
176 178
177 // increase the CrossSection by CrossSecFact 179 // increase the CrossSection by CrossSecFactor (default 1)
178 return (fCurrentSigma > 0.0) ? 1.0/(fCurrent 180 return (fCurrentSigma > 0.0) ? 1.0/(fCurrentSigma*fCrossSecFactor) : DBL_MAX;
179 } 181 }
180 182
181 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 184
183 G4VParticleChange* G4AnnihiToMuPair::PostStepD 185 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack,
184 186 const G4Step& aStep)
185 // 187 //
186 // generation of e+e- -> mu+mu- 188 // generation of e+e- -> mu+mu-
187 // 189 //
188 { 190 {
189 aParticleChange.Initialize(aTrack); 191 aParticleChange.Initialize(aTrack);
190 192
191 // current Positron energy and direction, re 193 // current Positron energy and direction, return if energy too low
192 const G4DynamicParticle *aDynamicPositron = 194 const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
193 const G4double Mele = CLHEP::electron_mass_c 195 const G4double Mele = CLHEP::electron_mass_c2;
194 G4double Epos = aDynamicPositron->GetTotalEn 196 G4double Epos = aDynamicPositron->GetTotalEnergy();
195 G4double xs = CrossSectionPerVolume(Epos, aT 197 G4double xs = CrossSectionPerVolume(Epos, aTrack.GetMaterial());
196 198
197 // test of cross section 199 // test of cross section
198 if(xs > 0.0 && fCurrentSigma*G4UniformRand() << 200 if(xs > 0.0 && fCurrentSigma*G4UniformRand() > xs)
199 return G4VDiscreteProcess::PostStepDoIt(aT << 201 {
200 } << 202 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
>> 203 }
201 204
202 const G4ThreeVector PosiDirection = aDynamic 205 const G4ThreeVector PosiDirection = aDynamicPositron->GetMomentumDirection();
203 G4double xi = fLowEnergyLimit/Epos; // xi is 206 G4double xi = fLowEnergyLimit/Epos; // xi is always less than 1,
204 // goes 207 // goes to 0 at high Epos
205 208
206 // generate cost; probability function 1+cos 209 // generate cost; probability function 1+cost**2 at high Epos
207 // 210 //
208 G4double cost; 211 G4double cost;
209 do { cost = 2.*G4UniformRand()-1.; } 212 do { cost = 2.*G4UniformRand()-1.; }
210 // Loop checking, 07-Aug-2015, Vladimir Ivan 213 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
211 while (2.*G4UniformRand() > 1.+xi+cost*cost* 214 while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) );
212 G4double sint = std::sqrt(1.-cost*cost); << 215 G4double sint = sqrt(1.-cost*cost);
213 216
214 // generate phi 217 // generate phi
215 // 218 //
216 G4double phi = 2.*CLHEP::pi*G4UniformRand(); 219 G4double phi = 2.*CLHEP::pi*G4UniformRand();
217 220
218 G4double Ecm = std::sqrt(0.5*Mele*(Epos+Me 221 G4double Ecm = std::sqrt(0.5*Mele*(Epos+Mele));
219 G4double Pcm = std::sqrt(Ecm*Ecm - fMass*f 222 G4double Pcm = std::sqrt(Ecm*Ecm - fMass*fMass);
220 G4double beta = std::sqrt((Epos-Mele)/(Epos 223 G4double beta = std::sqrt((Epos-Mele)/(Epos+Mele));
221 G4double gamma = Ecm/Mele; 224 G4double gamma = Ecm/Mele;
222 G4double Pt = Pcm*sint; 225 G4double Pt = Pcm*sint;
223 226
224 // energy and momentum of the muons in the L 227 // energy and momentum of the muons in the Lab
225 // 228 //
226 G4double EmuPlus = gamma*(Ecm + cost*beta* 229 G4double EmuPlus = gamma*(Ecm + cost*beta*Pcm);
227 G4double EmuMinus = gamma*(Ecm - cost*beta* 230 G4double EmuMinus = gamma*(Ecm - cost*beta*Pcm);
228 G4double PmuPlusZ = gamma*(beta*Ecm + cost* 231 G4double PmuPlusZ = gamma*(beta*Ecm + cost*Pcm);
229 G4double PmuMinusZ = gamma*(beta*Ecm - cost* 232 G4double PmuMinusZ = gamma*(beta*Ecm - cost*Pcm);
230 G4double PmuPlusX = Pt*std::cos(phi); 233 G4double PmuPlusX = Pt*std::cos(phi);
231 G4double PmuPlusY = Pt*std::sin(phi); 234 G4double PmuPlusY = Pt*std::sin(phi);
232 G4double PmuMinusX =-PmuPlusX; 235 G4double PmuMinusX =-PmuPlusX;
233 G4double PmuMinusY =-PmuPlusY; 236 G4double PmuMinusY =-PmuPlusY;
234 // absolute momenta 237 // absolute momenta
235 G4double PmuPlus = std::sqrt(Pt*Pt+PmuPlusZ 238 G4double PmuPlus = std::sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ );
236 G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinus 239 G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ);
237 240
238 // mu+ mu- directions for Positron in z-dire 241 // mu+ mu- directions for Positron in z-direction
239 // 242 //
240 G4ThreeVector MuPlusDirection(PmuPlusX / Pmu << 243 G4ThreeVector
241 G4ThreeVector MuMinusDirection(PmuMinusX / P << 244 MuPlusDirection(PmuPlusX/PmuPlus, PmuPlusY/PmuPlus, PmuPlusZ/PmuPlus);
>> 245 G4ThreeVector
>> 246 MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus);
242 247
243 // rotate to actual Positron direction 248 // rotate to actual Positron direction
244 // 249 //
245 MuPlusDirection.rotateUz(PosiDirection); 250 MuPlusDirection.rotateUz(PosiDirection);
246 MuMinusDirection.rotateUz(PosiDirection); 251 MuMinusDirection.rotateUz(PosiDirection);
247 252
248 aParticleChange.SetNumberOfSecondaries(2); 253 aParticleChange.SetNumberOfSecondaries(2);
249 254
250 // create G4DynamicParticle object for the p 255 // create G4DynamicParticle object for the particle1
251 auto aParticle1 = new G4DynamicParticle(part << 256 G4DynamicParticle* aParticle1 =
>> 257 new G4DynamicParticle(part1, MuPlusDirection, EmuPlus-fMass);
252 aParticleChange.AddSecondary(aParticle1); 258 aParticleChange.AddSecondary(aParticle1);
253 // create G4DynamicParticle object for the p 259 // create G4DynamicParticle object for the particle2
254 auto aParticle2 = new G4DynamicParticle(part << 260 G4DynamicParticle* aParticle2 =
>> 261 new G4DynamicParticle(part2, MuMinusDirection, EmuMinus-fMass);
255 aParticleChange.AddSecondary(aParticle2); 262 aParticleChange.AddSecondary(aParticle2);
256 263
257 // Kill the incident positron 264 // Kill the incident positron
258 // 265 //
259 aParticleChange.ProposeEnergy(0.); 266 aParticleChange.ProposeEnergy(0.);
260 aParticleChange.ProposeTrackStatus(fStopAndK 267 aParticleChange.ProposeTrackStatus(fStopAndKill);
261 268
262 return &aParticleChange; 269 return &aParticleChange;
263 } 270 }
264 271
265 //....oooOO0OOooo........oooOO0OOooo........oo 272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 273
267 void G4AnnihiToMuPair::PrintInfoDefinition() 274 void G4AnnihiToMuPair::PrintInfoDefinition()
268 { 275 {
269 G4String comments = fInfo + " annihilation, 276 G4String comments = fInfo + " annihilation, atomic e- at rest, SubType=";
270 G4cout << G4endl << GetProcessName() << ": << 277 G4cout << G4endl << GetProcessName() << ": " << comments
271 G4cout << " threshold at " << fLowEne << 278 << GetProcessSubType() << G4endl;
272 << " good description up to " << fHig << 279 G4cout << " threshold at " << fLowEnergyLimit/CLHEP::GeV << " GeV"
273 << G4endl; << 280 << " good description up to "
>> 281 << fHighEnergyLimit/CLHEP::TeV << " TeV for all Z." << G4endl;
274 } 282 }
275 283
276 //....oooOO0OOooo........oooOO0OOooo........oo 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 285