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