EIC code displayed by LXR master/​include/​Pythia8Plugins/​PythiaCascade.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-09-15 09:06:57

 // PythiaCascade.h is a part of the PYTHIA event generator.
 // Copyright (C) 2025 Torbjorn Sjostrand.
 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
 // Please respect the MCnet Guidelines, see GUIDELINES for details.
 // Author: Torbjorn Sjostrand.
 
 #ifndef Pythia8_PythiaCascade_H
 #define Pythia8_PythiaCascade_H
 
 #include "Pythia8/Pythia.h"
 
 namespace Pythia8 {
 
 //==========================================================================
 
 // Wrapper class for Pythia usage in cosmic ray or detector cascades.
 // Here it is assumed that the user wants to control the full
 // evolution.  The complete event record is then kept somewhere
 // separate from PYTHIA, and specific input on production and decay
 // generates the next subevent.
 
 // Intended flow:
 
 // - init sets up all generation, given a maximum energy. This may be
 //   time-consuming, less so if MPI initialization data can be reused.
 
 // - sigmaSetuphN prepares a possible collision for a given hadron by
 //   calculating the hadron-nucleon cross section (if possible).
 //   Moderately time-consuming.
 
 // - sigmahA uses the hN cross section calculated by sigaSetuphN and
 //   returns hadron-ion cross section for a collision with a specified
 //   nucleus.  It can be called several times to cover a mix of target
 //   nuclei, with minimal time usage.
 
 // - nextColl performs the hadron-nucleus collision, as a sequences of
 //   hadron-nucleon ones. Can be quite time-consuming.
 
 // - nextDecay can be used anytime to decay a particle. Each
 //   individual decay is rather fast, but there may be many of them.
 
 // - stat can be used at end of run to give a summary of error
 //   messages.  Negligible time usage.
 
 // - references to particleData() and rndm() can be used in the main
 // - program.
 
 // Note: when a hadron interacts with a medium particle, the latter is
 // added to the event record. This program uses these additional
 // status codes for target particles in the medium:
 //  181: the first (or only) target nucleon in a collision.
 //  182: nucleons from further subcollisions with the same target nucleus.
 
 // Warning: the large boosts involved at the higher cosmic ray
 // energies are not perfectly managed, which leads to non-negligible
 // energy-momentum non-conservation. The worst (sub)events are caught
 // and regenerated.
 
 //--------------------------------------------------------------------------
 
 class PythiaCascade {
 
 public:
 
   // Default constructor, all setup is done in init().
   PythiaCascade() = default;
 
   //--------------------------------------------------------------------------
 
   // Initialize PythiaCascade for a given maximal incoming energy.
 
   // Set eMax to the maximal incoming energy (in GeV) on fixed target.
 
   // Keep listFinal = false if you want to get back the full event
   // record, else only the "final" particles of the collision or decay
   // are returned.
 
   // Keep rapidDecays = false if you want to do every decay yourself,
   // else all particle types with a tau0 below smallTau0 will be
   // decayed immediately.  The tau0 units are mm/c, so default gives c
   // * tau0 = 1e-10 mm = 100 fm.  Note that time dilation effects are
   // not included, and they can be large.
 
   // The reuseMPI argument mimics MultipartonInteractions:reuseInit,
   // but streamlined into three main options;
   // 0 = current run is self-contained, so MPI is initialized from scratch,
   //     but MPI data is stored on initFile for future use;
   // 3 = read MPI initialization data from initFile, if it exists,
   //     and else initialize and save the MPI data on initFile;
   // -1 = read MPI initialization data and other run settings from
   //    a .cmnd file, if not empty.
   // In each case, the saved initialization must have been done with
   // the same or a larger eMax than the one now being used.
 
   bool init(double eMaxIn = 1e9, bool listFinalIn = false,
     bool rapidDecaysIn = false, double smallTau0In = 1e-10,
     int reuseMPI = 3, string initFile = "pythiaCascade.mpi") {
 
     // Store input for future usage.
     eMax        = eMaxIn;
     listFinal   = listFinalIn;
     rapidDecays = rapidDecaysIn;
     smallTau0   = smallTau0In;
 
     // Proton mass.
     mp      = pythiaMain.particleData.m0(2212);
 
     // Main Pythia object for managing the cascade evolution in a
     // nucleus. Can also do decays, but no hard processes.
     pythiaMain.readString("ProcessLevel:all = off");
     pythiaMain.readString("13:mayDecay  = on");
     pythiaMain.readString("211:mayDecay = on");
     pythiaMain.readString("321:mayDecay = on");
     pythiaMain.readString("130:mayDecay = on");
     pythiaMain.settings.flag("ParticleDecays:limitTau0", rapidDecays);
     pythiaMain.settings.parm("ParticleDecays:tau0Max", smallTau0);
 
     // Redure statistics printout to relevant ones.
     pythiaMain.readString("Stat:showProcessLevel = off");
     pythiaMain.readString("Stat:showPartonLevel = off");
 
     // Initialize. Return if failure.
     if (!pythiaMain.init()) return false;
 
     if ( reuseMPI < 0 ) {
       pythiaColl.readFile(initFile);
       initFile = "";
     }
 
     // Secondary Pythia object for performing individual collisions,
     // or decays. Variable incoming beam type and energy.
     pythiaColl.readString("Beams:allowVariableEnergy = on");
     pythiaColl.readString("Beams:allowIDAswitch = on");
 
     // Set up for fixed-target collisions.
     pythiaColl.readString("Beams:frameType = 3");
     pythiaColl.settings.parm("Beams:pzA", eMax);
     pythiaColl.readString("Beams:pzB = 0.");
 
     // Must use the soft and low-energy QCD processes.
     pythiaColl.readString("SoftQCD:all = on");
     pythiaColl.readString("LowEnergyQCD:all = on");
 
     // Primary (single) decay to be done by pythiaColl, to circumvent
     // limitTau0.
     pythiaColl.readString("13:mayDecay  = on");
     pythiaColl.readString("211:mayDecay = on");
     pythiaColl.readString("321:mayDecay = on");
     pythiaColl.readString("130:mayDecay = on");
 
     // Secondary decays to be done by pythiaMain, respecting limitTau0.
     pythiaColl.readString("HadronLevel:Decay = off");
 
     // Reduce printout and relax energy-momentum conservation.
     // (Unusually large errors unfortunate consequence of large boosts.)
     pythiaColl.readString("Print:quiet = on");
     pythiaColl.readString("Check:epTolErr = 0.01");
     pythiaColl.readString("Check:epTolWarn = 0.0001");
     pythiaColl.readString("Check:mTolErr = 0.01");
 
     // Redure statistics printout to relevant ones.
     pythiaColl.readString("Stat:showProcessLevel = off");
     pythiaColl.readString("Stat:showPartonLevel = off");
 
     // Reuse MPI initialization file if it exists; else create a new one.
     if (reuseMPI > 0)
       pythiaColl.readString("MultipartonInteractions:reuseInit = 3");
     else if (reuseMPI == 0)
       pythiaColl.readString("MultipartonInteractions:reuseInit = 1");
     if (reuseMPI >= 0)
       pythiaColl.settings.word("MultipartonInteractions:initFile", initFile);
 
     // Initialize.
     if (!pythiaColl.init()) return false;
     return true;
 
   }
 
   //--------------------------------------------------------------------------
 
   // Calculate the average number of collisions. Average number of
   // hadron-nucleon collisions in a hadron-nucleus one. If not in
   // table then interpolate, knowing that <n> - 1 propto A^{2/3}.
 
   double nCollAvg(int A) {
     for (int i = 0; i < nA; ++i) {
       if (A == tabA[i]) {
         return (sigmaNow < tabBorder) ? 1. + tabSlopeLo[i] * sigmaNow
           : 1. + tabOffset[i] + tabSlope[i] * sigmaNow;
       } else if (A < tabA[i]) {
         double nColl1 = (sigmaNow < tabBorder) ? tabSlopeLo[i - 1] * sigmaNow
           : tabOffset[i - 1] + tabSlope[i - 1] * sigmaNow;
         double nColl2 = (sigmaNow < tabBorder) ? tabSlopeLo[i] * sigmaNow
           : tabOffset[i] + tabSlope[i] * sigmaNow;
         double wt1 = (tabA[i] - A) / (tabA[i] - tabA[i - 1]);
         return 1. + wt1 * pow( A / tabA[i - 1], 2./3.) * nColl1
           + (1. - wt1) * pow( A / tabA[i], 2./3.) * nColl2;
       }
     }
 
     return numeric_limits<double>::quiet_NaN();
   }
 
   //--------------------------------------------------------------------------
 
   // Calculate the hadron-proton collision cross section.
   // Return false if not possible to find.
 
   bool sigmaSetuphN(int idNowIn, Vec4 pNowIn, double mNowIn) {
 
     // Cannot handle low-energy hadrons.
     if (pNowIn.e() - mNowIn < eKinMin) return false;
 
     // Cannot handle hadrons above maximum energy set at initialization.
     if (pNowIn.e() > eMax) {
       logger.ERROR_MSG("too high energy");
       return false;
     }
 
     // Save incoming quantities for reuse in later methods.
     idNow = idNowIn;
     pNow  = pNowIn;
     mNow  = mNowIn;
 
     // Calculate hadron-nucleon cross section. Check if cross section
     // vanishes.
     eCMNow = (pNow + Vec4(0, 0, 0, mp)).mCalc();
     sigmaNow = pythiaColl.getSigmaTotal(idNow, 2212, eCMNow, mNow, mp);
     if (sigmaNow <= 0.) {
       if (eCMNow - mNow - mp > eKinMin)
         logger.ERROR_MSG("vanishing cross section");
        return false;
     }
 
     // Done.
     return true;
 
   }
 
   //--------------------------------------------------------------------------
 
   // Calculate the hadron-nucleus cross section for a given nucleon
   // number A, using hN cross section from sigmaSetuphN. Interpolate
   // where not (yet) available.
 
   double sigmahA(int A) {
 
     // Restrict to allowed range 1 <= A <= 208.
     if (A < 1 || A > 208) {
       logger.ERROR_MSG("A is outside of valid range (1 <= A <= 208)");
       return 0.;
     }
 
     // Correction factor for number of h-nucleon collisions per
     // h-nucleus one.
     double sigmahA = A * sigmaNow / nCollAvg(A);
 
     // Done.
     return sigmahA;
 
   }
 
   //--------------------------------------------------------------------------
 
   // Generate a collision, and return the event record.
   // Input (Z, A) of nucleus, and optionally collision vertex.
 
   Event& nextColl(int Znow, int Anow, Vec4 vNow = Vec4() ) {
 
     // References to the two event records. Clear main event record.
     Event& eventMain = pythiaMain.event;
     Event& eventColl = pythiaColl.event;
     eventMain.clear();
 
     // Restrict to allowed range 1 <= A <= 208.
     if (Anow < 1 || Anow > 208) {
       logger.ERROR_MSG("A is outside of valid range (1 <= A <= 208)");
       return eventMain;
     }
 
     // Insert incoming particle in cleared main event record.
     eventMain.append(90,   -11, 0, 0, 1, 1, 0, 0, pNow, mNow);
     int iHad = eventMain.append(idNow, 12, 0, 0, 0, 0, 0, 0, pNow, mNow);
     eventMain[iHad].vProd(vNow);
 
     // Set up for collisions on a nucleus.
     int np      = Znow;
     int nn      = Anow - Znow;
     int sizeOld = 0;
     int sizeNew = 0;
     Vec4 dirNow = pNow / pNow.pAbs();
     Rndm& rndm  = pythiaMain.rndm;
 
     // Drop rate of geometric series. (Deuterium has corrected nCollAvg.)
     double probMore = 1. - 1. / nCollAvg(Anow);
 
     // Loop over varying number of hit nucleons in target nucleus.
     for (int iColl = 1; iColl <= Anow; ++iColl) {
       if (iColl > 1 && rndm.flat() > probMore) break;
 
       // Pick incoming projectile: trivial for first subcollision, else ...
       int iProj    = iHad;
       int procType = 0;
 
       // ... find highest-pLongitudinal particle from latest subcollision.
       if (iColl > 1) {
         iProj = 0;
         double pMax = 0.;
         for (int i = sizeOld; i < sizeNew; ++i)
         if ( eventMain[i].isFinal() && eventMain[i].isHadron()) {
           double pp = dot3(dirNow, eventMain[i].p());
           if (pp > pMax) {
             iProj = i;
             pMax  = pp;
           }
         }
 
         // No further subcollision if no particle with enough energy.
         if ( iProj == 0 || eventMain[iProj].e() - eventMain[iProj].m()
           < eKinMin) break;
 
         // Choose process; only SD or ND at perturbative energies.
         double eCMSub = (eventMain[iProj].p() + Vec4(0, 0, 0, mp)).mCalc();
         if (eCMSub > 10.) procType = (rndm.flat() < probSD) ? 4 : 1;
       }
 
       // Pick one p or n from target.
       int idProj = eventMain[iProj].id();
       bool doProton = rndm.flat() < (np / double(np + nn));
       if (doProton) np -= 1;
       else          nn -= 1;
       int idNuc = (doProton) ? 2212 : 2112;
 
       // Do a projectile-nucleon subcollision. Return empty event if failure.
       pythiaColl.setBeamIDs(idProj, idNuc);
       pythiaColl.setKinematics(eventMain[iProj].p(), Vec4());
       if (!pythiaColl.next(procType)) {
         eventMain.clear();
         return eventMain;
       }
 
       // Insert target nucleon. Mothers are (0,iProj) to mark who it
       // interacted with. Always use proton mass for simplicity.
       int statusNuc = (iColl == 1) ? -181 : -182;
       int iNuc = eventMain.append( idNuc, statusNuc, 0, iProj, 0, 0, 0, 0,
         0., 0., 0., mp, mp);
       eventMain[iNuc].vProdAdd(vNow);
 
       // Update full energy of the event with the proton mass.
       eventMain[0].e( eventMain[0].e() + mp);
       eventMain[0].m( eventMain[0].p().mCalc() );
 
       // Insert secondary produced particles (but skip intermediate partons)
       // into main event record and shift to correct production vertex.
       sizeOld = eventMain.size();
       for (int iSub = 3; iSub < eventColl.size(); ++iSub) {
         if (!eventColl[iSub].isFinal()) continue;
         int iNew = eventMain.append(eventColl[iSub]);
         eventMain[iNew].mothers(iNuc, iProj);
         eventMain[iNew].vProdAdd(vNow);
       }
       sizeNew = eventMain.size();
 
       // Update daughters of colliding hadrons and other history.
       eventMain[iProj].daughters(sizeOld, sizeNew - 1);
       eventMain[iNuc].daughters(sizeOld, sizeNew - 1);
       eventMain[iProj].statusNeg();
       eventMain[iProj].tau(0.);
 
     // End of loop over interactions in a nucleus.
     }
 
     // Optionally do decays of short-lived particles.
     if (rapidDecays) pythiaMain.moreDecays();
 
     // Optionally compress event record.
     if (listFinal) compress();
 
     // Return generated collision.
     return eventMain;
 
   }
 
   //--------------------------------------------------------------------------
 
   // Generate a particle decay, and return the event record.
   // You can allow sequential decays, if they occur rapidly enough.
 
   Event& nextDecay(int idNowIn, Vec4 pNowIn, double mNowIn,
     Vec4 vNow = Vec4() ) {
 
     // Save incoming quantities. (Not needed, but by analogy with
     // collisions.)
     idNow = idNowIn;
     pNow  = pNowIn;
     mNow  = mNowIn;
 
     // References to the event records. Clear them.
     Event& eventMain = pythiaMain.event;
     Event& eventColl = pythiaColl.event;
     eventMain.clear();
     eventColl.clear();
 
     // Insert incoming particle in cleared collision event record.
     eventColl.append(90,   -11, 0, 0, 1, 1, 0, 0, pNow, mNow);
     int iHad = eventColl.append(idNow, 12, 0, 0, 0, 0, 0, 0, pNow, mNow);
     eventColl[iHad].vProd(vNow);
 
     // Decay incoming particle. Return empty event if fail. Copy event record.
     if (!pythiaColl.moreDecays(iHad)) return eventMain;
     eventMain = eventColl;
 
     // Optionally do secondary decays of short-lived particles.
     if (rapidDecays) pythiaMain.moreDecays();
 
     // Optionally compress event record.
     if (listFinal) compress();
 
     // Return generated collision.
     return eventMain;
 
   }
 
   //--------------------------------------------------------------------------
 
   // Compress the event record by removing initial and intermediate
   // particles.  Keep line 0, since the += operator for event records
   // only adds from 1 on.
 
   void compress() {
 
     // Reference to the main event record. Original and new size.
     Event& eventMain = pythiaMain.event;
     int sizeOld = eventMain.size();
     int sizeNew = 1;
 
     // Loop through all particles and move up the final ones to the top.
     // Remove history information. Update event record size.
     for (int i = 1; i < sizeOld; ++i) if (eventMain[i].isFinal()) {
       eventMain[sizeNew] = eventMain[i];
       eventMain[sizeNew].mothers( 0, 0);
       eventMain[sizeNew].daughters( 0, 0);
       ++sizeNew;
     }
 
     // Shrink event record to new size.
     eventMain.popBack( sizeOld - sizeNew);
 
   }
 
   //--------------------------------------------------------------------------
 
   // Summary of aborts, errors and warnings.
 
   void stat() {
     pythiaMain.stat();
     pythiaColl.stat();
     logger.errorStatistics();
   }
 
   //--------------------------------------------------------------------------
 
   // Possibility to access particle data and random numbers from
   // pythiaMain.
 
   ParticleData& particleData() {return pythiaMain.particleData;}
 
   Rndm& rndm() {return pythiaMain.rndm;}
 
 //--------------------------------------------------------------------------
 
 private:
 
   // The Pythia instances used for decays and for collisions. Could
   // be made public, but cleaner to allow only limited access as
   // above.
   Pythia pythiaMain, pythiaColl;
 
   // Logger instance for errors in this class.
   Logger logger;
 
   // Save quantities.
   bool   listFinal, rapidDecays;
   int    idNow;
   double eMax, smallTau0, mp, mNow, eCMNow, sigmaNow;
   Vec4   pNow;
 
   // Fraction of single diffractive events beyond first collision in nucleus.
   static constexpr double probSD = 0.3;
 
   // Hadrons below the kinetic energy threshold cannot interact with medium.
   // (At least not using Pythia.) Used both in lab and in CM frame.
   static constexpr double eKinMin = 0.2;
 
   // Studied nuclei by A number, with offset and slope of <nColl>(sigma):
   // 1H, 2H, 4He, 9Be, 12C, 14N, 16O, 27Al, 40Ar, 56Fe, 63Cu, 84Kr,
   // 107Ag, 129Xe, 197Au, 208Pb.
   static constexpr int nA = 16;
   static constexpr double tabBorder = 20.;
   static const double tabA[nA];
   static const double tabOffset[nA];
   static const double tabSlope[nA];
   static const double tabSlopeLo[nA];
 
 };
 
 // Constant definitions.
 const double PythiaCascade::tabA[] = {1, 2, 4, 9, 12, 14, 16, 27, 40, 56,
   63, 84, 107, 129, 197, 208};
 const double PythiaCascade::tabOffset[] = {0., 0.03, 0.08, 0.15, 0.20, 0.20,
   0.20, 0.26, 0.30, 0.34, 0.40, 0.40, 0.40, 0.50, 0.50, 0.60};
 const double PythiaCascade::tabSlope[] = {0., 0.0016, 0.0033, 0.0075, 0.0092,
   0.0105, 0.012, 0.017, 0.022, 0.027, 0.028, 0.034, 0.040, 0.044, 0.055,
   0.055};
 const double PythiaCascade::tabSlopeLo[] = {0., 0.0031, 0.0073, 0.015, 0.0192,
   0.0205, 0.022, 0.03, 0.037, 0.044, 0.048, 0.054, 0.06, 0.069, 0.08, 0.085};
 
 //==========================================================================
 
 } // end namespace Pythia8
 
 #endif // end Pythia8_PythiaCascade_H

This page was automatically generated by the 2.3.7 LXR engine. The LXR team