StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
HardDiffraction.cc
1 // HardDiffraction.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Author: Christine O. Rasmussen.
7 
8 // Function definitions (not found in the header) for the
9 // HardDiffraction class.
10 
11 #include "Pythia8/HardDiffraction.h"
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // Constants: could be changed here if desired, but normally should not.
17 // These are of technical nature, as described for each.
18 
19 // Lower limit on PDF value in order to avoid division by zero.
20 const double HardDiffraction::TINYPDF = 1e-10;
21 
22 // Ficticious Pomeron mass to leave room for beam remnant
23 const double HardDiffraction::POMERONMASS = 1.;
24 const double HardDiffraction::RHOMASS = 0.77549;
25 const double HardDiffraction::PROTONMASS = 0.93827;
26 
27 // Safetymargin for diffractive masses
28 const double HardDiffraction::DIFFMASSMARGIN = 0.2;
29 //--------------------------------------------------------------------------
30 
31 void HardDiffraction::init(BeamParticle* beamAPtrIn,
32  BeamParticle* beamBPtrIn) {
33 
34  // The beam pointers may differ from originally set in PhysicsBase.
35  beamAPtr = beamAPtrIn;
36  beamBPtr = beamBPtrIn;
37 
38  // Set diffraction parameters.
39  pomFlux = mode("SigmaDiffractive:PomFlux");
40 
41  // Read out some properties of beams to allow shorthand.
42  idA = (beamAPtr != 0) ? beamAPtr->id() : 0;
43  idB = (beamBPtr != 0) ? beamBPtr->id() : 0;
44  mA = (beamAPtr != 0) ? beamAPtr->m() : 0.;
45  mB = (beamBPtr != 0) ? beamBPtr->m() : 0.;
46  isGammaA = (beamAPtr != 0) ? beamAPtr->isGamma() : false;
47  isGammaB = (beamBPtr != 0) ? beamBPtr->isGamma() : false;
48  isGammaGamma = (isGammaA && isGammaB);
49 
50  // Set up Pomeron flux constants.
51  rescale = parm("Diffraction:PomFluxRescale");
52  a0 = 1. + parm("SigmaDiffractive:PomFluxEpsilon");
53  ap = parm("SigmaDiffractive:PomFluxAlphaPrime");
54 
55  if (pomFlux == 1) {
56  double sigmaRefPomP = parm("Diffraction:sigmaRefPomP");
57  normPom = pow2(sigmaRefPomP) * 0.02;
58  b0 = 2.3;
59  } else if (pomFlux == 2) {
60  normPom = 1/2.3;
61  A1 = 6.38;
62  A2 = 0.424;
63  a1 = 8.;
64  a2 = 3.;
65  } else if (pomFlux == 3) {
66  double beta = 10.;
67  normPom = pow2(beta)/(16.*M_PI);
68  a1 = 4.7;
69  } else if (pomFlux == 4) {
70  double beta = 1.8;
71  normPom = 9. * pow2(beta) / (4. * pow2(M_PI));
72  A1 = 0.27;
73  a1 = 8.38;
74  A2 = 0.56;
75  a2 = 3.78;
76  A3 = 0.18;
77  a3 = 1.36;
78  } else if (pomFlux == 5) {
79  A1 = 0.9;
80  a1 = 4.6;
81  A2 = 0.1;
82  a2 = 0.6;
83  a0 = 1. + parm("SigmaDiffractive:MBRepsilon");
84  ap = parm("SigmaDiffractive:MBRalpha");
85  bool renormalize = flag("Diffraction:useMBRrenormalization");
86  double cflux = 0.858;
87  double m2min = parm("SigmaDiffractive:MBRm2Min");
88  double dyminSDflux = parm("SigmaDiffractive:MBRdyminSDflux");
89  double dymaxSD = log(infoPtr->eCM()*infoPtr->eCM() / m2min);
90  double nGap = 0.;
91  if (renormalize){
92  double step = (dymaxSD - dyminSDflux) / 1000.;
93  // Calculate the integral of the flux
94  // to renormalize the gap:
95  for (int i = 0; i < 1000; ++i) {
96  double dy = dyminSDflux + (i + 0.5) * step;
97  double f = exp(2.*(a0 - 1.)*dy) * ( (A1/(a1 + 2.*ap*dy))
98  + (A2/(a2 + 2. * ap*dy)) );
99  nGap += step * cflux * f;
100  }
101  }
102  if (nGap < 1.) nGap = 1.;
103  normPom = cflux/nGap;
104  } else if (pomFlux == 6 || pomFlux == 7) {
105  // Has fixed values of eps and alpha' to get normalisation correct
106  ap = 0.06;
107  b0 = 5.5;
108  if (pomFlux == 6) a0 = 1.1182;
109  else a0 = 1.1110;
110  double xNorm = 0.003;
111  double b = b0 + 2. * ap * log(1./xNorm);
112  double mMin = (isGammaA || isGammaB) ? RHOMASS : PROTONMASS;
113  double tmin = -pow(mMin * xNorm, 2.)/(1. - xNorm);
114  double tcut = -1.;
115  double fNorm = exp(log(1./xNorm) * ( 2.*a0 - 2.));
116  fNorm *= (exp(b*tmin) - exp(b*tcut))/b;
117  normPom = 1./fNorm;
118  }
119 
120  // Initialise Pomeron values to zero.
121  xPomA = tPomA = thetaPomA = 0.;
122  xPomB = tPomB = thetaPomB = 0.;
123 
124  // Calculate rescaling factor for Pomeron flux in photons.
125  sigTotRatio = 1.;
126  if (isGammaA || isGammaB) {
127  sigmaTotPtr->calc(22, 2212, infoPtr->eCM());
128  double sigGamP = sigmaTotPtr->sigmaTot();
129  sigmaTotPtr->calc(2212, 2212, infoPtr->eCM());
130  double sigPP = sigmaTotPtr->sigmaTot();
131  sigTotRatio = sigGamP / sigPP;
132  }
133 
134  // Done.
135 }
136 
137 //--------------------------------------------------------------------------
138 
139 bool HardDiffraction::isDiffractive( int iBeamIn, int partonIn,
140  double xIn, double Q2In, double xfIncIn) {
141 
142  // iBeam = 1 means A B -> A' + (Pom + B) -> A' X
143  // iBeam = 2 means A B -> (A + Pom) + B' -> X B'
144 
145  // Store incoming values.
146  iBeam = iBeamIn;
147  int parton = partonIn;
148  double x = xIn;
149  double Q2 = Q2In;
150  double xfInc = xfIncIn;
151  tmpPomPtr = (iBeam == 1) ? beamPomAPtr : beamPomBPtr;
152  usePomInPhoton = ((iBeam == 1 && isGammaA) || (iBeam == 2 && isGammaB))
153  ? true : false;
154 
155  // Return false if value of inclusive PDF is zero.
156  if (xfInc < TINYPDF) {
157  infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
158  "inclusive PDF is zero");
159  return false;
160  }
161 
162  // Generate an xNow = x_P according to dx_P / x_P.
163  double xNow = pow(x, rndmPtr->flat());
164 
165  // Find estimate of diffractive PDF based on x_P choice above.
166  // f_i(xP) = int_x^1 d(xP)/xP * xP f_{P/p}(xP) * x/xP f_{i/P}(x/xP, Q^2)
167  // = ln (1/x) < xP f_{P/p}(xP) * x/xP f_{i/P}(x/xP, Q^2) >
168  double xfEst = log(1./x) * xfPom(xNow) * tmpPomPtr->xf(parton, x/xNow, Q2);
169 
170  // Warn if the estimated function exceeds the inclusive PDF.
171  if (xfEst > xfInc) {
172  ostringstream msg;
173  msg << ", id = " << parton;
174  infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
175  "weight above unity", msg.str());
176  }
177  // Discard if estimate/inclusive PDF is less than random number.
178  if (xfEst < rndmPtr->flat() * xfInc) return false;
179 
180  // Make sure there is momentum left for beam remnant.
181  // Make Pomeron massless and let mMin be the mass of the subcollision
182  // particle: if Pp then mMin = mp, if Pgamma then mMin = mrho.
183  double mMin = (usePomInPhoton) ? RHOMASS : PROTONMASS;
184  double m2Diff = xNow * pow2( infoPtr->eCM());
185  double mDiff = sqrt(m2Diff);
186  double mDiffA = (iBeam == 1) ? 0. : mMin;
187  double mDiffB = (iBeam == 2) ? 0. : mMin;
188  double m2DiffA = mDiffA * mDiffA;
189  double m2DiffB = mDiffB * mDiffB;
190  double eDiff = (iBeam == 1)
191  ? 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff
192  : 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
193  if ( 1. - x / xNow < POMERONMASS / eDiff) {
194  infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
195  "No momentum left for beam remnant.");
196  return false;
197  }
198 
199  // Make sure that the diffractive mass is not too high.
200  double m3 = ((iBeam == 1 && isGammaA) || (iBeam == 2 && isGammaB))
201  ? RHOMASS : PROTONMASS;
202  double m4 = mDiff;
203  if (m3 + m4 + DIFFMASSMARGIN >= infoPtr->eCM()) {
204  infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
205  "Too high diffractive mass.");
206  return false;
207  }
208 
209  // The chosen xNow is accepted, now find t and theta.
210  double tNow = pickTNow(xNow);
211  double thetaNow = getThetaNow(xNow, tNow);
212 
213  // Set the chosen diffractive values, to be able to use them later.
214  if (iBeam == 1) {
215  xPomA = xNow;
216  tPomA = tNow;
217  thetaPomA = thetaNow;
218  } else {
219  xPomB = xNow;
220  tPomB = tNow;
221  thetaPomB = thetaNow;
222  }
223 
224  // Done.
225  return true;
226 
227 }
228 
229 //--------------------------------------------------------------------------
230 
231 // Return x*f_P/p( x), ie. Pomeron flux inside proton, integrated over t.
232 
233 double HardDiffraction::xfPom(double xIn) {
234 
235  // Setup t range.
236  pair<double, double> tLim = tRange(xIn);
237  double tMin = tLim.first;
238  double tMax = tLim.second;
239  double x = xIn;
240  double xFlux = 0.;
241  if (tMin > 0. || tMax > 0.) return 0.;
242 
243  // Schuler-Sjostrand Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
244  // flux = normPom * 1/x * exp(2t(2.3 + 0.25 * log(1/x)))
245  // => x * flux = normPom * exp(2t(2.3 + 0.25*log(1/x)))
246  if (pomFlux == 1) {
247  double b = b0 + ap * log(1./x);
248  xFlux = normPom/(2.*b) * ( exp(2.*b*tMax) - exp(2.*b*tMin));
249  }
250 
251  // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
252  // flux = normPom * (1/x) * (6.38 *exp(8*t)+ 0.424 * exp(3*t))
253  // => x * flux = normPom * (6.38 *exp(8*t)+ 0.424 * exp(3*t))
254  else if (pomFlux == 2) {
255  xFlux = normPom * (A1/a1 * (exp(a1*tMax) - exp(a1*tMin))
256  + A2/a2 * (exp(a2*tMax) - exp(a2*tMin)));
257  }
258 
259  // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
260  // flux = normPom * x^(1 - 2*alpha(t)) * exp(-R_N^2 * t)
261  // => x * flux = normPom * x^(2 - 2*alpha(t)) * exp(-R_N^2 * t)
262  else if (pomFlux == 3) {
263  double b = (a1 + 2. * ap * log(1./x));
264  xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.));
265  xFlux *= (exp(b*tMax) - exp(b*tMin))/b;
266  }
267 
268  // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
269  // flux = 9 beta^2(0)/(4 pi^2) x^(1 - 2*alpha(t)) F_1(t)^2 with
270  // F_1(t)^2 = 0.27 exp(8.38 t) + 0.56 exp(3.78 t) + 0.18 exp(1.36 t)
271  // = (4m_p^2-2.8t)^2/(4m_p^2-t)^2*(1/(1-t/0.7))^4
272  // => x * flux = 9 beta^2(0)/(4 pi^2) * x^(2 - 2*\alpha(t)) F_1(t)^2
273  else if (pomFlux == 4) {
274  double Q = 2. * ap * log(1./x);
275  xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.));
276  xFlux *= (A1/(Q + a1) * (exp((Q + a1)*tMax) - exp((Q + a1)*tMin))
277  + A2/(Q + a2) * (exp((Q + a2)*tMax) - exp((Q + a2)*tMin))
278  + A3/(Q + a3) * (exp((Q + a3)*tMax) - exp((Q + a3)*tMin)));
279  }
280 
281  // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
282  // flux = normPom * F_1(t)^2 * exp((2 * alpha(t) - 1)*log(1/x))
283  // F_1(t)^2 = 0.9 exp(4.6 t) + 0.1 exp(0.6 t)
284  // = (4m_p^2-2.8t)^2/(4m_p^2-t)^2*(1/(1-t/0.7))^4
285  // => x * flux = normPom * F_1(t)^2 * exp( 2*(alpha(t) -1)*log(1/x))
286  else if (pomFlux == 5) {
287  double Q = 2. * ap * log(1./x);
288  xFlux = normPom * exp(log(1./x) * ( 2.*a0 - 2.));
289  xFlux *= (A1/(Q + a1) * (exp((Q + a1)*tMax) - exp((Q + a1)*tMin))
290  + A2/(Q + a2) * (exp((Q + a2)*tMax) - exp((Q + a2)*tMin)));
291  }
292 
293  // H1 Fit A, B Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
294  // flux = normPom * exp(B_Pom*t)/x^(2*\alpha(t)-1)
295  // => x * flux = normPom * exp(B_Pom * t) / x^(2*\alpha(t)-2)
296  else if (pomFlux == 6 || pomFlux == 7) {
297  double b = b0 + 2. * ap * log(1./x);
298  xFlux = normPom * exp(log(1./x) * ( 2.*a0 - 2.));
299  xFlux *= (exp(b*tMax) - exp(b*tMin))/b;
300  }
301 
302  // Done
303  if (usePomInPhoton) return xFlux * rescale * sigTotRatio;
304  else return xFlux * rescale;
305 
306 }
307 
308 //--------------------------------------------------------------------------
309 
310 // Pick a t value according to different Pomeron flux parametrizations.
311 
312 double HardDiffraction::pickTNow(double xIn) {
313 
314  // Get kinematical limits for t. initial values.
315  pair<double, double> tLim = HardDiffraction::tRange(xIn);
316  double tMin = tLim.first;
317  double tMax = tLim.second;
318  double tTmp = 0.;
319  double rndm = rndmPtr->flat();
320 
321  // Schuler-Sjostrand Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
322  if (pomFlux == 1) {
323  double b = b0 + ap * log(1./xIn);
324  tTmp = log( rndm*exp(2.*b*tMin) + (1. - rndm)*exp(2.*b*tMax))/(2.*b);
325  }
326 
327  // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
328  else if (pomFlux == 2) {
329  double prob1 = A1/a1 * (exp(a1*tMax) - exp(a1*tMin));
330  double prob2 = A2/a2 * (exp(a2*tMax) - exp(a2*tMin));
331  prob1 /= (prob1 + prob2);
332  tTmp = (prob1 > rndmPtr->flat())
333  ? log( rndm * exp(a1*tMin) + (1. - rndm) * exp(a1*tMax))/a1
334  : log( rndm * exp(a2*tMin) + (1. - rndm) * exp(a2*tMax))/a2;
335  }
336 
337  // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
338  else if (pomFlux == 3) {
339  double b = (2. * ap * log(1./xIn) + a1);
340  tTmp = log( rndm * exp(b*tMin) + (1. - rndm) * exp(b*tMax))/b;
341  }
342 
343  // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
344  else if (pomFlux == 4) {
345  double b1 = 2. * ap * log(1./xIn) + a1;
346  double b2 = 2. * ap * log(1./xIn) + a2;
347  double b3 = 2. * ap * log(1./xIn) + a3;
348  double prob1 = A1/b1 * ( exp(b1*tMax) - exp(b1*tMin));
349  double prob2 = A2/b2 * ( exp(b2*tMax) - exp(b2*tMin));
350  double prob3 = A3/b3 * ( exp(b3*tMax) - exp(b3*tMin));
351  double rndmProb = (prob1 + prob2 + prob3) * rndmPtr->flat() ;
352  if (rndmProb < prob1)
353  tTmp = log( rndm * exp(b1*tMin) + (1. - rndm) * exp(b1*tMax))/b1;
354  else if ( rndmProb < prob1 + prob2)
355  tTmp = log( rndm * exp(b2*tMin) + (1. - rndm) * exp(b2*tMax))/b2;
356  else
357  tTmp = log( rndm * exp(b3*tMin) + (1. - rndm) * exp(b3*tMax))/b3;
358  }
359 
360  // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
361  else if (pomFlux == 5) {
362  double b1 = a1 + 2. * ap * log(1./xIn);
363  double b2 = a2 + 2. * ap * log(1./xIn);
364  double prob1 = A1/b1 * (exp(b1*tMax) - exp(b1*tMin));
365  double prob2 = A2/b2 * (exp(b2*tMax) - exp(b2*tMin));
366  prob1 /= (prob1 + prob2);
367  tTmp = (prob1 > rndmPtr->flat())
368  ? log( rndm * exp(b1*tMin) + (1. - rndm) * exp(b1*tMax))/b1
369  : log( rndm * exp(b2*tMin) + (1. - rndm) * exp(b2*tMax))/b2;
370  }
371 
372  // H1 Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
373  else if (pomFlux == 6 || pomFlux == 7){
374  double b = b0 + 2. * ap * log(1./xIn);
375  tTmp = log( rndm * exp(b*tMin) + (1. - rndm) * exp(b*tMax))/b;
376  }
377 
378  // Done.
379  return tTmp;
380 
381 }
382 
383 //--------------------------------------------------------------------------
384 
385 // Return x*f_P/p( x, t), ie. Pomeron flux inside proton differential in t.
386 
387 double HardDiffraction::xfPomWithT(double xIn, double tIn) {
388 
389  // Initial values.
390  double x = xIn;
391  double t = tIn;
392  double xFlux = 0.;
393 
394  // Schuler-Sjostrand Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
395  if (pomFlux == 1) {
396  double b = b0 + ap * log(1./x);
397  xFlux = normPom * exp( 2.*b*t);
398  }
399 
400  // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
401  else if (pomFlux == 2)
402  xFlux = normPom * (A1 * exp(a1*t) + A2 * exp(a2*t));
403 
404  // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
405  else if (pomFlux == 3) {
406  xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.))
407  * exp(t * (a1 + 2.*ap*log(1./x)));
408  }
409 
410  // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
411  else if (pomFlux == 4){
412  double sqrF1 = A1 * exp(a1*t) + A2 * exp(a2*t) + A3 * exp(a3*t);
413  xFlux = normPom * pow(x, 2. + 2. * (a0 + ap*t)) * sqrF1;
414  }
415 
416  // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
417  else if (pomFlux == 5) {
418  double sqrF1 = A1 * exp(a1*t) + A2 * exp(a2*t);
419  xFlux = normPom * sqrF1 * exp(log(1./x) * (-2. + a0 + ap*t));
420  }
421 
422  // H1 Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
423  else if (pomFlux == 6 || pomFlux == 7)
424  xFlux = normPom * exp(b0*t)/pow(x, 2. * (a0 + ap*t) - 2.);
425 
426  // Done
427  if (usePomInPhoton) return xFlux * rescale * sigTotRatio;
428  else return xFlux * rescale;
429 
430 }
431 
432 //--------------------------------------------------------------------------
433 
434 // Set up t range. See p. 113 of 6.4 manual.
435 
436 pair<double, double> HardDiffraction::tRange(double xIn) {
437 
438  // Set up diffractive masses.
439  // s1 = mA^2, s2 = mB^2,
440  // s3 = M^2 (= mA^2 if A scatteres elastically)
441  // s4 = M^2 (= mB^2 if B scatteres elastically)
442  double eCM = infoPtr->eCM();
443  s = eCM * eCM;
444  double M2 = xIn * s;
445  s1 = pow2(mA);
446  s2 = pow2(mB);
447  s3 = (iBeam == 1) ? s1 : M2;
448  s4 = (iBeam == 2) ? s2 : M2;
449 
450  // Error exit if too large diffractive mass.
451  if (sqrt(s3) + sqrt(s4) >= eCM) return make_pair( 1., 1.);
452 
453  // Calculate kinematics.
454  double lambda12 = sqrtpos(pow2(s - s1 - s2) - 4. * s1 * s2);
455  double lambda34 = sqrtpos(pow2(s - s3 - s4) - 4. * s3 * s4);
456  double tmp1 = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
457  double tmp2 = lambda12 * lambda34 / s;
458  double tmp3 = (s1 + s4 - s2 - s3) * (s1 * s4 - s2 * s3) / s
459  + (s3 - s1) * (s4 - s2);
460  double tMin = -0.5 * (tmp1 + tmp2);
461  double tMax = tmp3 / tMin;
462 
463  // Done.
464  return make_pair(tMin, tMax);
465 
466 }
467 
468 //--------------------------------------------------------------------------
469 
470 // Get the scattering angle from the chosen t.
471 
472 double HardDiffraction::getThetaNow( double xIn, double tIn) {
473 
474  // Set up diffractive masses.
475  // s1 = mA^2, s2 = mB^2,
476  // s3 = M^2 (= mA^2 if A scatteres elastically)
477  // s4 = M^2 (= mB^2 if B scatteres elastically)
478  double eCM = infoPtr->eCM();
479  s = eCM * eCM;
480  double M2 = xIn * s;
481  s1 = pow2(mA);
482  s2 = pow2(mB);
483  s3 = (iBeam == 1) ? s1 : M2;
484  s4 = (iBeam == 2) ? s2 : M2;
485 
486  // Find theta from the chosen t.
487  double lambda12 = sqrtpos(pow2(s - s1 - s2) - 4. * s1 * s2);
488  double lambda34 = sqrtpos(pow2(s - s3 - s4) - 4. * s3 * s4);
489  double tmp1 = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4)/s;
490  double tmp2 = lambda12 * lambda34 / s;
491  double tmp3 = (s1 + s4 - s2 - s3) * (s1 * s4 - s2 * s3) / s
492  + (s3 - s1) * (s4 - s2);
493  double cosTheta = min(1., max(-1., (tmp1 + 2. * tIn) / tmp2));
494  double sinTheta = 2. * sqrtpos( -(tmp3 + tmp1 * tIn + tIn * tIn) ) / tmp2;
495  double theta = asin( min(1., sinTheta));
496 
497  if (cosTheta < 0.) theta = M_PI - theta;
498 
499  // Done.
500  return theta;
501 
502 }
503 
504 //==========================================================================
505 
506 } // end namespace Pythia8
Generated by   doxygen 1.8.5