Line data Source code
1 : // HadronScatter.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 : #include "Pythia8/HadronScatter.h"
7 :
8 : namespace Pythia8 {
9 :
10 : //==========================================================================
11 :
12 : // The SigmaPartialWave class
13 : // Reads in tables of partial wave data to provide dSigma/dCos(theta)
14 : // The generic classes of process are:
15 : // process = 0 (pi-pi), 1 (pi-K), 2 (pi-N)
16 : // Subprocesses are defined (along with isospin coefficients) in:
17 : // setupSubprocesses();
18 : // Individual subprocesses are selected using:
19 : // setSubprocess(subprocess); or setSubprocess(PDG1, PDG2);
20 : // Internally, there are two std::map's, to convert between:
21 : // subprocess <==> PDG1, PDG2
22 : //
23 : // Data are read in from files:
24 : // Lines starting with a '#' are comments
25 : // Lines starting with 'set' provide options:
26 : // set eType [Wcm | Tlab] - energy bins in Wcm or Tlab
27 : // set eUnit [MeV | GeV] - energy unit
28 : // set input [eta,delta | Sn,delta | Tr,Ti | mod,phi ]
29 : // - format of columns in partial waves
30 : // set dUnit [deg | rad] - unit of phase shifts
31 : // set norm [0 | 1] - normalisation
32 : // Column headers give L,2I[,2J] (2J for e.g. piN)
33 : // Input types: Sn,delta -> Sn = 1 - eta^2
34 : // mod,phi -> amplitude T_L = |T_L| exp(i phi_L)
35 : // Normalisation: 0 -> dSigma/dOmega = 1 / k^2 |T_L|^2
36 : // 1 -> dSigma/dOmega = 16 / s |T_L|^2
37 : //
38 : // Internally data is stored as (J = 0 for spinless):
39 : // pwData[L * LSHIFT + 2I * ISHIFT + J][energy_bin_centre] = T
40 : // where the energy is Wcm in GeV.
41 : //
42 : // This is stored using std::map's, to take into account that not all
43 : // L,I,J states are always present (e.g. negligable contributions or
44 : // conservation rules) and that bin sizes are not fixed.
45 : //
46 : // Re[T] and Im[T] are interpolated between bins and extrapolated down to
47 : // threshold from the first two bins. Above energy_bin_centre of the final
48 : // bin, no extrapolation is done and the final bin value is always used.
49 : //
50 : // A simple scheme to provide correct distributions for cos(theta) at a
51 : // given CM energy is included. Efficiency is not too bad, but can likely
52 : // be greatly improved.
53 : //
54 : // For each subprocess, a grid in bins of Wcm and cos(theta) is setup with:
55 : // setupGrid();
56 : // The size of the grid is set by the constants:
57 : // const double SigmaPartialWave::WCMBIN;
58 : // const double SigmaPartialWave::CTBIN;
59 : // For each bin of (Wcm, ct), the maximum sigma elastic is found by
60 : // splitting this bin into subbins multiple times, controlled by:
61 : // const int SigmaPartialWave::SUBBIN;
62 : // const int SigmaPartialWave::ITER
63 : // With the final maxium sigma elastic given by this value multipled by
64 : // a safety factor:
65 : // const double SigmaPartialWave::GRIDSAFETY
66 : //
67 : // To pick values of cos(theta) for a given CM energy, a:
68 : // pickCosTheta(Wcm);
69 : // function is provided. The above grid is used as an overestimate, to
70 : // pick properly distributed values of cos(theta).
71 :
72 : //--------------------------------------------------------------------------
73 :
74 : // Constants
75 :
76 : // pwData[L * LSHIFT + 2I * ISHIFT + J]
77 : const int SigmaPartialWave::LSHIFT = 1000000;
78 : const int SigmaPartialWave::ISHIFT = 1000;
79 :
80 : // Convert GeV^-2 to mb
81 : const double SigmaPartialWave::CONVERT2MB = 0.389380;
82 :
83 : // Size of bin in Wcm and cos(theta)
84 : const double SigmaPartialWave::WCMBIN = 0.005;
85 : const double SigmaPartialWave::CTBIN = 0.2;
86 : // Number of subbins and iterations
87 : const int SigmaPartialWave::SUBBIN = 2;
88 : const int SigmaPartialWave::ITER = 2;
89 : // Safety value to add on to grid maxima
90 : const double SigmaPartialWave::MASSSAFETY = 0.001;
91 : const double SigmaPartialWave::GRIDSAFETY = 0.05;
92 :
93 :
94 : //--------------------------------------------------------------------------
95 :
96 : // Perform initialization and store pointers.
97 :
98 : bool SigmaPartialWave::init(int processIn, string xmlPath, string filename,
99 : Info *infoPtrIn, ParticleData *particleDataPtrIn,
100 : Rndm *rndmPtrIn) {
101 : // Store incoming pointers
102 0 : infoPtr = infoPtrIn;
103 0 : particleDataPtr = particleDataPtrIn;
104 0 : rndmPtr = rndmPtrIn;
105 :
106 : // Check incoming process is okay
107 0 : if (processIn < 0 || processIn > 2) {
108 0 : infoPtr->errorMsg("Error in SigmaPartialWave::init: "
109 : "unknown process");
110 0 : return false;
111 : }
112 0 : process = processIn;
113 :
114 : // Setup subprocesses and isospin coefficients
115 0 : setupSubprocesses();
116 0 : setSubprocess(0);
117 :
118 : // Read in partial-wave data
119 0 : if (!readFile(xmlPath, filename)) return false;
120 :
121 : // Setup vector for Legendre polynomials
122 0 : PlVec.resize(Lmax);
123 0 : if (Lmax > 0) PlVec[0] = 1.;
124 : // And derivatives if needed
125 0 : if (process == 2) {
126 0 : PlpVec.resize(Lmax);
127 0 : if (Lmax > 0) PlpVec[0] = 0.;
128 0 : if (Lmax > 1) PlpVec[1] = 1.;
129 : }
130 :
131 : // Setup grid for integration
132 0 : setupGrid();
133 :
134 0 : return true;
135 0 : }
136 :
137 :
138 : //--------------------------------------------------------------------------
139 :
140 : // Read input data file
141 :
142 : bool SigmaPartialWave::readFile(string xmlPath, string filename) {
143 : // Create full path and open file
144 0 : string fullPath = xmlPath + filename;
145 0 : ifstream ifs(fullPath.c_str());
146 0 : if (!ifs.good()) {
147 0 : infoPtr->errorMsg("Error in SigmaPartialWave::init: "
148 : "could not read data file");
149 0 : return false;
150 : }
151 :
152 : // Default unit settings
153 : int eType = 0; // 0 = Wcm, 1 = Tlab
154 : int eUnit = 0; // 0 = GeV, 1 = MeV
155 : int input = 0; // 0 = eta, delta; 1 = Sn, delta (Sn = 1 - eta^2);
156 : // 2 = Treal, Tim, 3 = mod, phi
157 : int dUnit = 0; // 0 = deg, 1 = rad
158 0 : norm = 0; // 0 = standard, 1 = sqrt(s) / sqrt(s - 4Mpi^2)
159 :
160 : // Once we have a header line, each column corresponds to
161 : // values of L, I and J.
162 0 : Lmax = Imax = 0;
163 0 : binMax = 0.;
164 0 : vector < int > Lvec, Ivec, Jvec;
165 :
166 : // Parse the file
167 0 : string line;
168 0 : while (ifs.good()) {
169 : // Get line, convert to lowercase and strip leading whitespace
170 0 : getline(ifs, line);
171 0 : for (unsigned int i = 0; i < line.length(); i++)
172 0 : line[i] = tolower(line[i]);
173 0 : string::size_type startPos = line.find_first_not_of(" ");
174 0 : if (startPos != string::npos) line = line.substr(startPos);
175 : // Skip blank lines and lines that start with '#'
176 0 : if (line.length() == 0 || line[0] == '#') continue;
177 :
178 : // Tokenise line on whitespace (spaces or tabs)
179 0 : string lineT = line;
180 0 : vector < string > token;
181 0 : while (true) {
182 0 : startPos = lineT.find_first_of(" ");
183 0 : token.push_back(lineT.substr(0, startPos));
184 0 : if (startPos == string::npos) break;
185 0 : startPos = lineT.find_first_not_of(" ", startPos + 1);
186 0 : if (startPos == string::npos) break;
187 0 : lineT = lineT.substr(startPos);
188 : }
189 :
190 : // Settings
191 0 : if (token[0] == "set") {
192 : bool badSetting = false;
193 :
194 : // eType
195 0 : if (token[1] == "etype") {
196 0 : if (token[2] == "wcm") eType = 0;
197 0 : else if (token[2] == "tlab") eType = 1;
198 : else badSetting = true;
199 :
200 : // eUnit
201 0 : } else if (token[1] == "eunit") {
202 0 : if (token[2] == "gev") eUnit = 0;
203 0 : else if (token[2] == "mev") eUnit = 1;
204 : else badSetting = true;
205 :
206 : // input
207 0 : } else if (token[1] == "input") {
208 0 : if (token[2] == "eta,delta") input = 0;
209 0 : else if (token[2] == "sn,delta") input = 1;
210 0 : else if (token[2] == "tr,ti") input = 2;
211 0 : else if (token[2] == "mod,phi") input = 3;
212 : else badSetting = true;
213 :
214 : // dUnit
215 0 : } else if (token[1] == "dunit") {
216 0 : if (token[2] == "deg") dUnit = 0;
217 0 : else if (token[2] == "rad") dUnit = 1;
218 : else badSetting = true;
219 :
220 : // norm
221 0 : } else if (token[1] == "norm") {
222 0 : if (token[2] == "0") norm = 0;
223 0 : else if (token[2] == "1") norm = 1;
224 : else badSetting = true;
225 : }
226 :
227 : // Bad setting
228 0 : if (badSetting) {
229 0 : infoPtr->errorMsg("Error in SigmaPartialWave::init: "
230 : "bad setting line");
231 0 : return false;
232 : }
233 0 : continue;
234 : }
235 :
236 : // Header line
237 0 : if (line.substr(0, 1).find_first_of("0123456789.") != 0) {
238 : // Clear current stored L,2I,2J values
239 0 : Lvec.clear(); Ivec.clear(); Jvec.clear();
240 :
241 : // Parse header
242 : bool badHeader = false;
243 0 : for (unsigned int i = 1; i < token.size(); i++) {
244 : // Extract L
245 0 : startPos = token[i].find_first_of(",");
246 0 : if (startPos == string::npos) { badHeader = true; break; }
247 0 : string Lstr = token[i].substr(0, startPos);
248 0 : token[i] = token[i].substr(startPos + 1);
249 : // Extract 2I
250 0 : string Istr;
251 0 : startPos = token[i].find_first_of(", ");
252 0 : if (startPos == string::npos) {
253 0 : Istr = token[i];
254 0 : token[i] = "";
255 : } else {
256 0 : Istr = token[i].substr(0, startPos);
257 0 : token[i] = token[i].substr(startPos + 1);
258 : }
259 : // Extract 2J
260 0 : string Jstr("0");
261 0 : if (token[i].length() != 0) Jstr = token[i];
262 :
263 : // Convert to integers and store
264 0 : int L, I, J;
265 0 : stringstream Lss(Lstr); Lss >> L;
266 0 : stringstream Iss(Istr); Iss >> I;
267 0 : stringstream Jss(Jstr); Jss >> J;
268 0 : if (Lss.fail() || Iss.fail() || Jss.fail())
269 0 : { badHeader = true; break; }
270 0 : Lvec.push_back(L);
271 0 : Ivec.push_back(I);
272 0 : Jvec.push_back(J);
273 0 : Lmax = max(Lmax, L);
274 0 : Imax = max(Imax, I);
275 0 : }
276 0 : if (badHeader) {
277 0 : infoPtr->errorMsg("Error in SigmaPartialWave::init: "
278 : "malformed header line");
279 0 : return false;
280 : }
281 :
282 : // Data line
283 0 : } else {
284 : bool badData = false;
285 :
286 : // Check there are the correct number of columns
287 0 : if (token.size() != 2 * Lvec.size() + 1) badData = true;
288 :
289 : // Extract energy
290 0 : double eNow = 0.;
291 0 : if (!badData) {
292 0 : stringstream eSS(token[0]);
293 0 : eSS >> eNow;
294 0 : if (eSS.fail()) badData = true;
295 : // Convert to GeV if needed
296 0 : if (eUnit == 1) eNow *= 1e-3;
297 : // Convert to Wcm if needed
298 0 : if (eType == 1) eNow = sqrt(2. * mB * eNow + pow2(mA + mB));
299 0 : binMax = max(binMax, eNow);
300 0 : }
301 :
302 : // Extract eta/phase shifts
303 0 : if (!badData) {
304 0 : for (unsigned int i = 1; i < token.size(); i += 2) {
305 : // L,2I,2J
306 0 : int LIJidx = (i - 1) / 2;
307 0 : int L = Lvec[LIJidx];
308 0 : int I = Ivec[LIJidx];
309 0 : int J = Jvec[LIJidx];
310 :
311 0 : double i1, i2;
312 0 : stringstream i1SS(token[i]); i1SS >> i1;
313 0 : stringstream i2SS(token[i + 1]); i2SS >> i2;
314 0 : if (i1SS.fail() || i2SS.fail()) { badData = true; break; }
315 :
316 : // Sn to eta
317 0 : if (input == 1) i1 = sqrt(1. - i1);
318 : // Degrees to radians
319 0 : if ((input == 0 || input == 1 || input == 3) &&
320 0 : dUnit == 0) i2 *= M_PI / 180.;
321 :
322 : // Convert to Treal and Timg
323 0 : complex T(0., 0.);
324 0 : if (input == 0 || input == 1) {
325 0 : T = (i1 * exp(2. * complex(0., 1.) * i2) - 1.) /
326 0 : 2. / complex(0., 1.);
327 0 : } else if (input == 2) {
328 0 : T = complex(i1, i2);
329 0 : } else if (input == 3) {
330 0 : T = i1 * exp(complex(0., 1.) * i2);
331 0 : }
332 :
333 : // Store
334 0 : pwData[L * LSHIFT + I * ISHIFT + J][eNow] = T;
335 0 : }
336 0 : }
337 0 : if (badData) {
338 0 : infoPtr->errorMsg("Error in SigmaPartialWave::init: "
339 : "malformed data line");
340 0 : return false;
341 : }
342 0 : }
343 0 : }
344 :
345 : // Make sure it was EOF that caused us to end
346 0 : if (!ifs.eof()) { ifs.close(); return false; }
347 :
348 : // Maximum values of L and I
349 0 : Lmax++; Imax++;
350 :
351 0 : return true;
352 0 : }
353 :
354 :
355 : //--------------------------------------------------------------------------
356 :
357 : // Setup isospin coefficients and subprocess mapping
358 :
359 : void SigmaPartialWave::setupSubprocesses() {
360 :
361 : // Setup isospin coefficients
362 0 : switch (process) {
363 : // pi-pi
364 : case 0:
365 : // Map subprocess to incoming
366 0 : subprocessMax = 6;
367 0 : sp2in[0] = pair < int, int > ( 211, 211);
368 0 : sp2in[1] = pair < int, int > ( 211, -211);
369 0 : sp2in[2] = pair < int, int > ( 211, 111);
370 0 : sp2in[3] = pair < int, int > ( 111, 111);
371 0 : sp2in[4] = pair < int, int > (-211, 111);
372 0 : sp2in[5] = pair < int, int > (-211, -211);
373 : // Incoming to subprocess
374 0 : for (int i = 0; i < subprocessMax; i++) in2sp[sp2in[i]] = i;
375 : // Isospin coefficients
376 0 : isoCoeff[0][0] = 0.; isoCoeff[0][2] = 0.; isoCoeff[0][4] = 1.;
377 0 : isoCoeff[1][0] = 1./3.; isoCoeff[1][2] = 1./2.; isoCoeff[1][4] = 1./6.;
378 0 : isoCoeff[2][0] = 0.; isoCoeff[2][2] = 1./2.; isoCoeff[2][4] = 1./2.;
379 0 : isoCoeff[3][0] = 1./3.; isoCoeff[3][2] = 0.; isoCoeff[3][4] = 2./3.;
380 0 : isoCoeff[4][0] = 0.; isoCoeff[4][2] = 1./2.; isoCoeff[4][4] = 1./2.;
381 0 : isoCoeff[5][0] = 0.; isoCoeff[5][2] = 0.; isoCoeff[5][4] = 1.;
382 :
383 0 : break;
384 :
385 : // pi-K and pi-N
386 : case 1: case 2:
387 0 : int id1, id2;
388 0 : if (process == 1) { id1 = 321; id2 = 311; }
389 0 : else { id1 = 2212; id2 = 2112; }
390 :
391 : // Map subprocess to incoming
392 0 : subprocessMax = 12;
393 0 : sp2in[0] = pair < int, int > ( 211, id1);
394 0 : sp2in[1] = pair < int, int > ( 211, id2);
395 0 : sp2in[2] = pair < int, int > ( 111, id1);
396 0 : sp2in[3] = pair < int, int > ( 111, id2);
397 0 : sp2in[4] = pair < int, int > (-211, id1);
398 0 : sp2in[5] = pair < int, int > (-211, id2);
399 : // Isospin coefficients
400 0 : isoCoeff[0][1] = 0.; isoCoeff[0][3] = 1.;
401 0 : isoCoeff[1][1] = 2. / 3.; isoCoeff[1][3] = 1. / 3.;
402 0 : isoCoeff[2][1] = 1. / 3.; isoCoeff[2][3] = 2. / 3.;
403 0 : isoCoeff[3][1] = 1. / 3.; isoCoeff[3][3] = 2. / 3.;
404 0 : isoCoeff[4][1] = 2. / 3.; isoCoeff[4][3] = 1. / 3.;
405 0 : isoCoeff[5][1] = 0.; isoCoeff[5][3] = 1.;
406 : // Antiparticles
407 0 : for (int i = 0; i < 6; i++) {
408 0 : id1 = ((sp2in[i].first == 111) ? +1 : -1) * sp2in[i].first;
409 0 : sp2in[i + 6] = pair < int, int > (id1, -sp2in[i].second);
410 0 : isoCoeff[i + 6] = isoCoeff[i];
411 : }
412 : // Map incoming to subprocess
413 0 : for (int i = 0; i < subprocessMax; i++) in2sp[sp2in[i]] = i;
414 :
415 : break;
416 0 : }
417 :
418 0 : return;
419 : }
420 :
421 :
422 : //--------------------------------------------------------------------------
423 :
424 : // Setup grids for integration
425 :
426 : void SigmaPartialWave::setupGrid() {
427 : // Reset sigma maximum
428 0 : sigElMax = 0.;
429 :
430 : // Go through each subprocess
431 0 : gridMax.resize(subprocessMax);
432 0 : gridNorm.resize(subprocessMax);
433 0 : for (int sp = 0; sp < subprocessMax; sp++) {
434 : // Setup subprocess
435 0 : setSubprocess(sp);
436 :
437 : // Bins in Wcm
438 0 : int nBin1 = int( (binMax - mA - mB) / WCMBIN );
439 0 : gridMax[subprocess].resize(nBin1);
440 0 : gridNorm[subprocess].resize(nBin1);
441 0 : for (int n1 = 0; n1 < nBin1; n1++) {
442 : // Bin lower and upper
443 0 : double bl1 = mA + mB + double(n1) * WCMBIN;
444 0 : double bu1 = bl1 + WCMBIN;
445 :
446 : // Bins in cos(theta)
447 : int nBin2 = int( 2. / CTBIN );
448 0 : gridMax[subprocess][n1].resize(nBin2);
449 0 : for (int n2 = 0; n2 < nBin2; n2++) {
450 : // Bin lower and upper
451 0 : double bl2 = -1. + double(n2) * CTBIN;
452 0 : double bu2 = bl2 + CTBIN;
453 :
454 : // Find maximum
455 0 : double maxSig = 0.;
456 : double bl3 = bl1, bu3 = bu1, bl4 = bl2, bu4 = bu2;
457 0 : for (int iter = 0; iter < ITER; iter++) {
458 : int i3Save = -1, i4Save = -1;
459 0 : double step3 = (bu3 - bl3) / double(SUBBIN);
460 0 : double step4 = (bu4 - bl4) / double(SUBBIN);
461 0 : for (int i3 = 0; i3 <= SUBBIN; i3++) {
462 0 : double Wcm = bl3 + double(i3) * step3;
463 0 : for (int i4 = 0; i4 <= SUBBIN; i4++) {
464 0 : double ct = bl4 + double(i4) * step4;
465 0 : double ds = dSigma(Wcm, ct);
466 0 : if (ds > maxSig) {
467 : i3Save = i3;
468 : i4Save = i4;
469 0 : maxSig = ds;
470 0 : }
471 : }
472 : }
473 : // Set new min/max
474 0 : if (i3Save == -1 && i4Save == -1) break;
475 0 : if (i3Save > -1) {
476 0 : bl3 = bl3 + ((i3Save == 0) ? 0. : i3Save - 1.) * step3;
477 0 : bu3 = bl3 + ((i3Save == SUBBIN) ? 1. : 2.) * step3;
478 0 : }
479 0 : if (i4Save > -1) {
480 0 : bl4 = bl4 + ((i4Save == 0) ? 0. : i4Save - 1.) * step4;
481 0 : bu4 = bl4 + ((i4Save == SUBBIN) ? 1. : 2.) * step4;
482 0 : }
483 0 : } // for (iter)
484 :
485 : // Save maximum value
486 0 : gridMax[subprocess][n1][n2] = maxSig * (1. + GRIDSAFETY);
487 0 : gridNorm[subprocess][n1] += maxSig * (1. + GRIDSAFETY) * CTBIN;
488 0 : sigElMax = max(sigElMax, maxSig);
489 :
490 0 : } // for (n2)
491 : } // for (n1)
492 : } // for (sp)
493 :
494 0 : return;
495 : }
496 :
497 :
498 : //--------------------------------------------------------------------------
499 :
500 : // Pick a cos(theta) value
501 :
502 : double SigmaPartialWave::pickCosTheta(double Wcm) {
503 : // Find grid bin in Wcm
504 0 : int WcmBin = int((Wcm - mA - mB) / WCMBIN);
505 0 : if (WcmBin < 0) WcmBin = 0;
506 0 : if (WcmBin >= int(gridMax[subprocess].size()))
507 0 : WcmBin = int(gridMax[subprocess].size() - 1);
508 :
509 : // Pick a value of cos(theta)
510 : double ct, wgt;
511 :
512 0 : do {
513 : // Sample from overestimate and inverse
514 0 : double y = rndmPtr->flat() * gridNorm[subprocess][WcmBin];
515 : double sum = 0.;
516 : int ctBin;
517 0 : for (ctBin = 0; ctBin < int(2. / CTBIN); ctBin++) {
518 0 : if (sum + CTBIN * gridMax[subprocess][WcmBin][ctBin] > y) break;
519 0 : sum += CTBIN * gridMax[subprocess][WcmBin][ctBin];
520 : }
521 :
522 : // Linear interpolation
523 0 : double x1 = -1. + CTBIN * double(ctBin);
524 : double y1 = sum;
525 0 : double x2 = x1 + CTBIN;
526 0 : double y2 = sum + CTBIN * gridMax[subprocess][WcmBin][ctBin];
527 0 : ct = (x2 - x1) / (y2 - y1) * (y - y1) + x1;
528 0 : wgt = dSigma(Wcm, ct) / gridMax[subprocess][WcmBin][ctBin];
529 0 : if (wgt >= 1.) {
530 0 : infoPtr->errorMsg("Warning in SigmaPartialWave::pickCosTheta: "
531 : "weight above unity");
532 0 : break;
533 : }
534 0 : } while (wgt <= rndmPtr->flat());
535 :
536 0 : return ct;
537 0 : }
538 :
539 : //--------------------------------------------------------------------------
540 :
541 : // Set subprocess
542 :
543 : bool SigmaPartialWave::setSubprocess(int spIn) {
544 0 : if (sp2in.find(spIn) == sp2in.end()) return false;
545 0 : subprocess = spIn;
546 0 : pair < int, int > in = sp2in[spIn];
547 0 : idA = in.first;
548 0 : mA = particleDataPtr->m0(idA);
549 0 : idB = in.second;
550 0 : mB = particleDataPtr->m0(idB);
551 : return true;
552 0 : }
553 :
554 : bool SigmaPartialWave::setSubprocess(int idAin, int idBin) {
555 0 : pair < int, int > in = pair < int, int > (idAin, idBin);
556 0 : if (in2sp.find(in) == in2sp.end()) {
557 : // Try the other way around as well
558 0 : swap(in.first, in.second);
559 0 : if (in2sp.find(in) == in2sp.end()) return false;
560 : }
561 0 : subprocess = in2sp[in];
562 0 : idA = idAin;
563 0 : mA = particleDataPtr->m0(idA);
564 0 : idB = idBin;
565 0 : mB = particleDataPtr->m0(idB);
566 0 : return true;
567 0 : }
568 :
569 : //--------------------------------------------------------------------------
570 :
571 : // Calculate: mode = 0 (sigma elastic), 1 (sigma total), 2 (dSigma/dcTheta)
572 :
573 : double SigmaPartialWave::sigma(int mode, double Wcm, double cTheta) {
574 : // Below threshold, return 0
575 0 : if (Wcm < (mA + mB + MASSSAFETY)) return 0.;
576 :
577 : // Return values
578 0 : complex amp[2] = { complex(0., 0.) };
579 : double sig = 0.;
580 :
581 : // Kinematic variables
582 0 : double s = pow2(Wcm);
583 0 : double k2 = (s - pow2(mB + mA)) * (s - pow2(mB - mA)) / 4. / s;
584 :
585 : // Precompute all required Pl and Pl' values
586 0 : double sTheta = 0.;
587 0 : if (mode == 2) {
588 0 : if (process == 2) sTheta = sqrt(1. - pow2(cTheta));
589 0 : legendreP(cTheta, ((process == 2) ? true : false));
590 0 : }
591 :
592 : // Loop over L
593 0 : for (int L = 0; L < Lmax; L++) {
594 :
595 : // Loop over J (only J = 0 for spinless)
596 0 : complex ampJ[2] = { complex(0., 0.) };
597 0 : int Jstart = (process != 2) ? 0 : 2 * L - 1;
598 0 : int Jend = (process != 2) ? 1 : 2 * L + 2;
599 0 : int Jstep = (process != 2) ? 1 : 2;
600 : int Jcount = 0;
601 0 : for (int J = Jstart; J < Jend; J += Jstep, Jcount++) {
602 :
603 : // Loop over isospin coefficients
604 0 : for (int I = 0; I < Imax; I++) {
605 0 : if (isoCoeff[subprocess][I] == 0.) continue;
606 :
607 : // Check wave exists
608 0 : int LIJ = L * LSHIFT + I * ISHIFT + J;
609 0 : if (pwData.find(LIJ) == pwData.end()) continue;
610 :
611 : // Extrapolation / interpolation (not for last bin)
612 0 : map < double, complex >::iterator it = pwData[LIJ].upper_bound(Wcm);
613 0 : if (it == pwData[LIJ].begin()) ++it;
614 : double ar, ai;
615 0 : if (it == pwData[LIJ].end()) {
616 0 : ar = (--it)->second.real();
617 0 : ai = it->second.imag();
618 0 : } else {
619 0 : double eA = it->first;
620 0 : complex ampA = (it--)->second;
621 0 : double eB = it->first;
622 0 : complex ampB = it->second;
623 :
624 0 : ar = (ampA.real() - ampB.real()) / (eA - eB) *
625 0 : (Wcm - eB) + ampB.real();
626 0 : ai = (ampA.imag() - ampB.imag()) / (eA - eB) *
627 0 : (Wcm - eB) + ampB.imag();
628 0 : }
629 :
630 : // Isospin sum
631 0 : ampJ[Jcount] += isoCoeff[subprocess][I] * complex(ar, ai);
632 0 : }
633 : }
634 :
635 : // Partial wave sum. Sigma elastic
636 0 : if (mode == 0) {
637 0 : if (process == 0 || process == 1) {
638 0 : sig += (2. * L + 1.) * (ampJ[0] * conj(ampJ[0])).real();
639 0 : } else if (process == 2) {
640 0 : sig += ( (L + 0.) * (ampJ[0] * conj(ampJ[0])).real() +
641 0 : (L + 1.) * (ampJ[1] * conj(ampJ[1])).real() );
642 0 : }
643 :
644 : // Sigma total
645 0 : } else if (mode == 1) {
646 0 : if (process == 0 || process == 1) {
647 0 : sig += (2. * L + 1.) * ampJ[0].imag();
648 0 : } else if (process == 2) {
649 0 : sig += ( (L + 0.) * ampJ[0].imag() + (L + 1.) * ampJ[1].imag() );
650 0 : }
651 :
652 : // dSigma
653 0 : } else if (mode == 2) {
654 0 : if (process == 0 || process == 1) {
655 0 : amp[0] += (2. * L + 1.) * ampJ[0] * PlVec[L];
656 0 : } else if (process == 2) {
657 0 : amp[0] += ((L + 0.) * ampJ[0] + double(L + 1.) * ampJ[1]) * PlVec[L];
658 0 : amp[1] += complex(0., 1.) * (ampJ[1] - ampJ[0]) * sTheta * PlpVec[L];
659 0 : }
660 : }
661 :
662 0 : } // for (L)
663 :
664 : // Normalisation and return
665 0 : if (mode == 0 || mode == 1) {
666 0 : if (norm == 0) sig *= 4. * M_PI / k2 * CONVERT2MB;
667 0 : else if (norm == 1) sig *= 64. * M_PI / s * CONVERT2MB;
668 :
669 0 : } else if (mode == 2) {
670 0 : sig = (amp[0] * conj(amp[0])).real() + (amp[1] * conj(amp[1])).real();
671 0 : if (norm == 0) sig *= 2. * M_PI / k2 * CONVERT2MB;
672 0 : else if (norm == 1) sig *= 32. * M_PI / s * CONVERT2MB;
673 : }
674 : // Half for identical
675 0 : return ((idA == idB) ? 0.5 : 1.) * sig;
676 0 : }
677 :
678 :
679 : //--------------------------------------------------------------------------
680 :
681 : // Bonnet's recursion formula for Legendre polynomials and derivatives
682 :
683 : void SigmaPartialWave::legendreP(double ct, bool deriv) {
684 0 : if (Lmax > 1) PlVec[1] = ct;
685 0 : for (int L = 2; L < Lmax; L++) {
686 0 : PlVec[L] = ( (2. * L - 1.) * ct * PlVec[L - 1] -
687 0 : (L - 1.) * PlVec[L - 2] ) / double(L);
688 0 : if (deriv)
689 0 : PlpVec[L] = ( (2. * L - 1.) * (PlVec[L - 1] + ct * PlpVec[L - 1]) -
690 0 : (L - 1.) * PlpVec[L - 2] ) / double(L);
691 : }
692 0 : return;
693 : }
694 :
695 :
696 : //==========================================================================
697 :
698 : // HadronScatter class
699 :
700 : //--------------------------------------------------------------------------
701 :
702 : // Perform initialization and store pointers.
703 :
704 : bool HadronScatter::init(Info* infoPtrIn, Settings& settings,
705 : Rndm *rndmPtrIn, ParticleData *particleDataPtr) {
706 : // Save incoming pointers
707 0 : infoPtr = infoPtrIn;
708 0 : rndmPtr = rndmPtrIn;
709 :
710 : // Main settings
711 0 : doHadronScatter = settings.flag("HadronScatter:scatter");
712 0 : afterDecay = settings.flag("HadronScatter:afterDecay");
713 0 : allowDecayProd = settings.flag("HadronScatter:allowDecayProd");
714 0 : scatterRepeat = settings.flag("HadronScatter:scatterRepeat");
715 : // Hadron selection
716 0 : hadronSelect = settings.mode("HadronScatter:hadronSelect");
717 0 : Npar = settings.parm("HadronScatter:N");
718 0 : kPar = settings.parm("HadronScatter:k");
719 0 : pPar = settings.parm("HadronScatter:p");
720 : // Scattering probability
721 0 : scatterProb = settings.mode("HadronScatter:scatterProb");
722 0 : jPar = settings.parm("HadronScatter:j");
723 0 : rMax = settings.parm("HadronScatter:rMax");
724 0 : rMax2 = rMax * rMax;
725 0 : doTile = settings.flag("HadronScatter:tile");
726 :
727 : // String fragmentation and MPI settings
728 0 : pTsigma = 2.0 * settings.parm("StringPT:sigma");
729 0 : pTsigma2 = pTsigma * pTsigma;
730 0 : double pT0ref = settings.parm("MultipartonInteractions:pT0ref");
731 0 : double eCMref = settings.parm("MultipartonInteractions:eCMref");
732 0 : double eCMpow = settings.parm("MultipartonInteractions:eCMpow");
733 0 : double eCMnow = infoPtr->eCM();
734 0 : pT0MPI = pT0ref * pow(eCMnow / eCMref, eCMpow);
735 :
736 : // Tiling
737 0 : double mp2 = particleDataPtr->m0(111) * particleDataPtr->m0(111);
738 0 : double eA = infoPtr->eA();
739 0 : double eB = infoPtr->eB();
740 0 : double pzA = sqrt(eA * eA - mp2);
741 0 : double pzB = -sqrt(eB * eB - mp2);
742 0 : yMax = 0.5 * log((eA + pzA) / (eA - pzA));
743 0 : yMin = 0.5 * log((eB + pzB) / (eB - pzB));
744 : // Size in y and phi
745 0 : if (doTile) {
746 0 : ytMax = int((yMax - yMin) / rMax);
747 0 : ytSize = (yMax - yMin) / double(ytMax);
748 0 : ptMax = int(2. * M_PI / rMax);
749 0 : ptSize = 2. * M_PI / double(ptMax);
750 0 : } else {
751 0 : ytMax = 1;
752 0 : ytSize = yMax - yMin;
753 0 : ptMax = 1;
754 0 : ptSize = 2. * M_PI;
755 : }
756 : // Initialise tiles
757 0 : tile.resize(ytMax);
758 0 : for (int yt = 0; yt < ytMax; yt++) tile[yt].resize(ptMax);
759 :
760 : // Find path to data files, i.e. xmldoc directory location.
761 : // Environment variable takes precedence, else use constructor input.
762 : // XXX - as in Pythia.cc, but not passed around in e.g. Info/Settings,
763 : // so redo here
764 0 : string xmlPath = "";
765 : const char* PYTHIA8DATA = "PYTHIA8DATA";
766 0 : char* envPath = getenv(PYTHIA8DATA);
767 0 : if (envPath != 0 && *envPath != '\0') {
768 : int i = 0;
769 0 : while (*(envPath+i) != '\0') xmlPath += *(envPath+(i++));
770 0 : } else xmlPath = "../xmldoc";
771 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
772 :
773 : // Hadron scattering partial wave cross sections
774 0 : if ( !sigmaPW[0].init(0, xmlPath, "pipi-Froggatt.dat",
775 0 : infoPtr, particleDataPtr, rndmPtr) ) return false;
776 0 : if ( !sigmaPW[1].init(1, xmlPath, "piK-Estabrooks.dat",
777 0 : infoPtr, particleDataPtr, rndmPtr) ) return false;
778 0 : if ( !sigmaPW[2].init(2, xmlPath, "piN-SAID-WI08.dat",
779 0 : infoPtr, particleDataPtr, rndmPtr) ) return false;
780 0 : sigElMax = 0.;
781 0 : sigElMax = max(sigElMax, sigmaPW[0].getSigmaElMax());
782 0 : sigElMax = max(sigElMax, sigmaPW[1].getSigmaElMax());
783 0 : sigElMax = max(sigElMax, sigmaPW[2].getSigmaElMax());
784 :
785 : // DEBUG
786 0 : debugOutput();
787 :
788 0 : return true;
789 0 : }
790 :
791 :
792 : //--------------------------------------------------------------------------
793 :
794 : // Debug output
795 :
796 : void HadronScatter::debugOutput() {
797 : // Print settings
798 0 : cout << "Hadron scattering:" << endl
799 0 : << " scatter = " << ((doHadronScatter) ? "on" : "off") << endl
800 0 : << " afterDecay = " << ((afterDecay) ? "on" : "off") << endl
801 0 : << " allowDecayProd = " << ((allowDecayProd) ? "on" : "off") << endl
802 0 : << " scatterRepeat = " << ((scatterRepeat) ? "on" : "off") << endl
803 0 : << " tile = " << ((doTile) ? "on" : "off") << endl
804 0 : << " yMin = " << yMin << endl
805 0 : << " yMax = " << yMax << endl
806 0 : << " ytMax = " << ytMax << endl
807 0 : << " ytSize = " << ytSize << endl
808 0 : << " ptMax = " << ptMax << endl
809 0 : << " ptSize = " << ptSize << endl
810 0 : << endl
811 0 : << " hadronSelect = " << hadronSelect << endl
812 0 : << " N = " << Npar << endl
813 0 : << " k = " << kPar << endl
814 0 : << " p = " << pPar << endl
815 0 : << endl
816 0 : << " scatterProb = " << scatterProb << endl
817 0 : << " j = " << jPar << endl
818 0 : << " rMax = " << rMax << endl
819 0 : << endl
820 0 : << " pTsigma = " << pTsigma2 << endl
821 0 : << " pT0MPI = " << pT0MPI << endl
822 0 : << endl
823 0 : << " sigElMax = " << sigElMax << endl << endl;
824 :
825 0 : return;
826 : }
827 :
828 :
829 : //--------------------------------------------------------------------------
830 :
831 : // Perform hadron scattering
832 :
833 : void HadronScatter::scatter(Event& event) {
834 : // Reset tiles
835 0 : for (int yt = 0; yt < ytMax; yt++)
836 0 : for (int pt = 0; pt < ptMax; pt++)
837 0 : tile[yt][pt].clear();
838 :
839 : // Generate list of hadrons which can take part in the scattering
840 0 : for (int i = 0; i < event.size(); i++)
841 0 : if (event[i].isFinal() && event[i].isHadron() && canScatter(event, i))
842 0 : tile[yTile(event, i)][pTile(event, i)].insert(HSIndex(i, i));
843 :
844 : // Generate all pairwise interaction probabilities
845 0 : vector < HadronScatterPair > scatterList;
846 : // For each tile and for each hadron in the tile do pairing
847 0 : for (int pt1 = 0; pt1 < ptMax; pt1++)
848 0 : for (int yt1 = 0; yt1 < ytMax; yt1++)
849 0 : for (set < HSIndex >::iterator si1 = tile[yt1][pt1].begin();
850 0 : si1 != tile[yt1][pt1].end(); si1++)
851 0 : tileIntProb(scatterList, event, *si1, yt1, pt1, false);
852 : // Sort by ordering measure (largest to smallest)
853 0 : sort(scatterList.rbegin(), scatterList.rend());
854 :
855 : // Reset list of things that have scattered
856 0 : if (scatterRepeat) scattered.clear();
857 :
858 : // Do scatterings
859 0 : while (scatterList.size() > 0) {
860 : // Check still valid and scatter
861 0 : HadronScatterPair &hsp = scatterList[0];
862 0 : if (!event[hsp.i1.second].isFinal() || !event[hsp.i2.second].isFinal()) {
863 0 : scatterList.erase(scatterList.begin());
864 0 : continue;
865 : }
866 : // Remove old entries in tiles and scatter
867 0 : tile[hsp.yt1][hsp.pt1].erase(hsp.i1);
868 0 : tile[hsp.yt2][hsp.pt2].erase(hsp.i2);
869 0 : hadronScatter(event, hsp);
870 : // Store new indices for later
871 0 : HSIndex iNew1 = hsp.i1, iNew2 = hsp.i2;
872 :
873 : // Check if hadrons can scatter again
874 : bool resort = false;
875 0 : if (scatterRepeat) {
876 : // Check for new scatters of iNew1 and iNew2
877 0 : HSIndex iNew = iNew1;
878 0 : for (int i = 0; i < 2; i++) {
879 0 : if (canScatter(event, iNew.second)) {
880 : // If both can scatter again, make sure they can't scatter
881 : // with each other
882 0 : if (i == 1) scattered.insert(HSIndex(min(iNew1.first, iNew2.first),
883 0 : max(iNew1.first, iNew2.first)));
884 0 : int yt = yTile(event, iNew.second);
885 0 : int pt = pTile(event, iNew.second);
886 0 : resort = tileIntProb(scatterList, event, iNew, yt, pt, true);
887 0 : tile[yt][pt].insert(iNew);
888 0 : }
889 0 : iNew = iNew2;
890 : } // for (i)
891 :
892 0 : } // if (scatterRepeat)
893 :
894 : // Remove the old entry and onto next
895 0 : scatterList.erase(scatterList.begin());
896 : // Resort list if anything added
897 0 : if (resort) sort(scatterList.rbegin(), scatterList.rend());
898 :
899 0 : } // while (scatterList.size() > 0)
900 :
901 : // Done.
902 : return;
903 0 : }
904 :
905 :
906 : //--------------------------------------------------------------------------
907 :
908 : // Criteria for if a hadron is allowed to scatter or not
909 :
910 : bool HadronScatter::canScatter(Event& event, int i) {
911 : // Probability to accept
912 : double p = 0.;
913 :
914 : // Pions, K+, K-, p+, pbar- only
915 0 : if (scatterProb == 1 || scatterProb == 2)
916 0 : if (event[i].idAbs() != 111 && event[i].idAbs() != 211 &&
917 0 : event[i].idAbs() != 321 && event[i].idAbs() != 2212)
918 0 : return false;
919 :
920 : // Probability
921 0 : switch (hadronSelect) {
922 : case 0:
923 0 : double t1 = exp( - event[i].pT2() / 2. / pTsigma2);
924 0 : double t2 = pow(pT0MPI, pPar) /
925 0 : pow(pT0MPI * pT0MPI + event[i].pT2(), pPar / 2.);
926 0 : p = Npar * t1 / ( (1 - kPar) * t1 + kPar * t2 );
927 : break;
928 : }
929 :
930 : // Return true/false
931 0 : if (p > rndmPtr->flat()) return true;
932 0 : else return false;
933 0 : }
934 :
935 :
936 : //--------------------------------------------------------------------------
937 :
938 : // Probability for scattering
939 : bool HadronScatter::doesScatter(Event& event, const HSIndex &i1,
940 : const HSIndex &i2) {
941 0 : Particle &p1 = event[i1.second];
942 0 : Particle &p2 = event[i2.second];
943 :
944 : // Can't come from the same decay
945 0 : if (!allowDecayProd && event[i1.first].mother1() ==
946 0 : event[i2.first].mother1() && event[event[i1.first].mother1()].isHadron())
947 0 : return false;
948 :
949 : // Check that the two hadrons have not already scattered
950 0 : if (scatterRepeat &&
951 0 : scattered.find(HSIndex(min(i1.first, i2.first), max(i1.first, i2.first)))
952 0 : != scattered.end()) return false;
953 :
954 : // K-K, p-p and K-p not allowed
955 0 : int id1 = min(p1.idAbs(), p2.idAbs());
956 0 : int id2 = max(p1.idAbs(), p2.idAbs());
957 0 : if (scatterProb == 1 || scatterProb == 2) {
958 0 : if ((id1 == 321 || id1 == 2212) && id1 == id2) return false;
959 0 : if (id1 == 321 && id2 == 2212) return false;
960 : }
961 :
962 : // Distance in y - phi space
963 0 : double dy = p1.y() - p2.y();
964 0 : double dp = abs(p1.phi() - p2.phi());
965 0 : if (dp > M_PI) dp = 2 * M_PI - dp;
966 0 : double dr2 = dy * dy + dp * dp;
967 0 : double p = max(0., 1. - dr2 / rMax2);
968 :
969 : // Additional factor depending on scatterProb
970 0 : if (scatterProb == 0 || scatterProb == 1) {
971 0 : p *= jPar;
972 :
973 : // Additional pair dependent probability
974 0 : } else if (scatterProb == 2) {
975 : // Wcm and which partial wave amplitude to use
976 0 : double Wcm = (p1.p() + p2.p()).mCalc();
977 : int pw = 0;
978 0 : if ((id1 == 111 || id1 == 211) && (id2 == 111 || id2 == 211)) pw = 0;
979 0 : else if ((id1 == 111 || id1 == 211) && id2 == 321) pw = 1;
980 0 : else if ((id1 == 111 || id1 == 211) && id2 == 2212) pw = 2;
981 : else {
982 0 : infoPtr->errorMsg("Error in HadronScatter::doesScatter:"
983 : "unknown subprocess");
984 : }
985 0 : if (!sigmaPW[pw].setSubprocess(p1.id(), p2.id())) {
986 0 : infoPtr->errorMsg("Error in HadronScatter::doesScatter:"
987 : "setSubprocess failed");
988 0 : } else {
989 0 : p *= 1 - exp(-jPar * sigmaPW[pw].sigmaEl(Wcm));
990 : }
991 0 : }
992 :
993 : // Return true/false
994 0 : return (p > rndmPtr->flat());
995 0 : }
996 :
997 :
998 : //--------------------------------------------------------------------------
999 :
1000 : // Calculate ordering measure
1001 :
1002 : double HadronScatter::measure(Event& event, int idx1, int idx2) {
1003 0 : Particle &p1 = event[idx1];
1004 0 : Particle &p2 = event[idx2];
1005 0 : return abs(p1.pT() / p1.mT() - p2.pT() / p2.mT());
1006 : }
1007 :
1008 :
1009 : //--------------------------------------------------------------------------
1010 :
1011 : // Scatter a pair
1012 :
1013 : bool HadronScatter::hadronScatter(Event& event, HadronScatterPair &hsp) {
1014 0 : bool doSwap = (0.5 < rndmPtr->flat()) ? true : false;
1015 0 : if (doSwap) swap(hsp.i1, hsp.i2);
1016 :
1017 0 : Particle &p1 = event[hsp.i1.second];
1018 0 : Particle &p2 = event[hsp.i2.second];
1019 :
1020 : // Pick theta and phi (always flat)
1021 0 : double ct = 0., phi = 2 * M_PI * rndmPtr->flat();
1022 :
1023 : // Flat for all flavours
1024 0 : if (scatterProb == 0 || scatterProb == 1) {
1025 0 : ct = 2. * rndmPtr->flat() - 1.;
1026 :
1027 : // pi-pi, pi-K and pi-p only using partial wave amplitudes
1028 0 : } else if (scatterProb == 2) {
1029 : // Wcm and which partial wave amplitude to use
1030 0 : int id1 = min(p1.idAbs(), p2.idAbs());
1031 0 : int id2 = max(p1.idAbs(), p2.idAbs());
1032 0 : double Wcm = (p1.p() + p2.p()).mCalc();
1033 : int pw = 0;
1034 0 : if ((id1 == 111 || id1 == 211) && (id2 == 111 || id2 == 211)) pw = 0;
1035 0 : else if ((id1 == 111 || id1 == 211) && id2 == 321) pw = 1;
1036 0 : else if ((id1 == 111 || id1 == 211) && id2 == 2212) pw = 2;
1037 : else {
1038 0 : infoPtr->errorMsg("Error in HadronScatter::hadronScatter:"
1039 : "unknown subprocess");
1040 : }
1041 0 : sigmaPW[pw].setSubprocess(p1.id(), p2.id());
1042 0 : ct = sigmaPW[pw].pickCosTheta(Wcm);
1043 0 : }
1044 :
1045 : // Rotation
1046 0 : RotBstMatrix sMat;
1047 0 : sMat.toCMframe(p1.p(), p2.p());
1048 0 : sMat.rot(acos(ct), phi);
1049 0 : sMat.fromCMframe(p1.p(), p2.p());
1050 0 : Vec4 v1 = p1.p(), v2 = p2.p();
1051 0 : v1.rotbst(sMat);
1052 0 : v2.rotbst(sMat);
1053 :
1054 : // Update event record
1055 0 : int iNew1 = event.copy(hsp.i1.second, 98);
1056 0 : event[iNew1].p(v1);
1057 0 : event[iNew1].e(event[iNew1].eCalc());
1058 0 : event[hsp.i1.second].statusNeg();
1059 0 : int iNew2 = event.copy(hsp.i2.second, 98);
1060 0 : event[iNew2].p(v2);
1061 0 : event[iNew2].e(event[iNew2].eCalc());
1062 0 : event[hsp.i2.second].statusNeg();
1063 :
1064 : // New indices in event record
1065 0 : hsp.i1.second = iNew1;
1066 0 : hsp.i2.second = iNew2;
1067 :
1068 0 : if (doSwap) swap(hsp.i1, hsp.i2);
1069 0 : return true;
1070 0 : }
1071 :
1072 :
1073 : //--------------------------------------------------------------------------
1074 :
1075 : // Calculate pair interaction probabilities in nearest-neighbour tiles
1076 : // (yt1, pt1) represents centre cell (8):
1077 : // 7 | 0 | 1
1078 : // -----------
1079 : // 6 | 8 | 2
1080 : // -----------
1081 : // 5 | 4 | 3
1082 :
1083 : bool HadronScatter::tileIntProb(vector < HadronScatterPair > &hsp,
1084 : Event &event, const HSIndex &i1, int yt1, int pt1, bool doAll) {
1085 : // Track if a new pair is added
1086 : bool pairAdded = false;
1087 :
1088 : // Special handling for pairing in own tile
1089 0 : if (!doAll) {
1090 0 : set < HSIndex >::iterator si2 = tile[yt1][pt1].find(i1);
1091 0 : while (++si2 != tile[yt1][pt1].end())
1092 : // Calculate if scatters
1093 0 : if (doesScatter(event, i1, *si2)) {
1094 0 : double m = measure(event, i1.second, si2->second);
1095 0 : hsp.push_back(HadronScatterPair(i1, yt1, pt1, *si2, yt1, pt1, m));
1096 0 : }
1097 0 : }
1098 :
1099 : // And the rest
1100 0 : int tileMax = (doAll) ? 9 : 4;
1101 0 : for (int tileNow = 0; tileNow < tileMax; tileNow++) {
1102 : // New (yt, pt) coordinate
1103 : int yt2 = yt1, pt2 = pt1;
1104 0 : switch (tileNow) {
1105 0 : case 0: ++pt2; break;
1106 0 : case 1: ++yt2; ++pt2; break;
1107 0 : case 2: ++yt2; break;
1108 0 : case 3: ++yt2; --pt2; break;
1109 0 : case 4: --pt2; break;
1110 0 : case 5: --yt2; --pt2; break;
1111 0 : case 6: --yt2; break;
1112 0 : case 7: --yt2; ++pt2; break;
1113 : }
1114 :
1115 : // Limit in rapidity tiles
1116 0 : if (yt2 < 0 || yt2 >= ytMax) continue;
1117 : // Wraparound for phi tiles
1118 0 : if (pt2 < 0) {
1119 0 : pt2 = ptMax - 1;
1120 0 : if (pt2 == pt1 || pt2 == pt1 + 1) continue;
1121 :
1122 0 : } else if (pt2 >= ptMax) {
1123 : pt2 = 0;
1124 0 : if (pt2 == pt1 || pt2 == pt1 - 1) continue;
1125 : }
1126 :
1127 : // Interaction probability
1128 0 : for (set < HSIndex >::iterator si2 = tile[yt2][pt2].begin();
1129 0 : si2 != tile[yt2][pt2].end(); si2++) {
1130 : // Calculate if scatters
1131 0 : if (doesScatter(event, i1, *si2)) {
1132 0 : double m = measure(event, i1.second, si2->second);
1133 0 : hsp.push_back(HadronScatterPair(i1, yt1, pt1, *si2, yt2, pt2, m));
1134 : pairAdded = true;
1135 0 : }
1136 : }
1137 0 : }
1138 :
1139 0 : return pairAdded;
1140 : }
1141 :
1142 : //==========================================================================
1143 :
1144 : } // namespace Pythia8
|
