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