Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:06:24

0001 // HeavyIons.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 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 setKinematics(double eCMIn) override;
0198   bool setKinematics(double eAIn, double eBIn) override;
0199   bool setKinematics(double, double, double, double, double, double) override;
0200   bool setKinematics(Vec4, Vec4) override;
0201   bool setKinematics();
0202 
0203   // Set beam IDs.
0204   bool setBeamIDs(int idAIn, int idBIn = 0) override;
0205 
0206   // Make sure the correct information is available irrespective of frame type.
0207   void unifyFrames();
0208 
0209   // Print the Angantyr banner.
0210   void banner() const;
0211 
0212   // Subcollisions for the current event.
0213   const SubCollisionSet& subCollisions() const {
0214     return subColls; }
0215 
0216   // Get the underlying subcollision model.
0217   const SubCollisionModel& subCollisionModel() const {
0218     return *collPtr.get(); }
0219 
0220   // Get the underlying impact parameter generator.
0221   const ImpactParameterGenerator impactParameterGenerator() const {
0222     return *bGenPtr.get(); }
0223 
0224   // Projectile nucleus configuration for the current event.
0225   const Nucleus& projectile() const {
0226     return proj; }
0227 
0228   // Target nucleus configuration for the current event.
0229   const Nucleus& target() const {
0230     return targ; }
0231 
0232   // The underlying projectile nucleus model.
0233   const NucleusModel& projectileModel() const {
0234     return *projPtr.get(); }
0235 
0236   // The underlying target nucleus model.
0237   const NucleusModel& targetModel() const {
0238     return *targPtr.get(); }
0239 
0240   // Hadronic cross sections used by the subcollision model.
0241   const SigmaTotal sigmaNN() const {
0242     return sigTotNN; }
0243 
0244 protected:
0245 
0246   virtual void onInitInfoPtr() override {
0247     registerSubObject(sigTotNN); }
0248 
0249   // Figure out what beams the user want.
0250   void setBeamKinematics(int idA, int idB);
0251 
0252   // Initiaize a specific Pythia object and optionally run a number
0253   // of events to get a handle of the cross section.
0254   bool init(PythiaObject sel, string name, int n = 0);
0255 
0256   // Setup an EventInfo object from a Pythia instance.
0257   EventInfo mkEventInfo(Pythia &, Info &, const SubCollision * coll = 0);
0258 
0259   // Generate events from the internal Pythia oblects;
0260   EventInfo getSignal(const SubCollision& coll);
0261   EventInfo getND() { return getMBIAS(0, 101); }
0262   EventInfo getND(const SubCollision& coll) { return getMBIAS(&coll, 101); }
0263   EventInfo getEl(const SubCollision& coll) { return getMBIAS(&coll, 102); }
0264   EventInfo getSDP(const SubCollision& coll) { return getMBIAS(&coll, 103); }
0265   EventInfo getSDT(const SubCollision& coll) { return getMBIAS(&coll, 104); }
0266   EventInfo getDD(const SubCollision& coll) { return getMBIAS(&coll, 105); }
0267   EventInfo getCD(const SubCollision& coll) { return getMBIAS(&coll, 106); }
0268   EventInfo getSDabsP(const SubCollision& coll)
0269     { return getSASD(&coll, 103); }
0270   EventInfo getSDabsT(const SubCollision& coll)
0271     { return getSASD(&coll, 104); }
0272   EventInfo getMBIAS(const SubCollision * coll, int procid);
0273   EventInfo getSASD(const SubCollision * coll, int procid);
0274 
0275   bool genAbs(SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0276   void addSASD(const SubCollisionSet& subCollsIn);
0277   bool addDD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0278   bool addSD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0279   void addSDsecond(const SubCollisionSet& subCollsIn);
0280   bool addCD(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0281   void addCDsecond(const SubCollisionSet& subCollsIn);
0282   bool addEL(const SubCollisionSet& subCollsIn, list<EventInfo>& subEventsIn);
0283   void addELsecond(const SubCollisionSet& subCollsIn);
0284 
0285   bool buildEvent(list<EventInfo>& subEventsIn);
0286 
0287   bool setupFullCollision(EventInfo& ei, const SubCollision& coll,
0288     Nucleon::Status projStatus, Nucleon::Status targStatus);
0289   bool isRemnant(const EventInfo& ei, int i, int past = 1 ) const {
0290     int statNow = ei.event[i].status()*past;
0291     if ( statNow == 63 ) return true;
0292     if ( statNow > 70 && statNow < 80 )
0293       return isRemnant(ei, ei.event[i].mother1(), -1);
0294     return false;
0295   }
0296   bool fixIsoSpin(EventInfo& ei);
0297   EventInfo& shiftEvent(EventInfo& ei);
0298   static int getBeam(Event& ev, int i);
0299 
0300   // Generate a single diffractive
0301   bool nextSASD(int proc);
0302 
0303   // Add a diffractive event to an exsisting one. Optionally connect
0304   // the colours of the added event to the original.
0305   bool addNucleonExcitation(EventInfo& orig, EventInfo& add,
0306                             bool colConnect = false);
0307 
0308   // Find the recoilers in the current event to conserve energy and
0309   // momentum in addNucleonExcitation.
0310   vector<int> findRecoilers(const Event& e, bool tside, int beam, int end,
0311                             const Vec4& pdiff, const Vec4& pbeam);
0312 
0313   // Add a sub-event to the final event record.
0314   void addSubEvent(Event& evnt, Event& sub);
0315   static void addJunctions(Event& evnt, Event& sub, int coloff);
0316 
0317   // Add a nucleus remnant to the given event. Possibly introducing
0318   // a new particle type.
0319   bool addNucleusRemnants();
0320 
0321 public:
0322 
0323   // Helper function to construct two transformations that would give
0324   // the vectors p1 and p2 the total four momentum of p1p + p2p.
0325   static bool
0326   getTransforms(Vec4 p1, Vec4 p2, const Vec4& p1p,
0327                 pair<RotBstMatrix,RotBstMatrix>& R12);
0328   static double mT2(const Vec4& p) { return p.pPos()*p.pNeg(); }
0329   static double mT(const Vec4& p) { return sqrt(max(mT2(p), 0.0)); }
0330 
0331 private:
0332 
0333   // Private UserHooks class to select a specific process.
0334   struct ProcessSelectorHook: public UserHooks {
0335 
0336     ProcessSelectorHook(): proc(0), b(-1.0) {}
0337 
0338     // Yes we can veto event after process-level selection.
0339     virtual bool canVetoProcessLevel() {
0340       return true;
0341     }
0342 
0343     // Veto any unwanted process.
0344     virtual bool doVetoProcessLevel(Event&) {
0345       return proc > 0 && infoPtr->code() != proc;
0346     }
0347 
0348     // Can set the overall impact parameter for the MPI treatment.
0349     virtual bool canSetImpactParameter() const {
0350       return b >= 0.0;
0351     }
0352 
0353     // Set the overall impact parameter for the MPI treatment.
0354     virtual double doSetImpactParameter() {
0355       return b;
0356     }
0357 
0358     // The wanted process;
0359     int proc;
0360 
0361     // The selected b-value.
0362     double b;
0363 
0364   };
0365 
0366   // Holder class to temporarily select a specific process
0367   struct HoldProcess {
0368 
0369     // Set the given process for the given hook object.
0370     HoldProcess(shared_ptr<ProcessSelectorHook> hook, int proc,
0371       double b = -1.0) : saveHook(hook), saveProc(hook->proc), saveB(hook->b) {
0372       hook->proc = proc;
0373       hook->b = b;
0374     }
0375 
0376     // Reset the process of the hook object given in the constructor.
0377     ~HoldProcess() {
0378       if ( saveHook ) {
0379         saveHook->proc = saveProc;
0380         saveHook->b = saveB;
0381       }
0382     }
0383 
0384     // The hook object.
0385     shared_ptr<ProcessSelectorHook> saveHook;
0386 
0387     // The previous process of the hook object.
0388     int saveProc;
0389 
0390     // The previous b-value of the hook object.
0391     double saveB;
0392 
0393   };
0394 
0395   // The process selector for standard minimum bias processes.
0396   shared_ptr<ProcessSelectorHook> selectMB;
0397 
0398   // The process selector for the SASD object.
0399   shared_ptr<ProcessSelectorHook> selectSASD;
0400 
0401 private:
0402 
0403   static const int MAXTRY = 999;
0404   static const int MAXEVSAVE = 999;
0405 
0406   // Flag set if there is a specific signal process specified beyond
0407   // minimum bias.
0408   bool hasSignal;
0409 
0410   // Whether to do hadronization and hadron level effects.
0411   bool doHadronLevel;
0412 
0413   // Flag to determine whether to do single diffractive test.
0414   bool doSDTest;
0415 
0416   // Flag to determine whether to do only Glauber modelling.
0417   bool glauberOnly;
0418 
0419   // All subcollisions in current collision.
0420   SubCollisionSet subColls;
0421 
0422   // The underlying SubCollisionModel for generating nucleon-nucleon
0423   // subcollisions.
0424   shared_ptr<SubCollisionModel> collPtr;
0425 
0426   // The impact parameter generator.
0427   shared_ptr<ImpactParameterGenerator> bGenPtr;
0428 
0429   // The projectile and target nuclei in the current collision.
0430   Nucleus proj;
0431   Nucleus targ;
0432 
0433   // The underlying nucleus model for generating nuclons inside the
0434   // projectile and target nucleus.
0435   shared_ptr<NucleusModel> projPtr;
0436   shared_ptr<NucleusModel> targPtr;
0437 
0438   // Flag to indicate whether variable energy is enabled.
0439   bool doVarECM;
0440 
0441   // Different choices in choosing recoilers when adding
0442   // diffractively excited nucleon.
0443   int recoilerMode;
0444 
0445   // Different choices for handling impact parameters.
0446   int bMode;
0447 
0448   // Critical internal error, abort the event.
0449   bool doAbort;
0450 
0451 };
0452 
0453 //==========================================================================
0454 
0455 } // end namespace Pythia8
0456 
0457 #endif // Pythia8_HeavyIons_H