source: drsdaq/DRS/DRS.cc@ 68

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

  Name:         DRS.cc

  Created by:   Stefan Ritt, Matthias Schneebeli

  Modified by:  Sebastian Commichau, May-November 2008
                commichau@phys.ethz.ch

  Modification: This library works with:
                - Concurrent Technolgies VME single board PC (VP 315) 
                - Struck VME controller (SIS 3100) => faster!

  Contents:     Library functions for DRS board CMC card - requires
                DRS version 2 or 3 

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

#include "DRS.h"

#define DEBUG 0

// Minimal FPGA firmware version required for this library
#define REQUIRED_FIRMWARE_VERSION_DRS2 5268
#define REQUIRED_FIRMWARE_VERSION_DRS3 6981

#ifdef CT_VME
#define MEM_SEGMENT 0XA000 // Size of the memory segment allocated by each DRS board for BLT
#define BLT_TIMEOUT 1000   // Timeout for BLT [msec]
#endif

int drs_kbhit()
{
  int n;
  
  ioctl(0, FIONREAD, &n);
  return (n > 0);
}

static inline int getch() { return getchar(); }
inline void Sleep(useconds_t x) { usleep(x * 1000); }

// VME addresses

/* Assuming following DIP Switch settings:
   
SW1-1: 1 (off)       Use geographical addressing (1=left, 21=right)
SW1-2: 1 (off)       \
SW1-3: 1 (off)        >  VME_WINSIZE = 8MB, subwindow = 1MB
SW1-4: 0 (on)        /
SW1-5: 0 (on)        Reserved
SW1-6: 0 (on)        Reserved
SW1-7: 0 (on)        Reserved
SW1-8: 0 (on)        \
SW2-1: 0 (on)         |
SW2-2: 0 (on)         |
SW2-3: 0 (on)         |
SW2-4: 0 (on)         >  VME_ADDR_OFFSET = 0
SW2-5: 0 (on)         |
SW2-6: 0 (on)         |
SW2-7: 0 (on)         |
SW2-8: 0 (on)        /

which gives 
VME base address = SlotNo * VME_WINSIZE + VME_ADDR_OFFSET
= SlotNo * 0x80'0000 */

#define GEVPC_BASE_ADDR           0x00000000
#define GEVPC_WINSIZE               0x800000
#define GEVPC_USER_FPGA   (GEVPC_WINSIZE*2/8)
#define PMC1_OFFSET                  0x00000
#define PMC2_OFFSET                  0x80000
#define PMC_CTRL_OFFSET              0x00000  // All registers 32 bit!
#define PMC_STATUS_OFFSET            0x10000
#define PMC_FIFO_OFFSET              0x20000
#define PMC_RAM_OFFSET               0x40000

// DRS registers
#define T_CTRL                     1
#define T_STATUS                   2
#define T_RAM                      3
#define T_FIFO                     4
#define REG_CTRL             0x00000    // 32-bit control reg
#define REG_DAC_OFS          0x00004
#define REG_DAC0             0x00004
#define REG_DAC1             0x00006
#define REG_DAC2             0x00008
#define REG_DAC3             0x0000A
#define REG_DAC4             0x0000C
#define REG_DAC5             0x0000E
#define REG_DAC6             0x00010
#define REG_DAC7             0x00012
#define REG_CHANNEL_CONFIG   0x00014
#define REG_CHANNEL_SPAN     0x00016
#define REG_FREQ_SET_HI      0x00018
#define REG_FREQ_SET_LO      0x0001A
#define REG_TRIG_DELAY       0x0001C
#define REG_CALIB_TIMING     0x0001E
#define REG_STATUS           0x00000
#define REG_RDAC_OFS         0x0000A
#define REG_RDAC0            0x00004
#define REG_RDAC1            0x00006
#define REG_RDAC2            0x00008
#define REG_RDAC3            0x0000A
#define REG_RDAC4            0x0000C
#define REG_RDAC5            0x0000E
#define REG_RDAC6            0x00010
#define REG_RDAC7            0x00012
#define REG_EVENTS_IN_FIFO   0x00014
#define REG_EVENT_COUNT      0x00016
#define REG_FREQ1            0x00018
#define REG_FREQ2            0x0001A
#define REG_TEMPERATURE      0x0001C
#define REG_TRIGGER_BUS      0x0001E
#define REG_SERIAL_CMC       0x00020
#define REG_VERSION_FW       0x00022

using namespace std;


DRS::DRS() : fNumberOfBoards(0),
#ifdef STRUCK_VME
             fVMEInterface(0),
#endif
             First_VME_Slot(0), Last_VME_Slot(7)
{ }


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

DRS::~DRS() {
  
  int i;
  
  for (i = 0; i < fNumberOfBoards; i++) {
    delete fBoard[i];
  }
  
#ifdef CT_VME 
  if (!CloseVME()) 
    printf("VME connection closed\n");
  
  if (!CloseCMEM())
    printf("CMEM closed\n");
#endif

#ifdef STRUCK_VME
  if (fVMEInterface != NULL)
    if (mvme_close(fVMEInterface));
        printf("VME connection closed\n");

#endif

}

/*------------------------------------------------------------------*/
#ifdef CT_VME 
int DRS::OpenVME() {

  // Open VME connection
  if ((ErrorCode = VME_Open()) != VME_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return ErrorCode;
}
#endif
/*------------------------------------------------------------------*/

void DRS::InitialScan() {

  int index = 0;
  
  unsigned short Firmware, Serial, Temperature;

#ifdef CT_VME 
  
  unsigned int BoardAddress;
    
  if (!OpenVME()) {
    
    printf("VME connection opened\n");

    if (!OpenCMEM())
      printf("CMEM opened\n");
    else return;
    
    // Set master mapping input information
    MasterMap.vmebus_address    = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE; // Init VME board base address (derived from
                                                                           // the slot number => geographical addressing)

    MasterMap.window_size       = GEVPC_WINSIZE;                           // VME address window size
    
    MasterMap.address_modifier  = VME_A32;  
    MasterMap.options           = 0;
    
    
    // Check all VME slave slots 
    for (index = First_VME_Slot; index <= Last_VME_Slot; index++) {
      
      
      MasterMap.vmebus_address  = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE; // Update VME board base address

      
      if (DEBUG)
        printf("Checking VME slot %d (base address: 0X%08X)\n",index,MasterMap.vmebus_address);
                  
      MasterMapVME(&MasterMapping[index]);
      
      // **************************** Check PMC1 ****************************
      BoardAddress  = GEVPC_USER_FPGA;   // UsrFPGA base address
      BoardAddress += PMC1_OFFSET;       // PMC1 offset
      
      // Read firmware
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_VERSION_FW, &Firmware);
      
      // Read serial number
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_SERIAL_CMC, &Serial);

      // Read temperature register to see if CMC card is present 
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_TEMPERATURE, &Temperature);
      
      if (Firmware > 2400 && Firmware < 20000) {
        
        if (Temperature == 0XFFFF) {
          if (DEBUG)
            printf("No CMC board in slot %d\n", index);
        } else {
          if (DEBUG) {
            printf("Found CMC board in slot %d:\n", index);
            printf(" Firmware: %d\n",Firmware);
            printf(" Board serial nr.: %d\n",Serial);
            printf(" Temperature register: %d\n",Temperature);
          }
          
          fBoard[fNumberOfBoards] = new DRSBoard(MasterMapping[index], MasterMap.vmebus_address, BoardAddress, (index-2) << 1);
          
          if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
            fNumberOfBoards++;
          else
            if (DEBUG) 
              printf("Error: wrong firmware version: board has %d, required is %d\n",
                     fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                     fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
        }
      }

           
      // **************************** Check PMC2 ****************************
      BoardAddress  = GEVPC_USER_FPGA;   // UsrFPGA base address
      BoardAddress += PMC2_OFFSET;       // PMC2 offset
      
      // Read firmware
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_VERSION_FW, &Firmware);
      
      // Read serial number
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_SERIAL_CMC, &Serial);
      
      // Read temperature register to see if CMC card is present 
      VME_ReadFastUShort(MasterMapping[index], BoardAddress + PMC_STATUS_OFFSET + REG_TEMPERATURE, &Temperature);

      if (Firmware > 2400 && Firmware < 20000) {
        
        if (Temperature == 0XFFFF) {
          printf("No CMC board in slot %d\n", index);
        } else {
          if (DEBUG) {
            printf("Found CMC board in slot %d:\n", index);
            printf(" Firmware: %d\n",Firmware);
            printf(" Board serial nr.: %d\n",Serial);
            printf(" Temperature register: %d\n",Temperature);
          }
                  
          fBoard[fNumberOfBoards] = new DRSBoard(MasterMapping[index], MasterMap.vmebus_address, BoardAddress, ((index-2) << 1) | 1);
          
          if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
            fNumberOfBoards++;
          else
            if (DEBUG) 
              printf("Error: wrong firmware version: board has %d, required is %d\n",
                     fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                     fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
        }
      }
      else
        MasterUnMapVME(MasterMapping[index]);
    }
    if (DEBUG) printf("Found %d DRS boards in VME crate\n", fNumberOfBoards);
  } 
#endif
#ifdef STRUCK_VME

  int i = 0;

  mvme_addr_t Address;
  
  if (mvme_open(&fVMEInterface, 0) == MVME_SUCCESS) {
    
    printf("VME connection opened\n");

    mvme_set_am(fVMEInterface, MVME_AM_A32);
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D16);
    
    // Check all VME slave slots 
    for (index = 2; index <= 21; index++) {
      
      // **************************** Check PMC1 ****************************
      Address = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE;  // VME board base address
      Address += GEVPC_USER_FPGA;                         // UsrFPGA base address
      Address += PMC1_OFFSET;                             // PMC1 offset

      // Read firmware
      i = mvme_read(fVMEInterface, &Firmware, Address + PMC_STATUS_OFFSET + REG_VERSION_FW, 2);
      
      if (i == 2) {
        
        // Read serial number
        mvme_read(fVMEInterface, &Serial, Address + PMC_STATUS_OFFSET + REG_SERIAL_CMC, 2);
        
        // Read temperature register to see if CMC card is present 
        mvme_read(fVMEInterface, &Temperature, Address + PMC_STATUS_OFFSET + REG_TEMPERATURE, 2);
        
        if (Firmware > 2400 && Firmware < 20000) {
          
          if (Temperature == 0xFFFF) {
            printf("No CMC board in slot %d\n", index);
          } else {
            if (DEBUG) {
              printf("Found CMC board in slot %d:\n", index);
              printf(" Firmware: %d\n",Firmware);
              printf(" Board serial nr.: %d\n",Serial);
              printf(" Temperature register: %d\n",Temperature);
            }
            
            fBoard[fNumberOfBoards] = new DRSBoard(fVMEInterface, Address, (index-2) << 1);
            
            if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
              fNumberOfBoards++;
            else
              if (DEBUG) 
                printf("Error: wrong firmware version: board has %d, required is %d\n",
                       fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                       fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
            
          }
        }
      }
      
      // **************************** Check PMC2 ****************************
      Address = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE;  // VME board base address
      Address += GEVPC_USER_FPGA;                         // UsrFPGA base address
      Address += PMC2_OFFSET;                             // PMC2 offset
      
      // Read firmware
      i = mvme_read(fVMEInterface, &Firmware, Address + PMC_STATUS_OFFSET + REG_VERSION_FW, 2);

      if (i == 2) {

        // Read serial number
        mvme_read(fVMEInterface, &Serial, Address + PMC_STATUS_OFFSET + REG_SERIAL_CMC, 2);
        
        // Read temperature register to see if CMC card is present 
        mvme_read(fVMEInterface, &Temperature, Address + PMC_STATUS_OFFSET + REG_TEMPERATURE, 2);
        
        if (Firmware > 2400 && Firmware < 20000) {
          
          if (Temperature == 0xFFFF) {
            printf("No CMC board in slot %d\n", index);
          } else {
            if (DEBUG) {
              printf("Found CMC board in slot %d:\n", index);
              printf(" Firmware: %d\n",Firmware);
              printf(" Board serial nr.: %d\n",Serial);
              printf(" Temperature register: %d\n",Temperature);
            }
            
            fBoard[fNumberOfBoards] = new DRSBoard(fVMEInterface, Address, ((index-2) << 1) | 1);
        
            if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
              fNumberOfBoards++;
            else
              if (DEBUG) 
                printf("Error: wrong firmware version: board has %d, required is %d\n",
                       fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                       fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
          }
        }
      }
    }
    if (DEBUG) printf("Found %d DRS boards in VME crate\n", fNumberOfBoards);
  }   
#endif
  else printf("Error: cannot access VME crate, check driver, power and connection\n");

}

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

#ifdef CT_VME 
int DRS::MasterMapVME(int *MMap) {
  
  // Do master mapping  
  if (ErrorCode = VME_MasterMap(&MasterMap, MMap)) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return(ErrorCode);
}
#endif

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

#ifdef CT_VME 
int DRS::MasterUnMapVME(int MMap) {
  
  // Delete master mapping
  if (ErrorCode = VME_MasterUnmap(MMap)) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return(ErrorCode);  
}
#endif

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

#ifdef CT_VME 
int DRS::CloseVME() {
  
  // Close VME connection
  if ((ErrorCode = VME_Close()) != VME_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return ErrorCode;
}
#endif

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

#ifdef CT_VME 
int DRS::OpenCMEM() {
  
  // Open CMEM package
  if ((ErrorCode = CMEM_Open()) != CMEM_RCC_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return ErrorCode;
}
#endif

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

#ifdef CT_VME 
int DRS::CloseCMEM() {
  
  // Close CMEM package
  if ((ErrorCode = CMEM_Close()) != CMEM_RCC_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return ErrorCode;
}
#endif

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

#ifdef CT_VME
DRSBoard::DRSBoard(int MasterMapping, unsigned int BaseAddress, unsigned int BoardAddress, int SlotNumber)
#else
DRSBoard::DRSBoard(MVME_INTERFACE *MVME_Interface, mvme_addr_t BaseAddress, int SlotNumber)
#endif
  :fDAC_COFSA(0)
  ,fDAC_COFSB(0)
  ,fDAC_DRA(0)
  ,fDAC_DSA(0)
  ,fDAC_TLEVEL(0)
  ,fDAC_ACALIB(0)
  ,fDAC_DSB(0)
  ,fDAC_DRB(0)
  ,fDAC_COFS(0)
  ,fDAC_ADCOFS(0)
  ,fDAC_CLKOFS(0)
  ,fDAC_ROFS_1(0)
  ,fDAC_ROFS_2(0)
  ,fDAC_INOFS(0)
  ,fDAC_BIAS(0)
  ,fRequiredFirmwareVersion(0)
  ,fFirmwareVersion(0)
  ,fChipVersion(0)
  ,fBoardVersion(0)
  ,fCMCSerialNumber(0)
  ,fCtrlBits(0)
  ,fNumberOfReadoutChannels(0)
  ,fExternalClockFrequency(0)
  ,fBaseAddress(BaseAddress)
#ifdef CT_VME
  ,fBoardAddress(BoardAddress)
  ,fMasterMapping(MasterMapping)
  ,fSlotNumber(SlotNumber)
#endif
#ifdef STRUCK_VME
  ,fVMEInterface(MVME_Interface)
  ,fSlotNumber(SlotNumber)
#endif
  ,fFrequency(0)
  ,fDominoMode(0)
  ,fReadoutMode(0)
  ,fTriggerEnable(0)
  ,fDelayedStart(0)
  ,fTriggerCell(0)
  ,fMaxChips(0)
  ,fResponseCalibration(0)
  ,fTimeData(0)
  ,fNumberOfTimeData(0)
  ,fDebug(0)
  ,fTriggerStartBin(0)
  ,kRotateWave(0)
{
  ConstructBoard();
}


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

DRSBoard::~DRSBoard()
{
  // Response Calibration
  delete fResponseCalibration;
  
  int i;
  // Time Calibration
  for (i = 0; i < fNumberOfTimeData; i++) {
    delete fTimeData[i];
  }
  delete[]fTimeData;
  
#ifdef CT_VME 
  if (!FreeSegmentCMEM(CMEM_SegIdentifier) && DEBUG)
    printf("Memory segment %d (board #%d) freed\n",CMEM_SegIdentifier,fCMCSerialNumber);
#endif
}

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

void DRSBoard::ConstructBoard()
{
  fDebug = 0;
  fDominoMode = 1;
  fReadoutMode = 0;
  fTriggerEnable = 0;
  fCtrlBits = 0;
  fNumberOfReadoutChannels = 10;
  fExternalClockFrequency = 1000. / 30.;
  
  for (int i = 0; i < kNumberOfChips * kNumberOfChannels; i++)
    fWaveTransferred[i] = false;
  
  strcpy(fCalibDirectory, ".");
  
#ifdef CT_VME
  if (DEBUG) {
    printf("Base address:  0X%08X\n",fBaseAddress);
    printf("Board address: 0X%08X\n",fBoardAddress);
    printf("0X%08X\n",fBaseAddress+fBoardAddress);
  }
#endif
  ReadSerialNumber();
  
  // Check for required firmware version 
  if (!HasCorrectFirmware())
    return;
  
  if (DEBUG)
    printf("Board version: %d\n",fBoardVersion);

  if (fBoardVersion == 1) {
    
    fDAC_COFSA  = 0;
    fDAC_COFSB  = 1;
    fDAC_DRA    = 2;
    fDAC_DSA    = 3;
    fDAC_TLEVEL = 4;
    fDAC_ACALIB = 5;
    fDAC_DSB    = 6;
    fDAC_DRB    = 7;
    
  } else if (fBoardVersion == 2 || fBoardVersion == 3) {
    
    fDAC_COFS   = 0;
    fDAC_DSA    = 1;
    fDAC_DSB    = 2;
    fDAC_TLEVEL = 3;
    fDAC_CLKOFS = 5;
    fDAC_ACALIB = 6;
    fDAC_ADCOFS = 7;
    
  } else if (fBoardVersion == 4) {
    
    fDAC_ROFS_1 = 0;
    fDAC_DSA    = 1;
    fDAC_DSB    = 2;
    fDAC_ROFS_2 = 3;
    fDAC_BIAS   = 4;
    fDAC_ADCOFS = 7;
    fDAC_INOFS  = 5;
    fDAC_ACALIB = 6;
  }
  
  // Response Calibration
  fResponseCalibration = new ResponseCalibration(this);
  
  // Time Calibration
  fTimeData = new DRSBoard::TimeData *[kNumberOfChips];
  fNumberOfTimeData = 0;

#ifdef CT_VME 
  // Allocate contiguous memory for BLT
  AllocateSegmentCMEM(MEM_SEGMENT,&CMEM_SegIdentifier);
  if (DEBUG)
    printf("Memory segment %d (board #%d) allocated\n",CMEM_SegIdentifier,fCMCSerialNumber);

  AssignPhysicalSegAddressCMEM(CMEM_SegIdentifier, &PCIAddress);
  AssignVirtualSegAddressCMEM(CMEM_SegIdentifier, &VirtualAddress);
  if (DEBUG)
    printf("Physical address: 0X%08X, virtual address: 0X%08X\n", (unsigned int)PCIAddress,(unsigned int)VirtualAddress);
#endif

}

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

void DRSBoard::PrintBinary32(unsigned int i) {

  int k;

  for (k=31;k>=0;k--)
    if ((i & (1 << k)) !=0) {
        if ((k)%8)
        printf("1");
        else
        printf("1 ");
    }
    else {
        if ((k)%8)
        printf("0");
        else
        printf("0 ");
    }
  printf("\n");
}

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

long int DRSBoard::GetMicroSeconds() {

  struct tm * timeinfo;
  time_t rawtime;
  
  struct timezone tz;
  struct timeval actual_time;   // Actual time 
  
  gettimeofday(&actual_time, &tz);
  
  time(&rawtime);
  
  timeinfo = gmtime(&rawtime);  // Get UTC (or GMT timezone).
  
  return (timeinfo->tm_hour*3600 + timeinfo->tm_min*60 + timeinfo->tm_sec)*1000000 + actual_time.tv_usec;

}


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

int DRSBoard::TestRead(unsigned int n, int type) {

  float delta;
  float mbps;

  int errors = 0;
  
  const int size = 0X10000; // bytes
  
#ifdef STRUCK_VME 

  long int ts1, ts2;
  
  unsigned int Address = fBaseAddress + PMC_RAM_OFFSET;

  int read = 0, i;

  unsigned char data[size];
    
  printf("**************************************************\n");

  if (type==0) {
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D32);
    mvme_set_blt(fVMEInterface, MVME_BLT_BLT32); 
    printf(" Mode: VMEbus A32/D32 DMA read [64 kB]\n");
  }
  else if (type==1) {
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D32);
    mvme_set_blt(fVMEInterface,MVME_BLT_MBLT64);
    printf(" Mode: VMEbus A32/D64 DMA read [64 kB]\n");
  }
  else if (type==2) {
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D64);
    mvme_set_blt(fVMEInterface, MVME_BLT_2EVME); 
    printf(" Mode: VMEbus A32/D64 2eVME read [64 kB]\n");
  }

  ts1 = GetMicroSeconds();

  for (unsigned int j=0;j<n;j++) {
    read = mvme_read(fVMEInterface, static_cast<unsigned char*>(data), Address, size);
    
    while (read != size) {

      errors++;
      printf("Only read %d out of %d, retry with %d: ", read, size, size-read);
      i = mvme_read(fVMEInterface, static_cast<unsigned char*>(data) + read/4, Address + read, size - read);
      printf("read %d\n", i);
      if (i == 0) {
        printf("Error reading VME\n");
        return read;
      }
      read += i;
    }
  }

  ts2 = GetMicroSeconds();

  delta = (ts2 - ts1)/1000000.;

  mbps  = n * read * 1.0/(delta * 1024. * 1024.);

  printf(" %d BLT(s) finished...\n",n);
  
  if (!errors)
    printf(" %d errors... success!\n",errors);
  else
    printf(" %d errors...\n",errors);
  
  printf(" Duration: %.3f s\n", delta);
  printf(" Rate: %.3f MB/s\n", mbps);

  printf("**************************************************\n");
  
#endif
#ifdef CT_VME 
  tstamp ts1, ts2;

  unsigned int ret;

  ret = ts_open(1, TS_DUMMY);
  if (ret)
  {
    rcc_error_print(stdout, ret);
    exit(-2);
  }

  printf("**************************************************\n");

  // Assign fields for BLT
  BLT_List.number_of_items                       = 1;
  BLT_List.list_of_items[0].vmebus_address       = fBaseAddress + fBoardAddress + PMC_RAM_OFFSET;
  BLT_List.list_of_items[0].system_iobus_address = PCIAddress;
  BLT_List.list_of_items[0].size_requested       = size;

  if (type==0) {
    BLT_List.list_of_items[0].control_word         = VME_DMA_D32R | VME_A32; 
    printf(" Mode: VMEbus A32/D32 DMA read [64 kB]\n");
  }
  else if (type==1) {
    BLT_List.list_of_items[0].control_word         = VME_DMA_D64R | VME_A32; 
    printf(" Mode: VMEbus A32/D64 DMA read [64 kB]\n");
  }
  else if (type==2) {
    printf("Warning: current interface does not support 2exx transfer mode!\n");
    return 0;
  }

  printf("  VMEbus address:   0X%08X\n", BLT_List.list_of_items[0].vmebus_address);   
  printf(" Contiguous buffer:\n");
  printf("  Physical address: 0X%08X\n", (unsigned int)PCIAddress); 
  printf("  Virtual address:  0X%08X\n", (unsigned int)VirtualAddress); 
  
  ts_clock(&ts1);

  for (unsigned int i=0;i<n;i++)
    // Perfom BLT
    if ((ErrorCode = VME_BlockTransfer(&BLT_List,BLT_TIMEOUT)) != VME_SUCCESS) {
      VME_ErrorString(ErrorCode,ErrorString);
      printf(" VME (BLT): %s\n",ErrorString);
      
      errors++;
      return -1;
    }   

  ts_clock(&ts2);
  delta = ts_duration(ts1, ts2);
  mbps = n * size * 1.0 / (delta * 1024. * 1024.);
    
  printf(" %d BLT(s) finished...\n",n);

  if (!errors)
    printf(" %d errors... success!\n",errors);
  else
    printf(" %d errors...\n",errors);

  printf(" Duration: %.3f s\n", delta);
  printf(" Rate: %.3f MB/s\n", mbps);

  ret = ts_close(TS_DUMMY);
  if (ret)
  {
    rcc_error_print(stdout, ret);
    exit(-2);
  }

  printf("**************************************************\n");
#endif

  return 0;
}

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

void DRSBoard::ReadSerialNumber()
{
  unsigned char buffer[2];
  char str[80];
  int number;
  
  Read(T_STATUS, buffer, REG_VERSION_FW, 2);
  fFirmwareVersion = (static_cast<int>(buffer[1]) << 8) + buffer[0];
  
  // Retrieve board serial number
  Read(T_STATUS, buffer, REG_SERIAL_CMC, 2);
  number = (static_cast<int>(buffer[1]) << 8) + buffer[0];
  fCMCSerialNumber = number;
  
  // Determine DRS chip number from serial number
  fChipVersion = number < 1000 ? 2 : 3;
  
  if (number == 0xFFFF) {
    printf("Found new mezzanine board. Please select DRS version (2/[3]): ");
    fgets(str, sizeof(str), stdin);
    if (atoi(str) == 2)
      fChipVersion = 2;
    else
      fChipVersion = 3;
  }
  
  // Retrieve firmware version
  if (fChipVersion == 2)
    fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS2;
  if (fChipVersion == 3)
    fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS3;
  
  
  // Determine board version from serial number
  if (number >= 1000)
    fBoardVersion = 4;
  else if (number >= 100)
    fBoardVersion = 3;
  else if (number > 0)
    fBoardVersion = 2;
  else {
    fBoardVersion = 4;
    fChipVersion = 3;
    fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS3;
  }
  
  //Fixme (SCC 03032008): change function FlashEEPROM accordingly!
  //fChipVersion = 2;
  //fBoardVersion = 3;
  //fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS2;
}

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

bool DRSBoard::HasCorrectFirmware()
{
  // Check firmware version 
  return (fFirmwareVersion >= fRequiredFirmwareVersion);
}

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

int DRSBoard::Write(int type, unsigned int addr, void *data, int size)
{

  // Generic write function

#ifdef CT_VME 
  if (size > MEM_SEGMENT)
    size = MEM_SEGMENT;
  
  if (type == T_CTRL)
    addr += PMC_CTRL_OFFSET;
  else if (type == T_STATUS)
    addr += PMC_STATUS_OFFSET;
  else if (type == T_RAM)
    addr += PMC_RAM_OFFSET;
  
  if (size == 1) {
    
    VME_WriteSafeUChar(fMasterMapping, fBoardAddress + addr, *(static_cast<unsigned char*>(data)));
    
  } else if (size == 2) {
    
    VME_WriteSafeUShort(fMasterMapping, fBoardAddress + addr, *(static_cast<unsigned short*>(data)));
    
  } else if (size == 4) {
    
    VME_WriteSafeUInt(fMasterMapping, fBoardAddress + addr, *(static_cast<unsigned int*>(data)));
    
  } else {
    
    // Copy contiguous block of memory starting from VirtualAddress
    memcpy((unsigned long*)((unsigned long)VirtualAddress),(unsigned long*)data,size);

    // Assign fields for BLT
    BLT_List.list_of_items[0].vmebus_address       = fBaseAddress + fBoardAddress + PMC_RAM_OFFSET;
    BLT_List.list_of_items[0].system_iobus_address = PCIAddress;
    BLT_List.list_of_items[0].size_requested       = MEM_SEGMENT; 
    BLT_List.list_of_items[0].control_word         = VME_DMA_D64W | VME_A32;
   
    BLT_List.number_of_items = 1;

    // Perfom BLT
    if ((ErrorCode = VME_BlockTransfer(&BLT_List,BLT_TIMEOUT)) != VME_SUCCESS) {
      VME_ErrorString(ErrorCode,ErrorString);
      printf("VME (BLT): %s\n",ErrorString);
    }   
  }

#endif
#ifdef STRUCK_VME  
  
  unsigned int base_addr;
  
  base_addr = fBaseAddress;
  
  if (type == T_CTRL)
    base_addr += PMC_CTRL_OFFSET;
  else if (type == T_STATUS)
    base_addr += PMC_STATUS_OFFSET;
  else if (type == T_RAM)
    base_addr += PMC_RAM_OFFSET;
  
  if (size == 1) {
    // 8-bit write access
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D8);
    mvme_write(fVMEInterface, base_addr + addr, static_cast<mvme_locaddr_t*>(data), size);
  } else if (size == 2) {
    // 16-bit write access
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D16);
    mvme_write(fVMEInterface, base_addr + addr, static_cast<mvme_locaddr_t*>(data), size);
  } else {
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D32);
    
    // As long as no block transfer is supported, do pseudo block transfer
    mvme_set_blt(fVMEInterface, MVME_BLT_NONE);
    
    mvme_write(fVMEInterface, base_addr + addr, static_cast<mvme_locaddr_t*>(data), size);
  }
  
#endif

  return size;
}

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

int DRSBoard::Read(int type, void *data, unsigned int addr, int size)
{

  // Generic read function
  
#ifdef CT_VME

  if (size > MEM_SEGMENT)
    size = MEM_SEGMENT;

  if (type == T_CTRL)
    addr += PMC_CTRL_OFFSET;
  else if (type == T_STATUS)
    addr += PMC_STATUS_OFFSET;
  else if (type == T_RAM) {
    addr += PMC_RAM_OFFSET;
  }
  else if (type == T_FIFO) {
    addr += PMC_FIFO_OFFSET;
  }

  if (size == 1) {
    
    VME_ReadFastUChar(fMasterMapping, fBoardAddress + addr, static_cast<unsigned char*>(data)); 
  
  } else if (size == 2) {
    
    VME_ReadFastUShort(fMasterMapping, fBoardAddress + addr, static_cast<unsigned short*>(data)); 
    
  } else if (size == 4) {
    
    VME_ReadFastUInt(fMasterMapping, fBoardAddress + addr, static_cast<unsigned int*>(data)); 

  } else {

    // Assign fields for BLT
    BLT_List.list_of_items[0].vmebus_address       = fBaseAddress + fBoardAddress + PMC_RAM_OFFSET;
    BLT_List.list_of_items[0].system_iobus_address = PCIAddress;
    BLT_List.list_of_items[0].size_requested       = MEM_SEGMENT; 
    BLT_List.list_of_items[0].control_word         = VME_DMA_D64R | VME_A32;
   
    BLT_List.number_of_items = 1;

    // Perfom BLT
    if ((ErrorCode = VME_BlockTransfer(&BLT_List,BLT_TIMEOUT)) != VME_SUCCESS) {
      VME_ErrorString(ErrorCode,ErrorString);
      printf("VME (BLT): %s\n",ErrorString);
    }   
    
    // Copy contiguous block of memory starting from VirtualAddress
    memcpy((unsigned long*)data,(unsigned long*)((unsigned long)VirtualAddress),size);

    // Do pseudo BLT
    /*for (i = 0; i < size; i += 4) 
      VME_ReadFastUInt(fMasterMapping, fBoardAddress + addr + i, (unsigned int*)((unsigned char*)data + i));
    
    // Decode data
    printf("pseudo BLT: %d %d\n",((*((unsigned char*)data + 1) & 0x0f) << 8) | *((unsigned char*)data), 
                                 ((*((unsigned char*)data + 2)) << 4) | ((*((unsigned char*)data + 1)) >> 4) );

    printf("pseudo BLT: %d %d\n",((*((unsigned char*)data + 5) & 0x0f) << 8) | *((unsigned char*)data + 4), 
                                 ((*((unsigned char*)data + 6)) << 4) | ((*((unsigned char*)data + 5)) >> 4) );
    */
  }
  return size;

#endif
#ifdef STRUCK_VME

  unsigned int base_addr;
  int n, i;
  
  base_addr = fBaseAddress;
  
  if (type == T_CTRL)
    base_addr += PMC_CTRL_OFFSET;
  else if (type == T_STATUS)
    base_addr += PMC_STATUS_OFFSET;
  else if (type == T_RAM)
    base_addr += PMC_RAM_OFFSET;
  else if (type == T_FIFO)
    base_addr += PMC_FIFO_OFFSET;
  
  mvme_set_dmode(fVMEInterface, MVME_DMODE_D32);
  
  n = 0;
  if (size == 1) {
    // 8-bit read access
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D8);
    n = mvme_read(fVMEInterface, static_cast<mvme_locaddr_t*>(data), base_addr + addr, size);
  } else if (size == 2) {
    // 16-bit read access
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D16);
    n = mvme_read(fVMEInterface, static_cast<mvme_locaddr_t*>(data), base_addr + addr, size);
  } else {
    mvme_set_dmode(fVMEInterface, MVME_DMODE_D32);
  
    //mvme_set_blt(fVMEInterface, MVME_BLT_NONE); // Pseudo block transfer
    mvme_set_blt(fVMEInterface, MVME_BLT_2EVME);  // Use 2eVME if implemented
    //mvme_set_blt(fVMEInterface, MVME_BLT_MBLT64);  // Use MBLT64 if implemented
    //mvme_set_blt(fVMEInterface, MVME_BLT_2ESST);  // Use 2eSST if implemented

    n = mvme_read(fVMEInterface, static_cast<mvme_locaddr_t*>(data), base_addr + addr, size);
    while (n != size) {
      printf("Only read %d out of %d, retry with %d: ", n, size, size-n);
      i = mvme_read(fVMEInterface, static_cast<mvme_locaddr_t*>(data)+n/4, base_addr + addr + n, size-n);
      printf("read %d\n", i);
      if (i == 0) {
        printf("Error reading VME\n");
        return n;
      }
      n += i;
    }
    
    //for (i = 0; i < size; i += 4)
    //mvme_read(fVMEInterface, (mvme_locaddr_t *)((char *)data+i), base_addr + addr+i, 4);
  }
  return n;

#endif

}


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

#ifdef CT_VME 
int DRSBoard::AllocateSegmentCMEM(unsigned int SegSize, int *CMEM_SegIdentifier) {
  
  if ((ErrorCode = CMEM_SegmentAllocate(SegSize, "DMA_BUFFER", CMEM_SegIdentifier)) != CMEM_RCC_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }

  return ErrorCode;
}
#endif

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

#ifdef CT_VME 
int DRSBoard::AssignPhysicalSegAddressCMEM(int CMEM_SegIdentifier, unsigned long* PCIAddress) {
  
  if ((ErrorCode = CMEM_SegmentPhysicalAddress(CMEM_SegIdentifier, PCIAddress)) != CMEM_RCC_SUCCESS) {
    
    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);

  }

  return ErrorCode;
}
#endif

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

#ifdef CT_VME 
int DRSBoard::AssignVirtualSegAddressCMEM(int CMEM_SegIdentifier, unsigned long* VirtualAddress) {

  if ((ErrorCode = CMEM_SegmentVirtualAddress(CMEM_SegIdentifier, VirtualAddress)) != CMEM_RCC_SUCCESS) {

    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);

  }

  return ErrorCode;
}
#endif

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

#ifdef CT_VME 
 int DRSBoard::FreeSegmentCMEM(int CMEM_SegIdentifier) {
  
  // Free memory segment
  if ((ErrorCode = CMEM_SegmentFree(CMEM_SegIdentifier)) != CMEM_RCC_SUCCESS) {

    VME_ErrorString(ErrorCode,ErrorString);
    
    printf("Error: %s\n",ErrorString);
  }
  
  return ErrorCode;
 }
#endif

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

void DRSBoard::SetLED(int state)
{
  // Set LED state
  if (state)
    fCtrlBits |= BIT_LED;
  else
    fCtrlBits &= ~BIT_LED;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
}

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

void DRSBoard::SetChannelConfig(int firstChannel, int lastChannel, int nConfigChannels)
{
  // Set number of channels
  unsigned short d;
  
  if (lastChannel < 0 || lastChannel > 10) {
    printf("Invalid number of channels: %d (must be between 0 and 10)\n", lastChannel);
    return;
  }
  
  if (fChipVersion == 2) {
    // Register must contain last channel to read out starting from 9
    d = 9 - lastChannel;
    Write(T_CTRL, REG_CHANNEL_CONFIG, &d, 2);
  } else if (fChipVersion == 3) {
    // Upper four bits of register must contain last channel to read out starting from 9
    d = (firstChannel << 4) | lastChannel;
    Write(T_CTRL, REG_CHANNEL_SPAN, &d, 2);
    
    // Set bit pattern for write shift register
    switch (nConfigChannels ) {
    case  1: d = 0x001; break;
    case  2: d = 0x041; break;
    case  3: d = 0x111; break;
    case  4: d = 0x249; break;
    case  6: d = 0x555; break;
    case 12: d = 0xFFF; break;
    default:
      printf("Invalid channel configuration\n");
    }
    Write(T_CTRL, REG_CHANNEL_CONFIG, &d, 2);
  }
  
  fNumberOfReadoutChannels = lastChannel - firstChannel + 1;
}

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

void DRSBoard::SetNumberOfChannels(int nChannels)
{
  SetChannelConfig(0, nChannels-1, 12);
}

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

int DRSBoard::SetDAC(unsigned char channel, double value)
{
  // Set DAC value
  unsigned short d;
  
  // Normalize to 2.5V for 16 bit 
  d = static_cast<unsigned short>(value / 2.5 * 0xFFFF + 0.5);
  
  Write(T_CTRL, REG_DAC_OFS + (channel * 2), &d, 2);
  return 1;
}

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

int DRSBoard::ReadDAC(unsigned char channel, double *value)
{
  // Readback DAC value from control register
  unsigned char buffer[2];
  
  // Map 0->1, 1->0, 2->3, 3->2, etc. 
  // ofs = channel + 1 - 2*(channel % 2);
  
  Read(T_CTRL, buffer, REG_DAC_OFS + (channel * 2), 2);
  
  // Normalize to 2.5V for 16 bit 
  *value = 2.5 * (buffer[0] + (buffer[1] << 8)) / 0xFFFF;
  
  return 1;
}

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

int DRSBoard::GetRegulationDAC(double *value)
{
  // Get DAC value from status register (-> freq. regulation)
  unsigned char buffer[2];
  
  if (fBoardVersion == 1)
    Read(T_STATUS, buffer, REG_RDAC3, 2);
  else if (fBoardVersion == 2 || fBoardVersion == 3 || fBoardVersion == 4)
    Read(T_STATUS, buffer, REG_RDAC1, 2);
  
  // Normalize to 2.5V for 16 bit 
  *value = 2.5 * (buffer[0] + (buffer[1] << 8)) / 0xFFFF;
  
  return 1;
}

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

int DRSBoard::StartDomino()
{
  // Start domino sampling
  fCtrlBits |= BIT_START_TRIG;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  fCtrlBits &= ~BIT_START_TRIG;
  
  return 1;
}

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

int DRSBoard::Reinit()
{
  // Stop domino sampling
  // Reset readout state machine
  // Reset FIFO counters
  fCtrlBits |= BIT_REINIT_TRIG;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  fCtrlBits &= ~BIT_REINIT_TRIG;
  
  return 1;
}

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

int DRSBoard::Init()
{
  // Reinitialize
  fCtrlBits |= BIT_REINIT_TRIG;        // Reset readout state machine
  fCtrlBits &= ~BIT_FREQ_AUTO_ADJ;     // Turn auto. freq regul. off
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  fCtrlBits &= ~BIT_REINIT_TRIG;
  
  // Set default DAC voltages
  
  if (fBoardVersion == 1) {
    // Set max. domino speed
    SetDAC(fDAC_DRA, 2.5);
    SetDAC(fDAC_DSA, 2.5);
    SetDAC(fDAC_DRB, 2.5);
    SetDAC(fDAC_DSB, 2.5);
    // Set readout offset
    SetDAC(fDAC_COFSA, 0.9);
    SetDAC(fDAC_COFSB, 0.9);
    SetDAC(fDAC_TLEVEL, 1.7);
  } else if (fBoardVersion == 2 || fBoardVersion == 3) {
    // Set max. domino speed
    SetDAC(fDAC_DSA, 2.5);
    SetDAC(fDAC_DSB, 2.5);
    
    // Set readout offset
    SetDAC(fDAC_COFS, 0.9);
    SetDAC(fDAC_TLEVEL, 1.7);
    SetDAC(fDAC_ADCOFS, 1.7);  // 1.7 for DC coupling, 1.25 for AC
    SetDAC(fDAC_CLKOFS, 1);
  } else if (fBoardVersion == 4) {
    // Set max. domino speed
    SetDAC(fDAC_DSA, 2.5);
    SetDAC(fDAC_DSB, 2.5);
    
    // Set readout offset
    SetDAC(fDAC_ROFS_1, 1.25);   // LVDS level
    SetDAC(fDAC_ROFS_2, 0.85);   // Linear range  0.1V ... 1.1V
    
    SetDAC(fDAC_ADCOFS, 1.25);
    SetDAC(fDAC_ACALIB, 0.5);
    SetDAC(fDAC_INOFS, 0.6);
    SetDAC(fDAC_BIAS, 0.70);     // A bit above the internal bias of 0.68V
  }
  
  // Set default number of channels per chip 
  SetChannelConfig(0, fNumberOfReadoutChannels-1, 12);
  SetDominoMode(fDominoMode);
  SetReadoutMode(fReadoutMode);
  EnableTrigger(fTriggerEnable);
  
  // Disable calibration 
  EnableAcal(0, 0);
  SetCalibTiming(0, 0);
  
  return 1;
}

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

int DRSBoard::SetDominoMode(unsigned char mode)
{
  // Set domino mode
  // mode == 0: single sweep
  // mode == 1: run continously
  //
  fDominoMode = mode;
  if (mode)
    fCtrlBits |= BIT_DMODE;
  else
    fCtrlBits &= ~BIT_DMODE;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::SetDominoActive(unsigned char mode)
{
  // Set domino activity
  // mode == 0: stop during readout
  // mode == 1: keep domino wave running
  //
  fDominoMode = mode;
  if (mode)
    fCtrlBits |= BIT_DACTIVE;
  else
    fCtrlBits &= ~BIT_DACTIVE;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::SetReadoutMode(unsigned char mode)
{
  // Set readout mode
  // mode == 0: start from first bin
  // mode == 1: start from domino stop
  //
  fReadoutMode = mode;
  if (mode)
    fCtrlBits |= BIT_READOUT_MODE;
  else
    fCtrlBits &= ~BIT_READOUT_MODE;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::SoftTrigger(void)
{
  // Send a software trigger
  fCtrlBits |= BIT_SOFT_TRIG;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  fCtrlBits &= ~BIT_SOFT_TRIG;
  
  return 1;
}

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

int DRSBoard::EnableTrigger(int mode)
{
  // Enable external trigger
  fTriggerEnable = mode;
  if (mode)
    fCtrlBits |= BIT_ENABLE_TRIGGER;
  else
    fCtrlBits &= ~BIT_ENABLE_TRIGGER;
  
  if (mode == 2)
    fCtrlBits |= BIT_TRIGGER_DELAYED;
  else
    fCtrlBits &= ~BIT_TRIGGER_DELAYED;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::SetDelayedStart(int flag)
{
  // Enable external trigger
  fDelayedStart = flag;
  if (flag)
    fCtrlBits |= BIT_DELAYED_START;
  else
    fCtrlBits &= ~BIT_DELAYED_START;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::IsBusy()
{
  // Get running flag
  unsigned int status;
  
  Read(T_STATUS, &status, REG_STATUS, 4);
  return (status & BIT_RUNNING) > 0;
}

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

int DRSBoard::IsNewFreq(unsigned char chipIndex)
{
  unsigned int status;
  
  Read(T_STATUS, &status, REG_STATUS, 4);
  if (chipIndex == 0)
    return (status & BIT_NEW_FREQ1) > 0;
  return (status & BIT_NEW_FREQ2) > 0;
}

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

int DRSBoard::ReadFrequency(unsigned char chipIndex, double *f)
{
  // Read domino sampling frequency
  unsigned char buffer[2];
  
  if (chipIndex == 0)
    Read(T_STATUS, buffer, REG_FREQ1, 2);
  else
    Read(T_STATUS, buffer, REG_FREQ2, 2);
  
  *f = (static_cast<unsigned int>(buffer[1]) << 8) + buffer[0];
  
  // Convert counts to frequency 
  if (*f != 0)
    *f = 1024 * 200 * (32.768E6 * 4) / (*f) / 1E9;
  
  return 1;
}

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

double DRSBoard::VoltToFreq(double volt)
{
  if (fChipVersion == 3) {
    if (volt <= 1.2001)
      return (volt - 0.6) / 0.2;
    else
      return 0.73/0.28 + sqrt((0.73/0.28)*(0.73/0.28)-2.2/0.14+volt/0.14);
  } else
    return (volt - 0.5) / 0.2;
}

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

double DRSBoard::FreqToVolt(double freq)
{
  if (fChipVersion == 3) {
    if (freq <= 3)
      return 0.6 + 0.2 * freq;
    else
      return 2.2 - 0.73 * freq + 0.14 * freq * freq;
  } else
    return 0.55 + 0.25 * freq;
}

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

int DRSBoard::SetFrequency(double demand)
{
  // Set domino sampling frequency
  double freq, voltage, delta_voltage;
  unsigned short target;
  int i, index, timeout;
  int dominoModeSave = fDominoMode;
  int triggerEnableSave = fTriggerEnable;
  
  SetDominoMode(1);
  EnableTrigger(0);
  EnableAcal(0, 0);
  
  fFrequency = demand;
  
  // Turn automatic adjustment off 
  fCtrlBits &= ~BIT_FREQ_AUTO_ADJ;
  
  // Disable external trigger 
  fCtrlBits &= ~BIT_ENABLE_TRIGGER;
  
  // Set start pulse length for future DRSBoard_domino_start() 
  if (fChipVersion == 2) {
    if (demand < 0.8)
      fCtrlBits |= BIT_LONG_START_PULSE;
    else
      fCtrlBits &= ~BIT_LONG_START_PULSE;
    
    Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  }
  
  // Stop any running domino wave 
  Reinit();
  
  // Estimate DAC setting 
  voltage = FreqToVolt(demand);
  
  SetDAC(fDAC_DSA, voltage);
  SetDAC(fDAC_DSB, voltage);
  
  // Wait until new DAC value has settled 
  Sleep(10);
  
  // Restart domino wave 
  StartDomino();
  
  target = static_cast<unsigned short>(1024 * 200 * (32.768E6 * 4) / demand / 1E9);
  
  // Iterate over both DRS chips 
  for (index = 0; index < 2; index++) {
    
    // Starting voltage 
    voltage = FreqToVolt(demand);
    
    for (i = 0; i < 100; i++) {
      
      // Wait until measurement finished 
      for (timeout = 0; timeout < 1000; timeout++)
        if (IsNewFreq(index))
          break;
      
      freq = 0;
      if (timeout == 1000)
        break;
      
      ReadFrequency(index, &freq);
      
      delta_voltage = FreqToVolt(demand) - FreqToVolt(freq);
      
      if (fDebug) {
        if (fabs(freq - demand) < 0.001)
          printf("CHIP #%d, iter%3d: %1.5lf(%05d) %7.5lf\n", index, i, voltage,
                 static_cast<int>(voltage / 2.5 * 65535 + 0.5), freq);
        else
          printf("CHIP #%d, iter%3d: %1.5lf(%05d) %7.5lf %+5d\n", index, i, voltage,
                 static_cast<int>(voltage / 2.5 * 65535 + 0.5), freq,
                 static_cast<int>(delta_voltage / 2.5 * 65535 + 0.5));
      }
      
      if (fabs(freq - demand) < 0.001)
        break;
      
      voltage += delta_voltage;
      if (voltage > 2.5)
        voltage = 2.5;
      if (voltage < 0)
        voltage = 0;
      
      if (freq == 0)
        break;
      
      if (index == 0)
        SetDAC(fDAC_DSA, voltage);
      else
        SetDAC(fDAC_DSB, voltage);
      
      Sleep(10);
    }
    if (i == 100 || freq == 0 || timeout == 1000) {
      printf("Board %d: could not set frequency of CHIP #%d to %1.3f GHz\n", GetCMCSerialNumber(), index, demand);
      return 0;
    }
  }
  
  SetDominoMode(dominoModeSave);
  EnableTrigger(triggerEnableSave);
  return 1;
}

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

int DRSBoard::RegulateFrequency(double demand)
{
  // Set frequency regulation
  unsigned short target, target_hi, target_lo;
  
  if (demand < 0.42 || demand > 5.2)
    return 0;
  
  fFrequency = demand;
  
  // First iterate DAC value from host 
  if (!SetFrequency(demand))
    return 0;
  
  // Convert frequency in GHz into counts for 200 cycles 
  target = static_cast<unsigned short>(1024 * 200 * (32.768E6 * 4) / demand / 1E9);
  target_hi = target + 6;
  target_lo = target - 6;
  Write(T_CTRL, REG_FREQ_SET_HI, &target_hi, 2);
  Write(T_CTRL, REG_FREQ_SET_LO, &target_lo, 2);
  
  // Turn on regulation 
  fCtrlBits |= BIT_FREQ_AUTO_ADJ;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  // Optional monitoring code ... 
#if 0
  do {
    double freq;
    unsigned short dac, cnt;
    
    ReadFrequency(0, &freq);
    
    if (fBoardVersion == 1)
      Read(T_STATUS, &dac, REG_RDAC3, 2);
    else if (fBoardVersion == 2 || fBoardVersion == 3)
      Read(T_STATUS, &dac, REG_RDAC1, 2);
    
    Read(T_STATUS, &cnt, REG_FREQ1, 2);
    
    if (cnt < 65535)
      printf("%5d %5d %5d %1.5lf\n", dac, target, cnt, freq);
    
    Sleep(500);
  } while (1);
#endif
  
  return 1;
}

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

void DRSBoard::RegisterTest()
{
  // Register test
#define N_REG 8
  
  int i, n, n_err;
  unsigned int buffer[N_REG], ret[N_REG];
  
  printf("**************************************************\n");

  n_err = 0;
  for (n = 0; n < 100; n++) {
    for (i = 0; i < N_REG; i++)
      buffer[i] = (rand() << 16) | rand();
    Write(T_CTRL, 0, buffer, sizeof(buffer));
    
    memset(ret, 0, sizeof(ret));
    i = Read(T_CTRL, ret, 0, sizeof(ret));
    
    while (i != sizeof(ret)) {
      printf(" Read error!\n");
      return;
    }
    
    for (i = 0; i < N_REG; i++)
      if (buffer[i] != ret[i]) {
        n_err++;
        printf(" Error Reg.%d: %08X - %08X\n", i, buffer[i], ret[i]);
      }
      else
        printf(" OK   : %08X\n", buffer[i]);
  }
  
  printf(" Register test: %d errors\n", n_err);
  printf("**************************************************\n");
}

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

int DRSBoard::RAMTest(int flag)
{
#define N_DWORDS (40*1024/4)
  
  // RAM test
  int i, j, n, l=0;
  unsigned int buffer[N_DWORDS], ret[N_DWORDS];
  time_t now;
  
  // Integrity test
  printf("**************************************************\n");
  printf(" Integrity:\n");
  printf("  Buffer size (bytes): %d (%1.1lfk)\n", static_cast<int>(sizeof(ret)), sizeof(ret) / 1024.0);
  if (flag & 1) {
    for (i = 0; i < N_DWORDS; i++)
      buffer[i] = (rand() | rand() << 16) & 0x00FFFFFF;      // Random 24-bit values
    //buffer[i] = i;
    
    // Reset FIFO 
    Reinit();
    
    Write(T_RAM, 0, buffer, sizeof(buffer));
    memset(ret, 0, sizeof(ret));
    
    Read(T_FIFO, ret, 0, sizeof(ret));
    Reinit();
    
    for (i = n = 0; i < N_DWORDS; i++) {
      if (buffer[i] != ret[i])
        n++;
    }
    
    printf("  %d errors\n", n);
  }
  
  // Speed test 
  if (flag & 2) {

    printf(" Speed:\n");

    // Read continously to determine speed 
    time(&now);
    while (now == time(NULL));
    time(&now);
    i = n = 0;
    do {
      memset(ret, 0, sizeof(ret));
      
      for (j = 0; j < 10; j++) {
        Read(T_RAM, ret, 0, sizeof(ret));
        i += sizeof(ret);
      }
      
      if (flag & 1) {
        for (j = 0; j < N_DWORDS; j++)
          if (buffer[j] != ret[j])
            n++;
      }
      
      if (now != time(NULL)) {
        if (flag & 1)
          printf("  %d read/s, %1.2lf MB/s, %d errors\n", static_cast<int>(i / sizeof(ret)), i / 1024.0 / 1024.0,
                 n);
        else
          printf("  %d read/s, %1.2lf MB/s\n", static_cast<int>(i / sizeof(ret)), i / 1024.0 / 1024.0);
        time(&now);
        i = 0;   l++;
      }
    } while (l<5);
  }
  
  printf("**************************************************\n");

  return 0;
}

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

void DRSBoard::SetVoltageOffset(double offset1, double offset2)
{
  if (fChipVersion == 3) {
    SetDAC(fDAC_ROFS_1, 0.95 - offset1);
    SetDAC(fDAC_ROFS_2, 0.95 - offset2);
  } else if (fChipVersion == 2)
    SetDAC(fDAC_COFS, 0.9 - offset1);
}

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

int DRSBoard::SetExternalClockFrequency(double frequencyMHz)
{
  // Set the frequency of the external clock
  fExternalClockFrequency = frequencyMHz;
  return 0;
}

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

double DRSBoard::GetExternalClockFrequency()
{
  // Return the frequency of the external clock
  return fExternalClockFrequency;
}

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

int DRSBoard::TransferWaves(int numberOfChannels)
{
  return TransferWaves(fWaveforms,numberOfChannels);
}

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

int DRSBoard::TransferWaves(unsigned char *p, int numberOfChannels)
{
  return TransferWaves(p, 0, numberOfChannels-1);
}

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

int DRSBoard::TransferWaves(int firstChannel, int lastChannel)
{
  return TransferWaves(fWaveforms, firstChannel, lastChannel);
}

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

int DRSBoard::TransferWaves(unsigned char *p, int firstChannel, int lastChannel)
{

  // Transfer all waveforms at once
  int i, n, offset, n_requested;
  
  if (lastChannel < 0 || lastChannel > kNumberOfChips*kNumberOfChannels) {
    printf("Error: Invalid channel index %d\n", lastChannel);
    return 0;
  }
  
  if (firstChannel < 0 || firstChannel > kNumberOfChips*kNumberOfChannels) {
    printf("Error: Invalid channel index %d\n", firstChannel);
    return 0;
  }
  
  firstChannel = 0;
  lastChannel = kNumberOfChips*kNumberOfChannels -1;
  
  
  n_requested = (lastChannel - firstChannel + 1) * sizeof(short int) * kNumberOfBins;
  offset = firstChannel * sizeof(short int) * kNumberOfBins;
  
  n = Read(T_RAM, p, offset, n_requested);
    
  // Fixme (SCC 28082008): this check is now obsolete!!!! 
  if (n != n_requested) {
    printf("Error: only %d bytes read\n", n);
    return  n;
  }
  
  // Remember which waveforms have been transferred
  for (i = firstChannel; i <= lastChannel; i++)
    fWaveTransferred[i] = true;
    
  //fNumberOfTransferredWaves = numberOfChannels;
  return n;
}

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

int DRSBoard::TransferWaves(unsigned long *p, int numberOfChannels)
{
  return TransferWaves(p, 0, numberOfChannels-1);
}

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

int DRSBoard::TransferWaves(unsigned long *p, int firstChannel, int lastChannel)
{

  // Transfer all waveforms at once
  int i, n, offset, n_requested;
  
  if (lastChannel < 0 || lastChannel > kNumberOfChips*kNumberOfChannels) {
    printf("Error: Invalid channel index %d\n", lastChannel);
    return 0;
  }
  
  if (firstChannel < 0 || firstChannel > kNumberOfChips*kNumberOfChannels) {
    printf("Error: Invalid channel index %d\n", firstChannel);
    return 0;
  }
  
  firstChannel = 0;
  lastChannel = kNumberOfChips*kNumberOfChannels - 1;
  
  
  n_requested = (lastChannel - firstChannel + 1) * sizeof(short int) * kNumberOfBins;
  offset = firstChannel * sizeof(short int) * kNumberOfBins;

  n = Read(T_RAM, p, offset, n_requested);

  if (n != n_requested) {
    printf("Error: only %d bytes read\n", n);
    return  n;
  }
      
  // Remember which waveforms have been transferred
  for (i = firstChannel; i <= lastChannel; i++)
    fWaveTransferred[i] = true;
    
  //fNumberOfTransferredWaves = numberOfChannels;
  return n;
}

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

int DRSBoard::DecodeWave(unsigned int chipIndex, unsigned char channel, unsigned short *waveform)
{
  
  return DecodeWave(fWaveforms , chipIndex, channel, waveform);
}

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

int DRSBoard::DecodeWave(unsigned char *waveforms, unsigned int chipIndex, unsigned char channel, unsigned short *waveform)
{

  // Get waveform
  int i, offset, ind;
  
  // Check valid parameters 
  if (channel > 9 || chipIndex > 1)
    return kWrongChannelOrChip;
  
  // Re-map channel 
  if (fBoardVersion == 1) {
    if (channel < 8)
      channel = 7 - channel;
    else
      channel = 16 - channel;
  } else {
    channel = channel;
  } 
  
  offset = (kNumberOfBins * 4) * channel;

  if (DEBUG)
    printf("offset = %d (0X%X)\n",offset,offset);

  for (i = 0; i < kNumberOfBins; i++) {

    ind = i * 4 + offset;

    if (chipIndex == 0)
      // Lower 12 bit
      waveform[i] = ((waveforms[ind + 1] & 0x0f) << 8) | waveforms[ind];
    else
      // Upper 12 bit
      waveform[i] = (waveforms[ind + 2] << 4) | (waveforms[ind + 1] >> 4);
  }
  
  return kSuccess;
}

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

int DRSBoard::DecodeWave(unsigned long *waveforms, unsigned int chipIndex, unsigned char channel, unsigned short *waveform) 
{

  // Get waveform
  int i, offset, ind;
  
  // Check valid parameters 
  if (channel > 9 || chipIndex > 1)
    return kWrongChannelOrChip;
  
  // Re-map channel 
  if (fBoardVersion == 1) {
    if (channel < 8)
      channel = 7 - channel;
    else
      channel = 16 - channel;
  } else {
    channel = channel;
  } 
  
  offset = kNumberOfBins * channel;

  for (i = 0; i < kNumberOfBins; i++) {
  
    ind = i + offset;

    // Lower 12 bit
    if (chipIndex == 0)
      waveform[i] = (unsigned short)(waveforms[ind] & 0X00000FFF);
    // Upper 12 bit
    else 
      waveform[i] = (unsigned short)(((unsigned long)(waveforms[ind] & 0X00FFF000)) >> 12);
  }
  
  return kSuccess;
}

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

int DRSBoard::GetWave(unsigned int chipIndex, unsigned char channel, short *waveform, bool responseCalib,
                      int triggerCell, bool adjustToClock, float threshold)
{

  int ret;
  
  ret = GetWave(fWaveforms,chipIndex,channel,waveform,responseCalib,triggerCell,adjustToClock,threshold);

  if (kRotateWave) RotateWave((int)GetTriggerCell(chipIndex),waveform);  

  return ret;
}

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

int DRSBoard::GetWave(unsigned int chipIndex, unsigned char channel, float *waveform, bool responseCalib, 
                      int triggerCell, bool adjustToClock, float threshold)
{
  int ret,i;
  short waveS[kNumberOfBins];

  ret = GetWave(fWaveforms,chipIndex,channel,waveS,responseCalib,triggerCell,adjustToClock,threshold);

  if (responseCalib)
    for (i = 0; i < kNumberOfBins; i++)
      waveform[i] = static_cast<float>(waveS[i] * GetPrecision());
  else {
    for (i = 0; i < kNumberOfBins; i++) {
      if (fBoardVersion==4)
        waveform[i] = static_cast<float>(waveS[i] / 4.095); // 12-bit corresponding to 1V
      else
        waveform[i] = static_cast<float>(waveS[i]);
    }
  }

  if (kRotateWave)
    //RotateWave((int)GetTriggerCell(waveform,chipIndex),waveform);  
    RotateWave((int)GetTriggerCell(chipIndex),waveform);  

  return ret;
}

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

void DRSBoard::RotateWave(int triggerCell, float *waveform) 
{
  int i;
  float buffer[kNumberOfBins];

  memcpy((float*)buffer,(float*)waveform,sizeof(buffer));

  for (i=0;i<kNumberOfBins;i++)
    waveform[i] = buffer[(i + triggerCell)%kNumberOfBins];
}

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

void DRSBoard::RotateWave(int triggerCell, short *waveform) 
{
  int i;
  short buffer[kNumberOfBins];

  memcpy((short*)buffer,(short*)waveform,sizeof(buffer));

  for (i=0;i<kNumberOfBins;i++)
    waveform[i] = buffer[(i + triggerCell)%kNumberOfBins];
}

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


int DRSBoard::GetWave(unsigned char *waveforms, unsigned int chipIndex, unsigned char channel, float *waveform, bool responseCalib,
                      int triggerCell, bool adjustToClock, float threshold)
{
  int ret,i;
  short waveS[kNumberOfBins];

  ret = GetWave(waveforms,chipIndex,channel,waveS,responseCalib,triggerCell,adjustToClock,threshold);

  if (responseCalib)
    for (i = 0; i < kNumberOfBins; i++)
      waveform[i] = static_cast<float>(waveS[i] * GetPrecision());
  else {
    for (i = 0; i < kNumberOfBins; i++) {
      if (fBoardVersion == 4)
        waveform[i] = static_cast<float>(waveS[i] / 4.095); // 12-bit corresponding to 1V
      else
        waveform[i] = static_cast<float>(waveS[i]);
    }
  }
  return ret;
}

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


int DRSBoard::GetWave(unsigned char *waveforms, unsigned int chipIndex, unsigned char channel, short *waveform, bool responseCalib,
                      int triggerCell, bool adjustToClock, float threshold)
{
  if (!fWaveTransferred[chipIndex*kNumberOfChannels+channel])
    return kWaveNotAvailable;
  unsigned short adcWaveform[kNumberOfBins];
  int ret = DecodeWave(waveforms, chipIndex, channel, adcWaveform);
  if (ret!=kSuccess)
    return ret;
  
  return CalibrateWaveform(chipIndex, channel, adcWaveform, waveform, responseCalib,
                           triggerCell, adjustToClock, threshold);
}

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

int DRSBoard::GetWave(unsigned long *waveforms, unsigned int chipIndex, unsigned char channel, short *waveform, bool responseCalib, 
                      int triggerCell, bool adjustToClock, float threshold)
{
  if (!fWaveTransferred[chipIndex*kNumberOfChannels+channel])
    return kWaveNotAvailable;
  unsigned short adcWaveform[kNumberOfBins];
  int ret = DecodeWave(waveforms, chipIndex, channel, adcWaveform);
  if (ret!=kSuccess)
    return ret;
  
  return CalibrateWaveform(chipIndex, channel, adcWaveform, waveform, responseCalib,
                           triggerCell, adjustToClock, threshold);
}

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

int DRSBoard::GetADCWave(unsigned int chipIndex, unsigned char channel, unsigned short *waveform)
{
  return DecodeWave(chipIndex, channel, waveform);
}

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

int DRSBoard::GetADCWave(unsigned char *waveforms,unsigned int chipIndex, unsigned char channel, unsigned short *waveform)
{
  return DecodeWave(waveforms, chipIndex, channel, waveform);
}

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

int DRSBoard::GetADCWave(unsigned long *waveforms,unsigned int chipIndex, unsigned char channel, unsigned short *waveform)
{
  return DecodeWave(waveforms, chipIndex, channel, waveform);
}

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

int DRSBoard::CalibrateWaveform(unsigned int chipIndex, unsigned char channel, unsigned short *adcWaveform,
                                short *waveform, bool responseCalib,
                                int triggerCell, bool // adjustToClock 
                                , float threshold)
{
  int j;
  
  // Calibrate waveform
  if (responseCalib) {
    if (!fResponseCalibration->Calibrate(chipIndex, channel % 10, adcWaveform, waveform, triggerCell, threshold))
      return kZeroSuppression; // Return immediately if below threshold
  } else {
    for (j = 0; j < kNumberOfBins; j++) {
      waveform[j] = adcWaveform[j];
    }
  }
  // Fix bad cells for single turn mode
  if (fDominoMode == 0 && triggerCell==-1) {
    waveform[0] = 2 * waveform[1] - waveform[2];
    short m1 = (waveform[kNumberOfBins - 5] + waveform[kNumberOfBins - 6]) / 2;
    short m2 = (waveform[kNumberOfBins - 6] + waveform[kNumberOfBins - 7]) / 2;
    waveform[kNumberOfBins - 4] = m1 - 1 * (m2 - m1);
    waveform[kNumberOfBins - 3] = m1 - 2 * (m2 - m1);
    waveform[kNumberOfBins - 2] = m1 - 3 * (m2 - m1);
    waveform[kNumberOfBins - 1] = m1 - 4 * (m2 - m1);
  }
  return kSuccess;
}

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

int DRSBoard::GetStretchedTime(float *time, float *measurement, int numberOfMeasurements, float period)
{
  int j;
  if (*time >= measurement[numberOfMeasurements - 1]) {
    *time -= measurement[numberOfMeasurements - 1];
    return 1;
  }
  if (*time < measurement[0]) {
    *time = *time - measurement[0] - (numberOfMeasurements - 1) * period / 2;
    return 1;
  }
  for (j = 0; j < numberOfMeasurements - 1; j++) {
    if (*time > measurement[j] && *time <= measurement[j + 1]) {
      *time =
        (period / 2) / (measurement[j + 1] - measurement[j]) * (*time - measurement[j + 1]) -
        (numberOfMeasurements - 2 - j) * period / 2;
      return 1;
    }
  }
  return 0;
}

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

int DRSBoard::GetTriggerCell(unsigned int chipIndex)
{

  return GetTriggerCell(fWaveforms,chipIndex);
}

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

int DRSBoard::GetTriggerCell(unsigned char *waveforms,unsigned int chipIndex)
{
 
  int j, triggerCell;
  bool calib = 0;
  unsigned short baseLevel = 1000;
  unsigned short triggerChannel[1024];
  if (!fWaveTransferred[chipIndex*kNumberOfChannels+8])
    return -1;
  
  GetADCWave(waveforms,chipIndex, 8, triggerChannel);
  //calib = fResponseCalibration->SubtractADCOffset(chipIndex, 8, triggerChannel, triggerChannel,baseLevel); // Changed 24/10/2008, SCC


  triggerCell = -1;
  for (j = 0; j < kNumberOfBins; j++) {
    if (calib) {
      if (triggerChannel[j] <= baseLevel+200
          && triggerChannel[(j + 1) % kNumberOfBins] > baseLevel+200) {
        triggerCell = j;
        break;
      }
    } else {
      if (triggerChannel[j] >= 2000
          && triggerChannel[(j + 1) % kNumberOfBins] < 2000) {
        triggerCell = j;
        break;
      }
    }
  }
  if (triggerCell == -1) {
    return kInvalidTriggerSignal;
  }
  fTriggerCell = triggerCell;
  return triggerCell;
}

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


int DRSBoard::GetTriggerCell(unsigned long *waveforms,unsigned int chipIndex)
{

  int j, triggerCell;
  bool calib = 0;
  unsigned short baseLevel = 1000;
  unsigned short triggerChannel[1024];
  if (!fWaveTransferred[chipIndex*kNumberOfChannels+8])
    return -1;
  
  GetADCWave(waveforms,chipIndex, 8, triggerChannel);
  //calib = fResponseCalibration->SubtractADCOffset(chipIndex, 8, triggerChannel, triggerChannel, baseLevel); // Changed 07/10/2008, SCC
    
  triggerCell = -1;
  for (j = 0; j < kNumberOfBins; j++) {
    if (calib) {

      if (triggerChannel[j] <= baseLevel+200
          && triggerChannel[(j + 1) % kNumberOfBins] > baseLevel+200) {
        triggerCell = j;
        break;
      }
    } else {

      if (triggerChannel[j] >= 2000
          && triggerChannel[(j + 1) % kNumberOfBins] < 2000) {
        triggerCell = j;
        break;
      }
    }
  }
  
  if (triggerCell == -1) {
    return kInvalidTriggerSignal;
  }
  fTriggerCell = triggerCell;
  return triggerCell;

}

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

int DRSBoard::GetTriggerCell(float *waveform)
{
  int j, triggerCell;
    
  triggerCell = -1;

  for (j = 0; j < kNumberOfBins; j++) 
    if ((waveform[(j + 1) % kNumberOfBins]-waveform[j % kNumberOfBins]) > 400.)
      triggerCell = j;
  
    
  if (triggerCell == -1) {
    return kInvalidTriggerSignal;
  }
  fTriggerCell = triggerCell;
  return triggerCell;

}

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

void DRSBoard::TestDAC(int channel)
{
  // Test DAC
  int status;
  
  do {
    status = SetDAC(channel, 0);
    Sleep(1000);
    status = SetDAC(channel, 0.5);
    Sleep(1000);
    status = SetDAC(channel, 1);
    Sleep(1000);
    status = SetDAC(channel, 1.5);
    Sleep(1000);
    status = SetDAC(channel, 2);
    Sleep(1000);
    status = SetDAC(channel, 2.5);
    Sleep(1000);
  } while (status);
}

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

void DRSBoard::MeasureSpeed()
{
  // Measure domino sampling speed
  FILE *f;
  double vdr, vds, freq;
  
  f = fopen("speed.txt", "wt");
  fprintf(f, "\t");
  printf("\t");
  for (vdr = 0.5; vdr <= 2.501; vdr += 0.05) {
    fprintf(f, "%1.2lf\t", vdr);
    printf("%1.2lf\t", vdr);
  }
  fprintf(f, "\n");
  printf("\n");
  
  for (vds = 0.5; vds <= 2.501; vds += 0.05) {
    fprintf(f, "%1.2lf\t", vds);
    printf("%1.2lf\t", vds);
    
    SetDAC(fDAC_DSA, vds);
    StartDomino();
    Sleep(1000);
    ReadFrequency(0, &freq);
    
    fprintf(f, "%1.3lf\t", freq);
    printf("%1.3lf\t", freq);
    
    fprintf(f, "\n");
    printf("\n");
    fflush(f);
  }
}

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

void DRSBoard::InteractSpeed()
{
  int status, i;
  double freq, vds;
  
  do {
    printf("DS: ");
    scanf("%lf", &vds);
    if (vds == 0)
      break;
    
    SetDAC(fDAC_DSA, vds);
    SetDAC(fDAC_DSB, vds);
    
    StartDomino();
    for (i = 0; i < 4; i++) {
      Sleep(1000);
      
      status = ReadFrequency(0, &freq);
      if (!status)
        break;
      printf("%1.6lf GHz\n", freq);
    }
    
    // Turn CMC_LED off 
    SetLED(0);
    
  } while (1);
}

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

void DRSBoard::MonitorFrequency()
{
  // Monitor domino sampling frequency
  int status;
  unsigned int data;
  double freq, dac;
  FILE *f;
  time_t now;
  char str[256];
  
  f = fopen("DRSBoard.log", "w");
  
  do {
    Sleep(1000);
    
    status = ReadFrequency(0, &freq);
    if (!status)
      break;
    
    data = 0;
    if (fBoardVersion == 1)
      Read(T_STATUS, &data, REG_RDAC3, 2);
    else if (fBoardVersion == 2 || fBoardVersion == 3)
      Read(T_STATUS, &data, REG_RDAC1, 2);
    
    dac = data / 65536.0 * 2.5;
    printf("%1.6lf GHz, %1.4lf V\n", freq, dac);
    time(&now);
    strcpy(str, ctime(&now) + 11);
    str[8] = 0;
    
    fprintf(f, "%s %1.6lf GHz, %1.4lf V\n", str, freq, dac);
    fflush(f);
    
  } while (!drs_kbhit());
  
  fclose(f);
}

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

unsigned int DRSBoard::GetCtrlReg()
{
  unsigned int status;
  
  Read(T_CTRL, &status, REG_CTRL, 4);
  return status;
}

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

unsigned int DRSBoard::GetStatusReg()
{
  unsigned int status;
  
  Read(T_STATUS, &status, REG_STATUS, 4);
  return status;
}

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

int DRSBoard::EnableTcal(int flag)
{
  // Enable clock channel
  if (flag)
    fCtrlBits |= BIT_TCAL_EN;
  else
    fCtrlBits &= ~BIT_TCAL_EN;
  
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  
  return 1;
}

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

int DRSBoard::EnableAcal(int mode, double voltage)
{
  double t1, t2;
  
  if (mode == 0) {
    // Turn calibration off 
    SetCalibTiming(0, 0);
    fCtrlBits &= ~BIT_ACAL_EN;
    Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  } else if (mode == 1) {
    // Static calibration 
    SetCalibVoltage(voltage);
    SetCalibTiming(0, 0);
    fCtrlBits |= BIT_ACAL_EN;
    Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  } else if (mode == 2) {
    // First part calibration:
    // stop domino wave after 1.2 revolutions,
    // turn on calibration voltage after 0.1 revolutions 
    
    // Ensure circulating domino wave 
    SetDominoMode(1);
    
    // Set calibration voltage but do not turn it on now 
    SetCalibVoltage(voltage);
    fCtrlBits &= ~BIT_ACAL_EN;
    Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
    
    // Calculate duration of DENABLE signal as 1.2 revolutions 
    t1 = 1 / fFrequency * 1024 * 1.2; // ns
    t1 = static_cast<int>((t1 - 30) / 30 + 1);  // 30 ns offset, 30 ns units, rounded up
    t2 = 1 / fFrequency * 1024 * 0.1; // ns
    t2 = static_cast<int>((t2 - 30) / 30 + 1);  // 30 ns offset, 30 ns units, rounded up
    SetCalibTiming(static_cast<int>(t1), static_cast<int>(t2));
    
  } else if (mode == 3) {
    // Second part calibration:
    // stop domino wave after 1.05 revolutions 
    
    // Ensure circulating domino wave 
    SetDominoMode(1);
    
    // Turn on and let settle calibration voltage 
    SetCalibVoltage(voltage);
    fCtrlBits |= BIT_ACAL_EN;
    Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
    
    // Calculate duration of DENABLE signal as 1.1 revolutions 
    t1 = 1 / fFrequency * 1024 * 1.05;        // ns
    t1 = static_cast<int>((t1 - 30) / 30 + 1);  // 30 ns offset, 30 ns units, rounded up
    SetCalibTiming(static_cast<int>(t1), 0);
  }
  
  return 1;
}

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

int DRSBoard::SetCalibTiming(int t_enable, int t_cal)
{
  unsigned short d;
  
  if (fChipVersion == 2) {
    d = t_cal | (t_enable << 8);
    Write(T_CTRL, REG_CALIB_TIMING, &d, 2);
  }
  
  if (fChipVersion == 3) {
    d = t_cal;
    Write(T_CTRL, REG_CALIB_TIMING, &d, 2);
  }
  
  return 1;
}

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

int DRSBoard::SetCalibVoltage(double value)
{
  // Set Calibration Voltage
  SetDAC(fDAC_ACALIB, value);
  return 1;
}

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

double DRSBoard::GetTemperature()
{
  // Read Out Temperature Sensor
  unsigned char buffer[2];
  unsigned short d;
  double temperature;
  
  Read(T_STATUS, buffer, REG_TEMPERATURE, 2);
  
  d = (static_cast<unsigned int>(buffer[1]) << 8) + buffer[0];
  temperature = ((d >> 3) & 0x0FFF) * 0.0625;
  
  return temperature;
}

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

int DRSBoard::GetTriggerBus()
{
  unsigned short status;
  
  Read(T_STATUS, &status, REG_TRIGGER_BUS, 2);
  return static_cast<int>(status);
}

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

int DRSBoard::FlashEEPROM(unsigned short serial_cmc)
{
  unsigned short dac;
  
  // Read current DAC register
  Read(T_CTRL, &dac, REG_DAC0, 2);
  
  // Put serial in DAC register
  Write(T_CTRL, REG_DAC0, &serial_cmc, 2);
  
  // Execute flash
  fCtrlBits |= BIT_FLASH_TRIG;
  Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
  fCtrlBits &= ~BIT_FLASH_TRIG;
  
  // Wait 6 ms per word
  Sleep(20);
  
  // Write back old DAC registers
  Write(T_CTRL, REG_DAC0, &dac, 2);
  
  // Read back serial number
  ReadSerialNumber();
  
  return 1;
}

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

int DRSBoard::GetTime(unsigned int chipIndex, int frequencyMHz, float *time,int triggerCell)
{
  int i,irot;
  DRSBoard::TimeData * init;
  DRSBoard::TimeData::FrequencyData * freq;
  
  init = GetTimeCalibration(chipIndex);
  
  if (init == NULL) {
    for (i = 0; i < kNumberOfBins; i++)
      time[i] = static_cast<float>(i / fFrequency);
    return 1;
  }
  freq = NULL;
  for (i = 0; i < init->fNumberOfFrequencies; i++) {
    if (init->fFrequency[i]->fFrequency == frequencyMHz) {
      freq = init->fFrequency[i];
      break;
    }
  }
  if (freq == NULL) {
    for (i = 0; i < kNumberOfBins; i++)
      time[i] = static_cast<float>(i / fFrequency);
    return 1;
  }
  for (i = 0; i < kNumberOfBins; i++) {
    irot = i;
    if (triggerCell>-1)
      irot = (triggerCell + i) % kNumberOfBins;
    if (triggerCell + i < kNumberOfBins)
      time[i] = static_cast<float>((freq->fBin[irot] - freq->fBin[triggerCell]) / fFrequency);
    else
      time[i] = static_cast<float>((freq->fBin[irot] - freq->fBin[triggerCell] + freq->fBin[kNumberOfBins - 1] -
                                    2 * freq->fBin[0] + freq->fBin[1]) / fFrequency);
  }
  return 1;
}

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

bool DRSBoard::InitTimeCalibration(unsigned int chipIndex)
{
  return GetTimeCalibration(chipIndex, true) != NULL;
}

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

DRSBoard::TimeData *DRSBoard::GetTimeCalibration(unsigned int chipIndex, bool reinit)
{
  int i, l, index;
  char *cstop;
  char fileName[500];
  char error[240];
  PMXML_NODE node, rootNode, mainNode;
  
  index = fNumberOfTimeData;
  for (i = 0; i < fNumberOfTimeData; i++) {
    if (fTimeData[i]->fChip == static_cast < int >(chipIndex)) {
      if (!reinit)
        return fTimeData[i];
      else {
        index = i;
        break;
      }
    }
  }
  
  fTimeData[index] = new DRSBoard::TimeData();
  DRSBoard::TimeData * init = fTimeData[index];
  
  init->fChip = chipIndex;
  
  for (i = 0; i < init->kMaxNumberOfFrequencies; i++) {
    if (i <= 499 || (i >= 501 && i <= 999) || (i >= 1001 && i <= 1499) || (i >= 1501 && i <= 1999) || (i >= 2001 && i <= 2499) || i >= 2501)
      continue;
    sprintf(fileName, "%s/board%d/TimeCalib_board%d_chip%d_%dMHz.xml", fCalibDirectory, fCMCSerialNumber,
            fCMCSerialNumber, chipIndex, i);
    rootNode = mxml_parse_file(fileName, error, sizeof(error));
    if (rootNode == NULL)
      continue;
    
    init->fFrequency[init->fNumberOfFrequencies] = new DRSBoard::TimeData::FrequencyData();
    init->fFrequency[init->fNumberOfFrequencies]->fFrequency = i;
    
    mainNode = mxml_find_node(rootNode, "/DRSTimeCalibration");
    
    for (l = 0; l < kNumberOfBins; l++) {
      node = mxml_subnode(mainNode, l + 2);
      init->fFrequency[init->fNumberOfFrequencies]->fBin[l] = strtod(mxml_get_value(node), &cstop);
    }
    mxml_free_tree(rootNode);
    init->fNumberOfFrequencies++;
  }
  if (init->fNumberOfFrequencies == 0) {
    printf("Board %d --> Could not find time calibration file\n", GetCMCSerialNumber());
  }
  
  if (index == fNumberOfTimeData)
    fNumberOfTimeData++;
  
  return fTimeData[index];
}

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

void DRSBoard::SetCalibrationDirectory(const char *calibrationDirectoryPath)
{
  strncpy(fCalibDirectory, calibrationDirectoryPath, strlen(calibrationDirectoryPath));
  fCalibDirectory[strlen(calibrationDirectoryPath)] = 0;
};

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

void DRSBoard::GetCalibrationDirectory(char *calibrationDirectoryPath)
{
  strncpy(calibrationDirectoryPath, fCalibDirectory, strlen(fCalibDirectory));
  calibrationDirectoryPath[strlen(fCalibDirectory)] = 0;
};

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

void DRSBoard::LinearRegression(double *x, double *y, int n, double *a, double *b)
{
  int i;
  double sx, sxx, sy, sxy;
  
  sx = sxx = sy = sxy = 0;
  for (i = 0; i < n; i++) {
    sx  += x[i];
    sxx += x[i] * x[i];
    sy  += y[i];
    sxy += x[i] * y[i];
  }
  
  *a = (n * sxy - sx*sy) / (n * sxx - sx * sx);
  *b = (sy - *a * sx) / n;
}

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

void ResponseCalibration::SetCalibrationParameters(int numberOfPointsLowVolt, int numberOfPoints,
                                                   int numberOfMode2Bins, int numberOfSamples,
                                                   int numberOfGridPoints, int numberOfXConstPoints,
                                                   int numberOfXConstGridPoints, double triggerFrequency,
                                                   int showStatistics)
{
  DeleteFields();
  InitFields(numberOfPointsLowVolt, numberOfPoints, numberOfMode2Bins, numberOfSamples, numberOfGridPoints,
             numberOfXConstPoints, numberOfXConstGridPoints, triggerFrequency, showStatistics);
}

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

void ResponseCalibration::ResetCalibration()
{
  int i;
  for (i = 0; i < kNumberOfChips; i++)
    fCalibrationData[i]->fRead = false;
  fCurrentPoint = 0;
  fCurrentLowVoltPoint = 0;
  fCurrentSample = 0;
  fCurrentFitChannel = 0;
  fCurrentFitBin = 0;
  fRecorded = false;
  fFitted = false;
  fOffset = false;
};

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

bool ResponseCalibration::WriteCalibration(unsigned int chipIndex)
{
  if (!fOffset)
    return false;
  if (fBoard->GetChipVersion() == 3)
    return WriteCalibrationV4(chipIndex);
  else
    return WriteCalibrationV3(chipIndex);
}

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

bool ResponseCalibration::WriteCalibrationV3(unsigned int chipIndex)
{
  if (!fOffset)
    return false;
  
  int ii, j, k;
  char str[1000];
  char strt[1000];
  short tempShort;
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  
  // Open File
  fBoard->GetCalibrationDirectory(strt);
  sprintf(str, "%s/board%d", strt, fBoard->GetCMCSerialNumber());
  if (MakeDir(str) == -1) {
    printf("Error: Cannot create directory \"%s\"\n", str);
    return false;
  }
  sprintf(str, "%s/board%d/ResponseCalib_board%d_chip%d_%dMHz.bin", strt, fBoard->GetCMCSerialNumber(),
          fBoard->GetCMCSerialNumber(), chipIndex, static_cast<int>(fBoard->GetFrequency() * 1000));
  fCalibFile = fopen(str, "wb");
  if (fCalibFile == NULL) {
    printf("Error: Cannot write to file \"%s\"\n", str);
    return false;
  }
  // Write File
  fwrite(&data->fNumberOfGridPoints, 1, 1, fCalibFile);
  tempShort = static_cast<short>(data->fStartTemperature) * 10;
  fwrite(&tempShort, 2, 1, fCalibFile);
  tempShort = static_cast<short>(data->fEndTemperature) * 10;
  fwrite(&tempShort, 2, 1, fCalibFile);
  fwrite(&data->fMin, 4, 1, fCalibFile);
  fwrite(&data->fMax, 4, 1, fCalibFile);
  fwrite(&data->fNumberOfLimitGroups, 1, 1, fCalibFile);
  
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    chn = data->fChannel[ii];
    for (j = 0; j < kNumberOfBins; j++) {
      fwrite(&chn->fLimitGroup[j], 1, 1, fCalibFile);
      fwrite(&chn->fLookUpOffset[j], 2, 1, fCalibFile);
      fwrite(&chn->fNumberOfLookUpPoints[j], 1, 1, fCalibFile);
      for (k = 0; k < chn->fNumberOfLookUpPoints[j]; k++) {
        fwrite(&chn->fLookUp[j][k], 1, 1, fCalibFile);
      }
      for (k = 0; k < data->fNumberOfGridPoints; k++) {
        fwrite(&chn->fData[j][k], 2, 1, fCalibFile);
      }
      fwrite(&chn->fOffsetADC[j], 2, 1, fCalibFile);
      fwrite(&chn->fOffset[j], 2, 1, fCalibFile);
    }
  }
  fclose(fCalibFile);
  
  printf("Calibration successfully written to\n\"%s\"\n", str);
  return true;
}

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

bool ResponseCalibration::WriteCalibrationV4(unsigned int chipIndex)
{
  if (!fOffset)
    return false;
  
  int ii, j;
  char str[1000];
  char strt[1000];
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  
  // Open File
  fBoard->GetCalibrationDirectory(strt);
  sprintf(str, "%s/board%d", strt, fBoard->GetCMCSerialNumber());
  if (MakeDir(str) == -1) {
    printf("Error: Cannot create directory \"%s\"\n", str);
    return false;
  }
  sprintf(str, "%s/board%d/ResponseCalib_board%d_chip%d_%dMHz.bin", strt, fBoard->GetCMCSerialNumber(),
          fBoard->GetCMCSerialNumber(), chipIndex, static_cast<int>(fBoard->GetFrequency() * 1000));
  fCalibFile = fopen(str, "wb");
  if (fCalibFile == NULL) {
    printf("Error: Cannot write to file \"%s\"\n", str);
    return false;
  }
  
  // Write File
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    chn = data->fChannel[ii];
    for (j = 0; j < kNumberOfBins; j++) {
      fwrite(&chn->fOffset[j], 2, 1, fCalibFile);
      fwrite(&chn->fGain[j], 2, 1, fCalibFile);
    }
  }
  fclose(fCalibFile);
  
  printf("Calibration successfully written to\n\"%s\"\n", str);
  return true;
}

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

void ResponseCalibration::CalibrationTrigger(int mode, double voltage)
{
  fBoard->Reinit();
  fBoard->EnableAcal(mode, voltage);
  fBoard->StartDomino();
  fBoard->SoftTrigger();
  while (fBoard->IsBusy()) {
  }
}

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

void ResponseCalibration::CalibrationStart(double voltage)
{
  fBoard->SetDominoMode(1);
  fBoard->EnableAcal(0, voltage);
  fBoard->StartDomino();
  fBoard->IsBusy();
  fBoard->IsBusy();
  fBoard->IsBusy();
}

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

bool ResponseCalibration::RecordCalibrationPoints(int chipNumber)
{
  if (!fInitialized)
    return true;
  if (fBoard->GetChipVersion() == 3)
    return RecordCalibrationPointsV4(chipNumber);
  else
    return RecordCalibrationPointsV3(chipNumber);
}

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

bool ResponseCalibration::RecordCalibrationPointsV3(int chipNumber)
{

  int j, k, ii;
  int notdone, nsample;
  double voltage;
  float mean;
  const double minVolt = 0.006;
  const double xpos[50] =
    { 0.010, 0.027, 0.052, 0.074, 0.096, 0.117, 0.136, 0.155, 0.173, 0.191, 0.208, 0.226, 0.243, 0.260,
      0.277, 0.294, 0.310,
      0.325, 0.342, 0.358, 0.374, 0.390, 0.406, 0.422, 0.439, 0.457, 0.477, 0.497, 0.520, 0.546, 0.577, 0.611,
      0.656, 0.710,
      0.772, 0.842, 0.916,
      0.995, 1.075, 1.157, 1.240, 1.323, 1.407, 1.490, 1.575, 1.659, 1.744, 1.829, 1.914, 2.000
    };
  
  // Initialisations
  if (fCurrentLowVoltPoint == 0) {
    fBoard->SetDAC(fBoard->fDAC_CLKOFS, 0);
    // Record Temperature
    fCalibrationData[chipNumber]->fStartTemperature = static_cast<float>(fBoard->GetTemperature());
  }
  // Record current Voltage
  if (fCurrentLowVoltPoint < fNumberOfPointsLowVolt)
    voltage = (xpos[0] - minVolt) * fCurrentLowVoltPoint / static_cast<double>(fNumberOfPointsLowVolt) + minVolt;
  else
    voltage = xpos[fCurrentPoint];
  fBoard->SetCalibVoltage(voltage);
  fResponseY[fCurrentPoint + fCurrentLowVoltPoint] = static_cast<float>(voltage) * 1000;
  
  // Loop Over Number Of Samples For Statistics
  for (j = 0; j < fNumberOfSamples; j++) {
    // Read Out Second Part of the Waveform
    CalibrationTrigger(3, voltage);
    fBoard->TransferWaves();
    for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
      fBoard->GetADCWave(chipNumber, ii, fWaveFormMode3[ii][j]);
    }
    // Read Out First Part of the Waveform
    CalibrationStart(voltage);
    CalibrationTrigger(2, voltage);
    fBoard->TransferWaves();
    for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
      fBoard->GetADCWave(chipNumber, ii, fWaveFormMode2[ii][j]);
    }
    CalibrationStart(voltage);
  }
  // Average Sample Points
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    for (k = 0; k < kNumberOfBins; k++) {
      fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] = 0;
      for (j = 0; j < fNumberOfSamples; j++) {
        fSampleUsed[j] = 1;
        if (k < fNumberOfMode2Bins)
          fSamples[j] = fWaveFormMode2[ii][j][k];
        else
          fSamples[j] = fWaveFormMode3[ii][j][k];
        fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] += fSamples[j];
      }
      mean = fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] / fNumberOfSamples;
      notdone = 1;
      nsample = fNumberOfSamples;
      while (notdone) {
        notdone = 0;
        for (j = 0; j < fNumberOfSamples; j++) {
          if (fSampleUsed[j] && abs(static_cast<int>(fSamples[j] - mean)) > 3) {
            notdone = 1;
            fSampleUsed[j] = 0;
            nsample--;
            fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] -= fSamples[j];
            mean = fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] / nsample;
          }
        }
      }
      fResponseX[ii][k][fCurrentPoint + fCurrentLowVoltPoint] = mean;
    }
  }
  if (fCurrentLowVoltPoint < fNumberOfPointsLowVolt)
    fCurrentLowVoltPoint++;
  else
    fCurrentPoint++;
  
  if (fCurrentPoint == fNumberOfPoints) {
    fCalibrationData[chipNumber]->fEndTemperature = static_cast<float>(fBoard->GetTemperature());
    fRecorded = true;
    fFitted = false;
    fOffset = false;
    fCalibrationData[chipNumber]->fRead = false;
    fCalibrationData[chipNumber]->fHasOffsetCalibration = false;
    fBoard->SetCalibVoltage(0.0);
    fBoard->EnableAcal(1, 0.0);
    fBoard->SetDAC(fBoard->fDAC_CLKOFS, 0.0);
    return true;
  }
  
  return false;
}

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

bool ResponseCalibration::RecordCalibrationPointsV4(int chipNumber)
{
  int i, j, k, n;
  double voltage, s, s2, average, sigma;
  
  if (fCurrentPoint == 0) {
    fBoard->SetDominoMode(1);
    fBoard->EnableAcal(1,0);
    fBoard->SoftTrigger();
    while (fBoard->IsBusy());
    fBoard->StartDomino();
    fCalibrationData[chipNumber]->fStartTemperature = static_cast<float>(fBoard->GetTemperature());
  }
  voltage = 1.0 * fCurrentPoint / (static_cast<double>(fNumberOfPoints) - 1) + 0.1;
  fBoard->SetCalibVoltage(voltage);
  Sleep(10);
  fBoard->SetCalibVoltage(voltage);
  Sleep(10);
  
  // One dummy cycle for unknown reasons
  fBoard->SoftTrigger();
  while (fBoard->IsBusy());
  fBoard->StartDomino();
  //Sleep(50);
  fBoard->TransferWaves();
  
  // Loop over number of samples for statistics
  for (i = 0; i < fNumberOfSamples; i++) {
#ifdef DEBUG_CALIB
    printf("%d   \r", fNumberOfSamples-i);
#endif
    fBoard->SoftTrigger();
    while (fBoard->IsBusy());
    fBoard->StartDomino();
    //Sleep(50);
    fBoard->TransferWaves();
    for (j = 0; j < kNumberOfCalibChannels; j++) {
      fBoard->GetADCWave(chipNumber, j, fWaveFormMode3[j][i]);
    }
  }
  
  for (i = 0; i < kNumberOfCalibChannels; i++) {
    for (k = 0; k < kNumberOfBins; k++) {
      s = s2 = 0;
      
      for (j=0 ; j<fNumberOfSamples ; j++) {
        s += fWaveFormMode3[i][j][k];
        s2 += fWaveFormMode3[i][j][k] * fWaveFormMode3[i][j][k];
      }
      n = fNumberOfSamples;
      average = s / n;
      sigma = sqrt( (n*s2 - s*s) / (n*(n - 1)) );
      
      fResponseX[i][k][fCurrentPoint] = static_cast<float>(average);
    }
  }
  
#ifdef DEBUG_CALIB
  for (j = 0; j < fNumberOfSamples; j++)
    printf("%d ", fWaveFormMode3[1][j][10]);
  
  s = s2 = 0;
  for (j = 0; j < fNumberOfSamples; j++) {
    s += fWaveFormMode3[i][j][k];
    s2 += fWaveFormMode3[i][j][k] * fWaveFormMode3[i][j][k];
  }
  n = fNumberOfSamples;
  average = s / n;
  sigma = sqrt( (n*s2 - s*s) / (n*(n - 1)) );
  
  printf("\n");
  printf("%1.2lf V: %6.1lf (%1.4lf)\n", voltage,
         fResponseX[1][10][fCurrentPoint],
         fResponseX[1][10][fCurrentPoint]/4096.0);
#endif
  
  fCurrentPoint++;
  if (fCurrentPoint == fNumberOfPoints) {
    fCalibrationData[chipNumber]->fEndTemperature = static_cast<float>(fBoard->GetTemperature());
    fRecorded = true;
    return true;
  }
  
  return false;
}

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

bool ResponseCalibration::FitCalibrationPoints(int chipNumber)
{
  if (!fRecorded || fFitted)
    return true;
  if (fBoard->GetChipVersion() == 3)
    return FitCalibrationPointsV4(chipNumber);
  else
    return FitCalibrationPointsV3(chipNumber);
}

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

bool ResponseCalibration::FitCalibrationPointsV3(int chipNumber)
{
  int i, j, k;
  float x1, x2, y1, y2;
  float uu;
  float yc, yr;
  float xminExt, xrangeExt;
  float xmin, xrange;
  float average, averageError, averageExt, averageErrorExt;
  unsigned short i0, i1;
  
  CalibrationData *data = fCalibrationData[chipNumber];
  CalibrationData::CalibrationDataChannel * chn = data->fChannel[fCurrentFitChannel];
  
  data->DeletePreCalculatedBSpline();
  
  if (fCurrentFitBin == 0 && fCurrentFitChannel == 0) {
    data->fNumberOfLimitGroups = 0;
    data->fMin = 100000;
    data->fMax = -100000;
    for (i = 0; i < kNumberOfCalibChannels; i++) {
      for (j = 0; j < kNumberOfBins; j++) {
        if (data->fMin > fResponseX[i][j][fNumberOfPointsLowVolt + fNumberOfPoints - 1])
          data->fMin = fResponseX[i][j][fNumberOfPointsLowVolt + fNumberOfPoints - 1];
        if (data->fMax < fResponseX[i][j][fNumberOfPointsLowVolt])
          data->fMax = fResponseX[i][j][fNumberOfPointsLowVolt];
      }
    }
  }

  // Low Volt
  i0 = static_cast<unsigned short>(fResponseX[fCurrentFitChannel][fCurrentFitBin][0]);
  i1 = static_cast<unsigned short>(fResponseX[fCurrentFitChannel][fCurrentFitBin][fNumberOfPointsLowVolt]) + 1;
  chn->fLookUpOffset[fCurrentFitBin] = i0;
  delete chn->fLookUp[fCurrentFitBin];
  if (i0 - i1 + 1 < 2) {
    chn->fNumberOfLookUpPoints[fCurrentFitBin] = 2;
    chn->fLookUp[fCurrentFitBin] = new unsigned char[2];
    chn->fLookUp[fCurrentFitBin][0] = 0;
    chn->fLookUp[fCurrentFitBin][1] = 0;
  } else {
    chn->fNumberOfLookUpPoints[fCurrentFitBin] = i0 - i1 + 1;
    chn->fLookUp[fCurrentFitBin] = new unsigned char[i0 - i1 + 1];
    for (i = 0; i < i0 - i1 + 1; i++) {
      for (j = 0; j < fNumberOfPointsLowVolt; j++) {
        if (i0 - i >= fResponseX[fCurrentFitChannel][fCurrentFitBin][j + 1]) {
          x1 = fResponseX[fCurrentFitChannel][fCurrentFitBin][j];
          x2 = fResponseX[fCurrentFitChannel][fCurrentFitBin][j + 1];
          y1 = fResponseY[j];
          y2 = fResponseY[j + 1];
          chn->fLookUp[fCurrentFitBin][i] = static_cast<unsigned char>(((y2 - y1) * (i0 - i - x1) / (x2 - x1) + y1)/fPrecision);
          break;
        }
      }
    }
  }
  

  // Copy Points
  for (i = 0; i < fNumberOfPoints; i++) {
    fPntX[0][i] = fResponseX[fCurrentFitChannel][fCurrentFitBin][fNumberOfPointsLowVolt + i];
    fPntY[0][i] = fResponseY[fNumberOfPointsLowVolt + i];
  }
  // Fit BSpline
  for (i = 0; i < fNumberOfPoints; i++) {
    fUValues[0][i] = static_cast<float>(1 - i / (fNumberOfPoints - 1.));
  }
  if (!Approx(fPntX[0], fUValues[0], fNumberOfPoints, fNumberOfGridPoints, fResX[fCurrentFitBin]))
    return true;
  if (!Approx(fPntY[0], fUValues[0], fNumberOfPoints, fNumberOfGridPoints, fRes[fCurrentFitBin]))
    return true;
  

  // X constant fit
  for (k = 0; k < fNumberOfXConstPoints - 2; k++) {
    fPntX[1][k + 1] =
      GetValue(fResX[fCurrentFitBin], static_cast<float>(1 - k / static_cast<float>(fNumberOfXConstPoints - 3)),
               fNumberOfGridPoints);
    fPntY[1][k + 1] =
      GetValue(fRes[fCurrentFitBin], static_cast<float>(1 - k / static_cast<float>(fNumberOfXConstPoints - 3)),
               fNumberOfGridPoints);
  }
  xmin = fPntX[1][fNumberOfXConstPoints - 2];
  xrange = fPntX[1][1] - xmin;
  
  for (i = 0; i < fNumberOfXConstPoints - 2; i++) {
    fUValues[1][i + 1] = (fPntX[1][i + 1] - xmin) / xrange;
  }
  
  if (!Approx
      (&fPntY[1][1], &fUValues[1][1], fNumberOfXConstPoints - 2, fNumberOfXConstGridPoints,
       chn->fTempData))
    return true;

  
  // Error statistics
  if (fShowStatistics) {
    for (i = 0; i < fNumberOfPoints; i++) {
      uu = (fPntX[0][i] - xmin) / xrange;
      yc = GetValue(chn->fTempData, uu, fNumberOfXConstGridPoints);
      yr = fPntY[0][i];
      fStatisticsApprox[i][fCurrentFitBin + fCurrentFitChannel * kNumberOfBins] = yc - yr;
    }
  }
  // Add min and max point
  chn->fLimitGroup[fCurrentFitBin] = 0;
  while (xmin - kBSplineXMinOffset > data->fMin + kBSplineXMinOffset * chn->fLimitGroup[fCurrentFitBin]) {
    chn->fLimitGroup[fCurrentFitBin]++;
  }
  if (data->fNumberOfLimitGroups <= chn->fLimitGroup[fCurrentFitBin])
    data->fNumberOfLimitGroups = chn->fLimitGroup[fCurrentFitBin] + 1;
  xminExt = data->fMin + kBSplineXMinOffset * chn->fLimitGroup[fCurrentFitBin];
  xrangeExt = data->fMax - xminExt;
  
  fPntX[1][0] = data->fMax;
  uu = (fPntX[1][0] - xmin) / xrange;
  fPntY[1][0] = GetValue(chn->fTempData, uu, fNumberOfXConstGridPoints);
  
  fPntX[1][fNumberOfXConstPoints - 1] = xminExt;
  uu = (fPntX[1][fNumberOfXConstPoints - 1] - xmin) / xrange;
  fPntY[1][fNumberOfXConstPoints - 1] = GetValue(chn->fTempData, uu, fNumberOfXConstGridPoints);
  
  for (i = 0; i < fNumberOfXConstPoints; i++) {
    fUValues[1][i] = (fPntX[1][i] - xminExt) / xrangeExt;
  }
  
  if (!Approx
      (fPntY[1], fUValues[1], fNumberOfXConstPoints, fNumberOfXConstGridPoints, chn->fTempData))
    return true;
  
  // Error statistics
  if (fShowStatistics) {
    for (i = 0; i < fNumberOfPoints; i++) {
      uu = (fPntX[0][i] - xminExt) / xrangeExt;
      yc = GetValue(chn->fTempData, uu, fNumberOfXConstGridPoints);
      yr = fPntY[0][i];
      fStatisticsApproxExt[i][fCurrentFitBin + fCurrentFitChannel * kNumberOfBins] = yc - yr;
    }
  }
  for (i = 0; i < fNumberOfXConstGridPoints; i++) {
    chn->fData[fCurrentFitBin][i] = static_cast<short>(chn->fTempData[i] / fPrecision);
  }
  

  // Write end of file
  fCurrentFitBin++;
  if (fCurrentFitBin == kNumberOfBins) {
    fCurrentFitChannel++;
    fCurrentFitBin = 0;
  }
  if (fCurrentFitChannel == kNumberOfCalibChannels) {
    if (fShowStatistics) {
      for (i = 0; i < fNumberOfPoints; i++) {
        average = 0;
        averageError = 0;
        averageExt = 0;
        averageErrorExt = 0;
        for (j = 0; j < kNumberOfCalibChannels * kNumberOfBins; j++) {
          average += fStatisticsApprox[i][j];
          averageError += fStatisticsApprox[i][j] * fStatisticsApprox[i][j];
          averageExt += fStatisticsApproxExt[i][j];
          averageErrorExt += fStatisticsApproxExt[i][j] * fStatisticsApproxExt[i][j];
        }
        average /= kNumberOfCalibChannels * kNumberOfBins;
        averageError =
          sqrt((averageError -
                average * average / kNumberOfCalibChannels * kNumberOfBins) / (kNumberOfCalibChannels *
                                                                               kNumberOfBins - 1));
        averageExt /= kNumberOfCalibChannels * kNumberOfBins;
        averageErrorExt =
          sqrt((averageErrorExt -
                averageExt * averageExt / kNumberOfCalibChannels * kNumberOfBins) /
               (kNumberOfCalibChannels * kNumberOfBins - 1));
        printf("Error at %3.1f V : % 2.3f +- % 2.3f ; % 2.3f +- % 2.3f\n", fPntY[0][i], average,
               averageError, averageExt, averageErrorExt);
      }
    }
    fFitted = true;
    fOffset = false;
    fCalibrationData[chipNumber]->fRead = true;
    fCalibrationData[chipNumber]->fHasOffsetCalibration = false;
    data->PreCalculateBSpline();
    return true;
  }
  return false;
}

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

bool ResponseCalibration::FitCalibrationPointsV4(int chipNumber)
{
  if (!fRecorded || fFitted)
    return true;
  int i;
  double par[2];
  
  CalibrationData *data = fCalibrationData[chipNumber];
  CalibrationData::CalibrationDataChannel * chn = data->fChannel[fCurrentFitChannel];
  
  if (fCurrentFitBin == 0 && fCurrentFitChannel == 0) {
    for (i = 0; i < fNumberOfPoints; i++)
      fWWFit[i] = 1;
  }
  
  for (i = 0; i < fNumberOfPoints; i++) {
    fXXFit[i] = 1.0 * i / (static_cast<double>(fNumberOfPoints) - 1) + 0.1;
    fYYFit[i] = fResponseX[fCurrentFitChannel][fCurrentFitBin][i];
    if (fCurrentFitBin == 10 && fCurrentFitChannel == 1) {
      fXXSave[i] = fXXFit[i];
      fYYSave[i] = fYYFit[i];
    }
  }
  
  DRSBoard::LinearRegression(fXXFit, fYYFit, fNumberOfPoints, &par[1], &par[0]);
  chn->fOffset[fCurrentFitBin] = static_cast<short>(par[0] + 0.5);
  chn->fGain[fCurrentFitBin] = static_cast<short>(par[1] + 0.5);
  
  if (fCurrentFitChannel == 1 && fCurrentFitBin == 10) {
#ifdef DEBUG_CALIB
    printf("gain:%d, offset:%d\n", chn->fGain[10], chn->fOffset[10]);
#endif
    for (i = 0; i < fNumberOfPoints; i++) {
      fXXSave[i] = fXXFit[i];
      fYYSave[i] = (fYYFit[i] - chn->fOffset[10]) / chn->fGain[10] - fXXFit[i];
    }
  }
  
  fCurrentFitBin++;
  if (fCurrentFitBin == kNumberOfBins) {
    fCurrentFitChannel++;
    fCurrentFitBin = 0;
  }
  if (fCurrentFitChannel == kNumberOfCalibChannels) {
    fFitted = true;
    fOffset = true;
    fCalibrationData[chipNumber]->fRead = true;
    fCalibrationData[chipNumber]->fHasOffsetCalibration = false;
    return true;
  }
  
  return false;
}

unsigned int millitime()
{
  struct timeval tv;
  
  gettimeofday(&tv, NULL);
  
  return tv.tv_sec * 1000 + tv.tv_usec / 1000;
}

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

bool ResponseCalibration::OffsetCalibration(int chipNumber)
{
  if (!fFitted || fOffset)
    return true;
  int k, ii, j;
  int t1, t2;
  float mean,error;
  CalibrationData *data = fCalibrationData[chipNumber];
  CalibrationData::CalibrationDataChannel * chn;
  
  if (fCurrentSample == 0) {
    data->fHasOffsetCalibration = false;
    fBoard->SetCalibVoltage(0.0);
    fBoard->EnableAcal(1, 0.0);
  }
  // Loop Over Number Of Samples For Statistics
  t1 = millitime();
  fBoard->SoftTrigger();
  while (fBoard->IsBusy()) {
  }
  fBoard->TransferWaves();
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    fBoard->GetADCWave(chipNumber, ii, fWaveFormOffsetADC[ii][fCurrentSample]);
    fBoard->CalibrateWaveform(chipNumber, ii, fWaveFormOffsetADC[ii][fCurrentSample], fWaveFormOffset[ii][fCurrentSample],
                              true, false, false, 0);
  }
  fBoard->StartDomino();
  fBoard->IsBusy();
  fBoard->IsBusy();
  fBoard->IsBusy();
  t2 = millitime();
  while (t2 - t1 < (1000 / fTriggerFrequency)) {
    t2 = millitime();
  }
  fCurrentSample++;
  
  if (fCurrentSample == fNumberOfSamples) {
    // Average Sample Points
    float *sample = new float[fNumberOfSamples];
    for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
      chn = data->fChannel[ii];
      for (k = 0; k < kNumberOfBins; k++) {
        for (j = 0; j < fNumberOfSamples; j++)
          sample[j] = static_cast<float>(fWaveFormOffset[ii][j][k]);
        Average(1, sample, fNumberOfSamples, mean, error, 2);
        chn->fOffset[k] = static_cast<short>(mean);
        for (j = 0; j < fNumberOfSamples; j++)
          sample[j] = fWaveFormOffsetADC[ii][j][k];
        Average(1, sample, fNumberOfSamples, mean, error, 2);
        chn->fOffsetADC[k] = static_cast<unsigned short>(mean);
      }
    }
    fOffset = true;
    fCalibrationData[chipNumber]->fHasOffsetCalibration = true;
    delete sample;
    return true;
  }
  
  return false;
}

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

void ResponseCalibration::InitFields(int numberOfPointsLowVolt, int numberOfPoints, int numberOfMode2Bins,
                                     int numberOfSamples, int numberOfGridPoints, int numberOfXConstPoints,
                                     int numberOfXConstGridPoints, double triggerFrequency,
                                     int showStatistics)
{
  int ii, j, i;
  fInitialized = true;
  fNumberOfPointsLowVolt = numberOfPointsLowVolt;
  fNumberOfPoints = numberOfPoints;
  fNumberOfMode2Bins = numberOfMode2Bins;
  fNumberOfSamples = numberOfSamples;
  fNumberOfGridPoints = numberOfGridPoints;
  fNumberOfXConstPoints = numberOfXConstPoints;
  fNumberOfXConstGridPoints = numberOfXConstGridPoints;
  fTriggerFrequency = triggerFrequency;
  fShowStatistics = showStatistics;
  fCurrentPoint = 0;
  fCurrentSample = 0;
  fCurrentFitChannel = 0;
  fCurrentFitBin = 0;
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    for (j = 0; j < kNumberOfBins; j++) {
      fResponseX[ii][j] = new float[fNumberOfPoints + fNumberOfPointsLowVolt];
    }
  }
  fResponseY = new float[fNumberOfPoints + fNumberOfPointsLowVolt];
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    fWaveFormMode3[ii] = new unsigned short *[fNumberOfSamples];
    fWaveFormMode2[ii] = new unsigned short *[fNumberOfSamples];
    fWaveFormOffset[ii] = new short *[fNumberOfSamples];
    fWaveFormOffsetADC[ii] = new unsigned short *[fNumberOfSamples];
    for (i = 0; i < fNumberOfSamples; i++) {
      fWaveFormMode3[ii][i] = new unsigned short[kNumberOfBins];
      fWaveFormMode2[ii][i] = new unsigned short[kNumberOfBins];
      fWaveFormOffset[ii][i] = new short[kNumberOfBins];
      fWaveFormOffsetADC[ii][i] = new unsigned short[kNumberOfBins];
    }
  }
  fSamples = new unsigned short[fNumberOfSamples];
  fSampleUsed = new int[fNumberOfSamples];
  
  for (j = 0; j < kNumberOfBins; j++) {
    fRes[j] = new float[fNumberOfGridPoints];
    fResX[j] = new float[fNumberOfGridPoints];
  }
  for (i = 0; i < 2; i++) {
    fPntX[i] = new float[fNumberOfPoints * (1 - i) + fNumberOfXConstPoints * i];
    fPntY[i] = new float[fNumberOfPoints * (1 - i) + fNumberOfXConstPoints * i];
    fUValues[i] = new float[fNumberOfPoints * (1 - i) + fNumberOfXConstPoints * i];
  }
  fXXFit = new double[fNumberOfPoints];
  fYYFit = new double[fNumberOfPoints];
  fWWFit = new double[fNumberOfPoints];
  fYYFitRes = new double[fNumberOfPoints];
  fYYSave = new double[fNumberOfPoints];
  fXXSave = new double[fNumberOfPoints];
  
  fStatisticsApprox = new float *[fNumberOfPoints];
  fStatisticsApproxExt = new float *[fNumberOfPoints];
  for (i = 0; i < fNumberOfPoints; i++) {
    fStatisticsApprox[i] = new float[kNumberOfCalibChannels * kNumberOfBins];
    fStatisticsApproxExt[i] = new float[kNumberOfCalibChannels * kNumberOfBins];
  }
  for (i = 0; i < kNumberOfChips; i++) {
    fCalibrationData[i] = new CalibrationData(numberOfXConstGridPoints);
  }
}

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

void ResponseCalibration::DeleteFields()
{
  if (!fInitialized)
    return;
  fInitialized = false;
  int ii, j, i;
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    for (j = 0; j < kNumberOfBins; j++) {
      delete fResponseX[ii][j];
    }
  }
  delete fResponseY;
  for (ii = 0; ii < kNumberOfCalibChannels; ii++) {
    for (i = 0; i < fNumberOfSamples; i++) {
      if (fWaveFormMode3[ii] != NULL)
        delete fWaveFormMode3[ii][i];
      if (fWaveFormMode2[ii] != NULL)
        delete fWaveFormMode2[ii][i];
      if (fWaveFormOffset[ii] != NULL)
        delete fWaveFormOffset[ii][i];
      if (fWaveFormOffsetADC[ii] != NULL)
        delete fWaveFormOffsetADC[ii][i];
    }
    delete fWaveFormMode3[ii];
    delete fWaveFormMode2[ii];
    delete fWaveFormOffset[ii];
    delete fWaveFormOffsetADC[ii];
  }
  delete fSamples;
  delete fSampleUsed;
  
  for (j = 0; j < kNumberOfBins; j++) {
    delete fRes[j];
    delete fResX[j];
  }
  for (i = 0; i < 2; i++) {
    delete fPntX[i];
    delete fPntY[i];
    delete fUValues[i];
  }
  delete fXXFit;
  delete fYYFit;
  delete fWWFit;
  delete fYYFitRes;
  delete fYYSave;
  delete fXXSave;
  
  for (i = 0; i < fNumberOfPoints; i++) {
    delete fStatisticsApprox[i];
    delete fStatisticsApproxExt[i];
  }
  delete fStatisticsApprox;
  delete fStatisticsApproxExt;
  for (i = 0; i < kNumberOfChips; i++)
    delete fCalibrationData[i];
}

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

double ResponseCalibration::GetTemperature(unsigned int chipIndex)
{
  if (fCalibrationData[chipIndex]==NULL)
    return 0;
  if (!fCalibrationData[chipIndex]->fRead)
    return 0;
  return (fCalibrationData[chipIndex]->fStartTemperature + fCalibrationData[chipIndex]->fEndTemperature) / 2;
}

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

bool ResponseCalibration::Calibrate(unsigned int chipIndex, unsigned int channel, unsigned short *adcWaveform,
                                    short *uWaveform, int triggerCell, float threshold)
{
  int i;
  int hasOffset;
  bool aboveThreshold;
  float wave, v;
  int j,irot;
  
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  
  if (channel > kNumberOfCalibChannels || data == NULL) {
    for (i = 0; i < kNumberOfBins; i++) {
      uWaveform[i] = adcWaveform[i];
    }
    return true;
  }
  if (!data->fRead) {
    for (i = 0; i < kNumberOfBins; i++) {
      uWaveform[i] = adcWaveform[i];
    }
    return true;
  }
  
  chn = data->fChannel[channel];
  
  hasOffset = data->fHasOffsetCalibration;
  aboveThreshold = (threshold == 0);   // If threshold equal zero, always return true
  
  // Calibrate
  for (i = 0; i < kNumberOfBins; i++) {
    if (fBoard->GetChipVersion() != 3) {
      irot = i;
      if (triggerCell > -1)
        irot = (triggerCell + i) % kNumberOfBins;
      if (adcWaveform[irot] > chn->fLookUpOffset[irot]) {
        uWaveform[i] =
          ((chn->fLookUp[irot][0] - chn->fLookUp[irot][1]) * (adcWaveform[irot] - chn->fLookUpOffset[irot]) +
           chn->fLookUp[irot][0]);
      } else if (adcWaveform[irot] <= chn->fLookUpOffset[irot]
                 && adcWaveform[irot] > chn->fLookUpOffset[irot] - chn->fNumberOfLookUpPoints[irot]) {
        uWaveform[i] = chn->fLookUp[irot][chn->fLookUpOffset[irot] - adcWaveform[irot]];
      } else {
        wave = 0;
        for (j = 0; j < kBSplineOrder; j++) {
          wave += chn->fData[irot][data->fBSplineOffsetLookUp[adcWaveform[irot]][chn->fLimitGroup[irot]] + j]
            * data->fBSplineLookUp[adcWaveform[irot]][chn->fLimitGroup[irot]][j];
        }
        uWaveform[i] = static_cast<short>(wave);
      }
      // Offset Calibration
      if (hasOffset)
        uWaveform[i] -= chn->fOffset[irot];
    } else {
      if (chn->fGain[i] > 3800) {
        //if (channel == 1 && i == 10)
        //   printf("gain:%d offset:%d value:%d\n", chn->fGain[i], chn->fOffset[i], adcWaveform[i]);
        v = static_cast<float>(adcWaveform[i] - chn->fOffset[i]) / chn->fGain[i];
        uWaveform[i] = static_cast<short>(v * 1000 / GetPrecision() + 0.5);
      } else
        uWaveform[i] = 0;
    }
    
    // Check for Threshold
    if (!aboveThreshold) {
      if (uWaveform[i] >= threshold)
        aboveThreshold = true;
    }
  }
  return aboveThreshold;
}

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

bool ResponseCalibration::SubtractADCOffset(unsigned int chipIndex, unsigned int channel, unsigned short *adcWaveform,
                                            unsigned short *adcCalibratedWaveform, unsigned short newBaseLevel)
{
  int i;
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  
  if (channel >= kNumberOfCalibChannels || data == NULL)
    return false;
  if (!data->fRead || !data->fHasOffsetCalibration)
    return false;
  
  chn = data->fChannel[channel];
  for (i = 0; i < kNumberOfBins; i++)
    adcCalibratedWaveform[i] = adcWaveform[i]-chn->fOffsetADC[i]+newBaseLevel;
  return true;
}


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

bool ResponseCalibration::ReadCalibration(unsigned int chipIndex)
{
  if (fBoard->GetChipVersion() == 3)
    return ReadCalibrationV4(chipIndex);
  else
    return ReadCalibrationV3(chipIndex);
}

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

bool ResponseCalibration::ReadCalibrationV3(unsigned int chipIndex)
{
  int k, l, m, num;
  unsigned char ng;
  short tempShort;
  char fileName[2000];
  FILE *fileHandle;
  char calibDir[1000];
  
  // Read Response Calibration
  delete fCalibrationData[chipIndex];
  fCalibrationData[chipIndex] = NULL;
  
  fBoard->GetCalibrationDirectory(calibDir);
  sprintf(fileName, "%s/board%d/ResponseCalib_board%d_chip%d_%dMHz.bin", calibDir,
          fBoard->GetCMCSerialNumber(), fBoard->GetCMCSerialNumber(), chipIndex,
          static_cast<int>(fBoard->GetFrequency() * 1000));
  
  fileHandle = fopen(fileName, "rb");
  if (fileHandle == NULL) {
    printf("Board %d --> Could not find response calibration file:\n", fBoard->GetCMCSerialNumber());
    printf("%s\n", fileName);
    return false;
  }
  // Number Of Grid Points
  num = fread(&ng, 1, 1, fileHandle);
  if (num != 1) {
    printf("Error while reading response calibration file '%s'\n", fileName);
    printf("   at 'NumberOfGridPoints'.\n");
    return false;
  }
  
  fCalibrationData[chipIndex] = new CalibrationData(ng);
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  data->fRead = true;
  data->fHasOffsetCalibration = 1;
  data->DeletePreCalculatedBSpline();
  fCalibrationValid[chipIndex] = true;
  
  // Start Temperature
  num = fread(&tempShort, 2, 1, fileHandle);
  if (num != 1) {
    printf("Error while reading response calibration file '%s'\n", fileName);
    printf("   at 'StartTemperature'.\n");
    return false;
  }
  data->fStartTemperature = static_cast<float>(tempShort) / 10;
  // End Temperature
  num = fread(&tempShort, 2, 1, fileHandle);
  if (num != 1) {
    printf("Error while reading response calibration file '%s'\n", fileName);
    printf("   at 'EndTemperature'.\n");
    return false;
  }
  data->fEndTemperature = static_cast<float>(tempShort) / 10;
  if (fBoard->GetChipVersion() != 3) {
    // Min
    num = fread(&data->fMin, 4, 1, fileHandle);
    if (num != 1) {
      printf("Error while reading response calibration file '%s'\n", fileName);
      printf("   at 'Min'.\n");
      return false;
    }
    // Max
    num = fread(&data->fMax, 4, 1, fileHandle);
    if (num != 1) {
      printf("Error while reading response calibration file '%s'\n", fileName);
      printf("   at 'Max'.\n");
      return false;
    }
    // Number Of Limit Groups
    num = fread(&data->fNumberOfLimitGroups, 1, 1, fileHandle);
    if (num != 1) {
      printf("Error while reading response calibration file '%s'\n", fileName);
      printf("   at 'NumberOfLimitGroups'.\n");
      return false;
    }
  }
  // Read channel
  for (k = 0; k < kNumberOfCalibChannels; k++) {
    chn = data->fChannel[k];
    for (l = 0; l < kNumberOfBins; l++) {
      if (fBoard->GetChipVersion() != 3) {
        // Range Group
        num = fread(&chn->fLimitGroup[l], 1, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'RangeGroup' of channel %d bin %d.\n", k, l);
          return false;
        }
        // Look Up Offset
        num = fread(&chn->fLookUpOffset[l], 2, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'LookUpOffset' of channel %d bin %d.\n", k, l);
          return false;
        }
        // Number Of Look Up Points
        num = fread(&chn->fNumberOfLookUpPoints[l], 1, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'NumberOfLookUpPoints' of channel %d bin %d.\n", k, l);
          return false;
        }
        // Look Up Points
        delete chn->fLookUp[l];
        chn->fLookUp[l] = new unsigned char[chn->fNumberOfLookUpPoints[l]];
        for (m = 0; m < chn->fNumberOfLookUpPoints[l]; m++) {
          num = fread(&chn->fLookUp[l][m], 1, 1, fileHandle);
          if (num != 1) {
            printf("Error while reading response calibration file '%s'\n", fileName);
            printf("   at 'LookUp %d' of channel %d bin %d.\n", m, k, l);
            return false;
          }
        }
        // Points
        for (m = 0; m < data->fNumberOfGridPoints; m++) {
          num = fread(&chn->fData[l][m], 2, 1, fileHandle);
          if (num != 1) {
            printf("Error while reading response calibration file '%s'\n", fileName);
            printf("   at 'Point %d' of channel %d bin %d.\n", m, k, l);
            return false;
          }
        }
        // ADC Offset
        num = fread(&chn->fOffsetADC[l], 2, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'ADC Offset' of channel %d bin %d.\n", k, l);
          return false;
        }
      }
      // Offset
      num = fread(&chn->fOffset[l], 2, 1, fileHandle);
      if (num != 1) {
        printf("Error while reading response calibration file '%s'\n", fileName);
        printf("   at 'Offset' of channel %d bin %d.\n", k, l);
        return false;
      }
      if (fBoard->GetChipVersion() == 3) {
        // Gain
        num = fread(&chn->fGain[l], 2, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'Gain' of channel %d bin %d.\n", k, l);
          return false;
        }
      }
    }
  }
  fclose(fileHandle);
  
  if (fBoard->GetChipVersion() != 3) {
    data->PreCalculateBSpline();
  }
  
  return true;
}

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

bool ResponseCalibration::ReadCalibrationV4(unsigned int chipIndex)
{
  int k, l, num;
  char fileName[2000];
  FILE *fileHandle;
  char calibDir[1000];
  
  // Read Response Calibration
  
  fBoard->GetCalibrationDirectory(calibDir);
  sprintf(fileName, "%s/board%d/ResponseCalib_board%d_chip%d_%dMHz.bin", calibDir,
          fBoard->GetCMCSerialNumber(), fBoard->GetCMCSerialNumber(), chipIndex,
          static_cast<int>(fBoard->GetFrequency() * 1000));
  
  fileHandle = fopen(fileName, "rb");
  if (fileHandle == NULL) {
    printf("Board %d --> Could not find response calibration file:\n", fBoard->GetCMCSerialNumber());
    printf("%s\n", fileName);
    return false;
  }
  
  if (fInitialized)
    delete fCalibrationData[chipIndex];
  fCalibrationData[chipIndex] = new CalibrationData(1);
  CalibrationData *data = fCalibrationData[chipIndex];
  CalibrationData::CalibrationDataChannel * chn;
  data->fRead = true;
  data->fHasOffsetCalibration = 1;
  fCalibrationValid[chipIndex] = true;
  data->fStartTemperature = 0;
  data->fEndTemperature = 0;
  
  // read channel
  for (k = 0; k < kNumberOfCalibChannels; k++) {
    chn = data->fChannel[k];
    for (l = 0; l < kNumberOfBins; l++) {
      // Offset
      num = fread(&chn->fOffset[l], 2, 1, fileHandle);
      if (num != 1) {
        printf("Error while reading response calibration file '%s'\n", fileName);
        printf("   at 'Offset' of channel %d bin %d.\n", k, l);
        return false;
      }
      if (fBoard->GetChipVersion() == 3) {
        // Gain
        num = fread(&chn->fGain[l], 2, 1, fileHandle);
        if (num != 1) {
          printf("Error while reading response calibration file '%s'\n", fileName);
          printf("   at 'Gain' of channel %d bin %d.\n", k, l);
          return false;
        }
      }
    }
  }
  
  fclose(fileHandle);
  return true;
}

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

float ResponseCalibration::GetValue(float *coefficients, float u, int n)
{
  int j, ii;
  float bsplines[4];
  ii = CalibrationData::CalculateBSpline(n, u, bsplines);
  
  float s = 0;
  for (j = 0; j < kBSplineOrder; j++) {
    s += coefficients[ii + j] * bsplines[j];
  }
  return s;
}

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

int ResponseCalibration::Approx(float *p, float *uu, int np, int nu, float *coef)
{
  int i, iu, j;
  
  const int mbloc = 50;
  int ip = 0;
  int ir = 0;
  int mt = 0;
  int ileft, irow;
  float bu[kBSplineOrder];
  float *matrix[kBSplineOrder + 2];
  for (i = 0; i < kBSplineOrder + 2; i++)
    matrix[i] = new float[mbloc + nu + 1];
  for (iu = kBSplineOrder - 1; iu < nu; iu++) {
    for (i = 0; i < np; i++) {
      if (1 <= uu[i])
        ileft = nu - 1;
      else if (uu[i] < 0)
        ileft = kBSplineOrder - 2;
      else
        ileft = kBSplineOrder - 1 + static_cast<int>(uu[i] * (nu - kBSplineOrder + 1));
      if (ileft != iu)
        continue;
      irow = ir + mt;
      mt++;
      CalibrationData::CalculateBSpline(nu, uu[i], bu);
      for (j = 0; j < kBSplineOrder; j++) {
        matrix[j][irow] = bu[j];
      }
      matrix[kBSplineOrder][irow] = p[i];
      if (mt < mbloc)
        continue;
      LeastSquaresAccumulation(matrix, kBSplineOrder, &ip, &ir, mt, iu - kBSplineOrder + 1);
      mt = 0;
    }
    if (mt == 0)
      continue;
    LeastSquaresAccumulation(matrix, kBSplineOrder, &ip, &ir, mt, iu - kBSplineOrder + 1);
    mt = 0;
  }
  if (!LeastSquaresSolving(matrix, kBSplineOrder, ip, ir, coef, nu)) {
    for (i = 0; i < kBSplineOrder + 2; i++)
      delete matrix[i];
    return 0;
  }
  
  for (i = 0; i < kBSplineOrder + 2; i++)
    delete matrix[i];
  return 1;
}

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

void ResponseCalibration::LeastSquaresAccumulation(float **matrix, int nb, int *ip, int *ir, int mt, int jt)
{
  int i, j, l, mu, k, kh;
  float rho;
  
  if (mt <= 0)
    return;
  if (jt != *ip) {
    if (jt > (*ir)) {
      for (i = 0; i < mt; i++) {
        for (j = 0; j < nb + 1; j++) {
          matrix[j][jt + mt - i] = matrix[j][(*ir) + mt - i];
        }
      }
      for (i = 0; i < jt - (*ir); i++) {
        for (j = 0; j < nb + 1; j++) {
          matrix[j][(*ir) + i] = 0;
        }
      }
      *ir = jt;
    }
    mu = min(nb - 1, (*ir) - (*ip) - 1);
    if (mu != 0) {
      for (l = 0; l < mu; l++) {
        k = min(l + 1, jt - (*ip));
        for (i = l + 1; i < nb; i++) {
          matrix[i - k][(*ip) + l + 1] = matrix[i][(*ip) + l + 1];
        }
        for (i = 0; i < k; i++) {
          matrix[nb - i - 1][(*ip) + l + 1] = 0;
        }
      }
    }
    *ip = jt;
  }
  kh = min(nb + 1, (*ir) + mt - (*ip));
  
  for (i = 0; i < kh; i++) {
    Housholder(i, max(i + 1, (*ir) - (*ip)), (*ir) + mt - (*ip), matrix, i, (*ip), &rho, matrix, i + 1,
               (*ip), 1, nb - i);
  }
  
  *ir = (*ip) + kh;
  if (kh < nb + 1)
    return;
  for (i = 0; i < nb; i++) {
    matrix[i][(*ir) - 1] = 0;
  }
}

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

int ResponseCalibration::LeastSquaresSolving(float **matrix, int nb, int ip, int ir, float *x, int n)
{
  int i, j, l, ii;
  float s, rsq;
  for (j = 0; j < n; j++) {
    x[j] = matrix[nb][j];
  }
  rsq = 0;
  if (n <= ir - 1) {
    for (j = n; j < ir; j++) {
      rsq += pow(matrix[nb][j], 2);
    }
  }
  
  for (ii = 0; ii < n; ii++) {
    i = n - ii - 1;
    s = 0;
    l = max(0, i - ip);
    if (i != n - 1) {
      for (j = 1; j < min(n - i, nb); j++) {
        s += matrix[j + l][i] * x[i + j];
      }
    }
    if (matrix[l][i] == 0) {
      printf("Error in LeastSquaresSolving.\n");
      return 0;
    }
    x[i] = (x[i] - s) / matrix[l][i];
  }
  return 1;
}

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

void ResponseCalibration::Housholder(int lpivot, int l1, int m, float **u, int iU1, int iU2, float *up,
                                     float **c, int iC1, int iC2, int ice, int ncv)
{
  int i, j, incr;
  float tol = static_cast<float>(1e-20);
  float tolb = static_cast<float>(1e-24);
  float cl, clinv, sm, b;
  
  if (lpivot < 0 || lpivot >= l1 || l1 > m - 1)
    return;
  cl = fabs(u[iU1][iU2 + lpivot]);
  
  // Construct the transformation
  for (j = l1 - 1; j < m; j++)
    cl = max(fabsf(u[iU1][iU2 + j]), cl);
  if (cl < tol)
    return;
  clinv = 1 / cl;
  sm = pow(u[iU1][iU2 + lpivot] * clinv, 2);
  for (j = l1; j < m; j++) {
    sm = sm + pow(u[iU1][iU2 + j] * clinv, 2);
  }
  cl *= sqrt(sm);
  if (u[iU1][iU2 + lpivot] > 0)
    cl = -cl;
  *up = u[iU1][iU2 + lpivot] - cl;
  u[iU1][iU2 + lpivot] = cl;
  
  if (ncv <= 0)
    return;
  b = (*up) * u[iU1][iU2 + lpivot];
  if (fabs(b) < tolb)
    return;
  if (b >= 0)
    return;
  b = 1 / b;
  incr = ice * (l1 - lpivot);
  for (j = 0; j < ncv; j++) {
    sm = c[iC1 + j][iC2 + lpivot] * (*up);
    for (i = l1; i < m; i++) {
      sm = sm + c[iC1 + j][iC2 + lpivot + incr + (i - l1) * ice] * u[iU1][iU2 + i];
    }
    if (sm == 0)
      continue;
    sm *= b;
    c[iC1 + j][iC2 + lpivot] = c[iC1 + j][iC2 + lpivot] + sm * (*up);
    for (i = l1; i < m; i++) {
      c[iC1 + j][iC2 + lpivot + incr + (i - l1) * ice] =
        c[iC1 + j][iC2 + lpivot + incr + (i - l1) * ice] + sm * u[iU1][iU2 + i];
    }
  }
}

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

int ResponseCalibration::MakeDir(const char *path)
{
  struct stat buf;
  if (stat(path, &buf)) 
    return mkdir(path, 0711);
  return 0;
}

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

ResponseCalibration::ResponseCalibration(DRSBoard *board)
  :fBoard(board)
   ,fPrecision(0.1) // mV
   ,fInitialized(false)
   ,fRecorded(false)
   ,fFitted(false)
   ,fOffset(false)
   ,fNumberOfPointsLowVolt(0)
   ,fNumberOfPoints(0)
   ,fNumberOfMode2Bins(0)
   ,fNumberOfSamples(0)
   ,fNumberOfGridPoints(0)
   ,fNumberOfXConstPoints(0)
   ,fNumberOfXConstGridPoints(0)
   ,fTriggerFrequency(0)
   ,fShowStatistics(0)
   ,fCalibFile(0)
   ,fCurrentLowVoltPoint(0)
   ,fCurrentPoint(0)
   ,fCurrentSample(0)
   ,fCurrentFitChannel(0)
   ,fCurrentFitBin(0)
   ,fResponseY(0)
   ,fSamples(0)
   ,fSampleUsed(0)
   ,fXXFit(0)
   ,fYYFit(0)
   ,fWWFit(0)
   ,fYYFitRes(0)
   ,fYYSave(0)
   ,fXXSave(0)
   ,fStatisticsApprox(0)
   ,fStatisticsApproxExt(0)
{
  int i;
  // Initializing the Calibration Class
  CalibrationData::fIntRevers[0] = 0;
  for (i = 1; i < 2 * kBSplineOrder - 2; i++) {
    CalibrationData::fIntRevers[i] = static_cast<float>(1.) / i;
  }
  for (i = 0; i < kNumberOfChips; i++) {
    fCalibrationData[i] = NULL;
  }
  // Initializing the Calibration Creation
  fCalibrationValid[0] = false;
  fCalibrationValid[1] = false;
}

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

ResponseCalibration::~ResponseCalibration()
{
  // Deleting the Calibration Creation
  DeleteFields();
}

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

float ResponseCalibration::CalibrationData::fIntRevers[2 * kBSplineOrder - 2];
ResponseCalibration::CalibrationData::CalibrationData(int numberOfGridPoints)
  :fRead(false)
   ,fNumberOfGridPoints(numberOfGridPoints)
   ,fHasOffsetCalibration(0)
   ,fStartTemperature(0)
   ,fEndTemperature(0)
   ,fMin(0)
   ,fMax(0)
   ,fNumberOfLimitGroups(0)
{
  int i;
  for (i = 0; i < kNumberOfCalibChannels; i++) {
    fChannel[i] = new CalibrationDataChannel(numberOfGridPoints);
  }
  for (i = 0; i < kNumberOfADCBins; i++) {
    fBSplineOffsetLookUp[i] = NULL;
    fBSplineLookUp[i] = NULL;
  }
};

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

void ResponseCalibration::CalibrationData::PreCalculateBSpline()
{
  int i, j;
  float uu;
  float xmin, xrange;
  int nk = fNumberOfGridPoints - kBSplineOrder + 1;
  for (i = 0; i < kNumberOfADCBins; i++) {
    fBSplineLookUp[i] = new float *[fNumberOfLimitGroups];
    fBSplineOffsetLookUp[i] = new int[fNumberOfLimitGroups];
    for (j = 0; j < fNumberOfLimitGroups; j++) {
      fBSplineLookUp[i][j] = new float[kBSplineOrder];
      xmin = fMin + j * kBSplineXMinOffset;
      xrange = fMax - xmin;
      uu = (i - xmin) / xrange;
      if (i < xmin) {
        uu = 0;
      }
      if (i - xmin > xrange) {
        uu = 1;
      }
      fBSplineOffsetLookUp[i][j] = static_cast<int>(uu * nk);
      CalculateBSpline(fNumberOfGridPoints, uu, fBSplineLookUp[i][j]);
    }
  }
}

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

void ResponseCalibration::CalibrationData::DeletePreCalculatedBSpline()
{
  int i, j;
  for (i = 0; i < kNumberOfADCBins; i++) {
    if (fBSplineLookUp[i]!=NULL) {
      for (j = 0; j < fNumberOfLimitGroups; j++)
        delete fBSplineLookUp[i][j];
    }
    delete fBSplineLookUp[i];
    delete fBSplineOffsetLookUp[i];
  }
}

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

ResponseCalibration::CalibrationData::~CalibrationData()
{
  int i, j;
  for (i = 0; i < kNumberOfCalibChannels; i++) {
    delete fChannel[i];
  }
  for (i = 0; i < kNumberOfADCBins; i++) {
    if (fBSplineLookUp[i]!=NULL) {
      for (j = 0; j < fNumberOfLimitGroups; j++) {
        delete fBSplineLookUp[i][j];
      }
    }
    delete fBSplineLookUp[i];
    delete fBSplineOffsetLookUp[i];
  }
};

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

int ResponseCalibration::CalibrationData::CalculateBSpline(int nGrid, float value, float *bsplines)
{
  int minimum;
  int maximum;
  float xl;
  
  int nk = nGrid - kBSplineOrder + 1;
  float vl = value * nk;
  int ivl = static_cast<int>(vl);
  
  if (1 <= value) {
    xl = vl - nk + 1;
    minimum = 1 - nk;
  } else if (value < 0) {
    xl = vl;
    minimum = 0;
  } else {
    xl = vl - ivl;
    minimum = -ivl;
  }
  maximum = nk + minimum;
  
  // printf("xl = %f\n",xl);
  float vm, vmprev;
  int jl, ju;
  int nb = 0;
  
  bsplines[0] = 1;
  for (int i = 0; i < kBSplineOrder - 1; i++) {
    vmprev = 0;
    for (int j = 0; j < nb + 1; j++) {
      jl = max(minimum, j - nb);
      ju = min(maximum, j + 1);
      vm = bsplines[j] * fIntRevers[ju - jl];
      bsplines[j] = vm * (ju - xl) + vmprev;
      vmprev = vm * (xl - jl);
    }
    nb++;
    bsplines[nb] = vmprev;
  }
  return -minimum;
}

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

void ResponseCalibration::Average(int method,float *points,int numberOfPoints,float &mean,float &error,float sigmaBoundary)
{
  // Methods :
  // 0 : Average
  // 1 : Average inside sigmaBoundary*sigma
  int i;
  float sum = 0;
  float sumSquare = 0;
  
  if (method == 0 || method == 1) {
    for (i = 0; i < numberOfPoints; i++) {
      sum += points[i];
      sumSquare += points[i]*points[i];
    }
    
    mean = sum / numberOfPoints;
    error = sqrt((sumSquare - sum * sum / numberOfPoints) / (numberOfPoints - 1));
  }
  if (method == 1) {
    int numberOfGoodPoints = numberOfPoints;
    bool found = true;
    bool *goodSample = new bool[numberOfGoodPoints];
    for (i = 0; i < numberOfGoodPoints; i++)
      goodSample[i] = true;
    
    while (found) {
      found = false;
      for (i = 0; i < numberOfPoints; i++) {
        if (goodSample[i] && fabs(points[i] - mean) > sigmaBoundary * error) {
          found = true;
          goodSample[i] = false;
          numberOfGoodPoints--;
          sum -= points[i];
          sumSquare -= points[i]*points[i];
          mean = sum/numberOfGoodPoints;
          error = sqrt((sumSquare - sum * sum / numberOfGoodPoints) / (numberOfGoodPoints - 1));
        }
      }
    }
    delete goodSample;
  }
}