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