StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaTotal.cc
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 //--------------------------------------------------------------------------
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  // Parameters for diffractive systems.
141  sigmaPomP = settings.parm("Diffraction:sigmaRefPomP");
142  mPomP = settings.parm("Diffraction:mRefPomP");
143  pPomP = settings.parm("Diffraction:mPowPomP");
144 
145 }
146 
147 //--------------------------------------------------------------------------
148 
149 // Function that calculates the relevant properties.
150 
151 bool SigmaTotal::calc( int idA, int idB, double eCM) {
152 
153  // Derived quantities.
154  alP2 = 2. * ALPHAPRIME;
155  s0 = 1. / ALPHAPRIME;
156 
157  // Reset everything to zero to begin with.
158  isCalc = false;
159  sigTot = sigEl = sigXB = sigAX = sigXX = sigND = bEl = s = bA = bB = 0.;
160 
161  // Order flavour of incoming hadrons: idAbsA < idAbsB (restore later).
162  int idAbsA = abs(idA);
163  int idAbsB = abs(idB);
164  bool swapped = false;
165  if (idAbsA > idAbsB) {
166  swap( idAbsA, idAbsB);
167  swapped = true;
168  }
169  double sameSign = (idA * idB > 0);
170 
171  // Find process number.
172  int iProc = -1;
173  if (idAbsA > 1000) {
174  iProc = (sameSign) ? 0 : 1;
175  } else if (idAbsA > 100 && idAbsB > 1000) {
176  iProc = (sameSign) ? 2 : 3;
177  if (idAbsA/10 == 11 || idAbsA/10 == 22) iProc = 4;
178  if (idAbsA > 300) iProc = 5;
179  if (idAbsA > 400) iProc = 6;
180  if (idAbsA > 900) iProc = 13;
181  } else if (idAbsA > 100) {
182  iProc = 7;
183  if (idAbsB > 300) iProc = 8;
184  if (idAbsB > 400) iProc = 9;
185  if (idAbsA > 300) iProc = 10;
186  if (idAbsA > 300 && idAbsB > 400) iProc = 11;
187  if (idAbsA > 400) iProc = 12;
188  }
189  if (iProc == -1) return false;
190 
191  // Primitive implementation of Pomeron + p.
192  if (iProc == 13) {
193  s = eCM*eCM;
194  sigTot = sigmaPomP * pow( eCM / mPomP, pPomP);
195  sigND = sigTot;
196  isCalc = true;
197  return true;
198  }
199 
200  // Find hadron masses and check that energy is enough.
201  // For mesons use the corresponding vector meson masses.
202  int idModA = (idAbsA > 1000) ? idAbsA : 10 * (idAbsA/10) + 3;
203  int idModB = (idAbsB > 1000) ? idAbsB : 10 * (idAbsB/10) + 3;
204  double mA = particleDataPtr->m0(idModA);
205  double mB = particleDataPtr->m0(idModB);
206  if (eCM < mA + mB + MMIN) {
207  infoPtr->errorMsg("Error in SigmaTotal::calc: too low energy");
208  return false;
209  }
210 
211  // Evaluate the total cross section.
212  s = eCM*eCM;
213  double sEps = pow( s, EPSILON);
214  double sEta = pow( s, ETA);
215  sigTot = X[iProc] * sEps + Y[iProc] * sEta;
216 
217  // Slope of hadron form factors.
218  int iHadA = IHADATABLE[iProc];
219  int iHadB = IHADBTABLE[iProc];
220  bA = BHAD[iHadA];
221  bB = BHAD[iHadB];
222 
223  // Elastic slope parameter and cross section.
224  bEl = 2.*bA + 2.*bB + 4.*sEps - 4.2;
225  sigEl = CONVERTEL * pow2(sigTot) / bEl;
226 
227  // Lookup coefficients for single and double diffraction.
228  int iSD = ISDTABLE[iProc];
229  int iDD = IDDTABLE[iProc];
230  double sum1, sum2, sum3, sum4;
231 
232  // Single diffractive scattering A + B -> X + B cross section.
233  mMinXBsave = mA + MMIN0;
234  double sMinXB = pow2(mMinXBsave);
235  mResXBsave = mA + MRES0;
236  double sResXB = pow2(mResXBsave);
237  double sRMavgXB = mResXBsave * mMinXBsave;
238  double sRMlogXB = log(1. + sResXB/sMinXB);
239  double sMaxXB = CSD[iSD][0] * s + CSD[iSD][1];
240  double BcorrXB = CSD[iSD][2] + CSD[iSD][3] / s;
241  sum1 = log( (2.*bB + alP2 * log(s/sMinXB))
242  / (2.*bB + alP2 * log(s/sMaxXB)) ) / alP2;
243  sum2 = CRES * sRMlogXB / (2.*bB + alP2 * log(s/sRMavgXB) + BcorrXB) ;
244  sigXB = CONVERTSD * X[iProc] * BETA0[iHadB] * max( 0., sum1 + sum2);
245 
246  // Single diffractive scattering A + B -> A + X cross section.
247  mMinAXsave = mB + MMIN0;
248  double sMinAX = pow2(mMinAXsave);
249  mResAXsave = mB + MRES0;
250  double sResAX = pow2(mResAXsave);
251  double sRMavgAX = mResAXsave * mMinAXsave;
252  double sRMlogAX = log(1. + sResAX/sMinAX);
253  double sMaxAX = CSD[iSD][4] * s + CSD[iSD][5];
254  double BcorrAX = CSD[iSD][6] + CSD[iSD][7] / s;
255  sum1 = log( (2.*bA + alP2 * log(s/sMinAX))
256  / (2.*bA + alP2 * log(s/sMaxAX)) ) / alP2;
257  sum2 = CRES * sRMlogAX / (2.*bA + alP2 * log(s/sRMavgAX) + BcorrAX) ;
258  sigAX = CONVERTSD * X[iProc] * BETA0[iHadA] * max( 0., sum1 + sum2);
259 
260  // Order single diffractive correctly.
261  if (swapped) {
262  swap( bB, bA);
263  swap( sigXB, sigAX);
264  swap( mMinXBsave, mMinAXsave);
265  swap( mResXBsave, mResAXsave);
266  }
267 
268  // Double diffractive scattering A + B -> X1 + X2 cross section.
269  double y0min = log( s * SPROTON / (sMinXB * sMinAX) ) ;
270  double sLog = log(s);
271  double Delta0 = CDD[iDD][0] + CDD[iDD][1] / sLog
272  + CDD[iDD][2] / pow2(sLog);
273  sum1 = (y0min * (log( max( 1e-10, y0min/Delta0) ) - 1.) + Delta0)/ alP2;
274  if (y0min < 0.) sum1 = 0.;
275  double sMaxXX = s * ( CDD[iDD][3] + CDD[iDD][4] / sLog
276  + CDD[iDD][5] / pow2(sLog) );
277  double sLogUp = log( max( 1.1, s * s0 / (sMinXB * sRMavgAX) ));
278  double sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgAX) ));
279  sum2 = CRES * log( sLogUp / sLogDn ) * sRMlogAX / alP2;
280  sLogUp = log( max( 1.1, s * s0 / (sMinAX * sRMavgXB) ));
281  sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgXB) ));
282  sum3 = CRES * log(sLogUp / sLogDn) * sRMlogXB / alP2;
283  double BcorrXX = CDD[iDD][6] + CDD[iDD][7] / eCM + CDD[iDD][8] / s;
284  sum4 = pow2(CRES) * sRMlogAX * sRMlogXB
285  / max( 0.1, alP2 * log( s * s0 / (sRMavgAX * sRMavgXB) ) + BcorrXX);
286  sigXX = CONVERTDD * X[iProc] * max( 0., sum1 + sum2 + sum3 + sum4);
287 
288  // Option with user-requested damping of diffractive cross sections.
289  if (doDampen) {
290  sigXB = sigXB * maxXBOwn / (sigXB + maxXBOwn);
291  sigAX = sigAX * maxAXOwn / (sigAX + maxAXOwn);
292  sigXX = sigXX * maxXXOwn / (sigXX + maxXXOwn);
293  }
294 
295  // Option with user-set values for total and partial cross sections.
296  // (Is not done earlier since want diffractive slopes anyway.)
297  double sigNDOwn = sigTotOwn - sigElOwn - sigXBOwn - sigAXOwn - sigXXOwn;
298  double sigElMax = sigEl;
299  if (setTotal && sigNDOwn > 0.) {
300  sigTot = sigTotOwn;
301  sigEl = sigElOwn;
302  sigXB = sigXBOwn;
303  sigAX = sigAXOwn;
304  sigXX = sigXXOwn;
305  sigElMax = sigEl;
306 
307  // Sub-option to set elastic parameters, including Coulomb contribution.
308  if (setElastic) {
309  bEl = bSlope;
310  sigEl = CONVERTEL * pow2(sigTot) * (1. + rho*rho) / bSlope;
311  sigElMax = 2. * (sigEl * exp(-bSlope * tAbsMin)
312  + alphaEM0 * alphaEM0 / (4. * CONVERTEL * tAbsMin) );
313  }
314  }
315 
316  // Inelastic nondiffractive by unitarity.
317  sigND = sigTot - sigEl - sigXB - sigAX - sigXX;
318  if (sigND < 0.) infoPtr->errorMsg("Error in SigmaTotal::init: "
319  "sigND < 0");
320  else if (sigND < 0.4 * sigTot) infoPtr->errorMsg("Warning in "
321  "SigmaTotal::init: sigND suspiciously low");
322 
323  // Upper estimate of elastic, including Coulomb term, where appropriate.
324  sigEl = sigElMax;
325 
326  // Done.
327  isCalc = true;
328  return true;
329 
330 }
331 
332 //==========================================================================
333 
334 } // end namespace Pythia8
Generated by   doxygen 1.8.5