1 | // HadronScatter.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 | #ifndef Pythia8_HadronScatter_H |
---|
7 | #define Pythia8_HadronScatter_H |
---|
8 | |
---|
9 | #include "Event.h" |
---|
10 | #include "Info.h" |
---|
11 | #include "ParticleData.h" |
---|
12 | #include "PythiaComplex.h" |
---|
13 | |
---|
14 | #include <algorithm> |
---|
15 | #include <map> |
---|
16 | #include <set> |
---|
17 | using std::map; |
---|
18 | using std::set; |
---|
19 | using std::sort; |
---|
20 | |
---|
21 | namespace Pythia8 { |
---|
22 | |
---|
23 | class SigmaPartialWave { |
---|
24 | public: |
---|
25 | // Initialisation |
---|
26 | bool init(int, string, string, Info *, ParticleData *, Rndm *); |
---|
27 | |
---|
28 | // Read data file |
---|
29 | bool readFile(string, string); |
---|
30 | |
---|
31 | // Set the subprocess/incoming particles |
---|
32 | bool setSubprocess(int); |
---|
33 | bool setSubprocess(int, int); |
---|
34 | |
---|
35 | // Return sigma total/elastic, dSigma/dCos(theta) |
---|
36 | double sigmaEl(double Wcm) { return sigma(0, Wcm); } |
---|
37 | double sigmaTot(double Wcm) { return sigma(1, Wcm); } |
---|
38 | double dSigma(double Wcm, double cTheta) { return sigma(2, Wcm, cTheta); } |
---|
39 | |
---|
40 | // Return a cos(theta) value |
---|
41 | double pickCosTheta(double); |
---|
42 | |
---|
43 | // Return maximum sigma elastic |
---|
44 | double getSigmaElMax() { return sigElMax; } |
---|
45 | |
---|
46 | private: |
---|
47 | // Pointers |
---|
48 | Info *infoPtr; |
---|
49 | ParticleData *particleDataPtr; |
---|
50 | Rndm *rndmPtr; |
---|
51 | |
---|
52 | // Constants |
---|
53 | static const int LSHIFT, ISHIFT, SUBBIN, ITER; |
---|
54 | static const double CONVERT2MB, WCMBIN, CTBIN, MASSSAFETY, GRIDSAFETY; |
---|
55 | static const complex BINEND; |
---|
56 | |
---|
57 | // Settings |
---|
58 | int process, subprocess, subprocessMax, norm; |
---|
59 | |
---|
60 | // Current incoming and maximum L/I values |
---|
61 | int idA, idB, Lmax, Imax; |
---|
62 | |
---|
63 | // Masses of incoming, last bin, maximum sigma elastic |
---|
64 | double mA, mB, binMax, sigElMax; |
---|
65 | |
---|
66 | // Map subprocess to incoming and vice versa: |
---|
67 | // sp2in[subprocess] = idA, idB |
---|
68 | // in2sp[idA, idB] = subprocess |
---|
69 | map < int, pair < int, int > > sp2in; |
---|
70 | map < pair < int, int >, int > in2sp; |
---|
71 | |
---|
72 | // Isospin coefficients isoCoeff[subprocess][2I] |
---|
73 | map < int, map < int, double > > isoCoeff; |
---|
74 | |
---|
75 | // Storage for partial wave data: |
---|
76 | // pwData[L * LSHIFT + 2I * ISHIFT + 2J][eNow] = T |
---|
77 | map < int, map < double, complex > > pwData; |
---|
78 | |
---|
79 | // Values of Pl and Pl' as computed by legendreP |
---|
80 | vector < double > PlVec, PlpVec; |
---|
81 | |
---|
82 | // Integration grid [subprocess][WcmBin][cosThetaBin] |
---|
83 | vector < vector < vector < double > > > gridMax; |
---|
84 | vector < vector < double > > gridNorm; |
---|
85 | |
---|
86 | // Setup subprocesses (including isospin coefficients) |
---|
87 | void setupSubprocesses(); |
---|
88 | |
---|
89 | // Setup grids for integration |
---|
90 | void setupGrid(); |
---|
91 | |
---|
92 | // Routine for calculating sigma total/elastic and dSigma/dCos(theta) |
---|
93 | double sigma(int, double, double = 0.); |
---|
94 | |
---|
95 | // Generate Legendre polynomials (and optionally derivatives) |
---|
96 | void legendreP(double, bool = false); |
---|
97 | |
---|
98 | }; |
---|
99 | |
---|
100 | |
---|
101 | //========================================================================== |
---|
102 | |
---|
103 | // HadronScatterPair class |
---|
104 | // Simple class to hold details of a pair of hadrons which will scatter. |
---|
105 | // Stores indices in event record and the measure used for ordering |
---|
106 | |
---|
107 | // Store a pair of indices |
---|
108 | typedef pair < int, int > HSIndex; |
---|
109 | |
---|
110 | class HadronScatterPair { |
---|
111 | public: |
---|
112 | // Constructor |
---|
113 | HadronScatterPair(const HSIndex &i1in, int yt1in, int pt1in, |
---|
114 | const HSIndex &i2in, int yt2in, int pt2in, |
---|
115 | double measureIn) : |
---|
116 | i1(i1in), yt1(yt1in), pt1(pt1in), |
---|
117 | i2(i2in), yt2(yt2in), pt2(pt2in), |
---|
118 | measure(measureIn) {} |
---|
119 | |
---|
120 | // Operator for sorting according to ordering measure |
---|
121 | bool operator<(const HadronScatterPair& in) const { |
---|
122 | return this->measure < in.measure; |
---|
123 | } |
---|
124 | |
---|
125 | // Indices into event record of hadrons to scatter |
---|
126 | HSIndex i1; |
---|
127 | int yt1, pt1; |
---|
128 | HSIndex i2; |
---|
129 | int yt2, pt2; |
---|
130 | // Ordering measure |
---|
131 | double measure; |
---|
132 | }; |
---|
133 | |
---|
134 | |
---|
135 | //========================================================================== |
---|
136 | |
---|
137 | // HadronScatter class |
---|
138 | |
---|
139 | class HadronScatter { |
---|
140 | |
---|
141 | public: |
---|
142 | |
---|
143 | // Constructor. |
---|
144 | HadronScatter() {} |
---|
145 | |
---|
146 | // Initialisation |
---|
147 | bool init(Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn, |
---|
148 | ParticleData *particleDataPtr); |
---|
149 | |
---|
150 | // Perform all hadron scatterings |
---|
151 | void scatter(Event&); |
---|
152 | |
---|
153 | private: |
---|
154 | |
---|
155 | // Pointer to various information on the generation. |
---|
156 | Info* infoPtr; |
---|
157 | Rndm* rndmPtr; |
---|
158 | |
---|
159 | // Settings |
---|
160 | bool doHadronScatter, afterDecay, allowDecayProd, |
---|
161 | scatterRepeat, doTile; |
---|
162 | int hadronSelect, scatterProb; |
---|
163 | double Npar, kPar, pPar, jPar, rMax, rMax2; |
---|
164 | double pTsigma, pTsigma2, pT0MPI; |
---|
165 | |
---|
166 | // Tiling |
---|
167 | int ytMax, ptMax; |
---|
168 | double yMin, yMax, ytSize, ptSize; |
---|
169 | vector < vector < set < HSIndex > > > tile; |
---|
170 | |
---|
171 | // Keep track of scattered pairs |
---|
172 | set < HSIndex > scattered; |
---|
173 | |
---|
174 | // Partial wave amplitudes |
---|
175 | SigmaPartialWave sigmaPW[3]; |
---|
176 | |
---|
177 | // Maximum sigma elastic |
---|
178 | double sigElMax; |
---|
179 | |
---|
180 | // Decide if a hadron can scatter |
---|
181 | bool canScatter(Event &, int); |
---|
182 | |
---|
183 | // Probability for a pair of hadrons to scatter |
---|
184 | bool doesScatter(Event &, const HSIndex &, const HSIndex &); |
---|
185 | |
---|
186 | // Calculate measure for ordering of scatterings |
---|
187 | double measure(Event &, int, int); |
---|
188 | |
---|
189 | // Perform a single hadron scattering |
---|
190 | bool hadronScatter(Event &, HadronScatterPair &); |
---|
191 | |
---|
192 | // Tiling functions |
---|
193 | bool tileIntProb(vector < HadronScatterPair > &, Event &, |
---|
194 | const HSIndex &, int, int, bool); |
---|
195 | int yTile(Event& event, int idx) { |
---|
196 | return (doTile) ? int((event[idx].y() - yMin) / ytSize) : 0; |
---|
197 | } |
---|
198 | int pTile(Event& event, int idx) { |
---|
199 | return (doTile) ? int((event[idx].phi() + M_PI) / ptSize) : 0; |
---|
200 | } |
---|
201 | |
---|
202 | // Debug |
---|
203 | void debugOutput(); |
---|
204 | }; |
---|
205 | |
---|
206 | //========================================================================== |
---|
207 | |
---|
208 | |
---|
209 | } // end namespace Pythia8 |
---|
210 | |
---|
211 | #endif // Pythia8_HadronScatter_H |
---|
212 | |
---|