Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
45 #include "G4CascadeCoalescence.hh"
46 #include "G4CascadeParameters.hh"
47 #include "G4CollisionOutput.hh"
49 #include "G4InuclNuclei.hh"
50 #include "G4InuclParticle.hh"
51 #include "G4ParticleLargerBeta.hh"
52 #include "G4ParticleLargerEkin.hh"
53 #include "G4ThreeVector.hh"
54 #include <vector>
55 #include <numeric>
56 #include <algorithm>
57 #include <iterator>
58 
59 
60 // Short names for lists and iterators for convenience
61 
62 typedef std::vector<G4InuclElementaryParticle> hadronList;
63 typedef hadronList::const_iterator hadronIter;
64 
65 // Maximum relative momentum (in Bertini GeV) for nucleons in each cluster type
66 
67 const G4double G4CascadeCoalescence::dpMaxDoublet = G4CascadeParameters::dpMaxDoublet();
68 
69 const G4double G4CascadeCoalescence::dpMaxTriplet = G4CascadeParameters::dpMaxTriplet();
70 
71 const G4double G4CascadeCoalescence::dpMaxAlpha = G4CascadeParameters::dpMaxAlpha();
72 
73 
74 // Constructor and Destructor
75 
77  : verboseLevel(verbose), thisFinalState(0), thisHadrons(0) {}
78 
80 
81 
82 // Final state particle list is modified directly
83 
85  if (verboseLevel)
86  G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
87 
88  thisFinalState = &finalState; // Save pointers for use in processing
89  thisHadrons = &finalState.getOutgoingParticles();
90 
91  if (verboseLevel>1) thisFinalState->printCollisionOutput(); // Before
92 
93  selectCandidates();
94  createNuclei();
95  removeNucleons();
96 
97  if (verboseLevel>1) thisFinalState->printCollisionOutput(); // After
98 }
99 
100 
101 // Scan list for possible nucleon clusters
102 
103 void G4CascadeCoalescence::selectCandidates() {
104  if (verboseLevel)
105  G4cout << " >>> G4CascadeCoalescence::selectCandidates()" << G4endl;
106 
107  triedClusters.clear();
108  allClusters.clear();
109  usedNucleons.clear();
110 
111  size_t nHad = thisHadrons->size();
112  for (size_t idx1=0; idx1<nHad; idx1++) {
113  if (!getHadron(idx1).nucleon()) continue;
114  for (size_t idx2=idx1+1; idx2<nHad; idx2++) {
115  if (!getHadron(idx2).nucleon()) continue;
116  for (size_t idx3=idx2+1; idx3<nHad; idx3++) {
117  if (!getHadron(idx3).nucleon()) continue;
118  for (size_t idx4=idx3+1; idx4<nHad; idx4++) {
119  if (!getHadron(idx4).nucleon()) continue;
120  tryClusters(idx1, idx2, idx3, idx4);
121  }
122  tryClusters(idx1, idx2, idx3); // If idx4 loop was empty
123  }
124  tryClusters(idx1, idx2); // If idx3 loop was empty
125  }
126  }
127 
128  // All potential candidates built; report statistics
129  if (verboseLevel>1) {
130  G4cout << " Found " << allClusters.size() << " candidates, using "
131  << usedNucleons.size() << " nucleons" << G4endl;
132  }
133 }
134 
135 
136 // Do combinatorics of current set of four, skip nucleons already used
137 
138 void G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2,
139  size_t idx3, size_t idx4) {
140  fillCluster(idx1,idx2,idx3,idx4);
141  if (clusterTried(thisCluster)) return;
142 
143  if (verboseLevel>1)
144  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
145  << " " << idx3 << " " << idx4 << G4endl;
146 
147  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
148 
149  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) &&
150  !nucleonUsed(idx3) && !nucleonUsed(idx4)) {
151  if (goodCluster(thisCluster)) {
152  allClusters.push_back(thisCluster);
153  usedNucleons.insert(idx1);
154  usedNucleons.insert(idx2);
155  usedNucleons.insert(idx3);
156  usedNucleons.insert(idx4);
157  return;
158  }
159  }
160 
161  tryClusters(idx2,idx3,idx4); // Try constituent subclusters
162  tryClusters(idx1,idx3,idx4);
163  tryClusters(idx1,idx2,idx4);
164  tryClusters(idx1,idx2,idx3);
165 }
166 
167 void
168 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2, size_t idx3) {
169  fillCluster(idx1,idx2,idx3);
170  if (clusterTried(thisCluster)) return;
171 
172  if (verboseLevel>1)
173  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
174  << " " << idx3 << G4endl;
175 
176  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
177 
178  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) && !nucleonUsed(idx3)) {
179  fillCluster(idx1,idx2,idx3);
180  if (goodCluster(thisCluster)) {
181  allClusters.push_back(thisCluster);
182  usedNucleons.insert(idx1);
183  usedNucleons.insert(idx2);
184  usedNucleons.insert(idx3);
185  return;
186  }
187  }
188 
189  tryClusters(idx2,idx3); // Try constituent subclusters
190  tryClusters(idx1,idx3);
191  tryClusters(idx1,idx2);
192 }
193 
194 void
195 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2) {
196  if (nucleonUsed(idx1) || nucleonUsed(idx2)) return;
197 
198  fillCluster(idx1,idx2);
199  if (clusterTried(thisCluster)) return;
200 
201  if (verboseLevel>1)
202  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
203  << G4endl;
204 
205  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
206 
207  fillCluster(idx1,idx2);
208  if (goodCluster(thisCluster)) {
209  allClusters.push_back(thisCluster);
210  usedNucleons.insert(idx1);
211  usedNucleons.insert(idx2);
212  }
213 }
214 
215 
216 // Process list of candidate clusters into light ions
217 
218 void G4CascadeCoalescence::createNuclei() {
219  if (verboseLevel)
220  G4cout << " >>> G4CascadeCoalescence::createNuclei()" << G4endl;
221 
222  usedNucleons.clear();
223 
224  for (size_t i=0; i<allClusters.size(); i++) {
225  if (verboseLevel>1) G4cout << " attempting candidate #" << i << G4endl;
226 
227  const ClusterCandidate& cand = allClusters[i];
228  if (makeLightIon(cand)) { // Success, put ion in output
229  thisFinalState->addOutgoingNucleus(thisLightIon);
230  usedNucleons.insert(cand.begin(), cand.end());
231  }
232  }
233 }
234 
235 
236 // Remove nucleons indexed in "usedNucleons" from output
237 
238 void G4CascadeCoalescence::removeNucleons() {
239  if (verboseLevel>1)
240  G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
241 
242  // Remove nucleons from output from last to first (to preserve indexing)
243  std::set<size_t>::reverse_iterator usedIter;
244  for (usedIter = usedNucleons.rbegin(); usedIter != usedNucleons.rend(); ++usedIter)
245  thisFinalState->removeOutgoingParticle(*usedIter);
246 
247  usedNucleons.clear();
248 }
249 
250 
251 // Compute momentum of whole cluster
252 
254 G4CascadeCoalescence::getClusterMomentum(const ClusterCandidate& aCluster) const {
255  static G4LorentzVector ptot;
256  ptot.set(0.,0.,0.,0.);
257  for (size_t i=0; i<aCluster.size(); i++)
258  ptot += getHadron(aCluster[i]).getMomentum();
259 
260  return ptot;
261 }
262 
263 
264 // Determine magnitude of largest momentum in CM frame
265 
266 G4double G4CascadeCoalescence::maxDeltaP(const ClusterCandidate& aCluster) const {
267  if (verboseLevel>1) reportArgs("maxDeltaP", aCluster);
268 
269  G4LorentzVector pcms = getClusterMomentum(aCluster);
270  G4ThreeVector boost = pcms.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 
288 G4int G4CascadeCoalescence::
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 
301 size_t G4CascadeCoalescence::
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 
340 bool G4CascadeCoalescence::allNucleons(const ClusterCandidate& clus) const {
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 
350 bool G4CascadeCoalescence::goodCluster(const ClusterCandidate& clus) const {
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 
372 bool G4CascadeCoalescence::makeLightIon(const ClusterCandidate& aCluster) {
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 
401 void G4CascadeCoalescence::reportArgs(const G4String& name,
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 
414 void G4CascadeCoalescence::reportResult(const G4String& name,
415  const G4InuclNuclei& nucl) const {
416  G4cout << " >>> G4CascadeCoalescence::" << name << G4endl << nucl << G4endl;
417 }