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) 2014 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 "Pythia8/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 // For now a neutron is treated like a proton.
30 
31 //--------------------------------------------------------------------------
32 
33 // Definitions of static variables.
34 // Note that a lot of parameters are hardcoded as const here, rather
35 // than being interfaced for public change, since any changes would
36 // have to be done in a globally consistent manner. Which basically
37 // means a rewrite/replacement of the whole class.
38 
39 // Minimum threshold below which no cross sections will be defined.
40 const double SigmaTotal::MMIN = 2.;
41 
42 // General constants in total cross section parametrization:
43 // sigmaTot = X * s^epsilon + Y * s^eta (pomeron + reggeon).
44 const double SigmaTotal::EPSILON = 0.0808;
45 const double SigmaTotal::ETA = -0.4525;
46 const double SigmaTotal::X[] = { 21.70, 21.70, 13.63, 13.63, 13.63,
47  10.01, 0.970, 8.56, 6.29, 0.609, 4.62, 0.447, 0.0434};
48 const double SigmaTotal::Y[] = { 56.08, 98.39, 27.56, 36.02, 31.79,
49  1.51, -0.146, 13.08, -0.62, -0.060, 0.030, -0.0028, 0.00028};
50 
51 // Type of the two incoming hadrons as function of the process number:
52 // = 0 : p/n ; = 1 : pi/rho/omega; = 2 : phi; = 3 : J/psi.
53 const int SigmaTotal::IHADATABLE[] = { 0, 0, 1, 1, 1, 2, 3, 1, 1,
54  1, 2, 2, 3};
55 const int SigmaTotal::IHADBTABLE[] = { 0, 0, 0, 0, 0, 0, 0, 1, 2,
56  3, 2, 3, 3};
57 
58 // Hadron-Pomeron coupling beta(t) = beta(0) * exp(b*t).
59 const double SigmaTotal::BETA0[] = { 4.658, 2.926, 2.149, 0.208};
60 const double SigmaTotal::BHAD[] = { 2.3, 1.4, 1.4, 0.23};
61 
62 // Pomeron trajectory alpha(t) = 1 + epsilon + alpha' * t
63 const double SigmaTotal::ALPHAPRIME = 0.25;
64 
65 // Conversion coefficients = 1/(16pi) * (mb <-> GeV^2) * (G_3P)^n,
66 // with n = 0 elastic, n = 1 single and n = 2 double diffractive.
67 const double SigmaTotal::CONVERTEL = 0.0510925;
68 const double SigmaTotal::CONVERTSD = 0.0336;
69 const double SigmaTotal::CONVERTDD = 0.0084;
70 
71 // Diffractive mass spectrum starts at m + MMIN0 and has a low-mass
72 // enhancement, factor cRes, up to around m + mRes0.
73 const double SigmaTotal::MMIN0 = 0.28;
74 const double SigmaTotal::CRES = 2.0;
75 const double SigmaTotal::MRES0 = 1.062;
76 
77 // Parameters and coefficients for single diffractive scattering.
78 const int SigmaTotal::ISDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
79  6, 7, 8, 9};
80 const double SigmaTotal::CSD[10][8] = {
81  { 0.213, 0.0, -0.47, 150., 0.213, 0.0, -0.47, 150., } ,
82  { 0.213, 0.0, -0.47, 150., 0.267, 0.0, -0.47, 100., } ,
83  { 0.213, 0.0, -0.47, 150., 0.232, 0.0, -0.47, 110., } ,
84  { 0.213, 7.0, -0.55, 800., 0.115, 0.0, -0.47, 110., } ,
85  { 0.267, 0.0, -0.46, 75., 0.267, 0.0, -0.46, 75., } ,
86  { 0.232, 0.0, -0.46, 85., 0.267, 0.0, -0.48, 100., } ,
87  { 0.115, 0.0, -0.50, 90., 0.267, 6.0, -0.56, 420., } ,
88  { 0.232, 0.0, -0.48, 110., 0.232, 0.0, -0.48, 110., } ,
89  { 0.115, 0.0, -0.52, 120., 0.232, 6.0, -0.56, 470., } ,
90  { 0.115, 5.5, -0.58, 570., 0.115, 5.5, -0.58, 570. } };
91 
92 // Parameters and coefficients for double diffractive scattering.
93 const int SigmaTotal::IDDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
94  6, 7, 8, 9};
95 const double SigmaTotal::CDD[10][9] = {
96  { 3.11, -7.34, 9.71, 0.068, -0.42, 1.31, -1.37, 35.0, 118., } ,
97  { 3.11, -7.10, 10.6, 0.073, -0.41, 1.17, -1.41, 31.6, 95., } ,
98  { 3.12, -7.43, 9.21, 0.067, -0.44, 1.41, -1.35, 36.5, 132., } ,
99  { 3.13, -8.18, -4.20, 0.056, -0.71, 3.12, -1.12, 55.2, 1298., } ,
100  { 3.11, -6.90, 11.4, 0.078, -0.40, 1.05, -1.40, 28.4, 78., } ,
101  { 3.11, -7.13, 10.0, 0.071, -0.41, 1.23, -1.34, 33.1, 105., } ,
102  { 3.12, -7.90, -1.49, 0.054, -0.64, 2.72, -1.13, 53.1, 995., } ,
103  { 3.11, -7.39, 8.22, 0.065, -0.44, 1.45, -1.36, 38.1, 148., } ,
104  { 3.18, -8.95, -3.37, 0.057, -0.76, 3.32, -1.12, 55.6, 1472., } ,
105  { 4.18, -29.2, 56.2, 0.074, -1.36, 6.67, -1.14, 116.2, 6532. } };
106 const double SigmaTotal::SPROTON = 0.880;
107 
108 // MBR parameters. Integration of MBR cross section.
109 const int SigmaTotal::NINTEG = 1000;
110 const int SigmaTotal::NINTEG2 = 40;
111 const double SigmaTotal::HBARC2 = 0.38938;
112 // MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV^-2.
113 const double SigmaTotal::FFA1 = 0.9;
114 const double SigmaTotal::FFA2 = 0.1;
115 const double SigmaTotal::FFB1 = 4.6;
116 const double SigmaTotal::FFB2 = 0.6;
117 
118 //--------------------------------------------------------------------------
119 
120 // Store pointer to Info and initialize data members.
121 
122 void SigmaTotal::init(Info* infoPtrIn, Settings& settings,
123  ParticleData* particleDataPtrIn) {
124 
125  // Store pointers.
126  infoPtr = infoPtrIn;
127  particleDataPtr = particleDataPtrIn;
128 
129  // Normalization of central diffractive cross section.
130  zeroAXB = settings.flag("SigmaTotal:zeroAXB");
131  sigAXB2TeV = settings.parm("SigmaTotal:sigmaAXB2TeV");
132 
133  // User-set values for cross sections.
134  setTotal = settings.flag("SigmaTotal:setOwn");
135  sigTotOwn = settings.parm("SigmaTotal:sigmaTot");
136  sigElOwn = settings.parm("SigmaTotal:sigmaEl");
137  sigXBOwn = settings.parm("SigmaTotal:sigmaXB");
138  sigAXOwn = settings.parm("SigmaTotal:sigmaAX");
139  sigXXOwn = settings.parm("SigmaTotal:sigmaXX");
140  sigAXBOwn = settings.parm("SigmaTotal:sigmaAXB");
141 
142  // User-set values to dampen diffractive cross sections.
143  doDampen = settings.flag("SigmaDiffractive:dampen");
144  maxXBOwn = settings.parm("SigmaDiffractive:maxXB");
145  maxAXOwn = settings.parm("SigmaDiffractive:maxAX");
146  maxXXOwn = settings.parm("SigmaDiffractive:maxXX");
147  maxAXBOwn = settings.parm("SigmaDiffractive:maxAXB");
148 
149  // User-set values for handling of elastic sacattering.
150  setElastic = settings.flag("SigmaElastic:setOwn");
151  bSlope = settings.parm("SigmaElastic:bSlope");
152  rho = settings.parm("SigmaElastic:rho");
153  lambda = settings.parm("SigmaElastic:lambda");
154  tAbsMin = settings.parm("SigmaElastic:tAbsMin");
155  alphaEM0 = settings.parm("StandardModel:alphaEM0");
156 
157  // Parameters for diffractive systems.
158  sigmaPomP = settings.parm("Diffraction:sigmaRefPomP");
159  mPomP = settings.parm("Diffraction:mRefPomP");
160  pPomP = settings.parm("Diffraction:mPowPomP");
161 
162  // Parameters for MBR model.
163  PomFlux = settings.mode("Diffraction:PomFlux");
164  MBReps = settings.parm("Diffraction:MBRepsilon");
165  MBRalpha = settings.parm("Diffraction:MBRalpha");
166  MBRbeta0 = settings.parm("Diffraction:MBRbeta0");
167  MBRsigma0 = settings.parm("Diffraction:MBRsigma0");
168  m2min = settings.parm("Diffraction:MBRm2Min");
169  dyminSDflux = settings.parm("Diffraction:MBRdyminSDflux");
170  dyminDDflux = settings.parm("Diffraction:MBRdyminDDflux");
171  dyminCDflux = settings.parm("Diffraction:MBRdyminCDflux");
172  dyminSD = settings.parm("Diffraction:MBRdyminSD");
173  dyminDD = settings.parm("Diffraction:MBRdyminDD");
174  dyminCD = settings.parm("Diffraction:MBRdyminCD");
175  dyminSigSD = settings.parm("Diffraction:MBRdyminSigSD");
176  dyminSigDD = settings.parm("Diffraction:MBRdyminSigDD");
177  dyminSigCD = settings.parm("Diffraction:MBRdyminSigCD");
178 
179 }
180 
181 //--------------------------------------------------------------------------
182 
183 // Function that calculates the relevant properties.
184 
185 bool SigmaTotal::calc( int idA, int idB, double eCM) {
186 
187  // Derived quantities.
188  alP2 = 2. * ALPHAPRIME;
189  s0 = 1. / ALPHAPRIME;
190 
191  // Reset everything to zero to begin with.
192  isCalc = false;
193  sigTot = sigEl = sigXB = sigAX = sigXX = sigAXB = sigND = bEl = s
194  = bA = bB = 0.;
195 
196  // Order flavour of incoming hadrons: idAbsA < idAbsB (restore later).
197  int idAbsA = abs(idA);
198  int idAbsB = abs(idB);
199  bool swapped = false;
200  if (idAbsA > idAbsB) {
201  swap( idAbsA, idAbsB);
202  swapped = true;
203  }
204  double sameSign = (idA * idB > 0);
205 
206  // Find process number.
207  int iProc = -1;
208  if (idAbsA > 1000) {
209  iProc = (sameSign) ? 0 : 1;
210  } else if (idAbsA > 100 && idAbsB > 1000) {
211  iProc = (sameSign) ? 2 : 3;
212  if (idAbsA/10 == 11 || idAbsA/10 == 22) iProc = 4;
213  if (idAbsA > 300) iProc = 5;
214  if (idAbsA > 400) iProc = 6;
215  if (idAbsA > 900) iProc = 13;
216  } else if (idAbsA > 100) {
217  iProc = 7;
218  if (idAbsB > 300) iProc = 8;
219  if (idAbsB > 400) iProc = 9;
220  if (idAbsA > 300) iProc = 10;
221  if (idAbsA > 300 && idAbsB > 400) iProc = 11;
222  if (idAbsA > 400) iProc = 12;
223  }
224  if (iProc == -1) return false;
225 
226  // Primitive implementation of Pomeron + p.
227  if (iProc == 13) {
228  s = eCM*eCM;
229  sigTot = sigmaPomP * pow( eCM / mPomP, pPomP);
230  sigND = sigTot;
231  isCalc = true;
232  return true;
233  }
234 
235  // Find hadron masses and check that energy is enough.
236  // For mesons use the corresponding vector meson masses.
237  int idModA = (idAbsA > 1000) ? idAbsA : 10 * (idAbsA/10) + 3;
238  int idModB = (idAbsB > 1000) ? idAbsB : 10 * (idAbsB/10) + 3;
239  double mA = particleDataPtr->m0(idModA);
240  double mB = particleDataPtr->m0(idModB);
241  if (eCM < mA + mB + MMIN) {
242  infoPtr->errorMsg("Error in SigmaTotal::calc: too low energy");
243  return false;
244  }
245 
246  // Evaluate the total cross section.
247  s = eCM*eCM;
248  double sEps = pow( s, EPSILON);
249  double sEta = pow( s, ETA);
250  sigTot = X[iProc] * sEps + Y[iProc] * sEta;
251 
252  // Slope of hadron form factors.
253  int iHadA = IHADATABLE[iProc];
254  int iHadB = IHADBTABLE[iProc];
255  bA = BHAD[iHadA];
256  bB = BHAD[iHadB];
257 
258  // Elastic slope parameter and cross section.
259  bEl = 2.*bA + 2.*bB + 4.*sEps - 4.2;
260  sigEl = CONVERTEL * pow2(sigTot) / bEl;
261 
262  // Lookup coefficients for single and double diffraction.
263  int iSD = ISDTABLE[iProc];
264  int iDD = IDDTABLE[iProc];
265  double sum1, sum2, sum3, sum4;
266 
267  // Single diffractive scattering A + B -> X + B cross section.
268  mMinXBsave = mA + MMIN0;
269  double sMinXB = pow2(mMinXBsave);
270  mResXBsave = mA + MRES0;
271  double sResXB = pow2(mResXBsave);
272  double sRMavgXB = mResXBsave * mMinXBsave;
273  double sRMlogXB = log(1. + sResXB/sMinXB);
274  double sMaxXB = CSD[iSD][0] * s + CSD[iSD][1];
275  double BcorrXB = CSD[iSD][2] + CSD[iSD][3] / s;
276  sum1 = log( (2.*bB + alP2 * log(s/sMinXB))
277  / (2.*bB + alP2 * log(s/sMaxXB)) ) / alP2;
278  sum2 = CRES * sRMlogXB / (2.*bB + alP2 * log(s/sRMavgXB) + BcorrXB) ;
279  sigXB = CONVERTSD * X[iProc] * BETA0[iHadB] * max( 0., sum1 + sum2);
280 
281  // Single diffractive scattering A + B -> A + X cross section.
282  mMinAXsave = mB + MMIN0;
283  double sMinAX = pow2(mMinAXsave);
284  mResAXsave = mB + MRES0;
285  double sResAX = pow2(mResAXsave);
286  double sRMavgAX = mResAXsave * mMinAXsave;
287  double sRMlogAX = log(1. + sResAX/sMinAX);
288  double sMaxAX = CSD[iSD][4] * s + CSD[iSD][5];
289  double BcorrAX = CSD[iSD][6] + CSD[iSD][7] / s;
290  sum1 = log( (2.*bA + alP2 * log(s/sMinAX))
291  / (2.*bA + alP2 * log(s/sMaxAX)) ) / alP2;
292  sum2 = CRES * sRMlogAX / (2.*bA + alP2 * log(s/sRMavgAX) + BcorrAX) ;
293  sigAX = CONVERTSD * X[iProc] * BETA0[iHadA] * max( 0., sum1 + sum2);
294 
295  // Order single diffractive correctly.
296  if (swapped) {
297  swap( bB, bA);
298  swap( sigXB, sigAX);
299  swap( mMinXBsave, mMinAXsave);
300  swap( mResXBsave, mResAXsave);
301  }
302 
303  // Double diffractive scattering A + B -> X1 + X2 cross section.
304  double y0min = log( s * SPROTON / (sMinXB * sMinAX) ) ;
305  double sLog = log(s);
306  double Delta0 = CDD[iDD][0] + CDD[iDD][1] / sLog
307  + CDD[iDD][2] / pow2(sLog);
308  sum1 = (y0min * (log( max( 1e-10, y0min/Delta0) ) - 1.) + Delta0)/ alP2;
309  if (y0min < 0.) sum1 = 0.;
310  double sMaxXX = s * ( CDD[iDD][3] + CDD[iDD][4] / sLog
311  + CDD[iDD][5] / pow2(sLog) );
312  double sLogUp = log( max( 1.1, s * s0 / (sMinXB * sRMavgAX) ));
313  double sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgAX) ));
314  sum2 = CRES * log( sLogUp / sLogDn ) * sRMlogAX / alP2;
315  sLogUp = log( max( 1.1, s * s0 / (sMinAX * sRMavgXB) ));
316  sLogDn = log( max( 1.1, s * s0 / (sMaxXX * sRMavgXB) ));
317  sum3 = CRES * log(sLogUp / sLogDn) * sRMlogXB / alP2;
318  double BcorrXX = CDD[iDD][6] + CDD[iDD][7] / eCM + CDD[iDD][8] / s;
319  sum4 = pow2(CRES) * sRMlogAX * sRMlogXB
320  / max( 0.1, alP2 * log( s * s0 / (sRMavgAX * sRMavgXB) ) + BcorrXX);
321  sigXX = CONVERTDD * X[iProc] * max( 0., sum1 + sum2 + sum3 + sum4);
322 
323  // Central diffractive scattering A + B -> A + X + B, only p and pbar.
324  mMinAXBsave = 1.;
325  if ( (idAbsA == 2212 || idAbsA == 2112)
326  && (idAbsB == 2212 || idAbsB == 2112) && !zeroAXB) {
327  double sMinAXB = pow2(mMinAXBsave);
328  double sRefAXB = pow2(2000.);
329  sigAXB = sigAXB2TeV * pow( log(0.06 * s / sMinAXB), 1.5 )
330  / pow( log(0.06 * sRefAXB / sMinAXB), 1.5 );
331  }
332 
333  // Option with user-requested damping of diffractive cross sections.
334  if (doDampen) {
335  sigXB = sigXB * maxXBOwn / (sigXB + maxXBOwn);
336  sigAX = sigAX * maxAXOwn / (sigAX + maxAXOwn);
337  sigXX = sigXX * maxXXOwn / (sigXX + maxXXOwn);
338  sigAXB = sigAXB * maxAXBOwn / (sigAXB + maxAXBOwn);
339  }
340 
341  // Calculate cross sections in MBR model.
342  if (PomFlux == 5) calcMBRxsecs(idA, idB, eCM);
343 
344  // Option with user-set values for total and partial cross sections.
345  // (Is not done earlier since want diffractive slopes anyway.)
346  double sigNDOwn = sigTotOwn - sigElOwn - sigXBOwn - sigAXOwn - sigXXOwn
347  - sigAXBOwn;
348  double sigElMax = sigEl;
349  if (setTotal && sigNDOwn > 0.) {
350  sigTot = sigTotOwn;
351  sigEl = sigElOwn;
352  sigXB = sigXBOwn;
353  sigAX = sigAXOwn;
354  sigXX = sigXXOwn;
355  sigAXB = sigAXBOwn;
356  sigElMax = sigEl;
357 
358  // Sub-option to set elastic parameters, including Coulomb contribution.
359  if (setElastic) {
360  bEl = bSlope;
361  sigEl = CONVERTEL * pow2(sigTot) * (1. + rho*rho) / bSlope;
362  sigElMax = 2. * (sigEl * exp(-bSlope * tAbsMin)
363  + alphaEM0 * alphaEM0 / (4. * CONVERTEL * tAbsMin) );
364  }
365  }
366 
367  // Inelastic nondiffractive by unitarity.
368  sigND = sigTot - sigEl - sigXB - sigAX - sigXX - sigAXB;
369  if (sigND < 0.) infoPtr->errorMsg("Error in SigmaTotal::init: "
370  "sigND < 0");
371  else if (sigND < 0.4 * sigTot) infoPtr->errorMsg("Warning in "
372  "SigmaTotal::init: sigND suspiciously low");
373 
374  // Upper estimate of elastic, including Coulomb term, where appropriate.
375  sigEl = sigElMax;
376 
377  // Done.
378  isCalc = true;
379  return true;
380 
381 }
382 
383 //--------------------------------------------------------------------------
384 
385 // Calculate parameters in the MBR model.
386 
387 bool SigmaTotal::calcMBRxsecs( int idA, int idB, double eCM) {
388 
389  // Local variables.
390  double sigtot, sigel, sigsd, sigdd, sigdpe;
391 
392  // MBR parameters locally.
393  double eps = MBReps;
394  double alph = MBRalpha;
395  double beta0gev = MBRbeta0;
396  double beta0mb = beta0gev * sqrt(HBARC2);
397  double sigma0mb = MBRsigma0;
398  double sigma0gev = sigma0mb/HBARC2;
399  double a1 = FFA1;
400  double a2 = FFA2;
401  double b1 = FFB1;
402  double b2 = FFB2;
403 
404  // Calculate total and elastic cross sections.
405  double ratio;
406  if (eCM <= 1800.0) {
407  double sign = (idA * idB > 0);
408  sigtot = 16.79 * pow(s, 0.104) + 60.81 * pow(s, -0.32)
409  - sign * 31.68 * pow(s, -0.54);
410  ratio = 0.100 * pow(s, 0.06) + 0.421 * pow(s, -0.52)
411  + sign * 0.160 * pow(s, -0.6);
412  } else {
413  double sigCDF = 80.03;
414  double sCDF = pow2(1800.);
415  double sF = pow2(22.);
416  sigtot = sigCDF + ( pow2( log(s / sF)) - pow2( log(sCDF / sF)) )
417  * M_PI / (3.7 / HBARC2);
418  ratio = 0.066 + 0.0119 * log(s);
419  }
420  sigel=sigtot*ratio;
421 
422  // Integrate SD, DD and DPE(CD) cross sections.
423  // Each cross section is obtained from the ratio of two integrals:
424  // the Regge cross section and the renormalized flux.
425  double cflux, csig, c1, step, f;
426  double dymin0 = 0.;
427 
428  // Calculate SD cross section.
429  double dymaxSD = log(s / m2min);
430  cflux = pow2(beta0gev) / (16. * M_PI);
431  csig = cflux * sigma0mb;
432 
433  // SD flux.
434  c1 = cflux;
435  double fluxsd = 0.;
436  step = (dymaxSD - dyminSDflux) / NINTEG;
437  for (int i = 0; i < NINTEG; ++i) {
438  double dy = dyminSDflux + (i + 0.5) * step;
439  f = exp(2. * eps * dy) * ( (a1 / (b1 + 2. * alph * dy))
440  + (a2 / (b2 + 2. * alph * dy)) );
441  f *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD));
442  fluxsd = fluxsd + step * c1 * f;
443  }
444  if (fluxsd < 1.) fluxsd = 1.;
445 
446  // Regge cross section.
447  c1 = csig * pow(s, eps);
448  sigsd = 0.;
449  sdpmax = 0.;
450  step = (dymaxSD - dymin0) / NINTEG;
451  for (int i = 0; i < NINTEG; ++i) {
452  double dy = dymin0 + (i + 0.5) * step;
453  f = exp(eps * dy) * ( (a1 / (b1 + 2. * alph * dy))
454  + (a2 / (b2 + 2. * alph * dy)) );
455  f *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD));
456  if (f > sdpmax) sdpmax = f;
457  sigsd = sigsd + step * c1 * f;
458  }
459  sdpmax *= 1.01;
460  sigsd /= fluxsd;
461 
462  // Calculate DD cross section.
463  // Note: dymaxDD = ln(s * s0 /mMin^4) with s0 = 1 GeV^2.
464  double dymaxDD = log(s / pow2(m2min));
465  cflux = sigma0gev / (16. * M_PI);
466  csig = cflux * sigma0mb;
467 
468  // DD flux.
469  c1 = cflux / (2. * alph);
470  double fluxdd = 0.;
471  step = (dymaxDD - dyminDDflux) / NINTEG;
472  for (int i = 0; i < NINTEG; ++i) {
473  double dy = dyminDDflux + (i + 0.5) * step;
474  f = (dymaxDD - dy) * exp(2. * eps * dy)
475  * ( exp(-2. * alph * dy * exp(-dy))
476  - exp(-2. * alph * dy * exp(dy)) ) / dy;
477  f *= 0.5 * (1. + erf( (dy - dyminDD) / dyminSigDD));
478  fluxdd = fluxdd + step * c1 * f;
479  }
480  if (fluxdd < 1.) fluxdd = 1.;
481 
482  // Regge cross section.
483  c1 = csig * pow(s, eps) / (2. * alph);
484  ddpmax = 0.;
485  sigdd = 0.;
486  step = (dymaxDD - dymin0) / NINTEG;
487  for (int i = 0; i < NINTEG; ++i) {
488  double dy = dymin0 + (i + 0.5) * step;
489  f = (dymaxDD - dy) * exp(eps * dy)
490  * ( exp(-2. * alph * dy * exp(-dy))
491  - exp(-2. * alph * dy * exp(dy)) ) / dy;
492  f *= 0.5 * (1. + erf( (dy - dyminDD) / dyminSigDD));
493  if (f > ddpmax) ddpmax = f;
494  sigdd = sigdd + step * c1 * f;
495  }
496  ddpmax *= 1.01;
497  sigdd /= fluxdd;
498 
499  // Calculate DPE (CD) cross section.
500  double dymaxCD = log(s / m2min);
501  cflux = pow4(beta0gev) / pow2(16. * M_PI);
502  csig = cflux * pow2(sigma0mb / beta0mb);
503  double dy1, dy2, f1, f2, step2;
504 
505  // DPE flux.
506  c1 = cflux;
507  double fluxdpe = 0.;
508  step = (dymaxCD - dyminCDflux) / NINTEG;
509  for (int i = 0; i < NINTEG; ++i) {
510  double dy = dyminCDflux + (i + 0.5) * step;
511  f = 0.;
512  step2 = (dy - dyminCDflux) / NINTEG2;
513  for (int j = 0; j < NINTEG2; ++j) {
514  double yc = -0.5 * (dy - dyminCDflux) + (j + 0.5) * step2;
515  dy1 = 0.5 * dy - yc;
516  dy2 = 0.5 * dy + yc;
517  f1 = exp(2. * eps * dy1) * ( (a1 / (b1 + 2. * alph * dy1))
518  + (a2 / (b2 + 2. * alph * dy1)) );
519  f2 = exp(2. * eps * dy2) * ( (a1 / (b1 + 2. * alph * dy2))
520  + (a2 / (b2 + 2. * alph * dy2)) );
521  f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminCD)
522  / (dyminSigCD / sqrt(2.))) );
523  f2 *= 0.5 * (1. + erf( (dy2 - 0.5 *dyminCD)
524  / (dyminSigCD / sqrt(2.))) );
525  f += f1 * f2 * step2;
526  }
527  fluxdpe += step * c1 * f;
528  }
529  if (fluxdpe < 1.) fluxdpe = 1.;
530 
531  // Regge cross section.
532  c1 = csig * pow(s, eps);
533  sigdpe = 0.;
534  dpepmax = 0;
535  step = (dymaxCD - dymin0) / NINTEG;
536  for (int i = 0; i < NINTEG; ++i) {
537  double dy = dymin0 + (i + 0.5) * step;
538  f = 0.;
539  step2 = (dy - dymin0) / NINTEG2;
540  for (int j = 0; j < NINTEG2; ++j) {
541  double yc = -0.5 * (dy - dymin0) + (j + 0.5) * step2;
542  dy1 = 0.5 * dy - yc;
543  dy2 = 0.5 * dy + yc;
544  f1 = exp(eps * dy1) * ( (a1 / (b1 + 2. * alph * dy1))
545  + (a2 / (b2 + 2. * alph * dy1)) );
546  f2 = exp(eps * dy2) * ( (a1 / (b1 + 2. * alph * dy2))
547  + (a2 / (b2 + 2. * alph * dy2)) );
548  f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminCD)
549  / (dyminSigCD / sqrt(2.))) );
550  f2 *= 0.5 * (1. + erf( (dy2 - 0.5 * dyminCD)
551  /(dyminSigCD / sqrt(2.))) );
552  f += f1 * f2 * step2;
553  }
554  sigdpe += step * c1 * f;
555  if ( f > dpepmax) dpepmax = f;
556  }
557  dpepmax *= 1.01;
558  sigdpe /= fluxdpe;
559 
560  // Diffraction done. Now calculate total inelastic cross section.
561  sigND = sigtot - (2. * sigsd + sigdd + sigel + sigdpe);
562  sigTot = sigtot;
563  sigEl = sigel;
564  sigAX = sigsd;
565  sigXB = sigsd;
566  sigXX = sigdd;
567  sigAXB = sigdpe;
568 
569  return true;
570 }
571 
572 //==========================================================================
573 
574 } // end namespace Pythia8
Generated by   doxygen 1.8.5