G4INCL::ClusteringModelIntercomparison Class Reference

Cluster coalescence algorithm used in the IAEA intercomparison. More...

#include <G4INCLClusteringModelIntercomparison.hh>

Inheritance diagram for G4INCL::ClusteringModelIntercomparison:

Public Member Functions
	ClusteringModelIntercomparison (Config const *const theConfig)
virtual	~ClusteringModelIntercomparison ()
virtual Cluster *	getCluster (Nucleus *, Particle *)
virtual G4bool	clusterCanEscape (Nucleus const *const, Cluster const *const)
Data Structures
class	SortedNucleonConfiguration
	Class for storing and comparing sorted nucleon configurations.

---

Detailed Description

Cluster coalescence algorithm used in the IAEA intercomparison.

Definition at line 58 of file G4INCLClusteringModelIntercomparison.hh.

---

Constructor & Destructor Documentation

G4INCL::ClusteringModelIntercomparison::ClusteringModelIntercomparison

(

Config const *const

theConfig

)

[inline]

Definition at line 60 of file G4INCLClusteringModelIntercomparison.hh.

References G4INCL::ParticleTable::clusterZMax, G4INCL::ParticleTable::clusterZMin, and G4INCL::ParticleTable::maxClusterMass.

                                                                   :
      theNucleus(NULL),
      selectedA(0),
      selectedZ(0),
      sqtot(0.),
      cascadingEnergyPool(0.),
      protonMass(ParticleTable::getRealMass(Proton)),
      neutronMass(ParticleTable::getRealMass(Neutron)),
      runningMaxClusterAlgorithmMass(theConfig->getClusterMaxMass()),
      nConsideredMax(0),
      nConsidered(0),
      consideredPartners(NULL),
      isInRunningConfiguration(NULL),
      maxMassConfigurationSkipping(ParticleTable::maxClusterMass)
    {
      // Set up the maximum charge and neutron number for clusters
      clusterZMaxAll = 0;
      clusterNMaxAll = 0;
      for(G4int A=0; A<=runningMaxClusterAlgorithmMass; ++A) {
        if(ParticleTable::clusterZMax[A]>clusterZMaxAll)
          clusterZMaxAll = ParticleTable::clusterZMax[A];
        if(A-ParticleTable::clusterZMin[A]>clusterNMaxAll)
          clusterNMaxAll = A-ParticleTable::clusterZMin[A];
      }
      std::fill(candidateConfiguration,
                candidateConfiguration + ParticleTable::maxClusterMass,
                static_cast<Particle*>(NULL));

      std::fill(runningEnergies,
                runningEnergies + ParticleTable::maxClusterMass,
                0.0);

      std::fill(runningPotentials,
                runningPotentials + ParticleTable::maxClusterMass,
                0.0);

      std::fill(runningConfiguration,
                runningConfiguration + ParticleTable::maxClusterMass,
                -1);

    }

virtual G4INCL::ClusteringModelIntercomparison::~ClusteringModelIntercomparison

(

)

[inline, virtual]

Definition at line 102 of file G4INCLClusteringModelIntercomparison.hh.

                                              {
      delete [] consideredPartners;
      delete [] isInRunningConfiguration;
    }

---

Member Function Documentation

G4bool G4INCL::ClusteringModelIntercomparison::clusterCanEscape	(	Nucleus const *	const,
		Cluster const *	const
	)			[virtual]

Determine whether cluster can escape or not.

Implements G4INCL::IClusteringModel.

Definition at line 344 of file G4INCLClusteringModelIntercomparison.cc.

References G4INCL::ThreeVector::dot(), G4INCL::Particle::getA(), G4INCL::Particle::getMomentum(), G4INCL::Particle::getPosition(), G4INCL::ThreeVector::mag2(), and CLHEP::detail::n.

                                                                                                          {
    // Forbid emission of the whole nucleus
    if(c->getA()>=n->getA())
      return false;

    // Check the escape angle of the cluster
    const ThreeVector &pos = c->getPosition();
    const ThreeVector &mom = c->getMomentum();
    const G4double cosEscapeAngle = pos.dot(mom) / std::sqrt(pos.mag2()*mom.mag2());
    if(cosEscapeAngle < limitCosEscapeAngle)
      return false;

    return true;
  }

Cluster * G4INCL::ClusteringModelIntercomparison::getCluster	(	Nucleus *	,
		Particle *
	)			[virtual]

Choose a cluster candidate to be produced. At this point we don't yet decide if it can pass through the Coulomb barrier or not.

Implements G4INCL::IClusteringModel.

Definition at line 50 of file G4INCLClusteringModelIntercomparison.cc.

References G4INCL::ParticleTable::clusterPhaseSpaceCut, G4INCL::ParticleTable::clusterPosFact2, G4INCL::ThreeVector::dot(), G4INCL::Particle::getA(), G4INCL::Config::getClusterMaxMass(), G4INCL::Store::getConfig(), G4INCL::Nucleus::getDensity(), G4INCL::Particle::getEnergy(), G4INCL::Particle::getID(), G4INCL::Particle::getMomentum(), G4INCL::NuclearDensity::getNuclearRadius(), G4INCL::Store::getParticles(), G4INCL::Particle::getPosition(), G4INCL::Particle::getPotentialEnergy(), G4INCL::Nucleus::getStore(), G4INCL::Nucleus::getUniverseRadius(), G4INCL::Particle::getZ(), G4INCL::ThreeVector::mag(), and G4INCL::Particle::setPosition().

                                                                                          {
    // Set the maximum clustering mass dynamically, based on the current nucleus
    const G4int maxClusterAlgorithmMass = nucleus->getStore()->getConfig()->getClusterMaxMass();
    runningMaxClusterAlgorithmMass = std::min(maxClusterAlgorithmMass, nucleus->getA()/2);

    // Nucleus too small?
    if(runningMaxClusterAlgorithmMass<=1)
      return NULL;

    theNucleus = nucleus;
    Particle *theLeadingParticle = particle;

    // Initialise sqtot to a large number
    sqtot = 50000.0;
    selectedA = 0;
    selectedZ = 0;

    // The distance parameter, known as h in publications.
    // Default value is 1 fm.
    const G4double transp = 1.0;

    const G4double rmaxws = theNucleus->getUniverseRadius();

    // Radius of the sphere where the leading particle is positioned.
    const G4double Rprime = theNucleus->getDensity()->getNuclearRadius() + transp;

    // Bring the leading particle back to the coalescence sphere
    const G4double pk = theLeadingParticle->getMomentum().mag();
    const G4double cospr = theLeadingParticle->getPosition().dot(theLeadingParticle->getMomentum())/(theNucleus->getUniverseRadius() * pk);
    const G4double arg = rmaxws*rmaxws - Rprime*Rprime;
    G4double translat;

    if(arg > 0.0) {
      // coalescence sphere smaller than Rmax
      const G4double cosmin = std::sqrt(arg)/rmaxws;
      if(cospr <= cosmin) {
        // there is an intersection with the coalescence sphere
        translat = rmaxws * cospr;
      } else {
        // no intersection with the coalescence sphere
        translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin));
      }
    } else {
      // coalescence sphere larger than Rmax
      translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr));
    }

    const ThreeVector oldLeadingParticlePosition = theLeadingParticle->getPosition();
    const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->getMomentum() * (translat/pk);
    const ThreeVector &leadingParticleMomentum = theLeadingParticle->getMomentum();
    theLeadingParticle->setPosition(leadingParticlePosition);

    // Initialise the array of considered nucleons
    const G4int theNucleusA = theNucleus->getA();
    if(nConsideredMax < theNucleusA) {
      delete [] consideredPartners;
      delete [] isInRunningConfiguration;
      nConsideredMax = 2*theNucleusA;
      consideredPartners = new Particle *[nConsideredMax];
      isInRunningConfiguration = new G4bool [nConsideredMax];
      std::fill(isInRunningConfiguration,
                isInRunningConfiguration + nConsideredMax,
                false);
    }

    // Select the subset of nucleons that will be considered in the
    // cluster production:
    cascadingEnergyPool = 0.;
    nConsidered = 0;
    const ParticleList particles = theNucleus->getStore()->getParticles();
    for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
      if (!(*i)->isNucleon()) continue; // Only nucleons are allowed in clusters
      if ((*i)->getID() == theLeadingParticle->getID()) continue; // Don't count the leading particle

      G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2();
      G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2();
      G4double size = space*momentum*ParticleTable::clusterPosFact2[runningMaxClusterAlgorithmMass];
      // Nucleons are accepted only if they are "close enough" in phase space
      // to the leading nucleon. The selected phase-space parameter corresponds
      // to the running maximum cluster mass.
      if(size < ParticleTable::clusterPhaseSpaceCut[runningMaxClusterAlgorithmMass]) {
        consideredPartners[nConsidered] = *i;
        // Keep trace of how much energy is carried by cascading nucleons. This
        // is used to stop the clustering algorithm as soon as possible.
        if(!(*i)->isTargetSpectator())
          cascadingEnergyPool += (*i)->getEnergy() - (*i)->getPotentialEnergy() - 931.3;
        nConsidered++;
        // Make sure we don't exceed the array size
// assert(nConsidered<=nConsideredMax);
      }
    }
    // Sort the list of considered partners so that we give priority
    // to participants. As soon as we encounter the first spectator in
    // the list we know that all the remaining nucleons will be
    // spectators too.
    std::partition(consideredPartners, consideredPartners+nConsidered, cascadingFirstPredicate);

    // Clear the sets of checked configurations
    // We stop caching two masses short of the max mass -- there seems to be a
    // performance hit above
    maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2;
    for(G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i) // no caching for A=1,2
      checkedConfigurations[i].clear();

    // Initialise position, momentum and energy of the running cluster
    // configuration
    runningPositions[1] = leadingParticlePosition;
    runningMomenta[1] = leadingParticleMomentum;
    runningEnergies[1] = theLeadingParticle->getEnergy();
    runningPotentials[1] = theLeadingParticle->getPotentialEnergy();

    // Make sure that all the elements of isInRunningConfiguration are false.
// assert(std::count(isInRunningConfiguration, isInRunningConfiguration+nConsidered, true)==0);

    // Start the cluster search!
    findClusterStartingFrom(1, theLeadingParticle->getZ());

    // Again, make sure that all the elements of isInRunningConfiguration have
    // been reset to false. This is a sanity check.
// assert(std::count(isInRunningConfiguration, isInRunningConfiguration+nConsidered, true)==0);

    Cluster *chosenCluster = 0;
    if(selectedA!=0) { // A cluster was found!
      candidateConfiguration[selectedA-1] = theLeadingParticle;
      chosenCluster =  new Cluster(candidateConfiguration,
          candidateConfiguration + selectedA);
    }

    // Restore the original position of the leading particle
    theLeadingParticle->setPosition(oldLeadingParticlePosition);

    return chosenCluster;
  }

---

The documentation for this class was generated from the following files:

• G4INCLClusteringModelIntercomparison.hh
• G4INCLClusteringModelIntercomparison.cc

---