| 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 // ------------ G4GammaConversionToMuo 27 // ------------ G4GammaConversionToMuons physics process ------
28 // by H.Burkhardt, S. Kelner and R. Ko 28 // by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
29 // 29 //
30 // 30 //
31 // 07-08-02: missprint in OR condition in DoIt 31 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
32 // 25-10-04: migrade to new interfaces of Part 32 // 25-10-04: migrade to new interfaces of ParticleChange (vi)
33 // ------------------------------------------- 33 // ---------------------------------------------------------------------------
34 34
35 #include "G4GammaConversionToMuons.hh" 35 #include "G4GammaConversionToMuons.hh"
36 36
37 #include "G4BetheHeitler5DModel.hh" 37 #include "G4BetheHeitler5DModel.hh"
38 #include "G4Electron.hh" 38 #include "G4Electron.hh"
39 #include "G4EmParameters.hh" 39 #include "G4EmParameters.hh"
40 #include "G4EmProcessSubType.hh" 40 #include "G4EmProcessSubType.hh"
41 #include "G4Exp.hh" 41 #include "G4Exp.hh"
42 #include "G4Gamma.hh" 42 #include "G4Gamma.hh"
43 #include "G4Log.hh" 43 #include "G4Log.hh"
44 #include "G4LossTableManager.hh" 44 #include "G4LossTableManager.hh"
45 #include "G4MuonMinus.hh" 45 #include "G4MuonMinus.hh"
46 #include "G4MuonPlus.hh" 46 #include "G4MuonPlus.hh"
47 #include "G4NistManager.hh" 47 #include "G4NistManager.hh"
48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
49 #include "G4Positron.hh" 49 #include "G4Positron.hh"
50 #include "G4ProductionCutsTable.hh" 50 #include "G4ProductionCutsTable.hh"
51 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
52 #include "G4UnitsTable.hh" 52 #include "G4UnitsTable.hh"
53 53
54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
55 55
56 static const G4double sqrte = std::sqrt(std::e 56 static const G4double sqrte = std::sqrt(std::exp(1.));
57 static const G4double PowSat = -0.88; 57 static const G4double PowSat = -0.88;
58 58
59 G4GammaConversionToMuons::G4GammaConversionToM 59 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
60 G4ProcessType type) 60 G4ProcessType type)
61 : G4VDiscreteProcess (processName, type), 61 : G4VDiscreteProcess (processName, type),
62 Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass() 62 Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
63 Rc(CLHEP::elm_coupling / Mmuon), 63 Rc(CLHEP::elm_coupling / Mmuon),
64 LimitEnergy(5. * Mmuon), 64 LimitEnergy(5. * Mmuon),
65 LowestEnergyLimit(2. * Mmuon), 65 LowestEnergyLimit(2. * Mmuon),
66 HighestEnergyLimit(1e12 * CLHEP::GeV), // 66 HighestEnergyLimit(1e12 * CLHEP::GeV), // ok to 1e12GeV, then LPM suppression
67 theGamma(G4Gamma::Gamma()), 67 theGamma(G4Gamma::Gamma()),
68 theMuonPlus(G4MuonPlus::MuonPlus()), 68 theMuonPlus(G4MuonPlus::MuonPlus()),
69 theMuonMinus(G4MuonMinus::MuonMinus()) 69 theMuonMinus(G4MuonMinus::MuonMinus())
70 { 70 {
71 SetProcessSubType(fGammaConversionToMuMu); 71 SetProcessSubType(fGammaConversionToMuMu);
72 fManager = G4LossTableManager::Instance(); 72 fManager = G4LossTableManager::Instance();
73 fManager->Register(this); 73 fManager->Register(this);
74 } 74 }
75 75
76 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
77 77
78 G4GammaConversionToMuons::~G4GammaConversionTo 78 G4GammaConversionToMuons::~G4GammaConversionToMuons()
79 { 79 {
80 fManager->DeRegister(this); 80 fManager->DeRegister(this);
81 } 81 }
82 82
83 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
84 84
85 G4bool G4GammaConversionToMuons::IsApplicable( 85 G4bool G4GammaConversionToMuons::IsApplicable(const G4ParticleDefinition& part)
86 { 86 {
87 return (&part == theGamma); 87 return (&part == theGamma);
88 } 88 }
89 89
90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 91
92 void G4GammaConversionToMuons::BuildPhysicsTab 92 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition& p)
93 { 93 {
94 Energy5DLimit = G4EmParameters::Instance()-> 94 Energy5DLimit = G4EmParameters::Instance()->MaxEnergyFor5DMuPair();
95 95
96 auto table = G4Material::GetMaterialTable(); 96 auto table = G4Material::GetMaterialTable();
97 std::size_t nelm = 0; 97 std::size_t nelm = 0;
98 for (auto const& mat : *table) { 98 for (auto const& mat : *table) {
99 std::size_t n = mat->GetNumberOfElements() 99 std::size_t n = mat->GetNumberOfElements();
100 nelm = std::max(nelm, n); 100 nelm = std::max(nelm, n);
101 } 101 }
102 temp.resize(nelm, 0); 102 temp.resize(nelm, 0);
103 103
104 if (Energy5DLimit > 0.0 && nullptr != f5Dmod 104 if (Energy5DLimit > 0.0 && nullptr != f5Dmodel) {
105 f5Dmodel = new G4BetheHeitler5DModel(); 105 f5Dmodel = new G4BetheHeitler5DModel();
106 f5Dmodel->SetLeptonPair(theMuonPlus, theMu 106 f5Dmodel->SetLeptonPair(theMuonPlus, theMuonMinus);
107 const std::size_t numElems = G4ProductionC 107 const std::size_t numElems = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
108 const G4DataVector cuts(numElems); 108 const G4DataVector cuts(numElems);
109 f5Dmodel->Initialise(&p, cuts); 109 f5Dmodel->Initialise(&p, cuts);
110 } 110 }
111 PrintInfoDefinition(); 111 PrintInfoDefinition();
112 } 112 }
113 113
114 //....oooOO0OOooo........oooOO0OOooo........oo 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 115
116 G4double G4GammaConversionToMuons::GetMeanFree 116 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack, G4double,
117 117 G4ForceCondition*)
118 // returns the photon mean free path in GEANT4 118 // returns the photon mean free path in GEANT4 internal units
119 { 119 {
120 const G4DynamicParticle* aDynamicGamma = aTr 120 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
121 G4double GammaEnergy = aDynamicGamma->GetKin 121 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
122 const G4Material* aMaterial = aTrack.GetMate 122 const G4Material* aMaterial = aTrack.GetMaterial();
123 return ComputeMeanFreePath(GammaEnergy, aMat 123 return ComputeMeanFreePath(GammaEnergy, aMaterial);
124 } 124 }
125 125
126 //....oooOO0OOooo........oooOO0OOooo........oo 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 127
128 G4double 128 G4double
129 G4GammaConversionToMuons::ComputeMeanFreePath( 129 G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
130 130 const G4Material* aMaterial)
131 131
132 // computes and returns the photon mean free p 132 // computes and returns the photon mean free path in GEANT4 internal units
133 { 133 {
134 if(GammaEnergy <= LowestEnergyLimit) { retur 134 if(GammaEnergy <= LowestEnergyLimit) { return DBL_MAX; }
135 const G4ElementVector* theElementVector = aM 135 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
136 const G4double* NbOfAtomsPerVolume = aMateri 136 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
137 137
138 G4double SIGMA = 0.0; 138 G4double SIGMA = 0.0;
139 G4double fact = 1.0; 139 G4double fact = 1.0;
140 G4double e = GammaEnergy; 140 G4double e = GammaEnergy;
141 // low energy approximation as in Bethe-Heit 141 // low energy approximation as in Bethe-Heitler model
142 if(e < LimitEnergy) { 142 if(e < LimitEnergy) {
143 G4double y = (e - LowestEnergyLimit)/(Limi 143 G4double y = (e - LowestEnergyLimit)/(LimitEnergy - LowestEnergyLimit);
144 fact = y*y; 144 fact = y*y;
145 e = LimitEnergy; 145 e = LimitEnergy;
146 } 146 }
147 147
148 for ( std::size_t i=0 ; i < aMaterial->GetNu 148 for ( std::size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
149 { 149 {
150 SIGMA += NbOfAtomsPerVolume[i] * fact * 150 SIGMA += NbOfAtomsPerVolume[i] * fact *
151 ComputeCrossSectionPerAtom(e, (*theEleme 151 ComputeCrossSectionPerAtom(e, (*theElementVector)[i]->GetZasInt());
152 } 152 }
153 return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX; 153 return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX;
154 } 154 }
155 155
156 //....oooOO0OOooo........oooOO0OOooo........oo 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 157
158 G4double G4GammaConversionToMuons::GetCrossSec 158 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
159 const G4Dyn 159 const G4DynamicParticle* aDynamicGamma,
160 const G4Ele 160 const G4Element* anElement)
161 161
162 // gives the total cross section per atom in G 162 // gives the total cross section per atom in GEANT4 internal units
163 { 163 {
164 return ComputeCrossSectionPerAtom(aDynamicG 164 return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(),
165 anElement 165 anElement->GetZasInt());
166 } 166 }
167 167
168 //....oooOO0OOooo........oooOO0OOooo........oo 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
169 169
170 G4double G4GammaConversionToMuons::ComputeCros 170 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
171 G4double Egam, G4int 171 G4double Egam, G4int Z)
172 172
173 // Calculates the microscopic cross section in 173 // Calculates the microscopic cross section in GEANT4 internal units.
174 // Total cross section parametrisation from H. 174 // Total cross section parametrisation from H.Burkhardt
175 // It gives a good description at any energy ( 175 // It gives a good description at any energy (from 0 to 10**21 eV)
176 { 176 {
177 if(Egam <= LowestEnergyLimit) { return 0.0; 177 if(Egam <= LowestEnergyLimit) { return 0.0; }
178 178
179 G4NistManager* nist = G4NistManager::Instanc 179 G4NistManager* nist = G4NistManager::Instance();
180 180
181 G4double PowThres, Ecor, B, Dn, Zthird, Winf 181 G4double PowThres, Ecor, B, Dn, Zthird, Winfty, WMedAppr, Wsatur, sigfac;
182 182
183 if (Z == 1) { // special case of Hydrogen 183 if (Z == 1) { // special case of Hydrogen
184 B = 202.4; 184 B = 202.4;
185 Dn = 1.49; 185 Dn = 1.49;
186 } 186 }
187 else { 187 else {
188 B = 183.; 188 B = 183.;
189 Dn = 1.54 * nist->GetA27(Z); 189 Dn = 1.54 * nist->GetA27(Z);
190 } 190 }
191 Zthird = 1. / nist->GetZ13(Z); // Z**(-1/3) 191 Zthird = 1. / nist->GetZ13(Z); // Z**(-1/3)
192 Winfty = B * Zthird * Mmuon / (Dn * electron 192 Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
193 WMedAppr = 1. / (4. * Dn * sqrte * Mmuon); 193 WMedAppr = 1. / (4. * Dn * sqrte * Mmuon);
194 Wsatur = Winfty / WMedAppr; 194 Wsatur = Winfty / WMedAppr;
195 sigfac = 4. * fine_structure_const * Z * Z * 195 sigfac = 4. * fine_structure_const * Z * Z * Rc * Rc;
196 PowThres = 1.479 + 0.00799 * Dn; 196 PowThres = 1.479 + 0.00799 * Dn;
197 Ecor = -18. + 4347. / (B * Zthird); 197 Ecor = -18. + 4347. / (B * Zthird);
198 198
199 G4double CorFuc = 1. + .04 * G4Log(1. + Ecor 199 G4double CorFuc = 1. + .04 * G4Log(1. + Ecor / Egam);
200 G4double Eg = 200 G4double Eg =
201 G4Exp(G4Log(1. - 4. * Mmuon / Egam) * PowT 201 G4Exp(G4Log(1. - 4. * Mmuon / Egam) * PowThres)
202 * G4Exp(G4Log(G4Exp(G4Log(Wsatur) * PowSat 202 * G4Exp(G4Log(G4Exp(G4Log(Wsatur) * PowSat) + G4Exp(G4Log(Egam) * PowSat)) / PowSat);
203 203
204 G4double CrossSection = 7. / 9. * sigfac * G 204 G4double CrossSection = 7. / 9. * sigfac * G4Log(1. + WMedAppr * CorFuc * Eg);
205 CrossSection *= CrossSecFactor; // increase 205 CrossSection *= CrossSecFactor; // increase the CrossSection by (by default 1)
206 return CrossSection; 206 return CrossSection;
207 } 207 }
208 208
209 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
210 210
211 void G4GammaConversionToMuons::SetCrossSecFact 211 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
212 // Set the factor to artificially increase the 212 // Set the factor to artificially increase the cross section
213 { 213 {
214 if (fac < 0.0) return; 214 if (fac < 0.0) return;
215 CrossSecFactor = fac; 215 CrossSecFactor = fac;
216 if (verboseLevel > 1) { 216 if (verboseLevel > 1) {
217 G4cout << "The cross section for GammaConv 217 G4cout << "The cross section for GammaConversionToMuons is artificially "
218 << "increased by the CrossSecFactor 218 << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
219 } 219 }
220 } 220 }
221 221
222 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
223 223
224 G4VParticleChange* G4GammaConversionToMuons::P 224 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
225 225 const G4Track& aTrack,
226 226 const G4Step& aStep)
227 // 227 //
228 // generation of gamma->mu+mu- 228 // generation of gamma->mu+mu-
229 // 229 //
230 { 230 {
231 aParticleChange.Initialize(aTrack); 231 aParticleChange.Initialize(aTrack);
232 const G4Material* aMaterial = aTrack.GetMate 232 const G4Material* aMaterial = aTrack.GetMaterial();
233 233
234 // current Gamma energy and direction, retur 234 // current Gamma energy and direction, return if energy too low
235 const G4DynamicParticle* aDynamicGamma = aTr 235 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
236 G4double Egam = aDynamicGamma->GetKineticEne 236 G4double Egam = aDynamicGamma->GetKineticEnergy();
237 if (Egam <= LowestEnergyLimit) { 237 if (Egam <= LowestEnergyLimit) {
238 return G4VDiscreteProcess::PostStepDoIt(aT 238 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
239 } 239 }
240 // 240 //
241 // Kill the incident photon 241 // Kill the incident photon
242 // 242 //
243 aParticleChange.ProposeMomentumDirection( 0. 243 aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
244 aParticleChange.ProposeEnergy( 0. ) ; 244 aParticleChange.ProposeEnergy( 0. ) ;
245 aParticleChange.ProposeTrackStatus( fStopAnd 245 aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
246 246
247 if (Egam <= Energy5DLimit) { 247 if (Egam <= Energy5DLimit) {
248 std::vector<G4DynamicParticle*> fvect; 248 std::vector<G4DynamicParticle*> fvect;
249 f5Dmodel->SampleSecondaries(&fvect, aTrack 249 f5Dmodel->SampleSecondaries(&fvect, aTrack.GetMaterialCutsCouple(),
250 aTrack.GetDynamicParticle(), 0.0, DBL_ 250 aTrack.GetDynamicParticle(), 0.0, DBL_MAX);
251 for(auto dp : fvect) { aParticleChange.Add 251 for(auto dp : fvect) { aParticleChange.AddSecondary(dp); }
252 return G4VDiscreteProcess::PostStepDoIt(aT 252 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
253 } 253 }
254 254
255 G4ParticleMomentum GammaDirection = aDynamic 255 G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
256 256
257 // select randomly one element constituting 257 // select randomly one element constituting the material
258 const G4Element* anElement = SelectRandomAto 258 const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
259 G4int Z = anElement->GetZasInt(); 259 G4int Z = anElement->GetZasInt();
260 G4NistManager* nist = G4NistManager::Instanc 260 G4NistManager* nist = G4NistManager::Instance();
261 261
262 G4double B, Dn; 262 G4double B, Dn;
263 G4double A027 = nist->GetA27(Z); 263 G4double A027 = nist->GetA27(Z);
264 264
265 if (Z == 1) { // special case of Hydrogen 265 if (Z == 1) { // special case of Hydrogen
266 B = 202.4; 266 B = 202.4;
267 Dn = 1.49; 267 Dn = 1.49;
268 } 268 }
269 else { 269 else {
270 B = 183.; 270 B = 183.;
271 Dn = 1.54 * A027; 271 Dn = 1.54 * A027;
272 } 272 }
273 G4double Zthird = 1. / nist->GetZ13(Z); // 273 G4double Zthird = 1. / nist->GetZ13(Z); // Z**(-1/3)
274 G4double Winfty = B * Zthird * Mmuon / (Dn * 274 G4double Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
275 275
276 G4double C1Num = 0.138 * A027; 276 G4double C1Num = 0.138 * A027;
277 G4double C1Num2 = C1Num * C1Num; 277 G4double C1Num2 = C1Num * C1Num;
278 G4double C2Term2 = electron_mass_c2 / (183. 278 G4double C2Term2 = electron_mass_c2 / (183. * Zthird * Mmuon);
279 279
280 G4double GammaMuonInv = Mmuon / Egam; 280 G4double GammaMuonInv = Mmuon / Egam;
281 281
282 // generate xPlus according to the different 282 // generate xPlus according to the differential cross section by rejection
283 G4double xmin = (Egam <= LimitEnergy) ? 0.5 283 G4double xmin = (Egam <= LimitEnergy) ? 0.5 : 0.5 - std::sqrt(0.25 - GammaMuonInv);
284 G4double xmax = 1. - xmin; 284 G4double xmax = 1. - xmin;
285 285
286 G4double Ds2 = (Dn * sqrte - 2.); 286 G4double Ds2 = (Dn * sqrte - 2.);
287 G4double sBZ = sqrte * B * Zthird / electron 287 G4double sBZ = sqrte * B * Zthird / electron_mass_c2;
288 G4double LogWmaxInv = 288 G4double LogWmaxInv =
289 1. / G4Log(Winfty * (1. + 2. * Ds2 * Gamma 289 1. / G4Log(Winfty * (1. + 2. * Ds2 * GammaMuonInv) / (1. + 2. * sBZ * Mmuon * GammaMuonInv));
290 G4double xPlus = 0.5; 290 G4double xPlus = 0.5;
291 G4double xMinus = 0.5; 291 G4double xMinus = 0.5;
292 G4double xPM = 0.25; 292 G4double xPM = 0.25;
293 293
294 G4int nn = 0; 294 G4int nn = 0;
295 const G4int nmax = 1000; 295 const G4int nmax = 1000;
296 296
297 // sampling for Egam > LimitEnergy 297 // sampling for Egam > LimitEnergy
298 if (xmin < 0.5) { 298 if (xmin < 0.5) {
299 G4double result, W; 299 G4double result, W;
300 do { 300 do {
301 xPlus = xmin + G4UniformRand() * (xmax - 301 xPlus = xmin + G4UniformRand() * (xmax - xmin);
302 xMinus = 1. - xPlus; 302 xMinus = 1. - xPlus;
303 xPM = xPlus * xMinus; 303 xPM = xPlus * xMinus;
304 G4double del = Mmuon * Mmuon / (2. * Ega 304 G4double del = Mmuon * Mmuon / (2. * Egam * xPM);
305 W = Winfty * (1. + Ds2 * del / Mmuon) / 305 W = Winfty * (1. + Ds2 * del / Mmuon) / (1. + sBZ * del);
306 G4double xxp = 1. - 4. / 3. * xPM; // t 306 G4double xxp = 1. - 4. / 3. * xPM; // the main xPlus dependence
307 result = (xxp > 0.) ? xxp * G4Log(W) * L 307 result = (xxp > 0.) ? xxp * G4Log(W) * LogWmaxInv : 0.0;
308 if (result > 1.) { 308 if (result > 1.) {
309 G4cout << "G4GammaConversionToMuons::P 309 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
310 << " in dSigxPlusGen, result=" 310 << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
311 } 311 }
312 ++nn; 312 ++nn;
313 if(nn >= nmax) { break; } 313 if(nn >= nmax) { break; }
314 } 314 }
315 // Loop checking, 07-Aug-2015, Vladimir Iv 315 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
316 while (G4UniformRand() > result); 316 while (G4UniformRand() > result);
317 } 317 }
318 318
319 // now generate the angular variables via th 319 // now generate the angular variables via the auxilary variables t,psi,rho
320 G4double t; 320 G4double t;
321 G4double psi; 321 G4double psi;
322 G4double rho; 322 G4double rho;
323 323
324 G4double a3 = (GammaMuonInv / (2. * xPM)); 324 G4double a3 = (GammaMuonInv / (2. * xPM));
325 G4double a33 = a3 * a3; 325 G4double a33 = a3 * a3;
326 G4double f1; 326 G4double f1;
327 G4double b1 = 1./(4.*C1Num2); 327 G4double b1 = 1./(4.*C1Num2);
328 G4double b3 = b1*b1*b1; 328 G4double b3 = b1*b1*b1;
329 G4double a21 = a33 + b1; 329 G4double a21 = a33 + b1;
330 330
331 G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G 331 G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3);
332 332
333 G4double thetaPlus,thetaMinus,phiHalf; // fi 333 G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
334 nn = 0; 334 nn = 0;
335 // t, psi, rho generation start (while angl 335 // t, psi, rho generation start (while angle < pi)
336 do { 336 do {
337 //generate t by the rejection method 337 //generate t by the rejection method
338 do { 338 do {
339 ++nn; 339 ++nn;
340 t=G4UniformRand(); 340 t=G4UniformRand();
341 G4double a34=a33/(t*t); 341 G4double a34=a33/(t*t);
342 G4double a22 = a34 + b1; 342 G4double a22 = a34 + b1;
343 if(std::abs(b1)<0.0001*a34) { 343 if(std::abs(b1)<0.0001*a34) {
344 // special case of a34=a22 because of 344 // special case of a34=a22 because of logarithm accuracy
345 f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a3 345 f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34);
346 } 346 }
347 else { 347 else {
348 f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1 348 f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3);
349 } 349 }
350 if (f1 < 0.0 || f1 > f1_max) { // shoul 350 if (f1 < 0.0 || f1 > f1_max) { // should never happend
351 G4cout << "G4GammaConversionToMuons::P 351 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
352 << "outside allowed range f1=" << f1 352 << "outside allowed range f1=" << f1
353 << " is set to zero, a34 = "<< a34 << 353 << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"."
354 << G4endl; 354 << G4endl;
355 f1 = 0.0; 355 f1 = 0.0;
356 } 356 }
357 if(nn > nmax) { break; } 357 if(nn > nmax) { break; }
358 // Loop checking, 07-Aug-2015, Vladimir 358 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
359 } while ( G4UniformRand()*f1_max > f1); 359 } while ( G4UniformRand()*f1_max > f1);
360 // generate psi by the rejection method 360 // generate psi by the rejection method
361 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t)) 361 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
362 // long version 362 // long version
363 G4double f2; 363 G4double f2;
364 do { 364 do {
365 ++nn; 365 ++nn;
366 psi=twopi*G4UniformRand(); 366 psi=twopi*G4UniformRand();
367 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+std::co 367 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+std::cos(2.*psi));
368 if(f2<0 || f2> f2_max) { // should never 368 if(f2<0 || f2> f2_max) { // should never happend
369 G4cout << "G4GammaConversionToMuons::P 369 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
370 << "outside allowed range f2=" 370 << "outside allowed range f2=" << f2 << " is set to zero" << G4endl;
371 f2 = 0.0; 371 f2 = 0.0;
372 } 372 }
373 if(nn >= nmax) { break; } 373 if(nn >= nmax) { break; }
374 // Loop checking, 07-Aug-2015, Vladimir 374 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
375 } while ( G4UniformRand()*f2_max > f2); 375 } while ( G4UniformRand()*f2_max > f2);
376 376
377 // generate rho by direct transformation 377 // generate rho by direct transformation
378 G4double C2Term1=GammaMuonInv/(2.*xPM*t); 378 G4double C2Term1=GammaMuonInv/(2.*xPM*t);
379 G4double C22 = C2Term1*C2Term1+C2Term2*C2T 379 G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2;
380 G4double C2=4.*C22*C22/std::sqrt(xPM); 380 G4double C2=4.*C22*C22/std::sqrt(xPM);
381 G4double rhomax=(1./t-1.)*1.9/A027; 381 G4double rhomax=(1./t-1.)*1.9/A027;
382 G4double beta=G4Log( (C2+rhomax*rhomax*rho 382 G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
383 rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4Uniform 383 rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
384 384
385 //now get from t and psi the kinematical v 385 //now get from t and psi the kinematical variables
386 G4double u=std::sqrt(1./t-1.); 386 G4double u=std::sqrt(1./t-1.);
387 G4double xiHalf=0.5*rho*std::cos(psi); 387 G4double xiHalf=0.5*rho*std::cos(psi);
388 phiHalf=0.5*rho/u*std::sin(psi); 388 phiHalf=0.5*rho/u*std::sin(psi);
389 389
390 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus; 390 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
391 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus; 391 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
392 392
393 // protection against infinite loop 393 // protection against infinite loop
394 if(nn > nmax) { 394 if(nn > nmax) {
395 if(std::abs(thetaPlus)>pi) { thetaPlus = 395 if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
396 if(std::abs(thetaMinus)>pi) { thetaMinus 396 if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
397 } 397 }
398 398
399 // Loop checking, 07-Aug-2015, Vladimir Iv 399 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
400 } while ( std::abs(thetaPlus)>pi || std::abs 400 } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
401 401
402 // now construct the vectors 402 // now construct the vectors
403 // azimuthal symmetry, take phi0 at random b 403 // azimuthal symmetry, take phi0 at random between 0 and 2 pi
404 G4double phi0=twopi*G4UniformRand(); 404 G4double phi0=twopi*G4UniformRand();
405 G4double EPlus=xPlus*Egam; 405 G4double EPlus=xPlus*Egam;
406 G4double EMinus=xMinus*Egam; 406 G4double EMinus=xMinus*Egam;
407 407
408 // mu+ mu- directions for gamma in z-directi 408 // mu+ mu- directions for gamma in z-direction
409 G4ThreeVector MuPlusDirection ( std::sin(th 409 G4ThreeVector MuPlusDirection ( std::sin(thetaPlus) *std::cos(phi0+phiHalf),
410 std::sin(thetaPlus) *std:: 410 std::sin(thetaPlus) *std::sin(phi0+phiHalf), std::cos(thetaPlus) );
411 G4ThreeVector MuMinusDirection (-std::sin(th 411 G4ThreeVector MuMinusDirection (-std::sin(thetaMinus)*std::cos(phi0-phiHalf),
412 -std::sin(thetaMinus) *std:: 412 -std::sin(thetaMinus) *std::sin(phi0-phiHalf), std::cos(thetaMinus) );
413 // rotate to actual gamma direction 413 // rotate to actual gamma direction
414 MuPlusDirection.rotateUz(GammaDirection); 414 MuPlusDirection.rotateUz(GammaDirection);
415 MuMinusDirection.rotateUz(GammaDirection); 415 MuMinusDirection.rotateUz(GammaDirection);
416 416
417 // create G4DynamicParticle object for the p 417 // create G4DynamicParticle object for the particle1
418 auto aParticle1 = new G4DynamicParticle(theM 418 auto aParticle1 = new G4DynamicParticle(theMuonPlus, MuPlusDirection, EPlus - Mmuon);
419 aParticleChange.AddSecondary(aParticle1); 419 aParticleChange.AddSecondary(aParticle1);
420 // create G4DynamicParticle object for the p 420 // create G4DynamicParticle object for the particle2
421 auto aParticle2 = new G4DynamicParticle(theM 421 auto aParticle2 = new G4DynamicParticle(theMuonMinus, MuMinusDirection, EMinus - Mmuon);
422 aParticleChange.AddSecondary(aParticle2); 422 aParticleChange.AddSecondary(aParticle2);
423 // Reset NbOfInteractionLengthLeft and retu 423 // Reset NbOfInteractionLengthLeft and return aParticleChange
424 return G4VDiscreteProcess::PostStepDoIt( aTr 424 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
425 } 425 }
426 426
427 //....oooOO0OOooo........oooOO0OOooo........oo 427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
428 428
429 const G4Element* G4GammaConversionToMuons::Sel 429 const G4Element* G4GammaConversionToMuons::SelectRandomAtom(
430 const G4DynamicParticle* aDynamicGamma, 430 const G4DynamicParticle* aDynamicGamma,
431 const G4Material* aMaterial) 431 const G4Material* aMaterial)
432 { 432 {
433 // select randomly 1 element within the mate 433 // select randomly 1 element within the material, invoked by PostStepDoIt
434 434
435 const std::size_t NumberOfElements = aM 435 const std::size_t NumberOfElements = aMaterial->GetNumberOfElements();
436 const G4ElementVector* theElementVector = aM 436 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
437 const G4Element* elm = (*theElementVector)[0 437 const G4Element* elm = (*theElementVector)[0];
438 438
439 if (NumberOfElements > 1) { 439 if (NumberOfElements > 1) {
440 G4double e = std::max(aDynamicGamma->GetKi 440 G4double e = std::max(aDynamicGamma->GetKineticEnergy(), LimitEnergy);
441 const G4double* natom = aMaterial->GetVecN 441 const G4double* natom = aMaterial->GetVecNbOfAtomsPerVolume();
442 442
443 G4double sum = 0.; 443 G4double sum = 0.;
444 for (std::size_t i=0; i<NumberOfElements; 444 for (std::size_t i=0; i<NumberOfElements; ++i) {
445 elm = (*theElementVector)[i]; 445 elm = (*theElementVector)[i];
446 sum += natom[i]*ComputeCrossSectionPerAt 446 sum += natom[i]*ComputeCrossSectionPerAtom(e, elm->GetZasInt());
447 temp[i] = sum; 447 temp[i] = sum;
448 } 448 }
449 sum *= G4UniformRand(); 449 sum *= G4UniformRand();
450 for (std::size_t i=0; i<NumberOfElements; 450 for (std::size_t i=0; i<NumberOfElements; ++i) {
451 if(sum <= temp[i]) { 451 if(sum <= temp[i]) {
452 elm = (*theElementVector)[i]; 452 elm = (*theElementVector)[i];
453 break; 453 break;
454 } 454 }
455 } 455 }
456 } 456 }
457 return elm; 457 return elm;
458 } 458 }
459 459
460 //....oooOO0OOooo........oooOO0OOooo........oo 460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
461 461
462 void G4GammaConversionToMuons::PrintInfoDefini 462 void G4GammaConversionToMuons::PrintInfoDefinition()
463 { 463 {
464 G4String comments = "gamma->mu+mu- Bethe Hei 464 G4String comments = "gamma->mu+mu- Bethe Heitler process, SubType= ";
465 G4cout << G4endl << GetProcessName() << ": 465 G4cout << G4endl << GetProcessName() << ": " << comments << GetProcessSubType() << G4endl;
466 G4cout << " good cross section parame 466 G4cout << " good cross section parametrization from "
467 << G4BestUnit(LowestEnergyLimit, "Ene 467 << G4BestUnit(LowestEnergyLimit, "Energy") << " to " << HighestEnergyLimit / GeV
468 << " GeV for all Z." << G4endl; 468 << " GeV for all Z." << G4endl;
469 G4cout << " cross section factor: " < 469 G4cout << " cross section factor: " << CrossSecFactor << G4endl;
470 } 470 }
471 471
472 //....oooOO0OOooo........oooOO0OOooo........oo 472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
473 473