]>
Commit | Line | Data |
---|---|---|
63ba5337 | 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 |