include/Geant4/G4INCLCascade.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 #ifndef G4INCLCascade_hh
0039 #define G4INCLCascade_hh 1
0040
0041 #include "G4INCLParticle.hh"
0042 #include "G4INCLNucleus.hh"
0043 #include "G4INCLIPropagationModel.hh"
0044 #include "G4INCLCascadeAction.hh"
0045 #include "G4INCLEventInfo.hh"
0046 #include "G4INCLGlobalInfo.hh"
0047 #include "G4INCLLogger.hh"
0048 #include "G4INCLConfig.hh"
0049 #include "G4INCLRootFinder.hh"
0050 #ifndef INCLXX_IN_GEANT4_MODE
0051  #include "G4INCLGeant4Compat.hh"
0052 #endif
0053
0054
0055 namespace G4INCL {
0056   class INCL {
0057     public:
0058       INCL(Config const * const config);
0059
0060       ~INCL();
0061
0062       /// \brief Dummy copy constructor to silence Coverity warning
0063       INCL(const INCL &rhs);
0064
0065       /// \brief Dummy assignment operator to silence Coverity warning
0066       INCL &operator=(const INCL &rhs);
0067
0068       G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S);
0069
0070       G4bool initializeTarget(const G4int A, const G4int Z, const G4int S, AnnihilationType theAType);
0071       G4double initUniverseRadiusForAntiprotonAtRest(const G4int A, const G4int Z);
0072       inline const EventInfo &processEvent() {
0073         return processEvent(
0074             theConfig->getProjectileSpecies(),
0075             theConfig->getProjectileKineticEnergy(),
0076             theConfig->getTargetA(),
0077             theConfig->getTargetZ(),
0078             theConfig->getTargetS()
0079             );
0080       }
0081       const EventInfo &processEvent(
0082           ParticleSpecies const &projectileSpecies,
0083           const G4double kineticEnergy,
0084           const G4int targetA,
0085           const G4int targetZ,
0086           const G4int targetS
0087           );
0088
0089       void finalizeGlobalInfo(Random::SeedVector const &initialSeeds);
0090       const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
0091
0092
0093     private:
0094       IPropagationModel *propagationModel;
0095       G4int theA, theZ, theS;
0096       G4bool targetInitSuccess;
0097       G4double maxImpactParameter;
0098       G4double maxUniverseRadius;
0099       G4double maxInteractionDistance;
0100       G4double fixedImpactParameter;
0101       CascadeAction *cascadeAction;
0102       Config const * const theConfig;
0103       Nucleus *nucleus;
0104       G4bool forceTransparent;
0105
0106       EventInfo theEventInfo;
0107       GlobalInfo theGlobalInfo;
0108
0109       /// \brief Remnant size below which cascade stops
0110       G4int minRemnantSize;
0111
0112       /// \brief Class to adjust remnant recoil
0113       class RecoilFunctor : public RootFunctor {
0114         public:
0115           /** \brief Prepare for calling the () operator and scaleParticleEnergies
0116            *
0117            * The constructor sets the private class members.
0118            */
0119           RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
0120             RootFunctor(0., 1E6),
0121             nucleus(n),
0122             outgoingParticles(n->getStore()->getOutgoingParticles()),
0123             theEventInfo(ei) {
0124               for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
0125                 particleMomenta.push_back((*p)->getMomentum());
0126                 particleKineticEnergies.push_back((*p)->getKineticEnergy());
0127               }
0128               ProjectileRemnant * const aPR = n->getProjectileRemnant();
0129               if(aPR && aPR->getA()>0) {
0130                 particleMomenta.push_back(aPR->getMomentum());
0131                 particleKineticEnergies.push_back(aPR->getKineticEnergy());
0132                 outgoingParticles.push_back(aPR);
0133               }
0134             }
0135           virtual ~RecoilFunctor() {}
0136
0137           /** \brief Compute the energy-conservation violation.
0138            *
0139            * \param x scale factor for the particle energies
0140            * \return the energy-conservation violation
0141            */
0142           G4double operator()(const G4double x) const {
0143             scaleParticleEnergies(x);
0144             return nucleus->getConservationBalance(theEventInfo,true).energy;
0145           }
0146
0147           /// \brief Clean up after root finding
0148           void cleanUp(const G4bool success) const {
0149             if(!success)
0150               scaleParticleEnergies(1.);
0151           }
0152
0153         private:
0154           /// \brief Pointer to the nucleus
0155           Nucleus *nucleus;
0156           /// \brief List of final-state particles.
0157           ParticleList outgoingParticles;
0158           // \brief Reference to the EventInfo object
0159           EventInfo const &theEventInfo;
0160           /// \brief Initial momenta of the outgoing particles
0161           std::list<ThreeVector> particleMomenta;
0162           /// \brief Initial kinetic energies of the outgoing particles
0163           std::list<G4double> particleKineticEnergies;
0164
0165           /** \brief Scale the kinetic energies of the outgoing particles.
0166            *
0167            * \param rescale scale factor
0168            */
0169           void scaleParticleEnergies(const G4double rescale) const {
0170             // Rescale the energies (and the momenta) of the outgoing particles.
0171             ThreeVector pBalance = nucleus->getIncomingMomentum();
0172             std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
0173             std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
0174             for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP, ++iE)
0175             {
0176               const G4double mass = (*i)->getMass();
0177               const G4double newKineticEnergy = (*iE) * rescale;
0178
0179               (*i)->setMomentum(*iP);
0180               (*i)->setEnergy(mass + newKineticEnergy);
0181               (*i)->adjustMomentumFromEnergy();
0182
0183               pBalance -= (*i)->getMomentum();
0184             }
0185
0186             nucleus->setMomentum(pBalance);
0187             const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ(),nucleus->getS()) + nucleus->getExcitationEnergy();
0188             const G4double pRem2 = pBalance.mag2();
0189             const G4double recoilEnergy = pRem2/
0190               (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
0191             nucleus->setEnergy(remnantMass + recoilEnergy);
0192           }
0193       };
0194
0195       /// \brief Class to adjust remnant recoil in the reaction CM system
0196       class RecoilCMFunctor : public RootFunctor {
0197         public:
0198           /** \brief Prepare for calling the () operator and scaleParticleEnergies
0199            *
0200            * The constructor sets the private class members.
0201            */
0202           RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
0203             RootFunctor(0., 1E6),
0204             nucleus(n),
0205             theIncomingMomentum(nucleus->getIncomingMomentum()),
0206             outgoingParticles(n->getStore()->getOutgoingParticles()),
0207             theEventInfo(ei) {
0208               if(theIncomingMomentum.mag() == 0.){ //change the condition
0209                 thePTBoostVector = {0., 0., 0.};
0210                 //INCL_WARN("PTBoostVector at rest is zero" << '\n');
0211               }
0212               else{
0213                 thePTBoostVector = nucleus->getIncomingMomentum()/(nucleus->getInitialEnergy()); //D
0214                 //INCL_WARN("PTBoostVector" << '\n');
0215               }
0216               for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
0217                 (*p)->boost(thePTBoostVector);
0218                 particleCMMomenta.push_back((*p)->getMomentum());
0219               }
0220               ProjectileRemnant * const aPR = n->getProjectileRemnant();
0221               if(aPR && aPR->getA()>0) {
0222                 aPR->boost(thePTBoostVector);
0223                 particleCMMomenta.push_back(aPR->getMomentum());
0224                 outgoingParticles.push_back(aPR);
0225               }
0226             }
0227           virtual ~RecoilCMFunctor() {}
0228
0229           /** \brief Compute the energy-conservation violation.
0230            *
0231            * \param x scale factor for the particle energies
0232            * \return the energy-conservation violation
0233            */
0234           G4double operator()(const G4double x) const {
0235             scaleParticleCMMomenta(x);
0236             return nucleus->getConservationBalance(theEventInfo,true).energy;
0237           }
0238
0239           /// \brief Clean up after root finding
0240           void cleanUp(const G4bool success) const {
0241             if(!success)
0242               scaleParticleCMMomenta(1.);
0243           }
0244
0245         private:
0246           /// \brief Pointer to the nucleus
0247           Nucleus *nucleus;
0248           /// \brief Projectile-target CM boost vector
0249           ThreeVector thePTBoostVector;
0250           /// \brief Incoming momentum
0251           ThreeVector theIncomingMomentum;
0252           /// \brief List of final-state particles.
0253           ParticleList outgoingParticles;
0254           /// \brief Reference to the EventInfo object
0255           EventInfo const &theEventInfo;
0256           /// \brief Initial CM momenta of the outgoing particles
0257           std::list<ThreeVector> particleCMMomenta;
0258
0259           /** \brief Scale the kinetic energies of the outgoing particles.
0260            *
0261            * \param rescale scale factor
0262            */
0263           void scaleParticleCMMomenta(const G4double rescale) const {
0264             // Rescale the CM momenta of the outgoing particles.
0265             ThreeVector remnantMomentum = theIncomingMomentum;
0266             std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
0267             for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP)
0268             {
0269               (*i)->setMomentum(*iP * rescale);
0270               (*i)->adjustEnergyFromMomentum();
0271               (*i)->boost(-thePTBoostVector);
0272
0273               remnantMomentum -= (*i)->getMomentum();
0274             }
0275
0276             nucleus->setMomentum(remnantMomentum);
0277             const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ(),nucleus->getS()) + nucleus->getExcitationEnergy();
0278             const G4double pRem2 = remnantMomentum.mag2();
0279             const G4double recoilEnergy = pRem2/
0280               (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
0281             nucleus->setEnergy(remnantMass + recoilEnergy);
0282           }
0283       };
0284
0285       /** \brief Rescale the energies of the outgoing particles.
0286        *
0287        * Allow for the remnant recoil energy by rescaling the energy (and
0288        * momenta) of the outgoing particles.
0289        */
0290       void rescaleOutgoingForRecoil();
0291
0292 #ifndef INCLXX_IN_GEANT4_MODE
0293       /** \brief Run global conservation checks
0294        *
0295        * Check that energy and momentum are correctly conserved. If not, issue
0296        * a warning.
0297        *
0298        * Also feeds the balance variables in theEventInfo.
0299        *
0300        * \param afterRecoil whether to take into account nuclear recoil
0301        */
0302       void globalConservationChecks(G4bool afterRecoil);
0303 #endif
0304
0305       /** \brief Stopping criterion for the cascade
0306        *
0307        * Returns true if the cascade should continue, and false if any of the
0308        * stopping criteria is satisfied.
0309        */
0310       G4bool continueCascade();
0311
0312       /** \brief Make a projectile pre-fragment out of geometrical spectators
0313        *
0314        * The projectile pre-fragment is assigned an excitation energy given
0315        * by \f$E_\mathrm{sp}-E_\mathrm{i,A}\f$, where \f$E_\mathrm{sp}\f$ is the
0316        * sum of the energies of the spectator particles, and \f$E_\mathrm{i,A}\f$
0317        * is the sum of the smallest \f$A\f$ particle energies initially present
0318        * in the projectile, \f$A\f$ being the mass of the projectile
0319        * pre-fragment. This is equivalent to assuming that the excitation
0320        * energy is given by the sum of the transitions of all excited
0321        * projectile components to the "holes" left by the participants.
0322        *
0323        * This method can modify the outgoing list and adds a projectile
0324        * pre-fragment.
0325        *
0326        * \return the number of dynamical spectators that were merged back in
0327        *         the projectile
0328        */
0329       G4int makeProjectileRemnant();
0330
0331       /** \brief Make a compound nucleus
0332        *
0333        * Selects the projectile components that can actually enter their
0334        * potential and puts them into the target nucleus. If the CN excitation
0335        * energy turns out to be negative, the event is considered a
0336        * transparent. This method modifies theEventInfo and theGlobalInfo.
0337        */
0338       void makeCompoundNucleus();
0339
0340       /// \brief Initialise the cascade
0341       G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
0342
0343       /// \brief The actual cascade loop
0344       void cascade();
0345
0346       /// \brief Finalise the cascade and clean up
0347       void postCascade(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy);
0348
0349       /** \brief Initialise the maximum interaction distance.
0350        *
0351        * Used in forced CN events.
0352        */
0353       void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
0354
0355       /** \brief Initialize the universe radius.
0356        *
0357        * Used for determining the energy-dependent size of the volume particles
0358        * live in.
0359        */
0360       void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
0361
0362       /// \brief Update global counters and other members of theGlobalInfo object
0363       void updateGlobalInfo();
0364
0365       /// \brief Fill probabilities and particle types from datafile and return probability sum for normalization
0366       G4double read_file(std::string filename, std::vector<G4double>& probabilities,
0367                          std::vector<std::vector<G4String>>& particle_types);
0368
0369       /// \brief This function will tell you the FS line number from the datafile based on your random value
0370       G4int findStringNumber(G4double rdm, std::vector<G4double> yields);
0371
0372       /// \brief Initialise the "cascade" for pbar on H1
0373       void preCascade_pbarH1(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
0374
0375       /// \brief Initialise the "cascade" for pbar on H2
0376       void preCascade_pbarH2(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
0377
0378       /// \brief Finalise the "cascade" and clean up for pbar on H1
0379       void postCascade_pbarH1(ParticleList const &outgoingParticles);
0380
0381       /// \brief Finalise the "cascade" and clean up for pbar on H2
0382       void postCascade_pbarH2(ParticleList const &outgoingParticles, ParticleList const &H2Particles);
0383   };
0384 }
0385
0386 #endif