| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // ------------ G4AnnihiToMuPair physi
29 // by H.Burkhardt, S. Kelner and R. Ko
30 // -------------------------------------------
31 //
32 //....oooOO0OOooo........oooOO0OOooo........oo
33 //
34 // 27.01.03 : first implementation (hbu)
35 // 04.02.03 : cosmetic simplifications (mma)
36 // 25.10.04 : migrade to new interfaces of Par
37 // 28.02.18 : cross section now including SSS
38 //
39 //....oooOO0OOooo........oooOO0OOooo........oo
40
41 #include "G4AnnihiToMuPair.hh"
42
43 #include "G4Exp.hh"
44 #include "G4LossTableManager.hh"
45 #include "G4Material.hh"
46 #include "G4MuonMinus.hh"
47 #include "G4MuonPlus.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4Positron.hh"
50 #include "G4Step.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4TauMinus.hh"
53 #include "G4TauPlus.hh"
54 #include "G4ios.hh"
55 #include "Randomize.hh"
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4Str
60 G4ProcessType type):G4VDiscreteProcess (pr
61 {
62 //e+ Energy threshold
63 if(processName == "AnnihiToTauPair") {
64 SetProcessSubType(fAnnihilationToTauTau);
65 part1 = G4TauPlus::TauPlus();
66 part2 = G4TauMinus::TauMinus();
67 fInfo = "e+e->tau+tau-";
68 } else {
69 SetProcessSubType(fAnnihilationToMuMu);
70 part1 = G4MuonPlus::MuonPlus();
71 part2 = G4MuonMinus::MuonMinus();
72 }
73 fMass = part1->GetPDGMass();
74 fLowEnergyLimit = 2. * fMass * fMass / CLHEP
75
76 // model is ok up to 1000 TeV due to neglect
77 fHighEnergyLimit = 1000. * TeV;
78
79 fCurrentSigma = 0.0;
80 fCrossSecFactor = 1.;
81 fManager = G4LossTableManager::Instance();
82 fManager->Register(this);
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oo
86
87 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empt
88 {
89 fManager->DeRegister(this);
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 G4bool G4AnnihiToMuPair::IsApplicable(const G4
95 {
96 return ( &particle == G4Positron::Positron()
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 void G4AnnihiToMuPair::BuildPhysicsTable(const
102 {
103 PrintInfoDefinition();
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 void G4AnnihiToMuPair::SetCrossSecFactor(G4dou
109 // Set the factor to artificially increase the
110 {
111 fCrossSecFactor = fac;
112 //G4cout << "The cross section for AnnihiToM
113 // << "increased by the CrossSecFactor
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117
118 G4double G4AnnihiToMuPair::ComputeCrossSection
119 // Calculates the microscopic cross section in
120 // It gives a good description from threshold
121 {
122 G4double rmuon = CLHEP::elm_coupling/fMass;
123 G4double sig0 = CLHEP::pi*rmuon*rmuon/3.;
124 const G4double pial = CLHEP::pi*CLHEP::fine_
125
126 if (e <= fLowEnergyLimit) return 0.0;
127
128 const G4double xi = fLowEnergyLimit/e;
129 const G4double piaxi = pial * std::sqrt(xi);
130 G4double sigma = sig0 * xi * (1. + xi*0.5);
131 //G4cout << "### xi= " << xi << " piaxi=" <<
132
133 // argument of the exponent below 0.1 or abo
134 // Sigma per electron * number of electrons
135 if(xi <= 1.0 - 100*piaxi*piaxi) {
136 sigma *= std::sqrt(1.0 - xi);
137 }
138 else if (xi >= 1.0 - 0.01 * piaxi * piaxi) {
139 sigma *= piaxi;
140 }
141 else {
142 sigma *= piaxi / (1. - G4Exp(-piaxi / std:
143 }
144 // G4cout << "### sigma= " << sigma << G4end
145 return sigma;
146 }
147
148 //....oooOO0OOooo........oooOO0OOooo........oo
149
150 G4double G4AnnihiToMuPair::ComputeCrossSection
151
152 {
153 return ComputeCrossSectionPerElectron(energy
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oo
157
158 G4double G4AnnihiToMuPair::CrossSectionPerVolu
159 const G4Material* aMaterial)
160 {
161 return ComputeCrossSectionPerElectron(energy
162 }
163
164 //....oooOO0OOooo........oooOO0OOooo........oo
165
166 G4double G4AnnihiToMuPair::GetMeanFreePath(con
167 G4d
168 // returns the positron mean free path in GEAN
169 {
170 const G4DynamicParticle* aDynamicPositron =
171 G4double energy = aDynamicPositron->GetTotal
172 const G4Material* aMaterial = aTrack.GetMate
173
174 // cross section before step
175 fCurrentSigma = CrossSectionPerVolume(energy
176
177 // increase the CrossSection by CrossSecFact
178 return (fCurrentSigma > 0.0) ? 1.0/(fCurrent
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oo
182
183 G4VParticleChange* G4AnnihiToMuPair::PostStepD
184
185 //
186 // generation of e+e- -> mu+mu-
187 //
188 {
189 aParticleChange.Initialize(aTrack);
190
191 // current Positron energy and direction, re
192 const G4DynamicParticle *aDynamicPositron =
193 const G4double Mele = CLHEP::electron_mass_c
194 G4double Epos = aDynamicPositron->GetTotalEn
195 G4double xs = CrossSectionPerVolume(Epos, aT
196
197 // test of cross section
198 if(xs > 0.0 && fCurrentSigma*G4UniformRand()
199 return G4VDiscreteProcess::PostStepDoIt(aT
200 }
201
202 const G4ThreeVector PosiDirection = aDynamic
203 G4double xi = fLowEnergyLimit/Epos; // xi is
204 // goes
205
206 // generate cost; probability function 1+cos
207 //
208 G4double cost;
209 do { cost = 2.*G4UniformRand()-1.; }
210 // Loop checking, 07-Aug-2015, Vladimir Ivan
211 while (2.*G4UniformRand() > 1.+xi+cost*cost*
212 G4double sint = std::sqrt(1.-cost*cost);
213
214 // generate phi
215 //
216 G4double phi = 2.*CLHEP::pi*G4UniformRand();
217
218 G4double Ecm = std::sqrt(0.5*Mele*(Epos+Me
219 G4double Pcm = std::sqrt(Ecm*Ecm - fMass*f
220 G4double beta = std::sqrt((Epos-Mele)/(Epos
221 G4double gamma = Ecm/Mele;
222 G4double Pt = Pcm*sint;
223
224 // energy and momentum of the muons in the L
225 //
226 G4double EmuPlus = gamma*(Ecm + cost*beta*
227 G4double EmuMinus = gamma*(Ecm - cost*beta*
228 G4double PmuPlusZ = gamma*(beta*Ecm + cost*
229 G4double PmuMinusZ = gamma*(beta*Ecm - cost*
230 G4double PmuPlusX = Pt*std::cos(phi);
231 G4double PmuPlusY = Pt*std::sin(phi);
232 G4double PmuMinusX =-PmuPlusX;
233 G4double PmuMinusY =-PmuPlusY;
234 // absolute momenta
235 G4double PmuPlus = std::sqrt(Pt*Pt+PmuPlusZ
236 G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinus
237
238 // mu+ mu- directions for Positron in z-dire
239 //
240 G4ThreeVector MuPlusDirection(PmuPlusX / Pmu
241 G4ThreeVector MuMinusDirection(PmuMinusX / P
242
243 // rotate to actual Positron direction
244 //
245 MuPlusDirection.rotateUz(PosiDirection);
246 MuMinusDirection.rotateUz(PosiDirection);
247
248 aParticleChange.SetNumberOfSecondaries(2);
249
250 // create G4DynamicParticle object for the p
251 auto aParticle1 = new G4DynamicParticle(part
252 aParticleChange.AddSecondary(aParticle1);
253 // create G4DynamicParticle object for the p
254 auto aParticle2 = new G4DynamicParticle(part
255 aParticleChange.AddSecondary(aParticle2);
256
257 // Kill the incident positron
258 //
259 aParticleChange.ProposeEnergy(0.);
260 aParticleChange.ProposeTrackStatus(fStopAndK
261
262 return &aParticleChange;
263 }
264
265 //....oooOO0OOooo........oooOO0OOooo........oo
266
267 void G4AnnihiToMuPair::PrintInfoDefinition()
268 {
269 G4String comments = fInfo + " annihilation,
270 G4cout << G4endl << GetProcessName() << ":
271 G4cout << " threshold at " << fLowEne
272 << " good description up to " << fHig
273 << G4endl;
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oo
277