Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland 31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1 35 35 36 #include "globals.hh" 36 #include "globals.hh" 37 37 38 #include "G4INCLClusteringModelIntercomparison 38 #include "G4INCLClusteringModelIntercomparison.hh" 39 #include "G4INCLCluster.hh" 39 #include "G4INCLCluster.hh" 40 #include "G4INCLRandom.hh" 40 #include "G4INCLRandom.hh" 41 #include "G4INCLHashing.hh" 41 #include "G4INCLHashing.hh" 42 #include <algorithm> 42 #include <algorithm> 43 43 44 namespace G4INCL { 44 namespace G4INCL { 45 45 46 const G4int ClusteringModelIntercomparison:: 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:: 47 const G4int ClusteringModelIntercomparison::clusterZMax[ParticleTable::maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8}; 48 48 49 const G4double ClusteringModelIntercompariso 49 const G4double ClusteringModelIntercomparison::clusterPosFact[ParticleTable::maxClusterMass+1] = {0.0, 1.0, 0.5, 50 0.33333, 0.25, 50 0.33333, 0.25, 51 0.2, 0.16667, 51 0.2, 0.16667, 52 0.14286, 0.125, 52 0.14286, 0.125, 53 0.11111, 0.1, 53 0.11111, 0.1, 54 0.09091, 0.083333}; 54 0.09091, 0.083333}; 55 55 56 const G4double ClusteringModelIntercompariso 56 const G4double ClusteringModelIntercomparison::clusterPosFact2[ParticleTable::maxClusterMass+1] = {0.0, 1.0, 0.25, 57 0.11111, 0.0625, 57 0.11111, 0.0625, 58 0.04, 0.0277778, 58 0.04, 0.0277778, 59 0.020408, 0.015625, 59 0.020408, 0.015625, 60 0.012346, 0.01, 60 0.012346, 0.01, 61 0.0082645, 0.0069444}; 61 0.0082645, 0.0069444}; 62 62 63 const G4double ClusteringModelIntercompariso 63 const G4double ClusteringModelIntercomparison::clusterPhaseSpaceCut[ParticleTable::maxClusterMass+1] = {0.0, 70000.0, 180000.0, 64 90000.0, 90000.0, 64 90000.0, 90000.0, 65 128941.0 ,145607.0, 65 128941.0 ,145607.0, 66 161365.0, 176389.0, 66 161365.0, 176389.0, 67 190798.0, 204681.0, 67 190798.0, 204681.0, 68 218109.0, 231135.0}; 68 218109.0, 231135.0}; 69 69 70 const G4double ClusteringModelIntercompariso 70 const G4double ClusteringModelIntercomparison::limitCosEscapeAngle = 0.7; 71 71 72 #ifdef INCL_DO_NOT_COMPILE 72 #ifdef INCL_DO_NOT_COMPILE 73 namespace { 73 namespace { 74 G4bool cascadingFirstPredicate(ConsideredP 74 G4bool cascadingFirstPredicate(ConsideredPartner const &aPartner) { 75 return !aPartner.isTargetSpectator; 75 return !aPartner.isTargetSpectator; 76 } 76 } 77 } 77 } 78 #endif 78 #endif 79 79 80 Cluster* ClusteringModelIntercomparison::get 80 Cluster* ClusteringModelIntercomparison::getCluster(Nucleus *nucleus, Particle *particle) { 81 // Set the maximum clustering mass dynamic 81 // Set the maximum clustering mass dynamically, based on the current nucleus 82 const G4int maxClusterAlgorithmMass = nucl 82 const G4int maxClusterAlgorithmMass = nucleus->getStore()->getConfig()->getClusterMaxMass(); 83 runningMaxClusterAlgorithmMass = std::min( 83 runningMaxClusterAlgorithmMass = std::min(maxClusterAlgorithmMass, nucleus->getA()/2); 84 84 85 // Nucleus too small? 85 // Nucleus too small? 86 if(runningMaxClusterAlgorithmMass<=1) 86 if(runningMaxClusterAlgorithmMass<=1) 87 return NULL; 87 return NULL; 88 88 89 theNucleus = nucleus; 89 theNucleus = nucleus; 90 Particle *theLeadingParticle = particle; 90 Particle *theLeadingParticle = particle; 91 91 92 // Initialise sqtot to a large number 92 // Initialise sqtot to a large number 93 sqtot = 50000.0; 93 sqtot = 50000.0; 94 selectedA = 0; 94 selectedA = 0; 95 selectedZ = 0; 95 selectedZ = 0; 96 96 97 // The distance parameter, known as h in p 97 // The distance parameter, known as h in publications. 98 // Default value is 1 fm. 98 // Default value is 1 fm. 99 const G4double transp = 1.0; 99 const G4double transp = 1.0; 100 100 101 const G4double rmaxws = theNucleus->getUni 101 const G4double rmaxws = theNucleus->getUniverseRadius(); 102 102 103 // Radius of the sphere where the leading 103 // Radius of the sphere where the leading particle is positioned. 104 const G4double Rprime = theNucleus->getDen 104 const G4double Rprime = theNucleus->getDensity()->getProtonNuclearRadius() + transp; 105 105 106 // Bring the leading particle back to the 106 // Bring the leading particle back to the coalescence sphere 107 const G4double pk = theLeadingParticle->ge 107 const G4double pk = theLeadingParticle->getMomentum().mag(); 108 const G4double cospr = theLeadingParticle- 108 const G4double cospr = theLeadingParticle->getPosition().dot(theLeadingParticle->getMomentum())/(theNucleus->getUniverseRadius() * pk); 109 const G4double arg = rmaxws*rmaxws - Rprim 109 const G4double arg = rmaxws*rmaxws - Rprime*Rprime; 110 G4double translat; 110 G4double translat; 111 111 112 if(arg > 0.0) { 112 if(arg > 0.0) { 113 // coalescence sphere smaller than Rmax 113 // coalescence sphere smaller than Rmax 114 const G4double cosmin = std::sqrt(arg)/r 114 const G4double cosmin = std::sqrt(arg)/rmaxws; 115 if(cospr <= cosmin) { 115 if(cospr <= cosmin) { 116 // there is an intersection with the c 116 // there is an intersection with the coalescence sphere 117 translat = rmaxws * cospr; 117 translat = rmaxws * cospr; 118 } else { 118 } else { 119 // no intersection with the coalescenc 119 // no intersection with the coalescence sphere 120 translat = rmaxws * (cospr - std::sqrt 120 translat = rmaxws * (cospr - std::sqrt(cospr*cospr - cosmin*cosmin)); 121 } 121 } 122 } else { 122 } else { 123 // coalescence sphere larger than Rmax 123 // coalescence sphere larger than Rmax 124 translat = rmaxws * cospr - std::sqrt(Rp 124 translat = rmaxws * cospr - std::sqrt(Rprime*Rprime - rmaxws*rmaxws*(1.0 - cospr*cospr)); 125 } 125 } 126 126 127 const ThreeVector oldLeadingParticlePositi 127 const ThreeVector oldLeadingParticlePosition = theLeadingParticle->getPosition(); 128 const ThreeVector leadingParticlePosition 128 const ThreeVector leadingParticlePosition = oldLeadingParticlePosition - theLeadingParticle->getMomentum() * (translat/pk); 129 const ThreeVector &leadingParticleMomentum 129 const ThreeVector &leadingParticleMomentum = theLeadingParticle->getMomentum(); 130 theLeadingParticle->setPosition(leadingPar 130 theLeadingParticle->setPosition(leadingParticlePosition); 131 131 132 // Initialise the array of considered nucl 132 // Initialise the array of considered nucleons 133 const G4int theNucleusA = theNucleus->getA 133 const G4int theNucleusA = theNucleus->getA(); 134 if(nConsideredMax < theNucleusA) { 134 if(nConsideredMax < theNucleusA) { 135 delete [] consideredPartners; 135 delete [] consideredPartners; 136 delete [] isInRunningConfiguration; 136 delete [] isInRunningConfiguration; 137 nConsideredMax = 2*theNucleusA; 137 nConsideredMax = 2*theNucleusA; 138 consideredPartners = new ConsideredPartn 138 consideredPartners = new ConsideredPartner[nConsideredMax]; 139 isInRunningConfiguration = new G4bool [n 139 isInRunningConfiguration = new G4bool [nConsideredMax]; 140 std::fill(isInRunningConfiguration, 140 std::fill(isInRunningConfiguration, 141 isInRunningConfiguration + nCo 141 isInRunningConfiguration + nConsideredMax, 142 false); 142 false); 143 } 143 } 144 144 145 // Select the subset of nucleons that will 145 // Select the subset of nucleons that will be considered in the 146 // cluster production: 146 // cluster production: 147 cascadingEnergyPool = 0.; 147 cascadingEnergyPool = 0.; 148 nConsidered = 0; 148 nConsidered = 0; 149 ParticleList const &particles = theNucleus 149 ParticleList const &particles = theNucleus->getStore()->getParticles(); 150 for(ParticleIter i=particles.begin(), e=pa 150 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) { 151 if (!(*i)->isNucleonorLambda()) continue << 151 if (!(*i)->isNucleon()) continue; // Only nucleons are allowed in clusters 152 if ((*i)->getID() == theLeadingParticle- 152 if ((*i)->getID() == theLeadingParticle->getID()) continue; // Don't count the leading particle 153 153 154 G4double space = ((*i)->getPosition() - 154 G4double space = ((*i)->getPosition() - leadingParticlePosition).mag2(); 155 G4double momentum = ((*i)->getMomentum() 155 G4double momentum = ((*i)->getMomentum() - leadingParticleMomentum).mag2(); 156 G4double size = space*momentum*clusterPo 156 G4double size = space*momentum*clusterPosFact2[runningMaxClusterAlgorithmMass]; 157 // Nucleons are accepted only if they ar 157 // Nucleons are accepted only if they are "close enough" in phase space 158 // to the leading nucleon. The selected 158 // to the leading nucleon. The selected phase-space parameter corresponds 159 // to the running maximum cluster mass. 159 // to the running maximum cluster mass. 160 if(size < clusterPhaseSpaceCut[runningMa 160 if(size < clusterPhaseSpaceCut[runningMaxClusterAlgorithmMass]) { 161 consideredPartners[nConsidered] = *i; << 161 consideredPartners[nConsidered] = *i; 162 // Keep trace of how much energy is ca 162 // Keep trace of how much energy is carried by cascading nucleons. This 163 // is used to stop the clustering algo 163 // is used to stop the clustering algorithm as soon as possible. 164 if(!(*i)->isTargetSpectator()) 164 if(!(*i)->isTargetSpectator()) 165 cascadingEnergyPool += consideredPar 165 cascadingEnergyPool += consideredPartners[nConsidered].energy - consideredPartners[nConsidered].potentialEnergy - 931.3; 166 nConsidered++; 166 nConsidered++; 167 // Make sure we don't exceed the array 167 // Make sure we don't exceed the array size 168 // assert(nConsidered<=nConsideredMax); 168 // assert(nConsidered<=nConsideredMax); 169 } 169 } 170 } 170 } 171 // Sort the list of considered partners so 171 // Sort the list of considered partners so that we give priority 172 // to participants. As soon as we encounte 172 // to participants. As soon as we encounter the first spectator in 173 // the list we know that all the remaining 173 // the list we know that all the remaining nucleons will be 174 // spectators too. 174 // spectators too. 175 // std::partition(consideredPartners, consi 175 // std::partition(consideredPartners, consideredPartners+nConsidered, cascadingFirstPredicate); 176 176 177 #ifndef INCL_CACHING_CLUSTERING_MODEL_INTERCOM 177 #ifndef INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None 178 // Clear the sets of checked configuration 178 // Clear the sets of checked configurations 179 // We stop caching two masses short of the 179 // We stop caching two masses short of the max mass -- there seems to be a 180 // performance hit above 180 // performance hit above 181 maxMassConfigurationSkipping = runningMaxC 181 maxMassConfigurationSkipping = runningMaxClusterAlgorithmMass-2; 182 for(G4int i=0; i<runningMaxClusterAlgorith 182 for(G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i) // no caching for A=1,2 183 checkedConfigurations[i].clear(); 183 checkedConfigurations[i].clear(); 184 #endif 184 #endif 185 185 186 // Initialise position, momentum and energ 186 // Initialise position, momentum and energy of the running cluster 187 // configuration 187 // configuration 188 runningPositions[1] = leadingParticlePosit 188 runningPositions[1] = leadingParticlePosition; 189 runningMomenta[1] = leadingParticleMomentu 189 runningMomenta[1] = leadingParticleMomentum; 190 runningEnergies[1] = theLeadingParticle->g 190 runningEnergies[1] = theLeadingParticle->getEnergy(); 191 runningPotentials[1] = theLeadingParticle- 191 runningPotentials[1] = theLeadingParticle->getPotentialEnergy(); 192 192 193 // Make sure that all the elements of isIn 193 // Make sure that all the elements of isInRunningConfiguration are false. 194 // assert(std::count(isInRunningConfiguration, 194 // assert(std::count(isInRunningConfiguration, isInRunningConfiguration+nConsidered, true)==0); 195 195 196 // Start the cluster search! 196 // Start the cluster search! 197 findClusterStartingFrom(1, theLeadingParti << 197 findClusterStartingFrom(1, theLeadingParticle->getZ()); 198 198 199 // Again, make sure that all the elements 199 // Again, make sure that all the elements of isInRunningConfiguration have 200 // been reset to false. This is a sanity c 200 // been reset to false. This is a sanity check. 201 // assert(std::count(isInRunningConfiguration, 201 // assert(std::count(isInRunningConfiguration, isInRunningConfiguration+nConsidered, true)==0); 202 202 203 Cluster *chosenCluster = 0; 203 Cluster *chosenCluster = 0; 204 if(selectedA!=0) { // A cluster was found! 204 if(selectedA!=0) { // A cluster was found! 205 candidateConfiguration[selectedA-1] = th 205 candidateConfiguration[selectedA-1] = theLeadingParticle; 206 chosenCluster = new Cluster(candidateCo 206 chosenCluster = new Cluster(candidateConfiguration, 207 candidateConfiguration + selectedA); 207 candidateConfiguration + selectedA); 208 } 208 } 209 209 210 // Restore the original position of the le 210 // Restore the original position of the leading particle 211 theLeadingParticle->setPosition(oldLeading 211 theLeadingParticle->setPosition(oldLeadingParticlePosition); 212 212 213 return chosenCluster; 213 return chosenCluster; 214 } 214 } 215 215 216 G4double ClusteringModelIntercomparison::get 216 G4double ClusteringModelIntercomparison::getPhaseSpace(const G4int oldA, ConsideredPartner const &p) { 217 const G4double psSpace = (p.position - run 217 const G4double psSpace = (p.position - runningPositions[oldA]).mag2(); 218 const G4double psMomentum = (p.momentum*ol 218 const G4double psMomentum = (p.momentum*oldA - runningMomenta[oldA]).mag2(); 219 return psSpace * psMomentum * clusterPosFa 219 return psSpace * psMomentum * clusterPosFact2[oldA + 1]; 220 } 220 } 221 221 222 void ClusteringModelIntercomparison::findClu << 222 void ClusteringModelIntercomparison::findClusterStartingFrom(const G4int oldA, const G4int oldZ) { 223 const G4int newA = oldA + 1; 223 const G4int newA = oldA + 1; 224 const G4int oldAMinusOne = oldA - 1; 224 const G4int oldAMinusOne = oldA - 1; 225 G4int newZ; 225 G4int newZ; 226 G4int newN; 226 G4int newN; 227 G4int newS; << 228 227 229 // Look up the phase-space cut 228 // Look up the phase-space cut 230 const G4double phaseSpaceCut = clusterPhas 229 const G4double phaseSpaceCut = clusterPhaseSpaceCut[newA]; 231 230 232 // Configuration caching enabled only for 231 // Configuration caching enabled only for a certain mass interval 233 const G4bool cachingEnabled = (newA<=maxMa 232 const G4bool cachingEnabled = (newA<=maxMassConfigurationSkipping && newA>=3); 234 233 235 // Set the pointer to the container of cac 234 // Set the pointer to the container of cached configurations 236 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 235 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask) 237 HashContainer *theHashContainer; 236 HashContainer *theHashContainer; 238 if(cachingEnabled) 237 if(cachingEnabled) 239 theHashContainer = &(checkedConfiguratio 238 theHashContainer = &(checkedConfigurations[oldA-2]); 240 else 239 else 241 theHashContainer = NULL; 240 theHashContainer = NULL; 242 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 241 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) 243 SortedNucleonConfigurationContainer *theCo 242 SortedNucleonConfigurationContainer *theConfigContainer; 244 if(cachingEnabled) 243 if(cachingEnabled) 245 theConfigContainer = &(checkedConfigurat 244 theConfigContainer = &(checkedConfigurations[oldA-2]); 246 else 245 else 247 theConfigContainer = NULL; 246 theConfigContainer = NULL; 248 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_I 247 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None) 249 #error Unrecognized INCL_CACHING_CLUSTERING_MO 248 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None. 250 #endif 249 #endif 251 250 252 // Minimum and maximum Z values for this m 251 // Minimum and maximum Z values for this mass 253 const G4int ZMinForNewA = clusterZMin[newA 252 const G4int ZMinForNewA = clusterZMin[newA]; 254 const G4int ZMaxForNewA = clusterZMax[newA 253 const G4int ZMaxForNewA = clusterZMax[newA]; 255 254 256 for(G4int i=0; i<nConsidered; ++i) { 255 for(G4int i=0; i<nConsidered; ++i) { 257 // Only accept particles that are not al 256 // Only accept particles that are not already part of the cluster 258 if(isInRunningConfiguration[i]) continue 257 if(isInRunningConfiguration[i]) continue; 259 258 260 ConsideredPartner const &candidateNucleo 259 ConsideredPartner const &candidateNucleon = consideredPartners[i]; 261 260 262 // Z and A of the new cluster 261 // Z and A of the new cluster 263 newZ = oldZ + candidateNucleon.Z; 262 newZ = oldZ + candidateNucleon.Z; 264 newS = oldS + candidateNucleon.S; << 265 newN = newA - newZ; 263 newN = newA - newZ; 266 264 267 // Skip this nucleon if we already have 265 // Skip this nucleon if we already have too many protons or neutrons 268 if(newZ > clusterZMaxAll || newN > clust << 266 if(newZ > clusterZMaxAll || newN > clusterNMaxAll) 269 continue; 267 continue; 270 268 271 // Compute the phase space factor for a 269 // Compute the phase space factor for a new cluster which 272 // consists of the previous running clus 270 // consists of the previous running cluster and the new 273 // candidate nucleon: 271 // candidate nucleon: 274 const G4double phaseSpace = getPhaseSpac 272 const G4double phaseSpace = getPhaseSpace(oldA, candidateNucleon); 275 if(phaseSpace > phaseSpaceCut) continue; 273 if(phaseSpace > phaseSpaceCut) continue; 276 274 277 // Store the candidate nucleon in the ru 275 // Store the candidate nucleon in the running configuration 278 runningConfiguration[oldAMinusOne] = i; 276 runningConfiguration[oldAMinusOne] = i; 279 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 277 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask) 280 Hashing::HashType configHash; 278 Hashing::HashType configHash; 281 HashIterator aHashIter; 279 HashIterator aHashIter; 282 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 280 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) 283 SortedNucleonConfiguration thisConfig; 281 SortedNucleonConfiguration thisConfig; 284 SortedNucleonConfigurationIterator thisC 282 SortedNucleonConfigurationIterator thisConfigIter; 285 #endif 283 #endif 286 if(cachingEnabled) { 284 if(cachingEnabled) { 287 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 285 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask) 288 configHash = Hashing::hashConfig(runni 286 configHash = Hashing::hashConfig(runningConfiguration, oldA); 289 aHashIter = theHashContainer->lower_bo 287 aHashIter = theHashContainer->lower_bound(configHash); 290 // If we have already checked this con 288 // If we have already checked this configuration, skip it 291 if(aHashIter!=theHashContainer->end() 289 if(aHashIter!=theHashContainer->end() 292 && !(configHash < *aHashIter)) 290 && !(configHash < *aHashIter)) 293 continue; 291 continue; 294 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 292 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) 295 thisConfig.fill(runningConfiguration,o 293 thisConfig.fill(runningConfiguration,oldA); 296 thisConfigIter = theConfigContainer->l 294 thisConfigIter = theConfigContainer->lower_bound(thisConfig); 297 // If we have already checked this con 295 // If we have already checked this configuration, skip it 298 if(thisConfigIter!=theConfigContainer- 296 if(thisConfigIter!=theConfigContainer->end() 299 && !(thisConfig < *thisConfigIter)) 297 && !(thisConfig < *thisConfigIter)) 300 continue; 298 continue; 301 #endif 299 #endif 302 } 300 } 303 301 304 // Sum of the total energies of the clus 302 // Sum of the total energies of the cluster components 305 runningEnergies[newA] = runningEnergies[ 303 runningEnergies[newA] = runningEnergies[oldA] + candidateNucleon.energy; 306 // Sum of the potential energies of the 304 // Sum of the potential energies of the cluster components 307 runningPotentials[newA] = runningPotenti 305 runningPotentials[newA] = runningPotentials[oldA] + candidateNucleon.potentialEnergy; 308 306 309 // Update the available cascading kineti 307 // Update the available cascading kinetic energy 310 G4double oldCascadingEnergyPool = cascad 308 G4double oldCascadingEnergyPool = cascadingEnergyPool; 311 if(!candidateNucleon.isTargetSpectator) 309 if(!candidateNucleon.isTargetSpectator) 312 cascadingEnergyPool -= candidateNucleo 310 cascadingEnergyPool -= candidateNucleon.energy - candidateNucleon.potentialEnergy - 931.3; 313 311 314 // Check an approximate Coulomb barrier. 312 // Check an approximate Coulomb barrier. If the cluster is below 315 // 0.5*barrier and the remaining availab 313 // 0.5*barrier and the remaining available energy from cascading nucleons 316 // will not bring it above, reject the c 314 // will not bring it above, reject the cluster. 317 const G4double halfB = 0.72 * newZ * 315 const G4double halfB = 0.72 * newZ * 318 theNucleus->getZ()/(theNucleus->getDen 316 theNucleus->getZ()/(theNucleus->getDensity()->getProtonNuclearRadius()+1.7); 319 const G4double tout = runningEnergies[ne 317 const G4double tout = runningEnergies[newA] - runningPotentials[newA] - 320 931.3*newA; 318 931.3*newA; 321 if(tout<=halfB && tout+cascadingEnergyPo 319 if(tout<=halfB && tout+cascadingEnergyPool<=halfB) { 322 cascadingEnergyPool = oldCascadingEner 320 cascadingEnergyPool = oldCascadingEnergyPool; 323 continue; 321 continue; 324 } 322 } 325 323 326 // Here the nucleon has passed all the t 324 // Here the nucleon has passed all the tests. Accept it in the cluster. 327 runningPositions[newA] = (runningPositio 325 runningPositions[newA] = (runningPositions[oldA] * oldA + candidateNucleon.position)*clusterPosFact[newA]; 328 runningMomenta[newA] = runningMomenta[ol 326 runningMomenta[newA] = runningMomenta[oldA] + candidateNucleon.momentum; 329 327 330 // Add the config to the container 328 // Add the config to the container 331 if(cachingEnabled) { 329 if(cachingEnabled) { 332 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 330 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask) 333 theHashContainer->insert(aHashIter, co 331 theHashContainer->insert(aHashIter, configHash); 334 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 332 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) 335 theConfigContainer->insert(thisConfigI 333 theConfigContainer->insert(thisConfigIter, thisConfig); 336 #endif 334 #endif 337 } 335 } 338 336 339 // Set the flag that reminds us that thi 337 // Set the flag that reminds us that this nucleon has already been taken 340 // in the running configuration 338 // in the running configuration 341 isInRunningConfiguration[i] = true; 339 isInRunningConfiguration[i] = true; 342 340 343 // Keep track of the best physical clust 341 // Keep track of the best physical cluster 344 if(newZ >= ZMinForNewA && newZ <= ZMaxFo 342 if(newZ >= ZMinForNewA && newZ <= ZMaxForNewA) { 345 // Note: sqc is real kinetic energy, n 343 // Note: sqc is real kinetic energy, not the square of the kinetic energy! 346 const G4double sqc = KinematicsUtils:: 344 const G4double sqc = KinematicsUtils::invariantMass(runningEnergies[newA], 347 runningMomenta[newA]); 345 runningMomenta[newA]); 348 const G4double sqct = (sqc - 2.*newZ*p << 346 const G4double sqct = (sqc - 2.*newZ*protonMass - 2.*(newA-newZ)*neutronMass 349 + ParticleTable::getRealMass(newA << 347 + ParticleTable::getRealMass(newA, newZ)) 350 *clusterPosFact[newA]; 348 *clusterPosFact[newA]; 351 349 352 if(sqct < sqtot) { 350 if(sqct < sqtot) { 353 // This is the best cluster we have 351 // This is the best cluster we have found so far. Store its 354 // kinematics. 352 // kinematics. 355 sqtot = sqct; 353 sqtot = sqct; 356 selectedA = newA; 354 selectedA = newA; 357 selectedZ = newZ; 355 selectedZ = newZ; 358 selectedS = newS; << 359 356 360 // Store the running configuration i 357 // Store the running configuration in a ParticleList 361 for(G4int j=0; j<oldA; ++j) 358 for(G4int j=0; j<oldA; ++j) 362 candidateConfiguration[j] = consid 359 candidateConfiguration[j] = consideredPartners[runningConfiguration[j]].particle; 363 360 364 // Sanity check on number of nucleon 361 // Sanity check on number of nucleons in running configuration 365 // assert(std::count(isInRunningConfiguration, 362 // assert(std::count(isInRunningConfiguration, isInRunningConfiguration+nConsidered, true)==selectedA-1); 366 } 363 } 367 } 364 } 368 365 369 // The method recursively calls itself f 366 // The method recursively calls itself for the next mass 370 if(newA < runningMaxClusterAlgorithmMass << 367 if(newA < runningMaxClusterAlgorithmMass && newA+1 < theNucleus->getA()) { 371 findClusterStartingFrom(newA, newZ, newS); << 368 findClusterStartingFrom(newA, newZ); 372 } 369 } 373 370 374 // Reset the running configuration flag 371 // Reset the running configuration flag and the cascading energy pool 375 isInRunningConfiguration[i] = false; 372 isInRunningConfiguration[i] = false; 376 cascadingEnergyPool = oldCascadingEnergy 373 cascadingEnergyPool = oldCascadingEnergyPool; 377 } 374 } 378 } 375 } 379 376 380 G4bool ClusteringModelIntercomparison::clust 377 G4bool ClusteringModelIntercomparison::clusterCanEscape(Nucleus const * const n, Cluster const * const c) { 381 // Forbid emission of the whole nucleus 378 // Forbid emission of the whole nucleus 382 if(c->getA()>=n->getA() || c->getS()>0) << 379 if(c->getA()>=n->getA()) 383 return false; 380 return false; 384 381 385 // Check the escape angle of the cluster 382 // Check the escape angle of the cluster 386 const ThreeVector &pos = c->getPosition(); 383 const ThreeVector &pos = c->getPosition(); 387 const ThreeVector &mom = c->getMomentum(); 384 const ThreeVector &mom = c->getMomentum(); 388 const G4double cosEscapeAngle = pos.dot(mo 385 const G4double cosEscapeAngle = pos.dot(mom) / std::sqrt(pos.mag2()*mom.mag2()); 389 if(cosEscapeAngle < limitCosEscapeAngle) 386 if(cosEscapeAngle < limitCosEscapeAngle) 390 return false; 387 return false; 391 388 392 return true; 389 return true; 393 } 390 } 394 391 395 } 392 } 396 393