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

  Name:         HVCalib.cc

  Created by:   Luisa Sabrina Stark Schneebeli, April 2009
                sabrina.stark@phys.ethz.ch

  Contents:     Class reading the calibration table for DAC to HV values

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


#include <stdio.h>
#include <iostream>
#include <fstream>
#include <math.h>

#include "HVCalib.h"

#include "ProcessIO.h" // Must be not in HVCalib.h to avoid problem with declaring class ProcessIO

using namespace std;


HVCalib::HVCalib(class ProcessIO *hvConfig) {

  char calibfile[80], dLine[6000];
  int idac, idacOld, istep, ihv1000, ihvOld;
  bool first = true;
  float hv1000,dac;
  double intpart;

  Config = hvConfig;

  iDACMin = Config->DACMin;
  fHVCalibSlope = Config->fHVCalibSlope;
  NDACValues = Config->DACMax - Config->DACMin;
  NHVValues = 10*NDACValues;

  HVArray = new double*** [Config->NumHVBoards];
  DACArray = new int*** [Config->NumHVBoards];
  for(int i=0; i<Config->NumHVBoards; i++) {
    HVArray[i] = new double** [NUM_CHAINS];
    DACArray[i] = new int** [NUM_CHAINS];
    for(int j=0; j<NUM_CHAINS; j++) {
      HVArray[i][j] = new double* [NUM_CHANNELS];
      DACArray[i][j] = new int* [NUM_CHANNELS];
      for(int k=0; k<NUM_CHANNELS; k++) {
	HVArray[i][j][k] = new double [NDACValues];
	DACArray[i][j][k] = new int [NHVValues];
      }
    }
  }

  for(int i=0; i<Config->NumHVBoards; i++){
    for(int j=0; j<NUM_CHAINS; j++){
      sprintf(calibfile,"Calib/%s_%c%d.txt",hvConfig->fHVBoard[i]->BoardName,65+j,1); 
      ifstream fin(calibfile);

      while(fin){
	fin >> dac;
	idacOld = 0;
	if ((int) dac >= Config->DACMin) {
	  if (first){
	    iDACMin = (int) dac;
	    first = false;
	  }
	  idac = (int)dac - iDACMin;
	  istep = idac - idacOld;
	  for(int k=0; k<NUM_CHANNELS; k++){
 	    if (idac>=NDACValues) break;

	    fin>>HVArray[i][j][k][idac];

	    for (int l=1; l<istep;l++){
	      HVArray[i][j][k][idac-istep+l] = HVArray[i][j][k][idac-istep] + l*(HVArray[i][j][k][idac]-HVArray[i][j][k][idac-istep])/istep;
	      //	      printf(" HVArray: %d  %f \n", idac-istep+l, HVArray[i][j][0][idac-istep+l]);
	      //	      if (idac > 11240-iDACMin)
	      //		return;
	    }
	  }
	  idacOld = idac;
	}
	else{
	  fin.getline(dLine,6000,'\n'); // Skip low voltage region
	}
      }
    }
  }

  for(int i=0; i<Config->NumHVBoards; i++){
    for(int j=0; j<NUM_CHAINS; j++){
      for(int k=0; k<NUM_CHANNELS; k++){
	ihvOld = 0;
	first = true;
	for (int l=0; l<NDACValues; l++){
	  if (l == 0){
	    ihvOld = 0;
	    DACArray[i][j][k][0] = iDACMin;	    
	    //	    printf("first :  DACArray[i][j][k][0] = %d, ihvOld = %d, i %d, j %d, k %d\n",  DACArray[i][j][k][0], ihvOld,i,j,k);
	  }
	  hv1000 = (HVArray[i][j][k][l] - HVArray[i][j][k][0]) * 1000;
	  if (modf(hv1000,&intpart) >= 0.5){
	    ihv1000 = (int)((hv1000>=0)?ceil(hv1000):floor(hv1000));
	  }
	  else ihv1000 = (int)((hv1000<0)?ceil(hv1000):floor(hv1000));

	  if (ihv1000 <= ihvOld) continue;
	  else if (ihv1000 > ihvOld){
	    for (int m=0; m<(ihv1000-ihvOld); m++){
	      DACArray[i][j][k][ihvOld+m] = l+iDACMin;
	      //printf("HV = %d, HVArr = %f, DACArray = %d, m = %d. \n",ihvOld+m,HVArray[i][j][k][l],DACArray[i][j][k][ihvOld+m],m);
	    //	    printf("DACArray[i][j][k][0] = %d\n",  DACArray[i][j][k][0]);	    
	    }
	    ihvOld = ihv1000;
	  }
	}	
      }
    }
  }
}

/////////////////////////////////////////////////////////////////////////

HVCalib::~HVCalib() {

  for(int i=0; i<Config->NumHVBoards; i++) {
    for(int j=0; j<NUM_CHAINS; j++) {
      for(int k=0; k<NUM_CHANNELS; k++) {
	delete[] HVArray[i][j][k];
	delete[] DACArray[i][j][k];
      }
      delete[] HVArray[i][j];
      delete[] DACArray[i][j];
    }
    delete[] HVArray[i];
    delete[] DACArray[i];
  }
  delete[] HVArray;
  delete[] DACArray;
}

/////////////////////////////////////////////////////////////////////////

double HVCalib::DACToHV(int dac, int board, int chain, int channel) {
  
  if (dac == 0) return 0;
  
  if (dac < iDACMin){
    return HVArray[board][chain][channel][0] + (dac - iDACMin)*fHVCalibSlope;
  }
  else
    return HVArray[board][chain][channel][dac-iDACMin];

}
/////////////////////////////////////////////////////////////////////////

int HVCalib::HVToDAC(double hv, int board, int chain, int channel) {

  if (hv == 0) return 0;
  if (hv < HVArray[board][chain][channel][0]){
        return iDACMin + (int)((hv - HVArray[board][chain][channel][0])/fHVCalibSlope);
  }
  else{
    int ihv = 0;
    double intpart;
    double hv1000 = (hv-HVArray[board][chain][channel][0])*1000.;
    if (modf(hv1000,&intpart) >= 0.5){
      ihv = (int) (hv1000>=0 ? ceil(hv1000) : floor(hv1000));
    }
    else {
      ihv = (int) (hv1000<0 ? ceil(hv1000) : floor(hv1000));
    }
    if (ihv >= NHVValues) ihv = NHVValues;
    return DACArray[board][chain][channel][ihv];
  }
}