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.
82 const G4int maxClusterAlgorithmMass = nucleus->getStore()->getConfig()->getClusterMaxMass();
83 runningMaxClusterAlgorithmMass =
std::min(maxClusterAlgorithmMass, nucleus->getA()/2);
86 if(runningMaxClusterAlgorithmMass<=1)
90 Particle *theLeadingParticle = particle;
107 const G4double pk = theLeadingParticle->getMomentum().mag();
108 const G4double cospr = theLeadingParticle->getPosition().dot(theLeadingParticle->getMomentum())/(theNucleus->
getUniverseRadius() * pk);
109 const G4double arg = rmaxws*rmaxws - Rprime*Rprime;
114 const G4double cosmin = std::sqrt(arg)/rmaxws;
115 if(cospr <= cosmin) {
117 translat = rmaxws * cospr;
120 translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin));
124 translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr));
127 const ThreeVector oldLeadingParticlePosition = theLeadingParticle->getPosition();
128 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->getMomentum() * (translat/pk);
129 const ThreeVector &leadingParticleMomentum = theLeadingParticle->getMomentum();
130 theLeadingParticle->setPosition(leadingParticlePosition);
133 const G4int theNucleusA = theNucleus->
getA();
134 if(nConsideredMax < theNucleusA) {
135 delete [] consideredPartners;
136 delete [] isInRunningConfiguration;
137 nConsideredMax = 2*theNucleusA;
138 consideredPartners =
new ConsideredPartner[nConsideredMax];
139 isInRunningConfiguration =
new G4bool [nConsideredMax];
140 std::fill(isInRunningConfiguration,
141 isInRunningConfiguration + nConsideredMax,
147 cascadingEnergyPool = 0.;
150 for(
ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
151 if (!(*i)->isNucleon())
continue;
152 if ((*i)->getID() == theLeadingParticle->getID())
continue;
154 G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2();
155 G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2();
156 G4double size = space*momentum*clusterPosFact2[runningMaxClusterAlgorithmMass];
160 if(size < clusterPhaseSpaceCut[runningMaxClusterAlgorithmMass]) {
161 consideredPartners[nConsidered] = *i;
164 if(!(*i)->isTargetSpectator())
165 cascadingEnergyPool += consideredPartners[nConsidered].
energy - consideredPartners[nConsidered].potentialEnergy - 931.3;
177 #ifndef INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None
181 maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2;
182 for(
G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i)
183 checkedConfigurations[i].clear();
188 runningPositions[1] = leadingParticlePosition;
189 runningMomenta[1] = leadingParticleMomentum;
190 runningEnergies[1] = theLeadingParticle->getEnergy();
191 runningPotentials[1] = theLeadingParticle->getPotentialEnergy();
197 findClusterStartingFrom(1, theLeadingParticle->getZ());
203 Cluster *chosenCluster = 0;
205 candidateConfiguration[selectedA-1] = theLeadingParticle;
206 chosenCluster =
new Cluster(candidateConfiguration,
207 candidateConfiguration + selectedA);
211 theLeadingParticle->setPosition(oldLeadingParticlePosition);
213 return chosenCluster;
G4int getA() const
Returns the baryon number.
ParticleList const & getParticles() const
G4double getProtonNuclearRadius() const
NuclearDensity const * getDensity() const
Getter for theDensity.
G4double energy(const ThreeVector &p, const G4double m)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double getUniverseRadius() const
Getter for theUniverseRadius.
ParticleList::const_iterator ParticleIter