G4CascadeCoalescence.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // G4CascadeCoalescence:  Factory model for final-state interactions to
00027 //   produce light ions from cascade nucleons.  The algorithm implemented
00028 //   here is descirbed in Section 2.3 of the LAQGSM documentation (p. 11-12)
00029 //   [http://lib-www.lanl.gov/la-pubs/00818645.pdf].
00030 //
00031 // The relative-momentum cut offs for each cluster type may be set with
00032 // environment variables:
00033 //      DPMAX_2CLUSTER          0.090 GeV/c for deuterons
00034 //      DPMAX_3CLUSTER          0.108 GeV/c for tritons, He-3
00035 //      DPMAX_4CLUSTER          0.115 GeV/c for alphas
00036 //
00037 // 20110917  Michael Kelsey
00038 // 20110920  M. Kelsey -- Use environment variables to set momentum cuts for
00039 //           tuning, replace polymorphic argument lists with use of
00040 //           "ClusterCandidate"
00041 // 20110922  M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
00042 // 20110927  M. Kelsey -- Bug fix; missing <iterator> header, strtof -> strtod
00043 // 20120822  M. Kelsey -- Move envvars to G4CascadeParameters.
00044 
00045 #include "G4CascadeCoalescence.hh"
00046 #include "G4CascadeParameters.hh"
00047 #include "G4CollisionOutput.hh"
00048 #include "G4InuclElementaryParticle.hh"
00049 #include "G4InuclNuclei.hh"
00050 #include "G4InuclParticle.hh"
00051 #include "G4ParticleLargerBeta.hh"
00052 #include "G4ParticleLargerEkin.hh"
00053 #include "G4ThreeVector.hh"
00054 #include <vector>
00055 #include <numeric>
00056 #include <algorithm>
00057 #include <iterator>
00058 
00059 
00060 // Short names for lists and iterators for convenience
00061 
00062 typedef std::vector<G4InuclElementaryParticle> hadronList;
00063 typedef hadronList::const_iterator hadronIter;
00064 
00065 // Maximum relative momentum (in Bertini GeV) for nucleons in each cluster type
00066 
00067 const G4double G4CascadeCoalescence::dpMaxDoublet = G4CascadeParameters::dpMaxDoublet();
00068 
00069 const G4double G4CascadeCoalescence::dpMaxTriplet = G4CascadeParameters::dpMaxTriplet();
00070 
00071 const G4double G4CascadeCoalescence::dpMaxAlpha = G4CascadeParameters::dpMaxAlpha();
00072 
00073 
00074 // Constructor and Destructor
00075 
00076 G4CascadeCoalescence::G4CascadeCoalescence(G4int verbose)
00077   : verboseLevel(verbose), thisFinalState(0), thisHadrons(0) {}
00078 
00079 G4CascadeCoalescence::~G4CascadeCoalescence() {}
00080 
00081 
00082 // Final state particle list is modified directly
00083 
00084 void G4CascadeCoalescence::FindClusters(G4CollisionOutput& finalState) {
00085   if (verboseLevel)
00086     G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
00087 
00088   thisFinalState = &finalState;         // Save pointers for use in processing
00089   thisHadrons = &finalState.getOutgoingParticles();
00090 
00091   if (verboseLevel>1) thisFinalState->printCollisionOutput();   // Before
00092 
00093   selectCandidates();
00094   createNuclei();
00095   removeNucleons();
00096 
00097   if (verboseLevel>1) thisFinalState->printCollisionOutput();   // After
00098 }
00099 
00100 
00101 // Scan list for possible nucleon clusters
00102 
00103 void G4CascadeCoalescence::selectCandidates() {
00104   if (verboseLevel)
00105     G4cout << " >>> G4CascadeCoalescence::selectCandidates()" << G4endl;
00106 
00107   triedClusters.clear();
00108   allClusters.clear();
00109   usedNucleons.clear();
00110 
00111   size_t nHad = thisHadrons->size();
00112   for (size_t idx1=0; idx1<nHad; idx1++) {
00113     if (!getHadron(idx1).nucleon()) continue;
00114     for (size_t idx2=idx1+1; idx2<nHad; idx2++) {
00115       if (!getHadron(idx2).nucleon()) continue;
00116       for (size_t idx3=idx2+1; idx3<nHad; idx3++) {
00117         if (!getHadron(idx3).nucleon()) continue;
00118         for (size_t idx4=idx3+1; idx4<nHad; idx4++) {
00119           if (!getHadron(idx4).nucleon()) continue;
00120           tryClusters(idx1, idx2, idx3, idx4);
00121         }
00122         tryClusters(idx1, idx2, idx3);          // If idx4 loop was empty
00123       }
00124       tryClusters(idx1, idx2);                  // If idx3 loop was empty
00125     }
00126   }
00127 
00128   // All potential candidates built; report statistics
00129   if (verboseLevel>1) {
00130     G4cout << " Found " << allClusters.size() << " candidates, using "
00131            << usedNucleons.size() << " nucleons" << G4endl;
00132   }
00133 }
00134 
00135 
00136 // Do combinatorics of current set of four, skip nucleons already used
00137 
00138 void G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2,
00139                                        size_t idx3, size_t idx4) {
00140   fillCluster(idx1,idx2,idx3,idx4);
00141   if (clusterTried(thisCluster)) return;
00142 
00143   if (verboseLevel>1)
00144     G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
00145            << " " << idx3 << " " << idx4 << G4endl;
00146 
00147   triedClusters.insert(clusterHash(thisCluster));  // Remember current attempt
00148 
00149   if (!nucleonUsed(idx1) && !nucleonUsed(idx2) &&
00150       !nucleonUsed(idx3) && !nucleonUsed(idx4)) {
00151     if (goodCluster(thisCluster)) {
00152       allClusters.push_back(thisCluster);
00153       usedNucleons.insert(idx1);
00154       usedNucleons.insert(idx2);
00155       usedNucleons.insert(idx3);
00156       usedNucleons.insert(idx4);
00157       return;
00158     }
00159   }
00160 
00161   tryClusters(idx2,idx3,idx4);          // Try constituent subclusters
00162   tryClusters(idx1,idx3,idx4);
00163   tryClusters(idx1,idx2,idx4);
00164   tryClusters(idx1,idx2,idx3);
00165 }
00166 
00167 void 
00168 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2, size_t idx3) {
00169   fillCluster(idx1,idx2,idx3);
00170   if (clusterTried(thisCluster)) return;
00171 
00172   if (verboseLevel>1)
00173     G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
00174            << " " << idx3 << G4endl;
00175 
00176   triedClusters.insert(clusterHash(thisCluster));  // Remember current attempt
00177 
00178   if (!nucleonUsed(idx1) && !nucleonUsed(idx2) && !nucleonUsed(idx3)) {
00179     fillCluster(idx1,idx2,idx3);
00180     if (goodCluster(thisCluster)) {
00181       allClusters.push_back(thisCluster);
00182       usedNucleons.insert(idx1);
00183       usedNucleons.insert(idx2);
00184       usedNucleons.insert(idx3);
00185       return;
00186     }
00187   }
00188 
00189   tryClusters(idx2,idx3);               // Try constituent subclusters
00190   tryClusters(idx1,idx3);
00191   tryClusters(idx1,idx2);
00192 }
00193 
00194 void 
00195 G4CascadeCoalescence::tryClusters(size_t idx1, size_t idx2) {
00196   if (nucleonUsed(idx1) || nucleonUsed(idx2)) return;
00197 
00198   fillCluster(idx1,idx2);
00199   if (clusterTried(thisCluster)) return;
00200 
00201   if (verboseLevel>1)
00202     G4cout << " >>> G4CascadeCoalescence::tryClusters " << idx1 << " " << idx2
00203            << G4endl;
00204 
00205   triedClusters.insert(clusterHash(thisCluster));  // Remember current attempt
00206 
00207   fillCluster(idx1,idx2);
00208   if (goodCluster(thisCluster)) {
00209     allClusters.push_back(thisCluster);
00210     usedNucleons.insert(idx1);
00211     usedNucleons.insert(idx2);
00212   }
00213 }
00214 
00215 
00216 // Process list of candidate clusters into light ions
00217 
00218 void G4CascadeCoalescence::createNuclei() {
00219   if (verboseLevel)
00220     G4cout << " >>> G4CascadeCoalescence::createNuclei()" << G4endl;
00221 
00222   usedNucleons.clear();
00223 
00224   for (size_t i=0; i<allClusters.size(); i++) {
00225     if (verboseLevel>1) G4cout << " attempting candidate #" << i << G4endl;
00226 
00227     const ClusterCandidate& cand = allClusters[i];
00228     if (makeLightIon(cand)) {                   // Success, put ion in output
00229       thisFinalState->addOutgoingNucleus(thisLightIon);
00230       usedNucleons.insert(cand.begin(), cand.end());
00231     }
00232   }
00233 }
00234 
00235 
00236 // Remove nucleons indexed in "usedNucleons" from output
00237 
00238 void G4CascadeCoalescence::removeNucleons() {
00239   if (verboseLevel>1)
00240     G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
00241 
00242   // Remove nucleons from output from last to first (to preserve indexing)
00243   std::set<size_t>::reverse_iterator usedIter;
00244   for (usedIter = usedNucleons.rbegin(); usedIter != usedNucleons.rend(); ++usedIter)
00245     thisFinalState->removeOutgoingParticle(*usedIter);
00246 
00247   usedNucleons.clear();
00248 }
00249 
00250 
00251 // Compute momentum of whole cluster
00252 
00253 G4LorentzVector 
00254 G4CascadeCoalescence::getClusterMomentum(const ClusterCandidate& aCluster) const {
00255   static G4LorentzVector ptot;
00256   ptot.set(0.,0.,0.,0.);
00257   for (size_t i=0; i<aCluster.size(); i++)
00258     ptot += getHadron(aCluster[i]).getMomentum();
00259 
00260   return ptot;
00261 }
00262 
00263 
00264 // Determine magnitude of largest momentum in CM frame
00265 
00266 G4double G4CascadeCoalescence::maxDeltaP(const ClusterCandidate& aCluster) const {
00267   if (verboseLevel>1) reportArgs("maxDeltaP", aCluster);
00268 
00269   G4LorentzVector pcms = getClusterMomentum(aCluster);
00270   G4ThreeVector boost = pcms.boostVector();
00271 
00272   G4double dp, maxDP = -1.;
00273   for (size_t i=0; i<aCluster.size(); i++) {
00274     const G4InuclElementaryParticle& nucl = getHadron(aCluster[i]);
00275 
00276     // NOTE:  getMomentum() returns by value, event kinematics are not changed
00277     if ((dp = nucl.getMomentum().boost(-boost).rho()) > maxDP) maxDP = dp;
00278   }
00279 
00280   if (verboseLevel>1) G4cout << " maxDP = " << maxDP << G4endl;
00281 
00282   return maxDP;
00283 }
00284 
00285 
00286 // Compute "cluster type code" as sum of nucleon codes
00287 
00288 G4int G4CascadeCoalescence::
00289 clusterType(const ClusterCandidate& aCluster) const {
00290   G4int type = 0;
00291   for (size_t i=0; i<aCluster.size(); i++) {
00292     const G4InuclElementaryParticle& had = getHadron(aCluster[i]);
00293     type += had.nucleon() ? had.type() : 0;
00294   }
00295 
00296   return type;
00297 }
00298 
00299 // Compute "cluster hash" as chain of indices
00300 
00301 size_t G4CascadeCoalescence::
00302 clusterHash(const ClusterCandidate& aCluster) const {
00303   G4int hash = 0;
00304   for (size_t i=0; i<aCluster.size(); i++) {
00305     hash = hash*1000 + aCluster[i];     // FIXME:  Use hadron list size instead?
00306   }
00307 
00308   return hash;
00309 }
00310 
00311 
00312 // Create cluster candidate with listed indices
00313 
00314 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2) {
00315   thisCluster.clear();
00316   thisCluster.push_back(idx1);
00317   thisCluster.push_back(idx2);
00318 }
00319 
00320 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3) {
00321   thisCluster.clear();
00322   thisCluster.push_back(idx1);
00323   thisCluster.push_back(idx2);
00324   thisCluster.push_back(idx3);
00325 }
00326 
00327 void G4CascadeCoalescence::fillCluster(size_t idx1, size_t idx2, size_t idx3,
00328                                       size_t idx4) {
00329   thisCluster.clear();
00330   thisCluster.push_back(idx1);
00331   thisCluster.push_back(idx2);
00332   thisCluster.push_back(idx3);
00333   thisCluster.push_back(idx4);
00334 }
00335 
00336 
00337 
00338 // Make sure all candidates in cluster are nucleons
00339 
00340 bool G4CascadeCoalescence::allNucleons(const ClusterCandidate& clus) const {
00341   bool result = true;
00342   for (size_t i=0; i<clus.size(); i++)
00343     result &= getHadron(clus[0]).nucleon();
00344   return result;
00345 }
00346 
00347 
00348 // Determine if collection of nucleons can form a light ion
00349 
00350 bool G4CascadeCoalescence::goodCluster(const ClusterCandidate& clus) const {
00351   if (verboseLevel>2) reportArgs("goodCluster?",clus);
00352 
00353   if (!allNucleons(clus)) return false;
00354 
00355   if (clus.size() == 2)                                 // Deuterons (pn)
00356     return (clusterType(clus) == 3 && maxDeltaP(clus) < dpMaxDoublet);
00357 
00358   if (clus.size() == 3)                                 // Tritons or He-3
00359     return ((clusterType(clus) == 4 || clusterType(clus) == 5)  // ppn OR pnn
00360             && maxDeltaP(clus) < dpMaxTriplet);
00361 
00362   if (clus.size() == 4)                                 // Alphas (ppnn)
00363     return (clusterType(clus) == 6 && maxDeltaP(clus) < dpMaxAlpha);
00364 
00365   return false;
00366 }
00367 
00368 
00369 
00370 // Convert candidate nucleon set into output nucleus
00371 
00372 bool G4CascadeCoalescence::makeLightIon(const ClusterCandidate& aCluster) {
00373   if (verboseLevel>1) reportArgs("makeLightIon",aCluster);
00374 
00375   thisLightIon.clear();         // Initialize nucleus buffer before filling
00376 
00377   if (aCluster.size()<2) return false;          // Sanity check
00378 
00379   G4int A = aCluster.size();
00380   G4int Z = -1;
00381 
00382   G4int type = clusterType(aCluster);
00383   if (A==2 && type==3) Z = 1;           // Deuteron (pn)
00384   if (A==3 && type==5) Z = 1;           // Triton (pnn)
00385   if (A==3 && type==4) Z = 2;           // He-3 (ppn)
00386   if (A==4 && type==6) Z = 2;           // He-4/alpha (ppnn)
00387 
00388   if (Z < 0) return false;              // Invalid cluster content
00389 
00390   // NOTE:  Four-momentum will not be conserved due to binding energy
00391   thisLightIon.fill(getClusterMomentum(aCluster), A, Z, 0.,
00392                     G4InuclParticle::Coalescence);
00393 
00394   if (verboseLevel>1) reportResult("makeLightIon output",thisLightIon);
00395   return true;
00396 }
00397 
00398 
00399 // Report cluster arguments for validation
00400 
00401 void G4CascadeCoalescence::reportArgs(const G4String& name,
00402                                       const ClusterCandidate& aCluster) const {
00403   G4cout << " >>> G4CascadeCoalescence::" << name << " ";
00404   std::copy(aCluster.begin(), aCluster.end(),
00405             std::ostream_iterator<size_t>(G4cout, " "));
00406   G4cout << G4endl;
00407 
00408   if (verboseLevel>2) {
00409     for (size_t i=0; i<aCluster.size(); i++)
00410       G4cout << getHadron(aCluster[i]) << G4endl;
00411   }
00412 }
00413 
00414 void G4CascadeCoalescence::reportResult(const G4String& name,
00415                                         const G4InuclNuclei& nucl) const {
00416   G4cout << " >>> G4CascadeCoalescence::" << name << G4endl << nucl << G4endl;
00417 }

Generated on Mon May 27 17:47:49 2013 for Geant4 by  doxygen 1.4.7