StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BoseEinstein.cc
1 // BoseEinstein.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for the BoseEinsten class.
7 
8 #include "Pythia8/BoseEinstein.h"
9 
10 namespace Pythia8 {
11 
12 //==========================================================================
13 
14 // The BoseEinstein class.
15 
16 //--------------------------------------------------------------------------
17 
18 // Constants: could be changed here if desired, but normally should not.
19 // These are of technical nature, as described for each.
20 
21 // Enumeration of id codes and table for particle species considered.
22 const int BoseEinstein::IDHADRON[9] = { 211, -211, 111, 321, -321,
23  130, 310, 221, 331 };
24 const int BoseEinstein::ITABLE[9] = { 0, 0, 0, 1, 1, 1, 1, 2, 3 };
25 
26 // Distance between table entries, normalized to min( 2*mass, QRef).
27 const double BoseEinstein::STEPSIZE = 0.05;
28 
29 // Skip shift for two extremely close particles, to avoid instabilities.
30 const double BoseEinstein::Q2MIN = 1e-8;
31 
32 // Parameters of energy compensation procedure: maximally allowed
33 // relative energy error, iterative stepsize, and number of iterations.
34 const double BoseEinstein::COMPRELERR = 1e-10;
35 const double BoseEinstein::COMPFACMAX = 1000.;
36 const int BoseEinstein::NCOMPSTEP = 10;
37 
38 //--------------------------------------------------------------------------
39 
40 // Find settings. Precalculate table used to find momentum shifts.
41 
42 bool BoseEinstein::init() {
43 
44  // Main flags.
45  doPion = flag("BoseEinstein:Pion");
46  doKaon = flag("BoseEinstein:Kaon");
47  doEta = flag("BoseEinstein:Eta");
48 
49  // Shape of Bose-Einstein enhancement/suppression.
50  lambda = parm("BoseEinstein:lambda");
51  QRef = parm("BoseEinstein:QRef");
52 
53  // Multiples and inverses (= "radii") of distance parameters in Q-space.
54  QRef2 = 2. * QRef;
55  QRef3 = 3. * QRef;
56  R2Ref = 1. / (QRef * QRef);
57  R2Ref2 = 1. / (QRef2 * QRef2);
58  R2Ref3 = 1. / (QRef3 * QRef3);
59 
60  // Masses of particles with Bose-Einstein implemented.
61  for (int iSpecies = 0; iSpecies < 9; ++iSpecies)
62  mHadron[iSpecies] = particleDataPtr->m0( IDHADRON[iSpecies] );
63 
64  // Pair pi, K, eta and eta' masses for use in tables.
65  mPair[0] = 2. * mHadron[0];
66  mPair[1] = 2. * mHadron[3];
67  mPair[2] = 2. * mHadron[7];
68  mPair[3] = 2. * mHadron[8];
69 
70  // Loop over the four required tables. Local variables.
71  double Qnow, Q2now, centerCorr;
72  for (int iTab = 0; iTab < 4; ++iTab) {
73  m2Pair[iTab] = mPair[iTab] * mPair[iTab];
74 
75  // Step size and number of steps in normal table.
76  deltaQ[iTab] = STEPSIZE * min(mPair[iTab], QRef);
77  nStep[iTab] = min( 199, 1 + int(3. * QRef / deltaQ[iTab]) );
78  maxQ[iTab] = (nStep[iTab] - 0.1) * deltaQ[iTab];
79  centerCorr = deltaQ[iTab] * deltaQ[iTab] / 12.;
80 
81  // Construct normal table recursively in Q space.
82  shift[iTab][0] = 0.;
83  for (int i = 1; i <= nStep[iTab]; ++i) {
84  Qnow = deltaQ[iTab] * (i - 0.5);
85  Q2now = Qnow * Qnow;
86  shift[iTab][i] = shift[iTab][i - 1] + exp(-Q2now * R2Ref)
87  * deltaQ[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
88  }
89 
90  // Step size and number of steps in compensation table.
91  deltaQ3[iTab] = STEPSIZE * min(mPair[iTab], QRef3);
92  nStep3[iTab] = min( 199, 1 + int(9. * QRef / deltaQ3[iTab]) );
93  maxQ3[iTab] = (nStep3[iTab] - 0.1) * deltaQ3[iTab];
94  centerCorr = deltaQ3[iTab] * deltaQ3[iTab] / 12.;
95 
96  // Construct compensation table recursively in Q space.
97  shift3[iTab][0] = 0.;
98  for (int i = 1; i <= nStep3[iTab]; ++i) {
99  Qnow = deltaQ3[iTab] * (i - 0.5);
100  Q2now = Qnow * Qnow;
101  shift3[iTab][i] = shift3[iTab][i - 1] + exp(-Q2now * R2Ref3)
102  * deltaQ3[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
103  }
104 
105  }
106 
107  // Done.
108  return true;
109 
110 }
111 
112 //--------------------------------------------------------------------------
113 
114 // Perform Bose-Einstein corrections on an event.
115 
116 bool BoseEinstein::shiftEvent( Event& event) {
117 
118  // Reset list of identical particles.
119  hadronBE.resize(0);
120 
121  // Loop over all hadron species with BE effects.
122  nStored[0] = 0;
123  for (int iSpecies = 0; iSpecies < 9; ++iSpecies) {
124  nStored[iSpecies + 1] = nStored[iSpecies];
125  if (!doPion && iSpecies <= 2) continue;
126  if (!doKaon && iSpecies >= 3 && iSpecies <= 6) continue;
127  if (!doEta && iSpecies >= 7) continue;
128 
129  // Properties of current hadron species.
130  int idNow = IDHADRON[ iSpecies ];
131  int iTab = ITABLE[ iSpecies ];
132 
133  // Loop through event record to store copies of current species.
134  for (int i = 0; i < event.size(); ++i)
135  if ( event[i].id() == idNow && event[i].isFinal() )
136  hadronBE.push_back(
137  BoseEinsteinHadron( idNow, i, event[i].p(), event[i].m() ) );
138  nStored[iSpecies + 1] = hadronBE.size();
139 
140  // Loop through pairs of identical particles and find shifts.
141  for (int i1 = nStored[iSpecies]; i1 < nStored[iSpecies+1] - 1; ++i1)
142  for (int i2 = i1 + 1; i2 < nStored[iSpecies+1]; ++i2)
143  shiftPair( i1, i2, iTab);
144  }
145 
146  // Must have at least two pairs to carry out compensation.
147  if (nStored[9] < 2) return true;
148 
149  // Shift momenta and recalculate energies.
150  double eSumOriginal = 0.;
151  double eSumShifted = 0.;
152  double eDiffByComp = 0.;
153  for (int i = 0; i < nStored[9]; ++i) {
154  eSumOriginal += hadronBE[i].p.e();
155  hadronBE[i].p += hadronBE[i].pShift;
156  hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
157  eSumShifted += hadronBE[i].p.e();
158  eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
159  / hadronBE[i].p.e();
160  }
161 
162  // Iterate compensation shift until convergence.
163  int iStep = 0;
164  while ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal
165  && abs(eSumShifted - eSumOriginal) < COMPFACMAX * abs(eDiffByComp)
166  && iStep < NCOMPSTEP ) {
167  ++iStep;
168  double compFac = (eSumOriginal - eSumShifted) / eDiffByComp;
169  eSumShifted = 0.;
170  eDiffByComp = 0.;
171  for (int i = 0; i < nStored[9]; ++i) {
172  hadronBE[i].p += compFac * hadronBE[i].pComp;
173  hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
174  eSumShifted += hadronBE[i].p.e();
175  eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
176  / hadronBE[i].p.e();
177  }
178  }
179 
180  // Error if no convergence, and then return without doing BE shift.
181  // However, not grave enough to kill event, so return true.
182  if ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal ) {
183  infoPtr->errorMsg("Warning in BoseEinstein::shiftEvent: "
184  "no consistent BE shift topology found, so skip BE");
185  return true;
186  }
187 
188  // Store new particle copies with shifted momenta.
189  for (int i = 0; i < nStored[9]; ++i) {
190  int iNew = event.copy( hadronBE[i].iPos, 99);
191  event[ iNew ].p( hadronBE[i].p );
192  }
193 
194  // Done.
195  return true;
196 
197 }
198 
199 //--------------------------------------------------------------------------
200 
201 // Calculate shift and (unnormalized) compensation for pair.
202 
203 void BoseEinstein::shiftPair( int i1, int i2, int iTab) {
204 
205  // Calculate old relative momentum.
206  double Q2old = m2(hadronBE[i1].p, hadronBE[i2].p) - m2Pair[iTab];
207  if (Q2old < Q2MIN) return;
208  double Qold = sqrt(Q2old);
209  double psFac = sqrt(Q2old + m2Pair[iTab]) / Q2old;
210 
211  // Calculate new relative momentum for normal shift.
212  double Qmove = 0.;
213  if (Qold < deltaQ[iTab]) Qmove = Qold / 3.;
214  else if (Qold < maxQ[iTab]) {
215  double realQbin = Qold / deltaQ[iTab];
216  int intQbin = int( realQbin );
217  double inter = (pow3(realQbin) - pow3(intQbin))
218  / (3 * intQbin * (intQbin + 1) + 1);
219  Qmove = ( shift[iTab][intQbin] + inter * (shift[iTab][intQbin + 1]
220  - shift[iTab][intQbin]) ) * psFac;
221  }
222  else Qmove = shift[iTab][nStep[iTab]] * psFac;
223  double Q2new = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove), 2. / 3.);
224 
225  // Calculate corresponding three-momentum shift.
226  double Q2Diff = Q2new - Q2old;
227  double p2DiffAbs = (hadronBE[i1].p - hadronBE[i2].p).pAbs2();
228  double p2AbsDiff = hadronBE[i1].p.pAbs2() - hadronBE[i2].p.pAbs2();
229  double eSum = hadronBE[i1].p.e() + hadronBE[i2].p.e();
230  double eDiff = hadronBE[i1].p.e() - hadronBE[i2].p.e();
231  double sumQ2E = Q2Diff + eSum * eSum;
232  double rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
233  double rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
234  double factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
235  + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
236 
237  // Add shifts to sum. (Energy component dummy.)
238  Vec4 pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
239  hadronBE[i1].pShift += pDiff;
240  hadronBE[i2].pShift -= pDiff;
241 
242  // Calculate new relative momentum for compensation shift.
243  double Qmove3 = 0.;
244  if (Qold < deltaQ3[iTab]) Qmove3 = Qold / 3.;
245  else if (Qold < maxQ3[iTab]) {
246  double realQbin = Qold / deltaQ3[iTab];
247  int intQbin = int( realQbin );
248  double inter = (pow3(realQbin) - pow3(intQbin))
249  / (3 * intQbin * (intQbin + 1) + 1);
250  Qmove3 = ( shift3[iTab][intQbin] + inter * (shift3[iTab][intQbin + 1]
251  - shift3[iTab][intQbin]) ) * psFac;
252  }
253  else Qmove3 = shift3[iTab][nStep3[iTab]] *psFac;
254  double Q2new3 = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove3), 2. / 3.);
255 
256  // Calculate corresponding three-momentum shift.
257  Q2Diff = Q2new3 - Q2old;
258  sumQ2E = Q2Diff + eSum * eSum;
259  rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
260  rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
261  factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
262  + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
263 
264  // Extra dampening factor to go from BE_3 to BE_32.
265  factor *= 1. - exp(-Q2old * R2Ref2);
266 
267  // Add shifts to sum. (Energy component dummy.)
268  pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
269  hadronBE[i1].pComp += pDiff;
270  hadronBE[i2].pComp -= pDiff;
271 
272 }
273 
274 //==========================================================================
275 
276 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5