Line data Source code
1 : // BoseEinstein.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 : // 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(Info* infoPtrIn, Settings& settings,
43 : ParticleData& particleData) {
44 :
45 : // Save pointer.
46 0 : infoPtr = infoPtrIn;
47 :
48 : // Main flags.
49 0 : doPion = settings.flag("BoseEinstein:Pion");
50 0 : doKaon = settings.flag("BoseEinstein:Kaon");
51 0 : doEta = settings.flag("BoseEinstein:Eta");
52 :
53 : // Shape of Bose-Einstein enhancement/suppression.
54 0 : lambda = settings.parm("BoseEinstein:lambda");
55 0 : QRef = settings.parm("BoseEinstein:QRef");
56 :
57 : // Multiples and inverses (= "radii") of distance parameters in Q-space.
58 0 : QRef2 = 2. * QRef;
59 0 : QRef3 = 3. * QRef;
60 0 : R2Ref = 1. / (QRef * QRef);
61 0 : R2Ref2 = 1. / (QRef2 * QRef2);
62 0 : R2Ref3 = 1. / (QRef3 * QRef3);
63 :
64 : // Masses of particles with Bose-Einstein implemented.
65 0 : for (int iSpecies = 0; iSpecies < 9; ++iSpecies)
66 0 : mHadron[iSpecies] = particleData.m0( IDHADRON[iSpecies] );
67 :
68 : // Pair pi, K, eta and eta' masses for use in tables.
69 0 : mPair[0] = 2. * mHadron[0];
70 0 : mPair[1] = 2. * mHadron[3];
71 0 : mPair[2] = 2. * mHadron[7];
72 0 : mPair[3] = 2. * mHadron[8];
73 :
74 : // Loop over the four required tables. Local variables.
75 : double Qnow, Q2now, centerCorr;
76 0 : for (int iTab = 0; iTab < 4; ++iTab) {
77 0 : m2Pair[iTab] = mPair[iTab] * mPair[iTab];
78 :
79 : // Step size and number of steps in normal table.
80 0 : deltaQ[iTab] = STEPSIZE * min(mPair[iTab], QRef);
81 0 : nStep[iTab] = min( 199, 1 + int(3. * QRef / deltaQ[iTab]) );
82 0 : maxQ[iTab] = (nStep[iTab] - 0.1) * deltaQ[iTab];
83 0 : centerCorr = deltaQ[iTab] * deltaQ[iTab] / 12.;
84 :
85 : // Construct normal table recursively in Q space.
86 0 : shift[iTab][0] = 0.;
87 0 : for (int i = 1; i <= nStep[iTab]; ++i) {
88 0 : Qnow = deltaQ[iTab] * (i - 0.5);
89 0 : Q2now = Qnow * Qnow;
90 0 : shift[iTab][i] = shift[iTab][i - 1] + exp(-Q2now * R2Ref)
91 0 : * deltaQ[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
92 : }
93 :
94 : // Step size and number of steps in compensation table.
95 0 : deltaQ3[iTab] = STEPSIZE * min(mPair[iTab], QRef3);
96 0 : nStep3[iTab] = min( 199, 1 + int(9. * QRef / deltaQ3[iTab]) );
97 0 : maxQ3[iTab] = (nStep3[iTab] - 0.1) * deltaQ3[iTab];
98 0 : centerCorr = deltaQ3[iTab] * deltaQ3[iTab] / 12.;
99 :
100 : // Construct compensation table recursively in Q space.
101 0 : shift3[iTab][0] = 0.;
102 0 : for (int i = 1; i <= nStep3[iTab]; ++i) {
103 0 : Qnow = deltaQ3[iTab] * (i - 0.5);
104 0 : Q2now = Qnow * Qnow;
105 0 : shift3[iTab][i] = shift3[iTab][i - 1] + exp(-Q2now * R2Ref3)
106 0 : * deltaQ3[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
107 : }
108 :
109 : }
110 :
111 : // Done.
112 0 : return true;
113 :
114 0 : }
115 :
116 : //--------------------------------------------------------------------------
117 :
118 : // Perform Bose-Einstein corrections on an event.
119 :
120 : bool BoseEinstein::shiftEvent( Event& event) {
121 :
122 : // Reset list of identical particles.
123 0 : hadronBE.resize(0);
124 :
125 : // Loop over all hadron species with BE effects.
126 0 : nStored[0] = 0;
127 0 : for (int iSpecies = 0; iSpecies < 9; ++iSpecies) {
128 0 : nStored[iSpecies + 1] = nStored[iSpecies];
129 0 : if (!doPion && iSpecies <= 2) continue;
130 0 : if (!doKaon && iSpecies >= 3 && iSpecies <= 6) continue;
131 0 : if (!doEta && iSpecies >= 7) continue;
132 :
133 : // Properties of current hadron species.
134 0 : int idNow = IDHADRON[ iSpecies ];
135 0 : int iTab = ITABLE[ iSpecies ];
136 :
137 : // Loop through event record to store copies of current species.
138 0 : for (int i = 0; i < event.size(); ++i)
139 0 : if ( event[i].id() == idNow && event[i].isFinal() )
140 0 : hadronBE.push_back(
141 0 : BoseEinsteinHadron( idNow, i, event[i].p(), event[i].m() ) );
142 0 : nStored[iSpecies + 1] = hadronBE.size();
143 :
144 : // Loop through pairs of identical particles and find shifts.
145 0 : for (int i1 = nStored[iSpecies]; i1 < nStored[iSpecies+1] - 1; ++i1)
146 0 : for (int i2 = i1 + 1; i2 < nStored[iSpecies+1]; ++i2)
147 0 : shiftPair( i1, i2, iTab);
148 0 : }
149 :
150 : // Must have at least two pairs to carry out compensation.
151 0 : if (nStored[9] < 2) return true;
152 :
153 : // Shift momenta and recalculate energies.
154 : double eSumOriginal = 0.;
155 : double eSumShifted = 0.;
156 : double eDiffByComp = 0.;
157 0 : for (int i = 0; i < nStored[9]; ++i) {
158 0 : eSumOriginal += hadronBE[i].p.e();
159 0 : hadronBE[i].p += hadronBE[i].pShift;
160 0 : hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
161 0 : eSumShifted += hadronBE[i].p.e();
162 0 : eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
163 0 : / hadronBE[i].p.e();
164 : }
165 :
166 : // Iterate compensation shift until convergence.
167 : int iStep = 0;
168 0 : while ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal
169 0 : && abs(eSumShifted - eSumOriginal) < COMPFACMAX * abs(eDiffByComp)
170 0 : && iStep < NCOMPSTEP ) {
171 0 : ++iStep;
172 0 : double compFac = (eSumOriginal - eSumShifted) / eDiffByComp;
173 : eSumShifted = 0.;
174 : eDiffByComp = 0.;
175 0 : for (int i = 0; i < nStored[9]; ++i) {
176 0 : hadronBE[i].p += compFac * hadronBE[i].pComp;
177 0 : hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
178 0 : eSumShifted += hadronBE[i].p.e();
179 0 : eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
180 0 : / hadronBE[i].p.e();
181 : }
182 : }
183 :
184 : // Error if no convergence, and then return without doing BE shift.
185 : // However, not grave enough to kill event, so return true.
186 0 : if ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal ) {
187 0 : infoPtr->errorMsg("Warning in BoseEinstein::shiftEvent: "
188 : "no consistent BE shift topology found, so skip BE");
189 0 : return true;
190 : }
191 :
192 : // Store new particle copies with shifted momenta.
193 0 : for (int i = 0; i < nStored[9]; ++i) {
194 0 : int iNew = event.copy( hadronBE[i].iPos, 99);
195 0 : event[ iNew ].p( hadronBE[i].p );
196 : }
197 :
198 : // Done.
199 0 : return true;
200 :
201 0 : }
202 :
203 : //--------------------------------------------------------------------------
204 :
205 : // Calculate shift and (unnormalized) compensation for pair.
206 :
207 : void BoseEinstein::shiftPair( int i1, int i2, int iTab) {
208 :
209 : // Calculate old relative momentum.
210 0 : double Q2old = m2(hadronBE[i1].p, hadronBE[i2].p) - m2Pair[iTab];
211 0 : if (Q2old < Q2MIN) return;
212 0 : double Qold = sqrt(Q2old);
213 0 : double psFac = sqrt(Q2old + m2Pair[iTab]) / Q2old;
214 :
215 : // Calculate new relative momentum for normal shift.
216 : double Qmove = 0.;
217 0 : if (Qold < deltaQ[iTab]) Qmove = Qold / 3.;
218 0 : else if (Qold < maxQ[iTab]) {
219 0 : double realQbin = Qold / deltaQ[iTab];
220 0 : int intQbin = int( realQbin );
221 0 : double inter = (pow3(realQbin) - pow3(intQbin))
222 0 : / (3 * intQbin * (intQbin + 1) + 1);
223 0 : Qmove = ( shift[iTab][intQbin] + inter * (shift[iTab][intQbin + 1]
224 0 : - shift[iTab][intQbin]) ) * psFac;
225 0 : }
226 0 : else Qmove = shift[iTab][nStep[iTab]] * psFac;
227 0 : double Q2new = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove), 2. / 3.);
228 :
229 : // Calculate corresponding three-momentum shift.
230 0 : double Q2Diff = Q2new - Q2old;
231 0 : double p2DiffAbs = (hadronBE[i1].p - hadronBE[i2].p).pAbs2();
232 0 : double p2AbsDiff = hadronBE[i1].p.pAbs2() - hadronBE[i2].p.pAbs2();
233 0 : double eSum = hadronBE[i1].p.e() + hadronBE[i2].p.e();
234 0 : double eDiff = hadronBE[i1].p.e() - hadronBE[i2].p.e();
235 0 : double sumQ2E = Q2Diff + eSum * eSum;
236 0 : double rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
237 0 : double rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
238 0 : double factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
239 0 : + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
240 :
241 : // Add shifts to sum. (Energy component dummy.)
242 0 : Vec4 pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
243 0 : hadronBE[i1].pShift += pDiff;
244 0 : hadronBE[i2].pShift -= pDiff;
245 :
246 : // Calculate new relative momentum for compensation shift.
247 : double Qmove3 = 0.;
248 0 : if (Qold < deltaQ3[iTab]) Qmove3 = Qold / 3.;
249 0 : else if (Qold < maxQ3[iTab]) {
250 0 : double realQbin = Qold / deltaQ3[iTab];
251 0 : int intQbin = int( realQbin );
252 0 : double inter = (pow3(realQbin) - pow3(intQbin))
253 0 : / (3 * intQbin * (intQbin + 1) + 1);
254 0 : Qmove3 = ( shift3[iTab][intQbin] + inter * (shift3[iTab][intQbin + 1]
255 0 : - shift3[iTab][intQbin]) ) * psFac;
256 0 : }
257 0 : else Qmove3 = shift3[iTab][nStep3[iTab]] *psFac;
258 0 : double Q2new3 = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove3), 2. / 3.);
259 :
260 : // Calculate corresponding three-momentum shift.
261 0 : Q2Diff = Q2new3 - Q2old;
262 0 : sumQ2E = Q2Diff + eSum * eSum;
263 0 : rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
264 0 : rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
265 0 : factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
266 0 : + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
267 :
268 : // Extra dampening factor to go from BE_3 to BE_32.
269 0 : factor *= 1. - exp(-Q2old * R2Ref2);
270 :
271 : // Add shifts to sum. (Energy component dummy.)
272 0 : pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
273 0 : hadronBE[i1].pComp += pDiff;
274 0 : hadronBE[i2].pComp -= pDiff;
275 :
276 0 : }
277 :
278 : //==========================================================================
279 :
280 : } // end namespace Pythia8
|
