ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLClusteringModelIntercomparison.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLClusteringModelIntercomparison.cc
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 
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 
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 
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 
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 
71 
72 #ifdef INCL_DO_NOT_COMPILE
73  namespace {
74  G4bool cascadingFirstPredicate(ConsideredPartner const &aPartner) {
75  return !aPartner.isTargetSpectator;
76  }
77  }
78 #endif
79 
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?
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;
140  std::fill(isInRunningConfiguration,
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();
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.
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())
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
182  for(G4int i=0; i<runningMaxClusterAlgorithmMass-2; ++i) // no caching for A=1,2
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,
208  }
209 
210  // Restore the original position of the leading particle
211  theLeadingParticle->setPosition(oldLeadingParticlePosition);
212 
213  return chosenCluster;
214  }
215 
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 
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 *
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!
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)
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 
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 }