Gabriele Cosmo
393 lines
16 KiB
C++
393 lines
16 KiB
C++
//
|
|
// ********************************************************************
|
|
// * 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"
|
|
|
|
#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"
|
|
#ifndef INCLXX_IN_GEANT4_MODE
|
|
#include "G4INCLGeant4Compat.hh"
|
|
#endif
|
|
|
|
|
|
namespace G4INCL {
|
|
class INCL {
|
|
public:
|
|
INCL(Config const * const config);
|
|
|
|
~INCL();
|
|
|
|
/// \brief Dummy copy constructor to silence Coverity warning
|
|
INCL(const INCL &rhs);
|
|
|
|
/// \brief Dummy assignment operator to silence Coverity warning
|
|
INCL &operator=(const INCL &rhs);
|
|
|
|
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S);
|
|
|
|
G4bool initializeTarget(const G4int A, const G4int Z, const G4int S, AnnihilationType theAType);
|
|
G4double initUniverseRadiusForAntiprotonAtRest(const G4int A, const G4int Z);
|
|
inline const EventInfo &processEvent() {
|
|
return processEvent(
|
|
theConfig->getProjectileSpecies(),
|
|
theConfig->getProjectileKineticEnergy(),
|
|
theConfig->getTargetA(),
|
|
theConfig->getTargetZ(),
|
|
theConfig->getTargetS()
|
|
);
|
|
}
|
|
const EventInfo &processEvent(
|
|
ParticleSpecies const &projectileSpecies,
|
|
const G4double kineticEnergy,
|
|
const G4int targetA,
|
|
const G4int targetZ,
|
|
const G4int targetS
|
|
);
|
|
|
|
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds);
|
|
const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
|
|
|
|
|
|
private:
|
|
IPropagationModel *propagationModel;
|
|
G4int theA, theZ, theS;
|
|
G4bool targetInitSuccess;
|
|
G4double maxImpactParameter;
|
|
G4double maxUniverseRadius;
|
|
G4double maxInteractionDistance;
|
|
G4double fixedImpactParameter;
|
|
CascadeAction *cascadeAction;
|
|
Config const * const theConfig;
|
|
Nucleus *nucleus;
|
|
G4bool forceTransparent;
|
|
|
|
EventInfo theEventInfo;
|
|
GlobalInfo theGlobalInfo;
|
|
|
|
/// \brief Remnant size below which cascade stops
|
|
G4int minRemnantSize;
|
|
|
|
/// \brief Class to adjust remnant recoil
|
|
class RecoilFunctor : public RootFunctor {
|
|
public:
|
|
/** \brief Prepare for calling the () operator and scaleParticleEnergies
|
|
*
|
|
* The constructor sets the private class members.
|
|
*/
|
|
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 || aPR->getA()<0)) {
|
|
particleMomenta.push_back(aPR->getMomentum());
|
|
particleKineticEnergies.push_back(aPR->getKineticEnergy());
|
|
outgoingParticles.push_back(aPR);
|
|
}
|
|
}
|
|
virtual ~RecoilFunctor() {}
|
|
|
|
/** \brief Compute the energy-conservation violation.
|
|
*
|
|
* \param x scale factor for the particle energies
|
|
* \return the energy-conservation violation
|
|
*/
|
|
G4double operator()(const G4double x) const {
|
|
scaleParticleEnergies(x);
|
|
return nucleus->getConservationBalance(theEventInfo,true).energy;
|
|
}
|
|
|
|
/// \brief Clean up after root finding
|
|
void cleanUp(const G4bool success) const {
|
|
if(!success)
|
|
scaleParticleEnergies(1.);
|
|
}
|
|
|
|
private:
|
|
/// \brief Pointer to the nucleus
|
|
Nucleus *nucleus;
|
|
/// \brief List of final-state particles.
|
|
ParticleList outgoingParticles;
|
|
// \brief Reference to the EventInfo object
|
|
EventInfo const &theEventInfo;
|
|
/// \brief Initial momenta of the outgoing particles
|
|
std::list<ThreeVector> particleMomenta;
|
|
/// \brief Initial kinetic energies of the outgoing particles
|
|
std::list<G4double> particleKineticEnergies;
|
|
|
|
/** \brief Scale the kinetic energies of the outgoing particles.
|
|
*
|
|
* \param rescale scale factor
|
|
*/
|
|
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->getS()) + nucleus->getExcitationEnergy();
|
|
const G4double pRem2 = pBalance.mag2();
|
|
const G4double recoilEnergy = pRem2/
|
|
(std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
|
|
nucleus->setEnergy(remnantMass + recoilEnergy);
|
|
}
|
|
};
|
|
|
|
/// \brief Class to adjust remnant recoil in the reaction CM system
|
|
class RecoilCMFunctor : public RootFunctor {
|
|
public:
|
|
/** \brief Prepare for calling the () operator and scaleParticleEnergies
|
|
*
|
|
* The constructor sets the private class members.
|
|
*/
|
|
RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
|
|
RootFunctor(0., 1E6),
|
|
nucleus(n),
|
|
theIncomingMomentum(nucleus->getIncomingMomentum()),
|
|
outgoingParticles(n->getStore()->getOutgoingParticles()),
|
|
theEventInfo(ei) {
|
|
if(theIncomingMomentum.mag() == 0.){ //change the condition
|
|
thePTBoostVector = {0., 0., 0.};
|
|
//INCL_WARN("PTBoostVector at rest is zero" << '\n');
|
|
}
|
|
else{
|
|
thePTBoostVector = nucleus->getIncomingMomentum()/(nucleus->getInitialEnergy()); //D
|
|
//INCL_WARN("PTBoostVector" << '\n');
|
|
}
|
|
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() {}
|
|
|
|
/** \brief Compute the energy-conservation violation.
|
|
*
|
|
* \param x scale factor for the particle energies
|
|
* \return the energy-conservation violation
|
|
*/
|
|
G4double operator()(const G4double x) const {
|
|
scaleParticleCMMomenta(x);
|
|
return nucleus->getConservationBalance(theEventInfo,true).energy;
|
|
}
|
|
|
|
/// \brief Clean up after root finding
|
|
void cleanUp(const G4bool success) const {
|
|
if(!success)
|
|
scaleParticleCMMomenta(1.);
|
|
}
|
|
|
|
private:
|
|
/// \brief Pointer to the nucleus
|
|
Nucleus *nucleus;
|
|
/// \brief Projectile-target CM boost vector
|
|
ThreeVector thePTBoostVector;
|
|
/// \brief Incoming momentum
|
|
ThreeVector theIncomingMomentum;
|
|
/// \brief List of final-state particles.
|
|
ParticleList outgoingParticles;
|
|
/// \brief Reference to the EventInfo object
|
|
EventInfo const &theEventInfo;
|
|
/// \brief Initial CM momenta of the outgoing particles
|
|
std::list<ThreeVector> particleCMMomenta;
|
|
|
|
/** \brief Scale the kinetic energies of the outgoing particles.
|
|
*
|
|
* \param rescale scale factor
|
|
*/
|
|
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->getS()) + nucleus->getExcitationEnergy();
|
|
const G4double pRem2 = remnantMomentum.mag2();
|
|
const G4double recoilEnergy = pRem2/
|
|
(std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
|
|
nucleus->setEnergy(remnantMass + recoilEnergy);
|
|
}
|
|
};
|
|
|
|
/** \brief Rescale the energies of the outgoing particles.
|
|
*
|
|
* Allow for the remnant recoil energy by rescaling the energy (and
|
|
* momenta) of the outgoing particles.
|
|
*/
|
|
void rescaleOutgoingForRecoil();
|
|
|
|
#ifndef INCLXX_IN_GEANT4_MODE
|
|
/** \brief Run global conservation checks
|
|
*
|
|
* Check that energy and momentum are correctly conserved. If not, issue
|
|
* a warning.
|
|
*
|
|
* Also feeds the balance variables in theEventInfo.
|
|
*
|
|
* \param afterRecoil whether to take into account nuclear recoil
|
|
*/
|
|
void globalConservationChecks(G4bool afterRecoil);
|
|
#endif
|
|
|
|
/** \brief Stopping criterion for the cascade
|
|
*
|
|
* Returns true if the cascade should continue, and false if any of the
|
|
* stopping criteria is satisfied.
|
|
*/
|
|
G4bool continueCascade();
|
|
|
|
/** \brief Make a projectile pre-fragment out of geometrical spectators
|
|
*
|
|
* The projectile pre-fragment is assigned an excitation energy given
|
|
* by \f$E_\mathrm{sp}-E_\mathrm{i,A}\f$, where \f$E_\mathrm{sp}\f$ is the
|
|
* sum of the energies of the spectator particles, and \f$E_\mathrm{i,A}\f$
|
|
* is the sum of the smallest \f$A\f$ particle energies initially present
|
|
* in the projectile, \f$A\f$ being the mass of the projectile
|
|
* pre-fragment. This is equivalent to assuming that the excitation
|
|
* energy is given by the sum of the transitions of all excited
|
|
* projectile components to the "holes" left by the participants.
|
|
*
|
|
* This method can modify the outgoing list and adds a projectile
|
|
* pre-fragment.
|
|
*
|
|
* \return the number of dynamical spectators that were merged back in
|
|
* the projectile
|
|
*/
|
|
G4int makeProjectileRemnant();
|
|
|
|
/** \brief Make a compound nucleus
|
|
*
|
|
* Selects the projectile components that can actually enter their
|
|
* potential and puts them into the target nucleus. If the CN excitation
|
|
* energy turns out to be negative, the event is considered a
|
|
* transparent. This method modifies theEventInfo and theGlobalInfo.
|
|
*/
|
|
void makeCompoundNucleus();
|
|
|
|
/// \brief Initialise the cascade
|
|
G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/// \brief The actual cascade loop
|
|
void cascade();
|
|
|
|
/// \brief Finalise the cascade and clean up
|
|
void postCascade(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/** \brief Initialise the maximum interaction distance.
|
|
*
|
|
* Used in forced CN events.
|
|
*/
|
|
void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
|
|
|
|
/** \brief Initialize the universe radius.
|
|
*
|
|
* Used for determining the energy-dependent size of the volume particles
|
|
* live in.
|
|
*/
|
|
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
|
|
|
|
/// \brief Update global counters and other members of theGlobalInfo object
|
|
void updateGlobalInfo();
|
|
|
|
/// \brief Fill probabilities and particle types from datafile and return probability sum for normalization
|
|
G4double read_file(std::string filename, std::vector<G4double>& probabilities,
|
|
std::vector<std::vector<G4String>>& particle_types);
|
|
|
|
/// \brief This function will tell you the FS line number from the datafile based on your random value
|
|
G4int findStringNumber(G4double rdm, std::vector<G4double> yields);
|
|
|
|
/// \brief Initialise the "cascade" for pbar on H1
|
|
void preCascade_pbarH1(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/// \brief Initialise the "cascade" for pbar on H2
|
|
void preCascade_pbarH2(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/// \brief Initialise the "cascade" for nbar on H1
|
|
void preCascade_nbarH1(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/// \brief Initialise the "cascade" for nbar on H2
|
|
void preCascade_nbarH2(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
|
|
|
|
/// \brief Finalise the "cascade" and clean up for pbar/nbar on H1
|
|
void postCascade_pbarH1(ParticleList const &outgoingParticles);
|
|
|
|
/// \brief Finalise the "cascade" and clean up for pbar/nbar on H2
|
|
void postCascade_pbarH2(ParticleList const &outgoingParticles, ParticleList const &H2Particles);
|
|
};
|
|
}
|
|
|
|
#endif
|