~blok/ANALYSIS/kin.definition


DEFINITION OF KINEMATICAL VARIABLES FOR COINCIDENCE STUDIES
-----------------------------------------------------------

First all vectors are calculated in the laboratory system, which has
z along the beam, x to the right, y down. 
The spectrometers deliver the absolute value of their momentum p, and
two (Cartesian) angles, theta (always larger than zero), which is the angle
of the projection of p on to the xz plane with the (downstream) beam,
and x'_tar, which is the angle of p with respect to the xz plane.
Because the y-axis points downwards, x'_tar>0 means p below the xz plane,
i.e., py>0.
Since SOS is on the negative x side, in calculating e.g. q, theta_SOS
should be taken negative. Thus the formula's to calculate q, ph (the
hadron momentum) and pm (the missing momentum) are:
 
      cqx =        - p_e*cos(xp_e_tar)*sin(theta_e)
      cqy =        - p_e*sin(xp_e_tar)
      cqz = gpbeam - p_e*cos(xp_e_tar)*cos(theta_e)
      
      phx =  p_h*cos(xp_h_tar)*sin(theta_h)
      phy =  p_h*sin(xp_h_tar)
      phz =  p_h*cos(xp_h_tar)*cos(theta_h)

      pmx = cqx - phx 
      pmy = cqy - phy
      pmz = cqz - phz

      pm_par =(pmx*cqx+pmz*cqz)/cqxzabs   with  cqxzabs=sqrt(cqx**2+cqz**2)
      pm_oop = pmy
      pm_perp=(pmz*cqx-pmx*cqz)/cqxzabs 
      if (theta_e.gt.0)pm_perp=-pm_perp

Since we want pm_perp > 0 when pm closer to downstream beam than q,
we need this sign in case e' is to the right, and thus q to the left,
of the beam.
Effectively the sign introduced by the if statement is qx/abs(qx).
In this way pm_perp equals pm_x' (the primed coordinate system,
which is used to calculate phi_pi, will be defined later) in case q
lies in the x-z plane (qy=0).
      
Now we want to calculate the angles theta_pq and phi_pq, which are the
spherical angles in a coordinate system that has the z axis parallel to q
and the x axis inside the q-z-plane, i.e.  z' = q, y' = q X z, x' = y' X q.
This gives phi the way it is usually defined, i.e. phi=0 for in-plane
particles closer to the downstream beamline than q, phi=90 is above the
horizontal plane when q points to the right, and below the horizontal
plane when q points to the left.

The new unit vectors are then given by:

      dummy=sqrt( (cqx**2 + cqy**2) * cqabs**2 )
      new_x_x = -cqx*cqz/dummy
      new_x_y = -cqy*cqz/dummy
      new_x_z = (cqx**2 + cqy**2)/dummy

      dummy   = sqrt(cqx**2 + cqy**2)
      new_y_x =  cqy/dummy
      new_y_y = -cqx/dummy
      new_y_z =  0.

      new_z_x = cqx/cqabs
      new_z_y = cqy/cqabs
      new_z_z = cqz/cqabs

              
The inner product of q and ph yields directly the angle theta_pq:
   
      dot = phx*cqx + phy*cqy + phz*cqz
      cthetapq = acos(dot/(cqabs*ph))

For calculating phi we need the full vector ph in the new system:
      p_new_x = phx*new_x_x + phy*new_x_y + phz*new_x_z
      p_new_y = phx*new_y_x + phy*new_y_y + phz*new_y_z
      p_new_z = phx*new_z_x + phy*new_z_y + phz*new_z_z

      cphipi = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
      if (p_new_y.lt.0.)  cphipi = 2*c_pi - cphipi

It is interesting to see now that the x component of pm in the new
system, pm_new_x, when q is 'in-plane (qy=0), is 

      pm_new_x = (-pmx*qz+pmz*qx)/cqxzabs*(cqx/abs(cqx))

which is just the formula for pm_perp with the inclusion of the sign
as given by the if clause.