g2t kinematics check

I was thinking that quite a bit of the info that we store in e.g. StPythiaEvent is redundant. I wrote a class that stored a much smaller set of persistent information (parton 4-vectors) and calculated all other derived quantities (s, t, u, cosTheta, …) on-the-fly.

As a check I compared the values from the g2t table that were stored in StPythiaEvent with the values I calculated (using the 4-vectors in StPythiaEvent). My calculated versions of the Mandelstam variables and cosTheta were identical to the g2t values, but I got rather different values for x1, x2, and hard pT. Here’s a plot:


The fact that there are ~no counts in the top-right and bottom-left quadrants makes it look like either the reco or the g2t values are neglecting initial state radiation. That’s not supposed to be the case.

I’m calculating x1 and x2 as

x1 = (pT1*exp(eta1) + pT2*exp(eta2)) / 200
x2 = (pT1*exp(-eta1) + pT2*exp(-eta2)) / 200
If you know of something else PYTHIA is doing that I’m overlooking, please let me know.

Update 11 July 2008

I did a more systematic investigation into the formulas that are used to calculate the hard scattering kinematics stored in the g2t_pythia table, specifically mand_s, mand_t, mand_u, hard_p, cos_th, bjor_1, and bjor_2. I tried a few alternatives (taking parton masses into account, etc.) and plotted the fractional deviation between my reco quantities and the ones stored in the table. If the maximum fractional deviation was less than 1E-6, I took that difference to be a result of floating-point imprecision and the formula to be correct. Here’s what I came up with:

mand_s = (p1+p2)^2
mand_t = -0.5 * s_hat * (1-cos_theta)
mand_u = -0.5 * s_hat * (1+cos_theta)
hard_p = sqrt(Q2)
Q2 = t_hat * u_hat / s_hat

cos_th: polar angle of parton 3 after boosting to c.m. of hard scattering and rotating so that p1 and p2 are traveling along the z-axis

bjor_1, bjor_2: I don’t have a complete definition of these in terms of the parton 4-vectors. It’s true that x1*x2 = s_hat/s, where s is exactly 40000 GeV^2, so given x1 it’s possible to calculate x2. PYTHIA evaluates x1 and x2 by boosting to the c.m. frame of the hard scattering, rotating so that the collision occurs along z-axis, and then boosting back along z. The trouble is determining the value of this final boost (it’s not just the magnitude of the initial boost). From Section 9.2 of the PYTHIA manual:

Since the initial-state radiation machinery assigns space-like virtualities to the incoming partons, the definitions of x in terms of energy fractions and in terms of momentum fractions no longer coincide, and so the interacting subsystem may receive a net longitudinal boost compared with naïve expectations, as part of the parton-shower machinery.

For completeness, note that the following formulas will not always give results that match the g2t_pythia quantities exactly. Differences are usually small (less than 0.5%), but can be larger especially in the case of heavy quark production.

t_hat = (p1-p3)^2 = (p2-p4)^2
u_hat = (p1-p4)^2 = (p2-p3)^2
Q2 = (m3^2 +m4^2)/2 + (t_hat * u_hat - m3^2 * m4^2)/s_hat

and of course, any aforementioned expressions for x1 and x2 in terms of parton 4-vectors are also unreliable.