| 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