source: HiSusy/trunk/Pythia8/pythia8170/src/SigmaTotal.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 21.6 KB
Line 
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
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
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.
39const double SigmaTotal::MMIN  = 2.; 
40
41// General constants in total cross section parametrization:
42// sigmaTot = X * s^epsilon + Y * s^eta (pomeron + reggeon).
43const double SigmaTotal::EPSILON = 0.0808;
44const double SigmaTotal::ETA     = -0.4525;
45const 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};
47const 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.
52const int SigmaTotal::IHADATABLE[] = { 0, 0, 1, 1, 1, 2, 3, 1, 1, 
53  1, 2, 2, 3};
54const 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).
58const double SigmaTotal::BETA0[] = { 4.658, 2.926, 2.149, 0.208}; 
59const double SigmaTotal::BHAD[]  = {   2.3,   1.4,   1.4,  0.23};
60
61// Pomeron trajectory alpha(t) = 1 + epsilon + alpha' * t
62const 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.
66const double SigmaTotal::CONVERTEL = 0.0510925;
67const double SigmaTotal::CONVERTSD = 0.0336; 
68const 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.
72const double SigmaTotal::MMIN0 = 0.28; 
73const double SigmaTotal::CRES  = 2.0; 
74const double SigmaTotal::MRES0 = 1.062;
75
76// Parameters and coefficients for single diffractive scattering.
77const int SigmaTotal::ISDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5, 
78  6, 7, 8, 9}; 
79const 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.
92const int SigmaTotal::IDDTABLE[] = { 0, 0, 1, 1, 1, 2, 3, 4, 5, 
93  6, 7, 8, 9}; 
94const 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.  } };
105const double SigmaTotal::SPROTON = 0.880;
106 
107// MBR parameters. Integration of MBR cross section.
108const int    SigmaTotal::NINTEG = 1000;
109const int    SigmaTotal::NINTEG2 = 40;
110const double SigmaTotal::HBARC2 = 0.38938;
111// MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV^-2.
112const double SigmaTotal::FFA1 = 0.9;
113const double SigmaTotal::FFA2 = 0.1;
114const double SigmaTotal::FFB1 = 4.6;
115const double SigmaTotal::FFB2 = 0.6;
116
117//--------------------------------------------------------------------------
118
119// Store pointer to Info and initialize data members.
120
121void 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
184bool 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
385bool 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
Note: See TracBrowser for help on using the repository browser.