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