include/Pythia8/MultipartonInteractions.h

0001 // MultipartonInteractions.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 main classes for multiparton interactions physics.
0007 // SigmaMultiparton stores allowed processes by in-flavour combination.
0008 // MultipartonInteractions: generates multiparton parton-parton interactions.
0009
0010 #ifndef Pythia8_MultipartonInteractions_H
0011 #define Pythia8_MultipartonInteractions_H
0012
0013 #include "Pythia8/Basics.h"
0014 #include "Pythia8/BeamParticle.h"
0015 #include "Pythia8/Event.h"
0016 #include "Pythia8/Info.h"
0017 #include "Pythia8/PartonSystems.h"
0018 #include "Pythia8/PartonVertex.h"
0019 #include "Pythia8/PhysicsBase.h"
0020 #include "Pythia8/PythiaStdlib.h"
0021 #include "Pythia8/Settings.h"
0022 #include "Pythia8/SigmaTotal.h"
0023 #include "Pythia8/SigmaProcess.h"
0024 #include "Pythia8/StandardModel.h"
0025 #include "Pythia8/UserHooks.h"
0026
0027 namespace Pythia8 {
0028
0029 //==========================================================================
0030
0031 // SigmaMultiparton is a helper class to MultipartonInteractions.
0032 // It packs pointers to the allowed processes for different
0033 // flavour combinations and levels of ambition.
0034
0035 class SigmaMultiparton {
0036
0037 public:
0038
0039   // Constructor.
0040   SigmaMultiparton() : nChan(), needMasses(), useNarrowBW3(), useNarrowBW4(),
0041     m3Fix(), m4Fix(), sHatMin(), sigmaT(), sigmaU(), sigmaTval(), sigmaUval(),
0042     sigmaTsum(), sigmaUsum(), pickOther(), pickedU(), particleDataPtr(),
0043     rndmPtr() {}
0044
0045   // Initialize list of processes.
0046   bool init(int inState, int processLevel, Info* infoPtr,
0047     BeamParticle* beamAPtr, BeamParticle* beamBPtr);
0048
0049   // Switch to new beam particle identities.
0050   void updateBeamIDs() {
0051     for (int i = 0; i < nChan; ++i) sigmaT[i]->updateBeamIDs();
0052     for (int i = 0; i < nChan; ++i) sigmaU[i]->updateBeamIDs(); }
0053
0054   // Calculate cross section summed over possibilities.
0055   double sigma( int id1, int id2, double x1, double x2, double sHat,
0056     double tHat, double uHat, double alpS, double alpEM,
0057     bool restore = false, bool pickOtherIn = false);
0058
0059   // Return whether the other, rare processes were selected.
0060   bool pickedOther() {return pickOther;}
0061
0062   // Return one subprocess, picked according to relative cross sections.
0063   SigmaProcessPtr sigmaSel();
0064   bool swapTU() {return pickedU;}
0065
0066   // Return code or name of a specified process, for statistics table.
0067   int    nProc() const {return nChan;}
0068   int    codeProc(int iProc) const {return sigmaT[iProc]->code();}
0069   string nameProc(int iProc) const {return sigmaT[iProc]->name();}
0070
0071 private:
0072
0073   // Constants: could only be changed in the code itself.
0074   static const double MASSMARGIN, OTHERFRAC;
0075
0076   // Number of processes. Some use massive matrix elements.
0077   int            nChan;
0078   vector<bool>   needMasses, useNarrowBW3, useNarrowBW4;
0079   vector<double> m3Fix, m4Fix, sHatMin;
0080
0081   // Vector of process list, one for t-channel and one for u-channel.
0082   vector<SigmaProcessPtr> sigmaT, sigmaU;
0083
0084   // Values of cross sections in process list above.
0085   vector<double> sigmaTval, sigmaUval;
0086   double         sigmaTsum, sigmaUsum;
0087   bool           pickOther, pickedU;
0088
0089   // Pointers to the particle data and random number generator.
0090   ParticleData*  particleDataPtr;
0091   Rndm*          rndmPtr;
0092
0093 };
0094
0095 //==========================================================================
0096
0097 // The MultipartonInteractions class contains the main methods for the
0098 // generation of multiparton parton-parton interactions in hadronic collisions.
0099
0100 class MultipartonInteractions : public PhysicsBase {
0101
0102 public:
0103
0104   // Constructor.
0105   MultipartonInteractions() : allowRescatter(), allowDoubleRes(), canVetoMPI(),
0106     doPartonVertex(), doVarEcm(), setAntiSame(), setAntiSameNow(),
0107     pTmaxMatch(), alphaSorder(), alphaEMorder(), alphaSnfmax(), bProfile(),
0108     processLevel(), bSelScale(), rescatterMode(), nQuarkIn(), nSample(),
0109     enhanceScreening(), pT0paramMode(), alphaSvalue(), Kfactor(), pT0Ref(),
0110     ecmRef(), ecmPow(), pTmin(), coreRadius(), coreFraction(), expPow(),
0111     ySepResc(), deltaYResc(), sigmaPomP(), mPomP(), pPomP(), sigmaPomPom(),
0112     mMaxPertDiff(), mMinPertDiff(), a1(),
0113     a0now(), a02now(), bstepNow(), a2max(), b2now(), enhanceBmax(),
0114     enhanceBnow(), id1Save(), id2Save(), pT2Save(), x1Save(), x2Save(),
0115     sHatSave(), tHatSave(), uHatSave(), alpSsave(), alpEMsave(), pT2FacSave(),
0116     pT2RenSave(), xPDF1nowSave(), xPDF2nowSave(), scaleLimitPTsave(),
0117     dSigmaDtSelSave(), vsc1(), vsc2(), hasPomeronBeams(), hasLowPow(),
0118     globalRecoilFSR(), iDiffSys(), nMaxGlobalRecoilFSR(), bSelHard(), eCM(),
0119     sCM(), pT0(), pT20(), pT2min(), pTmax(), pT2max(), pT20R(), pT20minR(),
0120     pT20maxR(), pT20min0maxR(), pT2maxmin(), sigmaND(), pT4dSigmaMax(),
0121     pT4dProbMax(), dSigmaApprox(), sigmaInt(), sudExpPT(), zeroIntCorr(),
0122     normOverlap(), nAvg(), kNow(), normPi(), bAvg(), bDiv(), probLowB(),
0123     radius2B(), radius2C(), fracA(), fracB(), fracC(), fracAhigh(),
0124     fracBhigh(), fracChigh(), fracABChigh(), expRev(), cDiv(), cMax(),
0125     enhanceBavg(), bIsSet(false), bSetInFirst(), isAtLowB(), pickOtherSel(),
0126     id1(), id2(), i1Sel(), i2Sel(), id1Sel(), id2Sel(), iPDFA(), nPDFA(1),
0127     bNow(), enhanceB(), pT2(), pT2shift(), pT2Ren(), pT2Fac(), x1(),
0128     x2(), xT(), xT2(), tau(), y(), sHat(), tHat(), uHat(), alpS(), alpEM(),
0129     xPDF1now(), xPDF2now(), dSigmaSum(), x1Sel(), x2Sel(), sHatSel(),
0130     tHatSel(), uHatSel(), iPDFAsave(), nStep(), iStepFrom(), iStepTo(),
0131     eCMsave(), eStepMin(), eStepMax(), eStepSize(), eStepMix(), eStepFrom(),
0132     eStepTo(), beamOffset(), mGmGmMin(), mGmGmMax(), hasGamma(),
0133     isGammaGamma(), isGammaHadron(), isHadronGamma(), partonVertexPtr(),
0134     sigma2Sel(), dSigmaDtSel() {}
0135
0136   // Initialize the generation process for given beams.
0137   bool init( bool doMPIinit, int iDiffSysIn,
0138     BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
0139     PartonVertexPtr partonVertexPtrIn, bool hasGammaIn = false);
0140
0141   // Special setup to allow switching between beam PDFs for MPI handling.
0142   void initSwitchID( const vector<int>& idAListIn) {
0143     idAList = idAListIn; nPDFA = idAList.size();
0144     mpis  = vector<MPIInterpolationInfo>(nPDFA);}
0145
0146     // Switch to new beam particle identities, and possibly PDFs.
0147   void setBeamID(int iPDFAin) { iPDFA = iPDFAin;
0148     sigma2gg.updateBeamIDs(); sigma2qg.updateBeamIDs();
0149     sigma2qqbarSame.updateBeamIDs(); sigma2qq.updateBeamIDs();
0150     setAntiSameNow = setAntiSame && particleDataPtr->hasAnti(infoPtr->idA())
0151       && particleDataPtr->hasAnti(infoPtr->idB());}
0152
0153   // Reset impact parameter choice and update the CM energy.
0154   void reset();
0155
0156   // Select first = hardest pT in nondiffractive process.
0157   void pTfirst();
0158
0159   // Set up kinematics for first = hardest pT in nondiffractive process.
0160   void setupFirstSys( Event& process);
0161
0162   // Find whether to limit maximum scale of emissions.
0163   // Provide sum pT / 2 as potential limit where relevant.
0164   bool limitPTmax( Event& event);
0165   double scaleLimitPT() const {return scaleLimitPTsave;}
0166
0167   // Prepare system for evolution.
0168   void prepare(Event& event, double pTscale = 1000., bool rehashB = false) {
0169     if (!bSetInFirst) overlapNext(event, pTscale, rehashB); }
0170
0171   // Select next pT in downwards evolution.
0172   double pTnext( double pTbegAll, double pTendAll, Event& event);
0173
0174   // Set up kinematics of acceptable interaction.
0175   bool scatter( Event& event);
0176
0177   // Set "empty" values to avoid query of undefined quantities.
0178   void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac
0179     = xPDF1now = xPDF2now = 0.; bIsSet = false;}
0180
0181   // Get some information on current interaction.
0182   double Q2Ren()      const {return pT2Ren;}
0183   double alphaSH()    const {return alpS;}
0184   double alphaEMH()   const {return alpEM;}
0185   double x1H()        const {return x1;}
0186   double x2H()        const {return x2;}
0187   double Q2Fac()      const {return pT2Fac;}
0188   double pdf1()       const {return (id1 == 21 ? 4./9. : 1.) * xPDF1now;}
0189   double pdf2()       const {return (id2 == 21 ? 4./9. : 1.) * xPDF2now;}
0190   double bMPI()       const {return (bIsSet) ? bNow : 0.;}
0191   double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
0192   double enhanceMPIavg() const {return enhanceBavg;}
0193
0194   // For x-dependent matter profile, return incoming valence/sea
0195   // decision from trial interactions.
0196   int    getVSC1()   const {return vsc1;}
0197   int    getVSC2()   const {return vsc2;}
0198
0199   // Set the offset wrt. to normal beam particle positions for hard diffraction
0200   // and for photon beam from lepton.
0201   int  getBeamOffset()       const {return beamOffset;}
0202   void setBeamOffset(int offsetIn) {beamOffset = offsetIn;}
0203
0204   // Update and print statistics on number of processes.
0205   // Note: currently only valid for nondiffractive systems, not diffraction??
0206   void accumulate() { int iBeg = (infoPtr->isNonDiffractive()) ? 0 : 1;
0207     for (int i = iBeg; i < infoPtr->nMPI(); ++i)
0208     ++nGen[ infoPtr->codeMPI(i) ];}
0209   void statistics(bool resetStat = false);
0210   void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin();
0211     iter != nGen.end(); ++iter) iter->second = 0; }
0212
0213 private:
0214
0215   // Constants: could only be changed in the code itself.
0216   static const bool   SHIFTFACSCALE, PREPICKRESCATTER;
0217   static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN,
0218                       EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX,
0219                       KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV, WTACCWARN,
0220                       SIGMAMBLIMIT;
0221
0222   // Initialization data, read from Settings.
0223   bool   allowRescatter, allowDoubleRes, canVetoMPI, doPartonVertex, doVarEcm,
0224          setAntiSame, setAntiSameNow, allowIDAswitch;
0225   int    pTmaxMatch, alphaSorder, alphaEMorder, alphaSnfmax, bProfile,
0226          processLevel, bSelScale, rescatterMode, nQuarkIn, nSample,
0227          enhanceScreening, pT0paramMode, reuseInit;
0228   double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius,
0229          coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP,
0230          sigmaPomPom, mMaxPertDiff, mMinPertDiff;
0231   string initFile;
0232
0233   // x-dependent matter profile:
0234   // Constants.
0235   static const int    XDEP_BBIN;
0236   static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC,
0237                       XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM;
0238
0239   // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast
0240   // integration. Only needed during initialisation.
0241   vector <double> sigmaIntWgt, sigmaSumWgt;
0242
0243   // a1:             value of a1 constant, taken from settings database.
0244   // a0now (a02now): tuned value of a0 (squared value).
0245   // bstepNow:       current size of bins in b.
0246   // a2max:          given an xMin, a maximal (squared) value of the
0247   //                 width, to be used in overestimate Omax(b).
0248   // enhanceBmax,    retain enhanceB as enhancement factor for the hardest
0249   // enhanceBnow:    interaction. Use enhanceBmax as overestimate for fastPT2,
0250   //                 and enhanceBnow to store the value for the current
0251   //                 interaction.
0252   double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow;
0253
0254   // Storage for trial interactions.
0255   int    id1Save, id2Save;
0256   double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave,
0257          alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave,
0258          xPDF2nowSave, scaleLimitPTsave;
0259   SigmaProcessPtr dSigmaDtSelSave;
0260
0261   // vsc1, vsc2:     for minimum-bias events with trial interaction, store
0262   //                 decision on whether hard interaction was valence or sea.
0263   int    vsc1, vsc2;
0264
0265   // Other initialization data.
0266   // Number of points in the interpolation of Sudakov exponents.
0267   static const int NSUDPTS = 50;
0268   friend class MPIInterpolationInfo;
0269   bool   hasPomeronBeams, hasLowPow, globalRecoilFSR;
0270   int    iDiffSys, nMaxGlobalRecoilFSR, bSelHard;
0271   double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
0272          pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
0273          pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[NSUDPTS + 1],
0274          zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
0275          probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
0276          fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax,
0277          enhanceBavg;
0278
0279   // Properties specific to current system.
0280   bool   bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
0281   int    id1, id2, i1Sel, i2Sel, id1Sel, id2Sel, iPDFA, nPDFA;
0282   vector<int> idAList;
0283   double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
0284          tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
0285          dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
0286
0287   // Local values for beam particle switch and mass interpolation.
0288   int    iPDFAsave, nStep, iStepFrom, iStepTo;
0289   double eCMsave, eStepMin, eStepMax, eStepSize, eStepMix, eStepFrom, eStepTo;
0290
0291   // Stored values for mass interpolation. First index projectile particle.
0292   struct MPIInterpolationInfo {
0293     int nStepSave;
0294     double eStepMinSave, eStepMaxSave, eStepSizeSave;
0295     vector<double> pT0Save, pT4dSigmaMaxSave, pT4dProbMaxSave,
0296          sigmaIntSave, zeroIntCorrSave, normOverlapSave, kNowSave,
0297          bAvgSave, bDivSave, probLowBSave,
0298          fracAhighSave, fracBhighSave, fracChighSave,
0299          fracABChighSave, cDivSave, cMaxSave;
0300     vector<array<double, MultipartonInteractions::NSUDPTS + 1> >  sudExpPTSave;
0301
0302     void init(int nStepIn);
0303   };
0304
0305   vector<MPIInterpolationInfo> mpis;
0306
0307   // Beam offset wrt. normal situation and other photon-related parameters.
0308   int    beamOffset;
0309   double mGmGmMin, mGmGmMax;
0310   bool   hasGamma, isGammaGamma, isGammaHadron, isHadronGamma;
0311
0312   // Pointer to assign space-time vertices during parton evolution.
0313   PartonVertexPtr  partonVertexPtr;
0314
0315   // Collections of parton-level 2 -> 2 cross sections. Selected one.
0316   SigmaMultiparton  sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
0317   SigmaMultiparton* sigma2Sel;
0318   SigmaProcessPtr   dSigmaDtSel;
0319
0320   // Statistics on generated 2 -> 2 processes.
0321   map<int, int>  nGen;
0322
0323   // alphaStrong and alphaEM calculations.
0324   AlphaStrong    alphaS;
0325   AlphaEM        alphaEM;
0326
0327   // Scattered partons.
0328   vector<int>    scatteredA, scatteredB;
0329
0330   // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
0331   void upperEnvelope();
0332
0333   // Integrate the parton-parton interaction cross section.
0334   void jetCrossSection();
0335
0336   // Read or write initialization data from/to settings/file, to save
0337   // startup time.
0338   bool saveMPIdata();
0339   bool loadMPIdata();
0340
0341   // Evaluate "Sudakov form factor" for not having a harder interaction.
0342   double sudakov(double pT2sud, double enhance = 1.);
0343
0344   // Do a quick evolution towards the next smaller pT.
0345   double fastPT2( double pT2beg);
0346
0347   // Calculate the actual cross section, either for the first interaction
0348   // (including at initialization) or for any subsequent in the sequence.
0349   double sigmaPT2scatter(bool isFirst = false, bool doSymmetrize = false);
0350
0351   // Find the partons that may rescatter.
0352   void findScatteredPartons( Event& event);
0353
0354   // Calculate the actual cross section for a rescattering.
0355   double sigmaPT2rescatter( Event& event);
0356
0357   // Calculate factor relating matter overlap and interaction rate.
0358   void overlapInit();
0359
0360   // Pick impact parameter and interaction rate enhancement,
0361   // either before the first interaction (for nondiffractive) or after it.
0362   void overlapFirst();
0363   void overlapNext(Event& event, double pTscale, bool rehashB);
0364
0365 };
0366
0367 //==========================================================================
0368
0369 } // end namespace Pythia8
0370
0371 #endif // Pythia8_MultipartonInteractions_H