52 parentmass(0.), theDaughterMasses(0)
59 G4int theNumberOfDaughters,
86 G4int theNumberOfDaughters,
96 parentmass(theParentMass),
105 G4int theNumberOfDaughters,
112 theNumberOfDaughters,
116 parentmass(theParentMass),
117 theDaughterMasses(masses)
137 G4cout <<
"G4GeneralPhaseSpaceDecay::DecayIt ";
155 G4cout <<
"G4GeneralPhaseSpaceDecay::DecayIt ";
174 delete parentparticle;
182 G4cout <<
"G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
183 G4cout <<
" create decay products in rest frame " <<
G4endl;
196 if ( theDaughterMasses )
198 daughtermass[0]= *(theDaughterMasses);
199 daughtermass[1] = *(theDaughterMasses+1);
213 delete parentparticle;
216 daughtermomentum =
Pmx(parentmass,daughtermass[0],daughtermass[1]);
218 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
223 G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum);
226 Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
232 G4cout <<
"G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
233 G4cout <<
" create decay products in rest frame " <<
G4endl;
247 for (
G4int index=0; index<3; index++)
249 if ( theDaughterMasses )
251 daughtermass[index]= *(theDaughterMasses+index);
255 sumofdaughtermass += daughtermass[index];
264 delete parentparticle;
270 G4double momentummax=0.0, momentumsum = 0.0;
287 energy = rd2*(parentmass - sumofdaughtermass);
288 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
289 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
290 momentumsum += daughtermomentum[0];
293 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
294 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
295 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
296 momentumsum += daughtermomentum[1];
299 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
300 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
301 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
302 momentumsum += daughtermomentum[2];
303 }
while (momentummax > momentumsum - momentummax );
307 G4cout <<
" daughter 0:" << daughtermomentum[0]/GeV <<
"[GeV/c]" <<
G4endl;
308 G4cout <<
" daughter 1:" << daughtermomentum[1]/GeV <<
"[GeV/c]" <<
G4endl;
309 G4cout <<
" daughter 2:" << daughtermomentum[2]/GeV <<
"[GeV/c]" <<
G4endl;
310 G4cout <<
" momentum sum:" << momentumsum/GeV <<
"[GeV/c]" <<
G4endl;
314 G4double costheta, sintheta, phi, sinphi, cosphi;
315 G4double costhetan, sinthetan, phin, sinphin, cosphin;
317 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
319 sinphi = std::sin(phi);
320 cosphi = std::cos(phi);
322 G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
327 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
328 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
330 sinphin = std::sin(phin);
331 cosphin = std::cos(phin);
333 direction2.
setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
334 direction2.
setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
335 direction2.
setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
336 Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.
mag2());
339 G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.
mag()))*(-1.0);
340 Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.
mag2() );
346 G4cout <<
"G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
347 G4cout <<
" create decay products in rest frame " <<
G4endl;
374 sumofdaughtermass += daughtermass[index];
384 G4int numberOfTry = 0;
397 if (rd[index1] < rd[index2]){
399 rd[index1] = rd[index2];
406 tmas = parentmass - sumofdaughtermass;
407 temp = sumofdaughtermass;
409 sm[index1] = rd[index1]*tmas + temp;
410 temp -= daughtermass[index1];
412 G4cout << index1 <<
" rundom number:" << rd[index1];
413 G4cout <<
" virtual mass:" << sm[index1]/GeV <<
"[GeV/c/c]" <<
G4endl;
421 daughtermomentum[index1]=
Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
424 G4cout <<
" momentum:" << daughtermomentum[index1]/GeV <<
"[GeV/c]" <<
G4endl;
428 daughtermomentum[index1]=
Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
429 if(daughtermomentum[index1] < 0.0) {
432 G4cout <<
"G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
433 G4cout <<
" can not calculate daughter momentum " <<
G4endl;
435 G4cout <<
" mass:" << parentmass/GeV <<
"[GeV/c/c]" <<
G4endl;
437 G4cout <<
" mass:" << daughtermass[index1]/GeV <<
"[GeV/c/c]" ;
438 G4cout <<
" mass:" << daughtermomentum[index1]/GeV <<
"[GeV/c]" <<
G4endl;
441 delete [] daughtermass;
442 delete [] daughtermomentum;
447 weight *= daughtermomentum[index1]/sm[index1];
450 G4cout <<
" momentum:" << daughtermomentum[index1]/GeV <<
"[GeV/c]" <<
G4endl;
459 if (numberOfTry++ >100) {
461 G4cout <<
"G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
462 G4cout <<
" can not determine Decay Kinematics " <<
G4endl;
465 delete [] daughtermass;
466 delete [] daughtermomentum;
471 G4cout <<
"Start calulation of daughters momentum vector "<<
G4endl;
480 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
482 direction.
setZ(costheta);
483 direction.
setY(sintheta*std::sin(phi));
484 direction.
setX(sintheta*std::cos(phi));
491 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
493 direction.
setZ(costheta);
494 direction.
setY(sintheta*std::sin(phi));
495 direction.
setX(sintheta*std::cos(phi));
498 beta = daughtermomentum[index1];
499 beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
506 p4.
boost( direction.
x()*beta, direction.
y()*beta, direction.
z()*beta);
517 direction.
setX(1.0); direction.
setY(0.0); direction.
setZ(0.0);
520 delete parentparticle;
525 G4cout <<
"G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
526 G4cout <<
" create decay products in rest frame " <<
G4endl;
530 delete [] daughterparticle;
531 delete [] daughtermomentum;
532 delete [] daughtermass;
G4DLLIMPORT std::ostream G4cout
HepLorentzVector & boost(double, double, double)
G4int PushProducts(G4DynamicParticle *aParticle)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4DecayProducts * OneBodyDecayIt()
static G4double Pmx(G4double e, G4double p1, G4double p2)
virtual ~G4GeneralPhaseSpaceDecay()
G4DecayProducts * ThreeBodyDecayIt()
G4DecayProducts * ManyBodyDecayIt()
G4DecayProducts * TwoBodyDecayIt()
virtual G4DecayProducts * DecayIt(G4double mass=0.0)
G4GeneralPhaseSpaceDecay(G4int Verbose=1)
G4double GetPDGMass() const
G4String ** daughters_name
G4int GetVerboseLevel() const
G4ParticleDefinition * parent
G4ParticleDefinition ** daughters