include/Pythia8/MultipartonInteractions.h

0001 // MultipartonInteractions.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 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(),
0108     alphaEMorder(), alphaSnfmax(), bProfile(), processLevel(), bSelScale(),
0109     rescatterMode(), nQuarkIn(), nSample(), enhanceScreening(), pT0paramMode(),
0110     alphaSvalue(), Kfactor(), pT0Ref(), ecmRef(), ecmPow(), pTmin(),
0111     coreRadius(), coreFraction(), expPow(), ySepResc(), deltaYResc(),
0112     sigmaPomP(), mPomP(), pPomP(), 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          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   bool   hasPomeronBeams, hasLowPow, globalRecoilFSR;
0267   int    iDiffSys, nMaxGlobalRecoilFSR, bSelHard;
0268   double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
0269          pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
0270          pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101],
0271          zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
0272          probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
0273          fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax,
0274          enhanceBavg;
0275
0276   // Properties specific to current system.
0277   bool   bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
0278   int    id1, id2, i1Sel, i2Sel, id1Sel, id2Sel, iPDFA, nPDFA;
0279   vector<int> idAList;
0280   double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
0281          tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
0282          dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
0283
0284   // Local values for beam particle switch and mass interpolation.
0285   int    iPDFAsave, nStep, iStepFrom, iStepTo;
0286   double eCMsave, eStepMin, eStepMax, eStepSize, eStepMix, eStepFrom, eStepTo;
0287
0288   // Stored values for mass interpolation. First index projectile particle.
0289   struct MPIInterpolationInfo {
0290     int nStepSave;
0291     double eStepMinSave, eStepMaxSave, eStepSizeSave;
0292     vector<double> pT0Save, pT4dSigmaMaxSave, pT4dProbMaxSave,
0293          sigmaIntSave, zeroIntCorrSave, normOverlapSave, kNowSave,
0294          bAvgSave, bDivSave, probLowBSave,
0295          fracAhighSave, fracBhighSave, fracChighSave,
0296          fracABChighSave, cDivSave, cMaxSave;
0297     vector<array<double, 101> >  sudExpPTSave;
0298
0299     void init(int nStepIn);
0300   };
0301
0302   vector<MPIInterpolationInfo> mpis;
0303
0304   // Beam offset wrt. normal situation and other photon-related parameters.
0305   int    beamOffset;
0306   double mGmGmMin, mGmGmMax;
0307   bool   hasGamma, isGammaGamma, isGammaHadron, isHadronGamma;
0308
0309   // Pointer to assign space-time vertices during parton evolution.
0310   PartonVertexPtr  partonVertexPtr;
0311
0312   // Collections of parton-level 2 -> 2 cross sections. Selected one.
0313   SigmaMultiparton  sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
0314   SigmaMultiparton* sigma2Sel;
0315   SigmaProcessPtr   dSigmaDtSel;
0316
0317   // Statistics on generated 2 -> 2 processes.
0318   map<int, int>  nGen;
0319
0320   // alphaStrong and alphaEM calculations.
0321   AlphaStrong    alphaS;
0322   AlphaEM        alphaEM;
0323
0324   // Scattered partons.
0325   vector<int>    scatteredA, scatteredB;
0326
0327   // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
0328   void upperEnvelope();
0329
0330   // Integrate the parton-parton interaction cross section.
0331   void jetCrossSection();
0332
0333   // Read or write initialization data from/to file, to save startup time.
0334   bool saveMPIdata();
0335   bool loadMPIdata();
0336
0337   // Evaluate "Sudakov form factor" for not having a harder interaction.
0338   double sudakov(double pT2sud, double enhance = 1.);
0339
0340   // Do a quick evolution towards the next smaller pT.
0341   double fastPT2( double pT2beg);
0342
0343   // Calculate the actual cross section, either for the first interaction
0344   // (including at initialization) or for any subsequent in the sequence.
0345   double sigmaPT2scatter(bool isFirst = false, bool doSymmetrize = false);
0346
0347   // Find the partons that may rescatter.
0348   void findScatteredPartons( Event& event);
0349
0350   // Calculate the actual cross section for a rescattering.
0351   double sigmaPT2rescatter( Event& event);
0352
0353   // Calculate factor relating matter overlap and interaction rate.
0354   void overlapInit();
0355
0356   // Pick impact parameter and interaction rate enhancement,
0357   // either before the first interaction (for nondiffractive) or after it.
0358   void overlapFirst();
0359   void overlapNext(Event& event, double pTscale, bool rehashB);
0360
0361 };
0362
0363 //==========================================================================
0364
0365 } // end namespace Pythia8
0366
0367 #endif // Pythia8_MultipartonInteractions_H