include/Geant4/G4INCLNucleus.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // INCL++ intra-nuclear cascade model
0027 // Alain Boudard, CEA-Saclay, France
0028 // Joseph Cugnon, University of Liege, Belgium
0029 // Jean-Christophe David, CEA-Saclay, France
0030 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
0031 // Sylvie Leray, CEA-Saclay, France
0032 // Davide Mancusi, CEA-Saclay, France
0033 //
0034 #define INCLXX_IN_GEANT4_MODE 1
0035
0036 #include "globals.hh"
0037
0038 /*
0039  * G4INCLNucleus.hh
0040  *
0041  *  \date Jun 5, 2009
0042  * \author Pekka Kaitaniemi
0043  */
0044
0045 #ifndef G4INCLNUCLEUS_HH_
0046 #define G4INCLNUCLEUS_HH_
0047
0048 #include <list>
0049 #include <string>
0050
0051 #include "G4INCLParticle.hh"
0052 #include "G4INCLEventInfo.hh"
0053 #include "G4INCLCluster.hh"
0054 #include "G4INCLFinalState.hh"
0055 #include "G4INCLStore.hh"
0056 #include "G4INCLGlobals.hh"
0057 #include "G4INCLParticleTable.hh"
0058 #include "G4INCLConfig.hh"
0059 #include "G4INCLConfigEnums.hh"
0060 #include "G4INCLCluster.hh"
0061 #include "G4INCLProjectileRemnant.hh"
0062
0063 namespace G4INCL {
0064
0065   enum AnnihilationType {Def=0, PType, NType, PTypeInFlight, NTypeInFlight, NbarPTypeInFlight, NbarNTypeInFlight};
0066
0067   class Nucleus : public Cluster {
0068   public:
0069     Nucleus(G4int mass, G4int charge, G4int strangess, Config const * const conf, const G4double universeRadius=-1., AnnihilationType AType=Def);
0070     virtual ~Nucleus();
0071
0072     /// \brief Dummy copy constructor to silence Coverity warning
0073     Nucleus(const Nucleus &rhs);
0074
0075     /// \brief Dummy assignment operator to silence Coverity warning
0076     Nucleus &operator=(const Nucleus &rhs);
0077
0078     AnnihilationType getAType() const;
0079     void setAType(AnnihilationType type);
0080
0081     /**
0082      * Call the Cluster method to generate the initial distribution of
0083      * particles. At the beginning all particles are assigned as spectators.
0084      */
0085     void initializeParticles();
0086
0087     /// \brief Insert a new particle (e.g. a projectile) in the nucleus.
0088     void insertParticle(Particle *p) {
0089       theZ += p->getZ();
0090       theA += p->getA();
0091       theS += p->getS();
0092       theStore->particleHasEntered(p);
0093       if(p->isNucleon()) {
0094         theNpInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
0095         theNnInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
0096       }
0097       if(p->isPion()) {
0098         theNpionplusInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
0099         theNpionminusInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
0100       }
0101       if(p->isKaon() || p->isAntiKaon()) {
0102         theNkaonplusInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
0103         theNkaonminusInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
0104       }
0105       if(p->isAntiNucleon()) {
0106         theNantiprotonInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
0107       }
0108       if(!p->isTargetSpectator()) theStore->getBook().incrementCascading();
0109     };
0110
0111     /**
0112      * Apply reaction final state information to the nucleus.
0113      */
0114     void applyFinalState(FinalState *);
0115
0116     G4int getInitialA() const { return theInitialA; };
0117     G4int getInitialZ() const { return theInitialZ; };
0118     G4int getInitialS() const { return theInitialS; };
0119
0120     /**
0121      * Propagate the particles one time step.
0122      *
0123      * @param step length of the time step
0124      */
0125     void propagateParticles(G4double step);
0126
0127     G4int getNumberOfEnteringProtons() const { return theNpInitial; };
0128     G4int getNumberOfEnteringNeutrons() const { return theNnInitial; };
0129     G4int getNumberOfEnteringPions() const { return theNpionplusInitial+theNpionminusInitial; };
0130     G4int getNumberOfEnteringKaons() const { return theNkaonplusInitial+theNkaonminusInitial; };
0131     G4int getNumberOfEnteringantiProtons() const { return theNantiprotonInitial; };
0132
0133     /** \brief Outgoing - incoming separation energies.
0134      *
0135      * Used by CDPP.
0136      */
0137     G4double computeSeparationEnergyBalance() const {
0138       G4double S = 0.0;
0139       ParticleList const &outgoing = theStore->getOutgoingParticles();
0140       for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
0141         const ParticleType t = (*i)->getType();
0142         switch(t) {
0143           case Proton:
0144           case Neutron:
0145           case DeltaPlusPlus:
0146           case DeltaPlus:
0147           case DeltaZero:
0148           case DeltaMinus:
0149           case Lambda:
0150           case PiPlus:
0151           case PiMinus:
0152           case KPlus:
0153           case KMinus:
0154           case KZero:
0155           case KZeroBar:
0156           case KShort:
0157           case KLong:
0158           case SigmaPlus:
0159           case SigmaZero:
0160           case SigmaMinus:
0161           case antiProton:
0162           //case antiNeutron:
0163           //case antiLambda:
0164             S += thePotential->getSeparationEnergy(*i);
0165             break;
0166           case Composite:
0167             S += (*i)->getZ() * thePotential->getSeparationEnergy(Proton)
0168               + ((*i)->getA() + (*i)->getS() - (*i)->getZ()) * thePotential->getSeparationEnergy(Neutron)
0169               - (*i)->getS() * thePotential->getSeparationEnergy(Lambda);
0170             break;
0171           default:
0172             break;
0173         }
0174       }
0175
0176       S -= theNpInitial * thePotential->getSeparationEnergy(Proton);
0177       S -= theNnInitial * thePotential->getSeparationEnergy(Neutron);
0178       S -= theNpionplusInitial*thePotential->getSeparationEnergy(PiPlus);;
0179       S -= theNkaonplusInitial*thePotential->getSeparationEnergy(KPlus);
0180       S -= theNpionminusInitial*thePotential->getSeparationEnergy(PiMinus);
0181       S -= theNkaonminusInitial*thePotential->getSeparationEnergy(KMinus);
0182       S -= theNantiprotonInitial*thePotential->getSeparationEnergy(antiProton);
0183       return S;
0184     }
0185
0186     /** \brief Force the decay of outgoing deltas.
0187      *
0188      * \return true if any delta was forced to decay.
0189      */
0190     G4bool decayOutgoingDeltas();
0191
0192     /** \brief Force the decay of deltas inside the nucleus.
0193      *
0194      * \return true if any delta was forced to decay.
0195      */
0196     G4bool decayInsideDeltas();
0197
0198     /** \brief Force the transformation of strange particles into a Lambda;
0199      *
0200      *  \return true if any strange particles was forced to absorb.
0201      */
0202     G4bool decayInsideStrangeParticles();
0203
0204     /** \brief Force the decay of outgoing PionResonances (eta/omega).
0205      *
0206      * \return true if any eta was forced to decay.
0207      */
0208     G4bool decayOutgoingPionResonances(G4double timeThreshold);
0209
0210     /** \brief Force the decay of outgoing Neutral Sigma.
0211      *
0212      * \return true if any Sigma was forced to decay.
0213      */
0214     G4bool decayOutgoingSigmaZero(G4double timeThreshold);
0215
0216     /** \brief Force the transformation of outgoing Neutral Kaon into propation eigenstate.
0217      *
0218      * \return true if any kaon was forced to decay.
0219      */
0220     G4bool decayOutgoingNeutralKaon();
0221
0222     /** \brief Force the decay of unstable outgoing clusters.
0223      *
0224      * \return true if any cluster was forced to decay.
0225      */
0226     G4bool decayOutgoingClusters();
0227
0228     /** \brief Force the phase-space decay of the Nucleus.
0229      *
0230      * Only applied if Z==0 or N==0.
0231      *
0232      * \return true if the nucleus was forced to decay.
0233      */
0234     G4bool decayMe();
0235
0236     /// \brief Force emission of all pions inside the nucleus.
0237     void emitInsidePions();
0238
0239     /// \brief Force emission of all strange particles inside the nucleus.
0240     void emitInsideStrangeParticles();
0241
0242     /// \brief Force emission of all Lambda (desexitation code with strangeness not implanted yet)
0243     G4int emitInsideLambda();
0244
0245     /// \brief Force emission of all Kaon inside the nucleus
0246     G4bool emitInsideKaon();
0247
0248     /** \brief Compute the recoil momentum and spin of the nucleus. */
0249     void computeRecoilKinematics();
0250
0251     /** \brief Compute the current center-of-mass position.
0252      *
0253      * \return the center-of-mass position vector [fm].
0254      */
0255     ThreeVector computeCenterOfMass() const;
0256
0257     /** \brief Compute the current total energy.
0258      *
0259      * \return the total energy [MeV]
0260      */
0261     G4double computeTotalEnergy() const;
0262
0263     /** \brief Compute the current excitation energy.
0264      *
0265      * \return the excitation energy [MeV]
0266      */
0267     G4double computeExcitationEnergy() const;
0268
0269     /** \brief Set the incoming angular-momentum vector. */
0270     void setIncomingAngularMomentum(const ThreeVector &j) {
0271       incomingAngularMomentum = j;
0272     }
0273
0274     /** \brief Get the incoming angular-momentum vector. */
0275     const ThreeVector &getIncomingAngularMomentum() const { return incomingAngularMomentum; }
0276
0277     /** \brief Set the incoming momentum vector. */
0278     void setIncomingMomentum(const ThreeVector &p) {
0279       incomingMomentum = p;
0280     }
0281
0282     /** \brief Get the incoming momentum vector. */
0283     const ThreeVector &getIncomingMomentum() const {
0284       return incomingMomentum;
0285     }
0286
0287     /** \brief Set the initial energy. */
0288     void setInitialEnergy(const G4double e) { initialEnergy = e; }
0289
0290     /** \brief Get the initial energy. */
0291     G4double getInitialEnergy() const { return initialEnergy; }
0292
0293     /** \brief Get the excitation energy of the nucleus.
0294      *
0295      * Method computeRecoilKinematics() should be called first.
0296      */
0297     G4double getExcitationEnergy() const { return theExcitationEnergy; }
0298
0299     ///\brief Returns true if the nucleus contains any deltas.
0300     inline G4bool containsDeltas() {
0301       ParticleList const &inside = theStore->getParticles();
0302       for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0303         if((*i)->isDelta()) return true;
0304       return false;
0305     }
0306
0307     ///\brief Returns true if the nucleus contains any anti Kaons.
0308     inline G4bool containsAntiKaon() {
0309       ParticleList const &inside = theStore->getParticles();
0310       for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0311         if((*i)->isAntiKaon()) return true;
0312       return false;
0313     }
0314
0315     ///\brief Returns true if the nucleus contains any Lambda.
0316     inline G4bool containsLambda() {
0317       ParticleList const &inside = theStore->getParticles();
0318       for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0319         if((*i)->isLambda()) return true;
0320       return false;
0321     }
0322
0323     ///\brief Returns true if the nucleus contains any Sigma.
0324     inline G4bool containsSigma() {
0325       ParticleList const &inside = theStore->getParticles();
0326       for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0327         if((*i)->isSigma()) return true;
0328       return false;
0329     }
0330
0331     ///\brief Returns true if the nucleus contains any Kaons.
0332     inline G4bool containsKaon() {
0333       ParticleList const &inside = theStore->getParticles();
0334       for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0335         if((*i)->isKaon()) return true;
0336       return false;
0337     }
0338
0339       ///\brief Returns true if the nucleus contains any etas.
0340       inline G4bool containsEtas() {
0341           ParticleList const &inside = theStore->getParticles();
0342           for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0343               if((*i)->isEta()) return true;
0344           return false;
0345       }
0346
0347       ///\brief Returns true if the nucleus contains any omegas.
0348       inline G4bool containsOmegas() {
0349           ParticleList const &inside = theStore->getParticles();
0350           for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
0351               if((*i)->isOmega()) return true;
0352           return false;
0353       }
0354
0355
0356     /**
0357      * Print the nucleus info
0358      */
0359     std::string print();
0360
0361     Store* getStore() const {return theStore; };
0362     void setStore(Store *str) {
0363       delete theStore;
0364       theStore = str;
0365     };
0366
0367     G4double getInitialInternalEnergy() const { return initialInternalEnergy; };
0368
0369     /** \brief Is the event transparent?
0370      *
0371      * To be called at the end of the cascade.
0372      **/
0373     G4bool isEventTransparent() const;
0374
0375     /** \brief Does the nucleus give a cascade remnant?
0376      *
0377      * To be called after computeRecoilKinematics().
0378      **/
0379     G4bool hasRemnant() const { return remnant; }
0380
0381     /**
0382      * Fill the event info which contains INCL output data
0383      */
0384     void fillEventInfo(EventInfo *eventInfo);
0385
0386     G4bool getTryCompoundNucleus() { return tryCN; }
0387
0388     /// \brief Get the transmission barrier
0389     G4double getTransmissionBarrier(Particle const * const p) {
0390       const G4double theTransmissionRadius = theDensity->getTransmissionRadius(p);
0391       const G4double theParticleZ = p->getZ();
0392       return PhysicalConstants::eSquared*(theZ-theParticleZ)*theParticleZ/theTransmissionRadius;
0393     }
0394
0395     /// \brief Struct for conservation laws
0396     struct ConservationBalance {
0397       ThreeVector momentum;
0398       G4double energy;
0399       G4int Z, A, S;
0400     };
0401
0402     /// \brief Compute charge, mass, energy and momentum balance
0403     ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const;
0404
0405     /// \brief Adjust the kinematics for complete-fusion events
0406     void useFusionKinematics();
0407
0408     /** \brief Get the maximum allowed radius for a given particle.
0409      *
0410      * Calls the NuclearDensity::getMaxRFromP() method for nucleons and deltas,
0411      * and the NuclearDensity::getTrasmissionRadius() method for pions.
0412      *
0413      * \param particle pointer to a particle
0414      * \return surface radius
0415      */
0416     G4double getSurfaceRadius(Particle const * const particle) const {
0417       if(particle->isNucleon() || particle->isLambda() || particle->isResonance()){
0418         const G4double pr = particle->getReflectionMomentum()/thePotential->getFermiMomentum(particle);
0419         if(pr>=1.)
0420           return getUniverseRadius();
0421         else
0422           return theDensity->getMaxRFromP(particle->getType(), pr);
0423         }
0424       else {
0425         // Temporarily set RPION = RMAX
0426         return getUniverseRadius();
0427         //return 0.5*(theDensity->getTransmissionRadius(particle)+getUniverseRadius());
0428       }
0429     }
0430
0431     /// \brief Getter for theUniverseRadius.
0432     G4double getUniverseRadius() const { return theUniverseRadius; }
0433
0434     /// \brief Setter for theUniverseRadius.
0435     void setUniverseRadius(const G4double universeRadius) { theUniverseRadius=universeRadius; }
0436
0437     /// \brief Is it a nucleus-nucleus collision?
0438     G4bool isNucleusNucleusCollision() const { return isNucleusNucleus; }
0439
0440     /// \brief Set a nucleus-nucleus collision
0441     void setNucleusNucleusCollision() { isNucleusNucleus=true; }
0442
0443     /// \brief Set a particle-nucleus collision
0444     void setParticleNucleusCollision() { isNucleusNucleus=false; }
0445
0446     /// \brief Set the projectile remnant
0447     void setProjectileRemnant(ProjectileRemnant * const c) {
0448       delete theProjectileRemnant;
0449       theProjectileRemnant = c;
0450     }
0451
0452     /// \brief Get the projectile remnant
0453     ProjectileRemnant *getProjectileRemnant() const { return theProjectileRemnant; }
0454
0455     /// \brief Delete the projectile remnant
0456     void deleteProjectileRemnant() {
0457       delete theProjectileRemnant;
0458       theProjectileRemnant = NULL;
0459     }
0460
0461     /** \brief Finalise the projectile remnant
0462      *
0463      * Complete the treatment of the projectile remnant. If it contains
0464      * nucleons, assign its excitation energy and spin. Move stuff to the
0465      * outgoing list, if appropriate.
0466      *
0467      * \param emissionTime the emission time of the projectile remnant
0468      */
0469     void finalizeProjectileRemnant(const G4double emissionTime);
0470
0471     /// \brief Update the particle potential energy.
0472     inline void updatePotentialEnergy(Particle *p) const {
0473       p->setPotentialEnergy(thePotential->computePotentialEnergy(p));
0474     }
0475
0476     /// \brief Setter for theDensity
0477     void setDensity(NuclearDensity const * const d) {
0478       theDensity=d;
0479       if(theParticleSampler)
0480         theParticleSampler->setDensity(theDensity);
0481     };
0482
0483     /// \brief Getter for theDensity
0484     NuclearDensity const *getDensity() const { return theDensity; };
0485
0486     /// \brief Getter for thePotential
0487     NuclearPotential::INuclearPotential const *getPotential() const { return thePotential; };
0488
0489     /// \brief Getter for theAnnihilationType
0490     AnnihilationType getAnnihilationType() const { return theAType; }; //D
0491
0492     /// \brief Setter for theAnnihilationType
0493     void setAnnihilationType(const AnnihilationType at){
0494       theAType = at;
0495     }; //D
0496
0497   private:
0498     /** \brief Compute the recoil kinematics for a 1-nucleon remnant.
0499      *
0500      * Puts the remnant nucleon on mass shell and tries to enforce approximate
0501      * energy conservation by modifying the masses of the outgoing particles.
0502      */
0503     void computeOneNucleonRecoilKinematics();
0504
0505   private:
0506
0507     G4int theInitialZ, theInitialA, theInitialS;
0508     /// \brief The number of entering protons
0509     G4int theNpInitial;
0510     /// \brief The number of entering neutrons
0511     G4int theNnInitial;
0512     /// \brief The number of entering pions
0513     G4int theNpionplusInitial;
0514     G4int theNpionminusInitial;
0515     /// \brief The number of entering kaons
0516     G4int theNkaonplusInitial;
0517     G4int theNkaonminusInitial;
0518     /// \brief The number of entering antiprotons
0519     G4int theNantiprotonInitial;
0520
0521     G4double initialInternalEnergy;
0522     ThreeVector incomingAngularMomentum, incomingMomentum;
0523     ThreeVector initialCenterOfMass;
0524     G4bool remnant;
0525
0526     G4double initialEnergy;
0527     Store *theStore;
0528     G4bool tryCN;
0529
0530     /// \brief The radius of the universe
0531     G4double theUniverseRadius;
0532
0533     /** \brief true if running a nucleus-nucleus collision
0534      *
0535      * Tells INCL whether to make a projectile-like pre-fragment or not.
0536      */
0537     G4bool isNucleusNucleus;
0538
0539     /** \brief Pointer to the quasi-projectile
0540      *
0541      * Owned by the Nucleus object.
0542      */
0543     ProjectileRemnant *theProjectileRemnant;
0544
0545     /// \brief Pointer to the NuclearDensity object
0546     NuclearDensity const *theDensity;
0547
0548     /// \brief Pointer to the NuclearPotential object
0549     NuclearPotential::INuclearPotential const *thePotential;
0550
0551     AnnihilationType theAType; //D same order as in the cc
0552
0553     INCL_DECLARE_ALLOCATION_POOL(Nucleus)
0554   };
0555
0556 }
0557
0558 #endif /* G4INCLNUCLEUS_HH_ */