60 G4int theNumberOfDaughters,
72 useGivenDaughterMass(false)
99 G4cout <<
"G4PhaseSpaceDecayChannel::DecayIt ";
119 G4cout <<
"G4PhaseSpaceDecayChannel::DecayIt ";
140 delete parentparticle;
149 G4cout <<
"G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
150 G4cout <<
" create decay products in rest frame " <<
G4endl;
166 G4double daughtermass[2], daughterwidth[2];
177 delete parentparticle;
180 G4bool withWidth = (daughterwidth[0]>1.0e-3*daughtermass[0])
181 || (daughterwidth[1]>1.0e-3*daughtermass[1]);
183 G4double sumofdaughterwidthsq = daughterwidth[0]*daughterwidth[0]+daughterwidth[1]*daughterwidth[1];
185 G4double maxDev = (parentmass - daughtermass[0] - daughtermass[1] )/std::sqrt(sumofdaughterwidthsq);
189 G4cout <<
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
190 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
196 G4Exception(
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
198 "Can not create decay products: sum of daughter mass is larger than parent mass");
203 while (dm1+dm2>parentmass){
204 dm1=
DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
205 dm2=
DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
207 daughtermass[0] = dm1;
208 daughtermass[1] = dm2;
215 if (parentmass < daughtermass[0] + daughtermass[1] ){
218 G4cout <<
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
219 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
228 G4Exception(
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
230 "Can not create decay products: sum of daughter mass is larger than parent mass");
235 G4double daughtermomentum =
Pmx(parentmass,daughtermass[0],daughtermass[1]);
238 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
240 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
243 G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[0]*daughtermass[0]) - daughtermass[0];
246 Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[1]*daughtermass[1]) - daughtermass[1];
252 G4cout <<
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
253 G4cout <<
" create decay products in rest frame " <<
G4endl;
269 G4double daughtermass[3], daughterwidth[3];
271 G4double sumofdaughterwidthsq = 0.0;
273 for (
G4int index=0; index<3; index++){
275 sumofdaughtermass += daughtermass[index];
277 sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
278 withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
288 delete parentparticle;
292 G4double maxDev = (parentmass - sumofdaughtermass )/std::sqrt(sumofdaughterwidthsq) ;
296 G4cout <<
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
297 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
304 G4Exception(
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
306 "Can not create decay products: sum of daughter mass is larger than parent mass");
312 while (dm1+dm2+dm3>parentmass){
313 dm1=
DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
314 dm2=
DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
315 dm3=
DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
317 daughtermass[0] = dm1;
318 daughtermass[1] = dm2;
319 daughtermass[2] = dm3;
320 sumofdaughtermass = dm1+dm2+dm3;
327 sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
330 if (sumofdaughtermass >parentmass) {
333 G4cout <<
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
334 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
344 G4Exception(
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
346 "Can not create decay products: sum of daughter mass is larger than parent mass");
356 G4double momentummax=0.0, momentumsum = 0.0;
358 const size_t MAX_LOOP=10000;
360 for (
size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
371 energy = rd2*(parentmass - sumofdaughtermass);
372 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
373 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
374 momentumsum += daughtermomentum[0];
376 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
377 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
378 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
379 momentumsum += daughtermomentum[1];
381 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
382 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
383 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
384 momentumsum += daughtermomentum[2];
385 if (momentummax <= momentumsum - momentummax )
break;
391 G4cout <<
" daughter 0:" << daughtermomentum[0]/
GeV <<
"[GeV/c]" <<
G4endl;
392 G4cout <<
" daughter 1:" << daughtermomentum[1]/
GeV <<
"[GeV/c]" <<
G4endl;
393 G4cout <<
" daughter 2:" << daughtermomentum[2]/
GeV <<
"[GeV/c]" <<
G4endl;
399 G4double costheta, sintheta, phi, sinphi, cosphi;
400 G4double costhetan, sinthetan, phin, sinphin, cosphin;
402 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
404 sinphi = std::sin(phi);
405 cosphi = std::cos(phi);
407 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
408 G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + daughtermass[0]*daughtermass[0]) - daughtermass[0];
411 Ekin, daughtermass[0]);
414 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
415 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
417 sinphin = std::sin(phin);
418 cosphin = std::cos(phin);
420 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
421 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
422 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
423 G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
424 Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2]) - daughtermass[2];
427 Ekin, daughtermass[2]);
430 pmom = (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
431 Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1]) - daughtermass[1];
435 Ekin, daughtermass[1]);
440 G4cout <<
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
441 G4cout <<
" create decay products in rest frame " <<
G4endl;
475 delete parentparticle;
487 sumofdaughtermass += daughtermass[index];
489 if (sumofdaughtermass >parentmass) {
492 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt "
493 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
501 G4Exception(
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
503 "Can not create decay products: sum of daughter mass is larger than parent mass");
514 G4int numberOfTry = 0;
515 const size_t MAX_LOOP=10000;
517 for (
size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
522 for(index =1; index < numberOfDaughters -1; index++) rd[index] =
G4UniformRand();
523 rd[ numberOfDaughters -1] = 0.0;
524 for(index =1; index < numberOfDaughters -1; index++) {
526 if (rd[index] < rd[index2]){
528 rd[index] = rd[index2];
535 tmas = parentmass - sumofdaughtermass;
536 temp = sumofdaughtermass;
538 sm[index] = rd[index]*tmas + temp;
539 temp -= daughtermass[index];
541 G4cout << index <<
" rundom number:" << rd[index];
550 for(index =0; index< numberOfDaughters-1 && smOK; index--) {
551 smOK = (sm[index]-daughtermass[index]- sm[index+1] >=0.);
555 index =numberOfDaughters-1;
556 daughtermomentum[index]=
Pmx( sm[index-1],daughtermass[index-1], sm[index]);
560 G4cout <<
" momentum:" << daughtermomentum[index]/
GeV <<
"[GeV/c]" <<
G4endl;
563 for(index =numberOfDaughters-2; index>=0; index--) {
565 daughtermomentum[index]=
Pmx( sm[index],daughtermass[index], sm[index +1]);
566 if(daughtermomentum[index] < 0.0) {
570 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
571 G4cout <<
" can not calculate daughter momentum " <<
G4endl;
575 G4cout <<
" mass:" << daughtermass[index]/
GeV <<
"[GeV/c/c]" ;
576 G4cout <<
" mass:" << daughtermomentum[index]/
GeV <<
"[GeV/c]" <<
G4endl;
583 delete [] daughtermass;
584 delete [] daughtermomentum;
586 G4Exception(
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
588 "Can not create decay products: sum of daughter mass is larger than parent mass");
594 weight *= daughtermomentum[index]/sm[index];
598 G4cout <<
" momentum:" << daughtermomentum[index]/
GeV <<
"[GeV/c]" <<
G4endl;
611 if (numberOfTry++ >100) {
614 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
615 G4cout <<
" can not determine Decay Kinematics " <<
G4endl;
624 G4Exception(
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
626 " Cannot decay : Decay Kinematics cannot be calculated ");
629 delete [] daughtermass;
630 delete [] daughtermomentum;
638 G4cout <<
"Start calulation of daughters momentum vector "<<
G4endl;
646 index = numberOfDaughters -2;
648 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
650 direction.setZ(costheta);
651 direction.setY(sintheta*std::sin(phi));
652 direction.setX(sintheta*std::cos(phi));
656 for (index = numberOfDaughters -3; index >= 0; index--) {
659 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
661 direction.setZ(costheta);
662 direction.setY(sintheta*std::sin(phi));
663 direction.setX(sintheta*std::cos(phi));
666 beta = daughtermomentum[index];
667 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
674 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
690 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
691 G4cout <<
" create decay products in rest frame " <<
G4endl;
696 delete [] daughterparticle;
697 delete [] daughtermomentum;
698 delete [] daughtermass;
729 return (parentMass >= sumOfDaughterMassMin);
735 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
736 if (ppp>0)
return std::sqrt(ppp);
virtual G4DecayProducts * DecayIt(G4double)
G4double * G4MT_daughters_mass
G4bool SetDaughterMasses(G4double masses[])
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
G4double givenDaughterMasses[MAX_N_DAUGHTERS]
G4ParticleDefinition * G4MT_parent
G4bool IsOKWithParentMass(G4double parentMass)
G4ParticleDefinition ** G4MT_daughters
G4double G4MT_parent_mass
virtual G4bool IsOKWithParentMass(G4double parentMass)
G4bool useGivenDaughterMass
G4bool SampleDaughterMasses()
G4DecayProducts * ManyBodyDecayIt()
G4DecayProducts * OneBodyDecayIt()
const G4String & GetParticleName() const
G4DecayProducts * ThreeBodyDecayIt()
virtual ~G4PhaseSpaceDecayChannel()
G4GLOB_DLL std::ostream G4cout
static const double twopi
G4DecayProducts * TwoBodyDecayIt()
G4LorentzVector Get4Momentum() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4ThreadLocal G4double current_parent_mass
void Set4Momentum(const G4LorentzVector &momentum)
G4int GetVerboseLevel() const
G4double * G4MT_daughters_width
G4double energy(const ThreeVector &p, const G4double m)
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
G4PhaseSpaceDecayChannel(G4int Verbose=1)
void SetMass(G4double mass)
G4String ** daughters_name
static G4double Pmx(G4double e, G4double p1, G4double p2)
CLHEP::HepLorentzVector G4LorentzVector