void elasticscatter(int ievent,int *kf){
  int alpha,i,j;
  double tevent,delt;
  double ptot[4],q[4],q0[4],b[4],sthet,cthet,theta;
  double ptot2,ptotdotq0,q02,q2,m2j,m2i,b2;
#ifdef EXTRA_TESTS
  double m2itest,m2jtest;
#endif
  double q0dotb,ptotdotb,normq,normb;
#ifdef EXTRA_TESTS
  double qdotbtest,pdotbtest;
#endif
  double delt_i,delt_j;
  double *p_i,*p_j,*r_i,*r_j;
  /* kf[0] is the number of outgoing particles */

  i=eventinfo[ievent].i;
  j=eventinfo[ievent].j;
  kf[0]=2;
  kf[1]=i;
  kf[2]=j;
  delt=eventinfo[ievent].delt;
  tevent=eventinfo[ievent].time;
  r_i=partinfo[i].r;
  r_j=partinfo[j].r;
  p_i=partinfo[i].p;
  p_j=partinfo[j].p;
  ptot2=0.0;
  q0dotb=0.0;
  ptotdotb=0.0;
  m2j=0.0;m2i=0.0;
  delt_j=(tevent+0.5*delt-r_j[0]);
  delt_i=(tevent-0.5*delt-r_i[0]);
  for(alpha=0;alpha<4;alpha++){
    m2i=m2i+g[alpha]*p_i[alpha]*p_i[alpha];
    m2j=m2j+g[alpha]*p_j[alpha]*p_j[alpha];
    ptot[alpha]=p_i[alpha]+p_j[alpha];
    ptot2=ptot2+g[alpha]*ptot[alpha]*ptot[alpha];
    q0[alpha]=p_j[alpha]-p_i[alpha];
    b[alpha]=r_j[alpha]-r_i[alpha];
    r_j[alpha]=r_j[alpha]+delt_j*p_j[alpha]/p_j[0];
    r_i[alpha]=r_i[alpha]+delt_i*p_i[alpha]/p_i[0];
    q0dotb=q0dotb+g[alpha]*q0[alpha]*b[alpha];
    ptotdotb=ptotdotb+g[alpha]*ptot[alpha]*b[alpha];
  }
#ifdef EXTRA_TESTS
  if(m2j<-0.001||m2i<-0.001) printf("Trouble in elast.. %g %g\n",m2j,m2i);
#endif
  ptotdotq0=m2j-m2i;
  q02=-ptot2+2.0*(m2i+m2j);
  q2=0.0;b2=0.0;
#ifdef EXTRA_TESTS
  qdotbtest=0.0;
  pdotbtest=0.0;
#endif
  for(alpha=0;alpha<4;alpha++){
    q[alpha]=q0[alpha]-ptotdotq0*ptot[alpha]/ptot2;
    q2=q2+g[alpha]*q[alpha]*q[alpha];
    b[alpha]=b[alpha]-(ptot[alpha]*ptotdotb/ptot2)
      -(q[alpha]*(q0dotb-ptotdotb*ptotdotq0/ptot2))
      /(q02-(ptotdotq0*ptotdotq0/ptot2));
    b2=b2+g[alpha]*b[alpha]*b[alpha];
#ifdef EXTRA_TESTS
    qdotbtest=qdotbtest+q[alpha]*b[alpha]*g[alpha];
    pdotbtest=pdotbtest+ptot[alpha]*b[alpha]*g[alpha];
#endif
  }

#ifdef EXTRA_TESTS
  if(fabs(qdotbtest)>0.001||fabs(pdotbtest)>0.001){
    printf("qdotbtest or pdotbtest failure: %g %g\n",qdotbtest,pdotbtest);
    exit(1);
  }
#endif

  normq=sqrt(fabs(-q2));
  normb=sqrt(fabs(-b2));
  theta=pi*ran2();
  cthet=cos(theta);
  sthet=sqrt(1.0-cthet*cthet);
  for(alpha=0;alpha<4;alpha++){
    q[alpha]=(cthet*q[alpha])+(normq*sthet*b[alpha]/normb);
    q0[alpha]=q[alpha]+ptot[alpha]*ptotdotq0/ptot2;
    p_j[alpha]=0.5*(ptot[alpha]+q0[alpha]);
    p_i[alpha]=0.5*(ptot[alpha]-q0[alpha]);
  }

#ifdef EXTRA_TESTS
  m2itest=m2i;
  m2jtest=m2j;
#endif
  m2i=0.0;m2j=0.0;
  for(alpha=0;alpha<4;alpha++){
    m2i=m2i+g[alpha]*p_i[alpha]*p_i[alpha];
    m2j=m2j+g[alpha]*p_j[alpha]*p_j[alpha];
  }
#ifdef EXTRA_TESTS
  if(fabs(m2itest-m2i)>0.001||fabs(m2jtest-m2j)>0.001){
    printf("i,j=%d,%d   ires_i,ires_j=%d,%d\n",i,j,partinfo[i].ires,
	   partinfo[j].ires);
    printf("m2i=%f,m2j=%f\n",m2i,m2j);
    printf("Failed scattering, m_i(before/after) %f %f m_j(before/after) %f %f\n",m2itest,m2i,m2jtest,m2j);
    exit(1);
  }
#endif
  partinfo[i].mass=sqrt(fabs(m2i));
  partinfo[j].mass=sqrt(fabs(m2j));
  //printf("Successful elastic scattering  %d %d  b=%f\n",
  //partinfo[i].ires,partinfo[j].ires,sqrt(fabs(-b2)));
}

/*---------------------------------------------------------------------------*/

void resonancemaker(int ievent,int *kf){
  int alpha,scatcode,i,j;
  double delt,tevent,delt_i,delt_j;
  double *p_i,*p_j,*p_k,*r_i,*r_j,*r_k,*b;

  i=eventinfo[ievent].i;
  j=eventinfo[ievent].j;
  scatcode=eventinfo[ievent].scatcode;
  tevent=eventinfo[ievent].time;
  delt=eventinfo[ievent].delt;

#ifdef EXTRA_TESTS
  if(i==j) {
    printf("i=j in resonance maker %d %d\n",i,j);
    exit(1);
  }
  if(scatcode>nres||scatcode<=0){
    printf("scatcode inappropriate in resmaker, =%g\n",scatcode);
    exit(1);
  }
#endif
  kf[0]=1;
  kf[1]=firstemptypart;
  firstemptypart=partinfo[firstemptypart].next;
  partinfo[kf[1]].next=firstpart;
  partinfo[kf[1]].ires=scatcode;
  firstpart=kf[1];
  
  p_i=partinfo[i].p;
  p_j=partinfo[j].p;
  p_k=partinfo[kf[1]].p;
  r_i=partinfo[i].r;
  r_j=partinfo[j].r;
  r_k=partinfo[kf[1]].r;
  delt_j=tevent+0.5*delt-r_j[0];
  delt_i=tevent-0.5*delt-r_i[0];
  for(alpha=0;alpha<4;alpha++){
    p_k[alpha]=p_i[alpha]+p_j[alpha];
    r_j[alpha]=r_j[alpha]+delt_j*p_j[alpha]/p_j[0];
    r_i[alpha]=r_i[alpha]+delt_i*p_i[alpha]/p_i[0];
    partinfo[kf[1]].b[alpha]=r_j[alpha]-r_i[alpha];
    r_k[alpha]=0.5*(r_i[alpha]+r_j[alpha]);
  }
  partinfo[kf[1]].mass=sqrt(fabs(p_k[0]*p_k[0]-
			    p_k[1]*p_k[1]-p_k[2]*p_k[2]-p_k[3]*p_k[3]));
  npart=npart+1;
#ifdef EXTRA_TESTS
  if(resinfo[scatcode].minmass>partinfo[kf[1]].mass){
    printf("Underneath Mass Min. in Res.Maker\n");
    printf("scatcode=%d  threshold:%g mass:%g\n",scatcode,resinfo[scatcode].minmass,partinfo[kf[1]].mass);
    printf("i,ires_i= %d %d    j ires_j= %d %d \n",i,partinfo[i].ires,j,partinfo[j].ires);
    exit(1);
  }
#endif
  //printf("+++++++++++ Resonance produced, scatcode=%d parts: %d\n",
	 //scatcode,kf[1]);
}

/*---------------------------------------------------------------------------*/

void resonancedecay(int ievent, int *kf){
  int ires[4],ibody,nbodies,alpha,i;
  double tevent,b[4];
  double u[4],mass[4];
  double p[4][4],kprime[4];
  double e1,e2,randy,branchtot,masstot;
  double q,weight,wmax,sthet,cthet,phi;
  double p3mag,kprimemax,p3max,kprimemax2,kprimemag2,e1prime,e2prime;
  double e12,u12[4],test,oldtime;
#ifdef EXTRA_TESTS
  double checker;
#endif
  int bptr,nextbptr;

  i=eventinfo[ievent].i;
  for(alpha=0;alpha<4;alpha++) b[alpha]=partinfo[i].b[alpha];
  ires[0]=-eventinfo[ievent].scatcode;
  tevent=eventinfo[ievent].time;

  /* Find the channel */
  
  randy=ran2();
#ifdef EXTRA_TESTS
  checker=partinfo[i].p[0]*partinfo[i].p[0];
  for(alpha=1;alpha<4;alpha++) checker=checker-
				 partinfo[i].p[alpha]*partinfo[i].p[alpha];
  if(checker<0){
    printf("in resdecay, decaying particle is time like\n");
    exit(1);
  }
  if(fabs(sqrt(fabs(checker))-partinfo[i].mass)>1.0){
    printf("masses don't compute in resdecay %g  %g\n",
	   sqrt(fabs(checker)),partinfo[i].mass);
    exit(1);
  }
#endif
  oldtime=partinfo[i].r[0];
  for(alpha=0;alpha<4;alpha++) {
    p[0][alpha]=partinfo[i].p[alpha];
    u[alpha]=p[0][alpha]/partinfo[i].mass;
    partinfo[i].r[alpha]=
      partinfo[i].r[alpha]+(u[alpha]/u[0])*(tevent-oldtime);
  }
  branchtot=0.0;
  bptr=resinfo[ires[0]].firstbranch;
  do{
    branchtot=branchtot+branchinfo[bptr].branching;
    nextbptr=branchinfo[bptr].nextbranch;
    if(nextbptr==-1||randy<branchtot) goto BPTRFIND;
    bptr=nextbptr;
  }while(bptr!=-1);
  printf("Should not be here !!!!\n");
  exit(1);
BPTRFIND:
  mass[0]=partinfo[i].mass;
  nbodies=branchinfo[bptr].nbodies;
  masstot=0.0;
  for(ibody=1;ibody<=nbodies;ibody++) {
    ires[ibody]=branchinfo[bptr].ires[ibody];
    do{
      mass[ibody]=resinfo[ires[ibody]].mass+
	0.5*resinfo[ires[ibody]].width*tan(pi*(ran2()-0.5));
    } while(mass[ibody]<resinfo[ires[ibody]].minmass);
    masstot=masstot+mass[ibody];
  }

  /* CANCEL DECAY ( -> ONE-BODY DECAY ) IF FORBIDDEN */
  if(masstot>mass[0]) {
    nbodies=1;
    ires[1]=partinfo[i].ires;;
    p[1][0]=mass[0];
    //printf("mass is %g\n",mass[0]);
    for(alpha=1;alpha<4;alpha++) p[1][alpha]=0.0;
  }

  /* TWO-BODY DECAYS */
  else if(nbodies==2){
    q=sqrt(triangle(mass[0],mass[1],mass[2]));
    cthet=1.0-2.0*ran2();
    sthet=sqrt(1.0-cthet*cthet);
    phi=2.0*pi*ran2();
    p[1][3]=q*cthet;
    p[1][1]=q*sthet*cos(phi);
    p[1][2]=q*sthet*sin(phi);
    p[2][3]=-p[1][3];
    p[2][2]=-p[1][2];
    p[2][1]=-p[1][1];
    p[1][0]=sqrt(mass[1]*mass[1]+p[1][1]*p[1][1]
		 +p[1][2]*p[1][2]+p[1][3]*p[1][3]);
    p[2][0]=sqrt(mass[2]*mass[2]+p[2][1]*p[2][1]
		 +p[2][2]*p[2][2]+p[2][3]*p[2][3]);
  }

  /* THREE-BODY DECAYS */
  else if(nbodies==3){
    //printf("Now attempting a 3-body decay\n");
    //printf("masses: %g -> %g %g %g\n",mass[0],mass[1],mass[2],mass[3]);
    kprimemax2=triangle(mass[0]-mass[3],mass[1],mass[2]);
    kprimemax=sqrt(kprimemax2);
    p3max=sqrt(triangle(mass[0],mass[1]+mass[2],mass[3]));
    e1=sqrt(pow(mass[1],2)+p3max*p3max);
    e2=sqrt(pow(mass[2],2)+p3max*p3max);
    wmax=p3max*(e1*e2/(mass[1]*mass[2]))*(mass[1]+mass[2])/(e1+e2);
    do{
    TRY_AGAIN:
      do{
	kprime[1]=kprimemax*(2.0*ran2()-1.0);
	kprime[2]=kprimemax*(2.0*ran2()-1.0);
	kprime[3]=kprimemax*(2.0*ran2()-1.0);
	kprimemag2=kprime[1]*kprime[1]+
	  kprime[2]*kprime[2]+kprime[3]*kprime[3];
      } while(kprimemag2>kprimemax2);
      e1prime=sqrt(kprimemag2+mass[1]*mass[1]);
      e2prime=sqrt(kprimemag2+mass[2]*mass[2]);
      if(e1prime+e2prime+mass[3]>mass[0]) goto TRY_AGAIN;
      p3mag=sqrt(triangle(mass[0],e1prime+e2prime,mass[3]));
      cthet=1.0-2.0*ran2();
      sthet=sqrt(1.0-cthet*cthet);
      phi=2.0*pi*ran2();
      p[3][3]=p3mag*cthet;
      p[3][1]=p3mag*sthet*cos(phi);
      p[3][2]=p3mag*sthet*sin(phi);
      p[3][0]=sqrt(p3mag*p3mag+mass[3]*mass[3]);
      e12=sqrt(pow(e1prime+e2prime,2)+p3mag*p3mag);
      for(alpha=1;alpha<=3;alpha++) u12[alpha]=-p[3][alpha]/e12;
      u12[0]=sqrt(1.0+u12[1]*u12[1]+u12[2]*u12[2]+u12[3]*u12[3]);
      kprime[0]=e1prime;
      lorentz1(u12,kprime,&p[1][0]);
      kprime[0]=e2prime;
      for(alpha=1;alpha<=3;alpha++) kprime[alpha]=-kprime[alpha];
      lorentz1(u12,kprime,&p[2][0]);
      weight=p3mag*(p[1][0]*p[2][0]/(e1prime*e2prime))
	*((e1prime+e2prime)/(p[1][0]*p[2][0]));
    } while(ran2()<weight/wmax);
  }
  
  /* ADD THE NEW PARTICLES TO THE LIST */
  npart=npart+nbodies;
  kf[0]=nbodies;
  for(ibody=1;ibody<=nbodies;ibody++){
    kf[ibody]=firstemptypart;
    firstemptypart=partinfo[firstemptypart].next;
    partinfo[kf[ibody]].next=firstpart;
    partinfo[kf[ibody]].ires=ires[ibody];
    firstpart=kf[ibody];
    partinfo[kf[ibody]].ires=ires[ibody];
    for(alpha=0;alpha<4;alpha++) {
      if(nbodies!=2) partinfo[kf[ibody]].r[alpha]=partinfo[i].r[alpha];
      if(nbodies==2&&ibody==1)
	partinfo[kf[ibody]].r[alpha]=partinfo[i].r[alpha]-0.5*b[alpha];
      if(nbodies==2&&ibody==2)
	partinfo[kf[ibody]].r[alpha]=partinfo[i].r[alpha]+0.5*b[alpha];
      partinfo[kf[ibody]].b[alpha]=0.0;
    }
  }
  if(nbodies==1) {
    lorentz1(u,&p[1][0],&(partinfo[kf[1]].p[0]));
    //printf("%%%%%%%%%%%%%%%%%  ONE-BODY DECAY %%%%%%%%%%%%%\n");
  }
  if(nbodies==2) {
    lorentz2( u,&(p[1][0]),partinfo[kf[1]].p,&(p[2][0]),
	     partinfo[kf[2]].p );
  }
  if(nbodies==3) {
    //printf("Doing 3-body decay\n");
#ifdef EXTRA_TESTS
    checker=u[0]*u[0]-u[1]*u[1]-u[2]*u[2]-u[3]*u[3]-1.0;
    if(fabs(checker)>0.0001){
      printf("usquared not unity\n");
      printf("%g %g %g %g\n",u[0],u[1],u[2],u[3]);
      exit(1);
    }
#endif
    lorentz3( u,&p[1][0],&partinfo[kf[1]].p[0],
	     &p[2][0],&partinfo[kf[2]].p[0],
	     &p[3][0],&partinfo[kf[3]].p[0] );
  }
  /* Test to see if particles are time-like */
  for(ibody=1;ibody<=nbodies;ibody++){
    test=p[ibody][0]*p[ibody][0];
    for(alpha=1;alpha<=3;alpha++) test=test-p[ibody][alpha]*p[ibody][alpha];
#ifdef EXTRA_TESTS
    if(test<0.0&&ires[ibody]>1) {
      printf("MAJOR FAULT, IMAGINARY MASS ibody=%d  nbodies=%d\n",ibody,nbodies);
      exit(1);
    }
#endif
    partinfo[kf[ibody]].mass=sqrt(fabs(test));
  }
  //printf("++++++++++++++++++Resonance decayed,nbodies=%d, parts: %d %d\n",
  //nbodies,kf[0],kf[1]);
}

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/