Geant4  10.00.p02
G4CascadeCoalescence.cc
Go to the documentation of this file.
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 // G4CascadeCoalescence: Factory model for final-state interactions to
27 // produce light ions from cascade nucleons. The algorithm implemented
28 // here is descirbed in Section 2.3 of the LAQGSM documentation (p. 11-12)
29 // [http://lib-www.lanl.gov/la-pubs/00818645.pdf].
30 //
31 // The relative-momentum cut offs for each cluster type may be set with
32 // environment variables:
33 // DPMAX_2CLUSTER 0.090 GeV/c for deuterons
34 // DPMAX_3CLUSTER 0.108 GeV/c for tritons, He-3
35 // DPMAX_4CLUSTER 0.115 GeV/c for alphas
36 //
37 // 20110917 Michael Kelsey
38 // 20110920 M. Kelsey -- Use environment variables to set momentum cuts for
39 // tuning, replace polymorphic argument lists with use of
40 // "ClusterCandidate"
41 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
42 // 20110927 M. Kelsey -- Bug fix; missing <iterator> header, strtof -> strtod
43 // 20120822 M. Kelsey -- Move envvars to G4CascadeParameters.
44 // 20130314 M. Kelsey -- Restore null initializer and if-block for _TLS_.
45 // 20130326 M. Kelsey -- Replace _TLS_ with mutable data member buffer.
46 
47 #include "G4CascadeCoalescence.hh"
48 #include "G4CascadeParameters.hh"
49 #include "G4CollisionOutput.hh"
51 #include "G4InuclNuclei.hh"
52 #include "G4InuclParticle.hh"
53 #include "G4ParticleLargerBeta.hh"
54 #include "G4ParticleLargerEkin.hh"
55 #include "G4ThreeVector.hh"
56 #include <vector>
57 #include <numeric>
58 #include <algorithm>
59 #include <iterator>
60 
61 
62 // Short names for lists and iterators for convenience
63 
64 typedef std::vector<G4InuclElementaryParticle> hadronList;
65 typedef hadronList::const_iterator hadronIter;
66 
67 // Maximum relative momentum (in Bertini GeV) for nucleons in each cluster type
68 
70 
72 
74 
75 
76 // Constructor and Destructor
77 
79  : verboseLevel(verbose), thisFinalState(0), thisHadrons(0) {}
80 
82 
83 
84 // Final state particle list is modified directly
85 
87  if (verboseLevel)
88  G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
89 
90  thisFinalState = &finalState; // Save pointers for use in processing
91  thisHadrons = &finalState.getOutgoingParticles();
92 
94 
96  createNuclei();
98 
100 }
101 
102 
103 // Scan list for possible nucleon clusters
104 
106  if (verboseLevel)
107  G4cout << " >>> G4CascadeCoalescence::selectCandidates()" << G4endl;
108 
109  triedClusters.clear();
110  allClusters.clear();
111  usedNucleons.clear();
112 
113  size_t nHad = thisHadrons->size();
114  for (size_t idx1=0; idx1<nHad; idx1++) {
115  if (!getHadron(idx1).nucleon()) continue;
116  for (size_t idx2=idx1+1; idx2<nHad; idx2++) {
117  if (!getHadron(idx2).nucleon()) continue;
118  for (size_t idx3=idx2+1; idx3<nHad; idx3++) {
119  if (!getHadron(idx3).nucleon()) continue;
120  for (size_t idx4=idx3+1; idx4<nHad; idx4++) {
121  if (!getHadron(idx4).nucleon()) continue;
122  tryClusters(idx1, idx2, idx3, idx4);
123  }
124  tryClusters(idx1, idx2, idx3); // If idx4 loop was empty
125  }
126  tryClusters(idx1, idx2); // If idx3 loop was empty
127  }
128  }
129 
130  // All potential candidates built; report statistics
131  if (verboseLevel>1) {
132  G4cout << " Found " << allClusters.size() << " candidates, using "
133  << usedNucleons.size() << " nucleons" << G4endl;
134  }
135 }
136 
137 
138 // Do combinatorics of current set of four, skip nucleons already used
139 
140 void G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2,
141  size_t idx3, size_t idx4) {
142  fillCluster(idx1,idx2,idx3,idx4);
143  if (clusterTried(thisCluster)) return;
144 
145  if (verboseLevel>1)
146  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
147  << " " << idx3 << " " << idx4 << G4endl;
148 
149  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
150 
151  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) &&
152  !nucleonUsed(idx3) && !nucleonUsed(idx4)) {
153  if (goodCluster(thisCluster)) {
154  allClusters.push_back(thisCluster);
155  usedNucleons.insert(idx1);
156  usedNucleons.insert(idx2);
157  usedNucleons.insert(idx3);
158  usedNucleons.insert(idx4);
159  return;
160  }
161  }
162 
163  tryClusters(idx2,idx3,idx4); // Try constituent subclusters
164  tryClusters(idx1,idx3,idx4);
165  tryClusters(idx1,idx2,idx4);
166  tryClusters(idx1,idx2,idx3);
167 }
168 
169 void
170 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2, size_t idx3) {
171  fillCluster(idx1,idx2,idx3);
172  if (clusterTried(thisCluster)) return;
173 
174  if (verboseLevel>1)
175  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
176  << " " << idx3 << G4endl;
177 
178  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
179 
180  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) && !nucleonUsed(idx3)) {
181  fillCluster(idx1,idx2,idx3);
182  if (goodCluster(thisCluster)) {
183  allClusters.push_back(thisCluster);
184  usedNucleons.insert(idx1);
185  usedNucleons.insert(idx2);
186  usedNucleons.insert(idx3);
187  return;
188  }
189  }
190 
191  tryClusters(idx2,idx3); // Try constituent subclusters
192  tryClusters(idx1,idx3);
193  tryClusters(idx1,idx2);
194 }
195 
196 void
197 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2) {
198  if (nucleonUsed(idx1) || nucleonUsed(idx2)) return;
199 
200  fillCluster(idx1,idx2);
201  if (clusterTried(thisCluster)) return;
202 
203  if (verboseLevel>1)
204  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
205  << G4endl;
206 
207  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
208 
209  fillCluster(idx1,idx2);
210  if (goodCluster(thisCluster)) {
211  allClusters.push_back(thisCluster);
212  usedNucleons.insert(idx1);
213  usedNucleons.insert(idx2);
214  }
215 }
216 
217 
218 // Process list of candidate clusters into light ions
219 
221  if (verboseLevel)
222  G4cout << " >>> G4CascadeCoalescence::createNuclei()" << G4endl;
223 
224  usedNucleons.clear();
225 
226  for (size_t i=0; i<allClusters.size(); i++) {
227  if (verboseLevel>1) G4cout << " attempting candidate #" << i << G4endl;
228 
229  const ClusterCandidate& cand = allClusters[i];
230  if (makeLightIon(cand)) { // Success, put ion in output
232  usedNucleons.insert(cand.begin(), cand.end());
233  }
234  }
235 }
236 
237 
238 // Remove nucleons indexed in "usedNucleons" from output
239 
241  if (verboseLevel>1)
242  G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
243 
244  // Remove nucleons from output from last to first (to preserve indexing)
245  std::set<size_t>::reverse_iterator usedIter;
246  for (usedIter = usedNucleons.rbegin(); usedIter != usedNucleons.rend(); ++usedIter)
248 
249  usedNucleons.clear();
250 }
251 
252 
253 // Compute momentum of whole cluster
254 
257  pCluster.set(0.,0.,0.,0.);
258  for (size_t i=0; i<aCluster.size(); i++)
259  pCluster += getHadron(aCluster[i]).getMomentum();
260 
261  return pCluster;
262 }
263 
264 
265 // Determine magnitude of largest momentum in CM frame
266 
268  if (verboseLevel>1) reportArgs("maxDeltaP", aCluster);
269 
270  G4ThreeVector boost = getClusterMomentum(aCluster).boostVector();
271 
272  G4double dp, maxDP = -1.;
273  for (size_t i=0; i<aCluster.size(); i++) {
274  const G4InuclElementaryParticle& nucl = getHadron(aCluster[i]);
275 
276  // NOTE: getMomentum() returns by value, event kinematics are not changed
277  if ((dp = nucl.getMomentum().boost(-boost).rho()) > maxDP) maxDP = dp;
278  }
279 
280  if (verboseLevel>1) G4cout << " maxDP = " << maxDP << G4endl;
281 
282  return maxDP;
283 }
284 
285 
286 // Compute "cluster type code" as sum of nucleon codes
287 
289 clusterType(const ClusterCandidate& aCluster) const {
290  G4int type = 0;
291  for (size_t i=0; i<aCluster.size(); i++) {
292  const G4InuclElementaryParticle& had = getHadron(aCluster[i]);
293  type += had.nucleon() ? had.type() : 0;
294  }
295 
296  return type;
297 }
298 
299 // Compute "cluster hash" as chain of indices
300 
302 clusterHash(const ClusterCandidate& aCluster) const {
303  G4int hash = 0;
304  for (size_t i=0; i<aCluster.size(); i++) {
305  hash = hash*1000 + aCluster[i]; // FIXME: Use hadron list size instead?
306  }
307 
308  return hash;
309 }
310 
311 
312 // Create cluster candidate with listed indices
313 
314 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2) {
315  thisCluster.clear();
316  thisCluster.push_back(idx1);
317  thisCluster.push_back(idx2);
318 }
319 
320 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3) {
321  thisCluster.clear();
322  thisCluster.push_back(idx1);
323  thisCluster.push_back(idx2);
324  thisCluster.push_back(idx3);
325 }
326 
327 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3,
328  size_t idx4) {
329  thisCluster.clear();
330  thisCluster.push_back(idx1);
331  thisCluster.push_back(idx2);
332  thisCluster.push_back(idx3);
333  thisCluster.push_back(idx4);
334 }
335 
336 
337 
338 // Make sure all candidates in cluster are nucleons
339 
341  bool result = true;
342  for (size_t i=0; i<clus.size(); i++)
343  result &= getHadron(clus[0]).nucleon();
344  return result;
345 }
346 
347 
348 // Determine if collection of nucleons can form a light ion
349 
351  if (verboseLevel>2) reportArgs("goodCluster?",clus);
352 
353  if (!allNucleons(clus)) return false;
354 
355  if (clus.size() == 2) // Deuterons (pn)
356  return (clusterType(clus) == 3 && maxDeltaP(clus) < dpMaxDoublet);
357 
358  if (clus.size() == 3) // Tritons or He-3
359  return ((clusterType(clus) == 4 || clusterType(clus) == 5) // ppn OR pnn
360  && maxDeltaP(clus) < dpMaxTriplet);
361 
362  if (clus.size() == 4) // Alphas (ppnn)
363  return (clusterType(clus) == 6 && maxDeltaP(clus) < dpMaxAlpha);
364 
365  return false;
366 }
367 
368 
369 
370 // Convert candidate nucleon set into output nucleus
371 
373  if (verboseLevel>1) reportArgs("makeLightIon",aCluster);
374 
375  thisLightIon.clear(); // Initialize nucleus buffer before filling
376 
377  if (aCluster.size()<2) return false; // Sanity check
378 
379  G4int A = aCluster.size();
380  G4int Z = -1;
381 
382  G4int type = clusterType(aCluster);
383  if (A==2 && type==3) Z = 1; // Deuteron (pn)
384  if (A==3 && type==5) Z = 1; // Triton (pnn)
385  if (A==3 && type==4) Z = 2; // He-3 (ppn)
386  if (A==4 && type==6) Z = 2; // He-4/alpha (ppnn)
387 
388  if (Z < 0) return false; // Invalid cluster content
389 
390  // NOTE: Four-momentum will not be conserved due to binding energy
391  thisLightIon.fill(getClusterMomentum(aCluster), A, Z, 0.,
393 
394  if (verboseLevel>1) reportResult("makeLightIon output",thisLightIon);
395  return true;
396 }
397 
398 
399 // Report cluster arguments for validation
400 
402  const ClusterCandidate& aCluster) const {
403  G4cout << " >>> G4CascadeCoalescence::" << name << " ";
404  std::copy(aCluster.begin(), aCluster.end(),
405  std::ostream_iterator<size_t>(G4cout, " "));
406  G4cout << G4endl;
407 
408  if (verboseLevel>2) {
409  for (size_t i=0; i<aCluster.size(); i++)
410  G4cout << getHadron(aCluster[i]) << G4endl;
411  }
412 }
413 
415  const G4InuclNuclei& nucl) const {
416  G4cout << " >>> G4CascadeCoalescence::" << name << G4endl << nucl << G4endl;
417 }
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
static const G4double dpMaxTriplet
CLHEP::Hep3Vector G4ThreeVector
static G4double dpMaxAlpha()
bool allNucleons(const ClusterCandidate &clus) const
G4LorentzVector getMomentum() const
size_t clusterHash(const ClusterCandidate &clus) const
G4String name
Definition: TRTMaterials.hh:40
std::vector< size_t > ClusterCandidate
void printCollisionOutput(std::ostream &os=G4cout) const
std::set< size_t > triedClusters
std::vector< G4InuclElementaryParticle > hadronList
G4double maxDeltaP(const ClusterCandidate &aCluster) const
bool makeLightIon(const ClusterCandidate &aCluster)
std::set< size_t > usedNucleons
int G4int
Definition: G4Types.hh:78
void reportArgs(const G4String &name, const ClusterCandidate &clus) const
bool goodCluster(const ClusterCandidate &clus) const
std::vector< ClusterCandidate > allClusters
void copy(std::vector< T > &main, const std::vector< T > &data)
Definition: DicomRun.hh:91
const std::vector< G4InuclElementaryParticle > * thisHadrons
G4CollisionOutput * thisFinalState
const G4InuclElementaryParticle & getHadron(size_t idx) const
G4GLOB_DLL std::ostream G4cout
G4bool nucleon(G4int ityp)
bool nucleonUsed(size_t idx) const
G4int clusterType(const ClusterCandidate &aCluster) const
void reportResult(const G4String &name, const G4InuclNuclei &nucl) const
bool clusterTried(const ClusterCandidate &clus) const
void removeOutgoingParticle(G4int index)
static const G4double dpMaxAlpha
G4CascadeCoalescence(G4int verbose=0)
static const G4double A[nN]
hadronList::const_iterator hadronIter
static G4double dpMaxTriplet()
void tryClusters(size_t idx1, size_t idx2)
static G4double dpMaxDoublet()
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
ClusterCandidate thisCluster
#define G4endl
Definition: G4ios.hh:61
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
void FindClusters(G4CollisionOutput &finalState)
double G4double
Definition: G4Types.hh:76
static const G4double dpMaxDoublet
static unsigned long FASTCALL hash(KEY s)
Definition: xmlparse.cc:5846
G4LorentzVector getClusterMomentum(const ClusterCandidate &aCluster) const
CLHEP::HepLorentzVector G4LorentzVector
void fillCluster(size_t idx1, size_t idx2)