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