00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
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
00061
00062 typedef std::vector<G4InuclElementaryParticle> hadronList;
00063 typedef hadronList::const_iterator hadronIter;
00064
00065
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
00075
00076 G4CascadeCoalescence::G4CascadeCoalescence(G4int verbose)
00077 : verboseLevel(verbose), thisFinalState(0), thisHadrons(0) {}
00078
00079 G4CascadeCoalescence::~G4CascadeCoalescence() {}
00080
00081
00082
00083
00084 void G4CascadeCoalescence::FindClusters(G4CollisionOutput& finalState) {
00085 if (verboseLevel)
00086 G4cout << " >>> G4CascadeCoalescence::FindClusters()" << G4endl;
00087
00088 thisFinalState = &finalState;
00089 thisHadrons = &finalState.getOutgoingParticles();
00090
00091 if (verboseLevel>1) thisFinalState->printCollisionOutput();
00092
00093 selectCandidates();
00094 createNuclei();
00095 removeNucleons();
00096
00097 if (verboseLevel>1) thisFinalState->printCollisionOutput();
00098 }
00099
00100
00101
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);
00123 }
00124 tryClusters(idx1, idx2);
00125 }
00126 }
00127
00128
00129 if (verboseLevel>1) {
00130 G4cout << " Found " << allClusters.size() << " candidates, using "
00131 << usedNucleons.size() << " nucleons" << G4endl;
00132 }
00133 }
00134
00135
00136
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));
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);
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));
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);
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));
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
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)) {
00229 thisFinalState->addOutgoingNucleus(thisLightIon);
00230 usedNucleons.insert(cand.begin(), cand.end());
00231 }
00232 }
00233 }
00234
00235
00236
00237
00238 void G4CascadeCoalescence::removeNucleons() {
00239 if (verboseLevel>1)
00240 G4cout << " >>> G4CascadeCoalescence::removeNucleons()" << G4endl;
00241
00242
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
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
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
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
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
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];
00306 }
00307
00308 return hash;
00309 }
00310
00311
00312
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
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
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)
00356 return (clusterType(clus) == 3 && maxDeltaP(clus) < dpMaxDoublet);
00357
00358 if (clus.size() == 3)
00359 return ((clusterType(clus) == 4 || clusterType(clus) == 5)
00360 && maxDeltaP(clus) < dpMaxTriplet);
00361
00362 if (clus.size() == 4)
00363 return (clusterType(clus) == 6 && maxDeltaP(clus) < dpMaxAlpha);
00364
00365 return false;
00366 }
00367
00368
00369
00370
00371
00372 bool G4CascadeCoalescence::makeLightIon(const ClusterCandidate& aCluster) {
00373 if (verboseLevel>1) reportArgs("makeLightIon",aCluster);
00374
00375 thisLightIon.clear();
00376
00377 if (aCluster.size()<2) return false;
00378
00379 G4int A = aCluster.size();
00380 G4int Z = -1;
00381
00382 G4int type = clusterType(aCluster);
00383 if (A==2 && type==3) Z = 1;
00384 if (A==3 && type==5) Z = 1;
00385 if (A==3 && type==4) Z = 2;
00386 if (A==4 && type==6) Z = 2;
00387
00388 if (Z < 0) return false;
00389
00390
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
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 }