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