00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4MuonRadiativeDecayChannelWithSpin.hh"
00046
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "Randomize.hh"
00050 #include "G4DecayProducts.hh"
00051 #include "G4LorentzVector.hh"
00052
00053 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin()
00054 : G4VDecayChannel(),
00055 parent_polarization()
00056 {
00057 }
00058
00059 G4MuonRadiativeDecayChannelWithSpin::
00060 G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName,
00061 G4double theBR)
00062 : G4VDecayChannel("Radiative Muon Decay",1),
00063 parent_polarization()
00064 {
00065
00066 if (theParentName == "mu+") {
00067 SetBR(theBR);
00068 SetParent("mu+");
00069 SetNumberOfDaughters(4);
00070 SetDaughter(0, "e+");
00071 SetDaughter(1, "gamma");
00072 SetDaughter(2, "nu_e");
00073 SetDaughter(3, "anti_nu_mu");
00074 } else if (theParentName == "mu-") {
00075 SetBR(theBR);
00076 SetParent("mu-");
00077 SetNumberOfDaughters(4);
00078 SetDaughter(0, "e-");
00079 SetDaughter(1, "gamma");
00080 SetDaughter(2, "anti_nu_e");
00081 SetDaughter(3, "nu_mu");
00082 } else {
00083 #ifdef G4VERBOSE
00084 if (GetVerboseLevel()>0) {
00085 G4cout << "G4RadiativeMuonDecayChannel:: constructor :";
00086 G4cout << " parent particle is not muon but ";
00087 G4cout << theParentName << G4endl;
00088 }
00089 #endif
00090 }
00091 }
00092
00093 G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin()
00094 {
00095 }
00096
00097 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(const G4MuonRadiativeDecayChannelWithSpin &right):
00098 G4VDecayChannel(right)
00099 {
00100 parent_polarization = right.parent_polarization;
00101 }
00102
00103 G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin & right)
00104 {
00105 if (this != &right) {
00106 kinematics_name = right.kinematics_name;
00107 verboseLevel = right.verboseLevel;
00108 rbranch = right.rbranch;
00109
00110
00111 parent_name = new G4String(*right.parent_name);
00112
00113
00114 ClearDaughtersName();
00115
00116
00117 numberOfDaughters = right.numberOfDaughters;
00118 if ( numberOfDaughters >0 ) {
00119 if (daughters_name !=0) ClearDaughtersName();
00120 daughters_name = new G4String*[numberOfDaughters];
00121
00122 for (G4int index=0; index < numberOfDaughters; index++) {
00123 daughters_name[index] = new G4String(*right.daughters_name[index]);
00124 }
00125 }
00126 parent_polarization = right.parent_polarization;
00127 }
00128 return *this;
00129 }
00130
00131
00132 G4DecayProducts *G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double)
00133 {
00134
00135 #ifdef G4VERBOSE
00136 if (GetVerboseLevel()>1)
00137 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
00138 #endif
00139
00140 if (parent == 0) FillParent();
00141 if (daughters == 0) FillDaughters();
00142
00143
00144 G4double parentmass = parent->GetPDGMass();
00145
00146 G4double EMMU = parentmass;
00147
00148
00149 G4double daughtermass[4];
00150 G4double sumofdaughtermass = 0.0;
00151 for (G4int index=0; index<4; index++){
00152 daughtermass[index] = daughters[index]->GetPDGMass();
00153 sumofdaughtermass += daughtermass[index];
00154 }
00155
00156 G4double EMASS = daughtermass[0];
00157
00158
00159 G4ThreeVector dummy;
00160 G4DynamicParticle * parentparticle =
00161 new G4DynamicParticle( parent, dummy, 0.0);
00162
00163 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00164 delete parentparticle;
00165
00166 G4int i = 0;
00167
00168 G4double eps = EMASS/EMMU;
00169
00170 G4double som0, Qsqr, x, y, xx, yy, zz;
00171 G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
00172
00173 do {
00174
00175
00176
00177 i++;
00178
00179
00180
00181 do {
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 x = G4UniformRand();
00193
00194 rn3dim(xx,yy,zz,x);
00195
00196 if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){
00197 G4cout << "Norm of x not correct" << G4endl;
00198 }
00199
00200 phiE = atan4(xx,yy);
00201 cthetaE = zz/x;
00202 G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x;
00203
00204
00205
00206
00207
00208
00209
00210
00216
00217
00218
00219
00220
00221 y = G4UniformRand();
00222
00223 rn3dim(xx,yy,zz,y);
00224
00225 if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){
00226 G4cout << " Norm of y not correct " << G4endl;
00227 }
00228
00229 phiG = atan4(xx,yy);
00230 cthetaG = zz/y;
00231 G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y;
00232
00233
00234
00235
00236
00237
00238
00239
00245
00246
00247
00248
00249
00254
00255
00256
00257
00258
00259 cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG);
00260
00264
00265
00266
00267 G4double term0 = eps*eps;
00268 G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps;
00269 G4double beta = std::sqrt( x*((1.0-eps)*(1.0-eps))*
00270 (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1;
00271 G4double delta = 1.0-beta*cthetaGE;
00272
00273 G4double term3 = y*(1.0-(eps*eps));
00274 G4double term6 = term1*delta*term3;
00275
00276 Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
00277
00278
00279
00280
00281
00282 } while ( Qsqr<0.0 || Qsqr>1.0 );
00283
00285
00286
00287
00288 G4double Pmu = -1.0;
00289 if (GetParentName() == "mu-")Pmu = +1.0;
00290
00291
00292
00293
00294
00295 som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
00296
00304
00305
00306
00308
00309
00310
00311
00312
00313
00314 } while (G4UniformRand()*177.0 > som0);
00315
00317
00318
00319
00320 G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps);
00321 G4double G = EMMU/2.*y*(1.-eps*eps);
00322
00323
00324
00325
00326 if(E < EMASS) E = EMASS;
00327
00328
00329 G4double daughtermomentum[4];
00330
00331 daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
00332
00333 G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
00334 G4double cphiE = std::cos(phiE);
00335 G4double sphiE = std::sin(phiE);
00336
00337
00338
00339 G4double px = sthetaE*cphiE;
00340 G4double py = sthetaE*sphiE;
00341 G4double pz = cthetaE;
00342
00343 G4ThreeVector direction0(px,py,pz);
00344
00345 direction0.rotateUz(parent_polarization);
00346
00347 G4DynamicParticle * daughterparticle0
00348 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
00349
00350 products->PushProducts(daughterparticle0);
00351
00352 daughtermomentum[1] = G;
00353
00354 G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG);
00355 G4double cphiG = std::cos(phiG);
00356 G4double sphiG = std::sin(phiG);
00357
00358
00359
00360 px = sthetaG*cphiG;
00361 py = sthetaG*sphiG;
00362 pz = cthetaG;
00363
00364 G4ThreeVector direction1(px,py,pz);
00365
00366 direction1.rotateUz(parent_polarization);
00367
00368 G4DynamicParticle * daughterparticle1
00369 = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1);
00370
00371 products->PushProducts(daughterparticle1);
00372
00373
00374
00375
00376 G4double energy2 = parentmass*(1.0 - (x+y)/2.0);
00377
00378 G4double vmass = std::sqrt((energy2-
00379 (daughtermomentum[0]+daughtermomentum[1]))*
00380 (energy2+
00381 (daughtermomentum[0]+daughtermomentum[1])));
00382 G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
00383 beta = -1.0 * std::min(beta,0.99);
00384
00385 G4double costhetan = 2.*G4UniformRand()-1.0;
00386 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00387 G4double phin = twopi*G4UniformRand()*rad;
00388 G4double sinphin = std::sin(phin);
00389 G4double cosphin = std::cos(phin);
00390
00391 G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
00392
00393 G4DynamicParticle * daughterparticle2
00394 = new G4DynamicParticle( daughters[2], direction2*(vmass/2.));
00395 G4DynamicParticle * daughterparticle3
00396 = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.));
00397
00398
00399
00400 G4ThreeVector direction34(direction0.x()+direction1.x(),
00401 direction0.y()+direction1.y(),
00402 direction0.z()+direction1.z());
00403 direction34 = direction34.unit();
00404
00405 G4LorentzVector p4 = daughterparticle2->Get4Momentum();
00406 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
00407 daughterparticle2->Set4Momentum(p4);
00408
00409 p4 = daughterparticle3->Get4Momentum();
00410 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
00411 daughterparticle3->Set4Momentum(p4);
00412
00413 products->PushProducts(daughterparticle2);
00414 products->PushProducts(daughterparticle3);
00415
00416 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
00417 daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
00418
00419
00420 #ifdef G4VERBOSE
00421 if (GetVerboseLevel()>1) {
00422 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
00423 G4cout << " create decay products in rest frame " <<G4endl;
00424 products->DumpInfo();
00425 }
00426 #endif
00427 return products;
00428 }
00429
00430 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu,
00431 G4double x,
00432 G4double y,
00433 G4double cthetaE,
00434 G4double cthetaG,
00435 G4double cthetaGE)
00436 {
00437 G4double mu = 105.65;
00438 G4double me = 0.511;
00439 G4double rho = 0.75;
00440 G4double del = 0.75;
00441 G4double eps = 0.0;
00442 G4double kap = 0.0;
00443 G4double ksi = 1.0;
00444
00445 G4double delta = 1-cthetaGE;
00446
00447
00448
00449 G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y)
00450 +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x));
00451 G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y))
00452 -2.0*(x*x)*(1.0-y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
00453 G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y))
00454 -(x*x*x)*y*(4.0+3.0*y));
00455 G4double f2s = 1.5*((x*x*x)*(y*y)*(2.0+y));
00456
00457 G4double f_1se = 12.0*(x*y*(1.0-y)+(x*x)*(2.0-3.0*y)
00458 -2.0*(x*x*x));
00459 G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y))
00460 +(x*x*x)*(2.0+3.0*y));
00461 G4double f1se = -3.0*(x*x*x)*y*(2.0+y);
00462 G4double f2se = 0.0;
00463
00464 G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y)
00465 -(x*x)*y);
00466 G4double f0sg = 6.0*(-x*(y*y)*(2.0-3.0*y)-(x*x)*y*(1.0-4.0*y)
00467 +(x*x*x)*y);
00468 G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y)
00469 -2.0*(x*x*x)*(y*y));
00470 G4double f2sg = 1.5*(x*x*x)*(y*y*y);
00471
00472 G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y)
00473 +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x));
00474 G4double f0v = 8.0*(-x*y*(3.0-y-(y*y))-(x*x)*(3.0-y-4.0*(y*y))
00475 +2.0*(x*x*x)*(1.0+2.0*y));
00476 G4double f1v = 2.0*((x*x)*y*(6.0-5.0*y-2.0*(y*y))
00477 -2.0*(x*x*x)*y*(4.0+3.0*y));
00478 G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y);
00479
00480 G4double f_1ve = 8.0*(x*y*(1.0-2.0*y)
00481 +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x));
00482 G4double f0ve = 4.0*(-(x*x)*(2.0-3.0*y-4.0*(y*y))
00483 +2.0*(x*x*x)*(2.0+3.0*y));
00484 G4double f1ve = -4.0*(x*x*x)*y*(2.0+y);
00485 G4double f2ve = 0.0;
00486
00487 G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y)
00488 -2.0*(x*x)*y);
00489 G4double f0vg = 4.0*(2.0*x*(y*y)*(1.0+y)-(x*x)*y*(1.0-4.0*y)
00490 +2.0*(x*x*x)*y);
00491 G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y)
00492 -4.0*(x*x*x)*(y*y));
00493 G4double f2vg = 2.0*(x*x*x)*(y*y*y);
00494
00495 G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y)
00496 +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x));
00497 G4double f0t = 4.0*(-x*y*(6.0+(y*y))
00498 -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
00499 G4double f1t = 2.0*((x*x)*y*(6.0-5.0*y+(y*y))
00500 -(x*x*x)*y*(4.0+3.0*y));
00501 G4double f2t = (x*x*x)*(y*y)*(2.0+y);
00502
00503 G4double f_1te = -8.0*(x*y*(1.0+3.0*y)+(x*x)*(2.0+3.0*y)
00504 +2.0*(x*x*x));
00505 G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y))
00506 +(x*x*x)*(2.0+3.0*y));
00507 G4double f1te = -2.0*(x*x*x)*y*(2.0+y);
00508 G4double f2te = 0.0;
00509
00510 G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y);
00511 G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y)
00512 +(x*x*x)*y);
00513 G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y);
00514 G4double f2tg = (x*x*x)*(y*y*y);
00515
00516 G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x));
00517 term = 1.0/term;
00518
00519 G4double nss = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s;
00520 G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v;
00521 G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t;
00522
00523 G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se;
00524 G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve;
00525 G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te;
00526
00527 G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg;
00528 G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg;
00529 G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg;
00530
00531 G4double term1 = nv;
00532 G4double term2 = 2.0*nss+nv-nt;
00533 G4double term3 = 2.0*nss-2.0*nv+nt;
00534
00535 G4double term1e = 1.0/3.0*(1.0-4.0/3.0*del);
00536 G4double term2e = 2.0*nse+5.0*nve-nte;
00537 G4double term3e = 2.0*nse-2.0*nve+nte;
00538
00539 G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del);
00540 G4double term2g = 2.0*nsg+5.0*nvg-ntg;
00541 G4double term3g = 2.0*nsg-2.0*nvg+ntg;
00542
00543 G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3;
00544 G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e)
00545 +cthetaG*(nvg-term1g*term2g+kap*term3g));
00546
00547 G4double som0 = (som00+som01)/y;
00548 som0 = fine_structure_const/8./(twopi*twopi*twopi)*som0;
00549
00550
00551
00552
00553 return som0;
00554 }