G4INCL::ClusterDecay::anonymous_namespace{G4INCLClusterDecay.cc} Namespace Reference

Functions
void	phaseSpaceDecay (Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
	Disassembles unbound nuclei using the phase-space model. More...
void	recursiveDecay (Cluster *const c, ParticleList *decayProducts)
	Recursively decay clusters. More...
void	threeBodyDecay (Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
	Carries out three-body decays. More...
void	twoBodyDecay (Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
	Carries out two-body decays. More...

Function Documentation

phaseSpaceDecay()

void G4INCL::ClusterDecay::anonymous_namespace{G4INCLClusterDecay.cc}::phaseSpaceDecay	(	Cluster *const	c,
		ClusterDecayType	theDecayMode,
		ParticleList *	decayProducts
	)

Disassembles unbound nuclei using the phase-space model.

This implementation uses the Kopylov algorithm, defined in namespace PhaseSpaceGenerator.

Definition at line 405 of file G4INCLClusterDecay.cc.

                                                                                                          {
        const G4int theA = c->getA();
        const G4int theZ = c->getZ();
        const G4int theL = (-1)*(c->getS());
        const ThreeVector mom(0.0, 0.0, 0.0);
        const ThreeVector pos = c->getPosition();
 
        G4int theZStep;
        
        ParticleType theEjectileType;
        switch(theDecayMode) {
          case ProtonUnbound:
            theZStep = 1;
            theEjectileType = Proton;
            break;
          case NeutronUnbound:
            theZStep = 0;
            theEjectileType = Neutron;
            break;
          case LambdaUnbound: // Will always completly decay. Append only if theA == 0 and/or theZ == 0
            theZStep = -99;
            if(theZ==0) theEjectileType = Neutron;
            else theEjectileType = Proton; 
            break;
          default:
            INCL_ERROR("Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << '\n'
                  << c->print());
            return;
        }
        
        // Find the daughter cluster (first cluster which is not
        // proton/neutron-unbound, in the sense of the table)
        G4int finalDaughterZ, finalDaughterA, finalDaughterL;
        if(theZ<ParticleTable::clusterTableZSize && theA<ParticleTable::clusterTableASize && theZStep != -99) {
          finalDaughterZ=theZ;
          finalDaughterA=theA;
          finalDaughterL=theL;
          while(finalDaughterA>0 && clusterDecayMode[finalDaughterL][finalDaughterZ][finalDaughterA]!=StableCluster) { /* Loop modified, 15.01.18, J. Hirtz */
            finalDaughterA--;
            finalDaughterZ -= theZStep;
          }
        } else {
          finalDaughterA = 1;
          if(theDecayMode==ProtonUnbound){
            finalDaughterZ = 1;
            finalDaughterL = 0;
          }
          else if(theDecayMode==NeutronUnbound){
            finalDaughterZ = 0;
            finalDaughterL = 0;
          }
          else {
            finalDaughterZ = 0;
            finalDaughterL = 1;
          }
        }
// assert(finalDaughterZ<=theZ && finalDaughterA<theA && finalDaughterA>0 && finalDaughterZ>=0 && finalDaughterL>=0);
 
        // Compute the available decay energy
        const G4int nLambda = theL-finalDaughterL;
        const G4int nSplits = theA-finalDaughterA-nLambda;
        // c->getMass() can possibly contain some excitation energy, too
        const G4double availableEnergy = c->getMass();
 
        // Save the boost vector for the cluster
        const ThreeVector boostVector = - c->boostVector();
 
        // Create a list of decay products
        ParticleList products;
        c->setA(finalDaughterA);
        c->setZ(finalDaughterZ);
        c->setS((-1)*finalDaughterL);
        c->setRealMass();
        c->setMomentum(ThreeVector());
        c->adjustEnergyFromMomentum();
        products.push_back(c);
        
        for(G4int j=0; j<nLambda; ++j) {
          Particle *ejectile = new Particle(Lambda, mom, pos);
          ejectile->setRealMass();
          products.push_back(ejectile);
        }
        for(G4int i=0; i<nSplits; ++i) {
          Particle *ejectile = new Particle(theEjectileType, mom, pos);
          ejectile->setRealMass();
          products.push_back(ejectile);
        }
 
        PhaseSpaceGenerator::generate(availableEnergy, products);
        products.boost(boostVector);
 
        // Copy decay products in the output list (but skip the residue)
        ParticleList::iterator productsIter = products.begin();
        std::advance(productsIter, 1);
        decayProducts->insert(decayProducts->end(), productsIter, products.end());
 
        c->setExcitationEnergy(0.);
      }

References G4INCL::Particle::adjustEnergyFromMomentum(), G4INCL::ParticleList::boost(), G4INCL::Particle::boostVector(), G4INCL::ClusterDecay::clusterDecayMode, G4INCL::ParticleTable::clusterTableASize, G4INCL::ParticleTable::clusterTableZSize, G4INCL::PhaseSpaceGenerator::generate(), G4INCL::Particle::getA(), G4INCL::Particle::getMass(), G4INCL::Particle::getPosition(), G4INCL::Particle::getS(), G4INCL::Particle::getZ(), INCL_ERROR, G4INCL::Lambda, G4INCL::ClusterDecay::LambdaUnbound, G4INCL::Neutron, G4INCL::ClusterDecay::NeutronUnbound, pos, G4INCL::Cluster::print(), G4INCL::Proton, G4INCL::ClusterDecay::ProtonUnbound, G4INCL::Cluster::setA(), G4INCL::Cluster::setExcitationEnergy(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setRealMass(), G4INCL::Cluster::setS(), G4INCL::Cluster::setZ(), and G4INCL::ClusterDecay::StableCluster.

Referenced by recursiveDecay().

recursiveDecay()

void G4INCL::ClusterDecay::anonymous_namespace{G4INCLClusterDecay.cc}::recursiveDecay	(	Cluster *const	c,
		ParticleList *	decayProducts
	)

Recursively decay clusters.

Parameters

c	cluster that should decay
decayProducts	decay products are appended to the end of this list

Definition at line 509 of file G4INCLClusterDecay.cc.

                                                                          {
        const G4int Z = c->getZ();
        const G4int A = c->getA();
        const G4int S = c->getS();
// assert(c->getExcitationEnergy()>-1.e-5);
        if(c->getExcitationEnergy()<0.)
          c->setExcitationEnergy(0.);
 
        if(Z<ParticleTable::clusterTableZSize && A<ParticleTable::clusterTableASize && (S*(-1))<ParticleTable::clusterTableSSize) {
          ClusterDecayType theDecayMode = clusterDecayMode[(S*(-1))][Z][A];
 
          switch(theDecayMode) {
            default:
              INCL_ERROR("Unrecognized cluster-decay mode: " << theDecayMode << '\n'
                    << c->print());
              return;
              break;
            case StableCluster:
              // For stable clusters, just return
              return;
              break;
            case ProtonDecay:
            case NeutronDecay:
            case LambdaDecay:
            case AlphaDecay:
              // Two-body decays
              twoBodyDecay(c, theDecayMode, decayProducts);
              break;
            case TwoProtonDecay:
            case TwoNeutronDecay:
              // Three-body decays
              threeBodyDecay(c, theDecayMode, decayProducts);
              break;
            case ProtonUnbound:
            case NeutronUnbound:
            case LambdaUnbound:
              // Phase-space decays
              phaseSpaceDecay(c, theDecayMode, decayProducts);
              break;
          }
 
          // Calls itself recursively in case the produced remnant is still unstable.
          // Sneaky, isn't it.
          recursiveDecay(c,decayProducts);
 
        } else {
          // The cluster is too large for our decay-mode table. Decompose it only
          // if Z==0 || Z==A.
          INCL_DEBUG("Cluster is outside the decay-mode table." << c->print() << '\n');
          if(Z==A) {
            INCL_DEBUG("Z==A, will decompose it in free protons." << '\n');
            phaseSpaceDecay(c, ProtonUnbound, decayProducts);
          } else if(Z==0) {
            INCL_DEBUG("Z==0, will decompose it in free neutrons." << '\n');
            phaseSpaceDecay(c, NeutronUnbound, decayProducts);
          }
        }
      }

References A, G4INCL::ClusterDecay::AlphaDecay, G4INCL::ClusterDecay::clusterDecayMode, G4INCL::ParticleTable::clusterTableASize, G4INCL::ParticleTable::clusterTableSSize, G4INCL::ParticleTable::clusterTableZSize, G4INCL::Particle::getA(), G4INCL::Cluster::getExcitationEnergy(), G4INCL::Particle::getS(), G4INCL::Particle::getZ(), INCL_DEBUG, INCL_ERROR, G4INCL::ClusterDecay::LambdaDecay, G4INCL::ClusterDecay::LambdaUnbound, G4INCL::ClusterDecay::NeutronDecay, G4INCL::ClusterDecay::NeutronUnbound, phaseSpaceDecay(), G4INCL::Cluster::print(), G4INCL::ClusterDecay::ProtonDecay, G4INCL::ClusterDecay::ProtonUnbound, recursiveDecay(), S(), G4INCL::Cluster::setExcitationEnergy(), G4INCL::ClusterDecay::StableCluster, threeBodyDecay(), twoBodyDecay(), G4INCL::ClusterDecay::TwoNeutronDecay, G4INCL::ClusterDecay::TwoProtonDecay, and Z.

Referenced by G4INCL::ClusterDecay::decay(), and recursiveDecay().

threeBodyDecay()

void G4INCL::ClusterDecay::anonymous_namespace{G4INCLClusterDecay.cc}::threeBodyDecay	(	Cluster *const	c,
		ClusterDecayType	theDecayMode,
		ParticleList *	decayProducts
	)

Carries out three-body decays.

Definition at line 137 of file G4INCLClusterDecay.cc.

                                                                                                         {
        Particle *decayParticle1 = 0;
        Particle *decayParticle2 = 0;
        const ThreeVector mom(0.0, 0.0, 0.0);
        const ThreeVector pos = c->getPosition();
 
        // Create the emitted particles
        switch(theDecayMode) {
          case TwoProtonDecay:
            decayParticle1 = new Particle(Proton, mom, pos);
            decayParticle2 = new Particle(Proton, mom, pos);
            break;
          case TwoNeutronDecay:
            decayParticle1 = new Particle(Neutron, mom, pos);
            decayParticle2 = new Particle(Neutron, mom, pos);
            break;
          default:
            INCL_ERROR("Unrecognized cluster-decay mode in three-body decay: " << theDecayMode << '\n'
                  << c->print());
            return;
        }
        decayParticle1->makeParticipant();
        decayParticle2->makeParticipant();
        decayParticle1->setNumberOfDecays(1);
        decayParticle2->setNumberOfDecays(1);
        decayParticle1->setRealMass();
        decayParticle2->setRealMass();
 
        // Save some variables of the mother cluster
        const G4double motherMass = c->getMass();
        const ThreeVector velocity = -c->boostVector();
 
        // Masses and charges of the daughter particle and of the decay products
        const G4int decayZ1 = decayParticle1->getZ();
        const G4int decayA1 = decayParticle1->getA();
        const G4int decayS1 = decayParticle1->getS();
        const G4int decayZ2 = decayParticle2->getZ();
        const G4int decayA2 = decayParticle2->getA();
        const G4int decayS2 = decayParticle2->getS();
        const G4int decayZ = decayZ1 + decayZ2;
        const G4int decayA = decayA1 + decayA2;
        const G4int decayS = decayS1 + decayS2;
        const G4int daughterZ = c->getZ() - decayZ;
        const G4int daughterA = c->getA() - decayA;
        const G4int daughterS = c->getS() - decayS;
        const G4double decayMass1 = decayParticle1->getMass();
        const G4double decayMass2 = decayParticle2->getMass();
        const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ,daughterS);
 
        // Q-values
        G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
// assert(qValue > -1e-5); // Q-value should be >0
        if(qValue<0.)
          qValue=0.;
        const G4double qValueB = qValue * Random::shoot();
 
        // The decay particles behave as if they had more mass until the second
        // decay
        const G4double decayMass = decayMass1 + decayMass2 + qValueB;
 
        /* Stage A: mother --> daughter + (decay1+decay2) */
 
        // The mother cluster becomes the daughter
        c->setZ(daughterZ);
        c->setA(daughterA);
        c->setS(daughterS);
        c->setMass(daughterMass);
        c->setExcitationEnergy(0.);
 
        // Decay kinematics in the mother rest frame
        const G4double pCMA = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
        const ThreeVector momentumA = Random::normVector(pCMA);
        c->setMomentum(momentumA);
        c->adjustEnergyFromMomentum();
        const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.mag2());
 
        /* Stage B: (decay1+decay2) --> decay1 + decay2 */
 
        // Decay kinematics in the (decay1+decay2) rest frame
        const G4double pCMB = KinematicsUtils::momentumInCM(decayMass, decayMass1, decayMass2);
        const ThreeVector momentumB = Random::normVector(pCMB);
        decayParticle1->setMomentum(momentumB);
        decayParticle2->setMomentum(-momentumB);
        decayParticle1->adjustEnergyFromMomentum();
        decayParticle2->adjustEnergyFromMomentum();
 
        // Boost decay1 and decay2 to the Stage-A decay frame
        decayParticle1->boost(decayBoostVector);
        decayParticle2->boost(decayBoostVector);
 
        // Boost all particles to the lab frame
        decayParticle1->boost(velocity);
        decayParticle2->boost(velocity);
        c->boost(velocity);
 
        // Add the decay particles to the list of decay products
        decayProducts->push_back(decayParticle1);
        decayProducts->push_back(decayParticle2);
      }

References G4INCL::Particle::adjustEnergyFromMomentum(), G4INCL::Cluster::boost(), G4INCL::Particle::boost(), G4INCL::Particle::boostVector(), G4INCL::Particle::getA(), G4INCL::Particle::getMass(), G4INCL::Particle::getPosition(), G4INCL::ParticleTable::getRealMass(), G4INCL::Particle::getS(), G4INCL::Particle::getZ(), INCL_ERROR, G4INCL::ThreeVector::mag2(), G4INCL::Particle::makeParticipant(), G4INCL::KinematicsUtils::momentumInCM(), G4INCL::Neutron, G4INCL::Random::normVector(), pos, G4INCL::Cluster::print(), G4INCL::Proton, G4INCL::Cluster::setA(), G4INCL::Cluster::setExcitationEnergy(), G4INCL::Particle::setMass(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setNumberOfDecays(), G4INCL::Particle::setRealMass(), G4INCL::Cluster::setS(), G4INCL::Cluster::setZ(), G4INCL::Random::shoot(), G4INCL::ClusterDecay::TwoNeutronDecay, and G4INCL::ClusterDecay::TwoProtonDecay.

Referenced by recursiveDecay().

twoBodyDecay()

void G4INCL::ClusterDecay::anonymous_namespace{G4INCLClusterDecay.cc}::twoBodyDecay	(	Cluster *const	c,
		ClusterDecayType	theDecayMode,
		ParticleList *	decayProducts
	)

Carries out two-body decays.

Definition at line 60 of file G4INCLClusterDecay.cc.

                                                                                                       {
        Particle *decayParticle = 0;
        const ThreeVector mom(0.0, 0.0, 0.0);
        const ThreeVector pos = c->getPosition();
 
        // Create the emitted particle
        switch(theDecayMode) {
          case ProtonDecay:
            decayParticle = new Particle(Proton, mom, pos);
            break;
          case NeutronDecay:
            decayParticle = new Particle(Neutron, mom, pos);
            break;
          case AlphaDecay:
            decayParticle = new Cluster(2,4,0,false);
            break;
          case LambdaDecay:
            decayParticle = new Particle(Lambda, mom, pos);
            break;
          default:
            INCL_ERROR("Unrecognized cluster-decay mode in two-body decay: " << theDecayMode << '\n'
                  << c->print());
            return;
        }
        decayParticle->makeParticipant();
        decayParticle->setNumberOfDecays(1);
        decayParticle->setPosition(c->getPosition());
        decayParticle->setEmissionTime(c->getEmissionTime());
        decayParticle->setRealMass();
 
        // Save some variables of the mother cluster
#ifdef INCLXX_IN_GEANT4_MODE
          if ((c->getZ() == 1) && (c->getA() == 2) && (c->getS() == -1)) { // no Mass for A=2,Z=1,S=-1 in Geant4
              c->setMass(2053.952);
              if (c->getEnergy() < 2053.952)  // Energy can be lower than the sum of p and Lambda masses (2053.952)...
                  c->setMomentum(c->getMomentum() * 0.) ;
              else
                  c->setMomentum(c->getMomentum() / (std::sqrt(c->getMomentum().mag2())/std::sqrt(c->getMomentum().mag2() - 2053.952*2053.952))) ;
          }
#endif
        G4double motherMass = c->getMass();
        const ThreeVector velocity = -c->boostVector();
 
        // Characteristics of the daughter particle
        const G4int daughterZ = c->getZ() - decayParticle->getZ();
        const G4int daughterA = c->getA() - decayParticle->getA();
        const G4int daughterS = c->getS() - decayParticle->getS();
        const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ,daughterS);
        
        // The mother cluster becomes the daughter
        c->setZ(daughterZ);
        c->setA(daughterA);
        c->setS(daughterS);
        c->setMass(daughterMass);
        c->setExcitationEnergy(0.);
 
        // Decay kinematics in the mother rest frame
        const G4double decayMass = decayParticle->getMass();
// assert(motherMass-daughterMass-decayMass>-1.e-5); // Q-value should be >0
        G4double pCM = 0.;
        if(motherMass-daughterMass-decayMass>0.)
          pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
        const ThreeVector momentum = Random::normVector(pCM);
        c->setMomentum(momentum);
        c->adjustEnergyFromMomentum();
        decayParticle->setMomentum(-momentum);
        decayParticle->adjustEnergyFromMomentum();
 
        // Boost to the lab frame
        decayParticle->boost(velocity);
        c->boost(velocity);
 
        // Add the decay particle to the list of decay products
        decayProducts->push_back(decayParticle);
      }

References G4INCL::Particle::adjustEnergyFromMomentum(), G4INCL::ClusterDecay::AlphaDecay, G4INCL::Cluster::boost(), G4INCL::Particle::boost(), G4INCL::Particle::boostVector(), G4INCL::Particle::getA(), G4INCL::Particle::getEmissionTime(), G4INCL::Particle::getEnergy(), G4INCL::Particle::getMass(), G4INCL::Particle::getMomentum(), G4INCL::Particle::getPosition(), G4INCL::ParticleTable::getRealMass(), G4INCL::Particle::getS(), G4INCL::Particle::getZ(), INCL_ERROR, G4INCL::Lambda, G4INCL::ClusterDecay::LambdaDecay, G4INCL::ThreeVector::mag2(), G4INCL::Particle::makeParticipant(), G4INCL::KinematicsUtils::momentumInCM(), G4INCL::Neutron, G4INCL::ClusterDecay::NeutronDecay, G4INCL::Random::normVector(), pos, G4INCL::Cluster::print(), G4INCL::Proton, G4INCL::ClusterDecay::ProtonDecay, G4INCL::Cluster::setA(), G4INCL::Particle::setEmissionTime(), G4INCL::Cluster::setExcitationEnergy(), G4INCL::Particle::setMass(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setNumberOfDecays(), G4INCL::Particle::setPosition(), G4INCL::Particle::setRealMass(), G4INCL::Cluster::setS(), and G4INCL::Cluster::setZ().

Referenced by recursiveDecay().