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 | // $Id: G4MuNuclearInteraction.cc,v 1.14 2010/09/08 08:59:29 vnivanch Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-03-ref-09 $ |
---|
28 | // |
---|
29 | // $Id: |
---|
30 | // -------------------------------------------------------------- |
---|
31 | // GEANT 4 class implementation file |
---|
32 | // |
---|
33 | // History: first implementation, based on object model of |
---|
34 | // 2nd December 1995, G.Cosmo |
---|
35 | // -------- G4MuNuclearInteraction physics process --------- |
---|
36 | // by Laszlo Urban, May 1998 |
---|
37 | // added simple model for hadronic vertex, J.P. Wellisch, November 1998 |
---|
38 | // -------------------------------------------------------------- |
---|
39 | // 26/10/1998: new corr.s from R.Kokoulin + cleanup , L.Urban |
---|
40 | // 23/01/2009 V.Ivanchenko Add deregistration |
---|
41 | // |
---|
42 | |
---|
43 | #include "G4MuNuclearInteraction.hh" |
---|
44 | #include "G4UnitsTable.hh" |
---|
45 | #include "G4HadronicProcessStore.hh" |
---|
46 | |
---|
47 | // static members ........ |
---|
48 | G4int G4MuNuclearInteraction::nzdat = 5 ; |
---|
49 | G4double G4MuNuclearInteraction::zdat[]={1.,4.,13.,29.,92.}; |
---|
50 | G4double G4MuNuclearInteraction::adat[]={1.01,9.01,26.98,63.55,238.03}; |
---|
51 | G4int G4MuNuclearInteraction::ntdat = 8 ; |
---|
52 | G4double G4MuNuclearInteraction::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8, |
---|
53 | 1.e9,1.e10}; |
---|
54 | G4int G4MuNuclearInteraction::NBIN = 1000 ; |
---|
55 | G4double G4MuNuclearInteraction::ya[1001]={0.}; |
---|
56 | G4double G4MuNuclearInteraction::proba[5][8][1001]={{{0.}}}; |
---|
57 | |
---|
58 | G4MuNuclearInteraction::G4MuNuclearInteraction(const G4String& processName) |
---|
59 | : G4VDiscreteProcess(processName), |
---|
60 | theMeanFreePathTable(0), |
---|
61 | theCrossSectionTable(0), |
---|
62 | LowestKineticEnergy (1.*GeV), |
---|
63 | HighestKineticEnergy (1000000.*TeV), |
---|
64 | TotBin(100), |
---|
65 | CutFixed ( 0.200*GeV), |
---|
66 | GramPerMole(g/mole), |
---|
67 | theMuonMinus ( G4MuonMinus::MuonMinus() ), |
---|
68 | theMuonPlus ( G4MuonPlus::MuonPlus() ), |
---|
69 | thePionZero (G4PionZero::PionZero() ) |
---|
70 | { |
---|
71 | SetProcessSubType(fHadronInelastic); |
---|
72 | G4HadronicProcessStore::Instance()->RegisterExtraProcess(this); |
---|
73 | } |
---|
74 | |
---|
75 | G4MuNuclearInteraction::~G4MuNuclearInteraction() |
---|
76 | { |
---|
77 | G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this); |
---|
78 | |
---|
79 | if (theMeanFreePathTable) { |
---|
80 | theMeanFreePathTable->clearAndDestroy(); |
---|
81 | delete theMeanFreePathTable; |
---|
82 | } |
---|
83 | if (theCrossSectionTable) { |
---|
84 | theCrossSectionTable->clearAndDestroy(); |
---|
85 | delete theCrossSectionTable; |
---|
86 | } |
---|
87 | |
---|
88 | if (&PartialSumSigma) { |
---|
89 | PartialSumSigma.clearAndDestroy(); |
---|
90 | } |
---|
91 | } |
---|
92 | |
---|
93 | void G4MuNuclearInteraction::SetPhysicsTableBining(G4double lowE, |
---|
94 | G4double highE, G4int nBins) |
---|
95 | { |
---|
96 | LowestKineticEnergy = lowE; HighestKineticEnergy = highE ; TotBin = nBins ; |
---|
97 | } |
---|
98 | |
---|
99 | |
---|
100 | G4bool G4MuNuclearInteraction::IsApplicable(const G4ParticleDefinition& particle) |
---|
101 | { |
---|
102 | return( (&particle == theMuonMinus)||(&particle == theMuonPlus)) ; |
---|
103 | } |
---|
104 | |
---|
105 | void G4MuNuclearInteraction::PreparePhysicsTable( |
---|
106 | const G4ParticleDefinition& aParticleType) |
---|
107 | { |
---|
108 | G4HadronicProcessStore::Instance() |
---|
109 | ->RegisterParticleForExtraProcess(this, &aParticleType); |
---|
110 | } |
---|
111 | |
---|
112 | void G4MuNuclearInteraction::BuildPhysicsTable( |
---|
113 | const G4ParticleDefinition& aParticleType) |
---|
114 | { |
---|
115 | G4HadronicProcessStore::Instance()->PrintInfo(&aParticleType); |
---|
116 | |
---|
117 | G4double LowEdgeEnergy , Value; |
---|
118 | G4PhysicsLogVector* ptrVector; |
---|
119 | |
---|
120 | if (theCrossSectionTable) { |
---|
121 | theCrossSectionTable->clearAndDestroy() ; |
---|
122 | delete theCrossSectionTable ; |
---|
123 | } |
---|
124 | |
---|
125 | // make tables for the sampling at initialization |
---|
126 | if (theMeanFreePathTable == 0) MakeSamplingTables(&aParticleType); |
---|
127 | |
---|
128 | theCrossSectionTable = new G4PhysicsTable (G4Element::GetNumberOfElements()); |
---|
129 | const G4ElementTable* theElementTable = G4Element::GetElementTable() ; |
---|
130 | G4double AtomicNumber,AtomicWeight ; |
---|
131 | |
---|
132 | for (size_t J=0; J < G4Element::GetNumberOfElements(); J++ ) |
---|
133 | { |
---|
134 | ptrVector = new G4PhysicsLogVector(LowestKineticEnergy, |
---|
135 | HighestKineticEnergy,TotBin) ; |
---|
136 | AtomicNumber = (*theElementTable )[J]->GetZ() ; |
---|
137 | AtomicWeight = (*theElementTable )[J]->GetA() ; |
---|
138 | |
---|
139 | for ( G4int i = 0 ; i < TotBin ; i++) |
---|
140 | { |
---|
141 | LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(i) ; |
---|
142 | Value = ComputeMicroscopicCrossSection(&aParticleType, |
---|
143 | LowEdgeEnergy, |
---|
144 | AtomicNumber,AtomicWeight) ; |
---|
145 | ptrVector->PutValue(i,Value) ; |
---|
146 | } |
---|
147 | |
---|
148 | theCrossSectionTable->insertAt( J , ptrVector ) ; |
---|
149 | } |
---|
150 | |
---|
151 | G4double FixedEnergy = (LowestKineticEnergy + HighestKineticEnergy)/2. ; |
---|
152 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ; |
---|
153 | if (theMeanFreePathTable) { |
---|
154 | theMeanFreePathTable->clearAndDestroy(); |
---|
155 | delete theMeanFreePathTable; |
---|
156 | } |
---|
157 | theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials()); |
---|
158 | |
---|
159 | PartialSumSigma.clearAndDestroy(); |
---|
160 | PartialSumSigma.resize(G4Material::GetNumberOfMaterials()); |
---|
161 | |
---|
162 | |
---|
163 | for (size_t K=0 ; K < G4Material::GetNumberOfMaterials(); K++ ) |
---|
164 | { |
---|
165 | ptrVector = new G4PhysicsLogVector(LowestKineticEnergy, |
---|
166 | HighestKineticEnergy, |
---|
167 | TotBin ) ; |
---|
168 | |
---|
169 | const G4Material* material= (*theMaterialTable)[K]; |
---|
170 | |
---|
171 | for ( G4int i = 0 ; i < TotBin ; i++ ) |
---|
172 | { |
---|
173 | LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ; |
---|
174 | Value = ComputeMeanFreePath( &aParticleType, LowEdgeEnergy, |
---|
175 | material ); |
---|
176 | |
---|
177 | ptrVector->PutValue( i , Value ) ; |
---|
178 | } |
---|
179 | |
---|
180 | theMeanFreePathTable->insertAt( K , ptrVector ); |
---|
181 | |
---|
182 | // Compute the PartialSumSigma table at a given fixed energy |
---|
183 | ComputePartialSumSigma( &aParticleType, FixedEnergy, material); |
---|
184 | } |
---|
185 | |
---|
186 | if (&aParticleType == theMuonPlus) PrintInfoDefinition(); |
---|
187 | } |
---|
188 | |
---|
189 | void G4MuNuclearInteraction::ComputePartialSumSigma( |
---|
190 | const G4ParticleDefinition* ParticleType, |
---|
191 | G4double KineticEnergy, |
---|
192 | const G4Material* aMaterial) |
---|
193 | |
---|
194 | // Build the table of cross section per element. The table is built for MATERIALS. |
---|
195 | // This table is used by DoIt to select randomly an element in the material. |
---|
196 | { |
---|
197 | G4int Imate = aMaterial->GetIndex(); |
---|
198 | G4int NbOfElements = aMaterial->GetNumberOfElements(); |
---|
199 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
200 | const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector(); |
---|
201 | |
---|
202 | PartialSumSigma[Imate] = new G4DataVector(); |
---|
203 | |
---|
204 | G4double SIGMA = 0. ; |
---|
205 | |
---|
206 | for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ ) |
---|
207 | { |
---|
208 | SIGMA += theAtomNumDensityVector[Ielem] * |
---|
209 | ComputeMicroscopicCrossSection( ParticleType, KineticEnergy, |
---|
210 | (*theElementVector)[Ielem]->GetZ(), |
---|
211 | (*theElementVector)[Ielem]->GetA()) ; |
---|
212 | PartialSumSigma[Imate]->push_back(SIGMA); |
---|
213 | } |
---|
214 | } |
---|
215 | |
---|
216 | G4double G4MuNuclearInteraction::ComputeMicroscopicCrossSection( |
---|
217 | const G4ParticleDefinition* ParticleType, |
---|
218 | G4double KineticEnergy, |
---|
219 | G4double AtomicNumber,G4double AtomicWeight) |
---|
220 | { |
---|
221 | static const G4double |
---|
222 | xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 }; |
---|
223 | static const G4double |
---|
224 | wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 }; |
---|
225 | static const G4double ak1=6.9 ; |
---|
226 | static const G4double ak2=1.0 ; |
---|
227 | |
---|
228 | G4double Mass,epmin,epmax,epln,ep,aaa,bbb,hhh,x ; |
---|
229 | G4int kkk ; |
---|
230 | |
---|
231 | Mass = ParticleType->GetPDGMass() ; |
---|
232 | |
---|
233 | G4double CrossSection = 0.0 ; |
---|
234 | |
---|
235 | if ( AtomicNumber < 1. ) return CrossSection; |
---|
236 | if ( KineticEnergy <= CutFixed ) return CrossSection; |
---|
237 | |
---|
238 | epmin = CutFixed ; |
---|
239 | epmax = KineticEnergy + Mass - 0.5*proton_mass_c2 ; |
---|
240 | if (epmax <= epmin ) return CrossSection; // NaN bug correction |
---|
241 | |
---|
242 | aaa = std::log(epmin) ; |
---|
243 | bbb = std::log(epmax) ; |
---|
244 | kkk = G4int((bbb-aaa)/ak1 +ak2) ; |
---|
245 | hhh = (bbb-aaa)/kkk ; |
---|
246 | |
---|
247 | for (G4int l=0 ; l<kkk; l++) |
---|
248 | { |
---|
249 | x = aaa + hhh*l ; |
---|
250 | for (G4int ll=0; ll<8; ll++) |
---|
251 | { |
---|
252 | epln=x+xgi[ll]*hhh ; |
---|
253 | ep = std::exp(epln) ; |
---|
254 | CrossSection += ep*wgi[ll]* ComputeDMicroscopicCrossSection(ParticleType, |
---|
255 | KineticEnergy, |
---|
256 | AtomicNumber,AtomicWeight, |
---|
257 | ep) ; |
---|
258 | } |
---|
259 | } |
---|
260 | CrossSection *= hhh ; |
---|
261 | |
---|
262 | if (CrossSection < 0.) CrossSection = 0.; |
---|
263 | |
---|
264 | return CrossSection; |
---|
265 | } |
---|
266 | |
---|
267 | G4double G4MuNuclearInteraction::ComputeDMicroscopicCrossSection( |
---|
268 | const G4ParticleDefinition* ParticleType, |
---|
269 | G4double KineticEnergy, |
---|
270 | G4double, G4double AtomicWeight, |
---|
271 | G4double epsilon) |
---|
272 | // Calculates the differential (D) microscopic cross section |
---|
273 | // using the cross section formula of R.P. Kokoulin (18/01/98) |
---|
274 | { |
---|
275 | const G4double alam2 = 0.400*GeV*GeV ; |
---|
276 | const G4double alam = 0.632456*GeV ; |
---|
277 | const G4double coeffn = fine_structure_const/pi ; |
---|
278 | |
---|
279 | G4double ep,a,aeff,sigph,v,v1,v2,mass2,up,down ; |
---|
280 | G4double ParticleMass = ParticleType->GetPDGMass() ; |
---|
281 | G4double TotalEnergy = KineticEnergy + ParticleMass ; |
---|
282 | |
---|
283 | G4double DCrossSection = 0. ; |
---|
284 | |
---|
285 | if((epsilon >= TotalEnergy - 0.5*proton_mass_c2) |
---|
286 | || |
---|
287 | (epsilon <= CutFixed)) |
---|
288 | return DCrossSection ; |
---|
289 | |
---|
290 | ep = epsilon/GeV ; |
---|
291 | a = AtomicWeight/(g/mole) ; |
---|
292 | |
---|
293 | aeff = 0.22*a+0.78*std::exp(0.89*std::log(a)) ; //shadowing |
---|
294 | sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn ; //!!!!!!!!!!! |
---|
295 | |
---|
296 | v=epsilon/TotalEnergy ; |
---|
297 | v1 = 1.-v ; |
---|
298 | v2 = v*v ; |
---|
299 | mass2 = ParticleMass*ParticleMass ; |
---|
300 | |
---|
301 | up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1)) ; |
---|
302 | down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam) ; |
---|
303 | |
---|
304 | DCrossSection = coeffn*aeff*sigph/epsilon* |
---|
305 | (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down)) ; |
---|
306 | |
---|
307 | if( DCrossSection < 0.) |
---|
308 | DCrossSection = 0. ; |
---|
309 | |
---|
310 | return DCrossSection ; |
---|
311 | } |
---|
312 | |
---|
313 | void G4MuNuclearInteraction::MakeSamplingTables( |
---|
314 | const G4ParticleDefinition* ParticleType) |
---|
315 | { |
---|
316 | |
---|
317 | G4int nbin; |
---|
318 | G4double AtomicNumber,AtomicWeight,KineticEnergy, |
---|
319 | TotalEnergy,Maxep ; |
---|
320 | |
---|
321 | G4double ParticleMass = ParticleType->GetPDGMass() ; |
---|
322 | |
---|
323 | for (G4int iz=0; iz<nzdat; iz++) |
---|
324 | { |
---|
325 | AtomicNumber = zdat[iz]; |
---|
326 | AtomicWeight = adat[iz]*GramPerMole ; |
---|
327 | |
---|
328 | for (G4int it=0; it<ntdat; it++) |
---|
329 | { |
---|
330 | KineticEnergy = tdat[it]; |
---|
331 | TotalEnergy = KineticEnergy + ParticleMass; |
---|
332 | Maxep = TotalEnergy - 0.5*proton_mass_c2 ; |
---|
333 | |
---|
334 | G4double CrossSection = 0.0 ; |
---|
335 | |
---|
336 | G4double c,y,ymin,ymax,dy,yy,dx,x,ep ; |
---|
337 | |
---|
338 | // calculate the differential cross section |
---|
339 | // numerical integration in |
---|
340 | // log ............... |
---|
341 | c = std::log(Maxep/CutFixed) ; |
---|
342 | ymin = -5. ; |
---|
343 | ymax = 0. ; |
---|
344 | dy = (ymax-ymin)/NBIN ; |
---|
345 | |
---|
346 | nbin=-1; |
---|
347 | |
---|
348 | y = ymin - 0.5*dy ; |
---|
349 | yy = ymin - dy ; |
---|
350 | for (G4int i=0 ; i<NBIN; i++) |
---|
351 | { |
---|
352 | y += dy ; |
---|
353 | x = std::exp(y) ; |
---|
354 | yy += dy ; |
---|
355 | dx = std::exp(yy+dy)-std::exp(yy) ; |
---|
356 | |
---|
357 | ep = CutFixed*std::exp(c*x) ; |
---|
358 | |
---|
359 | CrossSection += ep*dx*ComputeDMicroscopicCrossSection(ParticleType, |
---|
360 | KineticEnergy,AtomicNumber, |
---|
361 | AtomicWeight,ep) ; |
---|
362 | if(nbin<NBIN) |
---|
363 | { |
---|
364 | nbin += 1 ; |
---|
365 | ya[nbin]=y ; |
---|
366 | proba[iz][it][nbin] = CrossSection ; |
---|
367 | } |
---|
368 | } |
---|
369 | ya[NBIN]=0. ; |
---|
370 | |
---|
371 | if(CrossSection > 0.) |
---|
372 | { |
---|
373 | for(G4int ib=0; ib<=nbin; ib++) |
---|
374 | { |
---|
375 | proba[iz][it][ib] /= CrossSection ; |
---|
376 | } |
---|
377 | } |
---|
378 | } |
---|
379 | } |
---|
380 | } |
---|
381 | |
---|
382 | G4VParticleChange* G4MuNuclearInteraction::PostStepDoIt( |
---|
383 | const G4Track& trackData, |
---|
384 | const G4Step& stepData) |
---|
385 | { |
---|
386 | static const G4double Mass=theMuonPlus->GetPDGMass() ; |
---|
387 | static const G4double m0=0.2*GeV ; |
---|
388 | |
---|
389 | aParticleChange.Initialize(trackData); |
---|
390 | G4Material* aMaterial=trackData.GetMaterial() ; |
---|
391 | |
---|
392 | const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle(); |
---|
393 | |
---|
394 | G4double KineticEnergy = aDynamicParticle->GetKineticEnergy(); |
---|
395 | G4ParticleMomentum ParticleDirection = |
---|
396 | aDynamicParticle->GetMomentumDirection(); |
---|
397 | |
---|
398 | // limits of the energy sampling |
---|
399 | G4double TotalEnergy = KineticEnergy + Mass ; |
---|
400 | G4double epmin = CutFixed ; |
---|
401 | G4double epmax = TotalEnergy - 0.5*proton_mass_c2 ; |
---|
402 | // check against insufficient energy |
---|
403 | if (epmax <= epmin ) |
---|
404 | { |
---|
405 | aParticleChange.ProposeMomentumDirection( ParticleDirection ); |
---|
406 | aParticleChange.ProposeEnergy( KineticEnergy ); |
---|
407 | aParticleChange.ProposeLocalEnergyDeposit (0.); |
---|
408 | aParticleChange.SetNumberOfSecondaries(0); |
---|
409 | return G4VDiscreteProcess::PostStepDoIt(trackData,stepData); |
---|
410 | } |
---|
411 | |
---|
412 | // select randomly one element constituing the material |
---|
413 | G4Element* anElement = SelectRandomAtom(aMaterial); |
---|
414 | |
---|
415 | // sample energy of the secondary ( pi0) |
---|
416 | // sampling using tables |
---|
417 | G4double ep,x,y ; |
---|
418 | G4int iy ; |
---|
419 | |
---|
420 | // select sampling table ; |
---|
421 | G4double lnZ = std::log(anElement->GetZ()) ; |
---|
422 | G4double delmin = 1.e10 ; |
---|
423 | G4double del ; |
---|
424 | G4int izz=0,itt=0,NBINminus1 ; |
---|
425 | NBINminus1 = NBIN-1 ; |
---|
426 | for (G4int iz=0; iz<nzdat; iz++) |
---|
427 | { |
---|
428 | del = std::abs(lnZ-std::log(zdat[iz])) ; |
---|
429 | if(del<delmin) |
---|
430 | { |
---|
431 | delmin=del ; |
---|
432 | izz=iz ; |
---|
433 | } |
---|
434 | } |
---|
435 | delmin = 1.e10 ; |
---|
436 | for (G4int it=0; it<ntdat; it++) |
---|
437 | { |
---|
438 | del = std::abs(std::log(KineticEnergy)-std::log(tdat[it])) ; |
---|
439 | if(del<delmin) |
---|
440 | { |
---|
441 | delmin=del; |
---|
442 | itt=it ; |
---|
443 | } |
---|
444 | } |
---|
445 | |
---|
446 | //sample energy transfer according to the sampling table |
---|
447 | |
---|
448 | G4double r = G4UniformRand() ; |
---|
449 | |
---|
450 | iy = -1 ; |
---|
451 | do { |
---|
452 | iy += 1 ; |
---|
453 | } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ; |
---|
454 | |
---|
455 | //sampling is Done uniformly in y in the bin |
---|
456 | if( iy < NBIN ) |
---|
457 | y = ya[iy] + G4UniformRand() * ( ya[iy+1] - ya[iy] ) ; |
---|
458 | else |
---|
459 | y = ya[iy] ; |
---|
460 | |
---|
461 | x = std::exp(y) ; |
---|
462 | ep = epmin*std::exp(x*std::log(epmax/epmin)) ; |
---|
463 | |
---|
464 | // sample scattering angle of mu, but first t should be sampled. |
---|
465 | G4double yy = ep/TotalEnergy ; |
---|
466 | G4double tmin=Mass*Mass*yy*yy/(1.-yy) ; |
---|
467 | G4double tmax=2.*proton_mass_c2*ep ; |
---|
468 | G4double t1,t2,t,w1,w2,w3,y1,y2,y3,rej ; |
---|
469 | if(m0<ep) |
---|
470 | { |
---|
471 | t1=m0*m0; |
---|
472 | t2=ep*ep; |
---|
473 | } |
---|
474 | else |
---|
475 | { |
---|
476 | t1=ep*ep; |
---|
477 | t2=m0*m0; |
---|
478 | } |
---|
479 | |
---|
480 | w1=tmax*t1; |
---|
481 | w2=tmax+t1 ; |
---|
482 | w3=tmax*(tmin+t1)/(tmin*w2); |
---|
483 | y1=1.-yy; |
---|
484 | y2=0.5*yy*yy; |
---|
485 | y3=y1+y2; |
---|
486 | |
---|
487 | // now the sampling of t |
---|
488 | G4int ntry=0; |
---|
489 | do |
---|
490 | { |
---|
491 | ntry += 1 ; |
---|
492 | t=w1/(w2*std::exp(G4UniformRand()*std::log(w3))-tmax) ; |
---|
493 | rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2)); |
---|
494 | } while (G4UniformRand() > rej) ; |
---|
495 | |
---|
496 | // compute angle from t |
---|
497 | G4double sinth2,theta ; // sinth2 = std::sin(theta/2)*std::sin(theta/2) ! |
---|
498 | sinth2 = 0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin) ; |
---|
499 | theta = std::acos(1.-2.*sinth2) ; |
---|
500 | |
---|
501 | G4double phi=twopi*G4UniformRand() ; |
---|
502 | G4double sinth=std::sin(theta) ; |
---|
503 | G4double dirx=sinth*std::cos(phi) , diry=sinth*std::sin(phi) , dirz=std::cos(theta); |
---|
504 | |
---|
505 | G4ThreeVector finalDirection(dirx,diry,dirz) ; |
---|
506 | finalDirection.rotateUz(ParticleDirection) ; |
---|
507 | |
---|
508 | G4double NewKinEnergy = KineticEnergy - ep ; |
---|
509 | G4double finalMomentum=std::sqrt(NewKinEnergy* |
---|
510 | (NewKinEnergy+2.*Mass)) ; |
---|
511 | |
---|
512 | G4double Ef=NewKinEnergy+Mass ; |
---|
513 | G4double initMomentum=std::sqrt(KineticEnergy*(TotalEnergy+Mass)) ; |
---|
514 | |
---|
515 | // G4double Q2=2.*(TotalEnergy*Ef-initMomentum*finalMomentum*std::cos(theta)-Mass*Mass) ; |
---|
516 | |
---|
517 | aParticleChange.ProposeMomentumDirection( finalDirection ); |
---|
518 | aParticleChange.ProposeEnergy( NewKinEnergy ); |
---|
519 | |
---|
520 | G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy); |
---|
521 | G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef); |
---|
522 | G4LorentzVector momentumTransfere = primaryMomentum-fsMomentum; |
---|
523 | |
---|
524 | G4DynamicParticle* aGamma = |
---|
525 | new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfere); |
---|
526 | G4Track gammaTrack(aGamma, trackData.GetGlobalTime(), trackData.GetPosition() ); |
---|
527 | gammaTrack.SetStep(trackData.GetStep()); |
---|
528 | G4Nucleus theTarget(aMaterial); |
---|
529 | |
---|
530 | G4VParticleChange* aHadronicFS = theHadronicVertex.ApplyYourself(theTarget, gammaTrack); |
---|
531 | //delete aGamma; |
---|
532 | |
---|
533 | G4int numSecondaries = aHadronicFS->GetNumberOfSecondaries(); |
---|
534 | aParticleChange.SetNumberOfSecondaries(numSecondaries); |
---|
535 | |
---|
536 | // G4ParticleMomentum secondaryMomentum = G4ThreeVector(0.,0.,0.); |
---|
537 | for(G4int iSec=0; iSec<numSecondaries; iSec++) |
---|
538 | { |
---|
539 | //secondaryMomentum |
---|
540 | // = secondaryMomentum + aHadronicFS->GetSecondary(iSec)->GetMomentum(); |
---|
541 | aParticleChange.AddSecondary(aHadronicFS->GetSecondary(iSec)); |
---|
542 | } |
---|
543 | aHadronicFS->Clear(); |
---|
544 | |
---|
545 | return G4VDiscreteProcess::PostStepDoIt(trackData,stepData); |
---|
546 | } |
---|
547 | |
---|
548 | G4Element* G4MuNuclearInteraction::SelectRandomAtom(G4Material* aMaterial) const |
---|
549 | { |
---|
550 | // select randomly 1 element within the material |
---|
551 | |
---|
552 | G4int Index = aMaterial->GetIndex(); |
---|
553 | G4int NumberOfElements = aMaterial->GetNumberOfElements(); |
---|
554 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
555 | if(1 == NumberOfElements) return ((*theElementVector)[0]); |
---|
556 | |
---|
557 | G4double rval = G4UniformRand()*((*PartialSumSigma[Index])[NumberOfElements-1]); |
---|
558 | for ( G4int i=0; i < NumberOfElements; i++ ) { |
---|
559 | if (rval <= (*PartialSumSigma[Index])[i]) return ((*theElementVector)[i]); |
---|
560 | } |
---|
561 | G4cout << "G4MuNuclearInteraction WARNING !!! no element selected for '" |
---|
562 | << aMaterial->GetName() |
---|
563 | << " 1st element returned." << G4endl; |
---|
564 | return ((*theElementVector)[0]); |
---|
565 | } |
---|
566 | |
---|
567 | void G4MuNuclearInteraction::PrintInfoDefinition() |
---|
568 | { |
---|
569 | G4String comments = "cross sections from R. Kokoulin \n "; |
---|
570 | comments += " Good description up to 1000 PeV."; |
---|
571 | |
---|
572 | G4cout << G4endl << GetProcessName() << ": " << comments |
---|
573 | << "\n PhysicsTables from " << G4BestUnit(LowestKineticEnergy, |
---|
574 | "Energy") |
---|
575 | << " to " << G4BestUnit(HighestKineticEnergy,"Energy") |
---|
576 | << " in " << TotBin << " bins. \n"; |
---|
577 | |
---|
578 | G4cout << G4endl; |
---|
579 | } |
---|
580 | |
---|
581 | G4double G4MuNuclearInteraction::GetMeanFreePath(const G4Track& trackData, |
---|
582 | G4double, |
---|
583 | G4ForceCondition* condition) |
---|
584 | { |
---|
585 | const G4DynamicParticle* aDynamicParticle; |
---|
586 | G4Material* aMaterial; |
---|
587 | G4double MeanFreePath; |
---|
588 | G4bool isOutRange ; |
---|
589 | |
---|
590 | *condition = NotForced ; |
---|
591 | |
---|
592 | aDynamicParticle = trackData.GetDynamicParticle(); |
---|
593 | aMaterial = trackData.GetMaterial(); |
---|
594 | |
---|
595 | G4double KineticEnergy = aDynamicParticle->GetKineticEnergy(); |
---|
596 | |
---|
597 | if (KineticEnergy < LowestKineticEnergy) |
---|
598 | MeanFreePath = DBL_MAX ; |
---|
599 | else { |
---|
600 | if (KineticEnergy > HighestKineticEnergy) |
---|
601 | KineticEnergy = 0.99*HighestKineticEnergy ; |
---|
602 | MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())-> |
---|
603 | GetValue( KineticEnergy, isOutRange ); |
---|
604 | } |
---|
605 | return MeanFreePath; |
---|
606 | } |
---|
607 | |
---|
608 | G4double G4MuNuclearInteraction::ComputeMeanFreePath( |
---|
609 | const G4ParticleDefinition* ParticleType, |
---|
610 | G4double KineticEnergy, |
---|
611 | const G4Material* aMaterial) |
---|
612 | { |
---|
613 | const G4ElementVector* theElementVector = aMaterial->GetElementVector() ; |
---|
614 | const G4double* theAtomNumDensityVector = |
---|
615 | aMaterial->GetAtomicNumDensityVector(); |
---|
616 | |
---|
617 | G4double SIGMA = 0 ; |
---|
618 | |
---|
619 | for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ ) |
---|
620 | { |
---|
621 | SIGMA += theAtomNumDensityVector[i] * |
---|
622 | ComputeMicroscopicCrossSection( ParticleType, KineticEnergy, |
---|
623 | (*theElementVector)[i]->GetZ(), |
---|
624 | (*theElementVector)[i]->GetA()) ; |
---|
625 | } |
---|
626 | |
---|
627 | return SIGMA<=0.0 ? DBL_MAX : 1./SIGMA ; |
---|
628 | } |
---|
629 | |
---|
630 | |
---|