G4INCLCascade.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 G4INCLCascade_hh
#define G4INCLCascade_hh 1

#include "G4INCLParticle.hh"
#include "G4INCLNucleus.hh"
#include "G4INCLIPropagationModel.hh"
#include "G4INCLCascadeAction.hh"
#include "G4INCLEventInfo.hh"
#include "G4INCLGlobalInfo.hh"
#include "G4INCLLogger.hh"
#include "G4INCLConfig.hh"
#include "G4INCLRootFinder.hh"

namespace G4INCL {
  class INCL {
    public:
      INCL(Config const * const config);

      ~INCL();

      INCL(const INCL &rhs);

      INCL &operator=(const INCL &rhs);

      G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z);
      G4bool initializeTarget(const G4int A, const G4int Z);
      inline const EventInfo &processEvent() {
        return processEvent(
            theConfig->getProjectileSpecies(),
            theConfig->getProjectileKineticEnergy(),
            theConfig->getTargetA(),
            theConfig->getTargetZ()
            );
      }
      const EventInfo &processEvent(
          ParticleSpecies const &projectileSpecies,
          const G4double kineticEnergy,
          const G4int targetA,
          const G4int targetZ
          );

      void finalizeGlobalInfo();
      const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }

    private:
      IPropagationModel *propagationModel;
      G4int theA, theZ;
      G4bool targetInitSuccess;
      G4double maxImpactParameter;
      G4double maxUniverseRadius;
      G4double maxInteractionDistance;
      G4double fixedImpactParameter;
      CascadeAction *cascadeAction;
      Config const * const theConfig;
      Nucleus *nucleus;
      G4bool forceTransparent;

      EventInfo theEventInfo;
      GlobalInfo theGlobalInfo;

      G4int minRemnantSize;

      class RecoilFunctor : public RootFunctor {
        public:
          RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
            RootFunctor(0., 1E6),
            nucleus(n),
            outgoingParticles(n->getStore()->getOutgoingParticles()),
            theEventInfo(ei) {
              for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
                particleMomenta.push_back((*p)->getMomentum());
                particleKineticEnergies.push_back((*p)->getKineticEnergy());
              }
              ProjectileRemnant * const aPR = n->getProjectileRemnant();
              if(aPR && aPR->getA()>0) {
                particleMomenta.push_back(aPR->getMomentum());
                particleKineticEnergies.push_back(aPR->getKineticEnergy());
                outgoingParticles.push_back(aPR);
              }
            }
          virtual ~RecoilFunctor() {}

          G4double operator()(const G4double x) const {
            scaleParticleEnergies(x);
            return nucleus->getConservationBalance(theEventInfo,true).energy;
          }

          void cleanUp(const G4bool success) const {
            if(!success)
              scaleParticleEnergies(1.);
          }

        private:
          Nucleus *nucleus;
          ParticleList outgoingParticles;
          // \brief Reference to the EventInfo object
          EventInfo const &theEventInfo;
          std::list<ThreeVector> particleMomenta;
          std::list<G4double> particleKineticEnergies;

          void scaleParticleEnergies(const G4double rescale) const {
            // Rescale the energies (and the momenta) of the outgoing particles.
            ThreeVector pBalance = nucleus->getIncomingMomentum();
            std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
            std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
            for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP, ++iE)
            {
              const G4double mass = (*i)->getMass();
              const G4double newKineticEnergy = (*iE) * rescale;

              (*i)->setMomentum(*iP);
              (*i)->setEnergy(mass + newKineticEnergy);
              (*i)->adjustMomentumFromEnergy();

              pBalance -= (*i)->getMomentum();
            }

            nucleus->setMomentum(pBalance);
            const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
            const G4double pRem2 = pBalance.mag2();
            const G4double recoilEnergy = pRem2/
              (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
            nucleus->setEnergy(remnantMass + recoilEnergy);
          }
      };

      class RecoilCMFunctor : public RootFunctor {
        public:
          RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
            RootFunctor(0., 1E6),
            nucleus(n),
            theIncomingMomentum(nucleus->getIncomingMomentum()),
            outgoingParticles(n->getStore()->getOutgoingParticles()),
            theEventInfo(ei) {
              thePTBoostVector = nucleus->getIncomingMomentum()/nucleus->getInitialEnergy();
              for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
                (*p)->boost(thePTBoostVector);
                particleCMMomenta.push_back((*p)->getMomentum());
              }
              ProjectileRemnant * const aPR = n->getProjectileRemnant();
              if(aPR && aPR->getA()>0) {
                aPR->boost(thePTBoostVector);
                particleCMMomenta.push_back(aPR->getMomentum());
                outgoingParticles.push_back(aPR);
              }
            }
          virtual ~RecoilCMFunctor() {}

          G4double operator()(const G4double x) const {
            scaleParticleCMMomenta(x);
            return nucleus->getConservationBalance(theEventInfo,true).energy;
          }

          void cleanUp(const G4bool success) const {
            if(!success)
              scaleParticleCMMomenta(1.);
          }

        private:
          Nucleus *nucleus;
          ThreeVector thePTBoostVector;
          ThreeVector theIncomingMomentum;
          ParticleList outgoingParticles;
          // \brief Reference to the EventInfo object
          EventInfo const &theEventInfo;
          std::list<ThreeVector> particleCMMomenta;

          void scaleParticleCMMomenta(const G4double rescale) const {
            // Rescale the CM momenta of the outgoing particles.
            ThreeVector remnantMomentum = theIncomingMomentum;
            std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
            for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP)
            {
              (*i)->setMomentum(*iP * rescale);
              (*i)->adjustEnergyFromMomentum();
              (*i)->boost(-thePTBoostVector);

              remnantMomentum -= (*i)->getMomentum();
            }

            nucleus->setMomentum(remnantMomentum);
            const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
            const G4double pRem2 = remnantMomentum.mag2();
            const G4double recoilEnergy = pRem2/
              (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
            nucleus->setEnergy(remnantMass + recoilEnergy);
          }
      };

      void rescaleOutgoingForRecoil();

#ifndef INCLXX_IN_GEANT4_MODE

      void globalConservationChecks(G4bool afterRecoil);
#endif

      G4bool continueCascade();

      G4int makeProjectileRemnant();

      void makeCompoundNucleus();

      G4bool preCascade(ParticleSpecies const projectileSpecies, const G4double kineticEnergy);

      void cascade();

      void postCascade();

      void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);

      void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);

      void updateGlobalInfo();
  };
}

#endif