1 | // BeamParticle.h is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Header file for information on incoming beams. |
---|
7 | // ResolvedParton: an initiator or remnant in beam. |
---|
8 | // BeamParticle: contains partons, parton densities, etc. |
---|
9 | |
---|
10 | #ifndef Pythia8_BeamParticle_H |
---|
11 | #define Pythia8_BeamParticle_H |
---|
12 | |
---|
13 | #include "Basics.h" |
---|
14 | #include "Event.h" |
---|
15 | #include "FragmentationFlavZpT.h" |
---|
16 | #include "Info.h" |
---|
17 | #include "ParticleData.h" |
---|
18 | #include "PartonDistributions.h" |
---|
19 | #include "PythiaStdlib.h" |
---|
20 | #include "Settings.h" |
---|
21 | |
---|
22 | namespace Pythia8 { |
---|
23 | |
---|
24 | //========================================================================== |
---|
25 | |
---|
26 | // This class holds info on a parton resolved inside the incoming beam, |
---|
27 | // i.e. either an initiator (part of a hard or a multiparton interaction) |
---|
28 | // or a remnant (part of the beam remnant treatment). |
---|
29 | |
---|
30 | // The companion code is -1 from onset and for g, is -2 for an unmatched |
---|
31 | // sea quark, is >= 0 for a matched sea quark, with the number giving the |
---|
32 | // companion position, and is -3 for a valence quark. |
---|
33 | |
---|
34 | // Rescattering partons properly do not belong here, but bookkeeping is |
---|
35 | // simpler with them, so they are stored with companion code -10. |
---|
36 | |
---|
37 | class ResolvedParton { |
---|
38 | |
---|
39 | public: |
---|
40 | |
---|
41 | // Constructor. |
---|
42 | ResolvedParton( int iPosIn = 0, int idIn = 0, double xIn = 0., |
---|
43 | int companionIn = -1) : iPosRes(iPosIn), idRes(idIn), xRes(xIn), |
---|
44 | companionRes(companionIn), xqCompRes(0.), mRes(0.), factorRes(1.), |
---|
45 | colRes(0), acolRes(0) { } |
---|
46 | |
---|
47 | // Set info on initiator or remnant parton. |
---|
48 | void iPos( int iPosIn) {iPosRes = iPosIn;} |
---|
49 | void id( int idIn) {idRes = idIn;} |
---|
50 | void x( double xIn) {xRes = xIn;} |
---|
51 | void update( int iPosIn, int idIn, double xIn) {iPosRes = iPosIn; |
---|
52 | idRes = idIn; xRes = xIn;} |
---|
53 | void companion( int companionIn) {companionRes = companionIn;} |
---|
54 | void xqCompanion( double xqCompIn) {xqCompRes = xqCompIn;} |
---|
55 | void p(Vec4 pIn) {pRes = pIn;} |
---|
56 | void px(double pxIn) {pRes.px(pxIn);} |
---|
57 | void py(double pyIn) {pRes.py(pyIn);} |
---|
58 | void pz(double pzIn) {pRes.pz(pzIn);} |
---|
59 | void e(double eIn) {pRes.e(eIn);} |
---|
60 | void m(double mIn) {mRes = mIn;} |
---|
61 | void col(int colIn) {colRes = colIn;} |
---|
62 | void acol(int acolIn) {acolRes = acolIn;} |
---|
63 | void cols(int colIn = 0,int acolIn = 0) |
---|
64 | {colRes = colIn; acolRes = acolIn;} |
---|
65 | void scalePT( double factorIn) {pRes.px(factorIn * pRes.px()); |
---|
66 | pRes.py(factorIn * pRes.py()); factorRes *= factorIn;} |
---|
67 | void scaleX( double factorIn) {xRes *= factorIn;} |
---|
68 | |
---|
69 | // Get info on initiator or remnant parton. |
---|
70 | int iPos() const {return iPosRes;} |
---|
71 | int id() const {return idRes;} |
---|
72 | double x() const {return xRes;} |
---|
73 | int companion() const {return companionRes;} |
---|
74 | bool isValence() const {return (companionRes == -3);} |
---|
75 | bool isUnmatched() const {return (companionRes == -2);} |
---|
76 | bool isCompanion() const {return (companionRes >= 0);} |
---|
77 | bool isFromBeam() const {return (companionRes > -10);} |
---|
78 | double xqCompanion() const {return xqCompRes;} |
---|
79 | Vec4 p() const {return pRes;} |
---|
80 | double px() const {return pRes.px();} |
---|
81 | double py() const {return pRes.py();} |
---|
82 | double pz() const {return pRes.pz();} |
---|
83 | double e() const {return pRes.e();} |
---|
84 | double m() const {return mRes;} |
---|
85 | double pT() const {return pRes.pT();} |
---|
86 | double mT2() const {return (mRes >= 0.) |
---|
87 | ? mRes*mRes + pRes.pT2() : - mRes*mRes + pRes.pT2();} |
---|
88 | double pPos() const {return pRes.e() + pRes.pz();} |
---|
89 | double pNeg() const {return pRes.e() - pRes.pz();} |
---|
90 | int col() const {return colRes;} |
---|
91 | int acol() const {return acolRes;} |
---|
92 | double pTfactor() const {return factorRes;} |
---|
93 | |
---|
94 | private: |
---|
95 | |
---|
96 | // Properties of a resolved parton. |
---|
97 | int iPosRes, idRes; |
---|
98 | double xRes; |
---|
99 | // Companion code and distribution value, if any. |
---|
100 | int companionRes; |
---|
101 | double xqCompRes; |
---|
102 | // Four-momentum and mass; for remnant kinematics construction. |
---|
103 | Vec4 pRes; |
---|
104 | double mRes, factorRes; |
---|
105 | // Colour codes. |
---|
106 | int colRes, acolRes; |
---|
107 | |
---|
108 | }; |
---|
109 | |
---|
110 | //========================================================================== |
---|
111 | |
---|
112 | // This class holds info on a beam particle in the evolution of |
---|
113 | // initial-state radiation and multiparton interactions. |
---|
114 | |
---|
115 | class BeamParticle { |
---|
116 | |
---|
117 | public: |
---|
118 | |
---|
119 | // Constructor. |
---|
120 | BeamParticle() : nInit(0) {resolved.resize(0); Q2ValFracSav = -1.;} |
---|
121 | |
---|
122 | // Initialize data on a beam particle and save pointers. |
---|
123 | void init( int idIn, double pzIn, double eIn, double mIn, |
---|
124 | Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn, |
---|
125 | Rndm* rndmPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn, |
---|
126 | StringFlav* flavSelPtrIn); |
---|
127 | |
---|
128 | // For mesons like pi0 valence content varies from event to event. |
---|
129 | void newValenceContent(); |
---|
130 | |
---|
131 | // Set new pZ and E, but keep the rest the same. |
---|
132 | void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);} |
---|
133 | |
---|
134 | // Member functions for output. |
---|
135 | int id() const {return idBeam;} |
---|
136 | Vec4 p() const {return pBeam;} |
---|
137 | double px() const {return pBeam.px();} |
---|
138 | double py() const {return pBeam.py();} |
---|
139 | double pz() const {return pBeam.pz();} |
---|
140 | double e() const {return pBeam.e();} |
---|
141 | double m() const {return mBeam;} |
---|
142 | bool isLepton() const {return isLeptonBeam;} |
---|
143 | bool isUnresolved() const {return isUnresolvedBeam;} |
---|
144 | // As hadrons here we only count those we know how to handle remnants for. |
---|
145 | bool isHadron() const {return isHadronBeam;} |
---|
146 | bool isMeson() const {return isMesonBeam;} |
---|
147 | bool isBaryon() const {return isBaryonBeam;} |
---|
148 | |
---|
149 | // Maximum x remaining after previous MPI and ISR, plus safety margin. |
---|
150 | double xMax(int iSkip = -1); |
---|
151 | |
---|
152 | // Special hard-process parton distributions (can agree with standard ones). |
---|
153 | double xfHard(int idIn, double x, double Q2) |
---|
154 | {return pdfHardBeamPtr->xf(idIn, x, Q2);} |
---|
155 | |
---|
156 | // Standard parton distributions. |
---|
157 | double xf(int idIn, double x, double Q2) |
---|
158 | {return pdfBeamPtr->xf(idIn, x, Q2);} |
---|
159 | |
---|
160 | // Ditto, split into valence and sea parts (where gluon counts as sea). |
---|
161 | double xfVal(int idIn, double x, double Q2) |
---|
162 | {return pdfBeamPtr->xfVal(idIn, x, Q2);} |
---|
163 | double xfSea(int idIn, double x, double Q2) |
---|
164 | {return pdfBeamPtr->xfSea(idIn, x, Q2);} |
---|
165 | |
---|
166 | // Rescaled parton distributions, as needed for MPI and ISR. |
---|
167 | // For ISR also allow split valence/sea, and only return relevant part. |
---|
168 | double xfMPI(int idIn, double x, double Q2) |
---|
169 | {return xfModified(-1, idIn, x, Q2);} |
---|
170 | double xfISR(int indexMPI, int idIn, double x, double Q2) |
---|
171 | {return xfModified( indexMPI, idIn, x, Q2);} |
---|
172 | |
---|
173 | // Decide whether chosen quark is valence, sea or companion. |
---|
174 | int pickValSeaComp(); |
---|
175 | |
---|
176 | // Initialize kind of incoming beam particle. |
---|
177 | void initBeamKind(); |
---|
178 | |
---|
179 | // Overload index operator to access a resolved parton from the list. |
---|
180 | ResolvedParton& operator[](int i) {return resolved[i];} |
---|
181 | |
---|
182 | // Total number of partons extracted from beam, and initiators only. |
---|
183 | int size() const {return resolved.size();} |
---|
184 | int sizeInit() const {return nInit;} |
---|
185 | |
---|
186 | // Clear list of resolved partons. |
---|
187 | void clear() {resolved.resize(0); nInit = 0;} |
---|
188 | |
---|
189 | // Add a resolved parton to list. |
---|
190 | int append( int iPos, int idIn, double x, int companion = -1) |
---|
191 | {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) ); |
---|
192 | return resolved.size() - 1;} |
---|
193 | |
---|
194 | // Print extracted parton list; for debug mainly. |
---|
195 | void list(ostream& os = cout) const; |
---|
196 | |
---|
197 | // How many different flavours, and how many quarks of given flavour. |
---|
198 | int nValenceKinds() const {return nValKinds;} |
---|
199 | int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i) |
---|
200 | if (idIn == idVal[i]) return nVal[i]; return 0;} |
---|
201 | |
---|
202 | // Test whether a lepton is to be considered as unresolved. |
---|
203 | bool isUnresolvedLepton(); |
---|
204 | |
---|
205 | // Add extra remnant flavours to make valence and sea come out right. |
---|
206 | bool remnantFlavours(Event& event); |
---|
207 | |
---|
208 | // Correlate all initiators and remnants to make a colour singlet. |
---|
209 | bool remnantColours(Event& event, vector<int>& colFrom, |
---|
210 | vector<int>& colTo); |
---|
211 | |
---|
212 | // Pick unrescaled x of remnant parton (valence or sea). |
---|
213 | double xRemnant(int i); |
---|
214 | |
---|
215 | // Tell whether a junction has been resolved, and its junction colours. |
---|
216 | bool hasJunction() const {return hasJunctionBeam;} |
---|
217 | int junctionCol(int i) const {return junCol[i];} |
---|
218 | void junctionCol(int i, int col) {junCol[i] = col;} |
---|
219 | |
---|
220 | // For a diffractive system, decide whether to kick out gluon or quark. |
---|
221 | bool pickGluon(double mDiff); |
---|
222 | |
---|
223 | // Pick a valence quark at random, and provide the remaining flavour. |
---|
224 | int pickValence(); |
---|
225 | int pickRemnant() const {return idVal2;} |
---|
226 | |
---|
227 | // Share lightcone momentum between two remnants in a diffractive system. |
---|
228 | // At the same time generate a relative pT for the two. |
---|
229 | double zShare( double mDiff, double m1, double m2); |
---|
230 | double pxShare() const {return pxRel;} |
---|
231 | double pyShare() const {return pyRel;} |
---|
232 | |
---|
233 | private: |
---|
234 | |
---|
235 | // Constants: could only be changed in the code itself. |
---|
236 | static const double XMINUNRESOLVED; |
---|
237 | |
---|
238 | // Pointer to various information on the generation. |
---|
239 | Info* infoPtr; |
---|
240 | |
---|
241 | // Pointer to the particle data table. |
---|
242 | ParticleData* particleDataPtr; |
---|
243 | |
---|
244 | // Pointer to the random number generator. |
---|
245 | Rndm* rndmPtr; |
---|
246 | |
---|
247 | // Pointers to PDF sets. |
---|
248 | PDF* pdfBeamPtr; |
---|
249 | PDF* pdfHardBeamPtr; |
---|
250 | |
---|
251 | // Pointer to class for flavour generation. |
---|
252 | StringFlav* flavSelPtr; |
---|
253 | |
---|
254 | // Initialization data, normally only set once. |
---|
255 | bool allowJunction; |
---|
256 | int maxValQuark, companionPower; |
---|
257 | double valencePowerMeson, valencePowerUinP, valencePowerDinP, |
---|
258 | valenceDiqEnhance, pickQuarkNorm, pickQuarkPower, |
---|
259 | diffPrimKTwidth, diffLargeMassSuppress; |
---|
260 | |
---|
261 | // Basic properties of a beam particle. |
---|
262 | int idBeam, idBeamAbs; |
---|
263 | Vec4 pBeam; |
---|
264 | double mBeam; |
---|
265 | // Beam kind. Valence flavour content for hadrons. |
---|
266 | bool isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam, |
---|
267 | isBaryonBeam; |
---|
268 | int nValKinds, idVal[3], nVal[3]; |
---|
269 | |
---|
270 | // Current parton density, by valence, sea and companion. |
---|
271 | int idSave, iSkipSave, nValLeft[3]; |
---|
272 | double xqgTot, xqVal, xqgSea, xqCompSum; |
---|
273 | |
---|
274 | // The list of resolved partons. |
---|
275 | vector<ResolvedParton> resolved; |
---|
276 | |
---|
277 | // Status after all initiators have been accounted for. Junction content. |
---|
278 | int nInit; |
---|
279 | bool hasJunctionBeam; |
---|
280 | int junCol[3]; |
---|
281 | |
---|
282 | // Routine to calculate pdf's given previous interactions. |
---|
283 | double xfModified( int iSkip, int idIn, double x, double Q2); |
---|
284 | |
---|
285 | // Fraction of hadron momentum sitting in a valence quark distribution. |
---|
286 | double xValFrac(int j, double Q2); |
---|
287 | double Q2ValFracSav, uValInt, dValInt; |
---|
288 | |
---|
289 | // Fraction of hadron momentum sitting in a companion quark distribution. |
---|
290 | double xCompFrac(double xs); |
---|
291 | |
---|
292 | // Value of companion quark PDF, also given the sea quark x. |
---|
293 | double xCompDist(double xc, double xs); |
---|
294 | |
---|
295 | // Valence quark subdivision for diffractive systems. |
---|
296 | int idVal1, idVal2, idVal3; |
---|
297 | double zRel, pxRel, pyRel; |
---|
298 | |
---|
299 | }; |
---|
300 | |
---|
301 | //========================================================================== |
---|
302 | |
---|
303 | } // end namespace Pythia8 |
---|
304 | |
---|
305 | #endif // Pythia8_BeamParticle_H |
---|