34 #define INCLXX_IN_GEANT4_MODE 1
67 particle1(p1), particle2(NULL),
69 violationEFunctor(NULL)
76 particle1(p1), particle2(p2),
77 isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
78 violationEFunctor(NULL)
95 (*backupParticle1) = (*particle1);
101 (*backupParticle2) = (*particle2);
144 const short maxIterations=50;
146 if(pos2 < r*r)
return true;
148 while( pos2 >= r*r && iterations<maxIterations )
150 pos *= std::sqrt(r*r*0.9801/pos2);
154 if( iterations < maxIterations)
165 INCL_DEBUG(
"postInteraction: final state: " <<
'\n' << fs->
print() <<
'\n');
181 (*i)->makeParticipant();
182 (*i)->setOutOfWell();
184 INCL_DEBUG(
"Pion was created outside its potential well." <<
'\n'
192 INCL_DEBUG(
"Enforcing energy conservation: failed!" <<
'\n');
207 INCL_DEBUG(
"Enforcing energy conservation: success!" <<
'\n');
209 INCL_DEBUG(
"postInteraction after energy conservation: final state: " <<
'\n' << fs->
print() <<
'\n');
212 for(
ParticleIter i=modified.begin(), e=modified.end(); i!=e; ++i )
213 if((*i)->isDelta() &&
215 INCL_DEBUG(
"Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
216 (*i)->getMass() <<
'\n');
279 if((*i)->isOutOfWell())
continue;
282 if( !successBringParticlesInside ) {
283 INCL_ERROR(
"Failed to bring particle inside the nucleus!" <<
'\n');
289 if(!(*i)->isOutOfWell()) {
292 G4bool goesBackToSpectator =
false;
294 G4double threshold = (*i)->getPotentialEnergy();
295 if((*i)->getType()==
Proton)
297 if((*i)->getKineticEnergy() < threshold)
298 goesBackToSpectator =
true;
302 (*i)->thawPropagation();
305 if(goesBackToSpectator) {
306 INCL_DEBUG(
"The following particle goes back to spectator:" <<
'\n'
307 << (*i)->print() <<
'\n');
308 if(!(*i)->isTargetSpectator()) {
311 (*i)->makeTargetSpectator();
313 if((*i)->isTargetSpectator()) {
316 (*i)->makeParticipant();
321 for(
ParticleIter i=destroyed.begin(), e=destroyed.end(); i!=e; ++i )
322 if(!(*i)->isTargetSpectator())
329 (*particle1) = (*backupParticle1);
331 (*particle2) = (*backupParticle2);
350 if(manyBodyFinalState)
364 (*violationEFunctor)(theSolution.
x);
366 INCL_DEBUG(
"Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." <<
'\n');
380 finalParticles(modAndCre),
381 initialEnergy(totalEnergyBeforeInteraction),
384 shouldUseLocalEnergy(localE)
395 particleMomenta.clear();
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::vector<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(totalEnergyBeforeInteraction),
454 theParticle(aParticle),
455 theEnergy(theParticle->getEnergy()),
456 theMomentum(theParticle->getMomentum()),
464 setParticleEnergy(alpha);
465 return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
477 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
478 for(
G4int iterLocE=0;
482 G4double particleEnergy = energyThreshold + locE + alpha*(theEnergy-energyThreshold);
483 const G4double theMass2 = std::pow(particleEnergy,2.)-theMomentum.mag2();
486 theMass = std::sqrt(theMass2);
489 particleEnergy = energyThreshold;
491 theParticle->setMass(theMass);
492 theParticle->setEnergy(particleEnergy);
501 deltaLocE = std::abs(locE-locEOld);
508 setParticleEnergy(1.);
RootFunctor * violationEFunctor
ViolationEEnergyFunctor(Nucleus *const nucleus, Particle *const aParticle, const G4double totalEnergyBeforeInteraction, const G4bool localE)
Prepare for calling the () operator and setParticleEnergy.
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
G4bool isBlocked(ParticleList const &p, Nucleus const *const n)
Check Pauli blocking.
static const G4double locEAccuracy
Target accuracy in the determination of the local-energy Q-value.
G4bool isCDPPBlocked(ParticleList const &p, Nucleus const *const n)
Check CDPP blocking.
void cleanUp(const G4bool success) const
Clean up after root finding.
G4double getMass() const
Get the cached particle mass.
void transformToLocalEnergyFrame(Nucleus const *const n, Particle *const p)
G4int getAcceptedCollisions() const
void preInteractionLocalEnergy(Particle *const p)
Apply local-energy transformation, if appropriate.
G4ThreadLocal G4double minDeltaMass2
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 postInteraction(FinalState *)
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 delta production. ...
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)
virtual ~ViolationEMomentumFunctor()
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.
ThreeVector const & boostVector
Pointer to the boost vector.
ParticleList const & getCreatedParticles() const
void boost(const ThreeVector &b) const
Final state of an interaction.
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
Book & getBook()
Return the pointer to the Book object which keeps track of various counters.
G4bool getBackToSpectator() const
Get back-to-spectator.
ParticleList modifiedAndCreated
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()
void preInteractionBlocking()
Store the state of the particles before the interaction.
static G4ThreadLocal Particle * backupParticle1
void setTotalEnergyBeforeInteraction(G4double E)
ViolationEMomentumFunctor(Nucleus *const nucleus, ParticleList const &modAndCre, const G4double totalEnergyBeforeInteraction, ThreeVector const &boost, const G4bool localE)
Prepare for calling the () operator and scaleParticleMomenta.
G4double getTotalEnergyBeforeInteraction() const
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.
G4ThreadLocal G4double minDeltaMass
void makeNoEnergyConservation()
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4double energy(const ThreeVector &p, const G4double m)
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
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.
std::vector< ThreeVector > particleMomenta
CM particle momenta, as determined by the channel.
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.
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
AvatarType getType() const
ParticleList::const_iterator ParticleIter
Static root-finder algorithm.
G4double getLocalEnergy(Nucleus const *const n, Particle *const p)