PYTHIA8/pythia8170/src/BoseEinstein.cxx

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