33 #define INCLXX_IN_GEANT4_MODE 1
66 particle1(p1), particle2(NULL),
68 violationEFunctor(NULL)
75 particle1(p1), particle2(p2),
76 isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
77 violationEFunctor(NULL)
94 (*backupParticle1) = (*particle1);
100 (*backupParticle2) = (*particle2);
146 const short maxIterations=50;
148 if(pos2 < r*r)
return true;
150 while( pos2 >= r*r && iterations<maxIterations )
152 pos *= std::sqrt(r*r*0.99/pos2);
156 if( iterations < maxIterations)
167 INCL_DEBUG(
"postInteraction: final state: " << std::endl << fs->
print() << std::endl);
171 modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end());
175 for(
ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i )
184 for(
ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
186 (*i)->makeParticipant();
187 (*i)->setOutOfWell();
189 INCL_DEBUG(
"Pion was created outside its potential well." << std::endl
199 INCL_DEBUG(
"Enforcing energy conservation: failed!" << std::endl);
205 for(
ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
215 INCL_DEBUG(
"Enforcing energy conservation: success!" << std::endl);
217 INCL_DEBUG(
"postInteraction after energy conservation: final state: " << std::endl << fs->
print() << std::endl);
220 for(
ParticleIter i=modified.begin(), e=modified.end(); i!=e; ++i )
221 if((*i)->isDelta() &&
223 INCL_DEBUG(
"Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
224 (*i)->getMass() << std::endl);
230 for(
ParticleIter j=created.begin(), end=created.end(); j!=end; ++j )
252 for(
ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
274 for(
ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
287 for(
ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i )
290 if((*i)->isOutOfWell())
continue;
293 if( !successBringParticlesInside ) {
294 INCL_ERROR(
"Failed to bring particle inside the nucleus!" << std::endl);
299 for(
ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i ) {
300 if(!(*i)->isOutOfWell()) {
303 G4bool goesBackToSpectator =
false;
305 G4double threshold = (*i)->getPotentialEnergy();
306 if((*i)->getType()==
Proton)
308 if((*i)->getKineticEnergy() < threshold)
309 goesBackToSpectator =
true;
313 (*i)->thawPropagation();
316 if(goesBackToSpectator) {
317 INCL_DEBUG(
"The following particle goes back to spectator:" << std::endl
318 << (*i)->print() << std::endl);
319 if(!(*i)->isTargetSpectator()) {
322 (*i)->makeTargetSpectator();
324 if((*i)->isTargetSpectator()) {
327 (*i)->makeParticipant();
332 for(
ParticleIter i=destroyed.begin(), e=destroyed.end(); i!=e; ++i )
333 if(!(*i)->isTargetSpectator())
340 (*particle1) = (*backupParticle1);
342 (*particle2) = (*backupParticle2);
349 if(manyBodyFinalState)
352 Particle const *
const p = modified.front();
363 (*violationEFunctor)(theSolution.
x);
365 INCL_DEBUG(
"Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." << std::endl);
379 initialEnergy(finalState->getTotalEnergyBeforeInteraction()),
382 shouldUseLocalEnergy(localE)
399 scaleParticleMomenta(alpha);
402 for(
ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i)
403 deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy();
404 deltaE -= initialEnergy;
410 std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
411 for(
ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i, ++iP) {
412 (*i)->setMomentum((*iP)*alpha);
413 (*i)->adjustEnergyFromMomentum();
419 (*i)->setPotentialEnergy(0.);
426 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
427 for(
G4int iterLocE=0;
431 (*i)->setEnergy(energy + locE);
432 (*i)->adjustMomentumFromEnergy();
435 deltaLocE = std::abs(locE-locEOld);
443 scaleParticleMomenta(1.);
452 initialEnergy(finalState->getTotalEnergyBeforeInteraction()),
454 theParticle(finalState->getModifiedParticles().front()),
455 theEnergy(theParticle->getEnergy()),
456 theMomentum(theParticle->getMomentum()),
465 setParticleEnergy(alpha);
466 return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
478 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
479 for(
G4int iterLocE=0;
483 const G4double particleEnergy = energyThreshold + alpha*(theEnergy-energyThreshold);
484 const G4double theMass = std::sqrt(std::pow(particleEnergy,2.)-theMomentum.mag2());
485 theParticle->setMass(theMass);
486 theParticle->setEnergy(particleEnergy + locE);
487 theParticle->adjustMomentumFromEnergy();
496 deltaLocE = std::abs(locE-locEOld);
503 setParticleEnergy(1.);
RootFunctor * violationEFunctor
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
static const G4double locEAccuracy
Target accuracy in the determination of the local-energy Q-value.
void cleanUp(const G4bool success) const
Clean up after root finding.
G4double getMass() const
Get the cached particle mass.
static G4bool isCDPPBlocked(ParticleList const &p, Nucleus const *const n)
Check CDPP blocking.
void transformToLocalEnergyFrame(Nucleus const *const n, Particle *const p)
void preInteractionLocalEnergy(Particle *const p)
Apply local-energy transformation, if appropriate.
G4bool shouldUseLocalEnergy() const
true if the given avatar should use local energy
ParticleList const & getModifiedParticles() const
static const G4int maxIterLocE
Max number of iterations for the determination of the local-energy Q-value.
ParticleList finalParticles
List of final-state particles.
void boost(const ThreeVector &aBoostVector)
Boost the particle using a boost vector.
const G4INCL::ThreeVector & getMomentum() const
Get the momentum vector.
RootFunctor-derived object for enforcing energy conservation in pi-N.
Config const * getConfig()
Get the config object.
G4double getEnergy() const
Get the energy of the particle in MeV.
G4double getSurfaceRadius(Particle const *const particle) const
Get the maximum allowed radius for a given particle.
ThreeVector makeBoostVector(Particle const *const p1, Particle const *const p2)
ViolationEEnergyFunctor(Nucleus *const nucleus, FinalState const *const finalState, const G4bool localE)
Prepare for calling the () operator and scaleParticleMomenta.
G4double mag2() const
Get the square of the length.
G4double getPotentialEnergy() const
Get the particle potential energy.
static G4ThreadLocal Particle * backupParticle2
void rpCorrelate()
Make the particle follow a strict r-p correlation.
InteractionAvatar(G4double, G4INCL::Nucleus *, G4INCL::Particle *)
G4double operator()(const G4double x) const
Compute the energy-conservation violation.
SeedVector getSeeds()
Get the seeds of the current generator.
ParticleList const & getCreatedParticles() const
Final state of an interaction.
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
static G4bool isBlocked(ParticleList const &p, Nucleus const *const n)
Check Pauli blocking.
Book & getBook()
Return the pointer to the Book object which keeps track of various counters.
G4bool getBackToSpectator() const
Get back-to-spectator.
void setParticleEnergy(const G4double energy) const
Set the energy of the particle.
virtual void setPosition(const G4INCL::ThreeVector &position)
void addOutgoingParticle(Particle *p)
void incrementEnergyViolationInteraction()
FinalState * postInteraction(FinalState *)
void preInteractionBlocking()
Store the state of the particles before the interaction.
static G4ThreadLocal Particle * backupParticle1
void setTotalEnergyBeforeInteraction(G4double E)
void incrementCascading()
RootFunctor-derived object for enforcing energy conservation in N-N.
void decrementCascading()
G4double total(Particle const *const p1, Particle const *const p2)
std::string print() const
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
void makeNoEnergyConservation()
ViolationEMomentumFunctor(Nucleus *const nucleus, FinalState const *const finalState, ThreeVector const *const boost, const G4bool localE)
Prepare for calling the () operator and scaleParticleMomenta.
G4double energy(const ThreeVector &p, const G4double m)
std::list< ThreeVector > particleMomenta
CM particle momenta, as determined by the channel.
ThreeVector const * boostVector
Pointer to the boost vector.
ParticleList const & getDestroyedParticles() const
const G4INCL::ThreeVector & getPosition() const
Set the position vector.
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
G4bool enforceEnergyConservation(FinalState *const fs)
Enforce energy conservation.
void restoreParticles() const
Restore the state of both particles.
G4double operator()(const G4double x) const
Compute the energy-conservation violation.
G4double mag() const
Get the length of the vector.
G4ThreadLocal G4double effectiveDeltaDecayThreshold
G4bool bringParticleInside(Particle *const p)
void scaleParticleMomenta(const G4double alpha) const
Scale the momenta of the modified and created particles.
virtual ~InteractionAvatar()
static const G4double alpha
G4bool isPion() const
Is this a pion?
void cleanUp(const G4bool success) const
Clean up after root finding.
static const G4double pos
ParticleList::const_iterator ParticleIter
Static root-finder algorithm.
G4double getLocalEnergy(Nucleus const *const n, Particle *const p)