PYTHIA8/pythia8170/src/MiniStringFragmentation.cxx

   1 // MiniStringFragmentation.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 .
   7 // MiniStringFragmentation class
   8
   9 #include "MiniStringFragmentation.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The MiniStringFragmentation class.
  16
  17 //--------------------------------------------------------------------------
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20 // These are of technical nature, as described for each.
  21
  22 // Since diffractive by definition is > 1 particle, try hard.
  23 const int MiniStringFragmentation::NTRYDIFFRACTIVE = 200;
  24
  25 // After one-body fragmentation failed, try two-body once more.
  26 const int MiniStringFragmentation::NTRYLASTRESORT  = 100;
  27
  28 // Loop try to combine available endquarks to valid hadron.
  29 const int MiniStringFragmentation::NTRYFLAV        = 10;
  30
  31 //--------------------------------------------------------------------------
  32
  33 // Initialize and save pointers.
  34
  35 void MiniStringFragmentation::init(Info* infoPtrIn, Settings& settings,
  36    ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
  37    StringFlav* flavSelPtrIn, StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {
  38
  39   // Save pointers.
  40   infoPtr         = infoPtrIn;
  41   particleDataPtr = particleDataPtrIn;
  42   rndmPtr         = rndmPtrIn;
  43   flavSelPtr      = flavSelPtrIn;
  44   pTSelPtr        = pTSelPtrIn;
  45   zSelPtr         = zSelPtrIn;
  46
  47   // Initialize the MiniStringFragmentation class proper.
  48   nTryMass        = settings.mode("MiniStringFragmentation:nTry");
  49
  50   // Initialize the b parameter of the z spectrum, used when joining jets.
  51   bLund           = zSelPtr->bAreaLund();
  52
  53 }
  54
  55 //--------------------------------------------------------------------------
  56
  57 // Do the fragmentation: driver routine.
  58
  59 bool MiniStringFragmentation::fragment(int iSub, ColConfig& colConfig,
  60   Event& event, bool isDiff) {
  61
  62   // Read in info on system to be treated.
  63   iParton  = colConfig[iSub].iParton;
  64   flav1    = FlavContainer( event[ iParton.front() ].id() );
  65   flav2    = FlavContainer( event[ iParton.back() ].id() );
  66   pSum     = colConfig[iSub].pSum;
  67   mSum     = colConfig[iSub].mass;
  68   m2Sum    = mSum*mSum;
  69   isClosed = colConfig[iSub].isClosed;
  70
  71   // Do not want diffractive systems to easily collapse to one particle.
  72   int nTryFirst = (isDiff) ? NTRYDIFFRACTIVE : nTryMass;
  73
  74   // First try to produce two particles from the system.
  75   if (ministring2two( nTryFirst, event)) return true;
  76
  77   // If this fails, then form one hadron and shuffle momentum.
  78   if (ministring2one( iSub, colConfig, event)) return true;
  79
  80   // If also this fails, then try harder to produce two particles.
  81   if (ministring2two( NTRYLASTRESORT, event)) return true;
  82
  83   // Else complete failure.
  84   infoPtr->errorMsg("Error in MiniStringFragmentation::fragment: "
  85       "no 1- or 2-body state found above mass threshold");
  86   return false;
  87
  88 }
  89
  90 //--------------------------------------------------------------------------
  91
  92   // Attempt to produce two particles from the ministring.
  93
  94 bool MiniStringFragmentation::ministring2two( int nTry, Event& event) {
  95
  96   // Properties of the produced hadrons.
  97   int    idHad1  = 0;
  98   int    idHad2  = 0;
  99   double mHad1   = 0.;
 100   double mHad2   = 0.;
 101   double mHadSum = 0.;
 102
 103   // Allow a few attempts to find a particle pair with low enough masses.
 104   for (int iTry = 0; iTry < nTry; ++iTry) {
 105
 106     // For closed gluon loop need to pick an initial flavour.
 107     if (isClosed) do {
 108       int idStart = flavSelPtr->pickLightQ();
 109       FlavContainer flavStart(idStart, 1);
 110       flavStart = flavSelPtr->pick( flavStart);
 111       flav1 = flavSelPtr->pick( flavStart);
 112       flav2.anti(flav1);
 113     } while (flav1.id == 0 || flav1.nPop > 0);
 114
 115     // Create a new q qbar flavour to form two hadrons.
 116     // Start from a diquark, if any.
 117     do {
 118       FlavContainer flav3 =
 119         (flav1.isDiquark() || (!flav2.isDiquark() && rndmPtr->flat() < 0.5) )
 120         ? flavSelPtr->pick( flav1) : flavSelPtr->pick( flav2).anti();
 121       idHad1 = flavSelPtr->combine( flav1, flav3);
 122       idHad2 = flavSelPtr->combine( flav2, flav3.anti());
 123     } while (idHad1 == 0 || idHad2 == 0);
 124
 125     // Check whether the mass sum fits inside the available phase space.
 126     mHad1 = particleDataPtr->mass(idHad1);
 127     mHad2 = particleDataPtr->mass(idHad2);
 128     mHadSum = mHad1 + mHad2;
 129     if (mHadSum < mSum) break;
 130   }
 131   if (mHadSum >= mSum) return false;
 132
 133   // Define an effective two-parton string, by splitting intermediate
 134   // gluon momenta in proportion to their closeness to either endpoint.
 135   Vec4 pSum1 = event[ iParton.front() ].p();
 136   Vec4 pSum2 = event[ iParton.back() ].p();
 137   if (iParton.size() > 2) {
 138     Vec4 pEnd1 = pSum1;
 139     Vec4 pEnd2 = pSum2;
 140     Vec4 pEndSum = pEnd1 + pEnd2;
 141     for (int i = 1; i < int(iParton.size()) - 1 ; ++i) {
 142       Vec4 pNow = event[ iParton[i] ].p();
 143       double ratio = (pEnd2 * pNow) / (pEndSum * pNow);
 144       pSum1 += ratio * pNow;
 145       pSum2 += (1. - ratio) * pNow;
 146     }
 147   }
 148
 149   // Set up a string region based on the two effective endpoints.
 150   StringRegion region;
 151   region.setUp( pSum1, pSum2);
 152
 153   // Generate an isotropic decay in the ministring rest frame,
 154   // suppressed at large pT by a fragmentation pT Gaussian.
 155   double pAbs2 = 0.25 * ( pow2(m2Sum - mHad1*mHad1 - mHad2*mHad2)
 156     - pow2(2. * mHad1 * mHad2) ) / m2Sum;
 157   double pT2 = 0.;
 158   do {
 159     double cosTheta = rndmPtr->flat();
 160     pT2 = (1. - pow2(cosTheta)) * pAbs2;
 161   } while (pTSelPtr->suppressPT2(pT2) < rndmPtr->flat() );
 162
 163   // Construct the forward-backward asymmetry of the two particles.
 164   double mT21 = mHad1*mHad1 + pT2;
 165   double mT22 = mHad2*mHad2 + pT2;
 166   double lambda = sqrtpos( pow2(m2Sum  - mT21 - mT22) - 4. * mT21 * mT22 );
 167   double probReverse = 1. / (1. + exp( min( 50., bLund * lambda) ) );
 168
 169   // Construct kinematics, as viewed in the transverse rest frame.
 170   double xpz1 = 0.5 * lambda/ m2Sum;
 171   if (probReverse > rndmPtr->flat()) xpz1 = -xpz1;
 172   double xmDiff = (mT21 - mT22) / m2Sum;
 173   double xe1 = 0.5 * (1. + xmDiff);
 174   double xe2 = 0.5 * (1. - xmDiff );
 175
 176   // Distribute pT isotropically in angle.
 177   double phi = 2. * M_PI * rndmPtr->flat();
 178   double pT  = sqrt(pT2);
 179   double px  = pT * cos(phi);
 180   double py  = pT * sin(phi);
 181
 182   // Translate this into kinematics in the string frame.
 183   Vec4 pHad1 = region.pHad( xe1 + xpz1, xe1 - xpz1,  px,  py);
 184   Vec4 pHad2 = region.pHad( xe2 - xpz1, xe2 + xpz1, -px, -py);
 185
 186   // Add produced particles to the event record.
 187   int iFirst = event.append( idHad1, 82, iParton.front(), iParton.back(),
 188     0, 0, 0, 0, pHad1, mHad1);
 189   int iLast = event.append( idHad2, 82, iParton.front(), iParton.back(),
 190     0, 0, 0, 0, pHad2, mHad2);
 191
 192   // Set decay vertex when this is displaced.
 193   if (event[iParton.front()].hasVertex()) {
 194     Vec4 vDec = event[iParton.front()].vDec();
 195     event[iFirst].vProd( vDec );
 196     event[iLast].vProd( vDec );
 197   }
 198
 199   // Set lifetime of hadrons.
 200   event[iFirst].tau( event[iFirst].tau0() * rndmPtr->exp() );
 201   event[iLast].tau( event[iLast].tau0() * rndmPtr->exp() );
 202
 203   // Mark original partons as hadronized and set their daughter range.
 204   for (int i = 0; i < int(iParton.size()); ++i) {
 205     event[ iParton[i] ].statusNeg();
 206     event[ iParton[i] ].daughters(iFirst, iLast);
 207   }
 208
 209   // Successfully done.
 210   return true;
 211
 212 }
 213
 214 //--------------------------------------------------------------------------
 215
 216 // Attempt to produce one particle from a ministring.
 217 // Current algorithm: find the system with largest invariant mass
 218 // relative to the existing one, and boost that system appropriately.
 219 // Try more sophisticated alternatives later?? (Z0 mass shifted??)
 220 // Also, if problems, attempt several times to obtain closer mass match??
 221
 222 bool MiniStringFragmentation::ministring2one( int iSub,
 223   ColConfig& colConfig, Event& event) {
 224
 225   // Cannot handle qq + qbarqbar system.
 226   if (abs(flav1.id) > 100 && abs(flav2.id) > 100) return false;
 227
 228   // For closed gluon loop need to pick an initial flavour.
 229   if (isClosed) do {
 230     int idStart = flavSelPtr->pickLightQ();
 231     FlavContainer flavStart(idStart, 1);
 232     flav1 = flavSelPtr->pick( flavStart);
 233     flav2 = flav1.anti();
 234   } while (abs(flav1.id) > 100);
 235
 236   // Select hadron flavour from available quark flavours.
 237   int idHad = 0;
 238   for (int iTryFlav = 0; iTryFlav < NTRYFLAV; ++iTryFlav) {
 239     idHad = flavSelPtr->combine( flav1, flav2);
 240     if (idHad != 0) break;
 241   }
 242   if (idHad == 0) return false;
 243
 244   // Find mass.
 245   double mHad = particleDataPtr->mass(idHad);
 246
 247   // Find the untreated parton system which combines to the largest
 248   // squared mass above mimimum required.
 249   int iMax = -1;
 250   double deltaM2 = mHad*mHad - mSum*mSum;
 251   double delta2Max = 0.;
 252   for (int iRec = iSub + 1; iRec < colConfig.size(); ++iRec) {
 253     double delta2Rec = 2. * (pSum * colConfig[iRec].pSum) - deltaM2
 254       - 2. * mHad * colConfig[iRec].mass;
 255     if (delta2Rec > delta2Max) { iMax = iRec; delta2Max = delta2Rec;}
 256   }
 257   if (iMax == -1) return false;
 258
 259   // Construct kinematics of the hadron and recoiling system.
 260   Vec4& pRec     = colConfig[iMax].pSum;
 261   double mRec    = colConfig[iMax].mass;
 262   double vecProd = pSum * pRec;
 263   double coefOld = mSum*mSum + vecProd;
 264   double coefNew = mHad*mHad + vecProd;
 265   double coefRec = mRec*mRec + vecProd;
 266   double coefSum = coefOld + coefNew;
 267   double sHat    = coefOld + coefRec;
 268   double root    = sqrtpos( (pow2(coefSum) - 4. * sHat * mHad*mHad)
 269     / (pow2(vecProd) - pow2(mSum * mRec)) );
 270   double k2      = 0.5 * (coefOld * root - coefSum) / sHat;
 271   double k1      = (coefRec * k2 + 0.5 * deltaM2) / coefOld;
 272   Vec4 pHad      = (1. + k1) * pSum - k2 * pRec;
 273   Vec4 pRecNew   = (1. + k2) * pRec - k1 * pSum;
 274
 275   // Add the produced particle to the event record.
 276   int iHad = event.append( idHad, 81, iParton.front(), iParton.back(),
 277     0, 0, 0, 0, pHad, mHad);
 278
 279   // Set decay vertex when this is displaced.
 280   if (event[iParton.front()].hasVertex()) {
 281     Vec4 vDec = event[iParton.front()].vDec();
 282     event[iHad].vProd( vDec );
 283   }
 284
 285   // Set lifetime of hadron.
 286   event[iHad].tau( event[iHad].tau0() * rndmPtr->exp() );
 287
 288   // Mark original partons as hadronized and set their daughter range.
 289   for (int i = 0; i < int(iParton.size()); ++i) {
 290     event[ iParton[i] ].statusNeg();
 291     event[ iParton[i] ].daughters(iHad, iHad);
 292   }
 293
 294   // Copy down recoiling system, with boosted momentum. Update current partons.
 295   RotBstMatrix M;
 296   M.bst(pRec, pRecNew);
 297   for (int i = 0; i < colConfig[iMax].size(); ++i) {
 298     int iOld = colConfig[iMax].iParton[i];
 299     // Do not touch negative iOld = beginning of new junction leg.
 300     if (iOld >= 0) {
 301       int iNew = event.copy(iOld, 72);
 302       event[iNew].rotbst(M);
 303       colConfig[iMax].iParton[i] = iNew;
 304     }
 305   }
 306   colConfig[iMax].pSum = pRecNew;
 307   colConfig[iMax].isCollected = true;
 308
 309   // Successfully done.
 310   return true;
 311
 312 }
 313
 314 //==========================================================================
 315
 316 } // end namespace Pythia8