G4INCLStandardPropagationModel.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// INCL++ intra-nuclear cascade model
// Alain Boudard, CEA-Saclay, France
// Joseph Cugnon, University of Liege, Belgium
// Jean-Christophe David, CEA-Saclay, France
// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
// Sylvie Leray, CEA-Saclay, France
// Davide Mancusi, CEA-Saclay, France
//
#define INCLXX_IN_GEANT4_MODE 1
 
#include "globals.hh"
 
/*
 * StandardPropagationModel.cpp
 *
 *  \date 4 juin 2009
 * \author Pekka Kaitaniemi
 */
 
#include "G4INCLStandardPropagationModel.hh"
#include "G4INCLPbarAtrestEntryChannel.hh"
#include "G4INCLNbarAtrestEntryChannel.hh"
#include "G4INCLAntinucleiAtrestEntryChannel.hh"
#include "G4INCLSurfaceAvatar.hh"
#include "G4INCLBinaryCollisionAvatar.hh"
#include "G4INCLDecayAvatar.hh"
#include "G4INCLCrossSections.hh"
#include "G4INCLRandom.hh"
#include <iostream>
#include <algorithm>
#include "G4INCLLogger.hh"
#include "G4INCLGlobals.hh"
#include "G4INCLKinematicsUtils.hh"
#include "G4INCLCoulombDistortion.hh"
#include "G4INCLDeltaDecayChannel.hh"
#include "G4INCLSigmaZeroDecayChannel.hh"
#include "G4INCLPionResonanceDecayChannel.hh"
#include "G4INCLParticleEntryAvatar.hh"
#include "G4INCLIntersection.hh"
#include <vector>
 
namespace G4INCL {
 
    StandardPropagationModel::StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType, const G4double hTime)
      :theNucleus(0), maximumTime(70.0), currentTime(0.0),
      hadronizationTime(hTime),
      firstAvatar(true),
      theLocalEnergyType(localEnergyType),
      theLocalEnergyDeltaType(localEnergyDeltaType)
    {
    }
 
    StandardPropagationModel::~StandardPropagationModel()
    {
      delete theNucleus;
    }
 
    G4INCL::Nucleus* StandardPropagationModel::getNucleus()
    {
      return theNucleus;
    }
 
//D
 
    G4double StandardPropagationModel::shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
      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);
      }
    }
 
//D
 
    G4double StandardPropagationModel::shootAtrest(ParticleType const t, const G4double kineticEnergy) {
      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.;
    } 
       
//D
      
    G4double StandardPropagationModel::shootParticle(ParticleType const type, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
      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.;
      }
    }
 
    G4double StandardPropagationModel::shootComposite(ParticleSpecies const &species, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
      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();
    }
 
    G4double StandardPropagationModel::shootCompositeAtrest(ParticleSpecies const &species, const G4double kineticEnergy){
      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.;
      }
    }
 
    G4double StandardPropagationModel::getStoppingTime() {
      return maximumTime;
    }
 
    void StandardPropagationModel::setStoppingTime(G4double time) {
// assert(time>0.0);
      maximumTime = time;
    }
 
    G4double StandardPropagationModel::getCurrentTime() {
      return currentTime;
    }
 
    void StandardPropagationModel::setNucleus(G4INCL::Nucleus *nucleus)
    {
      theNucleus = nucleus;
    }
 
    void StandardPropagationModel::registerAvatar(G4INCL::IAvatar *anAvatar)
    {
      if(anAvatar) theNucleus->getStore()->add(anAvatar);
    }
 
    IAvatar *StandardPropagationModel::generateBinaryCollisionAvatar(Particle * const p1, Particle * const p2) {
      // 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);
    }
 
    G4double StandardPropagationModel::getReflectionTime(G4INCL::Particle const * const aParticle) {
      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;
    }
 
    G4double StandardPropagationModel::getTime(G4INCL::Particle const * const particleA,
                         G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
    {
      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;
    }
 
    void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
 
      // 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));
        }
      }
    }
 
    void StandardPropagationModel::generateCollisions(const ParticleList &particles) {
      // 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));
        }
      }
    }
 
    void StandardPropagationModel::generateCollisions(const ParticleList &particles, const ParticleList &except) {
 
      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));
        }
      }
 
    }
 
    void StandardPropagationModel::updateAvatars(const ParticleList &particles) {
 
      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
    }
 
    void StandardPropagationModel::generateAllAvatars() {
      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);
    }
 
#ifdef INCL_REGENERATE_AVATARS
    void StandardPropagationModel::generateAllAvatarsExceptUpdated(FinalState const * const fs) {
      ParticleList const &particles = theNucleus->getStore()->getParticles();
// assert(!particles.empty());
      for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
        G4double time = this->getReflectionTime(*i);
        if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
      }
      ParticleList except = fs->getModifiedParticles();
      ParticleList const &entering = fs->getEnteringParticles();
      except.insert(except.end(), entering.begin(), entering.end());
      generateCollisions(particles,except);
      generateDecays(particles);
    }
#endif
 
    void StandardPropagationModel::generateDecays(const ParticleList &particles) {
      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));
          }
        }
      }
    }
 
    G4INCL::IAvatar* StandardPropagationModel::propagate(FinalState const * const fs)
    {
      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;
    }
}