#include <math.h>
// *************************************************************
// * create hitsel object and time-of-flight array             *
// *************************************************************
inline likelihood::likelihood(float r,float z):
  fitquality(r,z,bon_dwallfit())
{
  nlike=0;
}

inline void likelihood::set_hits(hitsel *h)
{
  event_hits=h;
  nlike=0;
  if (event_hits!=NULL) event_hits->qsort();
}

// *************************************************************
// * calculate dodecahedran of search radius and all thirteen  *
// * qualitites                                                *
// *************************************************************
inline void likelihood::check_around(float *vertex,float *result,
		                     float radius,float *q,int &max)
{
  int i;

  orientation.set(result[2],result[3],result[4]);
  dod.surround(vertex+3,orientation,radius);
  for(max=i=0; i<13; i++)
    {
      if ((q[i]=quality(vertex+3*i))>q[max])
	{
	  max=i;
	  *result=verfit[3];
	  result[1]=like0;
	  result[2]=*dirfit;
	  result[3]=dirfit[1];
	  result[4]=dirfit[2];
	  result[5]=dirfit[3];
	  result[6]=dirfit[4];
	}
#ifdef DEBUG_TWO
      printf("%2d(%2d) (%8.2f,%8.2f,%8.2f,%6.1f): %8.6f(%8.6f) (%6.1f,%6.1f,%6.1f) %4.1f %8.5f\n",
	     i,max,vertex[3*i],vertex[3*i+1],vertex[3*i+2],verfit[3],q[i],
	     like0,dirfit[0]*57.29578,dirfit[1]*57.29578,
	     dirfit[2]*57.29578,acos(dirfit[3])*57.29578,dirfit[4]);
#endif
    }
}

// *************************************************************
// * return number of checked points for each iteration        *
// *************************************************************
inline char likelihood::ncheck(void)
{
  return(13);
}

// *************************************************************
// * calculate quality at interpolated point                   *
// *************************************************************
inline void likelihood::interpolate(float *vertex,
                                    float radius,float *q,float *inter)
{
  dod.interpolate(q,inter);
  vertex[39]=vertex[0]+radius*orientation.getx(inter);
  vertex[40]=vertex[1]+radius*orientation.gety(inter);
  vertex[41]=vertex[2]+radius*orientation.getz(inter);
  q[13]=quality(vertex+39);
#ifdef DEBUG_TWO
  printf("interpolate to (%8.2f,%8.2f,%8.2f): %8.6f\n",
         vertex[39],vertex[40],vertex[41],q[13]);
#endif
}

// *************************************************************
// * maximize criterion likelihood using bonsaifit object      *
// *************************************************************
inline void likelihood::maximize(bonsaifit *fit,searchgrid *grid)
{
  if (grid->nset()<1) return;
  //initializing grid search
  cang0=bangle();
  plusdang=bplusdevangle();
  minusdang=bminusdevangle();
#ifdef DEBUG
  printf("Initial search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	 cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
  if (grid->nset()<2)
    {
      if (grid->size(0)+grid->size(1)<1) return;
      fit->search(obon_minimum_radius(),grid,1);
    }
  else
    {
      if(bongrid()<clusgrid())
        {
          if (grid->size(0)+grid->size(3)<1) return;
          fit->search(obon_minimum_radius(),grid,3);
	}
      else
	{
      	  if (grid->size(0)+grid->size(4)<1) return;
          fit->search(obon_minimum_radius(),grid,4);
	}
    }

  //original search
  cang0=oangle();
  plusdang=oplusdevangle();
  minusdang=ominusdevangle();
#ifdef DEBUG
  printf("Coarse search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	 cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
  fit->search(obon_min_gdn_difference(),obon_gdn_fraction(),
	      obon_minimum_radius(),obon_stop_radius());
  if (walldist(fit->xfit(),fit->yfit(),fit->zfit())<bon_dwall())
    {
#ifdef DEBUG
      printf("dwall=%8.2f<%8.2f: Final search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	      walldist(fit->xfit(),fit->yfit(),fit->zfit()),
              bon_dwall(),cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
      fit->search(bon_last_min_gdn_difference(),
	              bon_last_gdn_fraction());
#ifdef DEBUG
      fit->print_branch_list();
#endif
      return;
    }
  // conditional search
  cang0=cangle();
  plusdang=cplusdevangle();
  minusdang=cminusdevangle();
#ifdef DEBUG
    printf("Fine search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	   cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
  fit->search(cbon_min_gdn_difference(),cbon_gdn_fraction(),
	        cbon_minimum_radius(),cbon_stop_radius());
  minusdang=cminusdevangle();
#ifdef DEBUG
    printf("Final search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	   cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
  fit->search(bon_last_min_gdn_difference(),
	      bon_last_gdn_fraction());
#ifdef DEBUG
  fit->print_branch_list();
#endif
}

// *************************************************************
// * maximize criterion likelihood using one starting point    *
// *************************************************************
inline void likelihood::maximize(bonsaifit *fit,float *point)
{
  cang0=cangle();
  plusdang=cplusdevangle();
  minusdang=cminusdevangle();
#ifdef DEBUG
    printf("Fine search: cos_theta=%6.2f+%6.2f-%6.2f\n",
	   cang0,sqrt(0.5/plusdang),sqrt(0.5/minusdang));
#endif
  fit->search(cbon_minimum_radius(),cbon_stop_radius(),point);
}

// *************************************************************
// * get last fitted light emission time                       *
// *************************************************************
inline float likelihood::get_zero(void)
{
  return(verfit[3]);
}

// *************************************************************
// * get last fitted likelihood                                  *
// *************************************************************
inline float likelihood::get_ll(void)
{
  return(like);
}

// *************************************************************
// * get last fitted likelihood                                  *
// *************************************************************
inline float likelihood::get_ll0(void)
{
  return(like0);
}

inline void likelihood::get_dir(float *dir)
{
  *dir=*dirfit;
  dir[1]=dirfit[1];
  dir[2]=dirfit[2];
  dir[3]=dirfit[3];
  dir[4]=dirfit[4];
}

// *************************************************************
// * get number of times, quality was called                   *
// *************************************************************
inline int likelihood::nfit(void)
{
  return(nlike);
}

// *************************************************************
// * return size of output memory segment                      *
// *************************************************************
inline int likelihood::nresult(void)
{
  return(7);
};

// *************************************************************
// * copy emission time, raw likelihood and direction fit to   *
// * memory segment r                                          *
// *************************************************************
inline void likelihood::get_result(float *r)
{
  *r=verfit[3];
  r[1]=like0;
  r[2]=*dirfit;
  r[3]=dirfit[1];
  r[4]=dirfit[2];
  r[5]=dirfit[3];
  r[6]=dirfit[4];
  return;
};

// *************************************************************
// * copy emission time, raw likelihood and direction fit from *
// *  memory segment r                                         *
// *************************************************************
inline void likelihood::set_result(float *r)
{
  verfit[3]=*r;
  like0=r[1];
  *dirfit=r[2];
  dirfit[1]=r[3];
  dirfit[2]=r[4];
  dirfit[3]=r[5];
  dirfit[4]=r[6];
  return;
};