/*
   Subroutine: abc_los.c

      This subroutine performs energy loss for particles traversing a material
   defined for object (i1).  The energy loss per distance (dE/dx) is estimated
   by the routine los_bbl, using the Bethe-Bloch equation with optional density
   correction.  The energy loss is then integrated (with step size dx set by
   inclos) over the total track step using a fourth-order Runge-Kutta scheme,
   and the final energy and momentum magnitude are adjusted accordingly.
*/

#include "abc.h"

void abc_los(int it,int i1,int i2)
{
   int i,io,im,ip,is,ic=2;
   double am,an,dn,mp,rl,qp;
   double rs,rm,pm,pe,dx,dy,dz,dr,de,e1,e2;
   double r1[4],r2[4],p1[4],p2[4];

/* Fill material data */
   io = i1; im = objmat[io];
   am = matatm[im]; an = matatn[im];
   dn = matden[im]; rl = matrad[im];
   if(am <= 0.0 || an <= 0.0 || dn <= 0.0 || rl <= 0.0) return;

/* Fill particle data and point arrays */
   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];

/* Determine step length, initial momentum, and direction */
   rm = sqrt((r2[1]-r1[1])*(r2[1]-r1[1]) + (r2[2]-r1[2])*(r2[2]-r1[2])
           + (r2[3]-r1[3])*(r2[3]-r1[3]));
   pm = sqrt((p2[1]*p2[1])+(p2[2]*p2[2])+(p2[3]*p2[3]));
   if(pm < minmom) return;
   dx = p2[1]/pm; dy = p2[2]/pm; dz = p2[3]/pm;

/* Integrate dE/dx */
   rs = 0.0; dr = inclos*rm;
   while(rs < rm)
   {
      /* Determine initial E,dE/dx and perform step */
      e1 = sqrt(pm*pm + mp*mp);
      de = los_bbl(ic,io,ip,e1);
      e2 = los_rk4(ic,io,ip,e1,de,rs,dr);
      pe = e2;
      /* Check for range out */
      if((pe*pe-mp*mp) < minmom*minmom)
      {
         pm = minmom;
         pe = sqrt(pm*pm + mp*mp);
         rm = rs;
      }
      /* Decrement momentum and increment step */
      else pm = sqrt(pe*pe - mp*mp);
      rs += dr;
   }

/* Determine new position,momentum */
   p2[0] = pe;    r2[0] = r1[0] + rm*pe/pm;
   p2[1] = pm*dx; r2[1] = r1[1] + rm*dx;
   p2[2] = pm*dy; r2[2] = r1[2] + rm*dy;
   p2[3] = pm*dz; r2[3] = r1[3] + rm*dz;

/* Adjust final point position,momentum */
   is = trkpnt[it] - 1;
   if(i1 == i2)
   {
      trktpo[it][is] = r2[0];
      trkxpo[it][is] = r2[1];
      trkypo[it][is] = r2[2];
      trkzpo[it][is] = r2[3];
   }
   trkemo[it][is] = p2[0];
   trkxmo[it][is] = p2[1];
   trkymo[it][is] = p2[2];
   trkzmo[it][is] = p2[3];
}

/*
   Subroutine: los_bbl.c

      This subroutine returns the energy loss dE/dx (GeV/m) for particle (ip)
   with energy (ep) in object (io), with the possible corrections:

         ic = 0  -->  simple 1/beta**2 dependence
         ic = 1  -->  Bethe-Bloch
         ic = 2  -->  Bethe-Bloch with density correction

   Formulas for the Bethe-Bloch and associated approximations are from
   [Leo,Techniques for Nuc. and Part. Phys. Exp.] and the Particle Physics
   Booklet [Phys.Rev. D54,22.1(1996)].
*/
double los_bbl(ic,io,ip,ep)
int ic,io,ip;
double ep;
{
   int im;
   double am,an,dn,rl,mp,qp;
   double b1,g1,b2,g2,e1,e2,m1,m2,mr,ti,tm,hw;
   double ci=0.0,cd=0.0,dedx=0.0;

/* Fill material data */
   im = objmat[io];
   am = matatm[im]; an = matatn[im];
   dn = matden[im]; rl = matrad[im];
   qp = prtchr[ip]; mp = prtmas[ip];
   if(am <= 0.0 || an <= 0.0 || dn <= 0.0 || rl <= 0.0) return(dedx);
   if(qp == 0.0 || mp <= 0.0 || mp >= ep) return(dedx);
/* Calculate beta,gamma */
   m1 = 1000.0*mp; m2 = m1*m1; mr = me/m1;   /* convert to MeV */
   e1 = 1000.0*ep; e2 = e1*e1;
   b2 = (e2-m2)/e2; b1 = sqrt(b2);
   g2 = e2/m2;      g1 = sqrt(g2);
/* Determine corrections */
   if(ic == 0)
   {
      ci = 7.0;
   }
   if(ic >= 1)
   {
      if(an < 13) ti = 12.0*an + 7.0;
      else        ti = 9.76*an + 58.8*pow(an,-0.19);
      ti = 0.000001*ti;
      tm = (2.0*me*b2*g2)/(1.0 + 2.0*g1*mr + mr*mr);
      ci = 0.5*log((2.0*me*b2*g2*tm)/(ti*ti));
   }
   if(ic >= 2)
   {
      hw = 0.000028816*sqrt(dn*an/am);
      cd = log(hw/ti) + log(b1*g1) - 0.5;
      cd = log(b1*g1) - 0.5;
   }
/* Calculate dE/dx */
   dedx = -0.1*kb*(dn*an/am)*(qp*qp/b2)*(ci - b2 - 0.5*cd);
   if(dedx > 0.0) dedx = -1000.0;
   return(dedx);
}


/*
   Subroutine: los_rk4.c

      This subroutine performs a single fourth-order Runge-Kutta step for a
   particle with energy (ep) and dE/dx (dedx) over a step of length (h) [see
   Numerical Recipes,(16.1)], returning the decremented energy.
*/
double los_rk4(ic,io,ip,ep,dedx,x,h)
int ic,io,ip;
double ep,dedx,x,h;
{
   double efin;
   double h6,hh,xh,dem,det,et;
   hh = h/2.0;
   h6 = h/6.0;
   xh = x + hh;

   et = ep + hh*dedx;
   det = los_bbl(ic,io,ip,et);
   et = ep + hh*det;
   dem = los_bbl(ic,io,ip,et);
   et = ep + h*dem;
   dem = det + dem;
   det = los_bbl(ic,io,ip,et);

   efin = ep + h6*(dedx + det + 2.0*dem);
   return(efin);
}