Geant4  10.01.p03
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 
68 // Constructor and Destructor
69 
71  : verboseLevel(verbose), thisFinalState(0), thisHadrons(0),
72  dpMaxDoublet(G4CascadeParameters::dpMaxDoublet()),
73  dpMaxTriplet(G4CascadeParameters::dpMaxTriplet()),
74  dpMaxAlpha(G4CascadeParameters::dpMaxAlpha()) {}
75 
77 
78 
79 // Final state particle list is modified directly
80 
82  if (verboseLevel)
83  G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
84 
85  thisFinalState = &finalState; // Save pointers for use in processing
86  thisHadrons = &finalState.getOutgoingParticles();
87 
89 
91  createNuclei();
93 
95 }
96 
97 
98 // Scan list for possible nucleon clusters
99 
101  if (verboseLevel)
102  G4cout << " >>> G4CascadeCoalescence::selectCandidates()" << G4endl;
103 
104  triedClusters.clear();
105  allClusters.clear();
106  usedNucleons.clear();
107 
108  size_t nHad = thisHadrons->size();
109  for (size_t idx1=0; idx1<nHad; idx1++) {
110  if (!getHadron(idx1).nucleon()) continue;
111  for (size_t idx2=idx1+1; idx2<nHad; idx2++) {
112  if (!getHadron(idx2).nucleon()) continue;
113  for (size_t idx3=idx2+1; idx3<nHad; idx3++) {
114  if (!getHadron(idx3).nucleon()) continue;
115  for (size_t idx4=idx3+1; idx4<nHad; idx4++) {
116  if (!getHadron(idx4).nucleon()) continue;
117  tryClusters(idx1, idx2, idx3, idx4);
118  }
119  tryClusters(idx1, idx2, idx3); // If idx4 loop was empty
120  }
121  tryClusters(idx1, idx2); // If idx3 loop was empty
122  }
123  }
124 
125  // All potential candidates built; report statistics
126  if (verboseLevel>1) {
127  G4cout << " Found " << allClusters.size() << " candidates, using "
128  << usedNucleons.size() << " nucleons" << G4endl;
129  }
130 }
131 
132 
133 // Do combinatorics of current set of four, skip nucleons already used
134 
135 void G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2,
136  size_t idx3, size_t idx4) {
137  fillCluster(idx1,idx2,idx3,idx4);
138  if (clusterTried(thisCluster)) return;
139 
140  if (verboseLevel>1)
141  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
142  << " " << idx3 << " " << idx4 << G4endl;
143 
144  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
145 
146  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) &&
147  !nucleonUsed(idx3) && !nucleonUsed(idx4)) {
148  if (goodCluster(thisCluster)) {
149  allClusters.push_back(thisCluster);
150  usedNucleons.insert(idx1);
151  usedNucleons.insert(idx2);
152  usedNucleons.insert(idx3);
153  usedNucleons.insert(idx4);
154  return;
155  }
156  }
157 
158  tryClusters(idx2,idx3,idx4); // Try constituent subclusters
159  tryClusters(idx1,idx3,idx4);
160  tryClusters(idx1,idx2,idx4);
161  tryClusters(idx1,idx2,idx3);
162 }
163 
164 void
165 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2, size_t idx3) {
166  fillCluster(idx1,idx2,idx3);
167  if (clusterTried(thisCluster)) return;
168 
169  if (verboseLevel>1)
170  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
171  << " " << idx3 << G4endl;
172 
173  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
174 
175  if (!nucleonUsed(idx1) && !nucleonUsed(idx2) && !nucleonUsed(idx3)) {
176  fillCluster(idx1,idx2,idx3);
177  if (goodCluster(thisCluster)) {
178  allClusters.push_back(thisCluster);
179  usedNucleons.insert(idx1);
180  usedNucleons.insert(idx2);
181  usedNucleons.insert(idx3);
182  return;
183  }
184  }
185 
186  tryClusters(idx2,idx3); // Try constituent subclusters
187  tryClusters(idx1,idx3);
188  tryClusters(idx1,idx2);
189 }
190 
191 void
192 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2) {
193  if (nucleonUsed(idx1) || nucleonUsed(idx2)) return;
194 
195  fillCluster(idx1,idx2);
196  if (clusterTried(thisCluster)) return;
197 
198  if (verboseLevel>1)
199  G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
200  << G4endl;
201 
202  triedClusters.insert(clusterHash(thisCluster)); // Remember current attempt
203 
204  fillCluster(idx1,idx2);
205  if (goodCluster(thisCluster)) {
206  allClusters.push_back(thisCluster);
207  usedNucleons.insert(idx1);
208  usedNucleons.insert(idx2);
209  }
210 }
211 
212 
213 // Process list of candidate clusters into light ions
214 
216  if (verboseLevel)
217  G4cout << " >>> G4CascadeCoalescence::createNuclei()" << G4endl;
218 
219  usedNucleons.clear();
220 
221  for (size_t i=0; i<allClusters.size(); i++) {
222  if (verboseLevel>1) G4cout << " attempting candidate #" << i << G4endl;
223 
224  const ClusterCandidate& cand = allClusters[i];
225  if (makeLightIon(cand)) { // Success, put ion in output
227  usedNucleons.insert(cand.begin(), cand.end());
228  }
229  }
230 }
231 
232 
233 // Remove nucleons indexed in "usedNucleons" from output
234 
236  if (verboseLevel>1)
237  G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
238 
239  // Remove nucleons from output from last to first (to preserve indexing)
240  std::set<size_t>::reverse_iterator usedIter;
241  for (usedIter = usedNucleons.rbegin(); usedIter != usedNucleons.rend(); ++usedIter)
243 
244  usedNucleons.clear();
245 }
246 
247 
248 // Compute momentum of whole cluster
249 
252  pCluster.set(0.,0.,0.,0.);
253  for (size_t i=0; i<aCluster.size(); i++)
254  pCluster += getHadron(aCluster[i]).getMomentum();
255 
256  return pCluster;
257 }
258 
259 
260 // Determine magnitude of largest momentum in CM frame
261 
263  if (verboseLevel>1) reportArgs("maxDeltaP", aCluster);
264 
265  G4ThreeVector boost = getClusterMomentum(aCluster).boostVector();
266 
267  G4double dp, maxDP = -1.;
268  for (size_t i=0; i<aCluster.size(); i++) {
269  const G4InuclElementaryParticle& nucl = getHadron(aCluster[i]);
270 
271  // NOTE: getMomentum() returns by value, event kinematics are not changed
272  if ((dp = nucl.getMomentum().boost(-boost).rho()) > maxDP) maxDP = dp;
273  }
274 
275  if (verboseLevel>1) G4cout << " maxDP = " << maxDP << G4endl;
276 
277  return maxDP;
278 }
279 
280 
281 // Compute "cluster type code" as sum of nucleon codes
282 
284 clusterType(const ClusterCandidate& aCluster) const {
285  G4int type = 0;
286  for (size_t i=0; i<aCluster.size(); i++) {
287  const G4InuclElementaryParticle& had = getHadron(aCluster[i]);
288  type += had.nucleon() ? had.type() : 0;
289  }
290 
291  return type;
292 }
293 
294 // Compute "cluster hash" as chain of indices
295 
297 clusterHash(const ClusterCandidate& aCluster) const {
298  G4int hash = 0;
299  for (size_t i=0; i<aCluster.size(); i++) {
300  hash = hash*1000 + aCluster[i]; // FIXME: Use hadron list size instead?
301  }
302 
303  return hash;
304 }
305 
306 
307 // Create cluster candidate with listed indices
308 
309 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2) {
310  thisCluster.clear();
311  thisCluster.push_back(idx1);
312  thisCluster.push_back(idx2);
313 }
314 
315 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3) {
316  thisCluster.clear();
317  thisCluster.push_back(idx1);
318  thisCluster.push_back(idx2);
319  thisCluster.push_back(idx3);
320 }
321 
322 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3,
323  size_t idx4) {
324  thisCluster.clear();
325  thisCluster.push_back(idx1);
326  thisCluster.push_back(idx2);
327  thisCluster.push_back(idx3);
328  thisCluster.push_back(idx4);
329 }
330 
331 
332 
333 // Make sure all candidates in cluster are nucleons
334 
336  bool result = true;
337  for (size_t i=0; i<clus.size(); i++)
338  result &= getHadron(clus[0]).nucleon();
339  return result;
340 }
341 
342 
343 // Determine if collection of nucleons can form a light ion
344 
346  if (verboseLevel>2) reportArgs("goodCluster?",clus);
347 
348  if (!allNucleons(clus)) return false;
349 
350  if (clus.size() == 2) // Deuterons (pn)
351  return (clusterType(clus) == 3 && maxDeltaP(clus) < dpMaxDoublet);
352 
353  if (clus.size() == 3) // Tritons or He-3
354  return ((clusterType(clus) == 4 || clusterType(clus) == 5) // ppn OR pnn
355  && maxDeltaP(clus) < dpMaxTriplet);
356 
357  if (clus.size() == 4) // Alphas (ppnn)
358  return (clusterType(clus) == 6 && maxDeltaP(clus) < dpMaxAlpha);
359 
360  return false;
361 }
362 
363 
364 
365 // Convert candidate nucleon set into output nucleus
366 
368  if (verboseLevel>1) reportArgs("makeLightIon",aCluster);
369 
370  thisLightIon.clear(); // Initialize nucleus buffer before filling
371 
372  if (aCluster.size()<2) return false; // Sanity check
373 
374  G4int A = aCluster.size();
375  G4int Z = -1;
376 
377  G4int type = clusterType(aCluster);
378  if (A==2 && type==3) Z = 1; // Deuteron (pn)
379  if (A==3 && type==5) Z = 1; // Triton (pnn)
380  if (A==3 && type==4) Z = 2; // He-3 (ppn)
381  if (A==4 && type==6) Z = 2; // He-4/alpha (ppnn)
382 
383  if (Z < 0) return false; // Invalid cluster content
384 
385  // NOTE: Four-momentum will not be conserved due to binding energy
386  thisLightIon.fill(getClusterMomentum(aCluster), A, Z, 0.,
388 
389  if (verboseLevel>1) reportResult("makeLightIon output",thisLightIon);
390  return true;
391 }
392 
393 
394 // Report cluster arguments for validation
395 
397  const ClusterCandidate& aCluster) const {
398  G4cout << " >>> G4CascadeCoalescence::" << name << " ";
399  std::copy(aCluster.begin(), aCluster.end(),
400  std::ostream_iterator<size_t>(G4cout, " "));
401  G4cout << G4endl;
402 
403  if (verboseLevel>2) {
404  for (size_t i=0; i<aCluster.size(); i++)
405  G4cout << getHadron(aCluster[i]) << G4endl;
406  }
407 }
408 
410  const G4InuclNuclei& nucl) const {
411  G4cout << " >>> G4CascadeCoalescence::" << name << G4endl << nucl << G4endl;
412 }
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
CLHEP::Hep3Vector G4ThreeVector
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)
G4CascadeCoalescence(G4int verbose=0)
static const G4double A[nN]
hadronList::const_iterator hadronIter
void tryClusters(size_t idx1, size_t idx2)
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 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)