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 | // GEANT 4 class header file |
---|
28 | // |
---|
29 | // History: |
---|
30 | // 01 August 2007 P.Gumplinger |
---|
31 | // Reference: TRIUMF TWIST Technotes TN-55: |
---|
32 | // Pierre Depommier - "Radiative MuonDecay" |
---|
33 | // |
---|
34 | // ------------------------------------------------------------ |
---|
35 | // |
---|
36 | // |
---|
37 | // |
---|
38 | |
---|
39 | #include "G4MuonRadiativeDecayChannelWithSpin.hh" |
---|
40 | |
---|
41 | #include "Randomize.hh" |
---|
42 | #include "G4DecayProducts.hh" |
---|
43 | #include "G4LorentzVector.hh" |
---|
44 | |
---|
45 | G4MuonRadiativeDecayChannelWithSpin:: |
---|
46 | G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName, |
---|
47 | G4double theBR) |
---|
48 | : G4VDecayChannel("Radiative Muon Decay",1) |
---|
49 | { |
---|
50 | // set names for daughter particles |
---|
51 | if (theParentName == "mu+") { |
---|
52 | SetBR(theBR); |
---|
53 | SetParent("mu+"); |
---|
54 | SetNumberOfDaughters(4); |
---|
55 | SetDaughter(0, "e+"); |
---|
56 | SetDaughter(1, "gamma"); |
---|
57 | SetDaughter(2, "nu_e"); |
---|
58 | SetDaughter(3, "anti_nu_mu"); |
---|
59 | } else if (theParentName == "mu-") { |
---|
60 | SetBR(theBR); |
---|
61 | SetParent("mu-"); |
---|
62 | SetNumberOfDaughters(4); |
---|
63 | SetDaughter(0, "e-"); |
---|
64 | SetDaughter(1, "gamma"); |
---|
65 | SetDaughter(2, "anti_nu_e"); |
---|
66 | SetDaughter(3, "nu_mu"); |
---|
67 | } else { |
---|
68 | #ifdef G4VERBOSE |
---|
69 | if (GetVerboseLevel()>0) { |
---|
70 | G4cout << "G4RadiativeMuonDecayChannel:: constructor :"; |
---|
71 | G4cout << " parent particle is not muon but "; |
---|
72 | G4cout << theParentName << G4endl; |
---|
73 | } |
---|
74 | #endif |
---|
75 | } |
---|
76 | EMMU = 0.*MeV; |
---|
77 | EMASS = 0.*MeV; |
---|
78 | } |
---|
79 | |
---|
80 | G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin() |
---|
81 | { |
---|
82 | } |
---|
83 | |
---|
84 | G4DecayProducts *G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double) |
---|
85 | { |
---|
86 | |
---|
87 | #ifdef G4VERBOSE |
---|
88 | if (GetVerboseLevel()>1) |
---|
89 | G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt "; |
---|
90 | #endif |
---|
91 | |
---|
92 | if (parent == 0) FillParent(); |
---|
93 | if (daughters == 0) FillDaughters(); |
---|
94 | |
---|
95 | // parent mass |
---|
96 | G4double parentmass = parent->GetPDGMass(); |
---|
97 | |
---|
98 | EMMU = parentmass; |
---|
99 | |
---|
100 | //daughters'mass |
---|
101 | G4double daughtermass[4]; |
---|
102 | G4double sumofdaughtermass = 0.0; |
---|
103 | for (G4int index=0; index<4; index++){ |
---|
104 | daughtermass[index] = daughters[index]->GetPDGMass(); |
---|
105 | sumofdaughtermass += daughtermass[index]; |
---|
106 | } |
---|
107 | |
---|
108 | EMASS = daughtermass[0]; |
---|
109 | |
---|
110 | //create parent G4DynamicParticle at rest |
---|
111 | G4ThreeVector dummy; |
---|
112 | G4DynamicParticle * parentparticle = |
---|
113 | new G4DynamicParticle( parent, dummy, 0.0); |
---|
114 | //create G4Decayproducts |
---|
115 | G4DecayProducts *products = new G4DecayProducts(*parentparticle); |
---|
116 | delete parentparticle; |
---|
117 | |
---|
118 | G4int i = 0; |
---|
119 | |
---|
120 | G4double eps = EMASS/EMMU; |
---|
121 | |
---|
122 | G4double som0, Qsqr, x, y, xx, yy, zz; |
---|
123 | G4double cthetaE, cthetaG, cthetaGE, phiE, phiG; |
---|
124 | |
---|
125 | do { |
---|
126 | |
---|
127 | // leap1: |
---|
128 | |
---|
129 | i++; |
---|
130 | |
---|
131 | // leap2: |
---|
132 | |
---|
133 | do { |
---|
134 | // |
---|
135 | //-------------------------------------------------------------------------- |
---|
136 | // Build two vectors of random length and random direction, for the |
---|
137 | // positron and the photon. |
---|
138 | // x/y is the length of the vector, xx, yy and zz the components, |
---|
139 | // phi is the azimutal angle, theta the polar angle. |
---|
140 | //-------------------------------------------------------------------------- |
---|
141 | // |
---|
142 | // For the positron |
---|
143 | // |
---|
144 | x = G4UniformRand(); |
---|
145 | |
---|
146 | rn3dim(xx,yy,zz,x); |
---|
147 | |
---|
148 | if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){ |
---|
149 | G4cout << "Norm of x not correct" << G4endl; |
---|
150 | } |
---|
151 | |
---|
152 | phiE = atan4(xx,yy); |
---|
153 | cthetaE = zz/x; |
---|
154 | G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x; |
---|
155 | // |
---|
156 | // What you get: |
---|
157 | // |
---|
158 | // x = positron energy |
---|
159 | // phiE = azimutal angle of positron momentum |
---|
160 | // cthetaE = cosine of polar angle of positron momentum |
---|
161 | // sthetaE = sine of polar angle of positron momentum |
---|
162 | // |
---|
163 | //// G4cout << " x, xx, yy, zz " << x << " " << xx << " " |
---|
164 | //// << yy << " " << zz << G4endl; |
---|
165 | //// G4cout << " phiE, cthetaE, sthetaE " << phiE << " " |
---|
166 | //// << cthetaE << " " |
---|
167 | //// << sthetaE << " " << G4endl; |
---|
168 | // |
---|
169 | //----------------------------------------------------------------------- |
---|
170 | // |
---|
171 | // For the photon |
---|
172 | // |
---|
173 | y = G4UniformRand(); |
---|
174 | |
---|
175 | rn3dim(xx,yy,zz,y); |
---|
176 | |
---|
177 | if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){ |
---|
178 | G4cout << " Norm of y not correct " << G4endl; |
---|
179 | } |
---|
180 | |
---|
181 | phiG = atan4(xx,yy); |
---|
182 | cthetaG = zz/y; |
---|
183 | G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y; |
---|
184 | // |
---|
185 | // What you get: |
---|
186 | // |
---|
187 | // y = photon energy |
---|
188 | // phiG = azimutal angle of photon momentum |
---|
189 | // cthetaG = cosine of polar angle of photon momentum |
---|
190 | // sthetaG = sine of polar angle of photon momentum |
---|
191 | // |
---|
192 | //// G4cout << " y, xx, yy, zz " << y << " " << xx << " " |
---|
193 | //// << yy << " " << zz << G4endl; |
---|
194 | //// G4cout << " phiG, cthetaG, sthetaG " << phiG << " " |
---|
195 | //// << cthetaG << " " |
---|
196 | //// << sthetaG << " " << G4endl; |
---|
197 | // |
---|
198 | //----------------------------------------------------------------------- |
---|
199 | // |
---|
200 | // Maybe certain restrictions on the kinematical variables: |
---|
201 | // |
---|
202 | //// if (cthetaE > 0.01)goto leap2; |
---|
203 | //// if (cthetaG > 0.01)goto leap2; |
---|
204 | //// if (std::fabs(x-0.5) > 0.5 )goto leap2; |
---|
205 | //// if (std::fabs(y-0.5) > 0.5 )goto leap2; |
---|
206 | // |
---|
207 | //----------------------------------------------------------------------- |
---|
208 | // |
---|
209 | // Calculate the angle between positron and photon (cosine) |
---|
210 | // |
---|
211 | cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG); |
---|
212 | // |
---|
213 | //// G4cout << x << " " << cthetaE << " " << sthetaE << " " |
---|
214 | //// << y << " " << cthetaG << " " << sthetaG << " " |
---|
215 | //// << cthetaGE |
---|
216 | // |
---|
217 | //----------------------------------------------------------------------- |
---|
218 | // |
---|
219 | G4double term0 = eps*eps; |
---|
220 | G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps; |
---|
221 | G4double beta = std::sqrt( x*((1.0-eps)*(1.0-eps))* |
---|
222 | (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1; |
---|
223 | G4double delta = 1.0-beta*cthetaGE; |
---|
224 | |
---|
225 | G4double term3 = y*(1.0-(eps*eps)); |
---|
226 | G4double term6 = term1*delta*term3; |
---|
227 | |
---|
228 | Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps)); |
---|
229 | // |
---|
230 | //----------------------------------------------------------------------- |
---|
231 | // |
---|
232 | // Check the kinematics. |
---|
233 | // |
---|
234 | } while ( Qsqr<0.0 || Qsqr>1.0 ); |
---|
235 | // |
---|
236 | //// G4cout << x << " " << y << " " << beta << " " << Qsqr << G4endl; |
---|
237 | // |
---|
238 | // Do the calculation for -1 muon polarization (i.e. mu+) |
---|
239 | // |
---|
240 | G4double Pmu = -1.0; |
---|
241 | if (GetParentName() == "mu-")Pmu = +1.0; |
---|
242 | // |
---|
243 | // and for Fronsdal |
---|
244 | // |
---|
245 | //----------------------------------------------------------------------- |
---|
246 | // |
---|
247 | som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE); |
---|
248 | // |
---|
249 | //// if(som0<0.0){ |
---|
250 | //// G4cout << " som0 < 0 in Fronsdal " << som0 |
---|
251 | //// << " at event " << i << G4endl; |
---|
252 | //// G4cout << Pmu << " " << x << " " << y << " " |
---|
253 | //// << cthetaE << " " << cthetaG << " " |
---|
254 | //// << cthetaGE << " " << som0 << G4endl; |
---|
255 | //// } |
---|
256 | // |
---|
257 | //----------------------------------------------------------------------- |
---|
258 | // |
---|
259 | //// G4cout << x << " " << y << " " << som0 << G4endl; |
---|
260 | // |
---|
261 | //---------------------------------------------------------------------- |
---|
262 | // |
---|
263 | // Sample the decay rate |
---|
264 | // |
---|
265 | |
---|
266 | } while (G4UniformRand()*250000.0 > som0); |
---|
267 | |
---|
268 | /// if(i<10000000)goto leap1: |
---|
269 | // |
---|
270 | //----------------------------------------------------------------------- |
---|
271 | // |
---|
272 | G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps); |
---|
273 | G4double G = EMMU/2.*y*(1.-eps*eps); |
---|
274 | // |
---|
275 | //----------------------------------------------------------------------- |
---|
276 | // |
---|
277 | |
---|
278 | if(E < EMASS) E = EMASS; |
---|
279 | |
---|
280 | // calculate daughter momentum |
---|
281 | G4double daughtermomentum[4]; |
---|
282 | |
---|
283 | daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS); |
---|
284 | |
---|
285 | G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE); |
---|
286 | G4double cphiE = std::cos(phiE); |
---|
287 | G4double sphiE = std::sin(phiE); |
---|
288 | |
---|
289 | //Coordinates of the decay positron with respect to the muon spin |
---|
290 | |
---|
291 | G4double px = sthetaE*cphiE; |
---|
292 | G4double py = sthetaE*sphiE; |
---|
293 | G4double pz = cthetaE; |
---|
294 | |
---|
295 | G4ThreeVector direction0(px,py,pz); |
---|
296 | |
---|
297 | direction0.rotateUz(parent_polarization); |
---|
298 | |
---|
299 | G4DynamicParticle * daughterparticle0 |
---|
300 | = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0); |
---|
301 | |
---|
302 | products->PushProducts(daughterparticle0); |
---|
303 | |
---|
304 | daughtermomentum[1] = G; |
---|
305 | |
---|
306 | G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG); |
---|
307 | G4double cphiG = std::cos(phiG); |
---|
308 | G4double sphiG = std::sin(phiG); |
---|
309 | |
---|
310 | //Coordinates of the decay gamma with respect to the muon spin |
---|
311 | |
---|
312 | px = sthetaG*cphiG; |
---|
313 | py = sthetaG*sphiG; |
---|
314 | pz = cthetaG; |
---|
315 | |
---|
316 | G4ThreeVector direction1(px,py,pz); |
---|
317 | |
---|
318 | direction1.rotateUz(parent_polarization); |
---|
319 | |
---|
320 | G4DynamicParticle * daughterparticle1 |
---|
321 | = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1); |
---|
322 | |
---|
323 | products->PushProducts(daughterparticle1); |
---|
324 | |
---|
325 | // daughter 3 ,4 (neutrinos) |
---|
326 | // create neutrinos in the C.M frame of two neutrinos |
---|
327 | |
---|
328 | G4double energy2 = parentmass*(1.0 - (x+y)/2.0); |
---|
329 | |
---|
330 | G4double vmass = std::sqrt((energy2- |
---|
331 | (daughtermomentum[0]+daughtermomentum[1]))* |
---|
332 | (energy2+ |
---|
333 | (daughtermomentum[0]+daughtermomentum[1]))); |
---|
334 | G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2; |
---|
335 | beta = -1.0 * std::min(beta,0.99); |
---|
336 | |
---|
337 | G4double costhetan = 2.*G4UniformRand()-1.0; |
---|
338 | G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); |
---|
339 | G4double phin = twopi*G4UniformRand()*rad; |
---|
340 | G4double sinphin = std::sin(phin); |
---|
341 | G4double cosphin = std::cos(phin); |
---|
342 | |
---|
343 | G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan); |
---|
344 | |
---|
345 | G4DynamicParticle * daughterparticle2 |
---|
346 | = new G4DynamicParticle( daughters[2], direction2*(vmass/2.)); |
---|
347 | G4DynamicParticle * daughterparticle3 |
---|
348 | = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.)); |
---|
349 | |
---|
350 | // boost to the muon rest frame |
---|
351 | |
---|
352 | G4ThreeVector direction34(direction0.x()+direction1.x(), |
---|
353 | direction0.y()+direction1.y(), |
---|
354 | direction0.z()+direction1.z()); |
---|
355 | direction34 = direction34.unit(); |
---|
356 | |
---|
357 | G4LorentzVector p4 = daughterparticle2->Get4Momentum(); |
---|
358 | p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta); |
---|
359 | daughterparticle2->Set4Momentum(p4); |
---|
360 | |
---|
361 | p4 = daughterparticle3->Get4Momentum(); |
---|
362 | p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta); |
---|
363 | daughterparticle3->Set4Momentum(p4); |
---|
364 | |
---|
365 | products->PushProducts(daughterparticle2); |
---|
366 | products->PushProducts(daughterparticle3); |
---|
367 | |
---|
368 | daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); |
---|
369 | daughtermomentum[3] = daughterparticle3->GetTotalMomentum(); |
---|
370 | |
---|
371 | // output message |
---|
372 | #ifdef G4VERBOSE |
---|
373 | if (GetVerboseLevel()>1) { |
---|
374 | G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt "; |
---|
375 | G4cout << " create decay products in rest frame " <<G4endl; |
---|
376 | products->DumpInfo(); |
---|
377 | } |
---|
378 | #endif |
---|
379 | return products; |
---|
380 | } |
---|
381 | |
---|
382 | G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu, |
---|
383 | G4double x, |
---|
384 | G4double y, |
---|
385 | G4double cthetaE, |
---|
386 | G4double cthetaG, |
---|
387 | G4double cthetaGE) |
---|
388 | { |
---|
389 | G4double mu = 105.65; |
---|
390 | G4double me = 0.511; |
---|
391 | G4double rho = 0.75; |
---|
392 | G4double del = 0.75; |
---|
393 | G4double eps = 0.0; |
---|
394 | G4double kap = 0.0; |
---|
395 | G4double ksi = 1.0; |
---|
396 | |
---|
397 | G4double delta = 1-cthetaGE; |
---|
398 | |
---|
399 | // Calculation of the functions f(x,y) |
---|
400 | |
---|
401 | G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y) |
---|
402 | +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x)); |
---|
403 | G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y)) |
---|
404 | -2.0*(x*x)*(1.0-y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y)); |
---|
405 | G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y)) |
---|
406 | -(x*x*x)*y*(4.0+3.0*y)); |
---|
407 | G4double f2s = 1.5*((x*x*x)*(y*y)*(2.0+y)); |
---|
408 | |
---|
409 | G4double f_1se = 12.0*(x*y*(1.0-x)+(x*x)*(2.0-3.0*y) |
---|
410 | -2.0*(x*x*x)); |
---|
411 | G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y)) |
---|
412 | +(x*x*x)*(2.0+3.0*y)); |
---|
413 | G4double f1se = -3.0*(x*x*x)*y*(2.0+y); |
---|
414 | G4double f2se = 0.0; |
---|
415 | |
---|
416 | G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y) |
---|
417 | -(x*x)*y); |
---|
418 | G4double f0sg = 6.0*(-x*(y*y)*(2.0-3.0*y)-(x*x)*y*(1.0-4.0*y) |
---|
419 | +(x*x*x)*y); |
---|
420 | G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y) |
---|
421 | -2.0*(x*x*x)*(y*y)); |
---|
422 | G4double f2sg = 1.5*(x*x*x)*(y*y*y); |
---|
423 | |
---|
424 | G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y) |
---|
425 | +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x)); |
---|
426 | G4double f0v = 8.0*(-x*y*(3.0-y-(y*y))-(x*x)*(3.0-y-4.0*(y*y)) |
---|
427 | +2.0*(x*x*x)*(1.0+2.0*y)); |
---|
428 | G4double f1v = 2.0*((x*x)*y*(6.0-5.0*y-2.0*(y*y)) |
---|
429 | -2.0*(x*x*x)*y*(4.0+3.0*y)); |
---|
430 | G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y); |
---|
431 | |
---|
432 | G4double f_1ve = 8.0*(x*y*(1.0-2.0*y) |
---|
433 | +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x)); |
---|
434 | G4double f0ve = 4.0*(-(x*x)*(2.0-3.0*y-4.0*(y*y)) |
---|
435 | +2.0*(x*x*x)*(2.0+3.0*y)); |
---|
436 | G4double f1ve = -4.0*(x*x*x)*y*(2.0+y); |
---|
437 | G4double f2ve = 0.0; |
---|
438 | |
---|
439 | G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y) |
---|
440 | -2.0*(x*x)*y); |
---|
441 | G4double f0vg = 4.0*(2.0*x*(y*y)*(1.0+y)-(x*x)*y*(1.0-4.0*y) |
---|
442 | +2.0*(x*x*x)*y); |
---|
443 | G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y) |
---|
444 | -4.0*(x*x*x)*(y*y)); |
---|
445 | G4double f2vg = 2.0*(x*x*x)*(y*y*y); |
---|
446 | |
---|
447 | G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y) |
---|
448 | +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x)); |
---|
449 | G4double f0t = 4.0*(-x*y*(6.0+(y*y)) |
---|
450 | -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y)); |
---|
451 | G4double f1t = 2.0*((x*x)*y*(6.0-5.0*y+(y*y)) |
---|
452 | -(x*x*x)*y*(4.0+3.0*y)); |
---|
453 | G4double f2t = (x*x*x)*(y*y)*(2.0+y); |
---|
454 | |
---|
455 | G4double f_1te = -8.0*(x*y*(1.0+3.0*y)+(x*x)*(2.0+3.0*y) |
---|
456 | +2.0*(x*x*x)); |
---|
457 | G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y)) |
---|
458 | +(x*x*x)*(2.0+3.0*y)); |
---|
459 | G4double f1te = -2.0*(x*x*x)*y*(2.0+y); |
---|
460 | G4double f2te = 0.0; |
---|
461 | |
---|
462 | G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y); |
---|
463 | G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y) |
---|
464 | +(x*x*x)*y); |
---|
465 | G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y); |
---|
466 | G4double f2tg = (x*x*x)*(y*y*y); |
---|
467 | |
---|
468 | G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x)); |
---|
469 | term = 1.0/term; |
---|
470 | |
---|
471 | G4double ns = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s; |
---|
472 | G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v; |
---|
473 | G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t; |
---|
474 | |
---|
475 | G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se; |
---|
476 | G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve; |
---|
477 | G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te; |
---|
478 | |
---|
479 | G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg; |
---|
480 | G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg; |
---|
481 | G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg; |
---|
482 | |
---|
483 | G4double term1 = nv; |
---|
484 | G4double term2 = 2.0*ns-nv-nt; |
---|
485 | G4double term3 = 2.0*ns-2.0*nv+nt; |
---|
486 | |
---|
487 | G4double term1e = 1.0/3.0*(1.0-4.0/3.0*del); |
---|
488 | G4double term2e = 2.0*nse+5.0*nve-nte; |
---|
489 | G4double term3e = 2.0*nse-2.0*nve+nte; |
---|
490 | |
---|
491 | G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del); |
---|
492 | G4double term2g = 2.0*nsg+5.0*nvg-ntg; |
---|
493 | G4double term3g = 2.0*nsg-2.0*nvg+ntg; |
---|
494 | |
---|
495 | G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3; |
---|
496 | G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e) |
---|
497 | +cthetaG*(nvg-term1g*term2g+kap*term3g)); |
---|
498 | G4double som0 = som00+som01; |
---|
499 | |
---|
500 | // G4cout << x << " " << y << " " << som00 << " " |
---|
501 | // << som01 << " " << som0 << G4endl; |
---|
502 | |
---|
503 | return som0; |
---|
504 | } |
---|