StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BeamParticle.h
1 // BeamParticle.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 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  // Initialize only the two pdf pointers.
129  void initPDFPtr(PDF* pdfInPtr, PDF* pdfHardInPtr) {
130  pdfBeamPtr = pdfInPtr; pdfHardBeamPtr = pdfHardInPtr; }
131 
132  // For mesons like pi0 valence content varies from event to event.
133  void newValenceContent();
134 
135  // Set new pZ and E, but keep the rest the same.
136  void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);}
137 
138  // Member functions for output.
139  int id() const {return idBeam;}
140  Vec4 p() const {return pBeam;}
141  double px() const {return pBeam.px();}
142  double py() const {return pBeam.py();}
143  double pz() const {return pBeam.pz();}
144  double e() const {return pBeam.e();}
145  double m() const {return mBeam;}
146  bool isLepton() const {return isLeptonBeam;}
147  bool isUnresolved() const {return isUnresolvedBeam;}
148  // As hadrons here we only count those we know how to handle remnants for.
149  bool isHadron() const {return isHadronBeam;}
150  bool isMeson() const {return isMesonBeam;}
151  bool isBaryon() const {return isBaryonBeam;}
152 
153  // Maximum x remaining after previous MPI and ISR, plus safety margin.
154  double xMax(int iSkip = -1);
155 
156  // Special hard-process parton distributions (can agree with standard ones).
157  double xfHard(int idIn, double x, double Q2)
158  {return pdfHardBeamPtr->xf(idIn, x, Q2);}
159 
160  // Standard parton distributions.
161  double xf(int idIn, double x, double Q2)
162  {return pdfBeamPtr->xf(idIn, x, Q2);}
163 
164  // Ditto, split into valence and sea parts (where gluon counts as sea).
165  double xfVal(int idIn, double x, double Q2)
166  {return pdfBeamPtr->xfVal(idIn, x, Q2);}
167  double xfSea(int idIn, double x, double Q2)
168  {return pdfBeamPtr->xfSea(idIn, x, Q2);}
169 
170  // Rescaled parton distributions, as needed for MPI and ISR.
171  // For ISR also allow split valence/sea, and only return relevant part.
172  double xfMPI(int idIn, double x, double Q2)
173  {return xfModified(-1, idIn, x, Q2);}
174  double xfISR(int indexMPI, int idIn, double x, double Q2)
175  {return xfModified( indexMPI, idIn, x, Q2);}
176 
177  // Decide whether chosen quark is valence, sea or companion.
178  int pickValSeaComp();
179 
180  // Initialize kind of incoming beam particle.
181  void initBeamKind();
182 
183  // Overload index operator to access a resolved parton from the list.
184  ResolvedParton& operator[](int i) {return resolved[i];}
185  const ResolvedParton& operator[](int i) const {return resolved[i];}
186 
187  // Total number of partons extracted from beam, and initiators only.
188  int size() const {return resolved.size();}
189  int sizeInit() const {return nInit;}
190 
191  // Clear list of resolved partons.
192  void clear() {resolved.resize(0); nInit = 0;}
193 
194  // Add a resolved parton to list.
195  int append( int iPos, int idIn, double x, int companion = -1)
196  {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) );
197  return resolved.size() - 1;}
198 
199  // Print extracted parton list; for debug mainly.
200  void list(ostream& os = cout) const;
201 
202  // How many different flavours, and how many quarks of given flavour.
203  int nValenceKinds() const {return nValKinds;}
204  int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i)
205  if (idIn == idVal[i]) return nVal[i]; return 0;}
206 
207  // Test whether a lepton is to be considered as unresolved.
208  bool isUnresolvedLepton();
209 
210  // Add extra remnant flavours to make valence and sea come out right.
211  bool remnantFlavours(Event& event);
212 
213  // Correlate all initiators and remnants to make a colour singlet.
214  bool remnantColours(Event& event, vector<int>& colFrom,
215  vector<int>& colTo);
216 
217  // Pick unrescaled x of remnant parton (valence or sea).
218  double xRemnant(int i);
219 
220  // Tell whether a junction has been resolved, and its junction colours.
221  bool hasJunction() const {return hasJunctionBeam;}
222  int junctionCol(int i) const {return junCol[i];}
223  void junctionCol(int i, int col) {junCol[i] = col;}
224 
225  // For a diffractive system, decide whether to kick out gluon or quark.
226  bool pickGluon(double mDiff);
227 
228  // Pick a valence quark at random, and provide the remaining flavour.
229  int pickValence();
230  int pickRemnant() const {return idVal2;}
231 
232  // Share lightcone momentum between two remnants in a diffractive system.
233  // At the same time generate a relative pT for the two.
234  double zShare( double mDiff, double m1, double m2);
235  double pxShare() const {return pxRel;}
236  double pyShare() const {return pyRel;}
237 
238 private:
239 
240  // Constants: could only be changed in the code itself.
241  static const double XMINUNRESOLVED;
242 
243  // Pointer to various information on the generation.
244  Info* infoPtr;
245 
246  // Pointer to the particle data table.
247  ParticleData* particleDataPtr;
248 
249  // Pointer to the random number generator.
250  Rndm* rndmPtr;
251 
252  // Pointers to PDF sets.
253  PDF* pdfBeamPtr;
254  PDF* pdfHardBeamPtr;
255 
256  // Pointer to class for flavour generation.
257  StringFlav* flavSelPtr;
258 
259  // Initialization data, normally only set once.
260  bool allowJunction;
261  int maxValQuark, companionPower;
262  double valencePowerMeson, valencePowerUinP, valencePowerDinP,
263  valenceDiqEnhance, pickQuarkNorm, pickQuarkPower,
264  diffPrimKTwidth, diffLargeMassSuppress;
265 
266  // Basic properties of a beam particle.
267  int idBeam, idBeamAbs;
268  Vec4 pBeam;
269  double mBeam;
270  // Beam kind. Valence flavour content for hadrons.
271  bool isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam,
272  isBaryonBeam;
273  int nValKinds, idVal[3], nVal[3];
274 
275  // Current parton density, by valence, sea and companion.
276  int idSave, iSkipSave, nValLeft[3];
277  double xqgTot, xqVal, xqgSea, xqCompSum;
278 
279  // The list of resolved partons.
280  vector<ResolvedParton> resolved;
281 
282  // Status after all initiators have been accounted for. Junction content.
283  int nInit;
284  bool hasJunctionBeam;
285  int junCol[3];
286 
287  // Routine to calculate pdf's given previous interactions.
288  double xfModified( int iSkip, int idIn, double x, double Q2);
289 
290  // Fraction of hadron momentum sitting in a valence quark distribution.
291  double xValFrac(int j, double Q2);
292  double Q2ValFracSav, uValInt, dValInt;
293 
294  // Fraction of hadron momentum sitting in a companion quark distribution.
295  double xCompFrac(double xs);
296 
297  // Value of companion quark PDF, also given the sea quark x.
298  double xCompDist(double xc, double xs);
299 
300  // Valence quark subdivision for diffractive systems.
301  int idVal1, idVal2, idVal3;
302  double zRel, pxRel, pyRel;
303 
304 };
305 
306 //==========================================================================
307 
308 } // end namespace Pythia8
309 
310 #endif // Pythia8_BeamParticle_H
Definition: AgUStep.h:26