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
00045
00046
00047
00048 #include "G4INCLInteractionAvatar.hh"
00049 #include "G4INCLKinematicsUtils.hh"
00050 #include "G4INCLCrossSections.hh"
00051 #include "G4INCLPauliBlocking.hh"
00052 #include "G4INCLRootFinder.hh"
00053 #include "G4INCLLogger.hh"
00054 #include "G4INCLConfigEnums.hh"
00055
00056
00057 namespace G4INCL {
00058
00059 const G4double InteractionAvatar::locEAccuracy = 1.E-4;
00060 const G4int InteractionAvatar::maxIterLocE = 50;
00061
00062 InteractionAvatar::InteractionAvatar(G4double time, G4INCL::Nucleus *n, G4INCL::Particle *p1)
00063 : IAvatar(time), theNucleus(n),
00064 particle1(p1), particle2(NULL), isPiN(false)
00065 {
00066 }
00067
00068 InteractionAvatar::InteractionAvatar(G4double time, G4INCL::Nucleus *n, G4INCL::Particle *p1,
00069 G4INCL::Particle *p2)
00070 : IAvatar(time), theNucleus(n),
00071 particle1(p1), particle2(p2),
00072 isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()))
00073 {
00074 }
00075
00076 InteractionAvatar::~InteractionAvatar() {
00077 }
00078
00079 void InteractionAvatar::preInteractionBlocking() {
00080 oldParticle1Type = particle1->getType();
00081 oldParticle1Energy = particle1->getEnergy();
00082 oldParticle1Potential = particle1->getPotentialEnergy();
00083 oldParticle1Momentum = particle1->getMomentum();
00084 oldParticle1Position = particle1->getPosition();
00085 oldParticle1Mass = particle1->getMass();
00086 oldParticle1Helicity = particle1->getHelicity();
00087
00088 if(particle2) {
00089 oldParticle2Type = particle2->getType();
00090 oldParticle2Energy = particle2->getEnergy();
00091 oldParticle2Potential = particle2->getPotentialEnergy();
00092 oldParticle2Momentum = particle2->getMomentum();
00093 oldParticle2Position = particle2->getPosition();
00094 oldParticle2Mass = particle2->getMass();
00095 oldParticle2Helicity = particle2->getHelicity();
00096 oldTotalEnergy = oldParticle1Energy + oldParticle2Energy
00097 - particle1->getPotentialEnergy() - particle2->getPotentialEnergy();
00098 oldXSec = CrossSections::total(particle1, particle2);
00099 } else {
00100 oldTotalEnergy = oldParticle1Energy - particle1->getPotentialEnergy();
00101 }
00102 }
00103
00104 void InteractionAvatar::preInteractionLocalEnergy(Particle * const p) {
00105 if(!theNucleus || p->isPion()) return;
00106
00107 if(shouldUseLocalEnergy())
00108 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p);
00109 }
00110
00111 void InteractionAvatar::preInteraction() {
00112 preInteractionBlocking();
00113
00114 preInteractionLocalEnergy(particle1);
00115
00116 if(particle2) {
00117 preInteractionLocalEnergy(particle2);
00118 if(!isPiN) {
00119 boostVector = KinematicsUtils::makeBoostVector(particle1, particle2);
00120 particle2->boost(boostVector);
00121 }
00122 } else {
00123 boostVector = particle1->getMomentum()/particle1->getEnergy();
00124 }
00125 if(!isPiN)
00126 particle1->boost(boostVector);
00127 }
00128
00129 G4bool InteractionAvatar::bringParticleInside(Particle * const p) {
00130 ThreeVector pos = p->getPosition();
00131 G4double pos2 = pos.mag2();
00132 const G4double r = theNucleus->getSurfaceRadius(p);
00133 short iterations=0;
00134 const short maxIterations=50;
00135
00136 if(pos2 < r*r) return true;
00137
00138 while( pos2 >= r*r && iterations<maxIterations )
00139 {
00140 pos *= std::sqrt(r*r*0.99/pos2);
00141 pos2 = pos.mag2();
00142 iterations++;
00143 }
00144 if( iterations < maxIterations)
00145 {
00146 DEBUG("Particle position vector length was : " << p->getPosition().mag() << ", rescaled to: " << pos.mag() << std::endl);
00147 p->setPosition(pos);
00148 return true;
00149 }
00150 else
00151 return false;
00152 }
00153
00154 FinalState *InteractionAvatar::postInteraction(FinalState *fs) {
00155 ParticleList modified = fs->getModifiedParticles();
00156 ParticleList modifiedAndCreated = modified;
00157 ParticleList created = fs->getCreatedParticles();
00158 modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end());
00159
00160 if(!isPiN) {
00161
00162 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i )
00163 (*i)->boost(-boostVector);
00164 }
00165
00166
00167 if(!theNucleus) return fs;
00168
00169
00170
00171 for( ParticleIter i = created.begin(); i != created.end(); ++i )
00172 if((*i)->isPion() && (*i)->getPosition().mag() > theNucleus->getSurfaceRadius(*i)) {
00173 (*i)->makeParticipant();
00174 (*i)->setOutOfWell();
00175 fs->addOutgoingParticle(*i);
00176 DEBUG("Pion was created outside its potential well." << std::endl
00177 << (*i)->print());
00178 }
00179
00180
00181 fs->setTotalEnergyBeforeInteraction(oldTotalEnergy);
00182 G4bool success = true;
00183 if(!isPiN || shouldUseLocalEnergy())
00184 success = enforceEnergyConservation(fs);
00185 if(!success) {
00186 DEBUG("Enforcing energy conservation: failed!" << std::endl);
00187
00188
00189 restoreParticles();
00190
00191
00192 for( ParticleIter i = created.begin(); i != created.end(); ++i )
00193 delete *i;
00194
00195 FinalState *fsBlocked = new FinalState;
00196 delete fs;
00197 fsBlocked->makeNoEnergyConservation();
00198 fsBlocked->setTotalEnergyBeforeInteraction(0.0);
00199
00200 return fsBlocked;
00201 }
00202 DEBUG("Enforcing energy conservation: success!" << std::endl);
00203
00204
00205 for( ParticleIter i = modified.begin(); i != modified.end(); ++i )
00206 if((*i)->isDelta() &&
00207 (*i)->getMass() < ParticleTable::effectiveDeltaDecayThreshold) {
00208 DEBUG("Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
00209 (*i)->getMass() << std::endl);
00210
00211
00212 restoreParticles();
00213
00214
00215 for( ParticleIter j = created.begin(); j != created.end(); ++j )
00216 delete *j;
00217
00218 FinalState *fsBlocked = new FinalState;
00219 delete fs;
00220 fsBlocked->makeNoEnergyConservation();
00221 fsBlocked->setTotalEnergyBeforeInteraction(0.0);
00222
00223 return fsBlocked;
00224 }
00225
00226
00227 G4bool isBlocked = Pauli::isBlocked(modifiedAndCreated, theNucleus);
00228
00229 if(isBlocked) {
00230 DEBUG("Pauli: Blocked!" << std::endl);
00231
00232
00233 restoreParticles();
00234
00235
00236 for( ParticleIter i = created.begin(); i != created.end(); ++i )
00237 delete *i;
00238
00239 FinalState *fsBlocked = new FinalState;
00240 delete fs;
00241 fsBlocked->makePauliBlocked();
00242 fsBlocked->setTotalEnergyBeforeInteraction(0.0);
00243
00244 return fsBlocked;
00245 }
00246 DEBUG("Pauli: Allowed!" << std::endl);
00247
00248
00249 G4bool isCDPPBlocked = Pauli::isCDPPBlocked(created, theNucleus);
00250
00251 if(isCDPPBlocked) {
00252 DEBUG("CDPP: Blocked!" << std::endl);
00253
00254
00255 restoreParticles();
00256
00257
00258 for( ParticleIter i = created.begin(); i != created.end(); ++i )
00259 delete *i;
00260
00261 FinalState *fsBlocked = new FinalState;
00262 delete fs;
00263 fsBlocked->makePauliBlocked();
00264 fsBlocked->setTotalEnergyBeforeInteraction(0.0);
00265
00266 return fsBlocked;
00267 }
00268 DEBUG("CDPP: Allowed!" << std::endl);
00269
00270
00271 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i )
00272 {
00273
00274 if((*i)->isOutOfWell()) continue;
00275
00276 const G4bool successBringParticlesInside = bringParticleInside(*i);
00277 if( !successBringParticlesInside ) {
00278 ERROR("Failed to bring particle inside the nucleus!" << std::endl);
00279 }
00280 }
00281
00282
00283 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i ) {
00284 if(!(*i)->isOutOfWell()) {
00285
00286
00287 G4bool goesBackToSpectator = false;
00288 if((*i)->isNucleon() && theNucleus->getStore()->getConfig()->getBackToSpectator()) {
00289 G4double threshold = (*i)->getPotentialEnergy();
00290 if((*i)->getType()==Proton)
00291 threshold += Math::twoThirds*theNucleus->getTransmissionBarrier(*i);
00292 if((*i)->getKineticEnergy() < threshold)
00293 goesBackToSpectator = true;
00294 }
00295
00296
00297 (*i)->thawPropagation();
00298
00299
00300 if(goesBackToSpectator) {
00301 DEBUG("The following particle goes back to spectator:" << std::endl
00302 << (*i)->print() << std::endl);
00303 if(!(*i)->isTargetSpectator()) {
00304 theNucleus->getStore()->getBook()->decrementCascading();
00305 }
00306 (*i)->makeTargetSpectator();
00307 } else {
00308 if((*i)->isTargetSpectator()) {
00309 theNucleus->getStore()->getBook()->incrementCascading();
00310 }
00311 (*i)->makeParticipant();
00312 }
00313 }
00314 }
00315 ParticleList destroyed = fs->getDestroyedParticles();
00316 for( ParticleIter i = destroyed.begin(); i != destroyed.end(); ++i )
00317 if(!(*i)->isTargetSpectator())
00318 theNucleus->getStore()->getBook()->decrementCascading();
00319
00320 return fs;
00321 }
00322
00323 void InteractionAvatar::restoreParticles() const {
00324 particle1->setType(oldParticle1Type);
00325 particle1->setEnergy(oldParticle1Energy);
00326 particle1->setPotentialEnergy(oldParticle1Potential);
00327 particle1->setMomentum(oldParticle1Momentum);
00328 particle1->setPosition(oldParticle1Position);
00329 particle1->setMass(oldParticle1Mass);
00330 particle1->setHelicity(oldParticle1Helicity);
00331
00332 if(particle2) {
00333 particle2->setType(oldParticle2Type);
00334 particle2->setEnergy(oldParticle2Energy);
00335 particle2->setPotentialEnergy(oldParticle2Potential);
00336 particle2->setMomentum(oldParticle2Momentum);
00337 particle2->setPosition(oldParticle2Position);
00338 particle2->setMass(oldParticle2Mass);
00339 particle2->setHelicity(oldParticle2Helicity);
00340 }
00341 }
00342
00343 G4bool InteractionAvatar::enforceEnergyConservation(FinalState * const fs) {
00344
00345 ParticleList modified = fs->getModifiedParticles();
00346 const G4bool manyBodyFinalState = (modified.size() + fs->getCreatedParticles().size() > 1);
00347 if(manyBodyFinalState)
00348 violationEFunctor = new ViolationEMomentumFunctor(theNucleus, fs, &boostVector, shouldUseLocalEnergy());
00349 else {
00350 Particle const * const p = modified.front();
00351
00352
00353 if(p->getMass() < ParticleTable::effectiveDeltaDecayThreshold)
00354 return false;
00355 violationEFunctor = new ViolationEEnergyFunctor(theNucleus, fs);
00356 }
00357
00358
00359 const G4bool success = RootFinder::solve(violationEFunctor, 1.0);
00360 if(success) {
00361 std::pair<G4double,G4double> theSolution = RootFinder::getSolution();
00362 (*violationEFunctor)(theSolution.first);
00363 } else {
00364 WARN("Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." << std::endl);
00365 }
00366 delete violationEFunctor;
00367 return success;
00368 }
00369
00370
00371
00372
00373
00374 InteractionAvatar::ViolationEMomentumFunctor::ViolationEMomentumFunctor(Nucleus * const nucleus, FinalState const * const finalState, ThreeVector const * const boost, const G4bool localE) :
00375 RootFunctor(0., 1E6),
00376 initialEnergy(finalState->getTotalEnergyBeforeInteraction()),
00377 theNucleus(nucleus),
00378 boostVector(boost),
00379 shouldUseLocalEnergy(localE)
00380 {
00381
00382 finalParticles = finalState->getModifiedParticles();
00383 ParticleList created = finalState->getCreatedParticles();
00384 finalParticles.splice(finalParticles.end(), created);
00385
00386
00387
00388 particleMomenta.clear();
00389 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i) {
00390 (*i)->boost(*boostVector);
00391 particleMomenta.push_back((*i)->getMomentum());
00392 }
00393 }
00394
00395 G4double InteractionAvatar::ViolationEMomentumFunctor::operator()(const G4double alpha) const {
00396 scaleParticleMomenta(alpha);
00397
00398 G4double deltaE = 0.0;
00399 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i)
00400 deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy();
00401 deltaE -= initialEnergy;
00402 return deltaE;
00403 }
00404
00405 void InteractionAvatar::ViolationEMomentumFunctor::scaleParticleMomenta(const G4double alpha) const {
00406
00407 std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
00408 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i, ++iP) {
00409 (*i)->setMomentum((*iP)*alpha);
00410 (*i)->adjustEnergyFromMomentum();
00411 (*i)->boost(-(*boostVector));
00412 if(theNucleus)
00413 theNucleus->updatePotentialEnergy(*i);
00414 else
00415 (*i)->setPotentialEnergy(0.);
00416
00417 if(shouldUseLocalEnergy && !(*i)->isPion()) {
00418
00419 const G4double energy = (*i)->getEnergy();
00420 G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, *i);
00421 G4double locEOld;
00422 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
00423 for(G4int iterLocE=0;
00424 deltaLocE>InteractionAvatar::locEAccuracy && iterLocE<InteractionAvatar::maxIterLocE;
00425 ++iterLocE) {
00426 locEOld = locE;
00427 (*i)->setEnergy(energy + locE);
00428 (*i)->adjustMomentumFromEnergy();
00429 theNucleus->updatePotentialEnergy(*i);
00430 locE = KinematicsUtils::getLocalEnergy(theNucleus, *i);
00431 deltaLocE = std::abs(locE-locEOld);
00432 }
00433 }
00434 }
00435 }
00436
00437 void InteractionAvatar::ViolationEMomentumFunctor::cleanUp(const G4bool success) const {
00438 if(!success)
00439 scaleParticleMomenta(1.);
00440 }
00441
00442
00443
00444
00445
00446 InteractionAvatar::ViolationEEnergyFunctor::ViolationEEnergyFunctor(Nucleus * const nucleus, FinalState const * const finalState) :
00447 RootFunctor(0., 1E6),
00448 initialEnergy(finalState->getTotalEnergyBeforeInteraction()),
00449 theNucleus(nucleus),
00450 theParticle(finalState->getModifiedParticles().front()),
00451 theEnergy(theParticle->getEnergy()),
00452 theMomentum(theParticle->getMomentum()),
00453 energyThreshold(KinematicsUtils::energy(theMomentum,ParticleTable::effectiveDeltaDecayThreshold))
00454 {
00455
00456
00457
00458 }
00459
00460 G4double InteractionAvatar::ViolationEEnergyFunctor::operator()(const G4double alpha) const {
00461 setParticleEnergy(alpha);
00462 return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
00463 }
00464
00465 void InteractionAvatar::ViolationEEnergyFunctor::setParticleEnergy(const G4double alpha) const {
00466
00467 G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle);
00468 G4double locEOld;
00469 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
00470 for(G4int iterLocE=0;
00471 deltaLocE>InteractionAvatar::locEAccuracy && iterLocE<InteractionAvatar::maxIterLocE;
00472 ++iterLocE) {
00473 locEOld = locE;
00474 const G4double particleEnergy = energyThreshold + alpha*(theEnergy-energyThreshold);
00475 const G4double theMass = std::sqrt(std::pow(particleEnergy,2.)-theMomentum.mag2());
00476 theParticle->setMass(theMass);
00477 theParticle->setEnergy(particleEnergy + locE);
00478 theParticle->adjustMomentumFromEnergy();
00479 theNucleus->updatePotentialEnergy(theParticle);
00480 locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle);
00481 deltaLocE = std::abs(locE-locEOld);
00482 }
00483
00484 }
00485
00486 void InteractionAvatar::ViolationEEnergyFunctor::cleanUp(const G4bool success) const {
00487 if(!success)
00488 setParticleEnergy(1.);
00489 }
00490
00491 }