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 #include "G4ParticleDefinition.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4DecayProducts.hh"
00041 #include "G4VDecayChannel.hh"
00042 #include "G4KL3DecayChannel.hh"
00043 #include "Randomize.hh"
00044 #include "G4LorentzVector.hh"
00045 #include "G4LorentzRotation.hh"
00046
00047 G4KL3DecayChannel::G4KL3DecayChannel()
00048 :G4VDecayChannel(),
00049 massK(0.0), pLambda(0.0), pXi0(0.0)
00050 {
00051 daughterM[idPi] = 0.0;
00052 daughterM[idLepton] = 0.0;
00053 daughterM[idNutrino] = 0.0;
00054 }
00055
00056
00057 G4KL3DecayChannel::G4KL3DecayChannel(
00058 const G4String& theParentName,
00059 G4double theBR,
00060 const G4String& thePionName,
00061 const G4String& theLeptonName,
00062 const G4String& theNutrinoName)
00063 :G4VDecayChannel("KL3 Decay",theParentName,
00064 theBR, 3,
00065 thePionName,theLeptonName,theNutrinoName)
00066 {
00067 static const G4String K_plus("kaon+");
00068 static const G4String K_minus("kaon-");
00069 static const G4String K_L("kaon0L");
00070 static const G4String Mu_plus("mu+");
00071 static const G4String Mu_minus("mu-");
00072 static const G4String E_plus("e+");
00073 static const G4String E_minus("e-");
00074
00075 massK = 0.0;
00076 daughterM[idPi] = 0.0;
00077 daughterM[idLepton] = 0.0;
00078 daughterM[idNutrino] = 0.0;
00079
00080
00081 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
00082 ((theParentName == K_minus)&&(theLeptonName == E_minus)) ) {
00083
00084 pLambda = 0.0286;
00085 pXi0 = -0.35;
00086 } else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
00087 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) ) {
00088
00089 pLambda = 0.033;
00090 pXi0 = -0.35;
00091 } else if ( (theParentName == K_L) &&
00092 ((theLeptonName == E_plus) ||(theLeptonName == E_minus)) ){
00093
00094 pLambda = 0.0300;
00095 pXi0 = -0.11;
00096 } else if ( (theParentName == K_L) &&
00097 ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus)) ){
00098
00099 pLambda = 0.034;
00100 pXi0 = -0.11;
00101 } else {
00102 #ifdef G4VERBOSE
00103 if (GetVerboseLevel()>2) {
00104 G4cout << "G4KL3DecayChannel:: constructor :";
00105 G4cout << "illegal arguments " << G4endl;;
00106 DumpInfo();
00107 }
00108 #endif
00109
00110 pLambda = 0.0300;
00111 pXi0 = -0.11;
00112 }
00113 }
00114
00115 G4KL3DecayChannel::~G4KL3DecayChannel()
00116 {
00117 }
00118
00119 G4KL3DecayChannel::G4KL3DecayChannel(const G4KL3DecayChannel &right):
00120 G4VDecayChannel(right),
00121 massK(right.massK),
00122 pLambda(right.pLambda),
00123 pXi0(right.pXi0)
00124 {
00125 daughterM[idPi] = right.daughterM[idPi];
00126 daughterM[idLepton] = right.daughterM[idLepton];
00127 daughterM[idNutrino] = right.daughterM[idNutrino];
00128 }
00129
00130 G4KL3DecayChannel & G4KL3DecayChannel::operator=(const G4KL3DecayChannel & right)
00131 {
00132 if (this != &right) {
00133 kinematics_name = right.kinematics_name;
00134 verboseLevel = right.verboseLevel;
00135 rbranch = right.rbranch;
00136
00137
00138 parent_name = new G4String(*right.parent_name);
00139
00140
00141 ClearDaughtersName();
00142
00143
00144 numberOfDaughters = right.numberOfDaughters;
00145 if ( numberOfDaughters >0 ) {
00146 if (daughters_name !=0) ClearDaughtersName();
00147 daughters_name = new G4String*[numberOfDaughters];
00148
00149 for (G4int index=0; index < numberOfDaughters; index++) {
00150 daughters_name[index] = new G4String(*right.daughters_name[index]);
00151 }
00152 }
00153 massK = right.massK;
00154 pLambda = right.pLambda;
00155 pXi0 = right.pXi0;
00156 daughterM[idPi] = right.daughterM[idPi];
00157 daughterM[idLepton] = right.daughterM[idLepton];
00158 daughterM[idNutrino] = right.daughterM[idNutrino];
00159 }
00160 return *this;
00161 }
00162
00163
00164 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4double)
00165 {
00166
00167
00168
00169 #ifdef G4VERBOSE
00170 if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
00171 #endif
00172
00173
00174 if (parent == 0) {
00175 FillParent();
00176 }
00177 massK = parent->GetPDGMass();
00178
00179
00180 if (daughters == 0) {
00181 FillDaughters();
00182 }
00183 daughterM[idPi] = daughters[idPi]->GetPDGMass();
00184 daughterM[idLepton] = daughters[idLepton]->GetPDGMass();
00185 daughterM[idNutrino] = daughters[idNutrino]->GetPDGMass();
00186
00187
00188
00189 G4double daughterP[3], daughterE[3];
00190 G4double w;
00191 G4double r;
00192 do {
00193 r = G4UniformRand();
00194 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
00195 w = DalitzDensity(daughterE[idPi],daughterE[idLepton],daughterE[idNutrino]);
00196 } while ( r > w);
00197
00198
00199 #ifdef G4VERBOSE
00200 if (GetVerboseLevel()>1) {
00201 G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]" <<G4endl;
00202 G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]" <<G4endl;
00203 G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]" <<G4endl;
00204 }
00205 #endif
00206
00207 G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
00208 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, *direction, 0.0);
00209 delete direction;
00210
00211
00212 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00213 delete parentparticle;
00214
00215
00216 G4double costheta, sintheta, phi, sinphi, cosphi;
00217 G4double costhetan, sinthetan, phin, sinphin, cosphin;
00218
00219
00220 costheta = 2.*G4UniformRand()-1.0;
00221 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00222 phi = twopi*G4UniformRand()*rad;
00223 sinphi = std::sin(phi);
00224 cosphi = std::cos(phi);
00225 direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
00226 G4ThreeVector momentum0 = (*direction)*daughterP[0];
00227 G4DynamicParticle * daughterparticle
00228 = new G4DynamicParticle( daughters[0], momentum0);
00229 products->PushProducts(daughterparticle);
00230
00231
00232 costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
00233 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00234 phin = twopi*G4UniformRand()*rad;
00235 sinphin = std::sin(phin);
00236 cosphin = std::cos(phin);
00237 direction->setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
00238 direction->setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
00239 direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
00240
00241 G4ThreeVector momentum2 = (*direction)*daughterP[2];
00242 daughterparticle = new G4DynamicParticle( daughters[2], momentum2);
00243 products->PushProducts(daughterparticle);
00244
00245
00246 G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
00247 daughterparticle =
00248 new G4DynamicParticle( daughters[1], momentum1);
00249 products->PushProducts(daughterparticle);
00250
00251 #ifdef G4VERBOSE
00252 if (GetVerboseLevel()>1) {
00253 G4cout << "G4KL3DecayChannel::DecayIt ";
00254 G4cout << " create decay products in rest frame " <<G4endl;
00255 G4cout << " decay products address=" << products << G4endl;
00256 products->DumpInfo();
00257 }
00258 #endif
00259 delete direction;
00260 return products;
00261 }
00262
00263 void G4KL3DecayChannel::PhaseSpace(G4double parentM,
00264 const G4double* M,
00265 G4double* E,
00266 G4double* P )
00267
00268 {
00269
00270
00271 G4double sumofdaughtermass = 0.0;
00272 G4int index;
00273 for (index=0; index<3; index++){
00274 sumofdaughtermass += M[index];
00275 }
00276
00277
00278
00279 G4double rd1, rd2, rd;
00280 G4double momentummax=0.0, momentumsum = 0.0;
00281 G4double energy;
00282
00283 do {
00284 rd1 = G4UniformRand();
00285 rd2 = G4UniformRand();
00286 if (rd2 > rd1) {
00287 rd = rd1;
00288 rd1 = rd2;
00289 rd2 = rd;
00290 }
00291 momentummax = 0.0;
00292 momentumsum = 0.0;
00293
00294 energy = rd2*(parentM - sumofdaughtermass);
00295 P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
00296 E[0] = energy;
00297 if ( P[0] >momentummax )momentummax = P[0];
00298 momentumsum += P[0];
00299
00300 energy = (1.-rd1)*(parentM - sumofdaughtermass);
00301 P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
00302 E[1] = energy;
00303 if ( P[1] >momentummax )momentummax = P[1];
00304 momentumsum += P[1];
00305
00306 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
00307 P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
00308 E[2] = energy;
00309 if ( P[2] >momentummax )momentummax = P[2];
00310 momentumsum += P[2];
00311 } while (momentummax > momentumsum - momentummax );
00312
00313 #ifdef G4VERBOSE
00314 if (GetVerboseLevel()>2) {
00315 G4cout << "G4KL3DecayChannel::PhaseSpace ";
00316 G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
00317 for (index=0; index<3; index++){
00318 G4cout << index << " : " << M[index]/GeV << "GeV/c/c ";
00319 G4cout << " : " << E[index]/GeV << "GeV ";
00320 G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
00321 }
00322 }
00323 #endif
00324 }
00325
00326
00327 G4double G4KL3DecayChannel::DalitzDensity(G4double Epi, G4double El, G4double Enu)
00328 {
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 G4double massPi = daughterM[idPi];
00346 G4double massL = daughterM[idLepton];
00347 G4double massNu = daughterM[idNutrino];
00348
00349
00350 Epi = Epi + massPi;
00351 El = El + massL;
00352 Enu = Enu + massNu;
00353
00354 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
00355 G4double E = Epi_max - Epi;
00356 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
00357
00358 G4double F = 1.0 + pLambda*q2/massPi/massPi;
00359 G4double Fmax = 1.0;
00360 if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
00361
00362 G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
00363
00364 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
00365 G4double coeffB = massL*massL*(Enu-E/2.0);
00366 G4double coeffC = massL*massL*E/4.0;
00367
00368 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
00369
00370 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
00371
00372 #ifdef G4VERBOSE
00373 if (GetVerboseLevel()>2) {
00374 G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl;
00375 G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
00376 G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
00377 G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
00378 G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
00379 G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC <<G4endl;
00380 G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
00381 }
00382 #endif
00383 return (Rho/RhoMax);
00384 }
00385
00386