]>
Commit | Line | Data |
---|---|---|
63ba5337 | 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 |