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