34 #define INCLXX_IN_GEANT4_MODE 1
46 const G4int ClusteringModelIntercomparison::clusterZMin[
ParticleTable::maxClusterMass+1] = {0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3};
47 const G4int ClusteringModelIntercomparison::clusterZMax[
ParticleTable::maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8};
61 0.0082645, 0.0069444};
70 const G4double ClusteringModelIntercomparison::limitCosEscapeAngle = 0.7;
72 #ifdef INCL_DO_NOT_COMPILE
74 G4bool cascadingFirstPredicate(ConsideredPartner
const &aPartner) {
75 return !aPartner.isTargetSpectator;
83 runningMaxClusterAlgorithmMass =
std::min(maxClusterAlgorithmMass, nucleus->
getA()/2);
86 if(runningMaxClusterAlgorithmMass<=1)
90 Particle *theLeadingParticle = particle;
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));
128 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->
getMomentum() * (translat/pk);
130 theLeadingParticle->
setPosition(leadingParticlePosition);
133 const G4int theNucleusA = theNucleus->
getA();
134 if(nConsideredMax < theNucleusA) {
135 delete [] consideredPartners;
136 delete [] isInRunningConfiguration;
137 nConsideredMax = 2*theNucleusA;
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();
197 findClusterStartingFrom(1, theLeadingParticle->
getZ());
205 candidateConfiguration[selectedA-1] = theLeadingParticle;
206 chosenCluster =
new Cluster(candidateConfiguration,
207 candidateConfiguration + selectedA);
211 theLeadingParticle->
setPosition(oldLeadingParticlePosition);
213 return chosenCluster;
218 const G4double psMomentum = (p.
momentum*oldA - runningMomenta[oldA]).mag2();
219 return psSpace * psMomentum * clusterPosFact2[oldA + 1];
222 void ClusteringModelIntercomparison::findClusterStartingFrom(
const G4int oldA,
const G4int oldZ) {
223 const G4int newA = oldA + 1;
224 const G4int oldAMinusOne = oldA - 1;
229 const G4double phaseSpaceCut = clusterPhaseSpaceCut[newA];
232 const G4bool cachingEnabled = (newA<=maxMassConfigurationSkipping && newA>=3);
235 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
236 HashContainer *theHashContainer;
238 theHashContainer = &(checkedConfigurations[oldA-2]);
240 theHashContainer = NULL;
241 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
242 SortedNucleonConfigurationContainer *theConfigContainer;
244 theConfigContainer = &(checkedConfigurations[oldA-2]);
246 theConfigContainer = NULL;
247 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None)
248 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None.
252 const G4int ZMinForNewA = clusterZMin[newA];
253 const G4int ZMaxForNewA = clusterZMax[newA];
255 for(
G4int i=0; i<nConsidered; ++i) {
257 if(isInRunningConfiguration[i])
continue;
259 ConsideredPartner
const &candidateNucleon = consideredPartners[i];
262 newZ = oldZ + candidateNucleon.
Z;
266 if(newZ > clusterZMaxAll || newN > clusterNMaxAll)
272 const G4double phaseSpace = getPhaseSpace(oldA, candidateNucleon);
273 if(phaseSpace > phaseSpaceCut)
continue;
276 runningConfiguration[oldAMinusOne] = i;
277 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
278 Hashing::HashType configHash;
279 HashIterator aHashIter;
280 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
281 SortedNucleonConfiguration thisConfig;
282 SortedNucleonConfigurationIterator thisConfigIter;
285 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
286 configHash = Hashing::hashConfig(runningConfiguration, oldA);
287 aHashIter = theHashContainer->lower_bound(configHash);
289 if(aHashIter!=theHashContainer->end()
290 && !(configHash < *aHashIter))
292 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
293 thisConfig.fill(runningConfiguration,oldA);
294 thisConfigIter = theConfigContainer->lower_bound(thisConfig);
296 if(thisConfigIter!=theConfigContainer->end()
297 && !(thisConfig < *thisConfigIter))
303 runningEnergies[newA] = runningEnergies[oldA] + candidateNucleon.energy;
305 runningPotentials[newA] = runningPotentials[oldA] + candidateNucleon.potentialEnergy;
308 G4double oldCascadingEnergyPool = cascadingEnergyPool;
309 if(!candidateNucleon.isTargetSpectator)
310 cascadingEnergyPool -= candidateNucleon.energy - candidateNucleon.potentialEnergy - 931.3;
315 const G4double halfB = 0.72 * newZ *
317 const G4double tout = runningEnergies[newA] - runningPotentials[newA] -
319 if(tout<=halfB && tout+cascadingEnergyPool<=halfB) {
320 cascadingEnergyPool = oldCascadingEnergyPool;
325 runningPositions[newA] = (runningPositions[oldA] * oldA + candidateNucleon.position)*clusterPosFact[newA];
326 runningMomenta[newA] = runningMomenta[oldA] + candidateNucleon.momentum;
330 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
331 theHashContainer->insert(aHashIter, configHash);
332 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
333 theConfigContainer->insert(thisConfigIter, thisConfig);
339 isInRunningConfiguration[i] =
true;
342 if(newZ >= ZMinForNewA && newZ <= ZMaxForNewA) {
345 runningMomenta[newA]);
346 const G4double sqct = (sqc - 2.*newZ*protonMass - 2.*(newA-newZ)*neutronMass
348 *clusterPosFact[newA];
358 for(
G4int j=0; j<oldA; ++j)
359 candidateConfiguration[j] = consideredPartners[runningConfiguration[j]].particle;
367 if(newA < runningMaxClusterAlgorithmMass && newA+1 < theNucleus->
getA()) {
368 findClusterStartingFrom(newA, newZ);
372 isInRunningConfiguration[i] =
false;
373 cascadingEnergyPool = oldCascadingEnergyPool;
386 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)
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.
G4double energy(const ThreeVector &p, const G4double m)
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.
static const G4double pos
ParticleList::const_iterator ParticleIter