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)
Public Member Functions inherited from G4INCL::IClusteringModel
	IClusteringModel ()
virtual	~IClusteringModel ()
virtual Cluster *	getCluster (Nucleus *, Particle *)=0
virtual G4bool	clusterCanEscape (Nucleus const *const, Cluster const *const)=0

Detailed Description

Cluster coalescence algorithm used in the IAEA intercomparison.

Definition at line 60 of file G4INCLClusteringModelIntercomparison.hh.

Constructor & Destructor Documentation

ClusteringModelIntercomparison()

G4INCL::ClusteringModelIntercomparison::ClusteringModelIntercomparison ( Config const *const  theConfig )

inline

Definition at line 62 of file G4INCLClusteringModelIntercomparison.hh.

                                                                   :
      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);
 
    }

~ClusteringModelIntercomparison()

virtual G4INCL::ClusteringModelIntercomparison::~ClusteringModelIntercomparison ( )

inlinevirtual

Definition at line 104 of file G4INCLClusteringModelIntercomparison.hh.

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

Member Function Documentation

clusterCanEscape()

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

virtual

Determine whether cluster can escape or not.

Implements G4INCL::IClusteringModel.

Definition at line 346 of file G4INCLClusteringModelIntercomparison.cc.

                                                                                                          {
    // 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;
  }

getCluster()

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 52 of file G4INCLClusteringModelIntercomparison.cc.

                                                                                          {
    // 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