33 #define INCLXX_IN_GEANT4_MODE 1
44 const G4double ClusteringModelIntercomparison::limitCosEscapeAngle = 0.7;
46 static G4bool cascadingFirstPredicate(Particle *aParticle) {
47 return !aParticle->isTargetSpectator();
53 runningMaxClusterAlgorithmMass = std::min(maxClusterAlgorithmMass, nucleus->
getA()/2);
56 if(runningMaxClusterAlgorithmMass<=1)
60 Particle *theLeadingParticle = particle;
79 const G4double arg = rmaxws*rmaxws - Rprime*Rprime;
84 const G4double cosmin = std::sqrt(arg)/rmaxws;
87 translat = rmaxws * cospr;
90 translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin));
94 translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr));
98 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->
getMomentum() * (translat/pk);
100 theLeadingParticle->
setPosition(leadingParticlePosition);
103 const G4int theNucleusA = theNucleus->
getA();
104 if(nConsideredMax < theNucleusA) {
105 delete [] consideredPartners;
106 delete [] isInRunningConfiguration;
107 nConsideredMax = 2*theNucleusA;
108 consideredPartners =
new Particle *[nConsideredMax];
109 isInRunningConfiguration =
new G4bool [nConsideredMax];
110 std::fill(isInRunningConfiguration,
111 isInRunningConfiguration + nConsideredMax,
117 cascadingEnergyPool = 0.;
120 for(
ParticleIter i = particles.begin(); i != particles.end(); ++i) {
121 if (!(*i)->isNucleon())
continue;
122 if ((*i)->getID() == theLeadingParticle->
getID())
continue;
124 G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2();
125 G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2();
131 consideredPartners[nConsidered] = *i;
134 if(!(*i)->isTargetSpectator())
135 cascadingEnergyPool += (*i)->getEnergy() - (*i)->getPotentialEnergy() - 931.3;
145 std::partition(consideredPartners, consideredPartners+nConsidered, cascadingFirstPredicate);
150 maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2;
151 for(
G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i)
152 checkedConfigurations[i].
clear();
156 runningPositions[1] = leadingParticlePosition;
157 runningMomenta[1] = leadingParticleMomentum;
158 runningEnergies[1] = theLeadingParticle->
getEnergy();
165 findClusterStartingFrom(1, theLeadingParticle->
getZ());
173 candidateConfiguration[selectedA-1] = theLeadingParticle;
174 chosenCluster =
new Cluster(candidateConfiguration,
175 candidateConfiguration + selectedA);
179 theLeadingParticle->
setPosition(oldLeadingParticlePosition);
181 return chosenCluster;
184 inline G4double ClusteringModelIntercomparison::getPhaseSpace(
const G4int oldA,
Particle const *
const p) {
190 void ClusteringModelIntercomparison::findClusterStartingFrom(
const G4int oldA,
const G4int oldZ) {
191 const G4int newA = oldA + 1;
192 const G4int oldAMinusOne = oldA - 1;
200 const G4bool cachingEnabled = (newA<=maxMassConfigurationSkipping && newA>=3);
203 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
204 HashContainer *theHashContainer;
206 theHashContainer = &(checkedConfigurations[oldA-2]);
208 theHashContainer = NULL;
209 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
210 SortedNucleonConfigurationContainer *theConfigContainer;
212 theConfigContainer = &(checkedConfigurations[oldA-2]);
214 theConfigContainer = NULL;
216 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask.
223 for(
G4int i=0; i<nConsidered; ++i) {
225 if(isInRunningConfiguration[i])
continue;
227 Particle *
const candidateNucleon = consideredPartners[i];
230 newZ = oldZ + candidateNucleon->
getZ();
234 if(newZ > clusterZMaxAll || newN > clusterNMaxAll)
240 const G4double phaseSpace = getPhaseSpace(oldA, candidateNucleon);
241 if(phaseSpace > phaseSpaceCut)
continue;
244 runningConfiguration[oldAMinusOne] = i;
245 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
246 Hashing::HashType configHash;
247 HashIterator aHashIter;
248 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
249 SortedNucleonConfiguration thisConfig;
250 SortedNucleonConfigurationIterator thisConfigIter;
253 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
254 configHash = Hashing::hashConfig(runningConfiguration, oldA);
255 aHashIter = theHashContainer->lower_bound(configHash);
257 if(aHashIter!=theHashContainer->end()
258 && !(configHash < *aHashIter))
260 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
261 thisConfig.fill(runningConfiguration,oldA);
262 thisConfigIter = theConfigContainer->lower_bound(thisConfig);
264 if(thisConfigIter!=theConfigContainer->end()
265 && !(thisConfig < *thisConfigIter))
271 runningEnergies[newA] = runningEnergies[oldA] + candidateNucleon->getEnergy();
273 runningPotentials[newA] = runningPotentials[oldA] + candidateNucleon->getPotentialEnergy();
276 G4double oldCascadingEnergyPool = cascadingEnergyPool;
277 if(!candidateNucleon->isTargetSpectator())
278 cascadingEnergyPool -= candidateNucleon->getEnergy() - candidateNucleon->getPotentialEnergy() - 931.3;
283 const G4double halfB = 0.72 * newZ *
285 const G4double tout = runningEnergies[newA] - runningPotentials[newA] -
287 if(tout<=halfB && tout+cascadingEnergyPool<=halfB) {
288 cascadingEnergyPool = oldCascadingEnergyPool;
294 runningMomenta[newA] = runningMomenta[oldA] + candidateNucleon->getMomentum();
298 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
299 theHashContainer->insert(aHashIter, configHash);
300 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
301 theConfigContainer->insert(thisConfigIter, thisConfig);
306 isInRunningConfiguration[i] =
true;
309 if(newZ >= ZMinForNewA && newZ <= ZMaxForNewA) {
312 runningMomenta[newA]);
313 const G4double sqct = (sqc - 2.*newZ*protonMass - 2.*(newA-newZ)*neutronMass
325 for(
G4int j=0; j<oldA; ++j)
326 candidateConfiguration[j] = consideredPartners[runningConfiguration[j]];
334 if(newA < runningMaxClusterAlgorithmMass && newA+1 < theNucleus->
getA()) {
335 findClusterStartingFrom(newA, newZ);
339 isInRunningConfiguration[i] =
false;
340 cascadingEnergyPool = oldCascadingEnergyPool;
353 if(cosEscapeAngle < limitCosEscapeAngle)