80 theStateToNucleus(undefined),
81 theProjectilePotential(0),
142 theDefinition(aDefinition),
143 theFormationTime(aFormationTime),
144 thePosition(aPosition),
145 the4Momentum(a4Momentum),
146 theFermi3Momentum(0),
147 theTotal4Momentum(a4Momentum),
149 theStateToNucleus(undefined),
150 theProjectilePotential(0),
171 if (theDecayTable != 0)
195 G4bool * theDaughterIsShortLived = 0;
201 for (index =
nChannels - 1; index >= 0; --index)
206 if (nDaughters == 2 || nDaughters == 3)
214 theDaughterIsShortLived =
new G4bool[nDaughters];
215 for (
G4int n = 0;
n < nDaughters; ++
n)
232 if ( !theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
252 else if ( !theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
274 else if ( theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
304 else if ( theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
331 G4int nShortLived = 0;
332 if ( theDaughterIsShortLived[0] )
336 if ( theDaughterIsShortLived[1] )
342 if ( theDaughterIsShortLived[2] )
348 if ( nShortLived == 0 )
357 else if ( nShortLived >= 1 )
375 G4double theMomRatio = theActualMom / thePoleMom;
387 delete [] theDaughterIsShortLived;
388 theDaughterIsShortLived = 0;
419 : theDefinition(
nucleon->GetDefinition()),
421 thePosition(aPosition),
422 the4Momentum(a4Momentum),
423 theFermi3Momentum(
nucleon->GetMomentum()),
426 theActualMass(
nucleon->GetDefinition()->GetPDGMass()),
430 theStateToNucleus(undefined),
431 theProjectilePotential(0),
473 return (
this == & right);
480 return (
this != & right);
499 G4cerr <<
"Error condition encountered in G4KineticTrack::Decay()"<<
G4endl;
500 G4cerr <<
" track has no particle definition associated."<<
G4endl;
506 G4cerr <<
"Error condition encountered in G4KineticTrack::Decay()"<<
G4endl;
507 G4cerr <<
" particle definition has no decay table associated."<<
G4endl;
516 if (theTotalActualWidth !=0)
521 G4bool isChannelBelowThreshold =
true;
522 const G4int maxNumberOfLoops = 10000;
523 G4int loopCounter = 0;
529 G4int theNumberOfDaughters;
543 theCumActualWidth[index] = theSumActualWidth;
553 if (r < theCumActualWidth[index])
562 delete [] theCumActualWidth;
566 G4cerr <<
"Error condition encountered in G4KineticTrack::Decay()"<<
G4endl;
567 G4cerr <<
" decay channel has 0x0 channel associated."<<
G4endl;
577 theDaughtersName1 =
"";
578 theDaughtersName2 =
"";
579 theDaughtersName3 =
"";
580 theDaughtersName4 =
"";
582 for (
G4int i=0; i < 4; ++i) masses[i]=0.;
583 G4int shortlivedDaughters[4];
584 G4int numberOfShortliveds(0);
586 for (
G4int aD=0; aD < theNumberOfDaughters ; ++aD)
592 shortlivedDaughters[numberOfShortliveds]=aD;
593 ++numberOfShortliveds;
599 switch (theNumberOfDaughters)
605 theDaughtersName2 =
"";
606 theDaughtersName3 =
"";
607 theDaughtersName4 =
"";
612 theDaughtersName3 =
"";
613 theDaughtersName4 =
"";
614 if ( numberOfShortliveds == 1)
616 G4double massmax=theParentMass - SumLongLivedMass;
618 masses[shortlivedDaughters[0]]= aSampler.
SampleMass(aDaughter,massmax);
619 }
else if ( numberOfShortliveds == 2) {
626 masses[shortlivedDaughters[
zero]]=aSampler.
SampleMass(aDaughter,massmax);
627 massmax=theParentMass - masses[shortlivedDaughters[
zero]];
628 aDaughter=theDecayChannel->
GetDaughter(shortlivedDaughters[one]);
629 masses[shortlivedDaughters[one]]=aSampler.
SampleMass(aDaughter,massmax);
636 theDaughtersName4 =
"";
637 if ( numberOfShortliveds == 1)
639 G4double massmax=theParentMass - SumLongLivedMass;
641 masses[shortlivedDaughters[0]]= aSampler.
SampleMass(aDaughter,massmax);
649 if ( numberOfShortliveds == 1)
651 G4double massmax=theParentMass - SumLongLivedMass;
653 masses[shortlivedDaughters[0]]= aSampler.
SampleMass(aDaughter,massmax);
655 if ( theNumberOfDaughters > 4 ) {
657 ed <<
"More than 4 decay daughters: kept only the first 4" <<
G4endl;
667 for (
G4int i=0; i < 4; ++i) sumDaughterMasses += masses[i];
668 if ( theParentMass - sumDaughterMasses > 0.0 ) isChannelBelowThreshold =
false;
670 }
while ( isChannelBelowThreshold && ++loopCounter < maxNumberOfLoops );
679 theNumberOfDaughters,
686 if(!theDecayProducts)
689 ed <<
"Error condition encountered: phase-space decay failed." <<
G4endl
691 <<
"\t the channel index is: "<< chosench <<
" of "<<
nChannels <<
"channels" <<
G4endl
692 <<
"\t " << theNumberOfDaughters <<
" daughter particles: "
693 << theDaughtersName1 <<
" " << theDaughtersName2 <<
" " << theDaughtersName3 <<
" "
694 << theDaughtersName4 <<
G4endl;
713 for (
G4int i=dEntries; i > 0; --i)
715 theDynamicParticle = theDecayProducts->
PopProducts();
719 momentumBalanceCMS += theDynamicParticle->
Get4Momentum();
721 energyMomentumBalance -= momentum;
727 theDecayProductList->push_back(aDaughter);
728 delete theDynamicParticle;
730 delete theDecayProducts;
731 if(std::getenv(
"DecayEnergyBalanceCheck"))
732 std::cout <<
"DEBUGGING energy balance in cms and lab, charge baryon balance : "
733 << momentumBalanceCMS <<
" "
734 <<energyMomentumBalance <<
" "
738 return theDecayProductList;
753 G4double result = (1. / (2 * mass)) *
754 std::sqrt(
std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
755 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
756 BrWig(gamma2, mass2, xmass);
766 G4double result = (1. / (2 * mass)) *
767 std::sqrt(
std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
768 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
769 BrWig(gamma2, mass2, xmass);
780 const G4double result = (1. / (2 * mass)) *
783 BrWig(gamma2, mass2, xmass);
797 const G4double theUpperLimit = mass - xmass;
798 const G4int nIterations = 100;
810 const G4int nIterations = 100;
812 if (theLowerLimit>=theUpperLimit)
return 0.0;
816 theLowerLimit, theUpperLimit, nIterations);
817 return theIntegralOverMass2;
823 const G4int nIterations = 100;
825 if (theLowerLimit>=theUpperLimit)
return 0.0;
829 theLowerLimit, theUpperLimit, nIterations);
830 return theIntegralOverMass2;
838 const G4int nIterations = 100;
840 if (theLowerLimit>=theUpperLimit)
return 0.0;
844 theLowerLimit, theUpperLimit, nIterations);
845 return theIntegralOverMass2;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
static G4ThreadLocal G4double G4KineticTrack_xmass1
static G4ThreadLocal G4double G4KineticTrack_Gmass
G4GLOB_DLL std::ostream G4cerr
static G4AntiKaonZero * AntiKaonZero()
G4DynamicParticle * PopProducts()
G4VDecayChannel * GetDecayChannel(G4int index) const
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
virtual G4DecayProducts * DecayIt(G4double mass=0.0)
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4KaonZero * KaonZero()
G4double GetFormationTime() const
G4double * theDaughterMass
G4double IntegrandFunction1(G4double xmass) const
G4double IntegrateCMMomentum2() const
G4KineticTrackVector * Decay()
G4KineticTrack & operator=(const G4KineticTrack &right)
void SetCreatorModelID(G4int id)
G4int GetCreatorModelID() const
G4LorentzVector the4Momentum
G4bool operator!=(const G4KineticTrack &right) const
G4ThreeVector thePosition
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4LorentzVector theTotal4Momentum
G4double * theActualWidth
G4int GetnChannels() const
const G4ThreeVector & GetPosition() const
const G4ParticleDefinition * GetDefinition() const
G4bool operator==(const G4KineticTrack &right) const
G4double theFormationTime
const G4LorentzVector & GetTrackingMomentum() const
G4double * theDaughterWidth
G4double theProjectilePotential
G4double IntegrateCMMomentum(const G4double lowerLimit) const
G4double IntegrandFunction3(G4double xmass) const
G4double IntegrandFunction4(G4double xmass) const
G4double IntegrandFunction2(G4double xmass) const
G4double EvaluateCMMomentum(const G4double mass, const G4double *m_ij) const
G4double BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
CascadeState theStateToNucleus
G4LorentzVector theFermi3Momentum
const G4LorentzVector & Get4Momentum() const
const G4ParticleDefinition * theDefinition
G4double GetActualMass() const
G4double EvaluateTotalActualWidth()
G4bool IsShortLived() const
G4double GetPDGMass() const
G4double GetPDGWidth() const
G4double GetPDGCharge() const
G4DecayTable * GetDecayTable() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
G4double GetMinimumMass(const G4ParticleDefinition *p) const
const G4String & GetParentName() const
G4int GetNumberOfDaughters() const
G4ParticleDefinition * GetDaughter(G4int anIndex)
const G4String & GetDaughterName(G4int anIndex) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool nucleon(G4int ityp)
static const G4LorentzVector zero(0., 0., 0., 0.)
void G4SwapObj(T *a, T *b)