Line data Source code
1 : // BeamParticle.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "Pythia8/Basics.h"
14 : #include "Pythia8/Event.h"
15 : #include "Pythia8/FragmentationFlavZpT.h"
16 : #include "Pythia8/Info.h"
17 : #include "Pythia8/ParticleData.h"
18 : #include "Pythia8/PartonDistributions.h"
19 : #include "Pythia8/PythiaStdlib.h"
20 : #include "Pythia8/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 0 : class ResolvedParton {
38 :
39 : public:
40 :
41 : // Constructor.
42 0 : ResolvedParton( int iPosIn = 0, int idIn = 0, double xIn = 0.,
43 0 : int companionIn = -1) : iPosRes(iPosIn), idRes(idIn), xRes(xIn),
44 0 : companionRes(companionIn), xqCompRes(0.), mRes(0.), factorRes(1.),
45 0 : colRes(0), acolRes(0) { }
46 :
47 : // Set info on initiator or remnant parton.
48 0 : void iPos( int iPosIn) {iPosRes = iPosIn;}
49 0 : void id( int idIn) {idRes = idIn;}
50 0 : void x( double xIn) {xRes = xIn;}
51 0 : void update( int iPosIn, int idIn, double xIn) {iPosRes = iPosIn;
52 0 : idRes = idIn; xRes = xIn;}
53 0 : void companion( int companionIn) {companionRes = companionIn;}
54 0 : void xqCompanion( double xqCompIn) {xqCompRes = xqCompIn;}
55 0 : void p(Vec4 pIn) {pRes = pIn;}
56 0 : void px(double pxIn) {pRes.px(pxIn);}
57 0 : void py(double pyIn) {pRes.py(pyIn);}
58 0 : void pz(double pzIn) {pRes.pz(pzIn);}
59 0 : void e(double eIn) {pRes.e(eIn);}
60 0 : void m(double mIn) {mRes = mIn;}
61 0 : void col(int colIn) {colRes = colIn;}
62 0 : void acol(int acolIn) {acolRes = acolIn;}
63 : void cols(int colIn = 0,int acolIn = 0)
64 0 : {colRes = colIn; acolRes = acolIn;}
65 0 : void scalePT( double factorIn) {pRes.px(factorIn * pRes.px());
66 0 : pRes.py(factorIn * pRes.py()); factorRes *= factorIn;}
67 0 : void scaleX( double factorIn) {xRes *= factorIn;}
68 :
69 : // Get info on initiator or remnant parton.
70 0 : int iPos() const {return iPosRes;}
71 0 : int id() const {return idRes;}
72 0 : double x() const {return xRes;}
73 0 : int companion() const {return companionRes;}
74 0 : bool isValence() const {return (companionRes == -3);}
75 0 : bool isUnmatched() const {return (companionRes == -2);}
76 0 : bool isCompanion() const {return (companionRes >= 0);}
77 0 : bool isFromBeam() const {return (companionRes > -10);}
78 0 : double xqCompanion() const {return xqCompRes;}
79 0 : Vec4 p() const {return pRes;}
80 0 : double px() const {return pRes.px();}
81 0 : double py() const {return pRes.py();}
82 0 : double pz() const {return pRes.pz();}
83 0 : double e() const {return pRes.e();}
84 0 : double m() const {return mRes;}
85 : double pT() const {return pRes.pT();}
86 0 : double mT2() const {return (mRes >= 0.)
87 0 : ? mRes*mRes + pRes.pT2() : - mRes*mRes + pRes.pT2();}
88 0 : double pPos() const {return pRes.e() + pRes.pz();}
89 0 : double pNeg() const {return pRes.e() - pRes.pz();}
90 0 : int col() const {return colRes;}
91 0 : int acol() const {return acolRes;}
92 0 : double pTfactor() const {return factorRes;}
93 0 : bool hasCol() const {return (idRes == 21 || (idRes > 0 && idRes < 9)
94 0 : || (-idRes > 1000 && -idRes < 10000 && (-idRes/10)%10 == 0));}
95 0 : bool hasAcol() const {return (idRes == 21 || (-idRes > 0 && -idRes < 9)
96 0 : || (idRes > 1000 && idRes < 10000 && (idRes/10)%10 == 0));}
97 :
98 : private:
99 :
100 : // Properties of a resolved parton.
101 : int iPosRes, idRes;
102 : double xRes;
103 : // Companion code and distribution value, if any.
104 : int companionRes;
105 : double xqCompRes;
106 : // Four-momentum and mass; for remnant kinematics construction.
107 : Vec4 pRes;
108 : double mRes, factorRes;
109 : // Colour codes.
110 : int colRes, acolRes;
111 :
112 : };
113 :
114 : //==========================================================================
115 :
116 : // This class holds info on a beam particle in the evolution of
117 : // initial-state radiation and multiparton interactions.
118 :
119 0 : class BeamParticle {
120 :
121 : public:
122 :
123 : // Constructor.
124 0 : BeamParticle() : nInit(0) {resolved.resize(0); Q2ValFracSav = -1.;}
125 :
126 : // Initialize data on a beam particle and save pointers.
127 : void init( int idIn, double pzIn, double eIn, double mIn,
128 : Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn,
129 : Rndm* rndmPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn,
130 : StringFlav* flavSelPtrIn);
131 :
132 : // Initialize only the two pdf pointers.
133 : void initPDFPtr(PDF* pdfInPtr, PDF* pdfHardInPtr) {
134 : pdfBeamPtr = pdfInPtr; pdfHardBeamPtr = pdfHardInPtr; }
135 :
136 : // For mesons like pi0 valence content varies from event to event.
137 : void newValenceContent();
138 :
139 : // Set new pZ and E, but keep the rest the same.
140 0 : void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);}
141 :
142 : // Member functions for output.
143 0 : int id() const {return idBeam;}
144 0 : Vec4 p() const {return pBeam;}
145 : double px() const {return pBeam.px();}
146 : double py() const {return pBeam.py();}
147 : double pz() const {return pBeam.pz();}
148 0 : double e() const {return pBeam.e();}
149 0 : double m() const {return mBeam;}
150 0 : bool isLepton() const {return isLeptonBeam;}
151 0 : bool isUnresolved() const {return isUnresolvedBeam;}
152 : // As hadrons here we only count those we know how to handle remnants for.
153 0 : bool isHadron() const {return isHadronBeam;}
154 : bool isMeson() const {return isMesonBeam;}
155 0 : bool isBaryon() const {return isBaryonBeam;}
156 :
157 : // Maximum x remaining after previous MPI and ISR, plus safety margin.
158 : double xMax(int iSkip = -1);
159 :
160 : // Special hard-process parton distributions (can agree with standard ones).
161 : double xfHard(int idIn, double x, double Q2)
162 0 : {return pdfHardBeamPtr->xf(idIn, x, Q2);}
163 :
164 : // Standard parton distributions.
165 : double xf(int idIn, double x, double Q2)
166 0 : {return pdfBeamPtr->xf(idIn, x, Q2);}
167 :
168 : // Ditto, split into valence and sea parts (where gluon counts as sea).
169 : double xfVal(int idIn, double x, double Q2)
170 0 : {return pdfBeamPtr->xfVal(idIn, x, Q2);}
171 : double xfSea(int idIn, double x, double Q2)
172 0 : {return pdfBeamPtr->xfSea(idIn, x, Q2);}
173 :
174 : // Rescaled parton distributions, as needed for MPI and ISR.
175 : // For ISR also allow split valence/sea, and only return relevant part.
176 : double xfMPI(int idIn, double x, double Q2)
177 0 : {return xfModified(-1, idIn, x, Q2);}
178 : double xfISR(int indexMPI, int idIn, double x, double Q2)
179 0 : {return xfModified( indexMPI, idIn, x, Q2);}
180 :
181 : // Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).
182 : bool insideBounds(double x, double Q2)
183 : {return pdfBeamPtr->insideBounds(x,Q2);}
184 :
185 : // Access the running alpha_s of a PDF set (LHAPDF6 only).
186 : double alphaS(double Q2) {return pdfBeamPtr->alphaS(Q2);}
187 :
188 : // Return quark masses used in the PDF fit (LHAPDF6 only).
189 : double mQuarkPDF(int idIn) {return pdfBeamPtr->mQuarkPDF(idIn);}
190 :
191 : // Decide whether chosen quark is valence, sea or companion.
192 : int pickValSeaComp();
193 :
194 : // Initialize kind of incoming beam particle.
195 : void initBeamKind();
196 :
197 : // Overload index operator to access a resolved parton from the list.
198 0 : ResolvedParton& operator[](int i) {return resolved[i];}
199 : const ResolvedParton& operator[](int i) const {return resolved[i];}
200 :
201 : // Total number of partons extracted from beam, and initiators only.
202 0 : int size() const {return resolved.size();}
203 0 : int sizeInit() const {return nInit;}
204 :
205 : // Clear list of resolved partons.
206 0 : void clear() {resolved.resize(0); nInit = 0;}
207 :
208 : // Add a resolved parton to list.
209 : int append( int iPos, int idIn, double x, int companion = -1)
210 0 : {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) );
211 0 : return resolved.size() - 1;}
212 :
213 : // Print extracted parton list; for debug mainly.
214 : void list(ostream& os = cout) const;
215 :
216 : // How many different flavours, and how many quarks of given flavour.
217 : int nValenceKinds() const {return nValKinds;}
218 0 : int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i)
219 0 : if (idIn == idVal[i]) return nVal[i]; return 0;}
220 :
221 : // Test whether a lepton is to be considered as unresolved.
222 : bool isUnresolvedLepton();
223 :
224 : // Add extra remnant flavours to make valence and sea come out right.
225 : bool remnantFlavours(Event& event, bool isDIS = false);
226 :
227 : // Correlate all initiators and remnants to make a colour singlet.
228 : bool remnantColours(Event& event, vector<int>& colFrom,
229 : vector<int>& colTo);
230 :
231 : // Pick unrescaled x of remnant parton (valence or sea).
232 : double xRemnant(int i);
233 :
234 : // Tell whether a junction has been resolved, and its junction colours.
235 : bool hasJunction() const {return hasJunctionBeam;}
236 : int junctionCol(int i) const {return junCol[i];}
237 : void junctionCol(int i, int col) {junCol[i] = col;}
238 :
239 : // For a diffractive system, decide whether to kick out gluon or quark.
240 : bool pickGluon(double mDiff);
241 :
242 : // Pick a valence quark at random, and provide the remaining flavour.
243 : int pickValence();
244 0 : int pickRemnant() const {return idVal2;}
245 :
246 : // Share lightcone momentum between two remnants in a diffractive system.
247 : // At the same time generate a relative pT for the two.
248 : double zShare( double mDiff, double m1, double m2);
249 0 : double pxShare() const {return pxRel;}
250 0 : double pyShare() const {return pyRel;}
251 :
252 : // Add extra remnant flavours to make valence and sea come out right.
253 : bool remnantFlavoursNew(Event& event);
254 :
255 : // Find the colour setup of the removed partons from the scatterings.
256 : void findColSetup(Event& event);
257 :
258 : // Set initial colours.
259 : void setInitialCol(Event & event);
260 :
261 : // Update colours.
262 : void updateCol(vector<pair<int,int> > colourChanges);
263 :
264 0 : vector<pair <int,int> > getColUpdates() {return colUpdates;}
265 :
266 : private:
267 :
268 : // Constants: could only be changed in the code itself.
269 : static const double XMINUNRESOLVED, POMERONMASS;
270 : static const int NMAX, NRANDOMTRIES;
271 :
272 : // Pointer to various information on the generation.
273 : Info* infoPtr;
274 :
275 : // Pointer to the particle data table.
276 : ParticleData* particleDataPtr;
277 :
278 : // Pointer to the random number generator.
279 : Rndm* rndmPtr;
280 :
281 : // Pointers to PDF sets.
282 : PDF* pdfBeamPtr;
283 : PDF* pdfHardBeamPtr;
284 :
285 : // Pointer to class for flavour generation.
286 : StringFlav* flavSelPtr;
287 :
288 : // Initialization data, normally only set once.
289 : bool allowJunction, beamJunction;
290 : int maxValQuark, companionPower;
291 : double valencePowerMeson, valencePowerUinP, valencePowerDinP,
292 : valenceDiqEnhance, pickQuarkNorm, pickQuarkPower,
293 : diffPrimKTwidth, diffLargeMassSuppress, beamSat, gluonPower,
294 : xGluonCutoff;
295 :
296 : // Basic properties of a beam particle.
297 : int idBeam, idBeamAbs;
298 : Vec4 pBeam;
299 : double mBeam;
300 : // Beam kind. Valence flavour content for hadrons.
301 : bool isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam,
302 : isBaryonBeam;
303 : int nValKinds, idVal[3], nVal[3];
304 :
305 : // Current parton density, by valence, sea and companion.
306 : int idSave, iSkipSave, nValLeft[3];
307 : double xqgTot, xqVal, xqgSea, xqCompSum;
308 :
309 : // The list of resolved partons.
310 : vector<ResolvedParton> resolved;
311 :
312 : // Status after all initiators have been accounted for. Junction content.
313 : int nInit;
314 : bool hasJunctionBeam;
315 : int junCol[3];
316 :
317 : // Variables for new colour reconnection;
318 : pair <int,int> colSetup;
319 : vector<int> acols, cols;
320 : vector<bool> usedCol,usedAcol;
321 : vector< pair<int,int> > colUpdates;
322 : int nJuncs, nAjuncs, nDiffJuncs;
323 : bool allowBeamJunctions;
324 :
325 : // Routine to calculate pdf's given previous interactions.
326 : double xfModified( int iSkip, int idIn, double x, double Q2);
327 :
328 : // Fraction of hadron momentum sitting in a valence quark distribution.
329 : double xValFrac(int j, double Q2);
330 : double Q2ValFracSav, uValInt, dValInt;
331 :
332 : // Fraction of hadron momentum sitting in a companion quark distribution.
333 : double xCompFrac(double xs);
334 :
335 : // Value of companion quark PDF, also given the sea quark x.
336 : double xCompDist(double xc, double xs);
337 :
338 : // Valence quark subdivision for diffractive systems.
339 : int idVal1, idVal2, idVal3;
340 : double zRel, pxRel, pyRel;
341 :
342 : // Update a single (anti-) colour of the event.
343 : void updateSingleCol(int oldCol, int newCol);
344 :
345 : // Find a single (anti-) colour in the beam,
346 : // if a beam remnant is set the new colour.
347 : int findSingleCol(Event& event, bool isAcol, bool useHardScatters);
348 :
349 : };
350 :
351 : //==========================================================================
352 :
353 : } // end namespace Pythia8
354 :
355 : #endif // Pythia8_BeamParticle_H
|
