33 #define INCLXX_IN_GEANT4_MODE 1
45 const G4int ClusteringModelIntercomparison::clusterZMin[
ParticleTable::maxClusterMass+1] = {0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3};
46 const G4int ClusteringModelIntercomparison::clusterZMax[
ParticleTable::maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8};
60 0.0082645, 0.0069444};
69 const G4double ClusteringModelIntercomparison::limitCosEscapeAngle = 0.7;
71 #ifndef INCLXX_IN_GEANT4_MODE
73 G4bool cascadingFirstPredicate(ConsideredPartner
const &aPartner) {
74 return !aPartner.isTargetSpectator;
82 runningMaxClusterAlgorithmMass =
std::min(maxClusterAlgorithmMass, nucleus->
getA()/2);
85 if(runningMaxClusterAlgorithmMass<=1)
89 Particle *theLeadingParticle = particle;
108 const G4double arg = rmaxws*rmaxws - Rprime*Rprime;
113 const G4double cosmin = std::sqrt(arg)/rmaxws;
114 if(cospr <= cosmin) {
116 translat = rmaxws * cospr;
119 translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin));
123 translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr));
127 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->
getMomentum() * (translat/pk);
129 theLeadingParticle->
setPosition(leadingParticlePosition);
132 const G4int theNucleusA = theNucleus->
getA();
133 if(nConsideredMax < theNucleusA) {
134 delete [] consideredPartners;
135 delete [] isInRunningConfiguration;
136 nConsideredMax = 2*theNucleusA;
138 isInRunningConfiguration =
new G4bool [nConsideredMax];
139 std::fill(isInRunningConfiguration,
140 isInRunningConfiguration + nConsideredMax,
146 cascadingEnergyPool = 0.;
149 for(
ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
150 if (!(*i)->isNucleon())
continue;
151 if ((*i)->getID() == theLeadingParticle->
getID())
continue;
153 G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2();
154 G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2();
155 G4double size = space*momentum*clusterPosFact2[runningMaxClusterAlgorithmMass];
159 if(size < clusterPhaseSpaceCut[runningMaxClusterAlgorithmMass]) {
160 consideredPartners[nConsidered] = *i;
163 if(!(*i)->isTargetSpectator())
164 cascadingEnergyPool += consideredPartners[nConsidered].
energy - consideredPartners[nConsidered].potentialEnergy - 931.3;
179 maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2;
180 for(
G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i)
181 checkedConfigurations[i].clear();
185 runningPositions[1] = leadingParticlePosition;
186 runningMomenta[1] = leadingParticleMomentum;
187 runningEnergies[1] = theLeadingParticle->
getEnergy();
194 findClusterStartingFrom(1, theLeadingParticle->
getZ());
202 candidateConfiguration[selectedA-1] = theLeadingParticle;
203 chosenCluster =
new Cluster(candidateConfiguration,
204 candidateConfiguration + selectedA);
208 theLeadingParticle->
setPosition(oldLeadingParticlePosition);
210 return chosenCluster;
215 const G4double psMomentum = (p.
momentum*oldA - runningMomenta[oldA]).mag2();
216 return psSpace * psMomentum * clusterPosFact2[oldA + 1];
219 void ClusteringModelIntercomparison::findClusterStartingFrom(
const G4int oldA,
const G4int oldZ) {
220 const G4int newA = oldA + 1;
221 const G4int oldAMinusOne = oldA - 1;
226 const G4double phaseSpaceCut = clusterPhaseSpaceCut[newA];
229 const G4bool cachingEnabled = (newA<=maxMassConfigurationSkipping && newA>=3);
232 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
233 HashContainer *theHashContainer;
235 theHashContainer = &(checkedConfigurations[oldA-2]);
237 theHashContainer = NULL;
238 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
239 SortedNucleonConfigurationContainer *theConfigContainer;
241 theConfigContainer = &(checkedConfigurations[oldA-2]);
243 theConfigContainer = NULL;
245 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask.
249 const G4int ZMinForNewA = clusterZMin[newA];
250 const G4int ZMaxForNewA = clusterZMax[newA];
252 for(
G4int i=0; i<nConsidered; ++i) {
254 if(isInRunningConfiguration[i])
continue;
256 ConsideredPartner
const &candidateNucleon = consideredPartners[i];
259 newZ = oldZ + candidateNucleon.
Z;
263 if(newZ > clusterZMaxAll || newN > clusterNMaxAll)
269 const G4double phaseSpace = getPhaseSpace(oldA, candidateNucleon);
270 if(phaseSpace > phaseSpaceCut)
continue;
273 runningConfiguration[oldAMinusOne] = i;
274 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
275 Hashing::HashType configHash;
276 HashIterator aHashIter;
277 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
278 SortedNucleonConfiguration thisConfig;
279 SortedNucleonConfigurationIterator thisConfigIter;
282 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
283 configHash = Hashing::hashConfig(runningConfiguration, oldA);
284 aHashIter = theHashContainer->lower_bound(configHash);
286 if(aHashIter!=theHashContainer->end()
287 && !(configHash < *aHashIter))
289 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
290 thisConfig.fill(runningConfiguration,oldA);
291 thisConfigIter = theConfigContainer->lower_bound(thisConfig);
293 if(thisConfigIter!=theConfigContainer->end()
294 && !(thisConfig < *thisConfigIter))
300 runningEnergies[newA] = runningEnergies[oldA] + candidateNucleon.energy;
302 runningPotentials[newA] = runningPotentials[oldA] + candidateNucleon.potentialEnergy;
305 G4double oldCascadingEnergyPool = cascadingEnergyPool;
306 if(!candidateNucleon.isTargetSpectator)
307 cascadingEnergyPool -= candidateNucleon.energy - candidateNucleon.potentialEnergy - 931.3;
312 const G4double halfB = 0.72 * newZ *
314 const G4double tout = runningEnergies[newA] - runningPotentials[newA] -
316 if(tout<=halfB && tout+cascadingEnergyPool<=halfB) {
317 cascadingEnergyPool = oldCascadingEnergyPool;
322 runningPositions[newA] = (runningPositions[oldA] * oldA + candidateNucleon.position)*clusterPosFact[newA];
323 runningMomenta[newA] = runningMomenta[oldA] + candidateNucleon.momentum;
327 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
328 theHashContainer->insert(aHashIter, configHash);
329 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
330 theConfigContainer->insert(thisConfigIter, thisConfig);
335 isInRunningConfiguration[i] =
true;
338 if(newZ >= ZMinForNewA && newZ <= ZMaxForNewA) {
341 runningMomenta[newA]);
342 const G4double sqct = (sqc - 2.*newZ*protonMass - 2.*(newA-newZ)*neutronMass
344 *clusterPosFact[newA];
354 for(
G4int j=0; j<oldA; ++j)
355 candidateConfiguration[j] = consideredPartners[runningConfiguration[j]].particle;
363 if(newA < runningMaxClusterAlgorithmMass && newA+1 < theNucleus->
getA()) {
364 findClusterStartingFrom(newA, newZ);
368 isInRunningConfiguration[i] =
false;
369 cascadingEnergyPool = oldCascadingEnergyPool;
382 if(cosEscapeAngle < limitCosEscapeAngle)
G4int getA() const
Returns the baryon number.
virtual G4bool clusterCanEscape(Nucleus const *const, Cluster const *const)
G4double dot(const ThreeVector &v) const
ParticleList const & getParticles() const
G4int getClusterMaxMass() const
Get the maximum mass for production of clusters.
const G4INCL::ThreeVector & getMomentum() const
Config const * getConfig()
G4double getEnergy() const
G4double getPotentialEnergy() const
Get the particle potential energy.
G4double getProtonNuclearRadius() const
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
double precision function energy(A, Z)
const G4int maxClusterMass
virtual void setPosition(const G4INCL::ThreeVector &position)
G4int getZ() const
Returns the charge number.
G4double invariantMass(const G4double E, const ThreeVector &p)
Container for the relevant information.
virtual Cluster * getCluster(Nucleus *, Particle *)
NuclearDensity const * getDensity() const
Getter for theDensity.
const G4INCL::ThreeVector & getPosition() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double getUniverseRadius() const
Getter for theUniverseRadius.
Functions for hashing a collection of NucleonItems.
ParticleList::const_iterator ParticleIter