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
00037
00038
00039
00040
00041
00042
00043
00044 #ifndef G4INCLNucleus_hh
00045 #define G4INCLNucleus_hh 1
00046
00047 #include "G4INCLGlobals.hh"
00048 #include "G4INCLLogger.hh"
00049 #include "G4INCLParticle.hh"
00050 #include "G4INCLIAvatar.hh"
00051 #include "G4INCLNucleus.hh"
00052 #include "G4INCLKinematicsUtils.hh"
00053 #include "G4INCLDecayAvatar.hh"
00054 #include "G4INCLCluster.hh"
00055 #include "G4INCLClusterDecay.hh"
00056 #include "G4INCLDeJongSpin.hh"
00057 #include "G4INCLNuclearPotentialEnergyIsospinSmooth.hh"
00058 #include "G4INCLNuclearPotentialEnergyIsospin.hh"
00059 #include "G4INCLNuclearPotentialIsospin.hh"
00060 #include "G4INCLNuclearPotentialConstant.hh"
00061 #include <iterator>
00062 #include <cstdlib>
00063 #include <sstream>
00064
00065
00066 namespace G4INCL {
00067
00068 Nucleus::Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius)
00069 : Cluster(charge,mass),
00070 theInitialZ(charge), theInitialA(mass),
00071 theNpInitial(0), theNnInitial(0),
00072 initialInternalEnergy(0.),
00073 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
00074 initialCenterOfMass(0.,0.,0.),
00075 remnant(true),
00076 blockedDelta(NULL),
00077 initialEnergy(0.),
00078 tryCN(false),
00079 forceTransparent(false),
00080 projectileZ(0),
00081 projectileA(0),
00082 theUniverseRadius(universeRadius),
00083 isNucleusNucleus(false),
00084 theProjectileRemnant(NULL),
00085 theDensity(NULL),
00086 thePotential(NULL)
00087 {
00088 PotentialType potentialType;
00089 G4bool pionPotential;
00090 if(conf) {
00091 potentialType = conf->getPotentialType();
00092 pionPotential = conf->getPionPotential();
00093 } else {
00094
00095 potentialType = IsospinPotential;
00096 pionPotential = true;
00097 }
00098 switch(potentialType) {
00099 case IsospinEnergySmoothPotential:
00100 thePotential = new NuclearPotential::NuclearPotentialEnergyIsospinSmooth(theA, theZ, pionPotential);
00101 break;
00102 case IsospinEnergyPotential:
00103 thePotential = new NuclearPotential::NuclearPotentialEnergyIsospin(theA, theZ, pionPotential);
00104 break;
00105 case IsospinPotential:
00106 thePotential = new NuclearPotential::NuclearPotentialIsospin(theA, theZ, pionPotential);
00107 break;
00108 case ConstantPotential:
00109 thePotential = new NuclearPotential::NuclearPotentialConstant(theA, theZ, pionPotential);
00110 break;
00111 default:
00112 FATAL("Unrecognized potential type at Nucleus creation." << std::endl);
00113 break;
00114 }
00115
00116 ParticleTable::setProtonSeparationEnergy(thePotential->getSeparationEnergy(Proton));
00117 ParticleTable::setNeutronSeparationEnergy(thePotential->getSeparationEnergy(Neutron));
00118
00119 theDensity = NuclearDensityFactory::createDensity(theA, theZ);
00120
00121 theParticleSampler->setPotential(thePotential);
00122 theParticleSampler->setDensity(theDensity);
00123
00124 if(theUniverseRadius<0)
00125 theUniverseRadius = theDensity->getMaximumRadius();
00126 theStore = new Store(conf);
00127 toBeUpdated.clear();
00128 }
00129
00130 Nucleus::~Nucleus() {
00131 delete theStore;
00132 delete thePotential;
00133
00134
00135 }
00136
00137 void Nucleus::initializeParticles() {
00138
00139 delete theProjectileRemnant;
00140 theProjectileRemnant = NULL;
00141
00142 Cluster::initializeParticles();
00143 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
00144 updatePotentialEnergy(*i);
00145 theStore->add(*i);
00146 }
00147 particles.clear();
00148 initialInternalEnergy = computeTotalEnergy();
00149 initialCenterOfMass = thePosition;
00150 }
00151
00152 std::string Nucleus::dump() {
00153 std::stringstream ss;
00154 ss <<"(list ;; List of participants " << std::endl;
00155 ParticleList participants = theStore->getParticipants();
00156 for(ParticleIter i = participants.begin(); i != participants.end(); ++i) {
00157 ss <<"(make-particle-avatar-map " << std::endl
00158 << (*i)->dump()
00159 << "(list ;; List of avatars in this particle" << std::endl
00160 << ")) ;; Close the list of avatars and the particle-avatar-map" << std::endl;
00161 }
00162 ss << ")" << std::endl;
00163 return ss.str();
00164 }
00165
00166 void Nucleus::applyFinalState(FinalState *finalstate) {
00167 justCreated.clear();
00168 toBeUpdated.clear();
00169 blockedDelta = NULL;
00170 G4double totalEnergy = 0.0;
00171
00172 FinalStateValidity const validity = finalstate->getValidity();
00173 if(validity == ValidFS) {
00174
00175 ParticleList const &created = finalstate->getCreatedParticles();
00176 for(ParticleIter iter = created.begin(); iter != created.end(); ++iter) {
00177 theStore->add((*iter));
00178 if(!(*iter)->isOutOfWell()) {
00179 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
00180 justCreated.push_back((*iter));
00181 }
00182 }
00183
00184 ParticleList const &deleted = finalstate->getDestroyedParticles();
00185 for(ParticleIter iter = deleted.begin(); iter != deleted.end(); ++iter) {
00186 theStore->particleHasBeenDestroyed((*iter)->getID());
00187 }
00188
00189 ParticleList const &modified = finalstate->getModifiedParticles();
00190 for(ParticleIter iter = modified.begin(); iter != modified.end(); ++iter) {
00191 theStore->particleHasBeenUpdated((*iter)->getID());
00192 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
00193 toBeUpdated.push_back((*iter));
00194 }
00195
00196 ParticleList const &out = finalstate->getOutgoingParticles();
00197 for(ParticleIter iter = out.begin(); iter != out.end(); ++iter) {
00198 if((*iter)->isCluster()) {
00199 Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
00200 ParticleList const components = clusterOut->getParticles();
00201 for(ParticleIter in = components.begin(); in != components.end(); ++in)
00202 theStore->particleHasBeenEjected((*in)->getID());
00203 } else {
00204 theStore->particleHasBeenEjected((*iter)->getID());
00205 }
00206 totalEnergy += (*iter)->getEnergy();
00207 theA -= (*iter)->getA();
00208 theZ -= (*iter)->getZ();
00209 theStore->addToOutgoing(*iter);
00210 (*iter)->setEmissionTime(theStore->getBook()->getCurrentTime());
00211 }
00212
00213 ParticleList const &entering = finalstate->getEnteringParticles();
00214 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
00215 insertParticle(*iter);
00216 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
00217 toBeUpdated.push_back((*iter));
00218 }
00219 } else if(validity == PauliBlockedFS) {
00220 blockedDelta = finalstate->getBlockedDelta();
00221 } else if(validity == ParticleBelowFermiFS) {
00222 DEBUG("A Particle is entering below the Fermi sea:" << std::endl << finalstate->print() << std::endl);
00223 tryCN = true;
00224 ParticleList const &entering = finalstate->getEnteringParticles();
00225 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
00226 insertParticle(*iter);
00227 }
00228 } else if(validity == ParticleBelowZeroFS) {
00229 DEBUG("A Particle is entering below zero energy:" << std::endl << finalstate->print() << std::endl);
00230 forceTransparent = true;
00231 ParticleList const &entering = finalstate->getEnteringParticles();
00232 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
00233 insertParticle(*iter);
00234 }
00235 }
00236
00237 if(validity==ValidFS &&
00238 std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
00239 ERROR("Energy nonconservation! Energy at the beginning of the event = "
00240 << finalstate->getTotalEnergyBeforeInteraction()
00241 <<" and after interaction = "
00242 << totalEnergy << std::endl
00243 << finalstate->print());
00244 }
00245 }
00246
00247 void Nucleus::propagateParticles(G4double ) {
00248 WARN("Useless Nucleus::propagateParticles -method called." << std::endl);
00249 }
00250
00251 G4double Nucleus::computeTotalEnergy() const {
00252 G4double totalEnergy = 0.0;
00253 ParticleList inside = theStore->getParticles();
00254 for(ParticleIter p=inside.begin(); p!=inside.end(); ++p) {
00255 if((*p)->isNucleon())
00256 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
00257 else if((*p)->isResonance())
00258 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
00259 else
00260 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
00261 }
00262 return totalEnergy;
00263 }
00264
00265 void Nucleus::computeRecoilKinematics() {
00266
00267
00268 if(theA==1) {
00269 emitInsidePions();
00270 computeOneNucleonRecoilKinematics();
00271 remnant=false;
00272 return;
00273 }
00274
00275
00276 theMomentum = incomingMomentum;
00277 theSpin = incomingAngularMomentum;
00278
00279 ParticleList outgoing = theStore->getOutgoingParticles();
00280 for(ParticleIter p=outgoing.begin(); p!=outgoing.end(); ++p)
00281 {
00282 theMomentum -= (*p)->getMomentum();
00283 theSpin -= (*p)->getAngularMomentum();
00284 }
00285
00286
00287 thePosition = computeCenterOfMass();
00288 theSpin -= (thePosition-initialCenterOfMass).vector(theMomentum);
00289
00290 setMass(ParticleTable::getTableMass(theA,theZ) + theExcitationEnergy);
00291 adjustEnergyFromMomentum();
00292 remnant=true;
00293 }
00294
00295 ThreeVector Nucleus::computeCenterOfMass() const {
00296 ThreeVector cm(0.,0.,0.);
00297 G4double totalMass = 0.0;
00298 ParticleList inside = theStore->getParticles();
00299 for(ParticleIter p=inside.begin(); p!=inside.end(); ++p) {
00300 const G4double mass = (*p)->getMass();
00301 cm += (*p)->getPosition() * mass;
00302 totalMass += mass;
00303 }
00304 cm /= totalMass;
00305 return cm;
00306 }
00307
00308 G4double Nucleus::computeExcitationEnergy() const {
00309 const G4double totalEnergy = computeTotalEnergy();
00310 const G4double separationEnergies = computeSeparationEnergyBalance();
00311
00312 return totalEnergy - initialInternalEnergy - separationEnergies;
00313 }
00314
00315 std::string Nucleus::print()
00316 {
00317 std::stringstream ss;
00318 ss << "Particles in the nucleus:" << std::endl
00319 << "Participants:" << std::endl;
00320 G4int counter = 1;
00321 ParticleList participants = theStore->getParticipants();
00322 for(ParticleIter p = participants.begin(); p != participants.end(); ++p) {
00323 ss << "index = " << counter << std::endl
00324 << (*p)->print();
00325 counter++;
00326 }
00327 ss <<"Spectators:" << std::endl;
00328 ParticleList spectators = theStore->getSpectators();
00329 for(ParticleIter p = spectators.begin(); p != spectators.end(); ++p)
00330 ss << (*p)->print();
00331 ss <<"Outgoing:" << std::endl;
00332 ParticleList outgoing = theStore->getOutgoingParticles();
00333 for(ParticleIter p = outgoing.begin(); p != outgoing.end(); ++p)
00334 ss << (*p)->print();
00335
00336 return ss.str();
00337 }
00338
00339 G4bool Nucleus::decayOutgoingDeltas() {
00340 ParticleList out = theStore->getOutgoingParticles();
00341 ParticleList deltas;
00342 for(ParticleIter i = out.begin(); i != out.end(); ++i) {
00343 if((*i)->isDelta()) deltas.push_back((*i));
00344 }
00345 if(deltas.empty()) return false;
00346
00347 for(ParticleIter i = deltas.begin(); i != deltas.end(); ++i) {
00348 DEBUG("Decay outgoing delta particle:" << std::endl
00349 << (*i)->print() << std::endl);
00350 const ThreeVector beta = -(*i)->boostVector();
00351 const G4double deltaMass = (*i)->getMass();
00352
00353
00354
00355 (*i)->setMomentum(ThreeVector());
00356 (*i)->setEnergy((*i)->getMass());
00357
00358
00359 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
00360 FinalState *fs = decay->getFinalState();
00361 Particle * const pion = fs->getCreatedParticles().front();
00362 Particle * const nucleon = fs->getModifiedParticles().front();
00363
00364
00365 const G4double decayMomentum = KinematicsUtils::momentumInCM(deltaMass,
00366 nucleon->getTableMass(),
00367 pion->getTableMass());
00368 ThreeVector newMomentum = pion->getMomentum();
00369 newMomentum *= decayMomentum / newMomentum.mag();
00370
00371 pion->setTableMass();
00372 pion->setMomentum(newMomentum);
00373 pion->adjustEnergyFromMomentum();
00374 pion->setEmissionTime(theStore->getBook()->getCurrentTime());
00375 pion->boost(beta);
00376
00377 nucleon->setTableMass();
00378 nucleon->setMomentum(-newMomentum);
00379 nucleon->adjustEnergyFromMomentum();
00380 nucleon->setEmissionTime(theStore->getBook()->getCurrentTime());
00381 nucleon->boost(beta);
00382
00383 theStore->addToOutgoing(pion);
00384
00385 delete fs;
00386 delete decay;
00387 }
00388
00389 return true;
00390 }
00391
00392 G4bool Nucleus::decayInsideDeltas() {
00393
00394
00395
00396
00397
00398
00399
00400 const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
00401 if(thePotential->hasPionPotential() && !unphysicalRemnant)
00402 return false;
00403
00404
00405 ParticleList inside = theStore->getParticles();
00406 ParticleList deltas;
00407 for(ParticleIter i = inside.begin(); i != inside.end(); ++i)
00408 if((*i)->isDelta()) deltas.push_back((*i));
00409
00410
00411 for(ParticleIter i = deltas.begin(); i != deltas.end(); ++i) {
00412 DEBUG("Decay inside delta particle:" << std::endl
00413 << (*i)->print() << std::endl);
00414
00415
00416
00417
00418 IAvatar *decay;
00419 if(unphysicalRemnant)
00420 decay = new DecayAvatar((*i), 0.0, NULL, true);
00421 else
00422 decay = new DecayAvatar((*i), 0.0, this, true);
00423 FinalState *fs = decay->getFinalState();
00424
00425
00426
00427
00428 if(fs->getValidity()==ValidFS) {
00429
00430 applyFinalState(fs);
00431 }
00432 delete fs;
00433 delete decay;
00434 }
00435
00436
00437 if(unphysicalRemnant) {
00438 DEBUG("Remnant is unphysical: Z=" << theZ << ", A=" << theA << std::endl);
00439 emitInsidePions();
00440 }
00441
00442 return true;
00443 }
00444
00445 G4bool Nucleus::decayOutgoingClusters() {
00446 ParticleList out = theStore->getOutgoingParticles();
00447 ParticleList clusters;
00448 for(ParticleIter i = out.begin(); i != out.end(); ++i) {
00449 if((*i)->isCluster()) clusters.push_back((*i));
00450 }
00451 if(clusters.empty()) return false;
00452
00453 for(ParticleIter i = clusters.begin(); i != clusters.end(); ++i) {
00454 Cluster *cluster = dynamic_cast<Cluster*>(*i);
00455 cluster->deleteParticles();
00456 ParticleList decayProducts = ClusterDecay::decay(cluster);
00457 for(ParticleIter j = decayProducts.begin(); j!=decayProducts.end(); ++j)
00458 theStore->addToOutgoing(*j);
00459 }
00460 return true;
00461 }
00462
00463 G4bool Nucleus::decayMe() {
00464
00465 if(theA<=1 || (theZ!=0 && theA!=theZ))
00466 return false;
00467
00468 ParticleList decayProducts = ClusterDecay::decay(this);
00469 for(ParticleIter j = decayProducts.begin(); j!=decayProducts.end(); ++j)
00470 theStore->addToOutgoing(*j);
00471
00472 return true;
00473 }
00474
00475 void Nucleus::emitInsidePions() {
00476
00477
00478
00479
00480 WARN("Forcing emissions of all pions in the nucleus." << std::endl);
00481
00482
00483 const G4double tinyPionEnergy = 0.1;
00484
00485
00486 ParticleList inside = theStore->getParticles();
00487 for(ParticleIter i = inside.begin(); i != inside.end(); ++i) {
00488 if((*i)->isPion()) {
00489 (*i)->setEmissionTime(theStore->getBook()->getCurrentTime());
00490
00491 const G4double theQValueCorrection = (*i)->getEmissionQValueCorrection(theA,theZ);
00492 const G4double kineticEnergyOutside = (*i)->getKineticEnergy() - (*i)->getPotentialEnergy() + theQValueCorrection;
00493 (*i)->setTableMass();
00494 if(kineticEnergyOutside > 0.0)
00495 (*i)->setEnergy((*i)->getMass()+kineticEnergyOutside);
00496 else
00497 (*i)->setEnergy((*i)->getMass()+tinyPionEnergy);
00498 (*i)->adjustMomentumFromEnergy();
00499 (*i)->setPotentialEnergy(0.);
00500 theZ -= (*i)->getZ();
00501 theStore->particleHasBeenEjected((*i)->getID());
00502 theStore->addToOutgoing(*i);
00503 }
00504 }
00505 }
00506
00507 G4bool Nucleus::isEventTransparent() const {
00508
00509
00510 if(forceTransparent)
00511 return true;
00512
00513 ParticleList const &pL = theStore->getOutgoingParticles();
00514 G4int outZ = 0, outA = 0;
00515
00516
00517
00518 for(ParticleIter p = pL.begin(); p != pL.end(); ++p ) {
00519 if( (*p)->getNumberOfCollisions() != 0 ) return false;
00520 if( (*p)->getNumberOfDecays() != 0 ) return false;
00521 outZ += (*p)->getZ();
00522 outA += (*p)->getA();
00523 }
00524
00525
00526 if(theProjectileRemnant) {
00527 outZ += theProjectileRemnant->getZ();
00528 outA += theProjectileRemnant->getA();
00529 }
00530
00531 if(outZ!=projectileZ || outA!=projectileA) return false;
00532
00533 return true;
00534
00535 }
00536
00537 void Nucleus::computeOneNucleonRecoilKinematics() {
00538
00539
00540
00541 ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << std::endl);
00542
00543
00544 theExcitationEnergy = 0.0;
00545
00546
00547 Particle *remN = theStore->getParticles().front();
00548 theA -= remN->getA();
00549 theZ -= remN->getZ();
00550 theStore->particleHasBeenEjected(remN->getID());
00551 theStore->addToOutgoing(remN);
00552 remN->setEmissionTime(theStore->getBook()->getCurrentTime());
00553
00554
00555 if(remN->isDelta()) {
00556 IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
00557 FinalState *fs = decay->getFinalState();
00558
00559 ParticleList created = fs->getCreatedParticles();
00560 for(ParticleIter j = created.begin(); j != created.end(); ++j)
00561 theStore->addToOutgoing(*j);
00562 delete fs;
00563 delete decay;
00564 }
00565
00566
00567 ParticleList outgoing = theStore->getOutgoingParticles();
00568 if(outgoing.size() == 2) {
00569
00570 DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << std::endl);
00571
00572
00573
00574 Particle *p1 = outgoing.front(), *p2 = outgoing.back();
00575 const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
00576
00577 p1->boost(aBoostVector);
00578 const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
00579 const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
00580 const G4double scale = pcm/(p1->getMomentum().mag());
00581
00582 p1->setMomentum(p1->getMomentum()*scale);
00583 p2->setMomentum(-p1->getMomentum());
00584 p1->adjustEnergyFromMomentum();
00585 p2->adjustEnergyFromMomentum();
00586
00587 p1->boost(-aBoostVector);
00588 p2->boost(-aBoostVector);
00589
00590 } else {
00591
00592 DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << std::endl);
00593
00594 const G4int maxIterations=8;
00595 G4double totalEnergy, energyScale;
00596 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
00597 ThreeVector totalMomentum, deltaP;
00598 std::vector<ThreeVector> minMomenta;
00599
00600
00601 minMomenta.reserve(outgoing.size());
00602
00603
00604 totalMomentum.setX(0.0);
00605 totalMomentum.setY(0.0);
00606 totalMomentum.setZ(0.0);
00607 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
00608 totalMomentum += (*i)->getMomentum();
00609
00610
00611 totalEnergy = 0.0;
00612 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
00613 totalEnergy += (*i)->getEnergy();
00614
00615
00616 for(G4int iterations=0; iterations < maxIterations; ++iterations) {
00617
00618
00619 oldOldOldVal = oldOldVal;
00620 oldOldVal = oldVal;
00621 oldVal = val;
00622
00623 if(iterations%2 == 0) {
00624 DEBUG("Momentum step" << std::endl);
00625
00626 deltaP = incomingMomentum - totalMomentum;
00627 G4double pOldTot = 0.0;
00628 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
00629 pOldTot += (*i)->getMomentum().mag();
00630 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
00631 const ThreeVector mom = (*i)->getMomentum();
00632 (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
00633 (*i)->adjustEnergyFromMomentum();
00634 }
00635 } else {
00636 DEBUG("Energy step" << std::endl);
00637
00638 energyScale = initialEnergy/totalEnergy;
00639 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
00640 const ThreeVector mom = (*i)->getMomentum();
00641 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
00642 if(pScale>0) {
00643 (*i)->setEnergy((*i)->getEnergy()*energyScale);
00644 (*i)->adjustMomentumFromEnergy();
00645 }
00646 }
00647 }
00648
00649
00650 totalMomentum.setX(0.0);
00651 totalMomentum.setY(0.0);
00652 totalMomentum.setZ(0.0);
00653 totalEnergy = 0.0;
00654 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
00655 totalMomentum += (*i)->getMomentum();
00656 totalEnergy += (*i)->getEnergy();
00657 }
00658
00659
00660 val = std::pow(totalEnergy - initialEnergy,2) +
00661 0.25*(totalMomentum - incomingMomentum).mag2();
00662 DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
00663
00664
00665 if(val < oldVal) {
00666 DEBUG("New minimum found, storing the particle momenta" << std::endl);
00667 minMomenta.clear();
00668 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
00669 minMomenta.push_back((*i)->getMomentum());
00670 }
00671
00672
00673 if(val > oldOldVal && oldVal > oldOldOldVal) {
00674 DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
00675 break;
00676 }
00677 }
00678
00679
00680
00681
00682
00683 DEBUG("Applying the solution" << std::endl);
00684 std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
00685 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i, ++v) {
00686 (*i)->setMomentum(*v);
00687 (*i)->adjustEnergyFromMomentum();
00688 DATABLOCK((*i)->print());
00689 }
00690
00691 }
00692
00693 }
00694
00695 void Nucleus::fillEventInfo(EventInfo *eventInfo) {
00696 eventInfo->nParticles = 0;
00697 G4bool isNucleonAbsorption = false;
00698
00699 G4bool isPionAbsorption = false;
00700
00701
00702 if(eventInfo->projectileType == PiPlus ||
00703 eventInfo->projectileType == PiMinus ||
00704 eventInfo->projectileType == PiZero) {
00705 isPionAbsorption = true;
00706 }
00707
00708
00709 eventInfo->forcedCompoundNucleus = tryCN;
00710
00711
00712 ParticleList outgoingParticles = getStore()->getOutgoingParticles();
00713
00714
00715
00716 if(outgoingParticles.size() == 0 &&
00717 (eventInfo->projectileType == Proton ||
00718 eventInfo->projectileType == Neutron)) {
00719 isNucleonAbsorption = true;
00720 }
00721
00722
00723 eventInfo->nRemnants = 0;
00724 eventInfo->history.clear();
00725
00726 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i ) {
00727
00728 if((*i)->isCluster()) {
00729 Cluster const * const c = dynamic_cast<Cluster *>(*i);
00730
00731 #ifdef INCLXX_IN_GEANT4_MODE
00732 if(!c)
00733 continue;
00734 #endif
00735 const G4double eStar = c->getExcitationEnergy();
00736 if(std::abs(eStar)>1E-10) {
00737 if(eStar<0.) {
00738 WARN("Negative excitation energy in outgoing cluster! EStar = " << eStar << std::endl);
00739 }
00740 eventInfo->ARem[eventInfo->nRemnants] = c->getA();
00741 eventInfo->ZRem[eventInfo->nRemnants] = c->getZ();
00742 eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
00743 ThreeVector remnantSpin = c->getSpin();
00744 Float_t remnantSpinMag;
00745 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) {
00746 remnantSpinMag = (G4int) (remnantSpin.mag()/PhysicalConstants::hc + 0.5);
00747 } else {
00748 remnantSpinMag = ((G4int) (remnantSpin.mag()/PhysicalConstants::hc)) + 0.5;
00749 }
00750 remnantSpin *= remnantSpinMag/remnantSpin.mag();
00751 eventInfo->JRem[eventInfo->nRemnants] = remnantSpinMag;
00752 eventInfo->jxRem[eventInfo->nRemnants] = remnantSpin.getX();
00753 eventInfo->jyRem[eventInfo->nRemnants] = remnantSpin.getY();
00754 eventInfo->jzRem[eventInfo->nRemnants] = remnantSpin.getZ();
00755 eventInfo->EKinRem[eventInfo->nRemnants] = c->getKineticEnergy();
00756 ThreeVector mom = c->getMomentum();
00757 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
00758 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
00759 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
00760 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
00761 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
00762 eventInfo->nRemnants++;
00763 continue;
00764 }
00765 }
00766
00767
00768
00769 if(isPionAbsorption) {
00770 if((*i)->isPion()) {
00771 isPionAbsorption = false;
00772 }
00773 }
00774 eventInfo->A[eventInfo->nParticles] = (*i)->getA();
00775 eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
00776 eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
00777 eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
00778 ThreeVector mom = (*i)->getMomentum();
00779 eventInfo->px[eventInfo->nParticles] = mom.getX();
00780 eventInfo->py[eventInfo->nParticles] = mom.getY();
00781 eventInfo->pz[eventInfo->nParticles] = mom.getZ();
00782 eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
00783 eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
00784 eventInfo->origin[eventInfo->nParticles] = -1;
00785 eventInfo->history.push_back("");
00786 eventInfo->nParticles++;
00787 }
00788 eventInfo->nucleonAbsorption = isNucleonAbsorption;
00789 eventInfo->pionAbsorption = isPionAbsorption;
00790 eventInfo->nCascadeParticles = eventInfo->nParticles;
00791
00792
00793 if(hasRemnant()) {
00794 eventInfo->ARem[eventInfo->nRemnants] = getA();
00795 eventInfo->ZRem[eventInfo->nRemnants] = getZ();
00796 eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
00797 if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
00798 WARN("Negative excitation energy! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << std::endl);
00799 }
00800 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) {
00801 eventInfo->JRem[eventInfo->nRemnants] = (G4int) (getSpin().mag()/PhysicalConstants::hc + 0.5);
00802 } else {
00803 eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (getSpin().mag()/PhysicalConstants::hc)) + 0.5;
00804 }
00805 eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
00806 ThreeVector mom = getMomentum();
00807 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
00808 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
00809 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
00810 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
00811 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
00812 eventInfo->nRemnants++;
00813 }
00814
00815
00816 eventInfo->nCollisions = getStore()->getBook()->getAcceptedCollisions();
00817 eventInfo->nBlockedCollisions = getStore()->getBook()->getBlockedCollisions();
00818 eventInfo->nDecays = getStore()->getBook()->getAcceptedDecays();
00819 eventInfo->nBlockedDecays = getStore()->getBook()->getBlockedDecays();
00820 eventInfo->firstCollisionTime = getStore()->getBook()->getFirstCollisionTime();
00821 eventInfo->firstCollisionXSec = getStore()->getBook()->getFirstCollisionXSec();
00822 eventInfo->nReflectionAvatars = getStore()->getBook()->getAvatars(SurfaceAvatarType);
00823 eventInfo->nCollisionAvatars = getStore()->getBook()->getAvatars(CollisionAvatarType);
00824 eventInfo->nDecayAvatars = getStore()->getBook()->getAvatars(DecayAvatarType);
00825 }
00826
00827 Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
00828 ConservationBalance theBalance;
00829
00830 theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
00831 theBalance.A = theEventInfo.Ap + theEventInfo.At;
00832
00833 theBalance.energy = getInitialEnergy();
00834 theBalance.momentum = getIncomingMomentum();
00835
00836
00837 ParticleList outgoingParticles = theStore->getOutgoingParticles();
00838 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i ) {
00839 theBalance.Z -= (*i)->getZ();
00840 theBalance.A -= (*i)->getA();
00841
00842
00843 theBalance.energy -= (*i)->getEnergy();
00844 theBalance.momentum -= (*i)->getMomentum();
00845 }
00846
00847
00848 if(hasRemnant()) {
00849 theBalance.Z -= getZ();
00850 theBalance.A -= getA();
00851 theBalance.energy -= ParticleTable::getTableMass(getA(),getZ()) +
00852 getExcitationEnergy();
00853 if(afterRecoil)
00854 theBalance.energy -= getKineticEnergy();
00855 theBalance.momentum -= getMomentum();
00856 }
00857
00858 return theBalance;
00859 }
00860
00861 void Nucleus::useFusionKinematics() {
00862 setEnergy(initialEnergy);
00863 setMomentum(incomingMomentum);
00864 setSpin(incomingAngularMomentum);
00865 theExcitationEnergy = std::sqrt(theEnergy*theEnergy-theMomentum.mag2()) - getTableMass();
00866 setMass(getTableMass() + theExcitationEnergy);
00867 }
00868
00869 void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
00870
00871 if(theProjectileRemnant->getA()>1) {
00872
00873 const G4double aMass = theProjectileRemnant->getInvariantMass();
00874 theProjectileRemnant->setMass(aMass);
00875
00876
00877 const G4double anExcitationEnergy = aMass
00878 - ParticleTable::getTableMass(theProjectileRemnant->getA(), theProjectileRemnant->getZ());
00879
00880
00881 theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
00882
00883
00884 theProjectileRemnant->setSpin(DeJongSpin::shoot(theProjectileRemnant->getNumberStoredComponents(), theProjectileRemnant->getA()));
00885
00886
00887 theProjectileRemnant->setEmissionTime(anEmissionTime);
00888
00889
00890 theStore->addToOutgoing(theProjectileRemnant);
00891
00892
00893 theProjectileRemnant = NULL;
00894 } else if(theProjectileRemnant->getA()==1) {
00895
00896 Particle *theNucleon = theProjectileRemnant->getParticles().front();
00897 theStore->addToOutgoing(theNucleon);
00898
00899 deleteProjectileRemnant();
00900 } else
00901 deleteProjectileRemnant();
00902 }
00903
00904 }
00905
00906 #endif