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

  Name:         histo.c
  Created by:   Emil Frlez
  updated by the whole PiBeta group
  Contents:     General purpos histogram booking/filling module

  $Revision: 1.16 $

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

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

/* standard includes */
#include <stdio.h>
#include <math.h>
#include "util.h"
#include <string.h>      /* WL */

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

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

#ifdef OS_LINUX
#define f2cFortran
#endif

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

#ifndef PI
#define PI 3.14159265359
#endif

/*-- Parameters ----------------------------------------------------*/
extern double csi_theta[240];
extern double csi_phi[240];
extern float adc_dsc[N_DSC*132];
extern float peak_dsc_tg;
extern int npienu = 99;
extern int npibeta1=0;
extern int npibeta2=0;
extern float tgt_dsc[132];
extern float del_tgt_dsc[132];
extern float b0_dsc[132];

/* for refined tdc offsets WL*/
typedef struct{
  INT channel;
  float time;
}TIM_HISTO;

#define N_DATA 350
extern INT n_data_histo;
extern TIM_HISTO tim_histo[N_DATA];
extern float time_deg_normal;
/*----*/

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

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

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

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


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

#define ADC_N_BINS   400
#define ADC_X_LOW      0
#define ADC_X_HIGH  2000

#define TDC_N_BINS   800
#define TDC_X_LOW   -200
#define TDC_X_HIGH   200

#define CSI_SUM_N_BINS   500
#define CSI_SUM_X_LOW      0
#define CSI_SUM_X_HIGH   160 /* if changed, also change gainmatch.c and refhistos.rz ! */

#define DISC_MIN 10.f
#define DISC_MIN1 3.f

#define Y_CV 1327.f

#define CHN_B0_TDC     400
#define CHN_DEGL_TDC   401
#define CHN_DEGH_TDC   402
#define CHN_T1_4_TDC   403
#define CHN_T5_8_TDC   404
#define CHN_TGTL_TDC   405
#define CHN_TGTH_TDC   406
#define CHN_AC0_3_TDC  407
#define CHN_AC4_7_TDC  408
#define CHN_PV_TDC     409
#define CHN_CMV_TDC    410

#define B0_GAIN  0.0130
#define DEG_GAIN 0.0472
#define TGT_GAIN 0.0748

#define csi_toff 23.3f
#define pv_toff  21.5f
#define b0_peak  528.f
#define deg_peak 2330.f
#define b0_w     250.f
#define deg_w    669.f

/* These are for tdc offset fit WL*/
#define TDC_N_BINS_HISTO        201
#define TDC_X_LOW_HISTO        -50
#define TDC_X_HIGH_HISTO        50
#define TDC_OFFSET_PROMPT     25000
#define TDC_OFFSET_PIENUH_G   25001
#define TDC_OFFSET_PIENUH_E   25002 
#define TDC_LIN_PIENUH_G_BASE      27000
#define TDC_LIN_PIENUH_E_BASE      28000
#define N_RUNS      10          /* number of runs to do tdcoffset WL*/

TDC_CALIBRATION_PARAM tdccalib_param;

int run_counter = 0;       /* count runs WL*/
INT   i,j,csi_ok;
char  str[80];



INT histo_booking(void)

{
  HBOOK2(TDC_OFFSET_PROMPT, "TDC_prompt ", 
         TDC_N_BINS_HISTO, (float)TDC_X_LOW_HISTO, (float)TDC_X_HIGH_HISTO, 
         N_TDC, (float)0, (float)N_TDC, 0.f);
  HBOOK2(TDC_OFFSET_PIENUH_G, "TDC_pienu ", 
         TDC_N_BINS_HISTO, (float)TDC_X_LOW_HISTO, (float)TDC_X_HIGH_HISTO, 
         N_TDC, (float)0, (float)N_TDC, 0.f);
  HBOOK2(TDC_OFFSET_PIENUH_E, "TDC_pienu_positron ", 
         TDC_N_BINS_HISTO, (float)TDC_X_LOW_HISTO, (float)TDC_X_HIGH_HISTO, 
         N_TDC, (float)0, (float)N_TDC, 0.f);

  return SUCCESS;
}

/*-- bor routine ---------------------------------------------------*/

INT histo_bor(INT run_number)
{
  if(run_counter == N_RUNS)
    {   
      HRESET(TDC_OFFSET_PROMPT, " ");
      HRESET(TDC_OFFSET_PIENUH_G, " ");
      HRESET(TDC_OFFSET_PIENUH_E, " ");
    }
 return SUCCESS;
}


/*-- event routine -------------------------------------------------*/

INT histo_filling(EVENT_HEADER *pheader, void *pevent)
{
float        *cadc;
float        *ttim;
float        *tpre;
float        *tpst;
float        *twid;
float        badc[12];
WORD         *tbit,*trig,*nhit,*Badc;
ASUM_BANK    *asum;
PVET_BANK    *pvet;
MWPC_BANK    *mwpc;
CLMP_BANK    *clmp;
BEAM_BANK    *beam;
TEMP_BANK    *temp;
TRAK_BANK    *trak;
double       track_phi, track_theta;
int tgt_n,deg_n;
int n_csi, nhit_csi;
float col1sum,col2sum;
float tgt_mt,deg_mt;
float tgt_sum,tgt_adcsum,deg_sum,deg_adcsum;

float        ratio12;
float        de_tgt;  
float        toffset;
float        best_pienu, best_pi0;

float mpi_inv[5][5];
float x,t,min,psq;
float tht,theta1,theta2,phi1,phi2;
float a12,a13,a23;
int i,j,i1,i2,ip;       

float x_cv_intersection, z_cv_intersection;

float adc_t0,adc_t1,adc_t2,adc_t3,adc_t4,adc_t5,adc_t6,
      adc_t7,adc_t8,adc_b0,adc_dg,adc_tg,adc_del_tg;
float peak_b0, tdc_b0;
float peak_dg, tdc_dg;
float peak_tg, tdc_tg;
int   i_b0_max, i_dg_max, i_tg_max;
float adc_dsc_gain[16] = {1.26,1.12,0.89,1.14,0.97,0.62,0.89,0.57,0.85,
                          1.0,1.0,1.0,1.0,1.0,1.0,1.0};
int i_count;  /* the only change WL */

  /* open banks of interest */
  if (bk_locate(pevent, "CADC", &cadc) == 0 ||
      bk_locate(pevent, "TTIM", &ttim) == 0 ||
      bk_locate(pevent, "TPRE", &tpre) == 0 ||
      bk_locate(pevent, "TPST", &tpst) == 0 ||  
      bk_locate(pevent, "TWID", &twid) == 0 ||
      bk_locate(pevent, "NHIT", &nhit) == 0 ||
      bk_locate(pevent, "PVET", &pvet) == 0 ||
      bk_locate(pevent, "MWPC", &mwpc) == 0 ||
      bk_locate(pevent, "ASUM", &asum) == 0 ||
      bk_locate(pevent, "CLMP", &clmp) == 0 ||
      bk_locate(pevent, "BEAM", &beam) == 0 ||
      bk_locate(pevent, "TBIT", &tbit) == 0 ||
      bk_locate(pevent, "TRIG", &trig) == 0 ||
      bk_locate(pevent, "BADC", &Badc) == 0 ||  
      bk_locate(pevent, "TRAK", &trak) == 0)
    /* return of one of these banks is not present */
    return ANA_CONTINUE;

  /* create temporary bank, to be filled  */
  bk_create(pevent, "TEMP", TID_STRUCT, &temp);
  temp->tmp1 = 0.f;
  temp->tmp2 = -1000.f;
  temp->tmp3 = 0.f;
  temp->tmp4 = 0.f;


 /* for tdc_offsets according to different trigger */      
 
 /*------prompt gamma -------------------*/
   if(tbit[4]==1 && pvet->pv_edep[0]<0.2f &&
      asum->csi_sum<200 &&
      mwpc->wc1_n_hits == 0 &&
      mwpc->wc2_n_hits == 0 
     )
     {
     /* fill histos with raw TDC value and channel */
     for(i_count=0; i_count<n_data_histo; i_count++)
       {
	 if(fabs((tim_histo[i_count].time-time_deg_normal))<10.f /*&& cadc[tim_histo[i_count].channel]>0*/)
       	   HF2(TDC_OFFSET_PROMPT, tim_histo[i_count].time, (float)(tim_histo[i_count].channel), 1.f);  
       }
     }

 /*-------- pienuH Gamma ----------------*/
 if(tbit[2]==1 && pvet->pv_edep[0]<0.2f &&
    asum->csi_sum<200 &&
    mwpc->wc1_n_hits == 0 &&
    mwpc->wc2_n_hits == 0 
   )
   {
    /* fill histos with raw TDC value and channel */
    for(i_count=0; i_count<n_data_histo; i_count++)
      {
	if((tim_histo[i_count].time-time_deg_normal) > 20.f && cadc[tim_histo[i_count].channel]>0)
	  {
	    HF2(TDC_OFFSET_PIENUH_G, tim_histo[i_count].time, (float)(tim_histo[i_count].channel), 1.f);
	  }
      }
   }
 /*---------------- pienuH e+ */
 if(tbit[2]==1 && 
    pvet->pv_edep[0]>0.2f &&
    pvet->pv_edep[0]<1.2  &&
    fabs(asum->csi_sum-70.)<10. &&
    mwpc->wc1_n_hits == 1 &&
    mwpc->wc2_n_hits == 1 
   )
   {
    /* fill histos with raw TDC value and channel */
    for(i_count=0; i_count<n_data_histo; i_count++)
      if((tim_histo[i_count].time-time_deg_normal) > 20.f && cadc[tim_histo[i_count].channel]>0)
	{
	  HF2(TDC_OFFSET_PIENUH_E, tim_histo[i_count].time, (float)(tim_histo[i_count].channel), 1.f);
	}
   }

  bk_close(pevent, temp+1); 
  return ANA_CONTINUE;
}

/*----------- eor routine ----------- */
static float yfloat[N_TDC][TDC_N_BINS_HISTO];
extern TDC_CALIBRATION_PARAM tdccalib_param;

INT histo_eor(INT run_number)
{

int    i, i_channel;
double x[TDC_N_BINS_HISTO], y[TDC_N_BINS_HISTO];
double param[3], error[3], chisqr, ymax;
char   str[80];
HNDLE  hDB, hkey;
FILE *f_final;
float tdc_offset[480];

  memset(tdc_offset, 0, sizeof(tdc_offset));

  run_counter++; 
  f_final = fopen("tdcoffset.dat", "w");

  fprintf(f_final, "[/Analyzer/Parameters/TDC calibration]\n");
  fprintf(f_final, "Correct Offsets = INT:2\n");
  fprintf(f_final, "Offset = FLOAT[480]:\n");
  
  cm_get_experiment_database(&hDB, NULL); 
  printf("run number = %d run_counter=%d\n", run_number, run_counter); 
 
  if (run_counter == N_RUNS)
    {
     run_counter = 0;
    /*-------------- for prompt ------------------------------------*/
    printf(" calculating tdcoffset fot trigger prompt\n");
    HUNPAK(TDC_OFFSET_PROMPT, yfloat[0][0], "HIST", 1);
    
    /*-analyze TDC histos ---*/

    for (i=0 ; i<TDC_N_BINS_HISTO ; i++)
      x[i] = (float) i/TDC_N_BINS_HISTO*(TDC_X_HIGH_HISTO-TDC_X_LOW_HISTO)+TDC_X_LOW_HISTO;

    /* loop through TDC channels */
    for (i_channel=CHN_B0_TDC; i_channel<=CHN_TGTH_TDC ; i_channel++)
      {
      /* move TDC values into array */

      ymax = 0;
      for (i=0 ; i<TDC_N_BINS_HISTO ; i++)
        {
        y[i] = yfloat[i_channel][i];
        if (y[i] > ymax && i>0)
          ymax = y[i];
        }
      if (ymax > 10)
        {
        chisqr = gauss_fit(x, y, TDC_N_BINS_HISTO, param, error, TRUE);
        if (chisqr != -1)
          {
	    tdc_offset[i_channel] =(float)param[1]; 
          }
        }
      }
    tdc_offset[CHN_T5_8_TDC] = 0.f;    /* can't get enough data */
      /*----------------- pienu gammas --------------------------------*/
    
    printf("calculating tdc_offsets for trigger PienuH gammas\n");
    HUNPAK(TDC_OFFSET_PIENUH_G, yfloat[0][0], "HIST", 1);

    /* analyze TDC histos*/

    /* loop through TDC channels */
    for (i_channel=0; i_channel<N_CSI ; i_channel++)
      {
      /* move TDC values into array */

      ymax = 0;
      for (i=0 ; i<TDC_N_BINS_HISTO ; i++)
        {
        y[i] = yfloat[i_channel][i];
        if (y[i] > ymax && i>0)
          ymax = y[i];
        }
      if (ymax > 10)
        {
        chisqr = gauss_fit(x, y, TDC_N_BINS_HISTO, param, error, TRUE);
        if (chisqr != -1)
          {
	    tdc_offset[i_channel] = (float)param[1];
          }
        }
      }
    /*--------------- pienu positrons ---------------------------------- */
    
    printf("calculating tdc_offsets for trigger PienuH positrons\n");
    HUNPAK(TDC_OFFSET_PIENUH_E, yfloat[0][0], "HIST", 1);

     /*analyze TDC histos */
    /* loop through TDC channels */
    for (i_channel=240; i_channel<280 ; i_channel++)  /* PVs */
      {
      /* move TDC values into array */

      ymax = 0;
      for (i=0 ; i<TDC_N_BINS_HISTO ; i++)
        {
        y[i] = yfloat[i_channel][i];
        if (y[i] > ymax && i>0)
          ymax = y[i];
        }
      if (ymax > 10)
        {
        chisqr = gauss_fit(x, y, TDC_N_BINS_HISTO, param, error, TRUE);
        if (chisqr != -1)
          {
	    tdc_offset[i_channel] = (float)param[1];  
          }
        }
      }
    /*------------ write final tdcoffsets into file tdcoffset.dat -----------*/
    for(i_channel=0; i_channel<480; i_channel++)
      fprintf(f_final, "[%d] %f\n", i_channel,(float)tdccalib_param.offset[i_channel]+tdc_offset[i_channel]);
 
  }
    fclose(f_final);
    return SUCCESS;
}