PYTHIA8/pythia8140/src/SigmaTotal.cxx

   1 // SigmaTotal.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 SigmaTotal class.
   7
   8 #include "SigmaTotal.h"
   9
  10 namespace Pythia8 {
  11
  12 //==========================================================================
  13
  14 // The SigmaTotal class.
  15
  16 // Formulae are taken from:
  17 // G.A. Schuler and T. Sjostrand, Phys. Rev. D49 (1994) 2257,
  18 //   Z. Phys. C73 (1997) 677
  19 // which borrows some total cross sections from
  20 // A. Donnachie and P.V. Landshoff, Phys. Lett. B296 (1992) 227.
  21
  22 // Implemented processes with their process number iProc:
  23 // =  0 : p + p;     =  1 : pbar + p;
  24 // =  2 : pi+ + p;   =  3 : pi- + p;     =  4 : pi0/rho0 + p;
  25 // =  5 : phi + p;   =  6 : J/psi + p;
  26 // =  7 : rho + rho; =  8 : rho + phi;   =  9 : rho + J/psi;
  27 // = 10 : phi + phi; = 11 : phi + J/psi; = 12 : J/psi + J/psi.
  28 // = 13 : Pom + p (preliminary).
  29
  30 //--------------------------------------------------------------------------
  31
  32 // Definitions of static variables.
  33 // Note that a lot of parameters are hardcoded as const here, rather
  34 // than being interfaced for public change, since any changes would
  35 // have to be done in a globally consistent manner. Which basically
  36 // means a rewrite/replacement of the whole class.
  37
  38 // Minimum threshold below which no cross sections will be defined.
  39 const double SigmaTotal::MMIN  = 2.;
  40
  41 // General constants in total cross section parametrization:
  42 // sigmaTot = X * s^epsilon + Y * s^eta (pomeron + reggeon).
  43 const double SigmaTotal::EPSILON = 0.0808;
  44 const double SigmaTotal::ETA     = -0.4525;
  45 const double SigmaTotal::X[] = { 21.70, 21.70, 13.63, 13.63, 13.63,
  46   10.01, 0.970, 8.56, 6.29, 0.609, 4.62, 0.447, 0.0434};
  47 const double SigmaTotal::Y[] = { 56.08, 98.39, 27.56, 36.02, 31.79,
  48   1.51, -0.146, 13.08, -0.62, -0.060, 0.030, -0.0028, 0.00028};
  49
  50 // Type of the two incoming hadrons as function of the process number:
  51 // = 0 : p/n ; = 1 : pi/rho/omega; = 2 : phi; = 3 : J/psi.
  52 const int SigmaTotal::IHADATABLE[] = { 0, 0, 1, 1, 1, 2, 3, 1, 1,
  53   1, 2, 2, 3};
  54 const int SigmaTotal::IHADBTABLE[] = { 0, 0, 0, 0, 0, 0, 0, 1, 2,
  55   3, 2, 3, 3};
  56
  57 // Hadron-Pomeron coupling beta(t) = beta(0) * exp(b*t).
  58 const double SigmaTotal::BETA0[] = { 4.658, 2.926, 2.149, 0.208};
  59 const double SigmaTotal::BHAD[]  = {   2.3,   1.4,   1.4,  0.23};
  60
  61 // Pomeron trajectory alpha(t) = 1 + epsilon + alpha' * t
  62 const double SigmaTotal::ALPHAPRIME = 0.25;
  63
  64 // Conversion coefficients = 1/(16pi) * (mb <-> GeV^2) * (G_3P)^n,
  65 // with n = 0 elastic, n = 1 single and n = 2 double diffractive.
  66 const double SigmaTotal::CONVERTEL = 0.0510925;
  67 const double SigmaTotal::CONVERTSD = 0.0336;
  68 const double SigmaTotal::CONVERTDD = 0.0084;
  69
  70 // Diffractive mass spectrum starts at m + MMIN0 and has a low-mass
  71 // enhancement, factor cRes, up to around m + mRes0.
  72 const double SigmaTotal::MMIN0 = 0.28;
  73 const double SigmaTotal::CRES  = 2.0;
  74 const double SigmaTotal::MRES0 = 1.062;
  75
  76 // Parameters and coefficients for single diffractive scattering.
  77 const int SigmaTotal::ISDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
  78   6, 7, 8, 9};
  79 const double SigmaTotal::CSD[10][8] = {
  80   { 0.213, 0.0, -0.47, 150., 0.213, 0.0, -0.47, 150., } ,
  81   { 0.213, 0.0, -0.47, 150., 0.267, 0.0, -0.47, 100., } ,
  82   { 0.213, 0.0, -0.47, 150., 0.232, 0.0, -0.47, 110., } ,
  83   { 0.213, 7.0, -0.55, 800., 0.115, 0.0, -0.47, 110., } ,
  84   { 0.267, 0.0, -0.46,  75., 0.267, 0.0, -0.46,  75., } ,
  85   { 0.232, 0.0, -0.46,  85., 0.267, 0.0, -0.48, 100., } ,
  86   { 0.115, 0.0, -0.50,  90., 0.267, 6.0, -0.56, 420., } ,
  87   { 0.232, 0.0, -0.48, 110., 0.232, 0.0, -0.48, 110., } ,
  88   { 0.115, 0.0, -0.52, 120., 0.232, 6.0, -0.56, 470., } ,
  89   { 0.115, 5.5, -0.58, 570., 0.115, 5.5, -0.58, 570.  } };
  90
  91 // Parameters and coefficients for double diffractive scattering.
  92 const int SigmaTotal::IDDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
  93   6, 7, 8, 9};
  94 const double SigmaTotal::CDD[10][9] = {
  95   { 3.11, -7.34,  9.71, 0.068, -0.42, 1.31, -1.37,  35.0,  118., } ,
  96   { 3.11, -7.10,  10.6, 0.073, -0.41, 1.17, -1.41,  31.6,   95., } ,
  97   { 3.12, -7.43,  9.21, 0.067, -0.44, 1.41, -1.35,  36.5,  132., } ,
  98   { 3.13, -8.18, -4.20, 0.056, -0.71, 3.12, -1.12,  55.2, 1298., } ,
  99   { 3.11, -6.90,  11.4, 0.078, -0.40, 1.05, -1.40,  28.4,   78., } ,
 100   { 3.11, -7.13,  10.0, 0.071, -0.41, 1.23, -1.34,  33.1,  105., } ,
 101   { 3.12, -7.90, -1.49, 0.054, -0.64, 2.72, -1.13,  53.1,  995., } ,
 102   { 3.11, -7.39,  8.22, 0.065, -0.44, 1.45, -1.36,  38.1,  148., } ,
 103   { 3.18, -8.95, -3.37, 0.057, -0.76, 3.32, -1.12,  55.6, 1472., } ,
 104   { 4.18, -29.2,  56.2, 0.074, -1.36, 6.67, -1.14, 116.2, 6532.  } };
 105 const double SigmaTotal::SPROTON = 0.880;
 106
 107 //--------------------------------------------------------------------------
 108
 109 // Store pointer to Info and initialize data members.
 110
 111 void SigmaTotal::init(Info* infoPtrIn, Settings& settings,
 112   ParticleData* particleDataPtrIn) {
 113
 114   // Store pointers.
 115   infoPtr         = infoPtrIn;
 116   particleDataPtr = particleDataPtrIn;
 117
 118   // User-set values for cross sections.
 119   setTotal   = settings.flag("SigmaTotal:setOwn");
 120   sigTotOwn  = settings.parm("SigmaTotal:sigmaTot");
 121   sigElOwn   = settings.parm("SigmaTotal:sigmaEl");
 122   sigXBOwn   = settings.parm("SigmaTotal:sigmaXB");
 123   sigAXOwn   = settings.parm("SigmaTotal:sigmaAX");
 124   sigXXOwn   = settings.parm("SigmaTotal:sigmaXX");
 125
 126   // User-set values to dampen diffractive cross sections.
 127   doDampen   = settings.flag("SigmaDiffractive:dampen");
 128   maxXBOwn   = settings.parm("SigmaDiffractive:maxXB");
 129   maxAXOwn   = settings.parm("SigmaDiffractive:maxAX");
 130   maxXXOwn   = settings.parm("SigmaDiffractive:maxXX");
 131
 132   // User-set values for handling of elastic sacattering.
 133   setElastic = settings.flag("SigmaElastic:setOwn");
 134   bSlope     = settings.parm("SigmaElastic:bSlope");
 135   rho        = settings.parm("SigmaElastic:rho");
 136   lambda     = settings.parm("SigmaElastic:lambda");
 137   tAbsMin    = settings.parm("SigmaElastic:tAbsMin");
 138   alphaEM0   = settings.parm("StandardModel:alphaEM0");
 139
 140   // Parameter for diffractive systems.
 141   sigmaPomP  = settings.parm("Diffraction:sigmaPomP");
 142
 143 }
 144
 145 //--------------------------------------------------------------------------
 146
 147 // Function that calculates the relevant properties.
 148
 149 bool SigmaTotal::calc( int idA, int idB, double eCM) {
 150
 151   // Derived quantities.
 152   alP2 = 2. * ALPHAPRIME;
 153   s0   = 1. / ALPHAPRIME;
 154
 155   // Reset everything to zero to begin with.
 156   isCalc = false;
 157   sigTot = sigEl = sigXB = sigAX = sigXX = sigND = bEl = s = bA = bB = 0.;
 158
 159   // Order flavour of incoming hadrons: idAbsA < idAbsB (restore later).
 160   int idAbsA = abs(idA);
 161   int idAbsB = abs(idB);
 162   bool swapped = false;
 163   if (idAbsA > idAbsB) {
 164     swap( idAbsA, idAbsB);
 165     swapped = true;
 166   }
 167   double sameSign = (idA * idB > 0);
 168
 169   // Find process number.
 170   int iProc                                       = -1;
 171   if (idAbsA > 1000) {
 172     iProc                                         = (sameSign) ? 0 : 1;
 173   } else if (idAbsA > 100 && idAbsB > 1000) {
 174     iProc                                         = (sameSign) ? 2 : 3;
 175     if (idAbsA/10 == 11 || idAbsA/10 == 22) iProc = 4;
 176     if (idAbsA > 300) iProc                       = 5;
 177     if (idAbsA > 400) iProc                       = 6;
 178     if (idAbsA > 900) iProc                       = 13;
 179   } else if (idAbsA > 100) {
 180     iProc                                         = 7;
 181     if (idAbsB > 300) iProc                       = 8;
 182     if (idAbsB > 400) iProc                       = 9;
 183     if (idAbsA > 300) iProc                       = 10;
 184     if (idAbsA > 300 && idAbsB > 400) iProc       = 11;
 185     if (idAbsA > 400) iProc                       = 12;
 186   }
 187   if (iProc == -1) return false;
 188
 189   // Primitive implementation of Pomeron + p.
 190   if (iProc == 13) {
 191     s      = eCM*eCM;
 192     sigTot = sigmaPomP;
 193     sigND  = sigTot;
 194     isCalc = true;
 195     return true;
 196   }
 197
 198   // Find hadron masses and check that energy is enough.
 199   // For mesons use the corresponding vector meson masses.
 200   int idModA = (idAbsA > 1000) ? idAbsA : 10 * (idAbsA/10) + 3;
 201   int idModB = (idAbsB > 1000) ? idAbsB : 10 * (idAbsB/10) + 3;
 202   double mA  = particleDataPtr->m0(idModA);
 203   double mB  = particleDataPtr->m0(idModB);
 204   if (eCM < mA + mB + MMIN) {
 205     infoPtr->errorMsg("Error in SigmaTotal::calc: too low energy");
 206     return false;
 207   }
 208
 209   // Evaluate the total cross section.
 210   s           = eCM*eCM;
 211   double sEps = pow( s, EPSILON);
 212   double sEta = pow( s, ETA);
 213   sigTot      = X[iProc] * sEps + Y[iProc] * sEta;
 214
 215   // Slope of hadron form factors.
 216   int iHadA = IHADATABLE[iProc];
 217   int iHadB = IHADBTABLE[iProc];
 218   bA        = BHAD[iHadA];
 219   bB        = BHAD[iHadB];
 220
 221   // Elastic slope parameter and cross section.
 222   bEl   = 2.*bA + 2.*bB + 4.*sEps - 4.2;
 223   sigEl = CONVERTEL * pow2(sigTot) / bEl;
 224
 225   // Lookup coefficients for single and double diffraction.
 226   int iSD = ISDTABLE[iProc];
 227   int iDD = IDDTABLE[iProc];
 228   double sum1, sum2, sum3, sum4;
 229
 230   // Single diffractive scattering A + B -> X + B cross section.
 231   mMinXBsave      = mA + MMIN0;
 232   double sMinXB   = pow2(mMinXBsave);
 233   mResXBsave      = mA + MRES0;
 234   double sResXB   = pow2(mResXBsave);
 235   double sRMavgXB = mResXBsave * mMinXBsave;
 236   double sRMlogXB = log(1. + sResXB/sMinXB);
 237   double sMaxXB   = CSD[iSD][0] * s + CSD[iSD][1];
 238   double BcorrXB  = CSD[iSD][2] + CSD[iSD][3] / s;
 239   sum1  = log( (2.*bB + alP2 * log(s/sMinXB))
 240     / (2.*bB + alP2 * log(s/sMaxXB)) ) / alP2;
 241   sum2  = CRES * sRMlogXB / (2.*bB + alP2 * log(s/sRMavgXB) + BcorrXB) ;
 242   sigXB = CONVERTSD * X[iProc] * BETA0[iHadB] * max( 0., sum1 + sum2);
 243
 244   // Single diffractive scattering A + B -> A + X cross section.
 245   mMinAXsave      = mB + MMIN0;
 246   double sMinAX   = pow2(mMinAXsave);
 247   mResAXsave      = mB + MRES0;
 248   double sResAX   = pow2(mResAXsave);
 249   double sRMavgAX = mResAXsave * mMinAXsave;
 250   double sRMlogAX = log(1. + sResAX/sMinAX);
 251   double sMaxAX   = CSD[iSD][4] * s + CSD[iSD][5];
 252   double BcorrAX  = CSD[iSD][6] + CSD[iSD][7] / s;
 253   sum1  = log( (2.*bA + alP2 * log(s/sMinAX))
 254     / (2.*bA + alP2 * log(s/sMaxAX)) ) / alP2;
 255   sum2  = CRES * sRMlogAX / (2.*bA + alP2 * log(s/sRMavgAX) + BcorrAX) ;
 256   sigAX = CONVERTSD * X[iProc] * BETA0[iHadA] * max( 0., sum1 + sum2);
 257
 258   // Order single diffractive correctly.
 259   if (swapped) {
 260     swap( bB, bA);
 261     swap( sigXB, sigAX);
 262     swap( mMinXBsave, mMinAXsave);
 263     swap( mResXBsave, mResAXsave);
 264    }
 265
 266   // Double diffractive scattering A + B -> X1 + X2 cross section.
 267   double y0min = log( s * SPROTON / (sMinXB * sMinAX) ) ;
 268   double sLog  = log(s);
 269   double Delta0 = CDD[iDD][0] + CDD[iDD][1] / sLog
 270     + CDD[iDD][2] / pow2(sLog);
 271   sum1 = (y0min * (log( max( 1e-10, y0min/Delta0) ) - 1.) + Delta0)/ alP2;
 272   if (y0min < 0.) sum1 = 0.;
 273   double sMaxXX = s * ( CDD[iDD][3] + CDD[iDD][4] / sLog
 274     + CDD[iDD][5] / pow2(sLog) );
 275   double sLogUp = log( max( 1.1, s * s0 / (sMinXB * sRMavgAX) ));
 276   double sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgAX) ));
 277   sum2   = CRES * log( sLogUp / sLogDn ) * sRMlogAX / alP2;
 278   sLogUp = log( max( 1.1, s * s0 / (sMinAX * sRMavgXB) ));
 279   sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgXB) ));
 280   sum3   = CRES * log(sLogUp / sLogDn) * sRMlogXB / alP2;
 281   double BcorrXX =  CDD[iDD][6] + CDD[iDD][7] / eCM + CDD[iDD][8] / s;
 282   sum4   = pow2(CRES) * sRMlogAX * sRMlogXB
 283     / max( 0.1, alP2 * log( s * s0 / (sRMavgAX * sRMavgXB) ) + BcorrXX);
 284   sigXX  = CONVERTDD * X[iProc] * max( 0., sum1 + sum2 + sum3 + sum4);
 285
 286   // Option with user-requested damping of diffractive cross sections.
 287   if (doDampen) {
 288     sigXB = sigXB * maxXBOwn / (sigXB + maxXBOwn);
 289     sigAX = sigAX * maxAXOwn / (sigAX + maxAXOwn);
 290     sigXX = sigXX * maxXXOwn / (sigXX + maxXXOwn);
 291   }
 292
 293   // Option with user-set values for total and partial cross sections.
 294   // (Is not done earlier since want diffractive slopes anyway.)
 295   double sigNDOwn = sigTotOwn - sigElOwn - sigXBOwn - sigAXOwn - sigXXOwn;
 296   double sigElMax = sigEl;
 297   if (setTotal && sigNDOwn > 0.) {
 298     sigTot   = sigTotOwn;
 299     sigEl    = sigElOwn;
 300     sigXB    = sigXBOwn;
 301     sigAX    = sigAXOwn;
 302     sigXX    = sigXXOwn;
 303     sigElMax = sigEl;
 304
 305     // Sub-option to set elastic parameters, including Coulomb contribution.
 306     if (setElastic) {
 307       bEl      = bSlope;
 308       sigEl    = CONVERTEL * pow2(sigTot) * (1. + rho*rho) / bSlope;
 309       sigElMax = 2. * (sigEl * exp(-bSlope * tAbsMin)
 310                + alphaEM0 * alphaEM0 / (4. * CONVERTEL * tAbsMin) );
 311     }
 312   }
 313
 314   // Inelastic nondiffractive by unitarity.
 315   sigND = sigTot - sigEl - sigXB - sigAX - sigXX;
 316   if (sigND < 0.) infoPtr->errorMsg("Error in SigmaTotal::init: "
 317     "sigND < 0");
 318   else if (sigND < 0.4 * sigTot) infoPtr->errorMsg("Warning in "
 319     "SigmaTotal::init: sigND suspiciously low");
 320
 321   // Upper estimate of elastic, including Coulomb term, where appropriate.
 322   sigEl = sigElMax;
 323
 324   // Done.
 325   isCalc = true;
 326   return true;
 327
 328 }
 329
 330 //==========================================================================
 331
 332 } // end namespace Pythia8