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