Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 09:08:24

0001 // HeavyIons.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This file contains the definition of the HeavyIons class which
0007 // Provides Pythia with infrastructure to combine several nucleon
0008 // collisions into a single heavy ion collision event. This file also
0009 // includes the definition of the Angantyr class which implements the
0010 // default model for heavy ion collisions in Pythia.
0011 
0012 #ifndef Pythia8_HeavyIons_H
0013 #define Pythia8_HeavyIons_H
0014 
0015 #include "Pythia8/HIInfo.h"
0016 #include "Pythia8/PhysicsBase.h"
0017 
0018 namespace Pythia8 {
0019 
0020 // Forward declaration of the Pythia class.
0021 class Pythia;
0022 
0023 //==========================================================================
0024 
0025 // HeavyIons contains several standard Pythia objects to allow for
0026 // the combination of different kinds of nucleon-nucleon events into
0027 // a heavy ion collision event. The actual model for doing this must
0028 // be implemented in a subclass overriding the the virtual'init()'
0029 // and 'next()' functions.
0030 
0031 class HeavyIons : public PhysicsBase {
0032 
0033 public:
0034 
0035   // The constructor needs a reference to the main Pythia object to
0036   // which it will belong. A HeavyIons object cannot belong to more
0037   // than one main Pythia object.
0038   HeavyIons(Pythia& mainPythiaIn)
0039     : mainPythiaPtr(&mainPythiaIn), HIHooksPtr(0),
0040       pythia(vector<Pythia*>(1, &mainPythiaIn)) {
0041   }
0042 
0043   // Destructor.
0044   virtual ~HeavyIons() {}
0045 
0046   // Virtual function to be implemented in a subclass. This will be
0047   // called in the beginning of the Pythia::init function if the mode
0048   // HeavyIon:mode is set non zero. The return value should be true
0049   // if this object is able to handle the requested collision
0050   // type. If false Pythia::init will set HeavyIon:mode to zero but
0051   // will try to cope anyway.
0052   virtual bool init() = 0;
0053 
0054   // Virtual function to be implemented in a subclass. This will be
0055   // called in the beginning of the Pythia::next function if
0056   // HeavyIon:mode is set non zero. After the call, Pythia::next will
0057   // return immediately with the return value of this function.
0058   virtual bool next() = 0;
0059 
0060   // Static function to allow Pythia to duplicate some setting names
0061   // to be used for secondary Pythia objects.
0062   static void addSpecialSettings(Settings& settings);
0063 
0064   // Return true if the beams in the Primary Pythia object contains
0065   // heavy ions.
0066   static bool isHeavyIon(Settings& settings);
0067 
0068   // Possibility to pass in pointer for special heavy ion user hooks.
0069   bool setHIUserHooksPtr(HIUserHooksPtr userHooksPtrIn) {
0070     HIHooksPtr = userHooksPtrIn; return true;
0071   }
0072 
0073   // Set beam kinematics.
0074   virtual bool setKinematics(double /*eCMIn*/) {
0075     loggerPtr->ERROR_MSG("method not implemented for this heavy ion model");
0076     return false; }
0077   virtual bool setKinematics(double /*eAIn*/, double /*eBIn*/) {
0078     loggerPtr->ERROR_MSG("method not implemented for this heavy ion model");
0079     return false; }
0080   virtual bool setKinematics(double, double, double, double, double, double) {
0081     loggerPtr->ERROR_MSG("method not implemented for this heavy ion model");
0082     return false; }
0083   virtual bool setKinematics(Vec4, Vec4) {
0084     loggerPtr->ERROR_MSG("method not implemented for this heavy ion model");
0085     return false; }
0086 
0087   // Set beam particles.
0088   virtual bool setBeamIDs(int /*idAIn*/, int /*idBIn*/ = 0) {
0089     loggerPtr->ERROR_MSG("method not implemented for this heavy ion model");
0090     return false; }
0091 
0092   // The HIInfo object contains information about the last generated heavy
0093   // ion event as well as overall statistics of the generated events.
0094   HIInfo hiInfo;
0095 
0096   // Print out statistics.
0097   virtual void stat();
0098 
0099 protected:
0100 
0101   // Update the cross section in the main Pythia Info object using
0102   // information in the hiInfo object.
0103   void updateInfo();
0104 
0105   // If subclasses has additional Pythia objects for generating
0106   // minimum bias nucleon collisions and the main Pythia object is
0107   // set up to generated a hard signal process, this function can be
0108   // used to clear all selected processes in a clone of the main
0109   // Pythia object.
0110   void clearProcessLevel(Pythia& pyt);
0111 
0112   // Duplicate setting on the form match: to settings on the form HImatch:
0113   static void setupSpecials(Settings& settings, string match);
0114 
0115   // Copy settings on the form HImatch: to the corresponding match:
0116   // in the given Pythia object.
0117   static void setupSpecials(Pythia& p, string match);
0118 
0119   // Save current beam configuration.
0120   int idProj, idTarg;
0121 
0122   // This is the pointer to the main Pythia object to which this
0123   // object is assigned.
0124   Pythia * mainPythiaPtr;
0125 
0126   // Object containing information on inclusive pp cross sections to
0127   // be used in Glauber calculations in subclasses.
0128   SigmaTotal sigTotNN;
0129 
0130   // Optional HIUserHooks object able to modify the behavior of the
0131   // HeavyIon model.
0132   HIUserHooksPtr HIHooksPtr;
0133 
0134   // The internal Pythia objects. Index zero will always correspond
0135   // to the mainPythiaPtr.
0136   vector<Pythia *> pythia;
0137 
0138   // The names associated with the secondary pythia objects.
0139   vector<string> pythiaNames;
0140 
0141   // The Info objects associated to the secondary the secondary
0142   // pythia objects.
0143   vector<Info*> info;
0144 
0145   // Helper class to gain access to the Info object in a pythia
0146   // instance.
0147   struct InfoGrabber : public UserHooks {
0148 
0149     // Only one function: return the info object.
0150     Info * getInfo() {
0151       return infoPtr;
0152     }
0153 
0154   };
0155 
0156 };
0157 
0158 //==========================================================================
0159 
0160 // The default HeavyIon model in Pythia.
0161 
0162 class Angantyr : public HeavyIons {
0163 
0164 public:
0165 
0166   // Enumerate the different internal Pythia objects.
0167   enum PythiaObject {
0168     HADRON = 0,   // For hadronization only.
0169     MBIAS = 1,    // Minimum Bias processed.
0170     SASD = 2,     // Single diffractive as one side of non-diffractive.
0171     SIGPP = 3,    // Optional object for signal processes (pp).
0172     SIGPN = 4,    // Optional object for signal processes (pn).
0173     SIGNP = 5,    // Optional object for signal processes (np).
0174     SIGNN = 6,    // Optional object for signal processes (nn).
0175     ALL = 7       // Indicates all objects.
0176   };
0177 
0178 public:
0179 
0180   // The constructor needs a reference to the main Pythia object to
0181   // which it will belong. A Angantyr object cannot belong to more
0182   // than one main Pythia object.
0183   Angantyr(Pythia& mainPythiaIn);
0184 
0185   virtual ~Angantyr();
0186 
0187   // Initialize Angantyr.
0188   virtual bool init() override;
0189 
0190   // Produce a collision involving heavy ions.
0191   virtual bool next() override;
0192 
0193   // Set UserHooks for specific (or ALL) internal Pythia objects.
0194   bool setUserHooksPtr(PythiaObject sel, UserHooksPtr userHooksPtrIn);
0195 
0196   // Set beam kinematics.
0197   bool setKinematicsCM();
0198   bool setKinematics(double eCMIn) override;
0199   bool setKinematics(double eAIn, double eBIn) override;
0200   bool setKinematics(double, double, double, double, double, double) override;
0201   bool setKinematics(Vec4, Vec4) override;
0202   bool setKinematics();
0203 
0204   // Set beam IDs.
0205   bool setBeamIDs(int idAIn, int idBIn = 0) override;
0206 
0207   // Make sure the correct information is available irrespective of frame type.
0208   void unifyFrames();
0209 
0210   // Print the Angantyr banner.
0211   void banner() const;
0212 
0213   // Subcollisions for the current event.
0214   const SubCollisionSet& subCollisions() const {
0215     return subColls; }
0216 
0217   // Get the underlying subcollision model.
0218   const SubCollisionModel& subCollisionModel() const {
0219     return *collPtr.get(); }
0220 
0221   SubCollisionModel* subCollPtr() {
0222     return collPtr.get();
0223   }
0224 
0225   // Get the underlying impact parameter generator.
0226   const ImpactParameterGenerator impactParameterGenerator() const {
0227     return *bGenPtr.get(); }
0228 
0229   // Projectile nucleus configuration for the current event.
0230   const Nucleus& projectile() const {
0231     return proj; }
0232 
0233   // Target nucleus configuration for the current event.
0234   const Nucleus& target() const {
0235     return targ; }
0236 
0237   // The underlying projectile nucleus model.
0238   const NucleusModel& projectileModel() const {
0239     return *projPtr.get(); }
0240 
0241   // The underlying target nucleus model.
0242   const NucleusModel& targetModel() const {
0243     return *targPtr.get(); }
0244 
0245   // Hadronic cross sections used by the subcollision model.
0246   const SigmaTotal sigmaNN() const {
0247     return sigTotNN; }
0248 
0249 protected:
0250 
0251   virtual void onInitInfoPtr() override {
0252     registerSubObject(sigTotNN); }
0253 
0254   // Figure out what beams the user want.
0255   void setBeamKinematics(int idA, int idB);
0256 
0257   // Initiaize a specific Pythia object and optionally run a number
0258   // of events to get a handle of the cross section.
0259   bool init(PythiaObject sel, string name, int n = 0);
0260 
0261   // Setup an EventInfo object from a Pythia instance.
0262   EventInfo mkEventInfo(Pythia &, Info &, const SubCollision * coll = 0);
0263 
0264   // Generate events from the internal Pythia oblects;
0265   EventInfo getSignal(const SubCollision& coll);
0266   EventInfo getND() { return getMBIAS(0, 101); }
0267   EventInfo getND(const SubCollision& coll) { return getMBIAS(&coll, 101); }
0268   EventInfo getEl(const SubCollision& coll) { return getMBIAS(&coll, 102); }
0269   EventInfo getSDP(const SubCollision& coll) { return getMBIAS(&coll, 103); }
0270   EventInfo getSDT(const SubCollision& coll) { return getMBIAS(&coll, 104); }
0271   EventInfo getDD(const SubCollision& coll) { return getMBIAS(&coll, 105); }
0272   EventInfo getCD(const SubCollision& coll) { return getMBIAS(&coll, 106); }
0273   EventInfo getSDabsP(const SubCollision& coll)
0274     { return getSASD(&coll, 103); }
0275   EventInfo getSDabsT(const SubCollision& coll)
0276     { return getSASD(&coll, 104); }
0277   EventInfo getMBIAS(const SubCollision * coll, int procid);
0278   EventInfo getSASD(const SubCollision * coll, int procid);
0279 
0280   bool genAbs(SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0281   void addSASD(const SubCollisionSet& subCollsIn);
0282   bool addDD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0283   bool addSD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0284   void addSDsecond(const SubCollisionSet& subCollsIn);
0285   bool addCD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0286   void addCDsecond(const SubCollisionSet& subCollsIn);
0287   bool addEL(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0288   void addELsecond(const SubCollisionSet& subCollsIn);
0289 
0290   void resetEvent();
0291   bool buildEvent(list<EventInfo>& subEventsIn);
0292 
0293   bool setupFullCollision(EventInfo& ei, const SubCollision& coll,
0294     Nucleon::Status projStatus, Nucleon::Status targStatus);
0295   bool isRemnant(const EventInfo& ei, int i, int past = 1 ) const {
0296     int statNow = ei.event[i].status()*past;
0297     if ( statNow == 63 ) return true;
0298     if ( statNow > 70 && statNow < 80 )
0299       return isRemnant(ei, ei.event[i].mother1(), -1);
0300     return false;
0301   }
0302   bool fixIsoSpin(EventInfo& ei);
0303   EventInfo& shiftEvent(EventInfo& ei);
0304   static int getBeam(Event& ev, int i);
0305 
0306   // Generate a single diffractive
0307   bool nextSASD(int proc);
0308 
0309   // Add a diffractive event to an exsisting one. Optionally connect
0310   // the colours of the added event to the original.
0311   bool addNucleonExcitation(EventInfo& orig, EventInfo& add,
0312                             bool colConnect = false);
0313 
0314   // Find the recoilers in the current event to conserve energy and
0315   // momentum in addNucleonExcitation.
0316   vector<int> findRecoilers(const Event& e, bool tside, int beam, int end,
0317                             const Vec4& pdiff, const Vec4& pbeam);
0318 
0319   // Add a sub-event to the final event record.
0320   void addSubEvent(Event& evnt, Event& sub);
0321   static void addJunctions(Event& evnt, Event& sub, int coloff);
0322 
0323   // Add a nucleus remnant to the given event. Possibly introducing
0324   // a new particle type.
0325   bool addNucleusRemnants();
0326 
0327 public:
0328 
0329   // Helper function to construct two transformations that would give
0330   // the vectors p1 and p2 the total four momentum of p1p + p2p.
0331   static bool
0332   getTransforms(Vec4 p1, Vec4 p2, const Vec4& p1p,
0333                 pair<RotBstMatrix,RotBstMatrix>& R12);
0334   static double mT2(const Vec4& p) { return p.pPos()*p.pNeg(); }
0335   static double mT(const Vec4& p) { return sqrt(max(mT2(p), 0.0)); }
0336 
0337 private:
0338 
0339   // Private UserHooks class to select a specific process.
0340   struct ProcessSelectorHook: public UserHooks {
0341 
0342     ProcessSelectorHook(): proc(0), b(-1.0) {}
0343 
0344     // Yes we can veto event after process-level selection.
0345     virtual bool canVetoProcessLevel() {
0346       return true;
0347     }
0348 
0349     // Veto any unwanted process.
0350     virtual bool doVetoProcessLevel(Event&) {
0351       return proc > 0 && infoPtr->code() != proc;
0352     }
0353 
0354     // Can set the overall impact parameter for the MPI treatment.
0355     virtual bool canSetImpactParameter() const {
0356       return b >= 0.0;
0357     }
0358 
0359     // Set the overall impact parameter for the MPI treatment.
0360     virtual double doSetImpactParameter() {
0361       return b;
0362     }
0363 
0364     // The wanted process;
0365     int proc;
0366 
0367     // The selected b-value.
0368     double b;
0369 
0370   };
0371 
0372   // Holder class to temporarily select a specific process
0373   struct HoldProcess {
0374 
0375     // Set the given process for the given hook object.
0376     HoldProcess(shared_ptr<ProcessSelectorHook> hook, int proc,
0377       double b = -1.0) : saveHook(hook), saveProc(hook->proc), saveB(hook->b) {
0378       hook->proc = proc;
0379       hook->b = b;
0380     }
0381 
0382     // Reset the process of the hook object given in the constructor.
0383     ~HoldProcess() {
0384       if ( saveHook ) {
0385         saveHook->proc = saveProc;
0386         saveHook->b = saveB;
0387       }
0388     }
0389 
0390     // The hook object.
0391     shared_ptr<ProcessSelectorHook> saveHook;
0392 
0393     // The previous process of the hook object.
0394     int saveProc;
0395 
0396     // The previous b-value of the hook object.
0397     double saveB;
0398 
0399   };
0400 
0401   // The process selector for standard minimum bias processes.
0402   shared_ptr<ProcessSelectorHook> selectMB;
0403 
0404   // The process selector for the SASD object.
0405   shared_ptr<ProcessSelectorHook> selectSASD;
0406 
0407 private:
0408 
0409   static const int MAXTRY = 999;
0410   static const int MAXEVSAVE = 999;
0411 
0412   // Flag set if there is a specific signal process specified beyond
0413   // minimum bias.
0414   bool hasSignal;
0415 
0416   // Whether to do hadronization and hadron level effects.
0417   bool doHadronLevel;
0418 
0419   // Flag to determine whether to do single diffractive test.
0420   bool doSDTest;
0421 
0422   // Flag to determine whether to do only Glauber modelling.
0423   bool glauberOnly;
0424 
0425   // All subcollisions in current collision.
0426   SubCollisionSet subColls;
0427 
0428   // The underlying SubCollisionModel for generating nucleon-nucleon
0429   // subcollisions.
0430   shared_ptr<SubCollisionModel> collPtr;
0431 
0432   // The impact parameter generator.
0433   shared_ptr<ImpactParameterGenerator> bGenPtr;
0434 
0435   // The projectile and target nuclei in the current collision.
0436   Nucleus proj;
0437   Nucleus targ;
0438 
0439   // The underlying nucleus model for generating nuclons inside the
0440   // projectile and target nucleus.
0441   shared_ptr<NucleusModel> projPtr;
0442   shared_ptr<NucleusModel> targPtr;
0443 
0444   // Flag to indicate whether variable energy is enabled.
0445   bool doVarECM;
0446 
0447   // Different choices in choosing recoilers when adding
0448   // diffractively excited nucleon.
0449   int recoilerMode;
0450 
0451   // Different choices for handling impact parameters.
0452   int bMode;
0453 
0454   // Critical internal error, abort the event.
0455   bool doAbort;
0456 
0457 };
0458 
0459 //==========================================================================
0460 
0461 } // end namespace Pythia8
0462 
0463 #endif // Pythia8_HeavyIons_H