1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4GammaConversionToMuons.cc,v 1.8 2010/10/26 14:15:40 vnivanch Exp $ |
---|
28 | // GEANT4 tag $Name: emhighenergy-V09-03-02 $ |
---|
29 | // |
---|
30 | // ------------ G4GammaConversionToMuons physics process ------ |
---|
31 | // by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002 |
---|
32 | // |
---|
33 | // |
---|
34 | // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ... |
---|
35 | // 25-10-04: migrade to new interfaces of ParticleChange (vi) |
---|
36 | // --------------------------------------------------------------------------- |
---|
37 | |
---|
38 | #include "G4GammaConversionToMuons.hh" |
---|
39 | #include "G4EnergyLossTables.hh" |
---|
40 | #include "G4UnitsTable.hh" |
---|
41 | #include "G4MuonPlus.hh" |
---|
42 | #include "G4MuonMinus.hh" |
---|
43 | |
---|
44 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
45 | |
---|
46 | using namespace std; |
---|
47 | |
---|
48 | G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName, |
---|
49 | G4ProcessType type):G4VDiscreteProcess (processName, type), |
---|
50 | LowestEnergyLimit (4*G4MuonPlus::MuonPlus()->GetPDGMass()), // 4*Mmuon |
---|
51 | HighestEnergyLimit(1e21*eV), // ok to 1e21eV=1e12GeV, then LPM suppression |
---|
52 | CrossSecFactor(1.) |
---|
53 | { |
---|
54 | SetProcessSubType(15); |
---|
55 | MeanFreePath = DBL_MAX; |
---|
56 | } |
---|
57 | |
---|
58 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
59 | |
---|
60 | // destructor |
---|
61 | |
---|
62 | G4GammaConversionToMuons::~G4GammaConversionToMuons() // (empty) destructor |
---|
63 | { } |
---|
64 | |
---|
65 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
66 | |
---|
67 | G4bool G4GammaConversionToMuons::IsApplicable( |
---|
68 | const G4ParticleDefinition& particle) |
---|
69 | { |
---|
70 | return ( &particle == G4Gamma::Gamma() ); |
---|
71 | } |
---|
72 | |
---|
73 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
74 | |
---|
75 | void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition&) |
---|
76 | // Build cross section and mean free path tables |
---|
77 | { //here no tables, just calling PrintInfoDefinition |
---|
78 | PrintInfoDefinition(); |
---|
79 | } |
---|
80 | |
---|
81 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
82 | |
---|
83 | G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack, |
---|
84 | G4double, G4ForceCondition*) |
---|
85 | |
---|
86 | // returns the photon mean free path in GEANT4 internal units |
---|
87 | // (MeanFreePath is a private member of the class) |
---|
88 | |
---|
89 | { |
---|
90 | const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); |
---|
91 | G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); |
---|
92 | G4Material* aMaterial = aTrack.GetMaterial(); |
---|
93 | |
---|
94 | if (GammaEnergy < LowestEnergyLimit) |
---|
95 | MeanFreePath = DBL_MAX; |
---|
96 | else |
---|
97 | MeanFreePath = ComputeMeanFreePath(GammaEnergy,aMaterial); |
---|
98 | |
---|
99 | return MeanFreePath; |
---|
100 | } |
---|
101 | |
---|
102 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
103 | |
---|
104 | G4double G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy, |
---|
105 | G4Material* aMaterial) |
---|
106 | |
---|
107 | // computes and returns the photon mean free path in GEANT4 internal units |
---|
108 | { |
---|
109 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
110 | const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); |
---|
111 | |
---|
112 | G4double SIGMA = 0 ; |
---|
113 | |
---|
114 | for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ ) |
---|
115 | { |
---|
116 | G4double AtomicZ = (*theElementVector)[i]->GetZ(); |
---|
117 | G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole); |
---|
118 | SIGMA += NbOfAtomsPerVolume[i] * |
---|
119 | ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA); |
---|
120 | } |
---|
121 | return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX; |
---|
122 | } |
---|
123 | |
---|
124 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
125 | |
---|
126 | G4double G4GammaConversionToMuons::GetCrossSectionPerAtom( |
---|
127 | const G4DynamicParticle* aDynamicGamma, |
---|
128 | G4Element* anElement) |
---|
129 | |
---|
130 | // gives the total cross section per atom in GEANT4 internal units |
---|
131 | { |
---|
132 | G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); |
---|
133 | G4double AtomicZ = anElement->GetZ(); |
---|
134 | G4double AtomicA = anElement->GetA()/(g/mole); |
---|
135 | G4double crossSection = |
---|
136 | ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA); |
---|
137 | return crossSection; |
---|
138 | } |
---|
139 | |
---|
140 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
141 | |
---|
142 | G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom( |
---|
143 | G4double Egam, G4double Z, G4double A) |
---|
144 | |
---|
145 | // Calculates the microscopic cross section in GEANT4 internal units. |
---|
146 | // Total cross section parametrisation from H.Burkhardt |
---|
147 | // It gives a good description at any energy (from 0 to 10**21 eV) |
---|
148 | { static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass(); |
---|
149 | static const G4double Mele=electron_mass_c2; |
---|
150 | static const G4double Rc=elm_coupling/Mmuon; // classical particle radius |
---|
151 | static const G4double sqrte=sqrt(exp(1.)); |
---|
152 | static const G4double PowSat=-0.88; |
---|
153 | |
---|
154 | static G4double CrossSection = 0.0 ; |
---|
155 | |
---|
156 | if ( A < 1. ) return 0; |
---|
157 | if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0 |
---|
158 | |
---|
159 | static G4double EgamLast=0,Zlast=0,PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr, |
---|
160 | Wsatur,sigfac; |
---|
161 | |
---|
162 | if(Zlast==Z && Egam==EgamLast) return CrossSection; // already calculated |
---|
163 | EgamLast=Egam; |
---|
164 | |
---|
165 | if(Zlast!=Z) // new element |
---|
166 | { Zlast=Z; |
---|
167 | if(Z==1) // special case of Hydrogen |
---|
168 | { B=202.4; |
---|
169 | Dn=1.49; |
---|
170 | } |
---|
171 | else |
---|
172 | { B=183.; |
---|
173 | Dn=1.54*pow(A,0.27); |
---|
174 | } |
---|
175 | Zthird=pow(Z,-1./3.); // Z**(-1/3) |
---|
176 | Winfty=B*Zthird*Mmuon/(Dn*Mele); |
---|
177 | WMedAppr=1./(4.*Dn*sqrte*Mmuon); |
---|
178 | Wsatur=Winfty/WMedAppr; |
---|
179 | sigfac=4.*fine_structure_const*Z*Z*Rc*Rc; |
---|
180 | PowThres=1.479+0.00799*Dn; |
---|
181 | Ecor=-18.+4347./(B*Zthird); |
---|
182 | } |
---|
183 | G4double CorFuc=1.+.04*log(1.+Ecor/Egam); |
---|
184 | G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+ |
---|
185 | pow(Egam,PowSat),1./PowSat); // threshold and saturation |
---|
186 | CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg); |
---|
187 | CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1) |
---|
188 | return CrossSection; |
---|
189 | } |
---|
190 | |
---|
191 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
192 | |
---|
193 | void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac) |
---|
194 | // Set the factor to artificially increase the cross section |
---|
195 | { CrossSecFactor=fac; |
---|
196 | G4cout << "The cross section for GammaConversionToMuons is artificially " |
---|
197 | << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl; |
---|
198 | } |
---|
199 | |
---|
200 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
201 | |
---|
202 | G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt( |
---|
203 | const G4Track& aTrack, |
---|
204 | const G4Step& aStep) |
---|
205 | // |
---|
206 | // generation of gamma->mu+mu- |
---|
207 | // |
---|
208 | { |
---|
209 | aParticleChange.Initialize(aTrack); |
---|
210 | G4Material* aMaterial = aTrack.GetMaterial(); |
---|
211 | |
---|
212 | static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass(); |
---|
213 | static const G4double Mele=electron_mass_c2; |
---|
214 | static const G4double sqrte=sqrt(exp(1.)); |
---|
215 | |
---|
216 | // current Gamma energy and direction, return if energy too low |
---|
217 | const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle(); |
---|
218 | G4double Egam = aDynamicGamma->GetKineticEnergy(); |
---|
219 | if (Egam < 4*Mmuon) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); |
---|
220 | G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection(); |
---|
221 | |
---|
222 | // select randomly one element constituting the material |
---|
223 | const G4Element& anElement = *SelectRandomAtom(aDynamicGamma, aMaterial); |
---|
224 | G4double Z = anElement.GetZ(); |
---|
225 | G4double A = anElement.GetA()/(g/mole); |
---|
226 | |
---|
227 | static G4double Zlast=0,B,Dn,Zthird,Winfty,A027,C1Num2,C2Term2; |
---|
228 | if(Zlast!=Z) // the element has changed |
---|
229 | { Zlast=Z; |
---|
230 | if(Z==1) // special case of Hydrogen |
---|
231 | { B=202.4; |
---|
232 | Dn=1.49; |
---|
233 | } |
---|
234 | else |
---|
235 | { B=183.; |
---|
236 | Dn=1.54*pow(A,0.27); |
---|
237 | } |
---|
238 | Zthird=pow(Z,-1./3.); // Z**(-1/3) |
---|
239 | Winfty=B*Zthird*Mmuon/(Dn*Mele); |
---|
240 | A027=pow(A,0.27); |
---|
241 | G4double C1Num=0.35*A027; |
---|
242 | C1Num2=C1Num*C1Num; |
---|
243 | C2Term2=Mele/(183.*Zthird*Mmuon); |
---|
244 | } |
---|
245 | |
---|
246 | G4double GammaMuonInv=Mmuon/Egam; |
---|
247 | G4double sqrtx=sqrt(.25-GammaMuonInv); |
---|
248 | G4double xmax=.5+sqrtx; |
---|
249 | G4double xmin=.5-sqrtx; |
---|
250 | |
---|
251 | // generate xPlus according to the differential cross section by rejection |
---|
252 | G4double Ds2=(Dn*sqrte-2.); |
---|
253 | G4double sBZ=sqrte*B*Zthird/Mele; |
---|
254 | G4double LogWmaxInv=1./log(Winfty*(1.+2.*Ds2*GammaMuonInv) |
---|
255 | /(1.+2.*sBZ*Mmuon*GammaMuonInv)); |
---|
256 | G4double xPlus,xMinus,xPM,result,W; |
---|
257 | do |
---|
258 | { xPlus=xmin+G4UniformRand()*(xmax-xmin); |
---|
259 | xMinus=1.-xPlus; |
---|
260 | xPM=xPlus*xMinus; |
---|
261 | G4double del=Mmuon*Mmuon/(2.*Egam*xPM); |
---|
262 | W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del); |
---|
263 | if(W<1.) W=1.; // to avoid negative cross section at xmin |
---|
264 | G4double xxp=1.-4./3.*xPM; // the main xPlus dependence |
---|
265 | result=xxp*log(W)*LogWmaxInv; |
---|
266 | if(result>1.) { |
---|
267 | G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" |
---|
268 | << " in dSigxPlusGen, result=" << result << " > 1" << G4endl; |
---|
269 | } |
---|
270 | } |
---|
271 | while (G4UniformRand() > result); |
---|
272 | |
---|
273 | // now generate the angular variables via the auxilary variables t,psi,rho |
---|
274 | G4double t; |
---|
275 | G4double psi; |
---|
276 | G4double rho; |
---|
277 | |
---|
278 | G4double thetaPlus,thetaMinus,phiHalf; // final angular variables |
---|
279 | |
---|
280 | do // t, psi, rho generation start (while angle < pi) |
---|
281 | { |
---|
282 | //generate t by the rejection method |
---|
283 | G4double C1=C1Num2* GammaMuonInv/xPM; |
---|
284 | G4double f1_max=(1.-xPM) / (1.+C1); |
---|
285 | G4double f1; // the probability density |
---|
286 | do |
---|
287 | { t=G4UniformRand(); |
---|
288 | f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t)); |
---|
289 | if(f1<0 || f1> f1_max) // should never happend |
---|
290 | { |
---|
291 | G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" |
---|
292 | << "outside allowed range f1=" << f1 << " is set to zero" |
---|
293 | << G4endl; |
---|
294 | f1 = 0.0; |
---|
295 | } |
---|
296 | } |
---|
297 | while ( G4UniformRand()*f1_max > f1); |
---|
298 | // generate psi by the rejection method |
---|
299 | G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t)); |
---|
300 | |
---|
301 | // long version |
---|
302 | G4double f2; |
---|
303 | do |
---|
304 | { psi=2.*pi*G4UniformRand(); |
---|
305 | f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi)); |
---|
306 | if(f2<0 || f2> f2_max) // should never happend |
---|
307 | { |
---|
308 | G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" |
---|
309 | << "outside allowed range f2=" << f2 << " is set to zero" |
---|
310 | << G4endl; |
---|
311 | f2 = 0.0; |
---|
312 | } |
---|
313 | } |
---|
314 | while ( G4UniformRand()*f2_max > f2); |
---|
315 | |
---|
316 | // generate rho by direct transformation |
---|
317 | G4double C2Term1=GammaMuonInv/(2.*xPM*t); |
---|
318 | G4double C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.); |
---|
319 | G4double rhomax=1.9/A027*(1./t-1.); |
---|
320 | G4double beta=log( (C2+pow(rhomax,4.))/C2 ); |
---|
321 | rho=pow(C2 *( exp(beta*G4UniformRand())-1. ) ,0.25); |
---|
322 | |
---|
323 | //now get from t and psi the kinematical variables |
---|
324 | G4double u=sqrt(1./t-1.); |
---|
325 | G4double xiHalf=0.5*rho*cos(psi); |
---|
326 | phiHalf=0.5*rho/u*sin(psi); |
---|
327 | |
---|
328 | thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus; |
---|
329 | thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus; |
---|
330 | |
---|
331 | } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi); |
---|
332 | |
---|
333 | // now construct the vectors |
---|
334 | // azimuthal symmetry, take phi0 at random between 0 and 2 pi |
---|
335 | G4double phi0=2.*pi*G4UniformRand(); |
---|
336 | G4double EPlus=xPlus*Egam; |
---|
337 | G4double EMinus=xMinus*Egam; |
---|
338 | |
---|
339 | // mu+ mu- directions for gamma in z-direction |
---|
340 | G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf), |
---|
341 | sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) ); |
---|
342 | G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf), |
---|
343 | -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) ); |
---|
344 | // rotate to actual gamma direction |
---|
345 | MuPlusDirection.rotateUz(GammaDirection); |
---|
346 | MuMinusDirection.rotateUz(GammaDirection); |
---|
347 | aParticleChange.SetNumberOfSecondaries(2); |
---|
348 | // create G4DynamicParticle object for the particle1 |
---|
349 | G4DynamicParticle* aParticle1= new G4DynamicParticle( |
---|
350 | G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon); |
---|
351 | aParticleChange.AddSecondary(aParticle1); |
---|
352 | // create G4DynamicParticle object for the particle2 |
---|
353 | G4DynamicParticle* aParticle2= new G4DynamicParticle( |
---|
354 | G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon); |
---|
355 | aParticleChange.AddSecondary(aParticle2); |
---|
356 | // |
---|
357 | // Kill the incident photon |
---|
358 | // |
---|
359 | aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ; |
---|
360 | aParticleChange.ProposeEnergy( 0. ) ; |
---|
361 | aParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
---|
362 | // Reset NbOfInteractionLengthLeft and return aParticleChange |
---|
363 | return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); |
---|
364 | } |
---|
365 | |
---|
366 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
367 | |
---|
368 | G4Element* G4GammaConversionToMuons::SelectRandomAtom( |
---|
369 | const G4DynamicParticle* aDynamicGamma, |
---|
370 | G4Material* aMaterial) |
---|
371 | { |
---|
372 | // select randomly 1 element within the material, invoked by PostStepDoIt |
---|
373 | |
---|
374 | const G4int NumberOfElements = aMaterial->GetNumberOfElements(); |
---|
375 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
376 | if (NumberOfElements == 1) return (*theElementVector)[0]; |
---|
377 | |
---|
378 | const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); |
---|
379 | |
---|
380 | G4double PartialSumSigma = 0. ; |
---|
381 | G4double rval = G4UniformRand()/MeanFreePath; |
---|
382 | |
---|
383 | |
---|
384 | for ( G4int i=0 ; i < NumberOfElements ; i++ ) |
---|
385 | { PartialSumSigma += NbOfAtomsPerVolume[i] * |
---|
386 | GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]); |
---|
387 | if (rval <= PartialSumSigma) return ((*theElementVector)[i]); |
---|
388 | } |
---|
389 | G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName() |
---|
390 | << "' has no elements, NULL pointer returned." << G4endl; |
---|
391 | return NULL; |
---|
392 | } |
---|
393 | |
---|
394 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... |
---|
395 | |
---|
396 | void G4GammaConversionToMuons::PrintInfoDefinition() |
---|
397 | { |
---|
398 | G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= "; |
---|
399 | G4cout << G4endl << GetProcessName() << ": " << comments |
---|
400 | << GetProcessSubType() << G4endl; |
---|
401 | G4cout << " good cross section parametrization from " |
---|
402 | << G4BestUnit(LowestEnergyLimit,"Energy") |
---|
403 | << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl; |
---|
404 | } |
---|
405 | |
---|
406 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|