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 #include "G4ParticleDefinition.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4DecayProducts.hh"
00047 #include "G4VDecayChannel.hh"
00048 #include "G4MuonDecayChannel.hh"
00049 #include "Randomize.hh"
00050 #include "G4LorentzVector.hh"
00051 #include "G4LorentzRotation.hh"
00052 #include "G4RotationMatrix.hh"
00053
00054 G4MuonDecayChannel::G4MuonDecayChannel()
00055 :G4VDecayChannel()
00056 {
00057 }
00058
00059 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName,
00060 G4double theBR)
00061 :G4VDecayChannel("Muon Decay",1)
00062 {
00063
00064 if (theParentName == "mu+") {
00065 SetBR(theBR);
00066 SetParent("mu+");
00067 SetNumberOfDaughters(3);
00068 SetDaughter(0, "e+");
00069 SetDaughter(1, "nu_e");
00070 SetDaughter(2, "anti_nu_mu");
00071 } else if (theParentName == "mu-") {
00072 SetBR(theBR);
00073 SetParent("mu-");
00074 SetNumberOfDaughters(3);
00075 SetDaughter(0, "e-");
00076 SetDaughter(1, "anti_nu_e");
00077 SetDaughter(2, "nu_mu");
00078 } else {
00079 #ifdef G4VERBOSE
00080 if (GetVerboseLevel()>0) {
00081 G4cout << "G4MuonDecayChannel:: constructor :";
00082 G4cout << " parent particle is not muon but ";
00083 G4cout << theParentName << G4endl;
00084 }
00085 #endif
00086 }
00087 }
00088
00089 G4MuonDecayChannel::G4MuonDecayChannel(const G4MuonDecayChannel &right):
00090 G4VDecayChannel(right)
00091 {
00092 }
00093
00094 G4MuonDecayChannel::~G4MuonDecayChannel()
00095 {
00096 }
00097
00098 G4MuonDecayChannel & G4MuonDecayChannel::operator=(const G4MuonDecayChannel & right)
00099 {
00100 if (this != &right) {
00101 kinematics_name = right.kinematics_name;
00102 verboseLevel = right.verboseLevel;
00103 rbranch = right.rbranch;
00104
00105
00106 parent_name = new G4String(*right.parent_name);
00107
00108
00109 ClearDaughtersName();
00110
00111
00112 numberOfDaughters = right.numberOfDaughters;
00113 if ( numberOfDaughters >0 ) {
00114 if (daughters_name !=0) ClearDaughtersName();
00115 daughters_name = new G4String*[numberOfDaughters];
00116
00117 for (G4int index=0; index < numberOfDaughters; index++) {
00118 daughters_name[index] = new G4String(*right.daughters_name[index]);
00119 }
00120 }
00121 }
00122 return *this;
00123 }
00124
00125 G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double)
00126 {
00127
00128
00129
00130 #ifdef G4VERBOSE
00131 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
00132 #endif
00133
00134 if (parent == 0) FillParent();
00135 if (daughters == 0) FillDaughters();
00136
00137
00138 G4double parentmass = parent->GetPDGMass();
00139
00140
00141 G4double daughtermass[3];
00142 G4double sumofdaughtermass = 0.0;
00143 for (G4int index=0; index<3; index++){
00144 daughtermass[index] = daughters[index]->GetPDGMass();
00145 sumofdaughtermass += daughtermass[index];
00146 }
00147
00148
00149 G4ThreeVector dummy;
00150 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00151
00152 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00153 delete parentparticle;
00154
00155
00156 G4double daughtermomentum[3];
00157
00158 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
00159 G4double x;
00160
00161 G4double Ee,Ene;
00162
00163 G4double gam;
00164 G4double EMax=parentmass/2-daughtermass[0];
00165
00166
00167
00168 do {
00169 Ee=G4UniformRand();
00170 do{
00171 x=xmax*G4UniformRand();
00172 gam=G4UniformRand();
00173 }while (gam >x*(1.-x));
00174 Ene=x;
00175 } while ( Ene < (1.-Ee));
00176 G4double Enm=(2.-Ee-Ene);
00177
00178
00179
00180
00181 G4double costheta,sintheta,rphi,rtheta,rpsi;
00182 costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
00183 sintheta=std::sqrt(1.-costheta*costheta);
00184
00185
00186 rphi=twopi*G4UniformRand()*rad;
00187 rtheta=(std::acos(2.*G4UniformRand()-1.));
00188 rpsi=twopi*G4UniformRand()*rad;
00189
00190 G4RotationMatrix rot;
00191 rot.set(rphi,rtheta,rpsi);
00192
00193
00194 daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
00195 G4ThreeVector direction0(0.0,0.0,1.0);
00196
00197 direction0 *= rot;
00198
00199 G4DynamicParticle * daughterparticle = new G4DynamicParticle ( daughters[0], direction0 * daughtermomentum[0]);
00200
00201 products->PushProducts(daughterparticle);
00202
00203
00204
00205 daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
00206 G4ThreeVector direction1(sintheta,0.0,costheta);
00207
00208 direction1 *= rot;
00209
00210 G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( daughters[1], direction1 * daughtermomentum[1]);
00211 products->PushProducts(daughterparticle1);
00212
00213
00214
00215 daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
00216 G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
00217
00218 direction2 *= rot;
00219
00220 G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( daughters[2],
00221 direction2 * daughtermomentum[2]);
00222 products->PushProducts(daughterparticle2);
00223
00224
00225
00226
00227
00228 #ifdef G4VERBOSE
00229 if (GetVerboseLevel()>1) {
00230 G4cout << "G4MuonDecayChannel::DecayIt ";
00231 G4cout << " create decay products in rest frame " <<G4endl;
00232 products->DumpInfo();
00233 }
00234 #endif
00235 return products;
00236 }
00237
00238
00239
00240
00241
00242