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

File indexing completed on 2025-12-15 10:23:35

 // VinciaQED.h is a part of the PYTHIA event generator.
 // Copyright (C) 2025 Peter Skands, 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.
 
 // This file contains the QED antenna-shower class and auxiliary
 // classes. Main author is Rob Verheyen.
 
 #ifndef Pythia8_VinciaQED_H
 #define Pythia8_VinciaQED_H
 
 // PYTHIA 8 headers.
 #include "Pythia8/BeamParticle.h"
 #include "Pythia8/Event.h"
 #include "Pythia8/StandardModel.h"
 #include "Pythia8/PartonSystems.h"
 
 // VINCIA headers.
 #include "Pythia8/VinciaCommon.h"
 #include "Pythia8/VinciaWeights.h"
 
 namespace Pythia8 {
 
 //==========================================================================
 
 // Class for QED emissions.
 
 class QEDemitElemental {
 
 public:
 
   // Friends for internal private members.
   friend class QEDemitSystem;
 
   // Constuctor.
   QEDemitElemental() : rndmPtr(nullptr), partonSystemsPtr(nullptr), q2Sav(0.),
     zetaSav(0.), phiSav(0.), sxjSav(0.), syjSav(0.), alpha(0.), c(0.),
     hasTrial(false), x(0), y(0), idx(0), idy(0), mx2(0.), my2(0.),
     ex(0.), ey(0.), m2Ant(0.), sAnt(0.), QQ(0.), isII(false), isIF(false),
     isFF(false), isRF(false), isIA(true), isDip(false), shh(0.),
     isInitPtr(false), isInit(false), verbose(1) {;}
 
   // Initialize the pointers.
   void initPtr(Rndm* rndmPtrIn, PartonSystems* partonSystemsPtrIn);
   // Initialize.
   void init(Event &event, int xIn, int yIn, double shhIn,
     double verboseIn);
   // Initialize.
   void init(Event &event, int xIn, vector<int> iRecoilIn, double shhIn,
     double verboseIn);
   // Generate a trial point.
   double generateTrial(Event &event, double q2Start, double q2Low,
     double alphaIn, double cIn);
 
 private:
 
   // Random pointer.
   Rndm* rndmPtr{};
 
   // Parton system pointer.
   PartonSystems* partonSystemsPtr{};
 
   // Trial variables.
   double q2Sav, zetaSav, phiSav, sxjSav, syjSav, alpha, c;
   bool hasTrial;
 
   // Particle indices.
   int x, y;
   // Recoiler indices.
   vector<int> iRecoil;
   // IDs.
   int idx, idy;
   // Number of spin states.
   int spinTypex, spinTypey;
   // Particle masses.
   double mx2, my2;
   // Particle energies.
   double ex, ey;
   // Antenna invariant mass.
   double m2Ant;
   // Antenna dot product.
   double sAnt;
   // The negative of the product of charges.
   double QQ;
 
   // Type switches.
   bool isII, isIF, isFF, isRF, isIA, isDip;
 
   // Hadronic invariant mass.
   double shh;
 
   // Initialization.
   bool isInitPtr, isInit;
   int verbose;
 
 };
 
 //==========================================================================
 
 // Base class for QED systems.
 
 class QEDsystem {
 
  public:
 
   // Constructor.
   QEDsystem() : infoPtr(nullptr), partonSystemsPtr(nullptr),
     particleDataPtr(nullptr), rndmPtr(nullptr), settingsPtr(nullptr),
     loggerPtr(nullptr), vinComPtr(nullptr), isInitPtr(false), iSys(-1),
     verbose(0), jNew(0), shat(0.) {;}
 
   // Destructor.
   virtual ~QEDsystem() = default;
 
   // Initialize pointers.
   void initPtr(Info* infoPtrIn, ParticleData* particleDataPtrIn,
     PartonSystems* partonSystemsPtrIn, Rndm* rndmPtrIn,
     Settings* settingsPtrIn, VinciaCommon* vinComPtrIn);
 
   // Initialise settings for current run.
   virtual void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
     int verboseIn) = 0;
   virtual void setVerbose(int verboseIn) { verbose = verboseIn; }
   // Prepare a parton system for evolution.
   virtual void prepare(const int iSysIn, Event &event, const double q2CutIn,
     const int scaleRegionIn, const vector<double> evolutionWindowsIn,
     AlphaEM alIn) = 0;
   // Build parton system.
   virtual void buildSystem(Event &event) = 0;
   // Generate a trial scale.
   virtual double q2Next(Event &event, double q2Start) = 0 ;
   // Generate kinematics and check for veto of branching.
   virtual bool acceptTrial(Event &event) = 0;
   // Update the envent.
   virtual void updateEvent(Event &event) = 0;
   // Update the parton systems.
   virtual void updatePartonSystems();
   // Print information about the system.
   virtual void print() = 0;
 
   // Methods to tell which type of brancher this is.
   virtual bool isSplitting() {return false;};
   virtual bool isInitial() {return false;};
 
  protected:
 
   // Pointers.
   Info* infoPtr{};
   PartonSystems* partonSystemsPtr{};
   ParticleData* particleDataPtr{};
   Rndm* rndmPtr{};
   Settings* settingsPtr{};
   Logger* loggerPtr{};
   VinciaCommon* vinComPtr{};
   bool isInitPtr;
 
   // Event system.
   int iSys;
   vector<Vec4> pNew;
 
   // Verbose setting.
   int verbose;
 
   // Information for partonSystems.
   int jNew;
   map<int,int> iReplace;
   double shat;
 
 };
 
 //==========================================================================
 
 // Class for a QED emission system.
 
 class QEDemitSystem : public QEDsystem {
 
 public:
 
   QEDemitSystem() : shh(-1.), cMat(0.), trialIsVec(false), beamAPtr(nullptr),
     beamBPtr(nullptr), qedMode(-1), qedModeMPI(-1),
     scaleRegion(0), emitBelowHad(false), q2Cut(-1.), isInit(false),
     TINYPDF(-1.), kMapTypeFinal(0) {;}
 
   // Initialise settings for current run.
   void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
     int verboseIn) override;
   // Prepare a parton system for photon emission evolution
   void prepare(const int iSysIn, Event &event, const double q2CutIn,
     const int scaleRegionIn, const vector<double> evolutionWindowsIn,
     AlphaEM alIn) override;
   // Set up antenna pairing for incoherent mode.
   void buildSystem(Event &event) override;
   // Generate a trial scale.
   double q2Next(Event &event, double q2Start) override;
   // Generate kinematics and check veto.
   bool acceptTrial(Event &event) override;
   // Update the envent
   void updateEvent(Event &event) override;
   // Branching type.
   bool isInitial() override {return eleTrial->isII || eleTrial->isIF;}
   // Print the QED emit internal system.
   void print() override;
 
   // Trial antenna function.
   double aTrial(QEDemitElemental* ele, double sxj, double syj, double sxy);
   // Physical antenna function.
   double aPhys (QEDemitElemental* ele, double sxj, double syj, double sxy);
   // Ratio between PDFs.
   double pdfRatio(bool isA, double eOld, double eNew, int id, double Qt2);
 
   private:
 
   // Event system.
   double shh;
 
   // Internal storage.
   vector<vector<QEDemitElemental> > eleMat;
   vector<int> iCoh;
   double cMat;
   vector<QEDemitElemental> eleVec;
 
   // AlphaEM.
   AlphaEM al;
 
   // Evolution window.
   vector<double> evolutionWindows;
 
   // Trial pointer.
   QEDemitElemental* eleTrial{};
   bool trialIsVec;
 
   // Beam pointers.
   BeamParticle* beamAPtr{};
   BeamParticle* beamBPtr{};
 
   // Settings.
   int qedMode, qedModeMPI;
   vector<bool> useSpinsQEDNow, useSpinsQED, useSpinsQEDHadDec;
   int scaleRegion, emitBelowHad, isHadronDecay;
   double q2Cut;
 
   // Initialization.
   bool isInit;
 
   // PDF check.
   double TINYPDF;
 
   // Recoil strategy.
   int kMapTypeFinal;
 
   // Global recoil momenta.
   Vec4 pRecSum;
   vector<Vec4> pRec;
   vector<int>  iRec;
 
 };
 
 //==========================================================================
 
 // Class for trial QED splittings.
 
 class QEDsplitElemental {
 
 public:
 
   // Friends for internal private members.
   friend class QEDsplitSystem;
 
   // Constructor.
   QEDsplitElemental(Event &event, int iPhotIn, int iSpecIn):
     iPhot(iPhotIn), iSpec(iSpecIn), ariWeight(0) {
     m2Ant = max(VinciaConstants::PICO,
       m2(event[iPhotIn], event[iSpecIn]));
     sAnt = max(VinciaConstants::PICO,
       2.*event[iPhotIn].p()*event[iSpecIn].p());
     m2Spec = max(0., event[iSpecIn].m2());
   }
 
   // Kallen function.
   double getKallen() {return m2Ant/(m2Ant - m2Spec);}
 
 private:
 
   // Internal members.
   int iPhot, iSpec;
   double m2Spec, m2Ant, sAnt;
   double ariWeight;
 };
 
 //==========================================================================
 
 // Class for a QED splitting system.
 
 class QEDsplitSystem : public QEDsystem {
 
 public:
 
   // Constructor.
   QEDsplitSystem() :
     totIdWeight(-1.), hasTrial(false),
     q2Trial(-1.), zTrial(-1.), phiTrial(-1.), idTrial(0), nQuark(-1),
     nLepton(-1), q2Max(-1.), q2Cut(-1.), scaleRegion(0),
     beamAPtr(nullptr), beamBPtr(nullptr), isInit(false), kMapTypeFinal(0) {;}
 
   // Initialize.
   void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
     int verboseIn) override;
   // Prepare list of final-state photons - with recoilers - for splittings.
   void prepare(const int iSysIn, Event &event, const double q2CutIn,
     const int scaleRegionIn, const vector<double> evolutionWindowsIn,
     AlphaEM alIn) override;
   // Build the splitting system.
   void buildSystem(Event &event) override;
   // Generate a scale for the system.
   double q2Next(Event &event, double q2Start) override;
   // Generate kinematics and check veto.
   bool acceptTrial(Event &event) override;
   // Update Event.
   void updateEvent(Event &event) override;
   // Branching type: isSplitting() = true.
   bool isSplitting() override { return true;}
   // Print the system.
   void print() override;
 
 private:
 
   // AlphaEM.
   AlphaEM al;
 
   // Evolution window.
   vector<double> evolutionWindows;
 
   // Weights for splitting IDs.
   vector<int> ids;
   vector<double> idWeights;
   double totIdWeight;
 
   // Antennae.
   vector<QEDsplitElemental> eleVec;
 
   // Trial variables.
   bool hasTrial;
   double q2Trial, zTrial, phiTrial, idTrial;
   QEDsplitElemental* eleTrial{};
 
   // Settings.
   int nQuark, nLepton;
   double q2Max, q2Cut;
   int scaleRegion;
 
   // Pointers.
   BeamParticle*  beamAPtr;
   BeamParticle*  beamBPtr;
 
   // Initialization.
   bool isInit;
 
   // Recoil strategy.
   int kMapTypeFinal;
 
 };
 
 //==========================================================================
 
 // Class for a QED conversion system.
 
 class QEDconvSystem : public QEDsystem {
 
 public:
 
   // Constructor.
   QEDconvSystem() : totIdWeight(-1.), maxIdWeight(-1.), shh(-1.), s(-1.),
     iA(-1), iB(-1), isAPhot(false), isBPhot(false), hasTrial(false),
     iPhotTrial(-1), iSpecTrial(-1), q2Trial(-1.), zTrial(-1.), phiTrial(-1.),
     idTrial(-1), nQuark(-1), q2Cut(-1.),
     scaleRegion(0), beamAPtr(nullptr), beamBPtr(nullptr),isInit(false),
     TINYPDF(-1.) {;}
 
   // Initialize.
   void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
     int verboseIn) override;
   // Prepare for backwards-evolution of photons.
   void prepare(const int iSysIn, Event &event, const double q2CutIn,
     const int scaleRegionIn, const vector<double> evolutionWindowsIn,
     AlphaEM alIn) override;
   // Build the system.
   void buildSystem(Event &event) override;
   // Generate a trial scale.
   double q2Next(Event &event, double q2Start) override;
   // Generate kinematics and check veto.
   bool acceptTrial(Event &event) override;
   // Update.
   void updateEvent(Event &event) override;
   // Branching type: isInitial() = true.
   bool isInitial() override {return true;};
   // Print.
   void print() override;
 
 private:
 
   // Trial pdf ratios for conversion.
   map<int,double> Rhat;
 
   // AlphaEM.
   AlphaEM al;
 
   // Evolution window.
   vector<double> evolutionWindows;
 
   // Weights for conversion IDs.
   vector<int> ids;
   vector<double> idWeights;
   double totIdWeight, maxIdWeight;
   double shh;
 
   // Antenna parameters.
   double s;
   int iA, iB;
   bool isAPhot, isBPhot;
 
   // Trial variables.
   bool hasTrial;
   int iPhotTrial, iSpecTrial;
   double q2Trial, zTrial, phiTrial, idTrial;
 
   // Settings.
   int nQuark;
   double q2Cut;
 
   // Scale Region determines what region and scale we are showering at.
   // Scale Region = 0 is above the hadronisation scale.
   // Scale Region = 1 is below the hadronisation scale.
   // Scale Region = 2 is below the hadronisation scale and contains remanants.
   int scaleRegion;
 
   // Pointers.
   BeamParticle*  beamAPtr;
   BeamParticle*  beamBPtr;
 
   // Initialization.
   bool isInit;
 
   // PDF check.
   double TINYPDF;
 
   // Global recoil momenta
   vector<Vec4> pRec;
   vector<int>  iRec;
 
 };
 
 //==========================================================================
 
 // Base class for Vincia's QED and EW shower modules.
 
 class VinciaModule {
 
 public:
 
   // Default constructor.
   VinciaModule(): verbose(-1), isInitPtr(false), isInitSav(false) {;};
 
   // Default destructor.
   virtual ~VinciaModule() = default;
 
   // Initialise pointers (called at construction time).
   virtual void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) = 0;
 
   //Some early initialisation (needed for EW shower).
   virtual void load() {;}
 
   // Initialise settings for current run (called as part of Pythia::init()).
   virtual void init(BeamParticle* beamAPtrIn = nullptr,
     BeamParticle* beamBPtrIn = nullptr) = 0;
   bool isInit() {return isInitSav;}
 
   // Select helicities for a system of particles.
   virtual bool polarise(vector<Particle>&) {return false;}
 
   // Prepare to shower a system.
   virtual bool prepare(int iSysIn, Event &event, int scaleRegionIn) = 0;
 
   // Update shower system each time something has changed in event.
   virtual void update(Event &event, int iSys) = 0;
 
   // Set verbosity level.
   virtual void setVerbose(int verboseIn) {verbose = verboseIn;}
 
   // Clear everything (optionally for specific system).
   virtual void clear(int iSys = -1) = 0;
 
   // Generate a trial scale.
   virtual double q2Next(Event &event, double q2Start, double q2End) = 0;
 
   // Which system does the winner belong to?
   virtual int sysWin() = 0;
 
   // Get information about the latest branching.
   virtual bool lastIsSplitting() = 0;
   virtual bool lastIsInitial() = 0;
   virtual bool lastIsResonanceDecay() {return false;}
 
   // Generate kinematics and check veto.
   virtual bool acceptTrial(Event &event) = 0;
 
   // Update event after branching accepted.
   virtual void updateEvent(Event &event) = 0;
 
   // Update partonSystems after branching accepted.
   virtual void updatePartonSystems(Event &event) = 0;
 
   // End scales.
   virtual double q2minColoured() = 0;
   virtual double q2min() = 0;
 
   // Get number of branchers / systems.
   virtual unsigned int nBranchers() = 0;
   virtual unsigned int nResDec() = 0;
 
   // Members.
   BeamParticle* beamAPtr{};
   BeamParticle* beamBPtr{};
   Info* infoPtr{};
   ParticleData* particleDataPtr{};
   Logger* loggerPtr{};
   PartonSystems* partonSystemsPtr{};
   Rndm* rndmPtr{};
   Settings* settingsPtr{};
   VinciaCommon* vinComPtr{};
 
  protected:
 
   int verbose;
   bool isInitPtr, isInitSav;
 
 };
 
 //==========================================================================
 
 // Class for performing QED showers.
 
 class VinciaQED : public VinciaModule {
 
 public:
 
   // Friends for internal private members.
   friend class VinciaFSR;
 
   // Constructor.
   VinciaQED() {;}
 
   // Initialise pointers (called at construction time).
   void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override;
   // Initialise settings for current run (called as part of Pythia::init()).
   void init(BeamParticle* beamAPtrIn = 0, BeamParticle* beamBPtrIn = 0)
     override;
   // Prepare to shower a system.
   bool prepare(int iSysIn, Event& event, int scaleRegionIn) override;
   // Update QED shower system(s) each time something has changed in event.
   void update(Event& event, int iSys) override;
   // Set or change verbosity level, and propagate to QED systems.
   virtual void setVerbose(int verboseIn) override {
     verbose = verboseIn;
     emptyQEDemitSystem.setVerbose(verboseIn);
     emptyQEDsplitSystem.setVerbose(verboseIn);
     emptyQEDconvSystem.setVerbose(verboseIn);
   }
   // Clear everything, or specific system.
   void clear(int iSys = -1) override {
     if (iSys < 0) {emitSystems.clear(); splitSystems.clear();
       convSystems.clear();}
     else {emitSystems.erase(iSys); splitSystems.erase(iSys);
       convSystems.erase(iSys);}
     qedTrialSysPtr = nullptr;}
 
   // Generate a trial scale.
   double q2Next(Event& event, double q2Start, double) override;
   // Return the system window.
   int  sysWin() override {return iSysTrial;}
   // Information about last branching.
   bool lastIsSplitting() override {
     if (qedTrialSysPtr != nullptr) return  qedTrialSysPtr->isSplitting();
     else return false;}
   bool lastIsInitial() override {
     if (qedTrialSysPtr != nullptr) return qedTrialSysPtr->isInitial();
     else return false;}
   // Generate kinematics and check veto.
   bool acceptTrial(Event &event) override;
   // Update event after branching accepted.
   void updateEvent(Event &event) override;
   // Update partonSystems after branching accepted.
   void updatePartonSystems(Event &event) override;
   // Return scales.
   double q2minColoured() override {return q2minColouredSav;}
   double q2min() override {return q2minSav;}
 
   // Getter for number of systems.
   unsigned int nBranchers() override {
     int sizeNow = emitSystems.size();
     sizeNow = max(sizeNow, (int)splitSystems.size());
     sizeNow = max(sizeNow, (int)convSystems.size());
     return sizeNow;}
 
   // This module does not implement resonance decays.
   unsigned int nResDec() override { return 0; }
 
 private:
 
   // Get Q2 for QED system.
   template <class T>
   void q2NextSystem(map<int, T>& QEDsystemList, Event& event, double q2Start);
 
   // Systems.
   QEDemitSystem emptyQEDemitSystem;
   QEDsplitSystem emptyQEDsplitSystem;
   QEDconvSystem emptyQEDconvSystem;
   map< int, QEDemitSystem> emitSystems;
   map< int, QEDsplitSystem> splitSystems;
   map< int, QEDconvSystem> convSystems;
 
   // Settings.
   bool doQED, doEmission;
   int  nGammaToLepton, nGammaToQuark;
   bool doConvertGamma, doConvertQuark;
 
   // Scales.
   double q2minSav, q2minColouredSav;
 
   // Trial information
   int iSysTrial;
   double q2Trial;
   QEDsystem* qedTrialSysPtr{};
 
   // AlphaEM.
   AlphaEM al;
 
   // Evolution windows
   vector<double> evolutionWindows;
 
 };
 
 //==========================================================================
 
 } // end namespace Pythia8
 
 #endif // Pythia8_VinciaQED_H

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