[1] | 1 | // SigmaTotal.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 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 | // MBR parameters. Integration of MBR cross section. |
---|
| 108 | const int SigmaTotal::NINTEG = 1000; |
---|
| 109 | const int SigmaTotal::NINTEG2 = 40; |
---|
| 110 | const double SigmaTotal::HBARC2 = 0.38938; |
---|
| 111 | // MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV^-2. |
---|
| 112 | const double SigmaTotal::FFA1 = 0.9; |
---|
| 113 | const double SigmaTotal::FFA2 = 0.1; |
---|
| 114 | const double SigmaTotal::FFB1 = 4.6; |
---|
| 115 | const double SigmaTotal::FFB2 = 0.6; |
---|
| 116 | |
---|
| 117 | //-------------------------------------------------------------------------- |
---|
| 118 | |
---|
| 119 | // Store pointer to Info and initialize data members. |
---|
| 120 | |
---|
| 121 | void SigmaTotal::init(Info* infoPtrIn, Settings& settings, |
---|
| 122 | ParticleData* particleDataPtrIn) { |
---|
| 123 | |
---|
| 124 | // Store pointers. |
---|
| 125 | infoPtr = infoPtrIn; |
---|
| 126 | particleDataPtr = particleDataPtrIn; |
---|
| 127 | |
---|
| 128 | // Normalization of central diffractive cross section. |
---|
| 129 | zeroAXB = settings.flag("SigmaTotal:zeroAXB"); |
---|
| 130 | sigAXB2TeV = settings.parm("SigmaTotal:sigmaAXB2TeV"); |
---|
| 131 | |
---|
| 132 | // User-set values for cross sections. |
---|
| 133 | setTotal = settings.flag("SigmaTotal:setOwn"); |
---|
| 134 | sigTotOwn = settings.parm("SigmaTotal:sigmaTot"); |
---|
| 135 | sigElOwn = settings.parm("SigmaTotal:sigmaEl"); |
---|
| 136 | sigXBOwn = settings.parm("SigmaTotal:sigmaXB"); |
---|
| 137 | sigAXOwn = settings.parm("SigmaTotal:sigmaAX"); |
---|
| 138 | sigXXOwn = settings.parm("SigmaTotal:sigmaXX"); |
---|
| 139 | sigAXBOwn = settings.parm("SigmaTotal:sigmaAXB"); |
---|
| 140 | |
---|
| 141 | // User-set values to dampen diffractive cross sections. |
---|
| 142 | doDampen = settings.flag("SigmaDiffractive:dampen"); |
---|
| 143 | maxXBOwn = settings.parm("SigmaDiffractive:maxXB"); |
---|
| 144 | maxAXOwn = settings.parm("SigmaDiffractive:maxAX"); |
---|
| 145 | maxXXOwn = settings.parm("SigmaDiffractive:maxXX"); |
---|
| 146 | maxAXBOwn = settings.parm("SigmaDiffractive:maxAXB"); |
---|
| 147 | |
---|
| 148 | // User-set values for handling of elastic sacattering. |
---|
| 149 | setElastic = settings.flag("SigmaElastic:setOwn"); |
---|
| 150 | bSlope = settings.parm("SigmaElastic:bSlope"); |
---|
| 151 | rho = settings.parm("SigmaElastic:rho"); |
---|
| 152 | lambda = settings.parm("SigmaElastic:lambda"); |
---|
| 153 | tAbsMin = settings.parm("SigmaElastic:tAbsMin"); |
---|
| 154 | alphaEM0 = settings.parm("StandardModel:alphaEM0"); |
---|
| 155 | |
---|
| 156 | // Parameters for diffractive systems. |
---|
| 157 | sigmaPomP = settings.parm("Diffraction:sigmaRefPomP"); |
---|
| 158 | mPomP = settings.parm("Diffraction:mRefPomP"); |
---|
| 159 | pPomP = settings.parm("Diffraction:mPowPomP"); |
---|
| 160 | |
---|
| 161 | // Parameters for MBR model. |
---|
| 162 | PomFlux = settings.mode("Diffraction:PomFlux"); |
---|
| 163 | MBReps = settings.parm("Diffraction:MBRepsilon"); |
---|
| 164 | MBRalpha = settings.parm("Diffraction:MBRalpha"); |
---|
| 165 | MBRbeta0 = settings.parm("Diffraction:MBRbeta0"); |
---|
| 166 | MBRsigma0 = settings.parm("Diffraction:MBRsigma0"); |
---|
| 167 | m2min = settings.parm("Diffraction:MBRm2Min"); |
---|
| 168 | dyminSDflux = settings.parm("Diffraction:MBRdyminSDflux"); |
---|
| 169 | dyminDDflux = settings.parm("Diffraction:MBRdyminDDflux"); |
---|
| 170 | dyminCDflux = settings.parm("Diffraction:MBRdyminCDflux"); |
---|
| 171 | dyminSD = settings.parm("Diffraction:MBRdyminSD"); |
---|
| 172 | dyminDD = settings.parm("Diffraction:MBRdyminDD"); |
---|
| 173 | dyminCD = settings.parm("Diffraction:MBRdyminCD"); |
---|
| 174 | dyminSigSD = settings.parm("Diffraction:MBRdyminSigSD"); |
---|
| 175 | dyminSigDD = settings.parm("Diffraction:MBRdyminSigDD"); |
---|
| 176 | dyminSigCD = settings.parm("Diffraction:MBRdyminSigCD"); |
---|
| 177 | |
---|
| 178 | } |
---|
| 179 | |
---|
| 180 | //-------------------------------------------------------------------------- |
---|
| 181 | |
---|
| 182 | // Function that calculates the relevant properties. |
---|
| 183 | |
---|
| 184 | bool SigmaTotal::calc( int idA, int idB, double eCM) { |
---|
| 185 | |
---|
| 186 | // Derived quantities. |
---|
| 187 | alP2 = 2. * ALPHAPRIME; |
---|
| 188 | s0 = 1. / ALPHAPRIME; |
---|
| 189 | |
---|
| 190 | // Reset everything to zero to begin with. |
---|
| 191 | isCalc = false; |
---|
| 192 | sigTot = sigEl = sigXB = sigAX = sigXX = sigAXB = sigND = bEl = s |
---|
| 193 | = bA = bB = 0.; |
---|
| 194 | |
---|
| 195 | // Order flavour of incoming hadrons: idAbsA < idAbsB (restore later). |
---|
| 196 | int idAbsA = abs(idA); |
---|
| 197 | int idAbsB = abs(idB); |
---|
| 198 | bool swapped = false; |
---|
| 199 | if (idAbsA > idAbsB) { |
---|
| 200 | swap( idAbsA, idAbsB); |
---|
| 201 | swapped = true; |
---|
| 202 | } |
---|
| 203 | double sameSign = (idA * idB > 0); |
---|
| 204 | |
---|
| 205 | // Find process number. |
---|
| 206 | int iProc = -1; |
---|
| 207 | if (idAbsA > 1000) { |
---|
| 208 | iProc = (sameSign) ? 0 : 1; |
---|
| 209 | } else if (idAbsA > 100 && idAbsB > 1000) { |
---|
| 210 | iProc = (sameSign) ? 2 : 3; |
---|
| 211 | if (idAbsA/10 == 11 || idAbsA/10 == 22) iProc = 4; |
---|
| 212 | if (idAbsA > 300) iProc = 5; |
---|
| 213 | if (idAbsA > 400) iProc = 6; |
---|
| 214 | if (idAbsA > 900) iProc = 13; |
---|
| 215 | } else if (idAbsA > 100) { |
---|
| 216 | iProc = 7; |
---|
| 217 | if (idAbsB > 300) iProc = 8; |
---|
| 218 | if (idAbsB > 400) iProc = 9; |
---|
| 219 | if (idAbsA > 300) iProc = 10; |
---|
| 220 | if (idAbsA > 300 && idAbsB > 400) iProc = 11; |
---|
| 221 | if (idAbsA > 400) iProc = 12; |
---|
| 222 | } |
---|
| 223 | if (iProc == -1) return false; |
---|
| 224 | |
---|
| 225 | // Primitive implementation of Pomeron + p. |
---|
| 226 | if (iProc == 13) { |
---|
| 227 | s = eCM*eCM; |
---|
| 228 | sigTot = sigmaPomP * pow( eCM / mPomP, pPomP); |
---|
| 229 | sigND = sigTot; |
---|
| 230 | isCalc = true; |
---|
| 231 | return true; |
---|
| 232 | } |
---|
| 233 | |
---|
| 234 | // Find hadron masses and check that energy is enough. |
---|
| 235 | // For mesons use the corresponding vector meson masses. |
---|
| 236 | int idModA = (idAbsA > 1000) ? idAbsA : 10 * (idAbsA/10) + 3; |
---|
| 237 | int idModB = (idAbsB > 1000) ? idAbsB : 10 * (idAbsB/10) + 3; |
---|
| 238 | double mA = particleDataPtr->m0(idModA); |
---|
| 239 | double mB = particleDataPtr->m0(idModB); |
---|
| 240 | if (eCM < mA + mB + MMIN) { |
---|
| 241 | infoPtr->errorMsg("Error in SigmaTotal::calc: too low energy"); |
---|
| 242 | return false; |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | // Evaluate the total cross section. |
---|
| 246 | s = eCM*eCM; |
---|
| 247 | double sEps = pow( s, EPSILON); |
---|
| 248 | double sEta = pow( s, ETA); |
---|
| 249 | sigTot = X[iProc] * sEps + Y[iProc] * sEta; |
---|
| 250 | |
---|
| 251 | // Slope of hadron form factors. |
---|
| 252 | int iHadA = IHADATABLE[iProc]; |
---|
| 253 | int iHadB = IHADBTABLE[iProc]; |
---|
| 254 | bA = BHAD[iHadA]; |
---|
| 255 | bB = BHAD[iHadB]; |
---|
| 256 | |
---|
| 257 | // Elastic slope parameter and cross section. |
---|
| 258 | bEl = 2.*bA + 2.*bB + 4.*sEps - 4.2; |
---|
| 259 | sigEl = CONVERTEL * pow2(sigTot) / bEl; |
---|
| 260 | |
---|
| 261 | // Lookup coefficients for single and double diffraction. |
---|
| 262 | int iSD = ISDTABLE[iProc]; |
---|
| 263 | int iDD = IDDTABLE[iProc]; |
---|
| 264 | double sum1, sum2, sum3, sum4; |
---|
| 265 | |
---|
| 266 | // Single diffractive scattering A + B -> X + B cross section. |
---|
| 267 | mMinXBsave = mA + MMIN0; |
---|
| 268 | double sMinXB = pow2(mMinXBsave); |
---|
| 269 | mResXBsave = mA + MRES0; |
---|
| 270 | double sResXB = pow2(mResXBsave); |
---|
| 271 | double sRMavgXB = mResXBsave * mMinXBsave; |
---|
| 272 | double sRMlogXB = log(1. + sResXB/sMinXB); |
---|
| 273 | double sMaxXB = CSD[iSD][0] * s + CSD[iSD][1]; |
---|
| 274 | double BcorrXB = CSD[iSD][2] + CSD[iSD][3] / s; |
---|
| 275 | sum1 = log( (2.*bB + alP2 * log(s/sMinXB)) |
---|
| 276 | / (2.*bB + alP2 * log(s/sMaxXB)) ) / alP2; |
---|
| 277 | sum2 = CRES * sRMlogXB / (2.*bB + alP2 * log(s/sRMavgXB) + BcorrXB) ; |
---|
| 278 | sigXB = CONVERTSD * X[iProc] * BETA0[iHadB] * max( 0., sum1 + sum2); |
---|
| 279 | |
---|
| 280 | // Single diffractive scattering A + B -> A + X cross section. |
---|
| 281 | mMinAXsave = mB + MMIN0; |
---|
| 282 | double sMinAX = pow2(mMinAXsave); |
---|
| 283 | mResAXsave = mB + MRES0; |
---|
| 284 | double sResAX = pow2(mResAXsave); |
---|
| 285 | double sRMavgAX = mResAXsave * mMinAXsave; |
---|
| 286 | double sRMlogAX = log(1. + sResAX/sMinAX); |
---|
| 287 | double sMaxAX = CSD[iSD][4] * s + CSD[iSD][5]; |
---|
| 288 | double BcorrAX = CSD[iSD][6] + CSD[iSD][7] / s; |
---|
| 289 | sum1 = log( (2.*bA + alP2 * log(s/sMinAX)) |
---|
| 290 | / (2.*bA + alP2 * log(s/sMaxAX)) ) / alP2; |
---|
| 291 | sum2 = CRES * sRMlogAX / (2.*bA + alP2 * log(s/sRMavgAX) + BcorrAX) ; |
---|
| 292 | sigAX = CONVERTSD * X[iProc] * BETA0[iHadA] * max( 0., sum1 + sum2); |
---|
| 293 | |
---|
| 294 | // Order single diffractive correctly. |
---|
| 295 | if (swapped) { |
---|
| 296 | swap( bB, bA); |
---|
| 297 | swap( sigXB, sigAX); |
---|
| 298 | swap( mMinXBsave, mMinAXsave); |
---|
| 299 | swap( mResXBsave, mResAXsave); |
---|
| 300 | } |
---|
| 301 | |
---|
| 302 | // Double diffractive scattering A + B -> X1 + X2 cross section. |
---|
| 303 | double y0min = log( s * SPROTON / (sMinXB * sMinAX) ) ; |
---|
| 304 | double sLog = log(s); |
---|
| 305 | double Delta0 = CDD[iDD][0] + CDD[iDD][1] / sLog |
---|
| 306 | + CDD[iDD][2] / pow2(sLog); |
---|
| 307 | sum1 = (y0min * (log( max( 1e-10, y0min/Delta0) ) - 1.) + Delta0)/ alP2; |
---|
| 308 | if (y0min < 0.) sum1 = 0.; |
---|
| 309 | double sMaxXX = s * ( CDD[iDD][3] + CDD[iDD][4] / sLog |
---|
| 310 | + CDD[iDD][5] / pow2(sLog) ); |
---|
| 311 | double sLogUp = log( max( 1.1, s * s0 / (sMinXB * sRMavgAX) )); |
---|
| 312 | double sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgAX) )); |
---|
| 313 | sum2 = CRES * log( sLogUp / sLogDn ) * sRMlogAX / alP2; |
---|
| 314 | sLogUp = log( max( 1.1, s * s0 / (sMinAX * sRMavgXB) )); |
---|
| 315 | sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgXB) )); |
---|
| 316 | sum3 = CRES * log(sLogUp / sLogDn) * sRMlogXB / alP2; |
---|
| 317 | double BcorrXX = CDD[iDD][6] + CDD[iDD][7] / eCM + CDD[iDD][8] / s; |
---|
| 318 | sum4 = pow2(CRES) * sRMlogAX * sRMlogXB |
---|
| 319 | / max( 0.1, alP2 * log( s * s0 / (sRMavgAX * sRMavgXB) ) + BcorrXX); |
---|
| 320 | sigXX = CONVERTDD * X[iProc] * max( 0., sum1 + sum2 + sum3 + sum4); |
---|
| 321 | |
---|
| 322 | // Central diffractive scattering A + B -> A + X + B, only p and pbar. |
---|
| 323 | mMinAXBsave = 1.; |
---|
| 324 | if (idAbsA == 2212 && idAbsB == 2212 && !zeroAXB) { |
---|
| 325 | double sMinAXB = pow2(mMinAXBsave); |
---|
| 326 | double sRefAXB = pow2(2000.); |
---|
| 327 | sigAXB = sigAXB2TeV * pow( log(0.06 * s / sMinAXB), 1.5 ) |
---|
| 328 | / pow( log(0.06 * sRefAXB / sMinAXB), 1.5 ); |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | // Option with user-requested damping of diffractive cross sections. |
---|
| 332 | if (doDampen) { |
---|
| 333 | sigXB = sigXB * maxXBOwn / (sigXB + maxXBOwn); |
---|
| 334 | sigAX = sigAX * maxAXOwn / (sigAX + maxAXOwn); |
---|
| 335 | sigXX = sigXX * maxXXOwn / (sigXX + maxXXOwn); |
---|
| 336 | sigAXB = sigAXB * maxAXBOwn / (sigAXB + maxAXBOwn); |
---|
| 337 | } |
---|
| 338 | |
---|
| 339 | // Calculate cross sections in MBR model. |
---|
| 340 | if (PomFlux == 5) calcMBRxsecs(idA, idB, eCM); |
---|
| 341 | |
---|
| 342 | // Option with user-set values for total and partial cross sections. |
---|
| 343 | // (Is not done earlier since want diffractive slopes anyway.) |
---|
| 344 | double sigNDOwn = sigTotOwn - sigElOwn - sigXBOwn - sigAXOwn - sigXXOwn |
---|
| 345 | - sigAXBOwn; |
---|
| 346 | double sigElMax = sigEl; |
---|
| 347 | if (setTotal && sigNDOwn > 0.) { |
---|
| 348 | sigTot = sigTotOwn; |
---|
| 349 | sigEl = sigElOwn; |
---|
| 350 | sigXB = sigXBOwn; |
---|
| 351 | sigAX = sigAXOwn; |
---|
| 352 | sigXX = sigXXOwn; |
---|
| 353 | sigAXB = sigAXBOwn; |
---|
| 354 | sigElMax = sigEl; |
---|
| 355 | |
---|
| 356 | // Sub-option to set elastic parameters, including Coulomb contribution. |
---|
| 357 | if (setElastic) { |
---|
| 358 | bEl = bSlope; |
---|
| 359 | sigEl = CONVERTEL * pow2(sigTot) * (1. + rho*rho) / bSlope; |
---|
| 360 | sigElMax = 2. * (sigEl * exp(-bSlope * tAbsMin) |
---|
| 361 | + alphaEM0 * alphaEM0 / (4. * CONVERTEL * tAbsMin) ); |
---|
| 362 | } |
---|
| 363 | } |
---|
| 364 | |
---|
| 365 | // Inelastic nondiffractive by unitarity. |
---|
| 366 | sigND = sigTot - sigEl - sigXB - sigAX - sigXX - sigAXB; |
---|
| 367 | if (sigND < 0.) infoPtr->errorMsg("Error in SigmaTotal::init: " |
---|
| 368 | "sigND < 0"); |
---|
| 369 | else if (sigND < 0.4 * sigTot) infoPtr->errorMsg("Warning in " |
---|
| 370 | "SigmaTotal::init: sigND suspiciously low"); |
---|
| 371 | |
---|
| 372 | // Upper estimate of elastic, including Coulomb term, where appropriate. |
---|
| 373 | sigEl = sigElMax; |
---|
| 374 | |
---|
| 375 | // Done. |
---|
| 376 | isCalc = true; |
---|
| 377 | return true; |
---|
| 378 | |
---|
| 379 | } |
---|
| 380 | |
---|
| 381 | //========================================================================== |
---|
| 382 | |
---|
| 383 | // Calculate parameters in the MBR model. |
---|
| 384 | |
---|
| 385 | bool SigmaTotal::calcMBRxsecs( int idA, int idB, double eCM) { |
---|
| 386 | |
---|
| 387 | // Local variables. |
---|
| 388 | double sigtot, sigel, sigsd, sigdd, sigdpe; |
---|
| 389 | |
---|
| 390 | // MBR parameters locally. |
---|
| 391 | double eps = MBReps; |
---|
| 392 | double alph = MBRalpha; |
---|
| 393 | double beta0gev = MBRbeta0; |
---|
| 394 | double beta0mb = beta0gev * sqrt(HBARC2); |
---|
| 395 | double sigma0mb = MBRsigma0; |
---|
| 396 | double sigma0gev = sigma0mb/HBARC2; |
---|
| 397 | double a1 = FFA1; |
---|
| 398 | double a2 = FFA2; |
---|
| 399 | double b1 = FFB1; |
---|
| 400 | double b2 = FFB2; |
---|
| 401 | |
---|
| 402 | // Calculate total and elastic cross sections. |
---|
| 403 | double ratio; |
---|
| 404 | if (eCM <= 1800.0) { |
---|
| 405 | double sign = (idA * idB > 0); |
---|
| 406 | sigtot = 16.79 * pow(s, 0.104) + 60.81 * pow(s, -0.32) |
---|
| 407 | - sign * 31.68 * pow(s, -0.54); |
---|
| 408 | ratio = 0.100 * pow(s, 0.06) + 0.421 * pow(s, -0.52) |
---|
| 409 | + sign * 0.160 * pow(s, -0.6); |
---|
| 410 | } else { |
---|
| 411 | double sigCDF = 80.03; |
---|
| 412 | double sCDF = pow2(1800.); |
---|
| 413 | double sF = pow2(22.); |
---|
| 414 | sigtot = sigCDF + ( pow2( log(s / sF)) - pow2( log(sCDF / sF)) ) |
---|
| 415 | * M_PI / (3.7 / HBARC2); |
---|
| 416 | ratio = 0.066 + 0.0119 * log(s); |
---|
| 417 | } |
---|
| 418 | sigel=sigtot*ratio; |
---|
| 419 | |
---|
| 420 | // Integrate SD, DD and DPE(CD) cross sections. |
---|
| 421 | // Each cross section is obtained from the ratio of two integrals: |
---|
| 422 | // the Regge cross section and the renormalized flux. |
---|
| 423 | double cflux, csig, c1, step, f; |
---|
| 424 | double dymin0 = 0.; |
---|
| 425 | |
---|
| 426 | // Calculate SD cross section. |
---|
| 427 | double dymaxSD = log(s / m2min); |
---|
| 428 | cflux = pow2(beta0gev) / (16. * M_PI); |
---|
| 429 | csig = cflux * sigma0mb; |
---|
| 430 | |
---|
| 431 | // SD flux. |
---|
| 432 | c1 = cflux; |
---|
| 433 | double fluxsd = 0.; |
---|
| 434 | step = (dymaxSD - dyminSDflux) / NINTEG; |
---|
| 435 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 436 | double dy = dyminSDflux + (i + 0.5) * step; |
---|
| 437 | f = exp(2. * eps * dy) * ( (a1 / (b1 + 2. * alph * dy)) |
---|
| 438 | + (a2 / (b2 + 2. * alph * dy)) ); |
---|
| 439 | f *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD)); |
---|
| 440 | fluxsd = fluxsd + step * c1 * f; |
---|
| 441 | } |
---|
| 442 | if (fluxsd < 1.) fluxsd = 1.; |
---|
| 443 | |
---|
| 444 | // Regge cross section. |
---|
| 445 | c1 = csig * pow(s, eps); |
---|
| 446 | sigsd = 0.; |
---|
| 447 | sdpmax = 0.; |
---|
| 448 | step = (dymaxSD - dymin0) / NINTEG; |
---|
| 449 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 450 | double dy = dymin0 + (i + 0.5) * step; |
---|
| 451 | f = exp(eps * dy) * ( (a1 / (b1 + 2. * alph * dy)) |
---|
| 452 | + (a2 / (b2 + 2. * alph * dy)) ); |
---|
| 453 | f *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD)); |
---|
| 454 | if (f > sdpmax) sdpmax = f; |
---|
| 455 | sigsd = sigsd + step * c1 * f; |
---|
| 456 | } |
---|
| 457 | sdpmax *= 1.01; |
---|
| 458 | sigsd /= fluxsd; |
---|
| 459 | |
---|
| 460 | // Calculate DD cross section. |
---|
| 461 | // Note: dymaxDD = ln(s * s0 /mMin^4) with s0 = 1 GeV^2. |
---|
| 462 | double dymaxDD = log(s / pow2(m2min)); |
---|
| 463 | cflux = sigma0gev / (16. * M_PI); |
---|
| 464 | csig = cflux * sigma0mb; |
---|
| 465 | |
---|
| 466 | // DD flux. |
---|
| 467 | c1 = cflux / (2. * alph); |
---|
| 468 | double fluxdd = 0.; |
---|
| 469 | step = (dymaxDD - dyminDDflux) / NINTEG; |
---|
| 470 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 471 | double dy = dyminDDflux + (i + 0.5) * step; |
---|
| 472 | f = (dymaxDD - dy) * exp(2. * eps * dy) |
---|
| 473 | * ( exp(-2. * alph * dy * exp(-dy)) |
---|
| 474 | - exp(-2. * alph * dy * exp(dy)) ) / dy; |
---|
| 475 | f *= 0.5 * (1. + erf( (dy - dyminDD) / dyminSigDD)); |
---|
| 476 | fluxdd = fluxdd + step * c1 * f; |
---|
| 477 | } |
---|
| 478 | if (fluxdd < 1.) fluxdd = 1.; |
---|
| 479 | |
---|
| 480 | // Regge cross section. |
---|
| 481 | c1 = csig * pow(s, eps) / (2. * alph); |
---|
| 482 | ddpmax = 0.; |
---|
| 483 | sigdd = 0.; |
---|
| 484 | step = (dymaxDD - dymin0) / NINTEG; |
---|
| 485 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 486 | double dy = dymin0 + (i + 0.5) * step; |
---|
| 487 | f = (dymaxDD - dy) * exp(eps * dy) |
---|
| 488 | * ( exp(-2. * alph * dy * exp(-dy)) |
---|
| 489 | - exp(-2. * alph * dy * exp(dy)) ) / dy; |
---|
| 490 | f *= 0.5 * (1. + erf( (dy - dyminDD) / dyminSigDD)); |
---|
| 491 | if (f > ddpmax) ddpmax = f; |
---|
| 492 | sigdd = sigdd + step * c1 * f; |
---|
| 493 | } |
---|
| 494 | ddpmax *= 1.01; |
---|
| 495 | sigdd /= fluxdd; |
---|
| 496 | |
---|
| 497 | // Calculate DPE (CD) cross section. |
---|
| 498 | double dymaxCD = log(s / m2min); |
---|
| 499 | cflux = pow4(beta0gev) / pow2(16. * M_PI); |
---|
| 500 | csig = cflux * pow2(sigma0mb / beta0mb); |
---|
| 501 | double dy1, dy2, f1, f2, step2; |
---|
| 502 | |
---|
| 503 | // DPE flux. |
---|
| 504 | c1 = cflux; |
---|
| 505 | double fluxdpe = 0.; |
---|
| 506 | step = (dymaxCD - dyminCDflux) / NINTEG; |
---|
| 507 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 508 | double dy = dyminCDflux + (i + 0.5) * step; |
---|
| 509 | f = 0.; |
---|
| 510 | step2 = (dy - dyminCDflux) / NINTEG2; |
---|
| 511 | for (int j = 0; j < NINTEG2; ++j) { |
---|
| 512 | double yc = -0.5 * (dy - dyminCDflux) + (j + 0.5) * step2; |
---|
| 513 | dy1 = 0.5 * dy - yc; |
---|
| 514 | dy2 = 0.5 * dy + yc; |
---|
| 515 | f1 = exp(2. * eps * dy1) * ( (a1 / (b1 + 2. * alph * dy1)) |
---|
| 516 | + (a2 / (b2 + 2. * alph * dy1)) ); |
---|
| 517 | f2 = exp(2. * eps * dy2) * ( (a1 / (b1 + 2. * alph * dy2)) |
---|
| 518 | + (a2 / (b2 + 2. * alph * dy2)) ); |
---|
| 519 | f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminCD) |
---|
| 520 | / (dyminSigCD / sqrt(2))) ); |
---|
| 521 | f2 *= 0.5 * (1. + erf( (dy2 - 0.5 *dyminCD) |
---|
| 522 | / (dyminSigCD / sqrt(2))) ); |
---|
| 523 | f += f1 * f2 * step2; |
---|
| 524 | } |
---|
| 525 | fluxdpe += step * c1 * f; |
---|
| 526 | } |
---|
| 527 | if (fluxdpe < 1.) fluxdpe = 1.; |
---|
| 528 | |
---|
| 529 | // Regge cross section. |
---|
| 530 | c1 = csig * pow(s, eps); |
---|
| 531 | sigdpe = 0.; |
---|
| 532 | dpepmax = 0; |
---|
| 533 | step = (dymaxCD - dymin0) / NINTEG; |
---|
| 534 | for (int i = 0; i < NINTEG; ++i) { |
---|
| 535 | double dy = dymin0 + (i + 0.5) * step; |
---|
| 536 | f = 0.; |
---|
| 537 | step2 = (dy - dymin0) / NINTEG2; |
---|
| 538 | for (int j = 0; j < NINTEG2; ++j) { |
---|
| 539 | double yc = -0.5 * (dy - dymin0) + (j + 0.5) * step2; |
---|
| 540 | dy1 = 0.5 * dy - yc; |
---|
| 541 | dy2 = 0.5 * dy + yc; |
---|
| 542 | f1 = exp(eps * dy1) * ( (a1 / (b1 + 2. * alph * dy1)) |
---|
| 543 | + (a2 / (b2 + 2. * alph * dy1)) ); |
---|
| 544 | f2 = exp(eps * dy2) * ( (a1 / (b1 + 2. * alph * dy2)) |
---|
| 545 | + (a2 / (b2 + 2. * alph * dy2)) ); |
---|
| 546 | f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminCD) |
---|
| 547 | / (dyminSigCD / sqrt(2))) ); |
---|
| 548 | f2 *= 0.5 * (1. + erf( (dy2 - 0.5 * dyminCD) |
---|
| 549 | /(dyminSigCD / sqrt(2))) ); |
---|
| 550 | f += f1 * f2 * step2; |
---|
| 551 | } |
---|
| 552 | sigdpe += step * c1 * f; |
---|
| 553 | if ( f > dpepmax) dpepmax = f; |
---|
| 554 | } |
---|
| 555 | dpepmax *= 1.01; |
---|
| 556 | sigdpe /= fluxdpe; |
---|
| 557 | |
---|
| 558 | // Diffraction done. Now calculate total inelastic cross section. |
---|
| 559 | sigND = sigtot - (2. * sigsd + sigdd + sigel + sigdpe); |
---|
| 560 | sigTot = sigtot; |
---|
| 561 | sigEl = sigel; |
---|
| 562 | sigAX = sigsd; |
---|
| 563 | sigXB = sigsd; |
---|
| 564 | sigXX = sigdd; |
---|
| 565 | sigAXB = sigdpe; |
---|
| 566 | |
---|
| 567 | return true; |
---|
| 568 | } |
---|
| 569 | |
---|
| 570 | //========================================================================== |
---|
| 571 | |
---|
| 572 | } // end namespace Pythia8 |
---|