#include <G4INCLInteractionAvatar.hh>

Inheritance diagram for G4INCL::InteractionAvatar:

Public Member Functions
	InteractionAvatar (G4double, G4INCL::Nucleus , G4INCL::Particle )

	InteractionAvatar (G4double, G4INCL::Nucleus , G4INCL::Particle , G4INCL::Particle *)

virtual	~InteractionAvatar ()

Public Member Functions inherited from G4INCL::IAvatar
	IAvatar ()

	IAvatar (G4double time)

virtual	~IAvatar ()

FinalState *	getFinalState ()

void	fillFinalState (FinalState *fs)

G4double	getTime () const

virtual ParticleList	getParticles () const =0

virtual std::string	dump () const =0

AvatarType	getType () const

G4bool	isACollision () const

G4bool	isADecay () const

void	setType (AvatarType t)

long	getID () const

std::string	toString ()

Static Public Member Functions
static void	deleteBackupParticles ()
	Release the memory allocated for the backup particles.

Static Public Attributes
static const G4double	locEAccuracy = 1.E-4
	Target accuracy in the determination of the local-energy Q-value.

static const G4int	maxIterLocE = 50
	Max number of iterations for the determination of the local-energy Q-value.

Protected Member Functions
virtual G4INCL::IChannel *	getChannel ()=0

G4bool	bringParticleInside (Particle *const p)

void	preInteractionLocalEnergy (Particle *const p)
	Apply local-energy transformation, if appropriate.

void	preInteractionBlocking ()
	Store the state of the particles before the interaction.

void	preInteraction ()

void	postInteraction (FinalState *)

void	restoreParticles () const
	Restore the state of both particles.

G4bool	shouldUseLocalEnergy () const
	true if the given avatar should use local energy

G4bool	enforceEnergyConservation (FinalState *const fs)
	Enforce energy conservation.

Protected Attributes
Nucleus *	theNucleus

Particle *	particle1

Particle *	particle2

ThreeVector	boostVector

G4double	oldTotalEnergy

G4double	oldXSec

G4bool	isPiN

G4double	weight

ParticleList	modified

ParticleList	created

ParticleList	modifiedAndCreated

ParticleList	Destroyed

ParticleList	ModifiedAndDestroyed

Protected Attributes inherited from G4INCL::IAvatar
G4double	theTime

Static Protected Attributes
static G4ThreadLocal Particle *	backupParticle1 = NULL

static G4ThreadLocal Particle *	backupParticle2 = NULL

Detailed Description

Definition at line 61 of file G4INCLInteractionAvatar.hh.

Constructor & Destructor Documentation

◆ InteractionAvatar() [1/2]

G4INCL::InteractionAvatar::InteractionAvatar	(	G4double	time,
		G4INCL::Nucleus *	n,
		G4INCL::Particle *	p1 )

Definition at line 66 of file G4INCLInteractionAvatar.cc.

    : IAvatar(time), theNucleus(n),
    particle1(p1), particle2(NULL),
    isPiN(false),
    weight(1.),
    violationEFunctor(NULL)
  {
  }

◆ InteractionAvatar() [2/2]

G4INCL::InteractionAvatar::InteractionAvatar	(	G4double	time,
		G4INCL::Nucleus *	n,
		G4INCL::Particle *	p1,
		G4INCL::Particle *	p2 )

Definition at line 75 of file G4INCLInteractionAvatar.cc.

    : IAvatar(time), theNucleus(n),
    particle1(p1), particle2(p2),
    isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
    weight(1.),
    violationEFunctor(NULL)
  {
  }

◆ ~InteractionAvatar()

G4INCL::InteractionAvatar::~InteractionAvatar ( )

virtual

Definition at line 85 of file G4INCLInteractionAvatar.cc.

85 {

86 }

Member Function Documentation

◆ bringParticleInside()

G4bool G4INCL::InteractionAvatar::bringParticleInside ( Particle *const p )

protected

Definition at line 138 of file G4INCLInteractionAvatar.cc.

                                                                  {
    if(!theNucleus)
      return false;
 
    ThreeVector pos = p->getPosition();
    p->rpCorrelate();
    G4double pos2 = pos.mag2();
    const G4double r = theNucleus->getSurfaceRadius(p);
    short iterations=0;
    const short maxIterations=50;
 
    if(pos2 < r*r) return true;
 
    while( pos2 >= r*r && iterations<maxIterations ) /* Loop checking, 10.07.2015, D.Mancusi */
    {
      pos *= std::sqrt(r*r*0.9801/pos2); // 0.9801 == 0.99*0.99
      pos2 = pos.mag2();
      iterations++;
    }
    if( iterations < maxIterations)
    {
      INCL_DEBUG("Particle position vector length was : " << p->getPosition().mag() << ", rescaled to: " << pos.mag() << '\n');
      p->setPosition(pos);
      return true;
    }
    else
      return false;
  }

Referenced by postInteraction().

◆ deleteBackupParticles()

void G4INCL::InteractionAvatar::deleteBackupParticles ( )

static

Release the memory allocated for the backup particles.

Definition at line 88 of file G4INCLInteractionAvatar.cc.

                                                {
    delete backupParticle1;
    if(backupParticle2)
      delete backupParticle2;
    backupParticle1 = NULL;
    backupParticle2 = NULL;
  }

Referenced by G4INCL::INCL::~INCL().

◆ enforceEnergyConservation()

G4bool G4INCL::InteractionAvatar::enforceEnergyConservation ( FinalState *const fs )

protected

Enforce energy conservation.

Final states generated by the channels might violate energy conservation because of different reasons (energy-dependent potentials, local energy...). This conservation law must therefore be enforced by hand. We do so by rescaling the momenta of the final-state particles in the CM frame. If this turns out to be impossible, this method returns false.

Returns: true if the algorithm succeeded

Definition at line 364 of file G4INCLInteractionAvatar.cc.

                                                                           {
    // Set up the violationE calculation
    const G4bool manyBodyFinalState = (modifiedAndCreated.size() > 1);
    
    if(manyBodyFinalState)
      violationEFunctor = new ViolationEMomentumFunctor(theNucleus, modifiedAndCreated, fs->getTotalEnergyBeforeInteraction(), boostVector, shouldUseLocalEnergy());
    else {
      if (modified.empty()) {
        Particle * const p1 = created.front(); //we destroy all nucleons during annihilation in NNbar case
         // The following condition is necessary for the functor to work
         // correctly. A similar condition exists in INCL4.6.
        if(p1->getMass() < ParticleTable::minDeltaMass)
          return false;
        violationEFunctor = new ViolationEEnergyFunctor(theNucleus, p1, fs->getTotalEnergyBeforeInteraction(), shouldUseLocalEnergy());
      }
      else{
        Particle * const p2 = modified.front(); // normal situation
         // The following condition is necessary for the functor to work
         // correctly. A similar condition exists in INCL4.6.
         if(p2->getMass() < ParticleTable::minDeltaMass)
           return false;
         violationEFunctor = new ViolationEEnergyFunctor(theNucleus, p2, fs->getTotalEnergyBeforeInteraction(), shouldUseLocalEnergy());
      }
    }
 
    // Apply the root-finding algorithm
    const RootFinder::Solution theSolution = RootFinder::solve(violationEFunctor, 1.0);
    if(theSolution.success) { // Apply the solution
      (*violationEFunctor)(theSolution.x);
    } else if(theNucleus){
      INCL_DEBUG("Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." << '\n');
      theNucleus->getStore()->getBook().incrementEnergyViolationInteraction();
    }
    delete violationEFunctor;
    violationEFunctor = NULL;
    return theSolution.success;
  }

Referenced by G4INCL::DecayAvatar::postInteraction(), and postInteraction().

◆ getChannel()

virtual G4INCL::IChannel * G4INCL::InteractionAvatar::getChannel ( )

protectedpure virtual

Implements G4INCL::IAvatar.

Implemented in G4INCL::BinaryCollisionAvatar, and G4INCL::DecayAvatar.

◆ postInteraction()

void G4INCL::InteractionAvatar::postInteraction ( FinalState * fs )

protectedvirtual

Implements G4INCL::IAvatar.

Definition at line 167 of file G4INCLInteractionAvatar.cc.

                                                        {
    INCL_DEBUG("postInteraction: final state: " << '\n' << fs->print() << '\n');
    modified = fs->getModifiedParticles();
    created = fs->getCreatedParticles();
    Destroyed = fs->getDestroyedParticles();
    modifiedAndCreated = modified;
    modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end());
    ModifiedAndDestroyed = modified;
    ModifiedAndDestroyed.insert(ModifiedAndDestroyed.end(), Destroyed.begin(), Destroyed.end());
 
    // Boost back to lab
    modifiedAndCreated.boost(-boostVector);
 
    // If there is no Nucleus, just return
    if(!theNucleus) return;
 
    // Mark pions and kaons that have been created outside their well (we will force them
    // to be emitted later).
    for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
      if(((*i)->isPion() || (*i)->isKaon() || (*i)->isAntiKaon()) && (*i)->getPosition().mag() > theNucleus->getSurfaceRadius(*i)) {
        (*i)->makeParticipant();
        (*i)->setOutOfWell();
        fs->addOutgoingParticle(*i);
        INCL_DEBUG("Pion was created outside its potential well." << '\n'
            << (*i)->print());
      }
 
    // Try to enforce energy conservation
    fs->setTotalEnergyBeforeInteraction(oldTotalEnergy);
    G4bool success = enforceEnergyConservation(fs);
    if(!success) {
      INCL_DEBUG("Enforcing energy conservation: failed!" << '\n');
 
      // Restore the state of the initial particles
      restoreParticles();
 
      // Delete newly created particles
      for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
        delete *i;
 
      fs->reset();
      fs->makeNoEnergyConservation();
      fs->setTotalEnergyBeforeInteraction(0.0);
 
      return; // Interaction is blocked. Return an empty final state.
    }
    INCL_DEBUG("Enforcing energy conservation: success!" << '\n');
 
    INCL_DEBUG("postInteraction after energy conservation: final state: " << '\n' << fs->print() << '\n');
 
    // Check that outgoing delta resonances can decay to pi-N
    for(ParticleIter i=modified.begin(), e=modified.end(); i!=e; ++i )
      if((*i)->isDelta() &&
          (*i)->getMass() < ParticleTable::minDeltaMass) {
        INCL_DEBUG("Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
            (*i)->getMass() << '\n');
 
        // Restore the state of the initial particles
        restoreParticles();
 
        // Delete newly created particles
        for(ParticleIter j=created.begin(), end=created.end(); j!=end; ++j )
          delete *j;
 
        fs->reset();
        fs->makeNoEnergyConservation();
        fs->setTotalEnergyBeforeInteraction(0.0);
 
        return; // Interaction is blocked. Return an empty final state.
      }
 
    INCL_DEBUG("Random seeds before Pauli blocking: " << Random::getSeeds() << '\n');
    // Test Pauli blocking
    G4bool isBlocked = Pauli::isBlocked(modifiedAndCreated, theNucleus);
 
    if(isBlocked) {
      INCL_DEBUG("Pauli: Blocked!" << '\n');
 
      // Restore the state of the initial particles
      restoreParticles();
 
      // Delete newly created particles
      for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
        delete *i;
 
      fs->reset();
      fs->makePauliBlocked();
      fs->setTotalEnergyBeforeInteraction(0.0);
 
      return; // Interaction is blocked. Return an empty final state.
    }
    INCL_DEBUG("Pauli: Allowed!" << '\n');
 
    // Test CDPP blocking
    G4bool isCDPPBlocked = Pauli::isCDPPBlocked(created, theNucleus);
 
    if(isCDPPBlocked) {
      INCL_DEBUG("CDPP: Blocked!" << '\n');
 
      // Restore the state of the initial particles
      restoreParticles();
 
      // Delete newly created particles
      for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
        delete *i;
 
      fs->reset();
      fs->makePauliBlocked();
      fs->setTotalEnergyBeforeInteraction(0.0);
 
      return; // Interaction is blocked. Return an empty final state.
    }
    INCL_DEBUG("CDPP: Allowed!" << '\n');
 
    // If all went well, try to bring particles inside the nucleus...
    for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i )
    {
      // ...except for pions beyond their surface radius.
      if((*i)->isOutOfWell()) continue;
 
      const G4bool successBringParticlesInside = bringParticleInside(*i);
      if( !successBringParticlesInside ) {
        INCL_ERROR("Failed to bring particle inside the nucleus!" << '\n');
      }
    }
 
    // Collision accepted!
    // Biasing of the final state
    std::vector<G4int> newBiasCollisionVector;
    newBiasCollisionVector = ModifiedAndDestroyed.getParticleListBiasVector();
    if(std::fabs(weight-1.) > 1E-6){
      newBiasCollisionVector.push_back(Particle::nextBiasedCollisionID);
      Particle::FillINCLBiasVector(1./weight);
      weight = 1.; // useless?
    }
    for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i ) {
      (*i)->setBiasCollisionVector(newBiasCollisionVector);
      if(!(*i)->isOutOfWell()) {
        // Decide if the particle should be made into a spectator
        // (Back to spectator)
        G4bool goesBackToSpectator = false;
        if((*i)->isNucleon() && theNucleus->getStore()->getConfig()->getBackToSpectator()) {
          G4double threshold = (*i)->getPotentialEnergy();
          if((*i)->getType()==Proton)
            threshold += Math::twoThirds*theNucleus->getTransmissionBarrier(*i);
          if((*i)->getKineticEnergy() < threshold)
            goesBackToSpectator = true;
        }
        // Thaw the particle propagation
        (*i)->thawPropagation();
 
        // Increment or decrement the participant counters
        if(goesBackToSpectator) {
          INCL_DEBUG("The following particle goes back to spectator:" << '\n'
              << (*i)->print() << '\n');
          if(!(*i)->isTargetSpectator()) {
            theNucleus->getStore()->getBook().decrementCascading();
          }
          (*i)->makeTargetSpectator();
        } else {
          if((*i)->isTargetSpectator()) {
            theNucleus->getStore()->getBook().incrementCascading();
          }
          (*i)->makeParticipant();
        }
      }
    }
    ParticleList destroyed = fs->getDestroyedParticles();
    for(ParticleIter i=destroyed.begin(), e=destroyed.end(); i!=e; ++i )
      if(!(*i)->isTargetSpectator())
        theNucleus->getStore()->getBook().decrementCascading();
    return;
  }

Referenced by G4INCL::BinaryCollisionAvatar::postInteraction(), and G4INCL::DecayAvatar::postInteraction().

◆ preInteraction()

void G4INCL::InteractionAvatar::preInteraction ( )

protectedvirtual

Implements G4INCL::IAvatar.

Definition at line 123 of file G4INCLInteractionAvatar.cc.

                                         {
    preInteractionBlocking();
 
    preInteractionLocalEnergy(particle1);
 
    if(particle2) {
      preInteractionLocalEnergy(particle2);
      boostVector = KinematicsUtils::makeBoostVector(particle1, particle2);
      particle2->boost(boostVector);
    } else {
      boostVector = particle1->getMomentum()/particle1->getEnergy();
    }
    particle1->boost(boostVector);
  }

Referenced by G4INCL::BinaryCollisionAvatar::preInteraction(), and G4INCL::DecayAvatar::preInteraction().

◆ preInteractionBlocking()

void G4INCL::InteractionAvatar::preInteractionBlocking ( )

protected

Store the state of the particles before the interaction.

If the interaction cannot be realised for any reason, we will need to restore the particle state as it was before. This is done by calling the restoreParticles() method.

Definition at line 96 of file G4INCLInteractionAvatar.cc.

                                                 {
    if(backupParticle1)
      (*backupParticle1) = (*particle1);
    else
      backupParticle1 = new Particle(*particle1);
 
    if(particle2) {
      if(backupParticle2)
        (*backupParticle2) = (*particle2);
      else
        backupParticle2 = new Particle(*particle2);
 
      oldTotalEnergy = particle1->getEnergy() + particle2->getEnergy()
        - particle1->getPotentialEnergy() - particle2->getPotentialEnergy();
      oldXSec = CrossSections::total(particle1, particle2);
    } else {
      oldTotalEnergy = particle1->getEnergy() - particle1->getPotentialEnergy();
    }
  }

Referenced by preInteraction().

◆ preInteractionLocalEnergy()

void G4INCL::InteractionAvatar::preInteractionLocalEnergy ( Particle *const p )

protected

Apply local-energy transformation, if appropriate.

Parameters

p	particle to apply the transformation to

Definition at line 116 of file G4INCLInteractionAvatar.cc.

                                                                      {
    if(!theNucleus || p->isMeson() || p->isPhoton() || p->isAntiNucleon()) return; // Local energy does not make any sense without a nucleus
 
    if(shouldUseLocalEnergy())
      KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p);
  }

Referenced by preInteraction().

◆ restoreParticles()

void G4INCL::InteractionAvatar::restoreParticles ( ) const

protected

Restore the state of both particles.

The state must first be stored by calling preInteractionBlocking().

Definition at line 341 of file G4INCLInteractionAvatar.cc.

                                                 {
    (*particle1) = (*backupParticle1);
    if(particle2)
      (*particle2) = (*backupParticle2);
  }

Referenced by G4INCL::BinaryCollisionAvatar::getChannel(), G4INCL::DecayAvatar::postInteraction(), and postInteraction().

◆ shouldUseLocalEnergy()

G4bool G4INCL::InteractionAvatar::shouldUseLocalEnergy ( ) const

protected

true if the given avatar should use local energy

Definition at line 347 of file G4INCLInteractionAvatar.cc.

                                                       {
    if(!theNucleus) return false;
    LocalEnergyType theLocalEnergyType;
    if(theNucleus->getStore()->getConfig()->getProjectileType()==antiProton ||
     theNucleus->getStore()->getConfig()->getProjectileType()==antiNeutron){
      return false;
    }
    if(getType()==DecayAvatarType || isPiN)
      theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyPiType();
    else
      theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyBBType();
 
    const G4bool firstAvatar = (theNucleus->getStore()->getBook().getAcceptedCollisions() == 0);
    return ((theLocalEnergyType == FirstCollisionLocalEnergy && firstAvatar) ||
            theLocalEnergyType == AlwaysLocalEnergy);
  }