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