00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037 #include "G4INCLClusteringModelIntercomparison.hh"
00038 #include "G4INCLCluster.hh"
00039 #include "G4INCLRandom.hh"
00040 #include "G4INCLHashing.hh"
00041 #include <algorithm>
00042
00043 namespace G4INCL {
00044 const G4double ClusteringModelIntercomparison::limitCosEscapeAngle = 0.7;
00045
00046 static G4bool cascadingFirstPredicate(Particle *aParticle) {
00047 return !aParticle->isTargetSpectator();
00048 }
00049
00050 Cluster* ClusteringModelIntercomparison::getCluster(Nucleus *nucleus, Particle *particle) {
00051
00052 const G4int maxClusterAlgorithmMass = nucleus->getStore()->getConfig()->getClusterMaxMass();
00053 runningMaxClusterAlgorithmMass = std::min(maxClusterAlgorithmMass, nucleus->getA()/2);
00054
00055
00056 if(runningMaxClusterAlgorithmMass<=1)
00057 return NULL;
00058
00059 theNucleus = nucleus;
00060 Particle *theLeadingParticle = particle;
00061
00062
00063 sqtot = 50000.0;
00064 selectedA = 0;
00065 selectedZ = 0;
00066
00067
00068
00069 const G4double transp = 1.0;
00070
00071 const G4double rmaxws = theNucleus->getUniverseRadius();
00072
00073
00074 const G4double Rprime = theNucleus->getDensity()->getNuclearRadius() + transp;
00075
00076
00077 const G4double pk = theLeadingParticle->getMomentum().mag();
00078 const G4double cospr = theLeadingParticle->getPosition().dot(theLeadingParticle->getMomentum())/(theNucleus->getUniverseRadius() * pk);
00079 const G4double arg = rmaxws*rmaxws - Rprime*Rprime;
00080 G4double translat;
00081
00082 if(arg > 0.0) {
00083
00084 const G4double cosmin = std::sqrt(arg)/rmaxws;
00085 if(cospr <= cosmin) {
00086
00087 translat = rmaxws * cospr;
00088 } else {
00089
00090 translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin));
00091 }
00092 } else {
00093
00094 translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr));
00095 }
00096
00097 const ThreeVector oldLeadingParticlePosition = theLeadingParticle->getPosition();
00098 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->getMomentum() * (translat/pk);
00099 const ThreeVector &leadingParticleMomentum = theLeadingParticle->getMomentum();
00100 theLeadingParticle->setPosition(leadingParticlePosition);
00101
00102
00103 const G4int theNucleusA = theNucleus->getA();
00104 if(nConsideredMax < theNucleusA) {
00105 delete [] consideredPartners;
00106 delete [] isInRunningConfiguration;
00107 nConsideredMax = 2*theNucleusA;
00108 consideredPartners = new Particle *[nConsideredMax];
00109 isInRunningConfiguration = new G4bool [nConsideredMax];
00110 std::fill(isInRunningConfiguration,
00111 isInRunningConfiguration + nConsideredMax,
00112 false);
00113 }
00114
00115
00116
00117 cascadingEnergyPool = 0.;
00118 nConsidered = 0;
00119 const ParticleList particles = theNucleus->getStore()->getParticles();
00120 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
00121 if (!(*i)->isNucleon()) continue;
00122 if ((*i)->getID() == theLeadingParticle->getID()) continue;
00123
00124 G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2();
00125 G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2();
00126 G4double size = space*momentum*ParticleTable::clusterPosFact2[runningMaxClusterAlgorithmMass];
00127
00128
00129
00130 if(size < ParticleTable::clusterPhaseSpaceCut[runningMaxClusterAlgorithmMass]) {
00131 consideredPartners[nConsidered] = *i;
00132
00133
00134 if(!(*i)->isTargetSpectator())
00135 cascadingEnergyPool += (*i)->getEnergy() - (*i)->getPotentialEnergy() - 931.3;
00136 nConsidered++;
00137
00138
00139 }
00140 }
00141
00142
00143
00144
00145 std::partition(consideredPartners, consideredPartners+nConsidered, cascadingFirstPredicate);
00146
00147
00148
00149
00150 maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2;
00151 for(G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i)
00152 checkedConfigurations[i].clear();
00153
00154
00155
00156 runningPositions[1] = leadingParticlePosition;
00157 runningMomenta[1] = leadingParticleMomentum;
00158 runningEnergies[1] = theLeadingParticle->getEnergy();
00159 runningPotentials[1] = theLeadingParticle->getPotentialEnergy();
00160
00161
00162
00163
00164
00165 findClusterStartingFrom(1, theLeadingParticle->getZ());
00166
00167
00168
00169
00170
00171 Cluster *chosenCluster = 0;
00172 if(selectedA!=0) {
00173 candidateConfiguration[selectedA-1] = theLeadingParticle;
00174 chosenCluster = new Cluster(candidateConfiguration,
00175 candidateConfiguration + selectedA);
00176 }
00177
00178
00179 theLeadingParticle->setPosition(oldLeadingParticlePosition);
00180
00181 return chosenCluster;
00182 }
00183
00184 inline G4double ClusteringModelIntercomparison::getPhaseSpace(const G4int oldA, Particle const * const p) {
00185 const G4double psSpace = (p->getPosition() - runningPositions[oldA]).mag2();
00186 const G4double psMomentum = (p->getMomentum()*oldA - runningMomenta[oldA]).mag2();
00187 return psSpace * psMomentum * ParticleTable::clusterPosFact2[oldA + 1];
00188 }
00189
00190 void ClusteringModelIntercomparison::findClusterStartingFrom(const G4int oldA, const G4int oldZ) {
00191 const G4int newA = oldA + 1;
00192 const G4int oldAMinusOne = oldA - 1;
00193 G4int newZ;
00194 G4int newN;
00195
00196
00197 const G4double phaseSpaceCut = ParticleTable::clusterPhaseSpaceCut[newA];
00198
00199
00200 const G4bool cachingEnabled = (newA<=maxMassConfigurationSkipping && newA>=3);
00201
00202
00203 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
00204 HashContainer *theHashContainer;
00205 if(cachingEnabled)
00206 theHashContainer = &(checkedConfigurations[oldA-2]);
00207 else
00208 theHashContainer = NULL;
00209 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
00210 SortedNucleonConfigurationContainer *theConfigContainer;
00211 if(cachingEnabled)
00212 theConfigContainer = &(checkedConfigurations[oldA-2]);
00213 else
00214 theConfigContainer = NULL;
00215 #else
00216 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask.
00217 #endif
00218
00219
00220 const G4int ZMinForNewA = ParticleTable::clusterZMin[newA];
00221 const G4int ZMaxForNewA = ParticleTable::clusterZMax[newA];
00222
00223 for(G4int i=0; i<nConsidered; ++i) {
00224
00225 if(isInRunningConfiguration[i]) continue;
00226
00227 Particle * const candidateNucleon = consideredPartners[i];
00228
00229
00230 newZ = oldZ + candidateNucleon->getZ();
00231 newN = newA - newZ;
00232
00233
00234 if(newZ > clusterZMaxAll || newN > clusterNMaxAll)
00235 continue;
00236
00237
00238
00239
00240 const G4double phaseSpace = getPhaseSpace(oldA, candidateNucleon);
00241 if(phaseSpace > phaseSpaceCut) continue;
00242
00243
00244 runningConfiguration[oldAMinusOne] = i;
00245 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
00246 Hashing::HashType configHash;
00247 HashIterator aHashIter;
00248 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
00249 SortedNucleonConfiguration thisConfig;
00250 SortedNucleonConfigurationIterator thisConfigIter;
00251 #endif
00252 if(cachingEnabled) {
00253 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
00254 configHash = Hashing::hashConfig(runningConfiguration, oldA);
00255 aHashIter = theHashContainer->lower_bound(configHash);
00256
00257 if(aHashIter!=theHashContainer->end()
00258 && !(configHash < *aHashIter))
00259 continue;
00260 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
00261 thisConfig.fill(runningConfiguration,oldA);
00262 thisConfigIter = theConfigContainer->lower_bound(thisConfig);
00263
00264 if(thisConfigIter!=theConfigContainer->end()
00265 && !(thisConfig < *thisConfigIter))
00266 continue;
00267 #endif
00268 }
00269
00270
00271 runningEnergies[newA] = runningEnergies[oldA] + candidateNucleon->getEnergy();
00272
00273 runningPotentials[newA] = runningPotentials[oldA] + candidateNucleon->getPotentialEnergy();
00274
00275
00276 G4double oldCascadingEnergyPool = cascadingEnergyPool;
00277 if(!candidateNucleon->isTargetSpectator())
00278 cascadingEnergyPool -= candidateNucleon->getEnergy() - candidateNucleon->getPotentialEnergy() - 931.3;
00279
00280
00281
00282
00283 const G4double halfB = 0.72 * newZ *
00284 theNucleus->getZ()/(theNucleus->getDensity()->getNuclearRadius()+1.7);
00285 const G4double tout = runningEnergies[newA] - runningPotentials[newA] -
00286 931.3*newA;
00287 if(tout<=halfB && tout+cascadingEnergyPool<=halfB) {
00288 cascadingEnergyPool = oldCascadingEnergyPool;
00289 continue;
00290 }
00291
00292
00293 runningPositions[newA] = (runningPositions[oldA] * oldA + candidateNucleon->getPosition())*ParticleTable::clusterPosFact[newA];
00294 runningMomenta[newA] = runningMomenta[oldA] + candidateNucleon->getMomentum();
00295
00296
00297 if(cachingEnabled)
00298 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
00299 theHashContainer->insert(aHashIter, configHash);
00300 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
00301 theConfigContainer->insert(thisConfigIter, thisConfig);
00302 #endif
00303
00304
00305
00306 isInRunningConfiguration[i] = true;
00307
00308
00309 if(newZ >= ZMinForNewA && newZ <= ZMaxForNewA) {
00310
00311 const G4double sqc = KinematicsUtils::invariantMass(runningEnergies[newA],
00312 runningMomenta[newA]);
00313 const G4double sqct = (sqc - 2.*newZ*protonMass - 2.*(newA-newZ)*neutronMass
00314 + ParticleTable::getRealMass(newA, newZ))
00315 *ParticleTable::clusterPosFact[newA];
00316
00317 if(sqct < sqtot) {
00318
00319
00320 sqtot = sqct;
00321 selectedA = newA;
00322 selectedZ = newZ;
00323
00324
00325 for(G4int j=0; j<oldA; ++j)
00326 candidateConfiguration[j] = consideredPartners[runningConfiguration[j]];
00327
00328
00329
00330 }
00331 }
00332
00333
00334 if(newA < runningMaxClusterAlgorithmMass && newA+1 < theNucleus->getA()) {
00335 findClusterStartingFrom(newA, newZ);
00336 }
00337
00338
00339 isInRunningConfiguration[i] = false;
00340 cascadingEnergyPool = oldCascadingEnergyPool;
00341 }
00342 }
00343
00344 G4bool ClusteringModelIntercomparison::clusterCanEscape(Nucleus const * const n, Cluster const * const c) {
00345
00346 if(c->getA()>=n->getA())
00347 return false;
00348
00349
00350 const ThreeVector &pos = c->getPosition();
00351 const ThreeVector &mom = c->getMomentum();
00352 const G4double cosEscapeAngle = pos.dot(mom) / std::sqrt(pos.mag2()*mom.mag2());
00353 if(cosEscapeAngle < limitCosEscapeAngle)
00354 return false;
00355
00356 return true;
00357 }
00358
00359 }