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 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00041 #include "G4INCLCascade.hh"
00042 #include "G4INCLRandom.hh"
00043 #include "G4INCLRanecu.hh"
00044 #include "G4INCLGeant4Random.hh"
00045 #include "G4INCLStandardPropagationModel.hh"
00046 #include "G4INCLParticleTable.hh"
00047 #include "G4INCLGlobalInfo.hh"
00048
00049 #include "G4INCLPauliBlocking.hh"
00050 #include "G4INCLIPauli.hh"
00051 #include "G4INCLPauliStrict.hh"
00052 #include "G4INCLPauliStandard.hh"
00053 #include "G4INCLPauliStrictStandard.hh"
00054 #include "G4INCLPauliGlobal.hh"
00055 #include "G4INCLCDPP.hh"
00056
00057 #include "G4INCLLogger.hh"
00058 #include "G4INCLGlobals.hh"
00059 #include "G4INCLNuclearDensityFactory.hh"
00060
00061 #include "G4INCLCoulombDistortion.hh"
00062 #include "G4INCLICoulomb.hh"
00063 #include "G4INCLCoulombNone.hh"
00064 #include "G4INCLCoulombNonRelativistic.hh"
00065
00066 #include "G4INCLClustering.hh"
00067 #include "G4INCLClusteringModelIntercomparison.hh"
00068 #include "G4INCLClusteringModelNone.hh"
00069
00070 #include "G4INCLIntersection.hh"
00071
00072 #include "G4INCLCrossSections.hh"
00073
00074 #include <cstring>
00075 #include <cstdlib>
00076 #include <numeric>
00077
00078 namespace G4INCL {
00079
00080 INCL::INCL(G4INCL::Config const * const config)
00081 :propagationModel(0), theA(208), theZ(82),
00082 targetInitSuccess(false),
00083 maxImpactParameter(0.),
00084 maxUniverseRadius(0.),
00085 maxInteractionDistance(0.),
00086 fixedImpactParameter(0.),
00087 theConfig(config),
00088 nucleus(NULL),
00089 minRemnantSize(4)
00090 {
00091
00092 G4INCL::Logger::setLoggerSlave(new G4INCL::LoggerSlave(theConfig->getLogFileName()));
00093 G4INCL::Logger::setVerbosityLevel(theConfig->getVerbosity());
00094
00095
00096
00097
00098 #ifdef INCLXX_IN_GEANT4_MODE
00099 G4INCL::Random::setGenerator(new G4INCL::Geant4RandomGenerator());
00100 #else
00101 G4INCL::Random::setGenerator(new G4INCL::Ranecu(theConfig->getRandomSeeds()));
00102 #endif // INCLXX_IN_GEANT4_MODE
00103
00104
00105 G4INCL::PauliType pauli = theConfig->getPauliType();
00106 if(pauli == G4INCL::StrictStatisticalPauli)
00107 G4INCL::Pauli::setBlocker(new G4INCL::PauliStrictStandard);
00108 else if(pauli == G4INCL::StatisticalPauli)
00109 G4INCL::Pauli::setBlocker(new G4INCL::PauliStandard);
00110 else if(pauli == G4INCL::StrictPauli)
00111 G4INCL::Pauli::setBlocker(new G4INCL::PauliStrict);
00112 else if(pauli == G4INCL::GlobalPauli)
00113 G4INCL::Pauli::setBlocker(new G4INCL::PauliGlobal);
00114 else if(pauli == G4INCL::NoPauli)
00115 G4INCL::Pauli::setBlocker(NULL);
00116
00117 if(theConfig->getCDPP())
00118 G4INCL::Pauli::setCDPP(new G4INCL::CDPP);
00119 else
00120 G4INCL::Pauli::setCDPP(NULL);
00121
00122
00123 G4INCL::CoulombType coulombType = theConfig->getCoulombType();
00124 if(coulombType == G4INCL::NonRelativisticCoulomb)
00125 G4INCL::CoulombDistortion::setCoulomb(new G4INCL::CoulombNonRelativistic);
00126 else
00127 G4INCL::CoulombDistortion::setCoulomb(new G4INCL::CoulombNone);
00128
00129
00130 G4INCL::ClusterAlgorithmType clusterAlgorithm = theConfig->getClusterAlgorithm();
00131 if(clusterAlgorithm == G4INCL::IntercomparisonClusterAlgorithm)
00132 G4INCL::Clustering::setClusteringModel(new G4INCL::ClusteringModelIntercomparison(theConfig));
00133 else
00134 G4INCL::Clustering::setClusteringModel(new G4INCL::ClusteringModelNone);
00135
00136
00137 G4INCL::ParticleTable::initialize(theConfig);
00138
00139
00140
00141
00142
00143
00144
00145
00146 propagationModel = new G4INCL::StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType());
00147 eventAction = new EventAction();
00148 propagationAction = new PropagationAction();
00149 avatarAction = new AvatarAction();
00150
00151 theGlobalInfo.cascadeModel = theConfig->getVersionID().c_str();
00152 theGlobalInfo.deexcitationModel = "none";
00153
00154 #ifndef INCLXX_IN_GEANT4_MODE
00155
00156 theGlobalInfo.At = theConfig->getTargetA();
00157 theGlobalInfo.Zt = theConfig->getTargetZ();
00158 const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
00159 theGlobalInfo.Ap = theSpecies.theA;
00160 theGlobalInfo.Zp = theSpecies.theZ;
00161 theGlobalInfo.Ep = theConfig->getProjectileKineticEnergy();
00162
00163 INFO(theConfig->echo() << std::endl);
00164 #endif
00165
00166 fixedImpactParameter = theConfig->getImpactParameter();
00167 }
00168
00169 INCL::~INCL() {
00170 G4INCL::Pauli::deleteBlockers();
00171 G4INCL::CoulombDistortion::deleteCoulomb();
00172 G4INCL::Random::deleteGenerator();
00173 G4INCL::Clustering::deleteClusteringModel();
00174 G4INCL::Logger::deleteLoggerSlave();
00175 G4INCL::NuclearDensityFactory::clearCache();
00176 delete avatarAction;
00177 delete propagationAction;
00178 delete eventAction;
00179 delete propagationModel;
00180 delete theConfig;
00181 }
00182
00183 G4bool INCL::prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z) {
00184 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
00185 ERROR("Unsupported target: A = " << A << " Z = " << Z << std::endl);
00186 ERROR("Target configuration rejected." << std::endl);
00187 return false;
00188 }
00189
00190
00191 initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
00192
00193
00194 theZ = Z;
00195 if(theConfig->isNaturalTarget())
00196 theA = ParticleTable::drawRandomNaturalIsotope(Z);
00197 else
00198 theA = A;
00199 initializeTarget(theA, theZ);
00200
00201
00202 maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
00203 initMaxInteractionDistance(projectileSpecies, kineticEnergy);
00204
00205
00206 theGlobalInfo.geometricCrossSection =
00207 Math::tenPi*std::pow(maxImpactParameter,2);
00208
00209
00210 if(projectileSpecies.theA > 0)
00211 minRemnantSize = std::min(theA, 4);
00212 else
00213 minRemnantSize = std::min(theA-1, 4);
00214
00215 return true;
00216 }
00217
00218 G4bool INCL::initializeTarget(const G4int A, const G4int Z) {
00219 delete nucleus;
00220
00221 nucleus = new Nucleus(A, Z, theConfig, maxUniverseRadius);
00222 nucleus->getStore()->getBook()->reset();
00223 nucleus->initializeParticles();
00224
00225 propagationModel->setNucleus(nucleus);
00226 return true;
00227 }
00228
00229 const EventInfo &INCL::processEvent(
00230 ParticleSpecies const &projectileSpecies,
00231 const G4double kineticEnergy,
00232 const G4int targetA,
00233 const G4int targetZ
00234 ) {
00235
00236 targetInitSuccess = prepareReaction(projectileSpecies, kineticEnergy, targetA, targetZ);
00237
00238 if(!targetInitSuccess) {
00239 WARN("Target initialisation failed for A=" << targetA << ", Z=" << targetZ << std::endl);
00240 theEventInfo.transparent=true;
00241 return theEventInfo;
00242 }
00243
00244 const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
00245 if(canRunCascade) {
00246 cascade();
00247 postCascade();
00248 }
00249 return theEventInfo;
00250 }
00251
00252 G4bool INCL::preCascade(ParticleSpecies const projectileSpecies, const G4double kineticEnergy) {
00253
00254 theEventInfo.reset();
00255
00256 EventInfo::eventNumber++;
00257
00258
00259 theGlobalInfo.nShots++;
00260
00261
00262 theEventInfo.projectileType = projectileSpecies.theType;
00263 theEventInfo.Ap = projectileSpecies.theA;
00264 theEventInfo.Zp = projectileSpecies.theZ;
00265 theEventInfo.Ep = kineticEnergy;
00266 theEventInfo.At = nucleus->getA();
00267 theEventInfo.Zt = nucleus->getZ();
00268
00269
00270 if(maxImpactParameter<=0.) {
00271
00272 theGlobalInfo.nTransparents++;
00273
00274
00275 theEventInfo.transparent = true;
00276
00277 return false;
00278 }
00279
00280
00281
00282 G4double impactParameter, phi;
00283 if(fixedImpactParameter<0.) {
00284 impactParameter = maxImpactParameter * std::sqrt(Random::shoot0());
00285 phi = Random::shoot() * Math::twoPi;
00286 } else {
00287 impactParameter = fixedImpactParameter;
00288 phi = 0.;
00289 }
00290
00291
00292 theEventInfo.impactParameter = impactParameter;
00293
00294 const G4double effectiveImpactParameter = propagationModel->shoot(projectileSpecies, kineticEnergy, impactParameter, phi);
00295 if(effectiveImpactParameter < 0.) {
00296
00297 theGlobalInfo.nTransparents++;
00298
00299
00300 theEventInfo.transparent = true;
00301
00302 return false;
00303 }
00304
00305
00306 theEventInfo.transparent = false;
00307 theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
00308
00309 return true;
00310 }
00311
00312 void INCL::cascade() {
00313 do {
00314
00315 propagationAction->beforePropagationAction(propagationModel);
00316
00317
00318
00319 G4INCL::IAvatar *avatar = propagationModel->propagate();
00320
00321
00322 propagationAction->afterPropagationAction(propagationModel, avatar);
00323
00324 if(avatar == 0) break;
00325
00326
00327 avatarAction->beforeAvatarAction(avatar, nucleus);
00328
00329
00330
00331
00332
00333
00334
00335
00336 G4INCL::FinalState *finalState = avatar->getFinalState();
00337
00338
00339 avatarAction->afterAvatarAction(avatar, nucleus, finalState);
00340
00341
00342 nucleus->applyFinalState(finalState);
00343
00344
00345 delete avatar;
00346 delete finalState;
00347 } while(continueCascade());
00348
00349 }
00350
00351 void INCL::postCascade() {
00352
00353 theEventInfo.stoppingTime = propagationModel->getCurrentTime();
00354
00355
00356 if(nucleus->getTryCompoundNucleus()) {
00357 DEBUG("Trying compound nucleus" << std::endl);
00358 makeCompoundNucleus();
00359
00360 #ifndef INCLXX_IN_GEANT4_MODE
00361 if(!theEventInfo.transparent) globalConservationChecks(true);
00362 #endif
00363 return;
00364 }
00365
00366 theEventInfo.transparent = nucleus->isEventTransparent();
00367
00368 if(theEventInfo.transparent) {
00369
00370 theGlobalInfo.nTransparents++;
00371 if(nucleus->isForcedTransparent())
00372 theGlobalInfo.nForcedTransparents++;
00373 ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
00374 if(projectileRemnant) {
00375
00376 projectileRemnant->deleteParticles();
00377 nucleus->deleteProjectileRemnant();
00378 nucleus->getStore()->clearIncoming();
00379 } else {
00380
00381 nucleus->getStore()->deleteIncoming();
00382 }
00383 } else {
00384
00385 theEventInfo.deltasInside = nucleus->containsDeltas();
00386
00387
00388 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
00389 theEventInfo.forcedDeltasInside = nucleus->decayInsideDeltas();
00390
00391
00392
00393
00394
00395
00396 G4INCL::CoulombDistortion::distortOut(nucleus->getStore()->getOutgoingParticles(), nucleus);
00397
00398
00399
00400 if(nucleus->getStore()->getOutgoingParticles().size()==0
00401 && nucleus->getProjectileRemnant()
00402 && nucleus->getProjectileRemnant()->getParticles().size()==0) {
00403
00404 DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << std::endl);
00405
00406 nucleus->useFusionKinematics();
00407 nucleus->deleteProjectileRemnant();
00408
00409 if(nucleus->getExcitationEnergy()<0.) {
00410
00411 WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << std::endl);
00412 theEventInfo.transparent = true;
00413 return;
00414 }
00415
00416 } else {
00417
00418
00419 nucleus->setExcitationEnergy(nucleus->computeExcitationEnergy());
00420
00421
00422
00423 theEventInfo.nUnmergedSpectators = makeProjectileRemnant();
00424
00425
00426 nucleus->computeRecoilKinematics();
00427
00428 #ifndef INCLXX_IN_GEANT4_MODE
00429
00430 globalConservationChecks(false);
00431 #endif
00432
00433
00434
00435 if(nucleus->hasRemnant()) rescaleOutgoingForRecoil();
00436
00437 }
00438
00439
00440 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
00441
00442 #ifndef INCLXX_IN_GEANT4_MODE
00443
00444 globalConservationChecks(true);
00445 #endif
00446
00447
00448 nucleus->fillEventInfo(&theEventInfo);
00449
00450
00451 if(theEventInfo.nucleonAbsorption) theGlobalInfo.nNucleonAbsorptions++;
00452 if(theEventInfo.pionAbsorption) theGlobalInfo.nPionAbsorptions++;
00453 }
00454 }
00455
00456 void INCL::makeCompoundNucleus() {
00457
00458 nucleus->getStore()->clearIncoming();
00459 nucleus->getStore()->clearOutgoing();
00460 nucleus->getProjectileRemnant()->reset();
00461 nucleus->setA(theEventInfo.At);
00462 nucleus->setZ(theEventInfo.Zt);
00463
00464
00465
00466
00467 ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
00468 ThreeVector theCNSpin = nucleus->getIncomingAngularMomentum();
00469 const G4double theTargetMass = ParticleTable::getTableMass(theEventInfo.At, theEventInfo.Zt);
00470 G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt;
00471 Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
00472 G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
00473
00474
00475 ParticleList initialProjectileComponents = theProjectileRemnant->getParticles();
00476 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
00477
00478 std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(), shuffleComponentsHelper);
00479
00480 G4bool success = true;
00481 G4bool atLeastOneNucleonEntering = false;
00482 for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(); p!=shuffledComponents.end(); ++p) {
00483
00484 Intersection intersectionUniverse(IntersectionFactory::getEarlierTrajectoryIntersection(
00485 (*p)->getPosition(),
00486 (*p)->getPropagationVelocity(),
00487 maxUniverseRadius));
00488 if(!intersectionUniverse.exists)
00489 continue;
00490
00491
00492 Intersection intersectionInteraction(IntersectionFactory::getEarlierTrajectoryIntersection(
00493 (*p)->getPosition(),
00494 (*p)->getPropagationVelocity(),
00495 maxInteractionDistance));
00496 if(intersectionInteraction.exists)
00497 atLeastOneNucleonEntering = true;
00498
00499
00500 ParticleEntryAvatar theAvatar(0.0, nucleus, *p);
00501 FinalState *fs = theAvatar.getFinalState();
00502 nucleus->applyFinalState(fs);
00503 FinalStateValidity validity = fs->getValidity();
00504 delete fs;
00505 switch(validity) {
00506 case ValidFS:
00507 case ParticleBelowFermiFS:
00508
00509 theCNA++;
00510 theCNZ += (*p)->getZ();
00511 break;
00512 case ParticleBelowZeroFS:
00513 case PauliBlockedFS:
00514 case NoEnergyConservationFS:
00515 default:
00516 success = false;
00517 break;
00518 }
00519 }
00520
00521 if(!success || !atLeastOneNucleonEntering) {
00522 DEBUG("No nucleon entering in forced CN, or some nucleons entering below zero, forcing a transparent" << std::endl);
00523 theEventInfo.transparent = true;
00524 theGlobalInfo.nTransparents++;
00525 theGlobalInfo.nForcedTransparents++;
00526 nucleus->getProjectileRemnant()->deleteParticles();
00527 nucleus->deleteProjectileRemnant();
00528 return;
00529 }
00530
00531
00532
00533
00534
00535 theCNEnergy -= theProjectileRemnant->getEnergy();
00536 theCNMomentum -= theProjectileRemnant->getMomentum();
00537
00538
00539 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
00540
00541
00542 ParticleList const &outgoing = nucleus->getStore()->getOutgoingParticles();
00543
00544 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
00545 theCNSpin -= (*i)->getAngularMomentum();
00546 }
00547
00548
00549 const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ);
00550 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
00551 if(theCNInvariantMassSquared<0.) {
00552
00553 theGlobalInfo.nTransparents++;
00554 theGlobalInfo.nForcedTransparents++;
00555 theEventInfo.transparent = true;
00556 return;
00557 }
00558 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
00559 if(theCNExcitationEnergy<0.) {
00560
00561 theGlobalInfo.nTransparents++;
00562 theGlobalInfo.nForcedTransparents++;
00563 theEventInfo.transparent = true;
00564 return;
00565 } else {
00566
00567 nucleus->setA(theCNA);
00568 nucleus->setZ(theCNZ);
00569 nucleus->setMomentum(theCNMomentum);
00570 nucleus->setEnergy(theCNEnergy);
00571 nucleus->setExcitationEnergy(theCNExcitationEnergy);
00572 nucleus->setMass(theCNMass+theCNExcitationEnergy);
00573 nucleus->setSpin(theCNSpin);
00574
00575
00576 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
00577
00578
00579 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
00580
00581
00582 nucleus->fillEventInfo(&theEventInfo);
00583 theGlobalInfo.nForcedCompoundNucleus++;
00584 }
00585 }
00586
00587 void INCL::rescaleOutgoingForRecoil() {
00588 RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
00589
00590
00591 const G4bool success = RootFinder::solve(&theRecoilFunctor, 1.0);
00592 if(success) {
00593 std::pair<G4double,G4double> theSolution = RootFinder::getSolution();
00594 theRecoilFunctor(theSolution.first);
00595 } else {
00596 WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << std::endl);
00597 }
00598
00599 }
00600
00601 #ifndef INCLXX_IN_GEANT4_MODE
00602 void INCL::globalConservationChecks(G4bool afterRecoil) {
00603 Nucleus::ConservationBalance theBalance = nucleus->getConservationBalance(theEventInfo,afterRecoil);
00604
00605
00606 const G4double pLongBalance = theBalance.momentum.getZ();
00607 const G4double pTransBalance = theBalance.momentum.perp();
00608 if(theBalance.Z != 0) {
00609 ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << std::endl);
00610 }
00611 if(theBalance.A != 0) {
00612 ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << std::endl);
00613 }
00614 G4double EThreshold, pLongThreshold, pTransThreshold;
00615 if(afterRecoil) {
00616
00617 EThreshold = 10.;
00618 pLongThreshold = 1.;
00619 pTransThreshold = 1.;
00620 } else {
00621
00622 EThreshold = 0.1;
00623 pLongThreshold = 0.1;
00624 pTransThreshold = 0.1;
00625 }
00626 if(std::abs(theBalance.energy)>EThreshold) {
00627 WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << std::endl);
00628 }
00629 if(std::abs(pLongBalance)>pLongThreshold) {
00630 WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << std::endl);
00631 }
00632 if(std::abs(pTransBalance)>pTransThreshold) {
00633 WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << std::endl);
00634 }
00635
00636
00637 theEventInfo.EBalance = theBalance.energy;
00638 theEventInfo.pLongBalance = pLongBalance;
00639 theEventInfo.pTransBalance = pTransBalance;
00640 }
00641 #endif
00642
00643 G4bool INCL::continueCascade() {
00644
00645 if(propagationModel->getCurrentTime() > propagationModel->getStoppingTime()) {
00646 DEBUG("Cascade time (" << propagationModel->getCurrentTime()
00647 << ") exceeded stopping time (" << propagationModel->getStoppingTime()
00648 << "), stopping cascade" << std::endl);
00649 return false;
00650 }
00651
00652 if(nucleus->getStore()->getBook()->getCascading()==0 &&
00653 nucleus->getStore()->getIncomingParticles().empty()) {
00654 DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << std::endl);
00655 return false;
00656 }
00657
00658 if(nucleus->getA() <= minRemnantSize) {
00659 DEBUG("Remnant size (" << nucleus->getA()
00660 << ") smaller than or equal to minimum (" << minRemnantSize
00661 << "), stopping cascade" << std::endl);
00662 return false;
00663 }
00664
00665
00666 if(nucleus->getTryCompoundNucleus()) {
00667 DEBUG("Trying to make a compound nucleus, stopping cascade" << std::endl);
00668 return false;
00669 }
00670 if(nucleus->isForcedTransparent()) {
00671 DEBUG("Forcing a transparent, stopping cascade" << std::endl);
00672 return false;
00673 }
00674
00675 return true;
00676 }
00677
00678 void INCL::finalizeGlobalInfo() {
00679 theGlobalInfo.nucleonAbsorptionCrossSection = theGlobalInfo.geometricCrossSection *
00680 ((G4double) theGlobalInfo.nNucleonAbsorptions) / ((G4double) theGlobalInfo.nShots);
00681 theGlobalInfo.pionAbsorptionCrossSection = theGlobalInfo.geometricCrossSection *
00682 ((G4double) theGlobalInfo.nPionAbsorptions) / ((G4double) theGlobalInfo.nShots);
00683 theGlobalInfo.reactionCrossSection = theGlobalInfo.geometricCrossSection *
00684 ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents)) /
00685 ((G4double) theGlobalInfo.nShots);
00686 theGlobalInfo.errorReactionCrossSection = theGlobalInfo.geometricCrossSection *
00687 std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents)) /
00688 ((G4double) theGlobalInfo.nShots);
00689 }
00690
00691 G4int INCL::makeProjectileRemnant() {
00692 G4int nUnmergedSpectators = 0;
00693
00694
00695 if(!nucleus->getProjectileRemnant())
00696 return 0;
00697
00698
00699 ParticleList geomSpectators(nucleus->getProjectileRemnant()->getParticles());
00700 ParticleList dynSpectators(nucleus->getStore()->extractDynamicalSpectators());
00701
00702
00703 if(dynSpectators.empty() && geomSpectators.empty()) {
00704 nucleus->deleteProjectileRemnant();
00705 return 0;
00706 } else if(geomSpectators.size()+dynSpectators.size()==1) {
00707 if(dynSpectators.empty()) {
00708
00709
00710 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
00711 Particle *theSpectator = geomSpectators.front();
00712 #endif
00713
00714 nucleus->moveProjectileRemnantComponentsToOutgoing();
00715 } else {
00716
00717
00718 nucleus->getStore()->addToOutgoing(dynSpectators.front());
00719 }
00720 nucleus->deleteProjectileRemnant();
00721 } else {
00722
00723 ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
00724
00725
00726 ParticleList rejected = theProjectileRemnant->addMostDynamicalSpectators(dynSpectators);
00727
00728 nUnmergedSpectators = rejected.size();
00729 nucleus->getStore()->addToOutgoing(rejected);
00730
00731
00732 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
00733
00734 }
00735
00736 return nUnmergedSpectators;
00737 }
00738
00739 void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
00740 if(projectileSpecies.theType != Composite) {
00741 maxInteractionDistance = 0.;
00742 return;
00743 }
00744
00745 const G4double projectileKineticEnergyPerNucleon = kineticEnergy/projectileSpecies.theA;
00746 const G4double r0 = NuclearDensityFactory::createDensity(theA, theZ)->getNuclearRadius();
00747
00748 maxInteractionDistance = r0 + CrossSections::interactionDistanceNN(projectileKineticEnergyPerNucleon);
00749 }
00750
00751 void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
00752 G4double rMax = 0.0;
00753 if(A==0) {
00754 IsotopicDistribution const &anIsotopicDistribution =
00755 ParticleTable::getNaturalIsotopicDistribution(Z);
00756 IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
00757 for(IsotopeIter i=theIsotopes.begin(); i!=theIsotopes.end(); ++i) {
00758 NuclearDensity *theDensity = NuclearDensityFactory::createDensity(i->theA,Z);
00759 if(!theDensity) {
00760 FATAL("NULL density in initUniverseRadius. "
00761 << "Projectile type=" << p.theType
00762 << ", A=" << p.theA
00763 << ", Z=" << p.theZ
00764 << ", kinE=" << kineticEnergy
00765 << ", target A=" << A
00766 << ", Z=" << Z
00767 << std::endl);
00768 }
00769 rMax = std::max(theDensity->getMaximumRadius(), rMax);
00770 }
00771 } else {
00772 NuclearDensity *theDensity = NuclearDensityFactory::createDensity(A,Z);
00773 if(!theDensity) {
00774 FATAL("NULL density in initUniverseRadius. "
00775 << "Projectile type=" << p.theType
00776 << ", A=" << p.theA
00777 << ", Z=" << p.theZ
00778 << ", kinE=" << kineticEnergy
00779 << ", target A=" << A
00780 << ", Z=" << Z
00781 << std::endl);
00782 }
00783 rMax = theDensity->getMaximumRadius();
00784 }
00785 if(p.theType==Composite) {
00786 maxUniverseRadius = rMax;
00787 } else if(p.theType==Proton || p.theType==Neutron) {
00788 const G4double interactionDistanceNN = CrossSections::interactionDistanceNN(kineticEnergy);
00789 if(interactionDistanceNN>CrossSections::interactionDistanceNN1GeV()) {
00790 maxUniverseRadius = rMax
00791 - CrossSections::interactionDistanceNN1GeV()
00792 + interactionDistanceNN;
00793 } else
00794 maxUniverseRadius = rMax;
00795 } else if(p.theType==PiPlus
00796 || p.theType==PiZero
00797 || p.theType==PiMinus) {
00798 const G4double interactionDistancePiN = CrossSections::interactionDistancePiN(kineticEnergy);
00799 if(interactionDistancePiN>CrossSections::interactionDistancePiN1GeV()) {
00800 maxUniverseRadius = rMax
00801 - CrossSections::interactionDistancePiN1GeV()
00802 + interactionDistancePiN;
00803 } else
00804 maxUniverseRadius = rMax;
00805 }
00806 }
00807
00808 }