G4INCL::StandardPropagationModel Class Reference

#include <G4INCLStandardPropagationModel.hh>

Inheritance diagram for G4INCL::StandardPropagationModel:

Public Member Functions
	StandardPropagationModel (LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType, const G4double hTime=0.0)
virtual	~StandardPropagationModel ()
G4double	getCurrentTime ()
void	setNucleus (G4INCL::Nucleus *nucleus)
G4INCL::Nucleus *	getNucleus ()
G4double	shoot (ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
G4double	shootParticle (ParticleType const t, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
G4double	shootComposite (ParticleSpecies const &s, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
G4double	shootAtrest (ParticleType const t, const G4double kineticEnergy)
G4double	shootCompositeAtrest (ParticleSpecies const &s, const G4double kineticEnergy)
void	setStoppingTime (G4double)
G4double	getStoppingTime ()
void	registerAvatar (G4INCL::IAvatar *anAvatar)
IAvatar *	generateBinaryCollisionAvatar (Particle *const p1, Particle *const p2)
	Generate a two-particle avatar.
G4double	getReflectionTime (G4INCL::Particle const *const aParticle)
	Get the reflection time.
G4double	getTime (G4INCL::Particle const *const particleA, G4INCL::Particle const *const particleB, G4double *minDistOfApproach) const
void	generateUpdatedCollisions (const ParticleList &updatedParticles, const ParticleList &particles)
	Generate and register collisions between a list of updated particles and all the other particles.
void	generateCollisions (const ParticleList &particles)
	Generate and register collisions among particles in a list, except between those in another list.
void	generateCollisions (const ParticleList &particles, const ParticleList &except)
	Generate and register collisions among particles in a list, except between those in another list.
void	generateDecays (const ParticleList &particles)
	Generate decays for particles that can decay.
void	updateAvatars (const ParticleList &particles)
void	generateAllAvatars ()
	(Re)Generate all possible avatars.
G4INCL::IAvatar *	propagate (FinalState const *const fs)
Public Member Functions inherited from G4INCL::IPropagationModel
	IPropagationModel ()
virtual	~IPropagationModel ()

Detailed Description

Standard INCL4 particle propagation and avatar prediction

This class implements the standard INCL4 avatar prediction and particle propagation logic. The main idea is to predict all collisions between particles and their reflections from the potential wall. After this we select the avatar with the smallest time, propagate all particles to their positions at that time and return the avatar to the INCL kernel

See also

G4INCL::Kernel.

The particle trajectories in this propagation model are straight lines and all particles are assumed to move with constant velocity.

Definition at line 69 of file G4INCLStandardPropagationModel.hh.

Constructor & Destructor Documentation

StandardPropagationModel()

G4INCL::StandardPropagationModel::StandardPropagationModel	(	LocalEnergyType	localEnergyType,
		LocalEnergyType	localEnergyDeltaType,
		const G4double	hTime = 0.0 )

Definition at line 69 of file G4INCLStandardPropagationModel.cc.

      :theNucleus(0), maximumTime(70.0), currentTime(0.0),
      hadronizationTime(hTime),
      firstAvatar(true),
      theLocalEnergyType(localEnergyType),
      theLocalEnergyDeltaType(localEnergyDeltaType)
    {
    }

~StandardPropagationModel()

G4INCL::StandardPropagationModel::~StandardPropagationModel ( )

virtual

Definition at line 78 of file G4INCLStandardPropagationModel.cc.

    {
      delete theNucleus;
    }

Member Function Documentation

generateAllAvatars()

void G4INCL::StandardPropagationModel::generateAllAvatars

(

)

(Re)Generate all possible avatars.

Definition at line 878 of file G4INCLStandardPropagationModel.cc.

                                                      {
      ParticleList const &particles = theNucleus->getStore()->getParticles();
// assert(!particles.empty());
      G4double time;
      for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
        time = this->getReflectionTime(*i);
        if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
      }
      generateCollisions(particles);
      generateDecays(particles);
    }

Referenced by shootAtrest(), shootComposite(), shootCompositeAtrest(), and shootParticle().

generateBinaryCollisionAvatar()

IAvatar * G4INCL::StandardPropagationModel::generateBinaryCollisionAvatar	(	Particle *const	p1,
		Particle *const	p2 )

Generate a two-particle avatar.

Generate a two-particle avatar, if all the appropriate conditions are met.

Definition at line 688 of file G4INCLStandardPropagationModel.cc.

                                                                                                             {
      // Is either particle a participant?
      if(!p1->isParticipant() && !p2->isParticipant() && p1->getParticipantType()==p2->getParticipantType()) return NULL;
 
      // Is it a pi-resonance collision (we don't treat them)?
      if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
        return NULL;
 
      // Is it a photon collision (we don't treat them)?
      if(p1->isPhoton() || p2->isPhoton())
        return NULL;
 
      // Will the avatar take place between now and the end of the cascade?
      G4double minDistOfApproachSquared = 0.0;
      G4double t = getTime(p1, p2, &minDistOfApproachSquared);
      if(t>maximumTime || t<currentTime+hadronizationTime) return NULL;
 
      // Local energy. Jump through some hoops to calculate the cross section
      // at the collision point, and clean up after yourself afterwards.
      G4bool hasLocalEnergy;
      if(p1->isPion() || p2->isPion())
        hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
              theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
            theLocalEnergyDeltaType == AlwaysLocalEnergy);
      else
        hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
              theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
            theLocalEnergyType == AlwaysLocalEnergy);
      const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isMeson() && !p1->isAntiNucleon());
      const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isMeson() && !p2->isAntiNucleon());
 
      if(p1HasLocalEnergy) {
        backupParticle1 = *p1;
        p1->propagate(t - currentTime);
        if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
          *p1 = backupParticle1;
          return NULL;
        }
        KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p1);
      } 
      if(p2HasLocalEnergy) {
        backupParticle2 = *p2;
        p2->propagate(t - currentTime);
        if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
          *p2 = backupParticle2;
          if(p1HasLocalEnergy) {
            *p1 = backupParticle1;
          }
          return NULL;
        }
        KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p2);
      }
 
      // Compute the total cross section
      const G4double totalCrossSection = CrossSections::total(p1, p2);
      const G4double squareTotalEnergyInCM = KinematicsUtils::squareTotalEnergyInCM(p1,p2);
 
      // Restore particles to their state before the local-energy tweak
      if(p1HasLocalEnergy) {
        *p1 = backupParticle1;
      }
      if(p2HasLocalEnergy) {
        *p2 = backupParticle2;
      }
 
      // Is the CM energy > cutNN? (no cutNN on the first collision)
      if(theNucleus->getStore()->getBook().getAcceptedCollisions()>0
          && p1->isNucleon() && p2->isNucleon()
          && squareTotalEnergyInCM < BinaryCollisionAvatar::getCutNNSquared()) return NULL;
 
      // Do the particles come close enough to each other?
      if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
 
      // Bomb out if the two collision partners are the same particle
// assert(p1->getID() != p2->getID());
 
      // Return a new avatar, then!
      return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
    }

Referenced by generateCollisions(), generateCollisions(), and generateUpdatedCollisions().

generateCollisions() [1/2]

void G4INCL::StandardPropagationModel::generateCollisions

(

const ParticleList &

particles

)

Generate and register collisions among particles in a list, except between those in another list.

This method generates all possible collisions among the particles. Each collision is generated only once.

Parameters

particles

list of particles

Definition at line 838 of file G4INCLStandardPropagationModel.cc.

                                                                                   {
      // Loop over all the particles
      for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1) {
        // Loop over the rest of the particles
        for(ParticleIter p2 = p1 + 1; p2 != particles.end(); ++p2) {
          registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
        }
      }
    }

Referenced by generateAllAvatars().

generateCollisions() [2/2]

void G4INCL::StandardPropagationModel::generateCollisions	(	const ParticleList &	particles,
		const ParticleList &	except )

Generate and register collisions among particles in a list, except between those in another list.

This method generates all possible collisions among the particles. Each collision is generated only once. The collision is NOT generated if BOTH collision partners belong to the except list.

You should pass an empty list as the except parameter if you want to generate all possible collisions among particles.

Parameters

particles	list of particles
except	list of excluded particles

Definition at line 848 of file G4INCLStandardPropagationModel.cc.

                                                                                                               {
 
      const G4bool haveExcept = !except.empty();
 
      // Loop over all the particles
      for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1)
      {
        // Loop over the rest of the particles
        ParticleIter p2 = p1;
        for(++p2; p2 != particles.end(); ++p2)
        {
          // Skip the collision if both particles must be excluded
          if(haveExcept && except.contains(*p1) && except.contains(*p2)) continue;
 
          registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
        }
      }
 
    }

generateDecays()

void G4INCL::StandardPropagationModel::generateDecays

(

const ParticleList &

particles

)

Generate decays for particles that can decay.

The list of particles given as an argument is allowed to contain also stable particles.

Parameters

particles

list of particles to (possibly) generate decays for

Definition at line 906 of file G4INCLStandardPropagationModel.cc.

                                                                               {
      for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
           if((*i)->isDelta()) {
             G4double decayTime = DeltaDecayChannel::computeDecayTime((*i)); // time in fm/c
             G4double time = currentTime + decayTime;
             if(time <= maximumTime) {
               registerAvatar(new DecayAvatar((*i), time, theNucleus));
             }
           }
           else if((*i)->getType() == SigmaZero) {
             G4double decayTime = SigmaZeroDecayChannel::computeDecayTime((*i)); // time in fm/c
             G4double time = currentTime + decayTime;
             if(time <= maximumTime) {
               registerAvatar(new DecayAvatar((*i), time, theNucleus));
             }
           }
        if((*i)->isOmega()) {
          G4double decayTimeOmega = PionResonanceDecayChannel::computeDecayTime((*i));
          G4double timeOmega = currentTime + decayTimeOmega;
          if(timeOmega <= maximumTime) {
            registerAvatar(new DecayAvatar((*i), timeOmega, theNucleus));
          }
        }
      }
    }

Referenced by generateAllAvatars(), and propagate().

generateUpdatedCollisions()

void G4INCL::StandardPropagationModel::generateUpdatedCollisions	(	const ParticleList &	updatedParticles,
		const ParticleList &	particles )

Generate and register collisions between a list of updated particles and all the other particles.

This method does not generate collisions among the particles in updatedParticles; in other words, it generates a collision between one of the updatedParticles and one of the particles ONLY IF the latter does not belong to updatedParticles.

If you intend to generate all possible collisions among particles in a list, use generateCollisions().

Parameters

updatedParticles	list of updated particles
particles	list of particles

Definition at line 819 of file G4INCLStandardPropagationModel.cc.

                                                                                                                                {
 
      // Loop over all the updated particles
      for(ParticleIter updated=updatedParticles.begin(), e=updatedParticles.end(); updated!=e; ++updated)
      {
        // Loop over all the particles
        for(ParticleIter particle=particles.begin(), end=particles.end(); particle!=end; ++particle)
        {
          /* Consider the generation of a collision avatar only if (*particle)
           * is not one of the updated particles.
           * The criterion makes sure that you don't generate avatars between
           * updated particles. */
          if(updatedParticles.contains(*particle)) continue;
 
          registerAvatar(generateBinaryCollisionAvatar(*particle,*updated));
        }
      }
    }

Referenced by updateAvatars().

getCurrentTime()

G4double G4INCL::StandardPropagationModel::getCurrentTime ( )

virtual

Returns the current global time of the system.

Implements G4INCL::IPropagationModel.

Definition at line 674 of file G4INCLStandardPropagationModel.cc.

                                                      {
      return currentTime;
    }

getNucleus()

G4INCL::Nucleus * G4INCL::StandardPropagationModel::getNucleus ( )

virtual

Get the nucleus.

Implements G4INCL::IPropagationModel.

Definition at line 83 of file G4INCLStandardPropagationModel.cc.

    {
      return theNucleus;
    }

getReflectionTime()

G4double G4INCL::StandardPropagationModel::getReflectionTime

(

G4INCL::Particle const *const

aParticle

)

Get the reflection time.

Returns the reflection time of a particle on the potential wall.

Parameters

aParticle

pointer to the particle

Definition at line 768 of file G4INCLStandardPropagationModel.cc.

                                                                                               {
      Intersection theIntersection(
          IntersectionFactory::getLaterTrajectoryIntersection(
            aParticle->getPosition(),
            aParticle->getPropagationVelocity(),
            theNucleus->getSurfaceRadius(aParticle)));
      G4double time;
      if(theIntersection.exists) {
        time = currentTime + theIntersection.time;
      } else {
        INCL_ERROR("Imaginary reflection time for particle: " << '\n'
          << aParticle->print());
        time = 10000.0;
      }
      return time;
    }

Referenced by generateAllAvatars(), and updateAvatars().

getStoppingTime()

G4double G4INCL::StandardPropagationModel::getStoppingTime ( )

virtual

Get the current stopping time.

Implements G4INCL::IPropagationModel.

Definition at line 665 of file G4INCLStandardPropagationModel.cc.

                                                       {
      return maximumTime;
    }

getTime()

G4double G4INCL::StandardPropagationModel::getTime	(	G4INCL::Particle const *const	particleA,
		G4INCL::Particle const *const	particleB,
		G4double *	minDistOfApproach ) const

Get the predicted time of the collision between two particles.

Definition at line 785 of file G4INCLStandardPropagationModel.cc.

    {
      G4double time;
 
      // When annihilation is forced for antinucleons below the threshold energy, set the time step at 0.001 (when to have the smallest avatar time)
      Config const *theConfig=theNucleus->getStore()->getConfig();
      if (((particleA->getType()==antiProton)  && (particleA->getKineticEnergy() <= theConfig->getAtrestThreshold()))   ||
          ((particleB->getType()==antiProton)  && (particleB->getKineticEnergy() <= theConfig->getAtrestThreshold()))   ||
          ((particleA->getType()==antiNeutron) && (particleA->getKineticEnergy() <= particleA->getPotentialEnergy())) ||
          ((particleB->getType()==antiNeutron) && (particleB->getKineticEnergy() <= particleB->getPotentialEnergy())) ||
          ((particleA->getType()==antiProton  && particleA->getEnergy() <= particleA->getINCLMass()))   ||
          ((particleB->getType()==antiProton  && particleB->getEnergy() <= particleB->getINCLMass()))   ||
          ((particleA->getType()==antiNeutron && particleA->getEnergy() <= particleA->getINCLMass())) ||
          ((particleB->getType()==antiNeutron && particleB->getEnergy() <= particleB->getINCLMass())))
      {
        return currentTime + 0.001;
      }
      G4INCL::ThreeVector t13 = particleA->getPropagationVelocity();
      t13 -= particleB->getPropagationVelocity();
      G4INCL::ThreeVector distance = particleA->getPosition();
      distance -= particleB->getPosition();
      const G4double t7 = t13.dot(distance);
      const G4double dt = t13.mag2();
      if(dt <= 1.0e-10) {
        (*minDistOfApproach) = 100000.0;
        return currentTime + 100000.0;
      } else {
        time = -t7/dt;
      }
      (*minDistOfApproach) = distance.mag2() + time * t7;
      return currentTime + time;
    }

Referenced by generateBinaryCollisionAvatar().

propagate()

G4INCL::IAvatar * G4INCL::StandardPropagationModel::propagate ( FinalState const *const fs )

virtual

Propagate all particles and return the first avatar.

Implements G4INCL::IPropagationModel.

Definition at line 932 of file G4INCLStandardPropagationModel.cc.

    {
      if(fs) {
        // We update only the information related to particles that were updated
        // by the previous avatar.
#ifdef INCL_REGENERATE_AVATARS
#warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
        if(!fs->getModifiedParticles().empty() || !fs->getEnteringParticles().empty() || !fs->getCreatedParticles().empty()) {
          // Regenerates the entire avatar list, skipping collisions between
          // updated particles
          theNucleus->getStore()->clearAvatars();
          theNucleus->getStore()->initialiseParticleAvatarConnections();
          generateAllAvatarsExceptUpdated(fs);
        }
#else
        ParticleList const &updatedParticles = fs->getModifiedParticles();
        if(fs->getValidity()==PauliBlockedFS) {
          // This final state might represents the outcome of a Pauli-blocked delta
          // decay
// assert(updatedParticles.empty() || (updatedParticles.size()==1 && updatedParticles.front()->isResonance()));
// assert(fs->getEnteringParticles().empty() && fs->getCreatedParticles().empty() && fs->getOutgoingParticles().empty() && fs->getDestroyedParticles().empty());
          generateDecays(updatedParticles);
        } else {
          ParticleList const &entering = fs->getEnteringParticles();
          generateDecays(updatedParticles);
          generateDecays(entering);
 
          ParticleList const &created = fs->getCreatedParticles();
          if(created.empty() && entering.empty())
            updateAvatars(updatedParticles);
          else {
            ParticleList updatedParticlesCopy = updatedParticles;
            updatedParticlesCopy.insert(updatedParticlesCopy.end(), entering.begin(), entering.end());
            updatedParticlesCopy.insert(updatedParticlesCopy.end(), created.begin(), created.end());
            updateAvatars(updatedParticlesCopy);
          }
        }
#endif
      }
 
      G4INCL::IAvatar *theAvatar = theNucleus->getStore()->findSmallestTime();
      if(theAvatar == 0) return 0; // Avatar list is empty
      //      theAvatar->dispose();
 
      if(theAvatar->getTime() < currentTime) {
        INCL_ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << '\n');
        return 0;
      } else if(theAvatar->getTime() > currentTime) {
        theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
 
        currentTime = theAvatar->getTime();
        theNucleus->getStore()->getBook().setCurrentTime(currentTime);
      }
 
      return theAvatar;
    }

registerAvatar()

void G4INCL::StandardPropagationModel::registerAvatar

(

G4INCL::IAvatar *

anAvatar

)

Add an avatar to the storage.

Definition at line 683 of file G4INCLStandardPropagationModel.cc.

    {
      if(anAvatar) theNucleus->getStore()->add(anAvatar);
    }

Referenced by generateAllAvatars(), generateCollisions(), generateCollisions(), generateDecays(), generateUpdatedCollisions(), and updateAvatars().

setNucleus()

void G4INCL::StandardPropagationModel::setNucleus ( G4INCL::Nucleus * nucleus )

virtual

Set the nucleus for this propagation model.

Implements G4INCL::IPropagationModel.

Definition at line 678 of file G4INCLStandardPropagationModel.cc.

    {
      theNucleus = nucleus;
    }

setStoppingTime()

void G4INCL::StandardPropagationModel::setStoppingTime ( G4double time )

virtual

Set the stopping time of the simulation.

Implements G4INCL::IPropagationModel.

Definition at line 669 of file G4INCLStandardPropagationModel.cc.

                                                                {
// assert(time>0.0);
      maximumTime = time;
    }

shoot()

G4double G4INCL::StandardPropagationModel::shoot ( ParticleSpecies const & projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi )

virtual

Implements G4INCL::IPropagationModel.

Definition at line 90 of file G4INCLStandardPropagationModel.cc.

                                                                                                                                                                       {
      if(projectileSpecies.theType==Composite || projectileSpecies.theType==antiComposite){
        if(theNucleus->getAnnihilationType()!=Def)
          return shootCompositeAtrest(projectileSpecies,kineticEnergy);
        else
          return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
      }
      else if((projectileSpecies.theType==antiProton || projectileSpecies.theType==antiNeutron) && theNucleus->getAnnihilationType()!=Def){
        return shootAtrest(projectileSpecies.theType, kineticEnergy);
      }
      else{
        return shootParticle(projectileSpecies.theType, kineticEnergy, impactParameter, phi);
      }
    }

shootAtrest()

G4double G4INCL::StandardPropagationModel::shootAtrest ( ParticleType const t, const G4double kineticEnergy )

virtual

Implements G4INCL::IPropagationModel.

Definition at line 107 of file G4INCLStandardPropagationModel.cc.

                                                                                                     {
      theNucleus->setParticleNucleusCollision(); 
      currentTime = 0.0;
 
      // Create final state particles
      const G4double projectileMass = ParticleTable::getTableParticleMass(t);
      G4double energy = kineticEnergy + projectileMass;
      G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
      ThreeVector momentum(0.0, 0.0, momentumZ);
      Particle *pb = new G4INCL::Particle(t, energy, momentum, ThreeVector());
      if (t == antiProton){
      PbarAtrestEntryChannel *obj = new PbarAtrestEntryChannel(theNucleus, pb); 
      ParticleList fslist = obj->makeMesonStar();
      const G4bool isProton = obj->ProtonIsTheVictim();
      delete pb;
 
      //set Stopping time according to highest meson energy of the star
      G4double temfin;
      G4double TLab;
      std::vector<G4double> energies;
      std::vector<G4double> projections;
      ThreeVector ab, cd;
 
      for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
        energies.push_back((*pit)->getKineticEnergy());
        ab = (*pit)->boostVector();
        cd = (*pit)->getPosition();
        projections.push_back(ab.dot(cd)); //projection length
      }// make vector of energies
      
      temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
      TLab = *max_element(energies.begin(), energies.end()); //choose max energy
 
      // energy-dependent stopping time above 2 AGeV
      if(TLab>2000.)
        temfin *= (5.8E4-TLab)/5.6E4;
 
      maximumTime = temfin;  
 
      // If the incoming particle is slow, use a larger stopping time
      const G4double rMax = theNucleus->getUniverseRadius();
      const G4double distance = 2.*rMax;
      const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
      const G4double traversalTime = distance / maxMesonVelocityProjection;
      if(maximumTime < traversalTime)
        maximumTime = traversalTime;
      INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
 
      // Fill in the relevant kinematic variables
      theNucleus->setIncomingAngularMomentum(G4INCL::ThreeVector(0., 0., 0.));
      theNucleus->setIncomingMomentum(G4INCL::ThreeVector(0., 0., 0.));
      if(isProton){
        theNucleus->setInitialEnergy(pb->getEnergy()
          + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ() + 1,theNucleus->getS()));
      }
      else{
        theNucleus->setInitialEnergy(pb->getEnergy()
          + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ(),theNucleus->getS()));
      }
      //kinetic energy excluded from the balance
 
      for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
        (*p)->makeProjectileSpectator();                                        
      }
      
      generateAllAvatars();
      firstAvatar = false;
 
      // Get the entry avatars for mesons
      IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
      delete obj;
      theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
      INCL_DEBUG("Avatars added" << '\n');
       
      } //end (t == antiProton)
      else if (t == antiNeutron){
        NbarAtrestEntryChannel *obj = new NbarAtrestEntryChannel(theNucleus, pb);
        ParticleList fslist = obj->makeMesonStar();
        const bool isProton = obj->ProtonIsTheVictim();
        delete pb;
 
        //set Stopping time according to highest meson energy of the star
        G4double temfin;
        G4double TLab;
        std::vector<double> energies;
        std::vector<double> projections;
        ThreeVector ab, cd;
 
        for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
          energies.push_back((*pit)->getKineticEnergy());
          ab = (*pit)->boostVector();
          cd = (*pit)->getPosition();
          projections.push_back(ab.dot(cd)); //projection length
        }// make vector of energies
      
        temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
        TLab = *max_element(energies.begin(), energies.end()); //choose max energy
 
        // energy-dependent stopping time above 2 AGeV
        if(TLab>2000.)
          temfin *= (5.8E4-TLab)/5.6E4;
 
        maximumTime = temfin;  
 
        // If the incoming particle is slow, use a larger stopping time
        const G4double rMax = theNucleus->getUniverseRadius();
        const G4double distance = 2.*rMax;
        const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
        const G4double traversalTime = distance / maxMesonVelocityProjection;
        if(maximumTime < traversalTime)
          maximumTime = traversalTime;
        INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
 
        // Fill in the relevant kinematic variables
        theNucleus->setIncomingAngularMomentum(G4INCL::ThreeVector(0., 0., 0.));
        theNucleus->setIncomingMomentum(G4INCL::ThreeVector(0., 0., 0.));
        if(isProton){
          theNucleus->setInitialEnergy(pb->getEnergy()
            + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ() + 1,theNucleus->getS()));
        }
        else{
          theNucleus->setInitialEnergy(pb->getEnergy()
            + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ(),theNucleus->getS()));
        }
        //kinetic energy excluded from the balance
 
        for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
          (*p)->makeProjectileSpectator();                                        
        }
      
        generateAllAvatars();
        firstAvatar = false;
 
        // Get the entry avatars for mesons
        IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
        delete obj;
        theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
        INCL_DEBUG("Avatars added" << '\n');
      } 
      
      return 99.;
    } 

Referenced by shoot().

shootComposite()

G4double G4INCL::StandardPropagationModel::shootComposite ( ParticleSpecies const & s, const G4double kineticEnergy, const G4double impactParameter, const G4double phi )

virtual

Implements G4INCL::IPropagationModel.

Definition at line 334 of file G4INCLStandardPropagationModel.cc.

                                                                                                                                                                      {
      theNucleus->setNucleusNucleusCollision();
      currentTime = 0.0;
 
      // Create the ProjectileRemnant object
      ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
 
      // Same stopping time as for nucleon-nucleus
      maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
 
      // If the incoming cluster is slow, use a larger stopping 
      G4double rms=0.;
      if(species.theType == Composite){
        rms = ParticleTable::getLargestNuclearRadius(pr->getA(), pr->getZ());
      }
      else if(species.theType == antiComposite){
        rms = ParticleTable::getLargestNuclearRadius(-(pr->getA()), -(pr->getZ()));
      }
      else {
        INCL_ERROR("a non-composite try to go through shootComposite : " << species.theType << '\n');
      }
      const G4double rMax = theNucleus->getUniverseRadius();
      const G4double distance = 2.*rMax + 2.725*rms;
      const G4double projectileVelocity = pr->boostVector().mag();
      const G4double traversalTime = distance / projectileVelocity;
      if(maximumTime < traversalTime)
        maximumTime = traversalTime;
      INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
      // If Coulomb is activated, do not process events with impact
      // parameter larger than the maximum impact parameter, taking into
      // account Coulomb distortion.
      if(impactParameter>CoulombDistortion::maxImpactParameter(pr,theNucleus)) {
        INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
        delete pr;
        return -1.;
      }
 
      // Position the cluster at the right impact parameter
      ThreeVector position(impactParameter * std::cos(phi),
          impactParameter * std::sin(phi),
          0.);
      pr->setPosition(position);
 
      // Fill in the relevant kinematic variables
      theNucleus->setIncomingAngularMomentum(pr->getAngularMomentum());
      theNucleus->setIncomingMomentum(pr->getMomentum());
      theNucleus->setInitialEnergy(pr->getEnergy()
          + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ(),theNucleus->getS()));
 
      generateAllAvatars();
      firstAvatar = false;
 
      // Get the entry avatars from Coulomb
      IAvatarList theAvatarList
        = CoulombDistortion::bringToSurface(pr, theNucleus);
 
      if(theAvatarList.empty()) {
        INCL_DEBUG("No ParticleEntryAvatar found, transparent event" << '\n');
        delete pr;
        return -1.;
      }
 
      /* Store the internal kinematics of the projectile remnant.
       *
       * Note that this is at variance with the Fortran version, which stores
       * the initial kinematics of the particles *after* putting the spectators
       * on mass shell, but *before* removing the same energy from the
       * participants. Due to the different code flow, doing so in the C++
       * version leads to wrong excitation energies for the forced compound
       * nucleus.
       */
      pr->storeComponents();
 
      // Tell the Nucleus about the ProjectileRemnant
      theNucleus->setProjectileRemnant(pr);
 
      // Register the ParticleEntryAvatars
      theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
 
      return pr->getTransversePosition().mag();
    }

Referenced by shoot().

shootCompositeAtrest()

G4double G4INCL::StandardPropagationModel::shootCompositeAtrest ( ParticleSpecies const & s, const G4double kineticEnergy )

virtual

Implements G4INCL::IPropagationModel.

Definition at line 417 of file G4INCLStandardPropagationModel.cc.

                                                                                                                       {
      if(theNucleus->getAnnihilationType()==PType || theNucleus->getAnnihilationType()==NType){
        INCL_DEBUG("Antideuteron annihilation Model B chosen, Annihilation of one antinucleon " << '\n');
        theNucleus->setParticleNucleusCollision();
        currentTime = 0.0;
        
        //Dummy Cluster to intialise the anticomposite and distribute the energy and positio
        Cluster *DummyC = new Cluster(-1,-2,0);
        DummyC->setTableMass();
        DummyC->initializeParticles();
        DummyC->internalBoostToCM();
        const G4double projectileMass = DummyC->getMass();
        const G4double energy = kineticEnergy + projectileMass;
        const G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
        const ThreeVector aBoostVector = ThreeVector(0.0, 0.0, momentumZ / energy);
        DummyC->boost(-aBoostVector);
        DummyC->makeProjectileSpectator();
 
        Particle *pb = new Particle(antiProton, 1, ThreeVector(), ThreeVector());
        Particle *nb= new Particle(antiNeutron, 1, ThreeVector(), ThreeVector());
        ParticleList Antis = DummyC->getParticles();
        for(ParticleIter i = Antis.begin(), e=Antis.end();i!=e;++i){
          if((*i)->getType()==antiProton){
            pb = (*i);
            DummyC->removeParticle((*i));
          }
          else if((*i)->getType()==antiNeutron){
            nb = (*i);
            DummyC->removeParticle((*i));
          }
        }
        delete DummyC;
        Config const *theConfig=theNucleus->getStore()->getConfig();
        if(nb->getKineticEnergy() <= theConfig->getnbAtrestThreshold() && (pb->getEnergy() >= pb->getMass())){
          INCL_DEBUG("Annihilation of the Antineutron " << '\n');
          NbarAtrestEntryChannel *obj = new NbarAtrestEntryChannel(theNucleus,  nb);
          ParticleList fslist = obj->makeMesonStar();
          const G4bool isProton = obj->ProtonIsTheVictim();
          //delete nb;
 
          //set Stopping time according to highest meson energy of the star
          G4double temfin;
          G4double TLab;
          std::vector<G4double> energies;
          std::vector<G4double> projections;
          ThreeVector ab, cd;
          for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
            energies.push_back((*pit)->getKineticEnergy());
            ab = (*pit)->boostVector();
            cd = (*pit)->getPosition();
            projections.push_back(ab.dot(cd)); //projection length
          }// make vector of energies
          temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
          TLab = *max_element(energies.begin(), energies.end()); //choose max energy
          if(TLab>2000.)
            temfin *= (5.8E4-TLab)/5.6E4;
          maximumTime = temfin;  
          // If the incoming particle is slow, use a larger stopping time
          const G4double rMax = theNucleus->getUniverseRadius();
          const G4double distance = 2.*rMax;
          const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
          const G4double traversalTime = distance / maxMesonVelocityProjection;
          if(maximumTime < traversalTime)
            maximumTime = traversalTime;
          INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
          // Fill in the relevant kinematic variables
          theNucleus->setIncomingAngularMomentum(pb->getAngularMomentum());
          theNucleus->setIncomingMomentum(pb->getMomentum());
          if(isProton){
            theNucleus->setInitialEnergy(nb->getMass() + pb->getEnergy()
              + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ() + 1,theNucleus->getS()));
          }
          else{
            theNucleus->setInitialEnergy(nb->getMass() + pb->getEnergy()
              + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ(),theNucleus->getS()));
          }
          // Reset the particle kinematics to the INCL values
          for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
            (*p)->makeProjectileSpectator();                                        
          }
          pb->makeProjectileSpectator();
        
          generateAllAvatars();
          firstAvatar = false;
 
          // Get the entry avatars for mesons
          IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
          delete obj;
          theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
          // Get the entry avatars from Coulomb and put them in the Store
          ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurfaceAbar(pb, theNucleus);
          if(theEntryAvatar) {
            theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
            INCL_DEBUG("Avatars added" << '\n');
            return pb->getTransversePosition().mag();
          } else {
            INCL_DEBUG("Antiproton is transparent, not entering the nucleus " << '\n');
            //Transparent event 
            theNucleus->getStore()->addToMissed(pb);
            delete nb;
            return 99.;
          }
 
        }
        else{
          INCL_DEBUG("Annihilation of the Antiproton " << '\n');
          PbarAtrestEntryChannel *obj = new PbarAtrestEntryChannel(theNucleus, pb);
          ParticleList fslist = obj->makeMesonStar();
          const G4bool isProton = obj->ProtonIsTheVictim();
          //delete pb;
        
          //set Stopping time according to highest meson energy of the star
          G4double temfin;
          G4double TLab;
          std::vector<G4double> energies;
          std::vector<G4double> projections;
          ThreeVector ab, cd;
          for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
            energies.push_back((*pit)->getKineticEnergy());
            ab = (*pit)->boostVector();
            cd = (*pit)->getPosition();
            projections.push_back(ab.dot(cd)); //projection length
          }// make vector of energies
          temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
          TLab = *max_element(energies.begin(), energies.end()); //choose max energy
          if(TLab>2000.)
            temfin *= (5.8E4-TLab)/5.6E4;
          maximumTime = temfin;  
          // If the incoming particle is slow, use a larger stopping time
          const G4double rMax = theNucleus->getUniverseRadius();
          const G4double distance = 2.*rMax;
          const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
          const G4double traversalTime = distance / maxMesonVelocityProjection;
          if(maximumTime < traversalTime)
            maximumTime = traversalTime;
          INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
 
 
          // Fill in the relevant kinematic variables
          theNucleus->setIncomingAngularMomentum(nb->getAngularMomentum());
          theNucleus->setIncomingMomentum(nb->getMomentum());
          if(isProton){
            theNucleus->setInitialEnergy(pb->getMass() + nb->getEnergy()
              + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ() + 1,theNucleus->getS()));
          }
          else{
            theNucleus->setInitialEnergy(pb->getMass() + nb->getEnergy()
              + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ(),theNucleus->getS()));
          }
          // Reset the particle kinematics to the INCL values
          for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
            (*p)->makeProjectileSpectator();                                        
          }
          nb->makeProjectileSpectator();
        
          generateAllAvatars();
          firstAvatar = false;
 
          // Get the entry avatars for mesons
          IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
          delete obj;
          theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
          // Get the entry avatars from Coulomb and put them in the Store
          ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurfaceAbar(nb, theNucleus);
          if(theEntryAvatar) {
            theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
            INCL_DEBUG("Avatars added" << '\n');
            return nb->getTransversePosition().mag();
          } else {
            INCL_DEBUG("Antineutron is transparent, not entering the nucleus " << '\n');
            theNucleus->setIncomingAngularMomentum(ThreeVector(0.,0.,0.));
            theNucleus->setIncomingMomentum(ThreeVector(0.,0.,0.));
            //Transparent event 
            theNucleus->getStore()->addToMissed(nb);
            delete pb;
            return 99.;
          }
          delete theConfig;
        }
      } else{
        theNucleus->setNucleusNucleusCollision();
        currentTime =0.0;
        maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
 
        ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
        INCL_DEBUG("Antideuteron annihilation Model A chosen, Annihilation of Antideuteron as a whole" << '\n');
        AntinucleiAtrestEntryChannel *obj = new AntinucleiAtrestEntryChannel(theNucleus, pr, ThreeVector(), ThreeVector());
        ParticleList fslist = obj->makeMesonStar();
        //set Stopping time according to highest meson energy of the star
        G4double temfin;
        G4double TLab;
        std::vector<G4double> energies;
        std::vector<G4double> projections;
        ThreeVector ab, cd;
 
        for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
          energies.push_back((*pit)->getKineticEnergy());
          ab = (*pit)->boostVector();
          cd = (*pit)->getPosition();
          projections.push_back(ab.dot(cd)); //projection length
        }// make vector of energies
      
        temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
        TLab = *max_element(energies.begin(), energies.end()); //choose max energy
 
        // energy-dependent stopping time above 2 AGeV
        if(TLab>2000.)
          temfin *= (5.8E4-TLab)/5.6E4;
 
        maximumTime = temfin;  
 
        // If the incoming particle is slow, use a larger stopping time
        const G4double rMax = theNucleus->getUniverseRadius();
        const G4double distance = 2.*rMax;
        const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
        const G4double traversalTime = distance / maxMesonVelocityProjection;
        if(maximumTime < traversalTime)
          maximumTime = traversalTime;
        INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
        // Fill in the relevant kinematic variables
        theNucleus->setIncomingAngularMomentum(G4INCL::ThreeVector(0., 0., 0.));
        theNucleus->setIncomingMomentum(G4INCL::ThreeVector(0.,0.,0.));
        if(theNucleus->getAnnihilationType()==DNbarNPbarNType)
          theNucleus->setInitialEnergy(pr->getMass() + ParticleTable::getTableMass(theNucleus->getA() + 2, theNucleus->getZ(), theNucleus->getS()));
        else if(theNucleus->getAnnihilationType()==DNbarPPbarPType)
          theNucleus->setInitialEnergy(pr->getMass() + ParticleTable::getTableMass(theNucleus->getA() + 2, theNucleus->getZ() + 2,theNucleus->getS()));
        else if(theNucleus->getAnnihilationType() == DNbarNPbarPType || theNucleus->getAnnihilationType()==DNbarPPbarNType)
          theNucleus->setInitialEnergy(pr->getMass() + ParticleTable::getTableMass(theNucleus->getA() + 2, theNucleus->getZ() +1, theNucleus->getS()));
 
        for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
            (*p)->makeProjectileSpectator();                                        
        }
      
        generateAllAvatars();
        firstAvatar = false;
 
        IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
        delete pr;
        delete obj;
        theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
        INCL_DEBUG("Avatars added" << '\n');
        return 99.;
      }
    }

Referenced by shoot().

shootParticle()

G4double G4INCL::StandardPropagationModel::shootParticle ( ParticleType const t, const G4double kineticEnergy, const G4double impactParameter, const G4double phi )

virtual

Implements G4INCL::IPropagationModel.

Definition at line 254 of file G4INCLStandardPropagationModel.cc.

                                                                                                                                                              {
      theNucleus->setParticleNucleusCollision();
      currentTime = 0.0;
 
      // Create the projectile particle
      const G4double projectileMass = ParticleTable::getTableParticleMass(type);
      G4double energy = kineticEnergy + projectileMass;
      G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
      ThreeVector momentum(0.0, 0.0, momentumZ);
      Particle *p= new G4INCL::Particle(type, energy, momentum, ThreeVector());
 
      G4double temfin;
      G4double TLab;
      if( p->isMeson()) {
        temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
        TLab = p->getKineticEnergy();
      } else {
        temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
        TLab = p->getKineticEnergy()/p->getA();
      }
 
      // energy-dependent stopping time above 2 AGeV
      if(TLab>2000.)
        temfin *= (5.8E4-TLab)/5.6E4;
 
      maximumTime = temfin;
 
      // If the incoming particle is slow, use a larger stopping time
      const G4double rMax = theNucleus->getUniverseRadius();
      const G4double distance = 2.*rMax;
      const G4double projectileVelocity = p->boostVector().mag();
      const G4double traversalTime = distance / projectileVelocity;
      if(maximumTime < traversalTime)
        maximumTime = traversalTime;
      INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
 
      // If the incoming particle is an antinucleon use a larger stopping time
      if( p->isAntiNucleon()) maximumTime *= 2.;
 
      // If Coulomb is activated, do not process events with impact
      // parameter larger than the maximum impact parameter, taking into
      // account Coulomb distortion.
      if(impactParameter>CoulombDistortion::maxImpactParameter(p->getSpecies(), kineticEnergy, theNucleus)) {
        INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
        delete p;
        return -1.;
      }
 
      ThreeVector position(impactParameter * std::cos(phi),
          impactParameter * std::sin(phi),
          0.);
      p->setPosition(position);
 
      // Fill in the relevant kinematic variables
      theNucleus->setIncomingAngularMomentum(p->getAngularMomentum());
      theNucleus->setIncomingMomentum(p->getMomentum());
      theNucleus->setInitialEnergy(p->getEnergy()
          + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ(),theNucleus->getS()));
 
      // Reset the particle kinematics to the INCL values
      p->setINCLMass();
      p->setEnergy(p->getMass() + kineticEnergy);
      p->adjustMomentumFromEnergy();
 
      p->makeProjectileSpectator();
      generateAllAvatars();
      firstAvatar = false;
 
      // Get the entry avatars from Coulomb and put them in the Store
      ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurface(p, theNucleus);
      if(theEntryAvatar) {
        theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
 
        return p->getTransversePosition().mag();
      } else {
        delete p;
        return -1.;
      }
    }

Referenced by shoot().

updateAvatars()

void G4INCL::StandardPropagationModel::updateAvatars

(

const ParticleList &

particles

)

Update all avatars related to a particle.

Definition at line 868 of file G4INCLStandardPropagationModel.cc.

                                                                              {
 
      for(ParticleIter iter=particles.begin(), e=particles.end(); iter!=e; ++iter) {
        G4double time = this->getReflectionTime(*iter);
        if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
      }
      ParticleList const &p = theNucleus->getStore()->getParticles();
      generateUpdatedCollisions(particles, p);             // Predict collisions with spectators and participants
    }

Referenced by propagate().

---

The documentation for this class was generated from the following files:

• G4INCLStandardPropagationModel.hh
• G4INCLStandardPropagationModel.cc