Possible Bug in ParticlePropagator.cc (3.1.2)

[Keith Pedersen]

Hello,

I have been looking through your modules all day to see how they work, and I found a potential bug in ParticlePropagator.cc. Allow me to quote from lines [140, 162].

140  if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
     {
142   // solve pt2*t^2 + 2*(px*x + py*y)*t + (fRadius2 - x*x - y*y) = 0
143   tmp = px*y - py*x;
144   discr2 = pt2*fRadius2 - tmp*tmp;

146   if(discr2 < 0)
      {
        // no solutions
        continue;
      }

152   tmp = px*x + py*y;
      discr = TMath::Sqrt(discr2);
      t1 = (-tmp + discr)/pt2;
      t2 = (-tmp - discr)/pt2;
156   t = (t1 < 0) ? t2 : t1;

158   z_t = z + pz*t;
      if(TMath::Abs(z_t) > fHalfLength)
      {
        t3 = (+fHalfLength - z) / pz;
        t4 = (-fHalfLength - z) / pz;
        t = (t3 < 0) ? t4 : t3;
      }

      x_t = x + px*t;
      y_t = y + py*t;
168   z_t = z + pz*t;

So we're dealing with particles travelling in a straight line. We want to use the vector equation (x2vec + p2vec*t)^2 = R2vec^2 (where all three vectors are in the transverse plane, R2vec is a spatial vector pointing from the beamline to the point where the particle leaves the cylinder, x2Vec is the particle's creation position, and p2Vec is its 2-momentum). This equation is clearly the logic behind line 142, although the sign is flipped in the C coefficient (A*t^2+B*t+C==0) of line 142.

What I don't understand is line 143 to 150. My interpretation of these lines is that discr2 = (B^2 - 4A*C)/4. If discr2 is negative, there are no real solutions. However, line 143 is not B, it is -Lz, which I do not understand. Additionally, even if line 143 were correct, tmp*tmp has the wrong sign in line 144 (also, while -pt2*(x^2 + y^2) is missing from line 144, this seems like a decent approximation for large R). Now, while line 152 uses the correct B, the sqrt still comes from the incorrect discr2.

Also, the logic in line 156 is a bit strange. If t1 is negative, then t2 will be MORE negative. This could lead to line 158 assigning the wrong z_t, producing an incorrect endcap punch-through test. Now, I know that t is not really time (t = c*time/E), but since c and E are manifestly positive, t < 0 corresponds to time < 0, which is always a bad solution.

In my interpretation, with the correct implementation of the quadratic equation, complex t are IMPOSSIBLE as long as pt > 0, so the sanity test really should be re-phrased. Furthermore, t should ALWAYS be assigned to the (+) solution in line 156 (t1), as |B|<sqrt(B^2 + pt2*(R^2 - x^2 - y^2)).

Am I completely missing the point here, or is this a bug that could be affecting neutral particles? Thanks for you help.

Keith Pedersen

[Pavel Demin]

Dear Keith,

Thank you for checking this part of the particle propagation code!

I agree that the equation in the commented line 142 is wrong. It should indeed read:

// solve pt2*t^2 + 2*(px*x + py*y)*t + (x*x + y*y - fRadius2) = 0

The lines 143 and 144 look correct. Here is how they are obtained:

a = pt*pt = (px*px + py*py)
b = 2*(px*x + py*y)
c = (x*x + y*y - R*R)

D*D/4 = b*b/4 - a*c = (px*x + py*y)*(px*x + py*y) - (px*px + py*py)*(x*x + y*y - R*R)

(px*x + py*y)*(px*x + py*y) = px*px*x*x + py*py*y*y + 2*px*py*x*y

(px*px + py*py)*(x*x + y*y - R*R) = px*px*x*x + px*px*y*y - px*px*R*R + py*py*x*x + py*py*y*y - py*py*R*R

D*D/4 = 2*px*py*x*y + px*px*R*R - px*px*y*y - py*py*x*x + py*py*R*R

D*D/4 = (px*px + py*py)*R*R - (px*px*y*y + py*py*x*x - 2*px*py*x*y)

D*D/4 = pt*pt*R*R - (px*y - py*x)*(px*y - py*x)

So, the D/4, t1 and t2 are all correct.

However, I agree that there is something wrong with the line 156. It should read:

t = (t2 < 0) ? t1 : t2;
if(t < 0)
{
  // no solutions
  continue;
}

Moreover, I'd say that the line 163 should be fixed in the same way:

t = (t4 < 0) ? t3 : t4;
if(t < 0)
{
  // no solutions
  continue;
}

Regards,

Pavel

[Pavel Demin]

After some thought, I'd say that the current version of the particle propagation code is OK and the corrections that I proposed in my previous comment aren't needed.

In the following I'll try to explain why.

The following expression on lines 123-126, filters out all the particles outside the the cylinder:

if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
{
  continue;
}

So, discr2 can't be negative.

Taking into account that

tmp = px*x + py*y

discr = sqrt(tmp*tmp - a*c)

a*c = pt*pt*(x*x + y*y - fRadius2) < 0

discr*discr >= tmp*tmp

t1*t2 = (-tmp + discr)*(-tmp - discr)/(pt*pt) = (tmp*tmp - discr*discr)/(pt*pt) <= 0

I can say that if t1 is not zero then t1 and t2 have always opposite signs.

So, the following two expressions are equivalent:

t = (t1 < 0) ? t2 : t1;

t = (t2 < 0) ? t1 : t2;

Similarly, we can see that t3 and t4 have opposite signs if they are not zero:

t3*t4 = (z*z - fHalfLength*fHalfLength)/(pz*pz) < 0

So, there is nothing to fix after all.

[Keith Pedersen]

Pavel,

Yes, I see that I didn't follow the math all the way to the end. Sorry to bother you.

Keith