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