- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / include / BeamParticle.h
1 // BeamParticle.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 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 multiple 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 multiple interactions.
114
115 class BeamParticle {
116
117 public:
118
119   // Constructor.
120   BeamParticle() : nInit(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 MI 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 MI and ISR.
167   // For ISR also allow split valence/sea, and only return relevant part.
168   double xfMI(int idIn, double x, double Q2) 
169     {return xfModified(-1, idIn, x, Q2);}
170   double xfISR(int indexMI, int idIn, double x, double Q2) 
171     {return xfModified( indexMI, 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   isLeptonBeam, isUnresolvedBeam, 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