| 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 #ifndef G4INCLClusteringModelIntercomparison_h 38 #ifndef G4INCLClusteringModelIntercomparison_hh
39 #define G4INCLClusteringModelIntercomparison_h 39 #define G4INCLClusteringModelIntercomparison_hh 1
40 40
41 #ifdef INCLXX_IN_GEANT4_MODE 41 #ifdef INCLXX_IN_GEANT4_MODE
42 #define INCL_CACHING_CLUSTERING_MODEL_INTERCOM 42 #define INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set 1
43 #endif // INCLXX_IN_GEANT4_MODE 43 #endif // INCLXX_IN_GEANT4_MODE
44 44
45 #include "G4INCLIClusteringModel.hh" 45 #include "G4INCLIClusteringModel.hh"
46 #include "G4INCLParticle.hh" 46 #include "G4INCLParticle.hh"
47 #include "G4INCLParticleTable.hh" 47 #include "G4INCLParticleTable.hh"
48 #include "G4INCLCluster.hh" 48 #include "G4INCLCluster.hh"
49 #include "G4INCLNucleus.hh" 49 #include "G4INCLNucleus.hh"
50 #include "G4INCLKinematicsUtils.hh" 50 #include "G4INCLKinematicsUtils.hh"
51 #include "G4INCLHashing.hh" 51 #include "G4INCLHashing.hh"
52 52
53 #include <set> 53 #include <set>
54 #include <algorithm> 54 #include <algorithm>
55 55
56 namespace G4INCL { 56 namespace G4INCL {
57 57
58 /** \brief Container for the relevant inform 58 /** \brief Container for the relevant information
59 * 59 *
60 * This struct contains all the information 60 * This struct contains all the information that is relevant for the
61 * clustering algorithm. It is probably more 61 * clustering algorithm. It is probably more compact than the Particles it
62 * feeds on, hopefully improving cache perfo 62 * feeds on, hopefully improving cache performance.
63 */ 63 */
64 struct ConsideredPartner { 64 struct ConsideredPartner {
65 Particle *particle; 65 Particle *particle;
66 G4bool isTargetSpectator; 66 G4bool isTargetSpectator;
67 G4int Z; 67 G4int Z;
68 G4int S; 68 G4int S;
69 ThreeVector position; 69 ThreeVector position;
70 ThreeVector momentum; 70 ThreeVector momentum;
71 G4double energy; 71 G4double energy;
72 G4double potentialEnergy; 72 G4double potentialEnergy;
73 73
74 ConsideredPartner() : 74 ConsideredPartner() :
75 particle(NULL), 75 particle(NULL),
76 isTargetSpectator(false), 76 isTargetSpectator(false),
77 Z(0), 77 Z(0),
78 S(0), 78 S(0),
79 energy(0.), 79 energy(0.),
80 potentialEnergy(0.) 80 potentialEnergy(0.)
81 {} 81 {}
82 82
83 ConsideredPartner(Particle * const p) : 83 ConsideredPartner(Particle * const p) :
84 particle(p), 84 particle(p),
85 isTargetSpectator(particle->isTargetSpec 85 isTargetSpectator(particle->isTargetSpectator()),
86 Z(particle->getZ()), 86 Z(particle->getZ()),
87 S(particle->getS()), 87 S(particle->getS()),
88 position(particle->getPosition()), 88 position(particle->getPosition()),
89 momentum(particle->getMomentum()), 89 momentum(particle->getMomentum()),
90 energy(particle->getEnergy()), 90 energy(particle->getEnergy()),
91 potentialEnergy(particle->getPotentialEn 91 potentialEnergy(particle->getPotentialEnergy())
92 {} 92 {}
93 }; 93 };
94 94
95 /// \brief Cluster coalescence algorithm use 95 /// \brief Cluster coalescence algorithm used in the IAEA intercomparison
96 class ClusteringModelIntercomparison : publi 96 class ClusteringModelIntercomparison : public IClusteringModel {
97 public: 97 public:
98 ClusteringModelIntercomparison(Config cons 98 ClusteringModelIntercomparison(Config const * const theConfig) :
99 theNucleus(NULL), 99 theNucleus(NULL),
100 selectedA(0), 100 selectedA(0),
101 selectedZ(0), 101 selectedZ(0),
102 selectedS(0), 102 selectedS(0),
103 sqtot(0.), 103 sqtot(0.),
104 cascadingEnergyPool(0.), 104 cascadingEnergyPool(0.),
105 protonMass(ParticleTable::getRealMass(Pr 105 protonMass(ParticleTable::getRealMass(Proton)),
106 neutronMass(ParticleTable::getRealMass(N 106 neutronMass(ParticleTable::getRealMass(Neutron)),
107 lambdaMass(ParticleTable::getRealMass(La 107 lambdaMass(ParticleTable::getRealMass(Lambda)),
108 runningMaxClusterAlgorithmMass(theConfig 108 runningMaxClusterAlgorithmMass(theConfig->getClusterMaxMass()),
109 nConsideredMax(0), 109 nConsideredMax(0),
110 nConsidered(0), 110 nConsidered(0),
111 consideredPartners(NULL), 111 consideredPartners(NULL),
112 isInRunningConfiguration(NULL), 112 isInRunningConfiguration(NULL),
113 maxMassConfigurationSkipping(ParticleTab 113 maxMassConfigurationSkipping(ParticleTable::maxClusterMass)
114 { 114 {
115 // Set up the maximum charge and neutron 115 // Set up the maximum charge and neutron number for clusters
116 clusterZMaxAll = 0; 116 clusterZMaxAll = 0;
117 clusterNMaxAll = 0; 117 clusterNMaxAll = 0;
118 for(G4int A=0; A<=runningMaxClusterAlgor 118 for(G4int A=0; A<=runningMaxClusterAlgorithmMass; ++A) {
119 if(clusterZMax[A]>clusterZMaxAll) 119 if(clusterZMax[A]>clusterZMaxAll)
120 clusterZMaxAll = clusterZMax[A]; 120 clusterZMaxAll = clusterZMax[A];
121 if(A-clusterZMin[A]>clusterNMaxAll) 121 if(A-clusterZMin[A]>clusterNMaxAll)
122 clusterNMaxAll = A-clusterZMin[A]; 122 clusterNMaxAll = A-clusterZMin[A];
123 } 123 }
124 std::fill(candidateConfiguration, 124 std::fill(candidateConfiguration,
125 candidateConfiguration + Parti 125 candidateConfiguration + ParticleTable::maxClusterMass,
126 static_cast<Particle*>(NULL)); 126 static_cast<Particle*>(NULL));
127 127
128 std::fill(runningEnergies, 128 std::fill(runningEnergies,
129 runningEnergies + ParticleTabl 129 runningEnergies + ParticleTable::maxClusterMass,
130 0.0); 130 0.0);
131 131
132 std::fill(runningPotentials, 132 std::fill(runningPotentials,
133 runningPotentials + ParticleTa 133 runningPotentials + ParticleTable::maxClusterMass,
134 0.0); 134 0.0);
135 135
136 std::fill(runningConfiguration, 136 std::fill(runningConfiguration,
137 runningConfiguration + Particl 137 runningConfiguration + ParticleTable::maxClusterMass,
138 -1); 138 -1);
139 139
140 } 140 }
141 141
142 virtual ~ClusteringModelIntercomparison() 142 virtual ~ClusteringModelIntercomparison() {
143 delete [] consideredPartners; 143 delete [] consideredPartners;
144 delete [] isInRunningConfiguration; 144 delete [] isInRunningConfiguration;
145 } 145 }
146 146
147 virtual Cluster* getCluster(Nucleus*, Part 147 virtual Cluster* getCluster(Nucleus*, Particle*);
148 virtual G4bool clusterCanEscape(Nucleus co 148 virtual G4bool clusterCanEscape(Nucleus const * const, Cluster const * const);
149 149
150 private: 150 private:
151 void findClusterStartingFrom(const G4int o 151 void findClusterStartingFrom(const G4int oldA, const G4int oldZ, const G4int oldS);
152 G4double getPhaseSpace(const G4int oldA, C 152 G4double getPhaseSpace(const G4int oldA, ConsideredPartner const &p);
153 153
154 Nucleus *theNucleus; 154 Nucleus *theNucleus;
155 155
156 G4double runningEnergies[ParticleTable::ma 156 G4double runningEnergies[ParticleTable::maxClusterMass+1];
157 ThreeVector runningMomenta[ParticleTable:: 157 ThreeVector runningMomenta[ParticleTable::maxClusterMass+1];
158 ThreeVector runningPositions[ParticleTable 158 ThreeVector runningPositions[ParticleTable::maxClusterMass+1];
159 G4double runningPotentials[ParticleTable:: 159 G4double runningPotentials[ParticleTable::maxClusterMass+1];
160 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 160 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
161 Hashing::NucleonItem runningConfiguration[ 161 Hashing::NucleonItem runningConfiguration[ParticleTable::maxClusterMass];
162 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 162 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) || defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None)
163 G4int runningConfiguration[ParticleTable:: 163 G4int runningConfiguration[ParticleTable::maxClusterMass];
164 #else 164 #else
165 #error Unrecognized INCL_CACHING_CLUSTERING_MO 165 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None.
166 #endif 166 #endif
167 167
168 G4int selectedA, selectedZ, selectedS; 168 G4int selectedA, selectedZ, selectedS;
169 G4double sqtot; 169 G4double sqtot;
170 170
171 G4int clusterZMaxAll, clusterNMaxAll; 171 G4int clusterZMaxAll, clusterNMaxAll;
172 172
173 G4double cascadingEnergyPool; 173 G4double cascadingEnergyPool;
174 174
175 /// \brief Lower limit of Z for cluster of 175 /// \brief Lower limit of Z for cluster of mass A
176 static const G4int clusterZMin[ParticleTab 176 static const G4int clusterZMin[ParticleTable::maxClusterMass+1];
177 /// \brief Upper limit of Z for cluster of 177 /// \brief Upper limit of Z for cluster of mass A
178 static const G4int clusterZMax[ParticleTab 178 static const G4int clusterZMax[ParticleTable::maxClusterMass+1];
179 179
180 /// \brief Precomputed factor 1.0/A 180 /// \brief Precomputed factor 1.0/A
181 static const G4double clusterPosFact[Parti 181 static const G4double clusterPosFact[ParticleTable::maxClusterMass+1];
182 182
183 /// \brief Precomputed factor (1.0/A)^2 183 /// \brief Precomputed factor (1.0/A)^2
184 static const G4double clusterPosFact2[Part 184 static const G4double clusterPosFact2[ParticleTable::maxClusterMass+1];
185 185
186 /// \brief Phase-space parameters for clus 186 /// \brief Phase-space parameters for cluster formation
187 static const G4double clusterPhaseSpaceCut 187 static const G4double clusterPhaseSpaceCut[ParticleTable::maxClusterMass+1];
188 188
189 static const G4double limitCosEscapeAngle; 189 static const G4double limitCosEscapeAngle;
190 190
191 const G4double protonMass; 191 const G4double protonMass;
192 const G4double neutronMass; 192 const G4double neutronMass;
193 const G4double lambdaMass; 193 const G4double lambdaMass;
194 194
195 G4int runningMaxClusterAlgorithmMass; 195 G4int runningMaxClusterAlgorithmMass;
196 196
197 G4int nConsideredMax; 197 G4int nConsideredMax;
198 G4int nConsidered; 198 G4int nConsidered;
199 199
200 /** \brief Array of considered cluster par 200 /** \brief Array of considered cluster partners
201 * 201 *
202 * A dynamical array of ConsideredPartner 202 * A dynamical array of ConsideredPartner objects is allocated on this
203 * variable and filled with pointers to nu 203 * variable and filled with pointers to nucleons which are eligible for
204 * clustering. We used to use a ParticleLi 204 * clustering. We used to use a ParticleList for this purpose, but this
205 * made it very cumbersome to check whethe 205 * made it very cumbersome to check whether nucleons had already been
206 * included in the running configuration. 206 * included in the running configuration. Using an array of Particle*
207 * coupled with a boolean mask (\see{isInR 207 * coupled with a boolean mask (\see{isInRunningConfiguration}) reduces the
208 * overhead by a large amount. Running ti 208 * overhead by a large amount. Running times for 1-GeV p+Pb208 went down
209 * by almost 30% (!). 209 * by almost 30% (!).
210 * 210 *
211 * Lesson learnt: when you need speed, not 211 * Lesson learnt: when you need speed, nothing beats a good ol' array.
212 */ 212 */
213 ConsideredPartner *consideredPartners; 213 ConsideredPartner *consideredPartners;
214 214
215 /** \brief Array of flags for nucleons in 215 /** \brief Array of flags for nucleons in the running configuration
216 * 216 *
217 * Clustering partners that are already us 217 * Clustering partners that are already used in the running cluster
218 * configuration are flagged as "true" in 218 * configuration are flagged as "true" in this array.
219 */ 219 */
220 G4bool *isInRunningConfiguration; 220 G4bool *isInRunningConfiguration;
221 221
222 /** \brief Best cluster configuration 222 /** \brief Best cluster configuration
223 * 223 *
224 * This array contains pointers to the nuc 224 * This array contains pointers to the nucleons which make up the best
225 * cluster configuration that has been fou 225 * cluster configuration that has been found so far.
226 */ 226 */
227 Particle *candidateConfiguration[ParticleT 227 Particle *candidateConfiguration[ParticleTable::maxClusterMass];
228 228
229 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE 229 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
230 typedef std::set<Hashing::HashType> HashCo 230 typedef std::set<Hashing::HashType> HashContainer;
231 typedef HashContainer::iterator HashIterat 231 typedef HashContainer::iterator HashIterator;
232 232
233 /// \brief Array of containers for configu 233 /// \brief Array of containers for configurations that have already been checked
234 HashContainer checkedConfigurations[Partic 234 HashContainer checkedConfigurations[ParticleTable::maxClusterMass-2];
235 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN 235 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
236 /** \brief Class for storing and comparing 236 /** \brief Class for storing and comparing sorted nucleon configurations
237 * 237 *
238 * This class is actually just a wrapper a 238 * This class is actually just a wrapper around an array of Particle*
239 * pointers. It provides a lexicographical 239 * pointers. It provides a lexicographical comparison operator
240 * (SortedNucleonConfiguration::operator<) 240 * (SortedNucleonConfiguration::operator<) for inclusion in std::set
241 * containers. 241 * containers.
242 */ 242 */
243 class SortedNucleonConfiguration { 243 class SortedNucleonConfiguration {
244 public: 244 public:
245 // Use Particle* as nucleon identifier 245 // Use Particle* as nucleon identifiers
246 typedef G4int NucleonItem; 246 typedef G4int NucleonItem;
247 247
248 /// \brief Constructor 248 /// \brief Constructor
249 SortedNucleonConfiguration() : theSize 249 SortedNucleonConfiguration() : theSize(0), nucleons(NULL) {}
250 250
251 /// \brief Copy constructor 251 /// \brief Copy constructor
252 SortedNucleonConfiguration(const Sorte 252 SortedNucleonConfiguration(const SortedNucleonConfiguration &rhs) :
253 theSize(rhs.theSize), 253 theSize(rhs.theSize),
254 nucleons(new NucleonItem[theSize]) 254 nucleons(new NucleonItem[theSize])
255 { 255 {
256 std::copy(rhs.nucleons, rhs.nucleons+t 256 std::copy(rhs.nucleons, rhs.nucleons+theSize, nucleons);
257 } 257 }
258 258
259 /// \brief Destructor 259 /// \brief Destructor
260 ~SortedNucleonConfiguration() { 260 ~SortedNucleonConfiguration() {
261 delete [] nucleons; 261 delete [] nucleons;
262 } 262 }
263 263
264 /// \brief Helper method for the assig 264 /// \brief Helper method for the assignment operator
265 void swap(SortedNucleonConfiguration & 265 void swap(SortedNucleonConfiguration &rhs) {
266 std::swap(theSize, rhs.theSize); 266 std::swap(theSize, rhs.theSize);
267 std::swap(nucleons, rhs.nucleons); 267 std::swap(nucleons, rhs.nucleons);
268 } 268 }
269 269
270 /// \brief Assignment operator 270 /// \brief Assignment operator
271 SortedNucleonConfiguration &operator=( 271 SortedNucleonConfiguration &operator=(const SortedNucleonConfiguration &rhs) {
272 SortedNucleonConfiguration tempConfi 272 SortedNucleonConfiguration tempConfig(rhs);
273 swap(tempConfig); 273 swap(tempConfig);
274 return *this; 274 return *this;
275 } 275 }
276 276
277 /** \brief Order operator for SortedNu 277 /** \brief Order operator for SortedNucleonConfiguration
278 * 278 *
279 * The comparison is done lexicographi 279 * The comparison is done lexicographically (i.e. from the first
280 * element to the last). 280 * element to the last).
281 */ 281 */
282 G4bool operator<(const SortedNucleonCo 282 G4bool operator<(const SortedNucleonConfiguration &rhs) const {
283 // assert(theSize==rhs.theSize); 283 // assert(theSize==rhs.theSize);
284 return std::lexicographical_compare( 284 return std::lexicographical_compare(nucleons, nucleons+theSize, rhs.nucleons, rhs.nucleons+theSize);
285 } 285 }
286 286
287 /// \brief Fill configuration with arr 287 /// \brief Fill configuration with array of NucleonItem
288 void fill(NucleonItem *config, size_t 288 void fill(NucleonItem *config, size_t n) {
289 theSize = n; 289 theSize = n;
290 nucleons = new NucleonItem[theSize]; 290 nucleons = new NucleonItem[theSize];
291 std::copy(config, config+theSize, nu 291 std::copy(config, config+theSize, nucleons);
292 std::sort(nucleons, nucleons+theSize 292 std::sort(nucleons, nucleons+theSize);
293 } 293 }
294 294
295 private: 295 private:
296 /// \brief Size of the array 296 /// \brief Size of the array
297 size_t theSize; 297 size_t theSize;
298 298
299 /// \brief The real array 299 /// \brief The real array
300 NucleonItem *nucleons; 300 NucleonItem *nucleons;
301 }; 301 };
302 302
303 typedef std::set<SortedNucleonConfiguratio 303 typedef std::set<SortedNucleonConfiguration> SortedNucleonConfigurationContainer;
304 typedef SortedNucleonConfigurationContaine 304 typedef SortedNucleonConfigurationContainer::iterator SortedNucleonConfigurationIterator;
305 305
306 /// \brief Array of containers for configu 306 /// \brief Array of containers for configurations that have already been checked
307 SortedNucleonConfigurationContainer checke 307 SortedNucleonConfigurationContainer checkedConfigurations[ParticleTable::maxClusterMass-2];
308 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_I 308 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None)
309 #error Unrecognized INCL_CACHING_CLUSTERING_MO 309 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None.
310 #endif 310 #endif
311 311
312 /** \brief Maximum mass for configuration 312 /** \brief Maximum mass for configuration storage
313 * 313 *
314 * Skipping configurations becomes ineffic 314 * Skipping configurations becomes inefficient above this mass.
315 */ 315 */
316 G4int maxMassConfigurationSkipping; 316 G4int maxMassConfigurationSkipping;
317 }; 317 };
318 318
319 } 319 }
320 320
321 #endif 321 #endif
322 322