| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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
39 #include "G4INCLCluster.hh"
40 #include "G4INCLRandom.hh"
41 #include "G4INCLHashing.hh"
42 #include <algorithm>
43
44 namespace G4INCL {
45
46 const G4int ClusteringModelIntercomparison::
47 const G4int ClusteringModelIntercomparison::
48
49 const G4double ClusteringModelIntercompariso
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 ClusteringModelIntercompariso
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 ClusteringModelIntercompariso
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 ClusteringModelIntercompariso
71
72 #ifdef INCL_DO_NOT_COMPILE
73 namespace {
74 G4bool cascadingFirstPredicate(ConsideredP
75 return !aPartner.isTargetSpectator;
76 }
77 }
78 #endif
79
80 Cluster* ClusteringModelIntercomparison::get
81 // Set the maximum clustering mass dynamic
82 const G4int maxClusterAlgorithmMass = nucl
83 runningMaxClusterAlgorithmMass = std::min(
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 p
98 // Default value is 1 fm.
99 const G4double transp = 1.0;
100
101 const G4double rmaxws = theNucleus->getUni
102
103 // Radius of the sphere where the leading
104 const G4double Rprime = theNucleus->getDen
105
106 // Bring the leading particle back to the
107 const G4double pk = theLeadingParticle->ge
108 const G4double cospr = theLeadingParticle-
109 const G4double arg = rmaxws*rmaxws - Rprim
110 G4double translat;
111
112 if(arg > 0.0) {
113 // coalescence sphere smaller than Rmax
114 const G4double cosmin = std::sqrt(arg)/r
115 if(cospr <= cosmin) {
116 // there is an intersection with the c
117 translat = rmaxws * cospr;
118 } else {
119 // no intersection with the coalescenc
120 translat = rmaxws * (cospr - std::sqrt
121 }
122 } else {
123 // coalescence sphere larger than Rmax
124 translat = rmaxws * cospr - std::sqrt(Rp
125 }
126
127 const ThreeVector oldLeadingParticlePositi
128 const ThreeVector leadingParticlePosition
129 const ThreeVector &leadingParticleMomentum
130 theLeadingParticle->setPosition(leadingPar
131
132 // Initialise the array of considered nucl
133 const G4int theNucleusA = theNucleus->getA
134 if(nConsideredMax < theNucleusA) {
135 delete [] consideredPartners;
136 delete [] isInRunningConfiguration;
137 nConsideredMax = 2*theNucleusA;
138 consideredPartners = new ConsideredPartn
139 isInRunningConfiguration = new G4bool [n
140 std::fill(isInRunningConfiguration,
141 isInRunningConfiguration + nCo
142 false);
143 }
144
145 // Select the subset of nucleons that will
146 // cluster production:
147 cascadingEnergyPool = 0.;
148 nConsidered = 0;
149 ParticleList const &particles = theNucleus
150 for(ParticleIter i=particles.begin(), e=pa
151 if (!(*i)->isNucleonorLambda()) continue
152 if ((*i)->getID() == theLeadingParticle-
153
154 G4double space = ((*i)->getPosition() -
155 G4double momentum = ((*i)->getMomentum()
156 G4double size = space*momentum*clusterPo
157 // Nucleons are accepted only if they ar
158 // to the leading nucleon. The selected
159 // to the running maximum cluster mass.
160 if(size < clusterPhaseSpaceCut[runningMa
161 consideredPartners[nConsidered] = *i;
162 // Keep trace of how much energy is ca
163 // is used to stop the clustering algo
164 if(!(*i)->isTargetSpectator())
165 cascadingEnergyPool += consideredPar
166 nConsidered++;
167 // Make sure we don't exceed the array
168 // assert(nConsidered<=nConsideredMax);
169 }
170 }
171 // Sort the list of considered partners so
172 // to participants. As soon as we encounte
173 // the list we know that all the remaining
174 // spectators too.
175 // std::partition(consideredPartners, consi
176
177 #ifndef INCL_CACHING_CLUSTERING_MODEL_INTERCOM
178 // Clear the sets of checked configuration
179 // We stop caching two masses short of the
180 // performance hit above
181 maxMassConfigurationSkipping = runningMaxC
182 for(G4int i=0; i<runningMaxClusterAlgorith
183 checkedConfigurations[i].clear();
184 #endif
185
186 // Initialise position, momentum and energ
187 // configuration
188 runningPositions[1] = leadingParticlePosit
189 runningMomenta[1] = leadingParticleMomentu
190 runningEnergies[1] = theLeadingParticle->g
191 runningPotentials[1] = theLeadingParticle-
192
193 // Make sure that all the elements of isIn
194 // assert(std::count(isInRunningConfiguration,
195
196 // Start the cluster search!
197 findClusterStartingFrom(1, theLeadingParti
198
199 // Again, make sure that all the elements
200 // been reset to false. This is a sanity c
201 // assert(std::count(isInRunningConfiguration,
202
203 Cluster *chosenCluster = 0;
204 if(selectedA!=0) { // A cluster was found!
205 candidateConfiguration[selectedA-1] = th
206 chosenCluster = new Cluster(candidateCo
207 candidateConfiguration + selectedA);
208 }
209
210 // Restore the original position of the le
211 theLeadingParticle->setPosition(oldLeading
212
213 return chosenCluster;
214 }
215
216 G4double ClusteringModelIntercomparison::get
217 const G4double psSpace = (p.position - run
218 const G4double psMomentum = (p.momentum*ol
219 return psSpace * psMomentum * clusterPosFa
220 }
221
222 void ClusteringModelIntercomparison::findClu
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 = clusterPhas
231
232 // Configuration caching enabled only for
233 const G4bool cachingEnabled = (newA<=maxMa
234
235 // Set the pointer to the container of cac
236 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE
237 HashContainer *theHashContainer;
238 if(cachingEnabled)
239 theHashContainer = &(checkedConfiguratio
240 else
241 theHashContainer = NULL;
242 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN
243 SortedNucleonConfigurationContainer *theCo
244 if(cachingEnabled)
245 theConfigContainer = &(checkedConfigurat
246 else
247 theConfigContainer = NULL;
248 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_I
249 #error Unrecognized INCL_CACHING_CLUSTERING_MO
250 #endif
251
252 // Minimum and maximum Z values for this m
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 al
258 if(isInRunningConfiguration[i]) continue
259
260 ConsideredPartner const &candidateNucleo
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
268 if(newZ > clusterZMaxAll || newN > clust
269 continue;
270
271 // Compute the phase space factor for a
272 // consists of the previous running clus
273 // candidate nucleon:
274 const G4double phaseSpace = getPhaseSpac
275 if(phaseSpace > phaseSpaceCut) continue;
276
277 // Store the candidate nucleon in the ru
278 runningConfiguration[oldAMinusOne] = i;
279 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE
280 Hashing::HashType configHash;
281 HashIterator aHashIter;
282 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN
283 SortedNucleonConfiguration thisConfig;
284 SortedNucleonConfigurationIterator thisC
285 #endif
286 if(cachingEnabled) {
287 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE
288 configHash = Hashing::hashConfig(runni
289 aHashIter = theHashContainer->lower_bo
290 // If we have already checked this con
291 if(aHashIter!=theHashContainer->end()
292 && !(configHash < *aHashIter))
293 continue;
294 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN
295 thisConfig.fill(runningConfiguration,o
296 thisConfigIter = theConfigContainer->l
297 // If we have already checked this con
298 if(thisConfigIter!=theConfigContainer-
299 && !(thisConfig < *thisConfigIter))
300 continue;
301 #endif
302 }
303
304 // Sum of the total energies of the clus
305 runningEnergies[newA] = runningEnergies[
306 // Sum of the potential energies of the
307 runningPotentials[newA] = runningPotenti
308
309 // Update the available cascading kineti
310 G4double oldCascadingEnergyPool = cascad
311 if(!candidateNucleon.isTargetSpectator)
312 cascadingEnergyPool -= candidateNucleo
313
314 // Check an approximate Coulomb barrier.
315 // 0.5*barrier and the remaining availab
316 // will not bring it above, reject the c
317 const G4double halfB = 0.72 * newZ *
318 theNucleus->getZ()/(theNucleus->getDen
319 const G4double tout = runningEnergies[ne
320 931.3*newA;
321 if(tout<=halfB && tout+cascadingEnergyPo
322 cascadingEnergyPool = oldCascadingEner
323 continue;
324 }
325
326 // Here the nucleon has passed all the t
327 runningPositions[newA] = (runningPositio
328 runningMomenta[newA] = runningMomenta[ol
329
330 // Add the config to the container
331 if(cachingEnabled) {
332 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTE
333 theHashContainer->insert(aHashIter, co
334 #elif defined(INCL_CACHING_CLUSTERING_MODEL_IN
335 theConfigContainer->insert(thisConfigI
336 #endif
337 }
338
339 // Set the flag that reminds us that thi
340 // in the running configuration
341 isInRunningConfiguration[i] = true;
342
343 // Keep track of the best physical clust
344 if(newZ >= ZMinForNewA && newZ <= ZMaxFo
345 // Note: sqc is real kinetic energy, n
346 const G4double sqc = KinematicsUtils::
347 runningMomenta[newA]);
348 const G4double sqct = (sqc - 2.*newZ*p
349 + ParticleTable::getRealMass(newA
350 *clusterPosFact[newA];
351
352 if(sqct < sqtot) {
353 // This is the best cluster we have
354 // kinematics.
355 sqtot = sqct;
356 selectedA = newA;
357 selectedZ = newZ;
358 selectedS = newS;
359
360 // Store the running configuration i
361 for(G4int j=0; j<oldA; ++j)
362 candidateConfiguration[j] = consid
363
364 // Sanity check on number of nucleon
365 // assert(std::count(isInRunningConfiguration,
366 }
367 }
368
369 // The method recursively calls itself f
370 if(newA < runningMaxClusterAlgorithmMass
371 findClusterStartingFrom(newA, newZ, newS);
372 }
373
374 // Reset the running configuration flag
375 isInRunningConfiguration[i] = false;
376 cascadingEnergyPool = oldCascadingEnergy
377 }
378 }
379
380 G4bool ClusteringModelIntercomparison::clust
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(mo
389 if(cosEscapeAngle < limitCosEscapeAngle)
390 return false;
391
392 return true;
393 }
394
395 }
396