Geant4  10.03.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 // 20130314 M. Kelsey -- Restore null initializer and if-block for _TLS_.
45 // 20130326 M. Kelsey -- Replace _TLS_ with mutable data member buffer.
46 // 20170406 M. Kelsey -- Remove recursive tryCluster() calls (redundant),
47 // and remove use of triedClusters registry.
48 
49 #include "G4CascadeCoalescence.hh"
50 #include "G4CascadeParameters.hh"
51 #include "G4CollisionOutput.hh"
53 #include "G4InuclNuclei.hh"
54 #include "G4InuclParticle.hh"
55 #include "G4ParticleLargerBeta.hh"
56 #include "G4ParticleLargerEkin.hh"
57 #include "G4ThreeVector.hh"
58 #include <vector>
59 #include <numeric>
60 #include <algorithm>
61 #include <iterator>
62 
63 
64 // Constructor and Destructor
65 
67  : verboseLevel(verbose), thisFinalState(0), thisHadrons(0),
68  dpMaxDoublet(G4CascadeParameters::dpMaxDoublet()),
69  dpMaxTriplet(G4CascadeParameters::dpMaxTriplet()),
70  dpMaxAlpha(G4CascadeParameters::dpMaxAlpha()) {}
71 
73 
74 
75 // Final state particle list is modified directly
76 
78  if (verboseLevel)
79  G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
80 
81  thisFinalState = &finalState; // Save pointers for use in processing
82  thisHadrons = &finalState.getOutgoingParticles();
83 
84  if (verboseLevel>1) thisFinalState->printCollisionOutput(); // Before
85 
86  selectCandidates();
87  createNuclei();
88  removeNucleons();
89 
90  if (verboseLevel>1) thisFinalState->printCollisionOutput(); // After
91 }
92 
93 
94 // Scan list for possible nucleon clusters
95 
96 void G4CascadeCoalescence::selectCandidates() {
97  if (verboseLevel)
98  G4cout << " >>> G4CascadeCoalescence::selectCandidates()" << G4endl;
99 
100  allClusters.clear();
101  usedNucleons.clear();
102 
103  size_t nHad = thisHadrons->size();
104  for (size_t idx1=0; idx1<nHad; idx1++) {
105  if (!getHadron(idx1).nucleon()) continue;
106  for (size_t idx2=idx1+1; idx2<nHad; idx2++) {
107  if (!getHadron(idx2).nucleon()) continue;
108  for (size_t idx3=idx2+1; idx3<nHad; idx3++) {
109  if (!getHadron(idx3).nucleon()) continue;
110  for (size_t idx4=idx3+1; idx4<nHad; idx4++) {
111  if (!getHadron(idx4).nucleon()) continue;
112  tryClusters(idx1, idx2, idx3, idx4);
113  }
114  tryClusters(idx1, idx2, idx3); // If idx4 loop was empty
115  }
116  tryClusters(idx1, idx2); // If idx3 loop was empty
117  }
118  }
119 
120  // All potential candidates built; report statistics
121  if (verboseLevel>1) {
122  G4cout << " Found " << allClusters.size() << " candidates, using "
123  << usedNucleons.size() << " nucleons" << G4endl;
124  }
125 }
126 
127 
128 // Do combinatorics of current set of four, skip nucleons already used
129 
130 void G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2,
131  size_t idx3, size_t idx4) {
132  if (nucleonUsed(idx1) || nucleonUsed(idx2) ||
133  nucleonUsed(idx3) || nucleonUsed(idx4)) return;
134 
135  fillCluster(idx1,idx2,idx3,idx4);
136  if (verboseLevel>1) reportArgs("tryClusters",thisCluster);
137 
138  if (goodCluster(thisCluster)) {
139  allClusters.push_back(thisCluster);
140  usedNucleons.insert(idx1);
141  usedNucleons.insert(idx2);
142  usedNucleons.insert(idx3);
143  usedNucleons.insert(idx4);
144  }
145 }
146 
147 void
148 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2, size_t idx3) {
149  if (nucleonUsed(idx1) || nucleonUsed(idx2) || nucleonUsed(idx3)) return;
150 
151  fillCluster(idx1,idx2,idx3);
152  if (verboseLevel>1) reportArgs("tryClusters",thisCluster);
153 
154  if (goodCluster(thisCluster)) {
155  allClusters.push_back(thisCluster);
156  usedNucleons.insert(idx1);
157  usedNucleons.insert(idx2);
158  usedNucleons.insert(idx3);
159  }
160 }
161 
162 void
163 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2) {
164  if (nucleonUsed(idx1) || nucleonUsed(idx2)) return;
165 
166  fillCluster(idx1,idx2);
167  if (verboseLevel>1) reportArgs("tryClusters",thisCluster);
168 
169  if (goodCluster(thisCluster)) {
170  allClusters.push_back(thisCluster);
171  usedNucleons.insert(idx1);
172  usedNucleons.insert(idx2);
173  }
174 }
175 
176 
177 // Process list of candidate clusters into light ions
178 
179 void G4CascadeCoalescence::createNuclei() {
180  if (verboseLevel)
181  G4cout << " >>> G4CascadeCoalescence::createNuclei()" << G4endl;
182 
183  usedNucleons.clear();
184 
185  for (size_t i=0; i<allClusters.size(); i++) {
186  if (verboseLevel>1) G4cout << " attempting candidate #" << i << G4endl;
187 
188  const ClusterCandidate& cand = allClusters[i];
189  if (makeLightIon(cand)) { // Success, put ion in output
190  thisFinalState->addOutgoingNucleus(thisLightIon);
191  usedNucleons.insert(cand.begin(), cand.end());
192  }
193  }
194 }
195 
196 
197 // Remove nucleons indexed in "usedNucleons" from output
198 
199 void G4CascadeCoalescence::removeNucleons() {
200  if (verboseLevel>1)
201  G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
202 
203  // Remove nucleons from output from last to first (to preserve indexing)
204  std::set<size_t>::reverse_iterator usedIter;
205  for (usedIter = usedNucleons.rbegin(); usedIter != usedNucleons.rend(); ++usedIter)
206  thisFinalState->removeOutgoingParticle(*usedIter);
207 
208  usedNucleons.clear();
209 }
210 
211 
212 // Compute momentum of whole cluster
213 
215 G4CascadeCoalescence::getClusterMomentum(const ClusterCandidate& aCluster) const {
216  pCluster.set(0.,0.,0.,0.);
217  for (size_t i=0; i<aCluster.size(); i++)
218  pCluster += getHadron(aCluster[i]).getMomentum();
219 
220  return pCluster;
221 }
222 
223 
224 // Determine magnitude of largest momentum in CM frame
225 
226 G4double G4CascadeCoalescence::maxDeltaP(const ClusterCandidate& aCluster) const {
227  if (verboseLevel>1) reportArgs("maxDeltaP", aCluster);
228 
229  G4ThreeVector boost = getClusterMomentum(aCluster).boostVector();
230 
231  G4double dp, maxDP = -1.;
232  for (size_t i=0; i<aCluster.size(); i++) {
233  const G4InuclElementaryParticle& nucl = getHadron(aCluster[i]);
234 
235  // NOTE: getMomentum() returns by value, event kinematics are not changed
236  if ((dp = nucl.getMomentum().boost(-boost).rho()) > maxDP) maxDP = dp;
237  }
238 
239  if (verboseLevel>1) G4cout << " maxDP = " << maxDP << G4endl;
240 
241  return maxDP;
242 }
243 
244 
245 // Compute "cluster type code" as sum of nucleon codes
246 
247 G4int G4CascadeCoalescence::
248 clusterType(const ClusterCandidate& aCluster) const {
249  G4int type = 0;
250  for (size_t i=0; i<aCluster.size(); i++) {
251  const G4InuclElementaryParticle& had = getHadron(aCluster[i]);
252  type += had.nucleon() ? had.type() : 0;
253  }
254 
255  return type;
256 }
257 
258 
259 // Create cluster candidate with listed indices
260 
261 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2) {
262  thisCluster.clear();
263  thisCluster.push_back(idx1);
264  thisCluster.push_back(idx2);
265 }
266 
267 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3) {
268  thisCluster.clear();
269  thisCluster.push_back(idx1);
270  thisCluster.push_back(idx2);
271  thisCluster.push_back(idx3);
272 }
273 
274 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3,
275  size_t idx4) {
276  thisCluster.clear();
277  thisCluster.push_back(idx1);
278  thisCluster.push_back(idx2);
279  thisCluster.push_back(idx3);
280  thisCluster.push_back(idx4);
281 }
282 
283 
284 
285 // Make sure all candidates in cluster are nucleons
286 
287 bool G4CascadeCoalescence::allNucleons(const ClusterCandidate& clus) const {
288  bool result = true;
289  for (size_t i=0; i<clus.size(); i++)
290  result &= getHadron(clus[0]).nucleon();
291  return result;
292 }
293 
294 
295 // Determine if collection of nucleons can form a light ion
296 
297 bool G4CascadeCoalescence::goodCluster(const ClusterCandidate& clus) const {
298  if (verboseLevel>2) reportArgs("goodCluster?",clus);
299 
300  if (!allNucleons(clus)) return false;
301 
302  if (clus.size() == 2) // Deuterons (pn)
303  return (clusterType(clus) == 3 && maxDeltaP(clus) < dpMaxDoublet);
304 
305  if (clus.size() == 3) // Tritons or He-3
306  return ((clusterType(clus) == 4 || clusterType(clus) == 5) // ppn OR pnn
307  && maxDeltaP(clus) < dpMaxTriplet);
308 
309  if (clus.size() == 4) // Alphas (ppnn)
310  return (clusterType(clus) == 6 && maxDeltaP(clus) < dpMaxAlpha);
311 
312  return false;
313 }
314 
315 
316 
317 // Convert candidate nucleon set into output nucleus
318 
319 bool G4CascadeCoalescence::makeLightIon(const ClusterCandidate& aCluster) {
320  if (verboseLevel>1) reportArgs("makeLightIon",aCluster);
321 
322  thisLightIon.clear(); // Initialize nucleus buffer before filling
323 
324  if (aCluster.size()<2) return false; // Sanity check
325 
326  G4int A = aCluster.size();
327  G4int Z = -1;
328 
329  G4int type = clusterType(aCluster);
330  if (A==2 && type==3) Z = 1; // Deuteron (pn)
331  if (A==3 && type==5) Z = 1; // Triton (pnn)
332  if (A==3 && type==4) Z = 2; // He-3 (ppn)
333  if (A==4 && type==6) Z = 2; // He-4/alpha (ppnn)
334 
335  if (Z < 0) return false; // Invalid cluster content
336 
337  // NOTE: Four-momentum will not be conserved due to binding energy
338  thisLightIon.fill(getClusterMomentum(aCluster), A, Z, 0.,
340 
341  if (verboseLevel>1) reportResult("makeLightIon output",thisLightIon);
342  return true;
343 }
344 
345 
346 // Report cluster arguments for validation
347 
348 void G4CascadeCoalescence::reportArgs(const G4String& name,
349  const ClusterCandidate& aCluster) const {
350  G4cout << " >>> G4CascadeCoalescence::" << name << " ";
351  std::copy(aCluster.begin(), aCluster.end(),
352  std::ostream_iterator<size_t>(G4cout, " "));
353  G4cout << G4endl;
354 
355  if (verboseLevel>2) {
356  for (size_t i=0; i<aCluster.size(); i++)
357  G4cout << getHadron(aCluster[i]) << G4endl;
358  }
359 }
360 
361 void G4CascadeCoalescence::reportResult(const G4String& name,
362  const G4InuclNuclei& nucl) const {
363  G4cout << " >>> G4CascadeCoalescence::" << name << G4endl << nucl << G4endl;
364 }
G4double G4ParticleHPJENDLHEData::G4double result
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
const XML_Char * name
Definition: expat.h:151
Hep3Vector boostVector() const
G4LorentzVector getMomentum() const
void printCollisionOutput(std::ostream &os=G4cout) const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4bool nucleon(G4int ityp)
HepLorentzVector & boost(double, double, double)
void removeOutgoingParticle(G4int index)
double rho() const
G4CascadeCoalescence(G4int verbose=0)
void set(double x, double y, double z, double t)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
#define G4endl
Definition: G4ios.hh:61
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
void FindClusters(G4CollisionOutput &finalState)
double G4double
Definition: G4Types.hh:76