G4INCLNucleus.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"

/*
 * G4INCLNucleus.hh
 *
 *  \date Jun 5, 2009
 * \author Pekka Kaitaniemi
 */

#ifndef G4INCLNUCLEUS_HH_
#define G4INCLNUCLEUS_HH_

#include <list>
#include <string>

#include "G4INCLParticle.hh"
#include "G4INCLEventInfo.hh"
#include "G4INCLCluster.hh"
#include "G4INCLFinalState.hh"
#include "G4INCLStore.hh"
#include "G4INCLGlobals.hh"
#include "G4INCLParticleTable.hh"
#include "G4INCLConfig.hh"
#include "G4INCLConfigEnums.hh"
#include "G4INCLCluster.hh"
#include "G4INCLProjectileRemnant.hh"

namespace G4INCL {

  class Nucleus : public Cluster {
  public:
    Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius=-1.);
    virtual ~Nucleus();

    void initializeParticles();

    void insertParticle(Particle *p) {
      theZ += p->getZ();
      theA += p->getA();
      theStore->particleHasEntered(p);
      if(p->isNucleon()) {
        theNpInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
        theNnInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
      }
      if(!p->isTargetSpectator()) theStore->getBook()->incrementCascading();
    };

    void applyFinalState(FinalState *);

    G4int getInitialA() const { return theInitialA; };
    G4int getInitialZ() const { return theInitialZ; };

    ParticleList const &getCreatedParticles() const { return justCreated; }

    ParticleList const &getUpdatedParticles() const { return toBeUpdated; }

    Particle *getBlockedDelta() const { return blockedDelta; }

    void propagateParticles(G4double step);

    G4int getNumberOfEnteringProtons() const { return theNpInitial; };
    G4int getNumberOfEnteringNeutrons() const { return theNnInitial; };

    G4double computeSeparationEnergyBalance() const {
      G4double S = 0.0;
      ParticleList outgoing = theStore->getOutgoingParticles();
      for(ParticleIter i = outgoing.begin(); i != outgoing.end(); ++i) {
        const ParticleType t = (*i)->getType();
        switch(t) {
          case Proton:
          case Neutron:
          case DeltaPlusPlus:
          case DeltaPlus:
          case DeltaZero:
          case DeltaMinus:
            S += thePotential->getSeparationEnergy(*i);
            break;
          case Composite:
            S += (*i)->getZ() * thePotential->getSeparationEnergy(Proton)
              + ((*i)->getA() - (*i)->getZ()) * thePotential->getSeparationEnergy(Neutron);
            break;
          case PiPlus:
            S += thePotential->getSeparationEnergy(Proton) - thePotential->getSeparationEnergy(Neutron);
            break;
          case PiMinus:
            S += thePotential->getSeparationEnergy(Neutron) - thePotential->getSeparationEnergy(Proton);
            break;
          default:
            break;
        }
      }

      S -= theNpInitial * thePotential->getSeparationEnergy(Proton);
      S -= theNnInitial * thePotential->getSeparationEnergy(Neutron);
      return S;
    }

    G4bool decayOutgoingDeltas();

    G4bool decayInsideDeltas();

    G4bool decayOutgoingClusters();

    G4bool decayMe();

    void emitInsidePions();

    void computeRecoilKinematics();

    ThreeVector computeCenterOfMass() const;

    G4double computeTotalEnergy() const;

    G4double computeExcitationEnergy() const;

    void setIncomingAngularMomentum(const ThreeVector &j) {
      incomingAngularMomentum = j;
    }

    const ThreeVector &getIncomingAngularMomentum() const { return incomingAngularMomentum; }

    void setIncomingMomentum(const ThreeVector &p) {
      incomingMomentum = p;
    }

    const ThreeVector &getIncomingMomentum() const {
      return incomingMomentum;
    }

    void setInitialEnergy(const G4double e) { initialEnergy = e; }

    G4double getInitialEnergy() const { return initialEnergy; }

    G4double getExcitationEnergy() const { return theExcitationEnergy; }

    inline G4bool containsDeltas() {
      ParticleList inside = theStore->getParticles();
      for(ParticleIter i=inside.begin(); i!=inside.end(); ++i)
        if((*i)->isDelta()) return true;
      return false;
    }

    std::string print();

    std::string dump();

    Store* getStore() const {return theStore; };
    void setStore(Store *s) {
      delete theStore;
      theStore = s;
    };

    G4double getInitialInternalEnergy() const { return initialInternalEnergy; };

    G4bool isEventTransparent() const;

    G4bool hasRemnant() const { return remnant; }

    //    void fillEventInfo(Results::EventInfo *eventInfo);
    void fillEventInfo(EventInfo *eventInfo);

    G4bool getTryCompoundNucleus() { return tryCN; }

    G4int getProjectileChargeNumber() const { return projectileZ; }

    G4int getProjectileMassNumber() const { return projectileA; }

    void setProjectileChargeNumber(G4int n) { projectileZ=n; }

    void setProjectileMassNumber(G4int n) { projectileA=n; }

    G4bool isForcedTransparent() { return forceTransparent; }

    G4double getTransmissionBarrier(Particle const * const p) {
      const G4double theTransmissionRadius = theDensity->getTransmissionRadius(p);
      const G4double theParticleZ = p->getZ();
      return PhysicalConstants::eSquared*(theZ-theParticleZ)*theParticleZ/theTransmissionRadius;
    }

    struct ConservationBalance {
      ThreeVector momentum;
      G4double energy;
      G4int Z, A;
    };

    ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const;

    void useFusionKinematics();

    G4double getSurfaceRadius(Particle const * const particle) const {
      if(particle->isPion())
        // Temporarily set RPION = RMAX
        return getUniverseRadius();
        //return 0.5*(theDensity->getTransmissionRadius(particle)+getUniverseRadius());
      else {
        const G4double pr = particle->getMomentum().mag()/thePotential->getFermiMomentum(particle);
        if(pr>=1.)
          return getUniverseRadius();
        else
          return theDensity->getMaxRFromP(pr);
      }
    }

    G4double getCoulombRadius(ParticleSpecies const &p) const {
      if(p.theType == Composite) {
        const G4int zp = p.theZ;
        const G4int ap = p.theA;
        G4double barr, radius = 0.;
        if(zp==1 && ap==2) { // d
          barr = 0.2565*Math::pow23((G4double)theA)-0.78;
          radius = ParticleTable::eSquared*zp*theZ/barr - 2.5;
        } else if((zp==1 || zp==2) && ap==3) { // t, He3
          barr = 0.5*(0.5009*Math::pow23((G4double)theA)-1.16);
          radius = ParticleTable::eSquared*theZ/barr - 0.5;
        } else if(zp==2 && ap==4) { // alpha
          barr = 0.5939*Math::pow23((G4double)theA)-1.64;
          radius = ParticleTable::eSquared*zp*theZ/barr - 0.5;
        } else if(zp>2) {
          radius = getUniverseRadius();
        }
        if(radius<=0.) {
          ERROR("Negative Coulomb radius! Using the universe radius" << std::endl);
          radius = getUniverseRadius();
        }
        return radius;
      } else
        return getUniverseRadius();
    }

    G4double getUniverseRadius() const { return theUniverseRadius; }

    void setUniverseRadius(const G4double universeRadius) { theUniverseRadius=universeRadius; }

    G4bool isNucleusNucleusCollision() const { return isNucleusNucleus; }

    void setNucleusNucleusCollision() { isNucleusNucleus=true; }

    void setParticleNucleusCollision() { isNucleusNucleus=false; }

    void setProjectileRemnant(ProjectileRemnant * const c) {
      delete theProjectileRemnant;
      theProjectileRemnant = c;
    }

    ProjectileRemnant *getProjectileRemnant() const { return theProjectileRemnant; }

    void deleteProjectileRemnant() {
      delete theProjectileRemnant;
      theProjectileRemnant = NULL;
    }

    void moveProjectileRemnantComponentsToOutgoing() {
      theStore->addToOutgoing(theProjectileRemnant->getParticles());
      theProjectileRemnant->clearParticles();
    }

    void finalizeProjectileRemnant(const G4double emissionTime);

    inline void updatePotentialEnergy(Particle *p) {
      p->setPotentialEnergy(thePotential->computePotentialEnergy(p));
    }

    NuclearDensity* getDensity() const { return theDensity; };

    NuclearPotential::INuclearPotential* getPotential() const { return thePotential; };

  private:
    void computeOneNucleonRecoilKinematics();

  private:
    G4int theInitialZ, theInitialA;
    G4int theNpInitial;
    G4int theNnInitial;
    G4double initialInternalEnergy;
    ThreeVector incomingAngularMomentum, incomingMomentum;
    ThreeVector initialCenterOfMass;
    G4bool remnant;

    ParticleList toBeUpdated;
    ParticleList justCreated;
    Particle *blockedDelta;
    G4double initialEnergy;
    Store *theStore;
    G4bool tryCN;
    G4bool forceTransparent;

    G4int projectileZ;
    G4int projectileA;

    G4double theUniverseRadius;

    G4bool isNucleusNucleus;

    ProjectileRemnant *theProjectileRemnant;

    NuclearDensity *theDensity;

    NuclearPotential::INuclearPotential *thePotential;

  };

}

#endif /* G4INCLNUCLEUS_HH_ */