#include <G4GeneralPhaseSpaceDecay.hh>
Inheritance diagram for G4GeneralPhaseSpaceDecay:
Public Member Functions | |
G4GeneralPhaseSpaceDecay (G4int Verbose=1) | |
G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="") | |
G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theParentMass, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="") | |
G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theParentMass, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2, const G4String &theDaughterName3, const G4double *masses) | |
virtual | ~G4GeneralPhaseSpaceDecay () |
G4double | GetParentMass () const |
void | SetParentMass (const G4double aParentMass) |
virtual G4DecayProducts * | DecayIt (G4double mass=0.0) |
Static Public Member Functions | |
static G4double | Pmx (G4double e, G4double p1, G4double p2) |
Protected Member Functions | |
G4DecayProducts * | OneBodyDecayIt () |
G4DecayProducts * | TwoBodyDecayIt () |
G4DecayProducts * | ThreeBodyDecayIt () |
G4DecayProducts * | ManyBodyDecayIt () |
Definition at line 45 of file G4GeneralPhaseSpaceDecay.hh.
G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay | ( | G4int | Verbose = 1 |
) |
Definition at line 50 of file G4GeneralPhaseSpaceDecay.cc.
References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().
00050 : 00051 G4VDecayChannel("Phase Space", Verbose), 00052 parentmass(0.), theDaughterMasses(0) 00053 { 00054 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl; 00055 }
G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay | ( | const G4String & | theParentName, | |
G4double | theBR, | |||
G4int | theNumberOfDaughters, | |||
const G4String & | theDaughterName1, | |||
const G4String & | theDaughterName2 = "" , |
|||
const G4String & | theDaughterName3 = "" | |||
) |
Definition at line 57 of file G4GeneralPhaseSpaceDecay.cc.
References G4cout, G4endl, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), and G4VDecayChannel::parent.
00062 : 00063 G4VDecayChannel("Phase Space", 00064 theParentName,theBR, 00065 theNumberOfDaughters, 00066 theDaughterName1, 00067 theDaughterName2, 00068 theDaughterName3), 00069 theDaughterMasses(0) 00070 { 00071 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl; 00072 00073 // Set the parent particle (resonance) mass to the (default) PDG vale 00074 if (parent != NULL) 00075 { 00076 parentmass = parent->GetPDGMass(); 00077 } else { 00078 parentmass=0.; 00079 } 00080 00081 }
G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay | ( | const G4String & | theParentName, | |
G4double | theParentMass, | |||
G4double | theBR, | |||
G4int | theNumberOfDaughters, | |||
const G4String & | theDaughterName1, | |||
const G4String & | theDaughterName2 = "" , |
|||
const G4String & | theDaughterName3 = "" | |||
) |
Definition at line 83 of file G4GeneralPhaseSpaceDecay.cc.
References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().
00089 : 00090 G4VDecayChannel("Phase Space", 00091 theParentName,theBR, 00092 theNumberOfDaughters, 00093 theDaughterName1, 00094 theDaughterName2, 00095 theDaughterName3), 00096 parentmass(theParentMass), 00097 theDaughterMasses(0) 00098 { 00099 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl; 00100 }
G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay | ( | const G4String & | theParentName, | |
G4double | theParentMass, | |||
G4double | theBR, | |||
G4int | theNumberOfDaughters, | |||
const G4String & | theDaughterName1, | |||
const G4String & | theDaughterName2, | |||
const G4String & | theDaughterName3, | |||
const G4double * | masses | |||
) |
Definition at line 102 of file G4GeneralPhaseSpaceDecay.cc.
References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().
00109 : 00110 G4VDecayChannel("Phase Space", 00111 theParentName,theBR, 00112 theNumberOfDaughters, 00113 theDaughterName1, 00114 theDaughterName2, 00115 theDaughterName3), 00116 parentmass(theParentMass), 00117 theDaughterMasses(masses) 00118 { 00119 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl; 00120 }
G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay | ( | ) | [virtual] |
G4DecayProducts * G4GeneralPhaseSpaceDecay::DecayIt | ( | G4double | mass = 0.0 |
) | [virtual] |
Implements G4VDecayChannel.
Reimplemented in G4NuclearDecayChannel.
Definition at line 126 of file G4GeneralPhaseSpaceDecay.cc.
References G4VDecayChannel::daughters, G4VDecayChannel::DumpInfo(), G4VDecayChannel::FillDaughters(), G4VDecayChannel::FillParent(), G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), ManyBodyDecayIt(), G4VDecayChannel::numberOfDaughters, OneBodyDecayIt(), G4VDecayChannel::parent, G4VDecayChannel::parent_name, ThreeBodyDecayIt(), and TwoBodyDecayIt().
Referenced by G4KineticTrack::Decay().
00127 { 00128 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt "; 00129 G4DecayProducts * products = NULL; 00130 00131 if (parent == NULL) FillParent(); 00132 if (daughters == NULL) FillDaughters(); 00133 00134 switch (numberOfDaughters){ 00135 case 0: 00136 if (GetVerboseLevel()>0) { 00137 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt "; 00138 G4cout << " daughters not defined " <<G4endl; 00139 } 00140 break; 00141 case 1: 00142 products = OneBodyDecayIt(); 00143 break; 00144 case 2: 00145 products = TwoBodyDecayIt(); 00146 break; 00147 case 3: 00148 products = ThreeBodyDecayIt(); 00149 break; 00150 default: 00151 products = ManyBodyDecayIt(); 00152 break; 00153 } 00154 if ((products == NULL) && (GetVerboseLevel()>0)) { 00155 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt "; 00156 G4cout << *parent_name << " can not decay " << G4endl; 00157 DumpInfo(); 00158 } 00159 return products; 00160 }
G4double G4GeneralPhaseSpaceDecay::GetParentMass | ( | ) | const [inline] |
G4DecayProducts * G4GeneralPhaseSpaceDecay::ManyBodyDecayIt | ( | ) | [protected] |
Definition at line 353 of file G4GeneralPhaseSpaceDecay.cc.
References G4VDecayChannel::daughters, G4VDecayChannel::daughters_name, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::numberOfDaughters, G4VDecayChannel::parent, G4VDecayChannel::parent_name, Pmx(), G4DecayProducts::PushProducts(), G4DynamicParticle::Set4Momentum(), and G4InuclParticleNames::sm.
Referenced by DecayIt().
00361 : 19/Apr/1995 00362 // 00363 { 00364 //return value 00365 G4DecayProducts *products; 00366 00367 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl; 00368 00369 //daughters'mass 00370 G4double *daughtermass = new G4double[numberOfDaughters]; 00371 G4double sumofdaughtermass = 0.0; 00372 for (G4int index=0; index<numberOfDaughters; index++){ 00373 daughtermass[index] = daughters[index]->GetPDGMass(); 00374 sumofdaughtermass += daughtermass[index]; 00375 } 00376 00377 //Calculate daughter momentum 00378 G4double *daughtermomentum = new G4double[numberOfDaughters]; 00379 G4ParticleMomentum direction; 00380 G4DynamicParticle **daughterparticle; 00381 G4double *sm = new G4double[numberOfDaughters]; 00382 G4double tmas; 00383 G4double weight = 1.0; 00384 G4int numberOfTry = 0; 00385 G4int index1; 00386 00387 do { 00388 //Generate rundom number in descending order 00389 G4double temp; 00390 G4double *rd = new G4double[numberOfDaughters]; 00391 rd[0] = 1.0; 00392 for(index1 =1; index1 < numberOfDaughters -1; index1++) 00393 rd[index1] = G4UniformRand(); 00394 rd[ numberOfDaughters -1] = 0.0; 00395 for(index1 =1; index1 < numberOfDaughters -1; index1++) { 00396 for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) { 00397 if (rd[index1] < rd[index2]){ 00398 temp = rd[index1]; 00399 rd[index1] = rd[index2]; 00400 rd[index2] = temp; 00401 } 00402 } 00403 } 00404 00405 //calcurate virtual mass 00406 tmas = parentmass - sumofdaughtermass; 00407 temp = sumofdaughtermass; 00408 for(index1 =0; index1 < numberOfDaughters; index1++) { 00409 sm[index1] = rd[index1]*tmas + temp; 00410 temp -= daughtermass[index1]; 00411 if (GetVerboseLevel()>1) { 00412 G4cout << index1 << " rundom number:" << rd[index1]; 00413 G4cout << " virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl; 00414 } 00415 } 00416 delete [] rd; 00417 00418 //Calculate daughter momentum 00419 weight = 1.0; 00420 index1 =numberOfDaughters-1; 00421 daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]); 00422 if (GetVerboseLevel()>1) { 00423 G4cout << " daughter " << index1 << ":" << *daughters_name[index1]; 00424 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl; 00425 } 00426 for(index1 =numberOfDaughters-2; index1>=0; index1--) { 00427 // calculate 00428 daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]); 00429 if(daughtermomentum[index1] < 0.0) { 00430 // !!! illegal momentum !!! 00431 if (GetVerboseLevel()>0) { 00432 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt "; 00433 G4cout << " can not calculate daughter momentum " <<G4endl; 00434 G4cout << " parent:" << *parent_name; 00435 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl; 00436 G4cout << " daughter " << index1 << ":" << *daughters_name[index1]; 00437 G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ; 00438 G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl; 00439 } 00440 delete [] sm; 00441 delete [] daughtermass; 00442 delete [] daughtermomentum; 00443 return NULL; // Error detection 00444 00445 } else { 00446 // calculate weight of this events 00447 weight *= daughtermomentum[index1]/sm[index1]; 00448 if (GetVerboseLevel()>1) { 00449 G4cout << " daughter " << index1 << ":" << *daughters_name[index1]; 00450 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl; 00451 } 00452 } 00453 } 00454 if (GetVerboseLevel()>1) { 00455 G4cout << " weight: " << weight <<G4endl; 00456 } 00457 00458 // exit if number of Try exceeds 100 00459 if (numberOfTry++ >100) { 00460 if (GetVerboseLevel()>0) { 00461 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: "; 00462 G4cout << " can not determine Decay Kinematics " << G4endl; 00463 } 00464 delete [] sm; 00465 delete [] daughtermass; 00466 delete [] daughtermomentum; 00467 return NULL; // Error detection 00468 } 00469 } while ( weight > G4UniformRand()); 00470 if (GetVerboseLevel()>1) { 00471 G4cout << "Start calulation of daughters momentum vector "<<G4endl; 00472 } 00473 00474 G4double costheta, sintheta, phi; 00475 G4double beta; 00476 daughterparticle = new G4DynamicParticle*[numberOfDaughters]; 00477 00478 index1 = numberOfDaughters -2; 00479 costheta = 2.*G4UniformRand()-1.0; 00480 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); 00481 phi = twopi*G4UniformRand()*rad; 00482 direction.setZ(costheta); 00483 direction.setY(sintheta*std::sin(phi)); 00484 direction.setX(sintheta*std::cos(phi)); 00485 daughterparticle[index1] = new G4DynamicParticle( daughters[index1], direction*daughtermomentum[index1] ); 00486 daughterparticle[index1+1] = new G4DynamicParticle( daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) ); 00487 00488 for (index1 = numberOfDaughters -3; index1 >= 0; index1--) { 00489 //calculate momentum direction 00490 costheta = 2.*G4UniformRand()-1.0; 00491 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); 00492 phi = twopi*G4UniformRand()*rad; 00493 direction.setZ(costheta); 00494 direction.setY(sintheta*std::sin(phi)); 00495 direction.setX(sintheta*std::cos(phi)); 00496 00497 // boost already created particles 00498 beta = daughtermomentum[index1]; 00499 beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] ); 00500 for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) { 00501 G4LorentzVector p4; 00502 // make G4LorentzVector for secondaries 00503 p4 = daughterparticle[index2]->Get4Momentum(); 00504 00505 // boost secondaries to new frame 00506 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta); 00507 00508 // change energy/momentum 00509 daughterparticle[index2]->Set4Momentum(p4); 00510 } 00511 //create daughter G4DynamicParticle 00512 daughterparticle[index1]= new G4DynamicParticle( daughters[index1], direction*(-1.0*daughtermomentum[index1])); 00513 } 00514 00515 //create G4Decayproducts 00516 G4DynamicParticle *parentparticle; 00517 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0); 00518 parentparticle = new G4DynamicParticle( parent, direction, 0.0); 00519 products = new G4DecayProducts(*parentparticle); 00520 delete parentparticle; 00521 for (index1 = 0; index1<numberOfDaughters; index1++) { 00522 products->PushProducts(daughterparticle[index1]); 00523 } 00524 if (GetVerboseLevel()>1) { 00525 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt "; 00526 G4cout << " create decay products in rest frame " << G4endl; 00527 products->DumpInfo(); 00528 } 00529 00530 delete [] daughterparticle; 00531 delete [] daughtermomentum; 00532 delete [] daughtermass; 00533 delete [] sm; 00534 00535 return products; 00536 }
G4DecayProducts * G4GeneralPhaseSpaceDecay::OneBodyDecayIt | ( | ) | [protected] |
Definition at line 162 of file G4GeneralPhaseSpaceDecay.cc.
References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, and G4DecayProducts::PushProducts().
Referenced by G4NuclearDecayChannel::DecayIt(), and DecayIt().
00163 { 00164 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl; 00165 00166 // G4double daughtermass = daughters[0]->GetPDGMass(); 00167 00168 //create parent G4DynamicParticle at rest 00169 G4ParticleMomentum dummy; 00170 G4DynamicParticle * parentparticle = new G4DynamicParticle(parent, dummy, 0.0); 00171 00172 //create G4Decayproducts 00173 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 00174 delete parentparticle; 00175 00176 //create daughter G4DynamicParticle at rest 00177 G4DynamicParticle * daughterparticle = new G4DynamicParticle(daughters[0], dummy, 0.0); 00178 products->PushProducts(daughterparticle); 00179 00180 if (GetVerboseLevel()>1) 00181 { 00182 G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt "; 00183 G4cout << " create decay products in rest frame " <<G4endl; 00184 products->DumpInfo(); 00185 } 00186 return products; 00187 }
Definition at line 111 of file G4GeneralPhaseSpaceDecay.hh.
Referenced by ManyBodyDecayIt(), and TwoBodyDecayIt().
00112 { 00113 // calculate momentum of daughter particles in two-body decay 00114 if (e-p1-p2 < 0 ) 00115 { 00116 G4HadronicException(__FILE__, __LINE__, "G4GeneralPhaseSpaceDecay::Pmx energy in cms > mass1+mass2"); 00117 } 00118 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e); 00119 if (ppp>0) return std::sqrt(ppp); 00120 else return -1.; 00121 }
void G4GeneralPhaseSpaceDecay::SetParentMass | ( | const G4double | aParentMass | ) | [inline] |
Definition at line 103 of file G4GeneralPhaseSpaceDecay.hh.
Referenced by G4NuclearDecayChannel::DecayIt().
G4DecayProducts * G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt | ( | ) | [protected] |
Definition at line 239 of file G4GeneralPhaseSpaceDecay.cc.
References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, and G4DecayProducts::PushProducts().
Referenced by DecayIt().
00241 { 00242 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl; 00243 00244 //daughters'mass 00245 G4double daughtermass[3]; 00246 G4double sumofdaughtermass = 0.0; 00247 for (G4int index=0; index<3; index++) 00248 { 00249 if ( theDaughterMasses ) 00250 { 00251 daughtermass[index]= *(theDaughterMasses+index); 00252 } else { 00253 daughtermass[index] = daughters[index]->GetPDGMass(); 00254 } 00255 sumofdaughtermass += daughtermass[index]; 00256 } 00257 00258 //create parent G4DynamicParticle at rest 00259 G4ParticleMomentum dummy; 00260 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 00261 00262 //create G4Decayproducts 00263 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 00264 delete parentparticle; 00265 00266 //calculate daughter momentum 00267 // Generate two 00268 G4double rd1, rd2, rd; 00269 G4double daughtermomentum[3]; 00270 G4double momentummax=0.0, momentumsum = 0.0; 00271 G4double energy; 00272 00273 do 00274 { 00275 rd1 = G4UniformRand(); 00276 rd2 = G4UniformRand(); 00277 if (rd2 > rd1) 00278 { 00279 rd = rd1; 00280 rd1 = rd2; 00281 rd2 = rd; 00282 } 00283 momentummax = 0.0; 00284 momentumsum = 0.0; 00285 // daughter 0 00286 00287 energy = rd2*(parentmass - sumofdaughtermass); 00288 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]); 00289 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0]; 00290 momentumsum += daughtermomentum[0]; 00291 00292 // daughter 1 00293 energy = (1.-rd1)*(parentmass - sumofdaughtermass); 00294 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]); 00295 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1]; 00296 momentumsum += daughtermomentum[1]; 00297 00298 // daughter 2 00299 energy = (rd1-rd2)*(parentmass - sumofdaughtermass); 00300 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]); 00301 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2]; 00302 momentumsum += daughtermomentum[2]; 00303 } while (momentummax > momentumsum - momentummax ); 00304 00305 // output message 00306 if (GetVerboseLevel()>1) { 00307 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl; 00308 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl; 00309 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl; 00310 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl; 00311 } 00312 00313 //create daughter G4DynamicParticle 00314 G4double costheta, sintheta, phi, sinphi, cosphi; 00315 G4double costhetan, sinthetan, phin, sinphin, cosphin; 00316 costheta = 2.*G4UniformRand()-1.0; 00317 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); 00318 phi = twopi*G4UniformRand()*rad; 00319 sinphi = std::sin(phi); 00320 cosphi = std::cos(phi); 00321 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta); 00322 G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]); 00323 G4DynamicParticle * daughterparticle 00324 = new G4DynamicParticle( daughters[0], Etotal, direction0*daughtermomentum[0]); 00325 products->PushProducts(daughterparticle); 00326 00327 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]); 00328 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); 00329 phin = twopi*G4UniformRand()*rad; 00330 sinphin = std::sin(phin); 00331 cosphin = std::cos(phin); 00332 G4ParticleMomentum direction2; 00333 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 00334 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 00335 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta); 00336 Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2()); 00337 daughterparticle = new G4DynamicParticle( daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag())); 00338 products->PushProducts(daughterparticle); 00339 G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0); 00340 Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() ); 00341 daughterparticle = 00342 new G4DynamicParticle(daughters[1], Etotal, mom); 00343 products->PushProducts(daughterparticle); 00344 00345 if (GetVerboseLevel()>1) { 00346 G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt "; 00347 G4cout << " create decay products in rest frame " <<G4endl; 00348 products->DumpInfo(); 00349 } 00350 return products; 00351 }
G4DecayProducts * G4GeneralPhaseSpaceDecay::TwoBodyDecayIt | ( | ) | [protected] |
Definition at line 189 of file G4GeneralPhaseSpaceDecay.cc.
References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, Pmx(), and G4DecayProducts::PushProducts().
Referenced by G4NuclearDecayChannel::DecayIt(), and DecayIt().
00190 { 00191 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl; 00192 00193 //daughters'mass 00194 G4double daughtermass[2]; 00195 G4double daughtermomentum; 00196 if ( theDaughterMasses ) 00197 { 00198 daughtermass[0]= *(theDaughterMasses); 00199 daughtermass[1] = *(theDaughterMasses+1); 00200 } else { 00201 daughtermass[0] = daughters[0]->GetPDGMass(); 00202 daughtermass[1] = daughters[1]->GetPDGMass(); 00203 } 00204 00205 // G4double sumofdaughtermass = daughtermass[0] + daughtermass[1]; 00206 00207 //create parent G4DynamicParticle at rest 00208 G4ParticleMomentum dummy; 00209 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 00210 00211 //create G4Decayproducts @@GF why dummy parentparticle? 00212 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 00213 delete parentparticle; 00214 00215 //calculate daughter momentum 00216 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]); 00217 G4double costheta = 2.*G4UniformRand()-1.0; 00218 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); 00219 G4double phi = twopi*G4UniformRand()*rad; 00220 G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); 00221 00222 //create daughter G4DynamicParticle 00223 G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum); 00224 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0],Etotal, direction*daughtermomentum); 00225 products->PushProducts(daughterparticle); 00226 Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum); 00227 daughterparticle = new G4DynamicParticle( daughters[1],Etotal, direction*(-1.0*daughtermomentum)); 00228 products->PushProducts(daughterparticle); 00229 00230 if (GetVerboseLevel()>1) 00231 { 00232 G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt "; 00233 G4cout << " create decay products in rest frame " <<G4endl; 00234 products->DumpInfo(); 00235 } 00236 return products; 00237 }