G4INCLCluster.hh

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
// Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
// Davide Mancusi, CEA
// Alain Boudard, CEA
// Sylvie Leray, CEA
// Joseph Cugnon, University of Liege
//
#define INCLXX_IN_GEANT4_MODE 1

#include "globals.hh"

#ifndef G4INCLCluster_hh
#define G4INCLCluster_hh 1

#include "G4INCLParticle.hh"
#include "G4INCLNuclearDensityFactory.hh"
#include "G4INCLParticleSampler.hh"

namespace G4INCL {

  class Cluster : public Particle {
  public:

    Cluster(const G4int Z, const G4int A, const G4bool createParticleSampler=true) :
      Particle(),
      theExcitationEnergy(0.),
      theSpin(0.,0.,0.),
      theParticleSampler(NULL)
    {
      setType(Composite);
      theZ = Z;
      theA = A;
      setINCLMass();
      if(createParticleSampler)
        theParticleSampler = new ParticleSampler(A,Z);
    }

    template<class Iterator>
      Cluster(Iterator begin, Iterator end) :
        Particle(),
        theExcitationEnergy(0.),
        theSpin(0.,0.,0.),
        theParticleSampler(NULL)
    {
      setType(Composite);
      for(Iterator i = begin; i != end; ++i) {
        addParticle(*i);
      }
      thePosition /= theA;
      setINCLMass();
      adjustMomentumFromEnergy();
    }

    virtual ~Cluster() {
      delete theParticleSampler;
    }

    Cluster(const Cluster &rhs) :
      Particle(rhs),
      theExcitationEnergy(rhs.theExcitationEnergy),
      theSpin(rhs.theSpin)
    {
      for(ParticleIter p=rhs.particles.begin(), e=rhs.particles.end(); p!=e; ++p) {
        particles.push_back(new Particle(**p));
      }
      if(rhs.theParticleSampler)
        theParticleSampler = new ParticleSampler(rhs.theA,rhs.theZ);
      else
        theParticleSampler = NULL;
    }

    Cluster &operator=(const Cluster &rhs) {
      Cluster temporaryCluster(rhs);
      Particle::operator=(temporaryCluster);
      swap(temporaryCluster);
      return *this;
    }

    void swap(Cluster &rhs) {
      Particle::swap(rhs);
      std::swap(theExcitationEnergy, rhs.theExcitationEnergy);
      std::swap(theSpin, rhs.theSpin);
      // std::swap is overloaded by std::list and guaranteed to operate in
      // constant time
      std::swap(particles, rhs.particles);
      std::swap(theParticleSampler, rhs.theParticleSampler);
    }

    ParticleSpecies getSpecies() const {
      return ParticleSpecies(theA, theZ);
    }

    void deleteParticles() {
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        delete (*p);
      }
      clearParticles();
    }

    void clearParticles() { particles.clear(); }

    void setZ(const G4int Z) { theZ = Z; }

    void setA(const G4int A) { theA = A; }

    G4double getExcitationEnergy() const { return theExcitationEnergy; }

    void setExcitationEnergy(const G4double e) { theExcitationEnergy=e; }

    inline virtual G4double getTableMass() const { return getRealMass(); }

    ParticleList const &getParticles() const { return particles; }

    void removeParticle(Particle * const p) { particles.remove(p); }

    void addParticle(Particle * const p) {
// assert(p->isNucleon());
      particles.push_back(p);
      theEnergy += p->getEnergy();
      thePotentialEnergy += p->getPotentialEnergy();
      theMomentum += p->getMomentum();
      thePosition += p->getPosition();
      theA += p->getA();
      theZ += p->getZ();
      nCollisions += p->getNumberOfCollisions();
    }

    void addParticles(ParticleList const &pL) {
      for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p)
        addParticle(*p);
    }

    ParticleList getParticleList() const { return particles; }

    std::string print() const {
      std::stringstream ss;
      ss << "Cluster (ID = " << ID << ") type = ";
      ss << ParticleTable::getName(theType);
      ss << std::endl
        << "   A = " << theA << std::endl
        << "   Z = " << theZ << std::endl
        << "   mass = " << getMass() << std::endl
        << "   energy = " << theEnergy << std::endl
        << "   momentum = "
        << theMomentum.print()
        << std::endl
        << "   position = "
        << thePosition.print()
        << std::endl
        << "Contains the following particles:"
        << std::endl;
      for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i)
        ss << (*i)->print();
      ss << std::endl;
      return ss.str();
    }

    virtual void initializeParticles();

    void internalBoostToCM() {

      // First compute the current CM position and total momentum
      ThreeVector theCMPosition, theTotalMomentum;
      G4double theTotalEnergy = 0.0;
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        theCMPosition += (*p)->getPosition();
        theTotalMomentum += (*p)->getMomentum();
        theTotalEnergy += (*p)->getEnergy();
      }
      theCMPosition /= theA;
// assert((unsigned int)theA==particles.size());

      // Now determine the CM velocity of the particles
      // commented out because currently unused, see below
      // ThreeVector betaCM = theTotalMomentum / theTotalEnergy;

      // The new particle positions and momenta are scaled by a factor of
      // \f$\sqrt{A/(A-1)}\f$, so that the resulting density distributions in
      // the CM have the same variance as the one we started with.
      const G4double rescaling = std::sqrt(((G4double)theA)/((G4double)(theA-1)));

      // Loop again to boost and reposition
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        // \bug{We should do the following, but the Fortran version actually
        // does not!
        // (*p)->boost(betaCM);
        // Here is what the Fortran version does:}
        (*p)->setMomentum(((*p)->getMomentum()-theTotalMomentum/theA)*rescaling);

        // Set the CM position of the particles
        (*p)->setPosition(((*p)->getPosition()-theCMPosition)*rescaling);
      }

      // Set the global cluster kinematic variables
      thePosition.setX(0.0);
      thePosition.setY(0.0);
      thePosition.setZ(0.0);
      theMomentum.setX(0.0);
      theMomentum.setY(0.0);
      theMomentum.setZ(0.0);
      theEnergy = getMass();

      INCL_DEBUG("Cluster boosted to internal CM:" << std::endl << print());

    }

    void putParticlesOffShell() {
      // Compute the dynamical potential
      const G4double theDynamicalPotential = computeDynamicalPotential();
      INCL_DEBUG("The dynamical potential is " << theDynamicalPotential << " MeV" << std::endl);

      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        const G4double energy = (*p)->getEnergy() - theDynamicalPotential;
        const ThreeVector &momentum = (*p)->getMomentum();
        // Here particles are put off-shell so that we can satisfy the energy-
        // and momentum-conservation laws
        (*p)->setEnergy(energy);
        (*p)->setMass(std::sqrt(energy*energy - momentum.mag2()));
      }
      INCL_DEBUG("Cluster components are now off shell:" << std::endl
            << print());
    }

    void setPosition(const ThreeVector &position) {
      ThreeVector shift(position-thePosition);
      Particle::setPosition(position);
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->setPosition((*p)->getPosition()+shift);
      }
    }

    void boost(const ThreeVector &aBoostVector) {
      Particle::boost(aBoostVector);
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->boost(aBoostVector);
        // Apply Lorentz contraction to the particle position
        (*p)->lorentzContract(aBoostVector,thePosition);
        (*p)->rpCorrelate();
      }

      INCL_DEBUG("Cluster was boosted with (bx,by,bz)=("
          << aBoostVector.getX() << ", " << aBoostVector.getY() << ", " << aBoostVector.getZ() << "):"
          << std::endl << print());

    }

    void freezeInternalMotion() {
      const ThreeVector &normMomentum = theMomentum / getMass();
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        const G4double pMass = (*p)->getMass();
        const ThreeVector frozenMomentum = normMomentum * pMass;
        const G4double frozenEnergy = std::sqrt(frozenMomentum.mag2()+pMass*pMass);
        (*p)->setFrozenMomentum(frozenMomentum);
        (*p)->setFrozenEnergy(frozenEnergy);
        (*p)->freezePropagation();
      }
    }

    virtual void rotate(const G4double angle, const ThreeVector &axis) {
      Particle::rotate(angle, axis);
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->rotate(angle, axis);
      }
    }

    virtual void makeProjectileSpectator() {
      Particle::makeProjectileSpectator();
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->makeProjectileSpectator();
      }
    }

    virtual void makeTargetSpectator() {
      Particle::makeTargetSpectator();
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->makeTargetSpectator();
      }
    }

    virtual void makeParticipant() {
      Particle::makeParticipant();
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        (*p)->makeParticipant();
      }
    }

    ThreeVector const &getSpin() const { return theSpin; }

    void setSpin(const ThreeVector &j) { theSpin = j; }

    G4INCL::ThreeVector getAngularMomentum() const {
      return Particle::getAngularMomentum() + getSpin();
    }

  private:
    G4double computeDynamicalPotential() {
      G4double theDynamicalPotential = 0.0;
      for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
        theDynamicalPotential += (*p)->getEnergy();
      }
      theDynamicalPotential -= getTableMass();
      theDynamicalPotential /= theA;

      return theDynamicalPotential;
    }

  protected:
    ParticleList particles;
    G4double theExcitationEnergy;
    ThreeVector theSpin;
    ParticleSampler *theParticleSampler;

  };

}

#endif