1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
|
// MultipartonInteractions.h is a part of the PYTHIA event generator.
// Copyright (C) 2012 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// This file contains the main classes for multiparton interactions physics.
// SigmaMultiparton stores allowed processes by in-flavour combination.
// MultipartonInteractions: generates multiparton parton-parton interactions.
#ifndef Pythia8_MultipartonInteractions_H
#define Pythia8_MultipartonInteractions_H
#include "Basics.h"
#include "BeamParticle.h"
#include "Event.h"
#include "Info.h"
#include "PartonSystems.h"
#include "PythiaStdlib.h"
#include "Settings.h"
#include "SigmaTotal.h"
#include "SigmaProcess.h"
#include "StandardModel.h"
#include "UserHooks.h"
namespace Pythia8 {
//==========================================================================
// SigmaMultiparton is a helper class to MultipartonInteractions.
// It packs pointers to the allowed processes for different
// flavour combinations and levels of ambition.
class SigmaMultiparton {
public:
// Constructor.
SigmaMultiparton() {}
// Destructor.
~SigmaMultiparton() {
for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i];
for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];}
// Initialize list of processes.
bool init(int inState, int processLevel, Info* infoPtr,
Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtrIn,
BeamParticle* beamAPtr, BeamParticle* beamBPtr, Couplings* couplingsPtr);
// Calculate cross section summed over possibilities.
double sigma( int id1, int id2, double x1, double x2, double sHat,
double tHat, double uHat, double alpS, double alpEM,
bool restore = false, bool pickOtherIn = false);
// Return whether the other, rare processes were selected.
bool pickedOther() {return pickOther;}
// Return one subprocess, picked according to relative cross sections.
SigmaProcess* sigmaSel();
bool swapTU() {return pickedU;}
// Return code or name of a specified process, for statistics table.
int nProc() const {return nChan;}
int codeProc(int iProc) const {return sigmaT[iProc]->code();}
string nameProc(int iProc) const {return sigmaT[iProc]->name();}
private:
// Constants: could only be changed in the code itself.
static const double MASSMARGIN, OTHERFRAC;
// Number of processes. Some use massive matrix elements.
int nChan;
vector<bool> needMasses;
vector<double> m3Fix, m4Fix, sHatMin;
// Vector of process list, one for t-channel and one for u-channel.
vector<SigmaProcess*> sigmaT, sigmaU;
// Values of cross sections in process list above.
vector<double> sigmaTval, sigmaUval;
double sigmaTsum, sigmaUsum;
bool pickOther, pickedU;
// Pointer to the random number generator.
Rndm* rndmPtr;
};
//==========================================================================
// The MultipartonInteractions class contains the main methods for the
// generation of multiparton parton-parton interactions in hadronic collisions.
class MultipartonInteractions {
public:
// Constructor.
MultipartonInteractions() : bIsSet(false) {}
// Initialize the generation process for given beams.
bool init( bool doMPIinit, int diffractiveModeIn, Info* infoPtrIn,
Settings& settings, ParticleData* particleDataPtr, Rndm* rndmPtrIn,
BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn,
ostream& os = cout);
// Reset impact parameter choice and update the CM energy.
void reset();
// Select first = hardest pT in minbias process.
void pTfirst();
// Set up kinematics for first = hardest pT in minbias process.
void setupFirstSys( Event& process);
// Find whether to limit maximum scale of emissions.
bool limitPTmax( Event& event);
// Prepare system for evolution.
void prepare(Event& event, double pTscale = 1000.) {
if (!bSetInFirst) overlapNext(event, pTscale); }
// Select next pT in downwards evolution.
double pTnext( double pTbegAll, double pTendAll, Event& event);
// Set up kinematics of acceptable interaction.
bool scatter( Event& event);
// Set "empty" values to avoid query of undefined quantities.
void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac
= xPDF1now = xPDF2now = 0.; bIsSet = false;}
// Get some information on current interaction.
double Q2Ren() const {return pT2Ren;}
double alphaSH() const {return alpS;}
double alphaEMH() const {return alpEM;}
double x1H() const {return x1;}
double x2H() const {return x2;}
double Q2Fac() const {return pT2Fac;}
double pdf1() const {return xPDF1now;}
double pdf2() const {return xPDF2now;}
double bMPI() const {return (bIsSet) ? bNow / bAvg : 0.;}
double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
// For x-dependent matter profile, return incoming valence/sea
// decision from trial interactions.
int getVSC1() const {return vsc1;}
int getVSC2() const {return vsc2;}
// Update and print statistics on number of processes.
void accumulate() { int iBeg = (infoPtr->isMinBias()) ? 0 : 1;
for (int i = iBeg; i < infoPtr->nMPI(); ++i)
++nGen[ infoPtr->codeMPI(i) ];}
void statistics(bool resetStat = false, ostream& os = cout);
void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin();
iter != nGen.end(); ++iter) iter->second = 0; }
private:
// Constants: could only be changed in the code itself.
static const bool SHIFTFACSCALE, PREPICKRESCATTER;
static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN,
EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX,
KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV;
// Initialization data, read from Settings.
bool allowRescatter, allowDoubleRes, canVetoMPI;
int pTmaxMatch, alphaSorder, alphaEMorder, bProfile, processLevel,
rescatterMode, nQuarkIn, nSample, enhanceScreening;
double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius,
coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP,
mMaxPertDiff, mMinPertDiff;
// x-dependent matter profile:
// Constants.
static const int XDEP_BBIN;
static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC,
XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM;
// Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast
// integration. Only needed during initialisation.
vector <double> sigmaIntWgt, sigmaSumWgt;
// a1: value of a1 constant, taken from settings database.
// a0now (a02now): tuned value of a0 (squared value).
// bstepNow: current size of bins in b.
// a2max: given an xMin, a maximal (squared) value of the
// width, to be used in overestimate Omax(b).
// enhanceBmax, retain enhanceB as enhancement factor for the hardest
// enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2,
// and enhanceBnow to store the value for the current
// interaction.
double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow;
// Storage for trial interactions.
int id1Save, id2Save;
double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave,
alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave,
xPDF2nowSave;
SigmaProcess *dSigmaDtSelSave;
// vsc1, vsc2: for minimum-bias events with trial interaction, store
// decision on whether hard interaction was valence or sea.
int vsc1, vsc2;
// Other initialization data.
bool hasBaryonBeams, hasLowPow, globalRecoilFSR;
int diffractiveMode, nMaxGlobalRecoilFSR;
double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101],
zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax;
// Properties specific to current system.
bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel;
double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
// Stored values for mass interpolation for diffractive systems.
int nStep, iStepFrom, iStepTo;
double eCMsave, eStepSize, eStepSave, eStepFrom, eStepTo, pT0Save[5],
pT4dSigmaMaxSave[5], pT4dProbMaxSave[5], sigmaIntSave[5],
sudExpPTSave[5][101], zeroIntCorrSave[5], normOverlapSave[5],
kNowSave[5], bAvgSave[5], bDivSave[5], probLowBSave[5],
fracAhighSave[5], fracBhighSave[5], fracChighSave[5],
fracABChighSave[5], cDivSave[5], cMaxSave[5];
// Pointer to various information on the generation.
Info* infoPtr;
// Pointer to the random number generator.
Rndm* rndmPtr;
// Pointers to the two incoming beams.
BeamParticle* beamAPtr;
BeamParticle* beamBPtr;
// Pointers to Standard Model couplings.
Couplings* couplingsPtr;
// Pointer to information on subcollision parton locations.
PartonSystems* partonSystemsPtr;
// Pointer to total cross section parametrization.
SigmaTotal* sigmaTotPtr;
// Pointer to user hooks.
UserHooks* userHooksPtr;
// Collections of parton-level 2 -> 2 cross sections. Selected one.
SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
SigmaMultiparton* sigma2Sel;
SigmaProcess* dSigmaDtSel;
// Statistics on generated 2 -> 2 processes.
map<int, int> nGen;
// alphaStrong and alphaEM calculations.
AlphaStrong alphaS;
AlphaEM alphaEM;
// Scattered partons.
vector<int> scatteredA, scatteredB;
// Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
void upperEnvelope();
// Integrate the parton-parton interaction cross section.
void jetCrossSection();
// Evaluate "Sudakov form factor" for not having a harder interaction.
double sudakov(double pT2sud, double enhance = 1.);
// Do a quick evolution towards the next smaller pT.
double fastPT2( double pT2beg);
// Calculate the actual cross section, either for the first interaction
// (including at initialization) or for any subsequent in the sequence.
double sigmaPT2scatter(bool isFirst = false);
// Find the partons that may rescatter.
void findScatteredPartons( Event& event);
// Calculate the actual cross section for a rescattering.
double sigmaPT2rescatter( Event& event);
// Calculate factor relating matter overlap and interaction rate.
void overlapInit();
// Pick impact parameter and interaction rate enhancement,
// either before the first interaction (for minbias) or after it.
void overlapFirst();
void overlapNext(Event& event, double pTscale);
};
//==========================================================================
} // end namespace Pythia8
#endif // Pythia8_MultipartonInteractions_H
|