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