/*
   Subroutine: abc_mlt.c

      This subroutine performs multiple scattering on particle trajectories.
   First, the rms plane scattering angle (t0) is calculated from the particle
   momentum and material (i1) radiation length [see Particle Physics Booklet,
   Phys.Rev. D50,1173(1994)] for two orthogonal directions, from which a 
   random scatter angle is derived.  The angle is then rotated back to lab
   coordinates, where particle direction is thus altered.
*/

#include "abc.h"

void abc_mlt(int it,int i1,int i2)
{
   int i,ip,is;
   double mp,rl,qp,t0,tx,ty,ts,ps,pt,pp,pq;
   double bx,by,bz,bm,rx,ry,rz,rm,px,py,pz,pm,pxp,pyp,pzp;
   double the1,phi1,the2,phi2,thes,phis,thed,phid,rh,ru,rv;
   double r1[4],r2[4],p1[4],p2[4];

/* Fill particle data and point arrays */
   rl = 0.01*matrad[objmat[i1]]; if(rl <= 0.0) return;
   ip = trkprt[it]; mp = prtmas[ip]; qp = prtchr[ip];
   is = trkpnt[it] - 2;
   r1[0] = trktpo[it][is]; p1[0] = trkemo[it][is];
   r1[1] = trkxpo[it][is]; p1[1] = trkxmo[it][is];
   r1[2] = trkypo[it][is]; p1[2] = trkymo[it][is];
   r1[3] = trkzpo[it][is]; p1[3] = trkzmo[it][is];
   is = trkpnt[it] - 1;
   r2[0] = trktpo[it][is]; p2[0] = trkemo[it][is];
   r2[1] = trkxpo[it][is]; p2[1] = trkxmo[it][is];
   r2[2] = trkypo[it][is]; p2[2] = trkymo[it][is];
   r2[3] = trkzpo[it][is]; p2[3] = trkzmo[it][is];

/* Calculate rms plane scattering angle (t0) */
   if(p1[0] <= 0.0) { fprintf(stderr,"abc_mlt: E <= 0 \n"); return; }
   bx = p1[1]/p1[0]; by = p1[2]/p1[0]; bz = p1[3]/p1[0];
   bm = sqrt(bx*bx + by*by + bz*bz);
   rx = r2[1]-r1[1]; ry = r2[2]-r1[2]; rz = r2[3]-r1[3];
   rm = sqrt(rx*rx + ry*ry + rz*rz);
   px = p1[1]; py = p1[2]; pz = p1[3];
   pm = sqrt(px*px + py*py + pz*pz);
   t0 = ((0.0136*qp)/(bm*pm))*sqrt(rm/rl)*(1.0+0.038*log(rm/rl));

/* Choose random space scattering angles (thes,phis) */
   tx = uti_rnn(0.0,t0);
   ty = uti_rnn(0.0,t0);
   thes = sqrt(tx*tx + ty*ty);
   phis = -pi + uti_rnu()*2.0*pi;

/* Add scatter to momentum */
   px = p2[1]; py = p2[2]; pz = p2[3];
   pm = sqrt(px*px + py*py + pz*pz);
   if(pm <= 0.0) { fprintf(stderr,"abc_mlt: p <= 0 \n"); return; }
   /* Momentum spherical angles (the1,phi1) */
   the1 = acos(pz/pm);
   phi1 = atan2(py,px);
   /* New spherical angles (the1+thed,phi1+phid) */
   rh = pm*sin(thes);
   ru = rh*sin(phis); phid = asin(ru/pm);
   rv = rh*cos(phis); thed = asin(rv/pm);
   px = pm*sin(the1+thed)*cos(phi1+phid);
   py = pm*sin(the1+thed)*sin(phi1+phid);
   pz = pm*cos(the1+thed);
   p2[1] = px; p2[2] = py; p2[3] = pz;

/* Adjust final point momentum */
   is = trkpnt[it] - 1;
   trkxmo[it][is] = p2[1];
   trkymo[it][is] = p2[2];
   trkzmo[it][is] = p2[3];
}