/*
   Subroutine: abc_dec.c

      This subroutine generates particle decays according to the following
   input modes (id), equivalent to generator flag:

         id > +0  -->  force decay through mode id
         id = +0  -->  normal decay (via lifetime, random mode)
         id = -1  -->  force decay (random mode)
         id = -2  -->  suppress decay

   One particle (trivial case), two-particle, and three-particle decay modes
   are possible.
*/

#include "abc.h"

int abc_dec(int it,int id)
{
   int ichs=0,maxchs=100;
   int i,ir,is,ip,im,in,nm,np,nf,i1,i2,i3;
   double tlif,tdel,pdec,pran,pmod,psum,ptot;
   double r0[4], m[4],bx,by,bz,at,ap,ea,eb,ec,e1,e2,e3,r1,r2,r3;
   double p0[4],p1[4],p2[4],p3[4];
   double c0[4],c1[4],c2[4],c3[4];
   double d0[4],d1[4],d2[4],d3[4];

/* Fill zeroth array */
   if(id < -1) return(0);                                /* deactivate decay */
   ip = trkprt[it]; nm = prtmdn[ip]; m[0] = prtmas[ip];
   is = trkpnt[it] - 1;
   r0[0] = trktpo[it][is]; p0[0] = trkemo[it][is];
   r0[1] = trkxpo[it][is]; p0[1] = trkxmo[it][is];
   r0[2] = trkypo[it][is]; p0[2] = trkymo[it][is];
   r0[3] = trkzpo[it][is]; p0[3] = trkzmo[it][is];
   for(i=0;i<4;i++) { p1[i] = 0.0; p2[i] = 0.0; p3[i] = 0.0; }

/* Test for decay */
   if(id == 0 && is > 0)
   {
      ir = 0; tlif = 0.0; tdel = 0.0;
      tdel = trktpo[it][is] - trktpo[it][is-1];             /* lab ct2 - ct1 */
      if(m[0] > 0.0) tlif = prtlif[ip]*p0[0]/m[0];   /* lab life c*tau*gamma */
      if(tlif <= 0.0)
      {
         fprintf(stderr,"abc_dec: zero mass/life %d \n",ip); return(0);
      }
      pdec = exp(-tdel/tlif); pran = uti_rnu();       /* Poisson probability */
      if(pran > pdec) ir = 1;
      if(trctag == 10 && ip == 9)
      fprintf(stderr,"abc_dec: %d %lf %lf %lf %lf \n",ir,tlif,tdel,pdec,pran);
   }
   else ir = 1;
   if(ir == 0) return(0);

/* Choose decay mode */
   if(id > 0) im = id;
   else
   {
      ptot = 0.0; for(i=1;i<=nm;i++) ptot += prtmdb[ip][i];
      im = 0; pmod = ptot*uti_rnu();
      psum = 0.0; for(i=1;i<=nm;i++)
      {
         if(pmod >= psum) im = i; psum = psum + prtmdb[ip][i];
      }
      if(im < 1 || im > nm)
      {
         fprintf(stderr,"  Error: no decay mode chosen \n"); return(0);
      }
   }

/* Determine final state particles and cf beta */
   np = prtmdp[ip][im]; nf = 0; i1 = 0; i2 = 0; i3 = 0;
   if(np > 9999) { nf++; i3 = np/10000; np -= i3*10000; m[3] = prtmas[i3]; }
   if(np >   99) { nf++; i2 = np/100;   np -= i2*100;   m[2] = prtmas[i2]; }
   if(np >    0) { nf++; i1 = np/1;     np -= i1*1;     m[1] = prtmas[i1]; }
   bx = p0[1]/p0[0]; by = p0[2]/p0[0]; bz = p0[3]/p0[0];
   uti_bst(+bx,+by,+bz,p0,c0);

   if(trctag == 10)
   {
      fprintf(stderr,"abc_dec: %d %d %d %d %d \n",im,nf,i1,i2,i3);
      fprintf(stderr,"         %lf %lf %lf %lf \n",p0[0],p0[1],p0[2],p0[3]);
   }

/* Generate 1-body final state (non-random) */
   if(nf == 1)
   {
      for(i=0;i<4;i++) p1[i] = p0[i];
   }

/* Generate 2-body final state */
   if(nf == 2)
   {
   /* Particle 1 */
      at = acos(-1.0 + uti_rnu()*2.0);
      ap = -pi + uti_rnu()*2.0*pi;
      c1[0] = (c0[0]*c0[0] + m[1]*m[1] - m[2]*m[2])/(2.0*c0[0]);
      c1[1] = sqrt(c1[0]*c1[0] - m[1]*m[1])*sin(at)*cos(ap);
      c1[2] = sqrt(c1[0]*c1[0] - m[1]*m[1])*sin(at)*sin(ap);
      c1[3] = sqrt(c1[0]*c1[0] - m[1]*m[1])*cos(at);
      uti_bst(-bx,-by,-bz,c1,p1);
   /* Particle 2 */
      at = pi - at;
      ap = ap - pi; if(ap < -pi) ap = ap + 2.0*pi;
      c2[0] = (c0[0]*c0[0] + m[2]*m[2] - m[1]*m[1])/(2.0*c0[0]);
      c2[1] = sqrt(c2[0]*c2[0] - m[2]*m[2])*sin(at)*cos(ap);
      c2[2] = sqrt(c2[0]*c2[0] - m[2]*m[2])*sin(at)*sin(ap);
      c2[3] = sqrt(c2[0]*c2[0] - m[2]*m[2])*cos(at);
      uti_bst(-bx,-by,-bz,c2,p2);
   }

/* Generate 3-body final state */
   if(nf == 3)
   {
   /* Choose particle cm energies */
      choose:
      ichs++; if(ichs > maxchs)
      {
         fprintf(stderr,"abc_dec: %d exceeded \n",maxchs); return(0);
      }
      e1 = m[1] + uti_rnu()*(c0[0]-m[2]-m[3]);
      e2 = m[2] + uti_rnu()*(c0[0]-m[1]-m[3]);
      e3 = c0[0] - e1 - e2;
      if(e3 < m[3]) goto choose;
      r1 = sqrt(e1*e1 - m[1]*m[1]);
      r2 = sqrt(e2*e2 - m[2]*m[2]);
      r3 = sqrt(e3*e3 - m[3]*m[3]);
      if(r1 > (r2+r3)) goto choose;
      if(r2 > (r3+r1)) goto choose;
      if(r3 > (r1+r2)) goto choose;
   /* Calculate production plane momentum */
      d1[1] = r1;
      d2[1] = (r3*r3 - r1*r1 - r2*r2)/(2.0*r1);
      d3[1] = (r2*r2 - r3*r3 - r1*r1)/(2.0*r1);
      d1[2] = 0.0;
      d2[2] = +sqrt(r2*r2 - d2[1]*d2[1]);
      d3[2] = -sqrt(r3*r3 - d3[1]*d3[1]);
      d1[3] = 0.0;
      d2[3] = 0.0;
      d3[3] = 0.0;
   /* Choose production plane angles, rotate to cm, boost to lab */
      ea = -180.0 + uti_rnu()*360.0;
      eb = (180.0/pi)*acos(-1.0 + uti_rnu()*2.0);
      ec = -180.0 + uti_rnu()*360.0;
      uti_rot(+1,d1[1],d1[2],d1[3],0.,0.,0.,ea,eb,ec,&c1[1],&c1[2],&c1[3]);
      uti_rot(+1,d2[1],d2[2],d2[3],0.,0.,0.,ea,eb,ec,&c2[1],&c2[2],&c2[3]);
      uti_rot(+1,d3[1],d3[2],d3[3],0.,0.,0.,ea,eb,ec,&c3[1],&c3[2],&c3[3]);
      c1[0] = e1; uti_bst(-bx,-by,-bz,c1,p1);
      c2[0] = e2; uti_bst(-bx,-by,-bz,c2,p2);
      c3[0] = e3; uti_bst(-bx,-by,-bz,c3,p3);
   }

/* Fill new track arrays */
   if(nf > 0)
   {
      trknum++; in = trknum - 1;
      trktpo[in][0] = r0[0]; trkemo[in][0] = p1[0]; trkpnt[in] = 1;
      trkxpo[in][0] = r0[1]; trkxmo[in][0] = p1[1]; trkprt[in] = i1;
      trkypo[in][0] = r0[2]; trkymo[in][0] = p1[2];
      trkzpo[in][0] = r0[3]; trkzmo[in][0] = p1[3];
   }
   if(nf > 1)
   {
      trknum++; in = trknum - 1;
      trktpo[in][0] = r0[0]; trkemo[in][0] = p2[0]; trkpnt[in] = 1;
      trkxpo[in][0] = r0[1]; trkxmo[in][0] = p2[1]; trkprt[in] = i2;
      trkypo[in][0] = r0[2]; trkymo[in][0] = p2[2];
      trkzpo[in][0] = r0[3]; trkzmo[in][0] = p2[3];
   }
   if(nf > 2)
   {
     trknum++; in = trknum - 1;
      trktpo[in][0] = r0[0]; trkemo[in][0] = p3[0]; trkpnt[in] = 1;
      trkxpo[in][0] = r0[1]; trkxmo[in][0] = p3[1]; trkprt[in] = i3;
      trkypo[in][0] = r0[2]; trkymo[in][0] = p3[2];
      trkzpo[in][0] = r0[3]; trkzmo[in][0] = p3[3];
   }

   if(trctag == 10 && ip == 9)
   {
      fprintf(stderr,"         %lf %lf %lf %lf \n",p1[0],p1[1],p1[2],p1[3]);
      fprintf(stderr,"         %lf %lf %lf %lf \n",p2[0],p2[1],p2[2],p2[3]);
      fprintf(stderr,"         %lf %lf %lf %lf \n",p3[0],p3[1],p2[2],p2[3]);
   }
   return(1);
}