/********************************************************************\

  Name:         histo.c
  Created by:   Andre Lebedev

  Contents:     Produce files needed to do secondary TDC correction


\********************************************************************/
                                                        
/*-- Include files -------------------------------------------------*/

/* standard includes */
#include <stdio.h>
#include <math.h>

/* midas includes */
#include "midas.h"
#include "experim.h"
#include "analyzer.h"

/* cernlib includes */
#ifdef OS_WINNT
#define VISUAL_CPLUSPLUS
#endif

#ifndef sqr
#define sqr(x) ((x)*(x))
#endif

#ifdef OS_LINUX
#define f2cFortran
#endif

#include <cfortran.h>
#include <hbook.h>

/*-- Parameters ----------------------------------------------------*/

#ifndef ADC_TDC_LINEAR
#define ADC_TDC_LINEAR
#define N_MAX 130
#define MY_N_ADC 324
#endif

struct {
  float tdc[MY_N_ADC][N_MAX];
  float adc[MY_N_ADC][N_MAX];
  float tdc_err[MY_N_ADC][N_MAX];
  float adc_err[MY_N_ADC][N_MAX];
  INT n_entries[MY_N_ADC][N_MAX];
  INT valid[MY_N_ADC];
} adc_tdc_param;

extern ADC_CALIBRATION_PARAM adccalib_param; 

/*-- Static parameters ---------------------------------------------*/

/*-- Module declaration --------------------------------------------*/

INT histo_filling(EVENT_HEADER*,void*);
INT histo_booking(void);
INT histo_eor(INT run_number);

ANA_MODULE histo_module = {
  "General histo",               /* module name           */  
  "Stefan Ritt",                 /* author                */
  histo_filling,                 /* event routine         */
  NULL,                          /* BOR routine           */
  histo_eor,                     /* EOR routine           */
  histo_booking,                 /* init routine          */
  NULL,                          /* exit routine          */
  NULL,                          /* parameter structure   */
  0,                             /* structure size        */
  NULL,                          /* initial parameters    */
};

/*-- init routine --------------------------------------------------*/

INT histo_booking(void)
{
INT i, j;
char str[50];

  for(j=0;j<MY_N_ADC;j++)
    for(i=0;i<N_MAX;i++)
      {
      adc_tdc_param.tdc[j][i]=0.f;
      adc_tdc_param.adc[j][i]=0.f;
      adc_tdc_param.tdc_err[j][i]=0.f;
      adc_tdc_param.adc_err[j][i]=0.f;
      adc_tdc_param.n_entries[j][i]=0;
      adc_tdc_param.valid[j]=0;
      }

  /* Book 2-D histograms for each trigger */
  for (i=0; i<8; i++)
    {
    sprintf(str, "CsI ADC vs. TDC, TRIG=%d", 1<<i);
    HBOOK2(1990+i, str, 40, -20.f, 20.f, 100, 0.f, 3000.f, 0.f);
    sprintf(str, "PV ADC vs. TDC, TRIG=%d", 1<<i);
    HBOOK2(1980+i, str, 40, -20.f, 20.f, 100, 0.f, 3000.f, 0.f);
    }

  return SUCCESS;
}
                                 
/*-- event routine -------------------------------------------------*/
                                                                                       
INT histo_filling(EVENT_HEADER *pheader, void *pevent)
{
float        *cadc, m_adc;
float        *ttim, m_tdc;
float        ADC_MIN, BIN_0, BIN_N;
WORD         *trig;
INT i, j, n;

  /* open banks of interest */
  if (bk_locate(pevent, "CADC", &cadc) == 0 ||
      bk_locate(pevent, "TTIM", &ttim) == 0 ||
      bk_locate(pevent, "TRIG", &trig) == 0)
    /* return of one of these banks is not present */
    return 1;

  /* Fill in 2-D ADC vs. TDC histos for each trigger, for all  */
  for(i=0; i<8; i++)
    if((trig[0]) & (1<<i))
      for(j=0; j<240; j++)
	if(cadc[j] > 10.)
	  {
	  if(j<40)
	    HF2(1980+i, ttim[240+j], 
		(cadc[240+j]/adccalib_param.mev_per_channel_pv), 1.f);
	  HF2(1990+i, ttim[j], 
	      (cadc[j]/adccalib_param.mev_per_channel_csi), 1.f);
	  }

  /* Calculate the average of TDC and ADC value in each bin */
  for (i=0; i<MY_N_ADC; i++){

    /**********************************************************\
     One should adjust the numbers below if the numbers of
     ADC and TDC channels change.
     \*********************************************************/

    m_adc = -1.0;  /* Initialize default source of data as "dummy noise" */
    m_tdc = -1000.;

    if (i<240) /* CsI crystals */
      {
      ADC_MIN = 105.0;
      BIN_0 = 6.0;
      BIN_N = 35.0;
      m_adc=cadc[i] / adccalib_param.mev_per_channel_csi;
      m_tdc=ttim[i];
      if ((trig[0] & 8) == 0) /* Select pienu trigger for CsI */
         m_adc = -1.;
      }
    else if (i < CHN_PVETO+N_PVETO)
      {
      ADC_MIN = 35.0;
      BIN_0 = 2.0;
      BIN_N = 35.0;
      m_adc=cadc[i] / adccalib_param.mev_per_channel_pv;
      m_tdc=ttim[i];
      if ((trig[0] & 16) == 0) /* Select pienu: hi:PS' trigger for PV's */
         m_adc = -1.;
      }
    else if(i == 300) /* Beam Counters */
      {
      ADC_MIN = 5.0;
      BIN_0 = 1.0;
      BIN_N = 18.0;
      m_adc=cadc[392];
      m_tdc=ttim[400]-ttim[328];
      if(m_tdc < -8 || m_tdc>8)
	i++;
      if ((trig[0] & 128) == 0)
	m_adc = -1;
      }
    else if(i == 302 || i == 303) /* Degrader */
      {
      ADC_MIN = 5.0;
      BIN_0 = 4.0;
      BIN_N = 18.0;
      m_adc=cadc[393];
      m_tdc=ttim[i+99]-ttim[328];  /* DEG:L=TDC401, DEG:H=TDC402 */
      if ((trig[0] & 128) == 0)
	m_adc = -1;
      }
    else if(i >= 306 && i <= 314) /* Individual target counters */
      {
      cadc[i] /= adccalib_param.mev_per_channel_tgt; 
      ADC_MIN = 5.0;
      BIN_0 = 2.0;
      BIN_N = 20.0;
      m_adc=cadc[i];
      m_tdc=ttim[i];
      if ((trig[0] & 128) == 0)
	m_adc = -1;
      }
    else if(i >= 317 && i <= 320) /* Sum Target */
      {
      ADC_MIN = 5.0;
      BIN_0 = 2.0;
      BIN_N = 20.0;
      m_adc=cadc[78+i];
      m_tdc=ttim[i+86]-ttim[328]; /* sumT1:4=TDC403, sum T5:8=TDC404, 
			   T:L=TDC405, T:H=TDC406 */
      if ((trig[0] & 128) == 0)
	m_adc = -1;
      }
    else if(i == 321 || i == 322) /* Sum Active collimators */
      {
      ADC_MIN = 5.0;
      BIN_0 = 4.0;
      BIN_N = 18.0;
      m_adc = cadc[i+76];
      m_tdc = ttim[i+86]; /* AC0:3=TDC407, AC4:7=TDC408 */
      }

    /* Impose limits on acceptability of an event */
    if (m_adc >= ADC_MIN && 
	m_tdc >= -11. &&
	m_tdc <= 18.)
      {
      j = (INT) ((float) 0.5 * (-fabs((BIN_N - BIN_0 * (2 * N_MAX + 1)) /
                (BIN_N - BIN_0)) + 
	        sqrt(sqr((BIN_N - BIN_0 * (2 * N_MAX + 1)) / (BIN_N - BIN_0)) +
                8 * N_MAX * (m_adc - ADC_MIN) / (BIN_N - BIN_0))));

      if (j >= N_MAX)  /* ADC value too high */
	return ANA_CONTINUE ;

      adc_tdc_param.valid[i] = 1;
      n = adc_tdc_param.n_entries[i][j];
      adc_tdc_param.adc[i][j] *= (float) n/(n+1);
      adc_tdc_param.adc[i][j] += (float) m_adc/(n+1);
      adc_tdc_param.tdc[i][j] *= (float) n/(n+1);
      adc_tdc_param.tdc[i][j] += (float) m_tdc/(n+1);
      adc_tdc_param.adc_err[i][j] *= (float) n/(n+1);
      adc_tdc_param.adc_err[i][j] += (float) sqr(m_adc)/(n+1);
      adc_tdc_param.tdc_err[i][j] *= (float) n/(n+1);
      adc_tdc_param.tdc_err[i][j] += (float) sqr(m_tdc)/(n+1);
      adc_tdc_param.n_entries[i][j]++;
      }
  }
  return ANA_CONTINUE;
}

/* ----------- EOR routine --------------------- */

INT histo_eor(INT run_number)
{

INT i, j, n;
float adc, tdc, adc_err, tdc_err;
char str[80];
FILE *fi;

   for(i=0;i<MY_N_ADC;i++)
     {
     if(adc_tdc_param.valid[i] == 1)
     {
     sprintf(str, "tdc_data/tdc%03d.dat",i);
     fi=fopen(str,"w");
     if(fi==NULL)
        cm_msg(MERROR,"adc_tdc_eor", "Cannot open %s", str);
     for(j=0;j<N_MAX;j++)
       {
       n = adc_tdc_param.n_entries[i][j];
      if(n != 0)
	 {
	 adc = adc_tdc_param.adc[i][j];
	 tdc = adc_tdc_param.tdc[i][j];
	 if(n > 1)
	   {
	   adc_err = adc_tdc_param.adc_err[i][j] - sqr(adc);
	   if(adc_err < 0.005)
	     adc_err = 1.0;
           adc_err = sqrt((float)adc_err/(n-1));
	   tdc_err = adc_tdc_param.tdc_err[i][j] - sqr(tdc);
	   if(tdc_err < 0.003)
	     tdc_err = 1.0;
           tdc_err=sqrt((float)tdc_err/(n-1));
	   }
	 else
	   {
	   adc_err=5.f;
	   tdc_err=3.f;
	   }
	 fprintf (fi, "%4.2f\t%4.2f\t%3.3f\t%3.3f\n", tdc, adc, 
		 tdc_err, adc_err);
	 }
       }
     fclose(fi);
     }
     else
       printf("File \"tdc%03d.dat\" contains no data; it wasn't created.\n",i);
     }
   printf("Data files of TDC vs ADC values were written to 
          \"tdc_data/tdcXXX.dat\".\n");
return SUCCESS;
}