]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8175/src/SigmaTotal.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
CID 21256: Uninitialized pointer field (UNINIT_CTOR)
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8175 / src / SigmaTotal.cxx
CommitLineData
c6b60c38 1// SigmaTotal.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2013 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
10namespace 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.
40const double SigmaTotal::MMIN = 2.;
41
42// General constants in total cross section parametrization:
43// sigmaTot = X * s^epsilon + Y * s^eta (pomeron + reggeon).
44const double SigmaTotal::EPSILON = 0.0808;
45const double SigmaTotal::ETA = -0.4525;
46const 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};
48const 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.
53const int SigmaTotal::IHADATABLE[] = { 0, 0, 1, 1, 1, 2, 3, 1, 1,
54 1, 2, 2, 3};
55const 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).
59const double SigmaTotal::BETA0[] = { 4.658, 2.926, 2.149, 0.208};
60const double SigmaTotal::BHAD[] = { 2.3, 1.4, 1.4, 0.23};
61
62// Pomeron trajectory alpha(t) = 1 + epsilon + alpha' * t
63const 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.
67const double SigmaTotal::CONVERTEL = 0.0510925;
68const double SigmaTotal::CONVERTSD = 0.0336;
69const 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.
73const double SigmaTotal::MMIN0 = 0.28;
74const double SigmaTotal::CRES = 2.0;
75const double SigmaTotal::MRES0 = 1.062;
76
77// Parameters and coefficients for single diffractive scattering.
78const int SigmaTotal::ISDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
79 6, 7, 8, 9};
80const 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.
93const int SigmaTotal::IDDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5,
94 6, 7, 8, 9};
95const 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. } };
106const double SigmaTotal::SPROTON = 0.880;
107
108// MBR parameters. Integration of MBR cross section.
109const int SigmaTotal::NINTEG = 1000;
110const int SigmaTotal::NINTEG2 = 40;
111const double SigmaTotal::HBARC2 = 0.38938;
112// MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV^-2.
113const double SigmaTotal::FFA1 = 0.9;
114const double SigmaTotal::FFA2 = 0.1;
115const double SigmaTotal::FFB1 = 4.6;
116const double SigmaTotal::FFB2 = 0.6;
117
118//--------------------------------------------------------------------------
119
120// Store pointer to Info and initialize data members.
121
122void 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
185bool 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
387bool 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