59 G4int theNumberOfDaughters,
71 current_parent_mass(0.0)
91 if (parentMass >0.0) current_parent_mass = parentMass;
98 G4cout <<
"G4PhaseSpaceDecayChannel::DecayIt ";
104 products = OneBodyDecayIt();
107 products = TwoBodyDecayIt();
110 products = ThreeBodyDecayIt();
113 products = ManyBodyDecayIt();
118 G4cout <<
"G4PhaseSpaceDecayChannel::DecayIt ";
137 delete parentparticle;
145 G4cout <<
"G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
146 G4cout <<
" create decay products in rest frame " <<
G4endl;
159 G4double parentmass = current_parent_mass;
172 delete parentparticle;
175 daughtermomentum =
Pmx(parentmass,daughtermass[0],daughtermass[1]);
176 if (daughtermomentum <0.0) {
179 G4cerr <<
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
180 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
186 G4Exception(
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
188 "Can not create decay products: sum of daughter mass is larger than parent mass");
193 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
195 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
205 G4cout <<
"G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
206 G4cout <<
" create decay products in rest frame " <<
G4endl;
220 G4double parentmass = current_parent_mass;
224 for (
G4int index=0; index<3; index++){
226 sumofdaughtermass += daughtermass[index];
236 delete parentparticle;
238 if (sumofdaughtermass >parentmass) {
241 G4cerr <<
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
242 <<
"sum of daughter mass is larger than parent mass" <<
G4endl;
249 G4Exception(
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
251 "Can not create decay products: sum of daughter mass is larger than parent mass");
261 G4double momentummax=0.0, momentumsum = 0.0;
275 energy = rd2*(parentmass - sumofdaughtermass);
276 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
277 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
278 momentumsum += daughtermomentum[0];
280 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
281 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
282 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
283 momentumsum += daughtermomentum[1];
285 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
286 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
287 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
288 momentumsum += daughtermomentum[2];
289 }
while (momentummax > momentumsum - momentummax );
293 G4cout <<
" daughter 0:" << daughtermomentum[0]/GeV <<
"[GeV/c]" <<
G4endl;
294 G4cout <<
" daughter 1:" << daughtermomentum[1]/GeV <<
"[GeV/c]" <<
G4endl;
295 G4cout <<
" daughter 2:" << daughtermomentum[2]/GeV <<
"[GeV/c]" <<
G4endl;
296 G4cout <<
" momentum sum:" << momentumsum/GeV <<
"[GeV/c]" <<
G4endl;
301 G4double costheta, sintheta, phi, sinphi, cosphi;
302 G4double costhetan, sinthetan, phin, sinphin, cosphin;
304 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
306 sinphi = std::sin(phi);
307 cosphi = std::cos(phi);
308 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
313 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
314 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
316 sinphin = std::sin(phin);
317 cosphin = std::cos(phin);
319 direction2.
setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
320 direction2.
setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
321 direction2.
setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
328 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.
mag()))*(-1.0)
334 G4cout <<
"G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
335 G4cout <<
" create decay products in rest frame " <<
G4endl;
362 G4double parentmass = current_parent_mass;
368 sumofdaughtermass += daughtermass[index];
378 G4int numberOfTry = 0;
389 if (rd[index] < rd[index2]){
391 rd[index] = rd[index2];
398 tmas = parentmass - sumofdaughtermass;
399 temp = sumofdaughtermass;
401 sm[index] = rd[index]*tmas + temp;
402 temp -= daughtermass[index];
404 G4cout << index <<
" rundom number:" << rd[index];
405 G4cout <<
" virtual mass:" <<
sm[index]/GeV <<
"[GeV/c/c]" <<
G4endl;
413 daughtermomentum[index]=
Pmx( sm[index-1],daughtermass[index-1], sm[index]);
417 G4cout <<
" momentum:" << daughtermomentum[index]/GeV <<
"[GeV/c]" <<
G4endl;
422 daughtermomentum[index]=
Pmx( sm[index],daughtermass[index], sm[index +1]);
423 if(daughtermomentum[index] < 0.0) {
427 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
428 G4cout <<
" can not calculate daughter momentum " <<
G4endl;
430 G4cout <<
" mass:" << parentmass/GeV <<
"[GeV/c/c]" <<
G4endl;
432 G4cout <<
" mass:" << daughtermass[index]/GeV <<
"[GeV/c/c]" ;
433 G4cout <<
" mass:" << daughtermomentum[index]/GeV <<
"[GeV/c]" <<
G4endl;
437 delete [] daughtermass;
438 delete [] daughtermomentum;
440 G4Exception(
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
442 "Can not create decay products: sum of daughter mass is larger than parent mass");
448 weight *= daughtermomentum[index]/
sm[index];
452 G4cout <<
" momentum:" << daughtermomentum[index]/GeV <<
"[GeV/c]" <<
G4endl;
465 if (numberOfTry++ >100) {
468 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
469 G4cout <<
" can not determine Decay Kinematics " <<
G4endl;
478 G4Exception(
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
480 " Cannot decay : Decay Kinematics cannot be calculated ");
483 delete [] daughtermass;
484 delete [] daughtermomentum;
491 G4cout <<
"Start calulation of daughters momentum vector "<<
G4endl;
501 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
503 direction.
setZ(costheta);
504 direction.
setY(sintheta*std::sin(phi));
505 direction.
setX(sintheta*std::cos(phi));
512 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
514 direction.
setZ(costheta);
515 direction.
setY(sintheta*std::sin(phi));
516 direction.
setX(sintheta*std::cos(phi));
519 beta = daughtermomentum[index];
520 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
527 p4.
boost( direction.
x()*beta, direction.
y()*beta, direction.
z()*beta);
538 direction.
setX(1.0); direction.
setY(0.0); direction.
setZ(0.0);
541 delete parentparticle;
548 G4cout <<
"G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
549 G4cout <<
" create decay products in rest frame " <<
G4endl;
554 delete [] daughterparticle;
555 delete [] daughtermomentum;
556 delete [] daughtermass;
566 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
567 if (ppp>0)
return std::sqrt(ppp);
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
HepLorentzVector & boost(double, double, double)
G4int PushProducts(G4DynamicParticle *aParticle)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
const G4String & GetParticleName() const
virtual ~G4PhaseSpaceDecayChannel()
static G4double Pmx(G4double e, G4double p1, G4double p2)
virtual G4DecayProducts * DecayIt(G4double)
G4PhaseSpaceDecayChannel(G4int Verbose=1)
G4String ** daughters_name
G4int GetVerboseLevel() const
G4double * daughters_mass
G4ParticleDefinition * parent
G4ParticleDefinition ** daughters
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)