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 #include "G4MuonDecayChannelWithSpin.hh"
00043
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "Randomize.hh"
00047
00048 #include "G4DecayProducts.hh"
00049 #include "G4LorentzVector.hh"
00050
00051 G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin()
00052 : G4MuonDecayChannel(),
00053 parent_polarization(),
00054 EMMU( 0.*MeV),
00055 EMASS( 0.*MeV)
00056 {
00057 }
00058
00059 G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName,
00060 G4double theBR)
00061 : G4MuonDecayChannel(theParentName,theBR),
00062 parent_polarization(),
00063 EMMU( 0.*MeV),
00064 EMASS( 0.*MeV)
00065 {
00066 }
00067
00068 G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin()
00069 {
00070 }
00071
00072 G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4MuonDecayChannelWithSpin &right):
00073 G4MuonDecayChannel(right)
00074 {
00075 parent_polarization = right.parent_polarization;
00076 EMMU = right.EMMU;
00077 EMASS = right.EMASS;
00078 }
00079
00080 G4MuonDecayChannelWithSpin & G4MuonDecayChannelWithSpin::operator=(const G4MuonDecayChannelWithSpin & right)
00081 {
00082 if (this != &right) {
00083 kinematics_name = right.kinematics_name;
00084 verboseLevel = right.verboseLevel;
00085 rbranch = right.rbranch;
00086
00087
00088 parent_name = new G4String(*right.parent_name);
00089
00090
00091 ClearDaughtersName();
00092
00093
00094 numberOfDaughters = right.numberOfDaughters;
00095 if ( numberOfDaughters >0 ) {
00096 if (daughters_name !=0) ClearDaughtersName();
00097 daughters_name = new G4String*[numberOfDaughters];
00098
00099 for (G4int index=0; index < numberOfDaughters; index++) {
00100 daughters_name[index] = new G4String(*right.daughters_name[index]);
00101 }
00102 }
00103 parent_polarization = right.parent_polarization;
00104 EMMU = right.EMMU;
00105 EMASS = right.EMASS;
00106 }
00107 return *this;
00108 }
00109
00110
00111 G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double)
00112 {
00113
00114
00115
00116
00117 #ifdef G4VERBOSE
00118 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
00119 #endif
00120
00121 if (parent == 0) FillParent();
00122 if (daughters == 0) FillDaughters();
00123
00124
00125 G4double parentmass = parent->GetPDGMass();
00126
00127 EMMU = parentmass;
00128
00129
00130 G4double daughtermass[3];
00131 G4double sumofdaughtermass = 0.0;
00132 for (G4int index=0; index<3; index++){
00133 daughtermass[index] = daughters[index]->GetPDGMass();
00134 sumofdaughtermass += daughtermass[index];
00135 }
00136
00137 EMASS = daughtermass[0];
00138
00139
00140 G4ThreeVector dummy;
00141 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00142
00143 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00144 delete parentparticle;
00145
00146
00147
00148 G4double michel_rho = 0.75;
00149 G4double michel_delta = 0.75;
00150 G4double michel_xsi = 1.00;
00151 G4double michel_eta = 0.00;
00152
00153 G4double rndm, x, ctheta;
00154
00155 G4double FG;
00156 G4double FG_max = 2.00;
00157
00158 G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
00159 G4double x0 = EMASS/W_mue;
00160
00161 G4double x0_squared = x0*x0;
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 do{
00172
00173
00174
00175 rndm = G4UniformRand();
00176
00177 x = x0 + rndm*(1.-x0);
00178
00179 G4double x_squared = x*x;
00180
00181 G4double F_IS, F_AS, G_IS, G_AS;
00182
00183 F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
00184 F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
00185
00186 G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
00187 G_IS = G_IS + michel_eta*(1.-x)*x0;
00188
00189 G_AS = 3.*(michel_xsi-1.)*(1.-x);
00190 G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
00191 G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
00192
00193 F_IS = F_IS + G_IS;
00194 F_AS = F_AS + G_AS;
00195
00196
00197
00198 G4double R_IS = F_c(x,x0);
00199
00200 G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
00201
00202
00203
00204 G4double R_AS = F_theta(x,x0);
00205
00206 rndm = G4UniformRand();
00207
00208 ctheta = 2.*rndm-1.;
00209
00210 G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
00211
00212 FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
00213
00214 if(FG>FG_max){
00215 G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
00216 FG_max = FG;
00217 }
00218
00219 rndm = G4UniformRand();
00220
00221 }while(FG<rndm*FG_max);
00222
00223 G4double energy = x * W_mue;
00224
00225 rndm = G4UniformRand();
00226
00227 G4double phi = twopi * rndm;
00228
00229 if(energy < EMASS) energy = EMASS;
00230
00231
00232 G4double daughtermomentum[3];
00233
00234 daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
00235
00236 G4double stheta = std::sqrt(1.-ctheta*ctheta);
00237 G4double cphi = std::cos(phi);
00238 G4double sphi = std::sin(phi);
00239
00240
00241
00242 G4double px = stheta*cphi;
00243 G4double py = stheta*sphi;
00244 G4double pz = ctheta;
00245
00246 G4ThreeVector direction0(px,py,pz);
00247
00248 direction0.rotateUz(parent_polarization);
00249
00250 G4DynamicParticle * daughterparticle0
00251 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
00252
00253 products->PushProducts(daughterparticle0);
00254
00255
00256
00257
00258 G4double energy2 = parentmass*(1.0 - x/2.0);
00259 G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
00260 G4double beta = -1.0*daughtermomentum[0]/energy2;
00261 G4double costhetan = 2.*G4UniformRand()-1.0;
00262 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00263 G4double phin = twopi*G4UniformRand()*rad;
00264 G4double sinphin = std::sin(phin);
00265 G4double cosphin = std::cos(phin);
00266
00267 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
00268 G4DynamicParticle * daughterparticle1
00269 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
00270 G4DynamicParticle * daughterparticle2
00271 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
00272
00273
00274 G4LorentzVector p4;
00275 p4 = daughterparticle1->Get4Momentum();
00276 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00277 daughterparticle1->Set4Momentum(p4);
00278 p4 = daughterparticle2->Get4Momentum();
00279 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00280 daughterparticle2->Set4Momentum(p4);
00281 products->PushProducts(daughterparticle1);
00282 products->PushProducts(daughterparticle2);
00283 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
00284 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
00285
00286
00287 #ifdef G4VERBOSE
00288 if (GetVerboseLevel()>1) {
00289 G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
00290 G4cout << " create decay products in rest frame " <<G4endl;
00291 products->DumpInfo();
00292 }
00293 #endif
00294 return products;
00295 }
00296
00297 G4double G4MuonDecayChannelWithSpin::R_c(G4double x){
00298
00299 G4int n_max = (int)(100.*x);
00300
00301 if(n_max<10)n_max=10;
00302
00303 G4double L2 = 0.0;
00304
00305 for(G4int n=1; n<=n_max; n++){
00306 L2 += std::pow(x,n)/(n*n);
00307 }
00308
00309 G4double omega = std::log(EMMU/EMASS);
00310
00311 G4double r_c;
00312
00313 r_c = 2.*L2-(pi*pi/3.)-2.;
00314 r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
00315 r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
00316 r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
00317
00318 return r_c;
00319 }