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