/********************************************************************

  Name:         DRS.cpp
  Created by:   Stefan Ritt, Matthias Schneebeli

  Contents:     Library functions for DRS mezzanine and USB boards

  $Id: DRS.cpp 13351 2009-04-28 11:12:54Z ritt@PSI.CH $

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

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <assert.h>
#include <algorithm>
#include <sys/stat.h>
#include "strlcpy.h"

#ifdef _MSC_VER
#pragma warning(disable:4996)
#   include <windows.h>
#   include <direct.h>
#else
#   include <unistd.h>
#   include <sys/time.h>
inline void Sleep(useconds_t x)
{
   usleep(x * 1000);
}
#endif

#ifdef _MSC_VER
#include <conio.h>
#define drs_kbhit() kbhit()
#else
#include <sys/ioctl.h>
int drs_kbhit()
{
   int n;

   ioctl(0, FIONREAD, &n);
   return (n > 0);
}
static inline int getch()
{
   return getchar();
}
#endif

#include <DRS.h>

#ifdef _MSC_VER
extern "C" {
#endif

#include <mxml.h>

#ifdef _MSC_VER
}
#endif

/*---- minimal FPGA firmvare version required for this library -----*/
const int REQUIRED_FIRMWARE_VERSION_DRS2 = 5268;
const int REQUIRED_FIRMWARE_VERSION_DRS3 = 6981;
const int REQUIRED_FIRMWARE_VERSION_DRS4 = 13191;

/*---- VME addresses -----------------------------------------------*/
#ifdef HAVE_VME
/* 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)        reserverd
   SW1-6: 0 (on)        reserverd
   SW1-7: 0 (on)        reserverd
   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
#endif                          // HAVE_VME
/*---- USB addresses -----------------------------------------------*/
#define USB_TIMEOUT                     1000    // one second
#ifdef HAVE_USB
#define USB_CTRL_OFFSET                 0x00    /* all registers 32 bit */
#define USB_STATUS_OFFSET               0x40
#define USB_RAM_OFFSET                  0x80
#define USB_CMD_IDENT                      0    // Query identification
#define USB_CMD_ADDR                       1    // Address cycle
#define USB_CMD_READ                       2    // "VME" read <addr><size>
#define USB_CMD_WRITE                      3    // "VME" write <addr><size>
#define USB_CMD_READ12                     4    // 12-bit read <LSB><MSB>
#define USB_CMD_WRITE12                    5    // 12-bit write <LSB><MSB>

#define USB2_CMD_READ                      1
#define USB2_CMD_WRITE                     2
#define USB2_CTRL_OFFSET             0x00000    /* all registers 32 bit */
#define USB2_STATUS_OFFSET           0x10000
#define USB2_FIFO_OFFSET             0x20000
#define USB2_RAM_OFFSET              0x40000
#endif                          // HAVE_USB

/*---- Register addresses ------------------------------------------*/

#ifndef T_CTRL
#define T_CTRL                             1
#define T_STATUS                           2
#define T_RAM                              3
#define T_FIFO                             4
#endif

#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    // low byte
#define REG_CONFIG                   0x00014    // high byte
#define REG_CHANNEL_SPAN             0x00016
#define REG_FREQ_SET_HI              0x00018    // DRS2
#define REG_FREQ_SET_LO              0x0001A    // DRS2
#define REG_TRG_DELAY                0x00018    // DRS4
#define REG_FREQ_SET                 0x0001A    // DRS4
#define REG_TRIG_DELAY               0x0001C
#define REG_LMK_MSB                  0x0001C    // DRS4 Mezz
#define REG_CALIB_TIMING             0x0001E    // DRS2
#define REG_EEPROM_PAGE              0x0001E    // DRS4
#define REG_LMK_LSB                  0x0001E    // DRS4 Mezz

#define REG_MAGIC                    0x00000
#define REG_BOARD_TYPE               0x00002
#define REG_STATUS                   0x00004
#define REG_RDAC_OFS                 0x0000E
#define REG_RDAC0                    0x00008
#define REG_STOP_CELL0               0x00008
#define REG_RDAC1                    0x0000A
#define REG_STOP_CELL1               0x0000A
#define REG_RDAC2                    0x0000C
#define REG_STOP_CELL2               0x0000C
#define REG_RDAC3                    0x0000E
#define REG_STOP_CELL3               0x0000E
#define REG_RDAC4                    0x00000
#define REG_RDAC5                    0x00002
#define REG_RDAC6                    0x00014
#define REG_RDAC7                    0x00016
#define REG_EVENTS_IN_FIFO           0x00018
#define REG_EVENT_COUNT              0x0001A
#define REG_FREQ1                    0x0001C
#define REG_FREQ2                    0x0001E
#define REG_TEMPERATURE              0x00020
#define REG_TRIGGER_BUS              0x00022
#define REG_SERIAL_BOARD             0x00024
#define REG_VERSION_FW               0x00026

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

using namespace std;

#ifdef HAVE_USB
#define USB2_BUFFER_SIZE (1024*1024+10)
unsigned char static *usb2_buffer = NULL;
#endif

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

DRS::DRS()
:  fNumberOfBoards(0)
#ifdef HAVE_VME
    , fVmeInterface(0)
#endif
{
#ifdef HAVE_USB
   MUSB_INTERFACE *usb_interface;
#endif

#if defined(HAVE_VME) || defined(HAVE_USB)
   int index = 0, i = 0;
#endif

   memset(fError, 0, sizeof(fError));

#ifdef HAVE_VME
   unsigned short type, fw, magic, serial, temperature;
   mvme_addr_t addr;

   if (mvme_open(&fVmeInterface, 0) == MVME_SUCCESS) {

      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 */
         addr = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE;        // VME board base address
         addr += GEVPC_USER_FPGA;       // UsrFPGA base address
         addr += PMC1_OFFSET;   // PMC1 offset

         mvme_set_dmode(fVmeInterface, MVME_DMODE_D16);
         i = mvme_read(fVmeInterface, &magic, addr + PMC_STATUS_OFFSET + REG_MAGIC, 2);
         if (i == 2) {
            if (magic != 0xC0DE) {
               printf("Found old firmware, please upgrade immediately!\n");
               fBoard[fNumberOfBoards] = new DRSBoard(fVmeInterface, addr, (index - 2) << 1);
               fNumberOfBoards++;
            } else {

               /* read board type */
               mvme_read(fVmeInterface, &type, addr + PMC_STATUS_OFFSET + REG_BOARD_TYPE, 2);
               type &= 0xFF;
               if (type == 2 || type == 3 || type == 4) {    // DRS2 or DRS3 or DRS4

                  /* read firmware number */
                  mvme_read(fVmeInterface, &fw, addr + PMC_STATUS_OFFSET + REG_VERSION_FW, 2);

                  /* read serial number */
                  mvme_read(fVmeInterface, &serial, addr + PMC_STATUS_OFFSET + REG_SERIAL_BOARD, 2);

                  /* read temperature register to see if CMC card is present */
                  mvme_read(fVmeInterface, &temperature, addr + PMC_STATUS_OFFSET + REG_TEMPERATURE, 2);

                  /* LED blinking */
#if 0
                  do {
                     data = 0x00040000;
                     mvme_write(fVmeInterface, addr + PMC_CTRL_OFFSET + REG_CTRL, &data, sizeof(data));
                     mvme_write(fVmeInterface, addr + PMC2_OFFSET + PMC_CTRL_OFFSET + REG_CTRL, &data,
                                sizeof(data));

                     Sleep(500);

                     data = 0x00000000;
                     mvme_write(fVmeInterface, addr + PMC_CTRL_OFFSET + REG_CTRL, &data, sizeof(data));
                     mvme_write(fVmeInterface, addr + PMC2_OFFSET + PMC_CTRL_OFFSET + REG_CTRL, data,
                                sizeof(data));

                     Sleep(500);

                  } while (1);
#endif

                  if (temperature == 0xFFFF) {
                     //printf("slot %d, fw %d, no CMC board in upper slot\n", index, fw);
                  } else {
                     //printf("slot %d, fw %d, CMC serial %d in upper slot\n", index, fw, serial);

                     fBoard[fNumberOfBoards] = new DRSBoard(fVmeInterface, addr, (index - 2) << 1);
                     if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
                        fNumberOfBoards++;
                     else
                        sprintf(fError, "Wrong firmware version: board has %d, required is %d\n",
                                fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                                fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
                  }
               }
            }
         }

         /* check PMC2 */
         addr = GEVPC_BASE_ADDR + index * GEVPC_WINSIZE;        // VME board base address
         addr += GEVPC_USER_FPGA;       // UsrFPGA base address
         addr += PMC2_OFFSET;   // PMC2 offset

         mvme_set_dmode(fVmeInterface, MVME_DMODE_D16);
         i = mvme_read(fVmeInterface, &fw, addr + PMC_STATUS_OFFSET + REG_MAGIC, 2);
         if (i == 2) {
            if (magic != 0xC0DE) {
               printf("Found old firmware, please upgrade immediately!\n");
               fBoard[fNumberOfBoards] = new DRSBoard(fVmeInterface, addr, (index - 2) << 1 | 1);
               fNumberOfBoards++;
            } else {

               /* read board type */
               mvme_read(fVmeInterface, &type, addr + PMC_STATUS_OFFSET + REG_BOARD_TYPE, 2);
               type &= 0xFF;
               if (type == 2 || type == 3 || type == 4) {    // DRS2 or DRS3 or DRS4

                  /* read firmware number */
                  mvme_read(fVmeInterface, &fw, addr + PMC_STATUS_OFFSET + REG_VERSION_FW, 2);

                  /* read serial number */
                  mvme_read(fVmeInterface, &serial, addr + PMC_STATUS_OFFSET + REG_SERIAL_BOARD, 2);

                  /* read temperature register to see if CMC card is present */
                  mvme_read(fVmeInterface, &temperature, addr + PMC_STATUS_OFFSET + REG_TEMPERATURE, 2);

                  if (temperature == 0xFFFF) {
                     //printf("slot %d, fw %d, no CMC board in lower slot\n", index, fw);
                  } else {
                     //printf("slot %d, fw %d, CMC serial %d in lower slot\n", index, fw, serial);

                     fBoard[fNumberOfBoards] = new DRSBoard(fVmeInterface, addr, ((index - 2) << 1) | 1);
                     if (fBoard[fNumberOfBoards]->HasCorrectFirmware())
                        fNumberOfBoards++;
                     else
                        sprintf(fError, "Wrong firmware version: board has %d, required is %d\n",
                                fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                                fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
                  }
               }
            }
         }
      }
   } else
      printf("Cannot access VME crate, check driver, power and connection\n");
#endif                          // HAVE_VME

#ifdef HAVE_USB
   unsigned char buffer[512];
   int found, one_found;

   one_found = 0;
   for (index = 0; index < 127; index++) {
      found = 0;

      /* check for USB-Mezzanine test board */
      if (musb_open(&usb_interface, 0x10C4, 0x1175, index, 1, 0) == MUSB_SUCCESS) {

         /* check ID */
         buffer[0] = USB_CMD_IDENT;
         musb_write(usb_interface, 2, buffer, 1, USB_TIMEOUT);

         i = musb_read(usb_interface, 1, (char *) buffer, sizeof(buffer), USB_TIMEOUT);
         if (strcmp((char *) buffer, "USB_MEZZ2 V1.0") != 0) {
            /* no USB-Mezzanine board found */
            musb_close(usb_interface);
         } else {
            usb_interface->usb_type = 1;        // USB 1.1
            fBoard[fNumberOfBoards] = new DRSBoard(usb_interface);
            if (fBoard[fNumberOfBoards]->HasCorrectFirmware()) {
               fNumberOfBoards++;
               found = 1;
               one_found = 1;
            } else
               sprintf(fError, "Wrong firmware version: board has %d, required is %d\n",
                       fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                       fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
         }
      }

      /* check for DRS4 evaluation board */
      if (musb_open(&usb_interface, 0x04B4, 0x1175, index, 1, 0) == MUSB_SUCCESS) {

         /* check ID */
         struct usb_device_descriptor d;
         usb_get_descriptor(usb_interface->dev, USB_DT_DEVICE, 0, &d, sizeof(d));
         if (d.bcdDevice != 1) {
            /* no DRS evaluation board found */
            musb_close(usb_interface);
         } else {

            /* drain any data from Cy7C68013 FIFO if FPGA startup caused erratic write */
            do {
               i = musb_read(usb_interface, 8, buffer, sizeof(buffer), 100);
               if (i > 0)
                  printf("%d bytes stuck in buffer\n", i);
            } while (i > 0);

            usb_interface->usb_type = 2;        // USB 2.0
            fBoard[fNumberOfBoards] = new DRSBoard(usb_interface);
            if (fBoard[fNumberOfBoards]->HasCorrectFirmware()) {
               fNumberOfBoards++;
               found = 1;
               one_found = 1;
            } else {
               sprintf(fError, "Wrong firmware version: board has %d, required is %d\n",
                      fBoard[fNumberOfBoards]->GetFirmwareVersion(),
                      fBoard[fNumberOfBoards]->GetRequiredFirmwareVersion());
            }
         }
      }

      if (!found) {
         if (!one_found)
            printf("USB successfully scanned, but no boards found\n");
         break;
      }
   }
#endif                          // HAVE_USB

   return;
}

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

DRS::~DRS()
{
   int i;
   for (i = 0; i < fNumberOfBoards; i++) {
      delete fBoard[i];
   }

#ifdef HAVE_USB
   if (usb2_buffer) {
      free(usb2_buffer);
      usb2_buffer = NULL;
   }
#endif

#ifdef HAVE_VME
   mvme_close(fVmeInterface);
#endif
}

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

bool DRS::GetError(char *str, int size)
{
   if (fError[0])
      strlcpy(str, fError, size);

   return fError[0] > 0;
}

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

#ifdef HAVE_USB
DRSBoard::DRSBoard(MUSB_INTERFACE * musb_interface)
:  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)
    , fDRSType(0)
    , fBoardType(0)
    , fRequiredFirmwareVersion(0)
    , fFirmwareVersion(0)
    , fBoardSerialNumber(0)
    , fCtrlBits(0)
    , fNumberOfReadoutChannels(0)
    , fExternalClockFrequency(0)
#ifdef HAVE_USB
    , fUsbInterface(musb_interface)
#endif
#ifdef HAVE_VME
    , fVmeInterface(0)
    , fBaseAddress(0)
    , fSlotNumber(0)
#endif
    , fFrequency(0)
    , fDominoMode(0)
    , fDominoActive(0)
    , fChannelConfig(0)
    , fWSRLoop(0)
    , fReadoutMode(0)
    , fTriggerEnable1(0)
    , fTriggerEnable2(0)
    , fTriggerSource(0)
    , fDelayedStart(0)
    , fCommonMode(0.8)
    , fAcalMode(0)
    , fAcalVolt(0)
    , fTcalMode(0)
    , fTcalLevel(0)
    , fMaxChips(0)
    , fResponseCalibration(0)
    , fCellCalibrationValid(false)
    , fTimingCalibrationValid(false)
    , fTimeData(0)
    , fNumberOfTimeData(0)
    , fDebug(0)
    , fTriggerStartBin(0)
{
   if (musb_interface->usb_type == 1)
      fTransport = TR_USB;
   else
      fTransport = TR_USB2;
   memset(fTriggerCell, 0, sizeof(fTriggerCell));
   ConstructBoard();
}

#endif

#ifdef HAVE_VME
/*------------------------------------------------------------------*/

DRSBoard::DRSBoard(MVME_INTERFACE * mvme_interface, mvme_addr_t base_address, int slot_number)
: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)
, fDRSType(0)
, fBoardType(0)
, fRequiredFirmwareVersion(0)
, fFirmwareVersion(0)
, fBoardSerialNumber(0)
, fTransport(TR_VME)
, fCtrlBits(0)
, fNumberOfReadoutChannels(0)
, fExternalClockFrequency(0)
#ifdef HAVE_USB
, fUsbInterface(0)
#endif
#ifdef HAVE_VME
, fVmeInterface(mvme_interface)
, fBaseAddress(base_address)
, fSlotNumber(slot_number)
#endif
, fFrequency(0)
, fRefClock(0)
, fDominoMode(1)
, fDominoActive(1)
, fChannelConfig(0)
, fWSRLoop(1)
, fReadoutMode(0)
, fTriggerEnable1(0)
, fTriggerEnable2(0)
, fTriggerSource(0)
, fDelayedStart(0)
, fCommonMode(0.8)
, fMaxChips(0)
, fResponseCalibration(0)
, fTimeData(0)
, fNumberOfTimeData(0)
, fDebug(0)
, fTriggerStartBin(0)
{
   ConstructBoard();
}

#endif

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

DRSBoard::~DRSBoard()
{
   int i;
#ifdef HAVE_USB
   if (fTransport == TR_USB || fTransport == TR_USB2)
      musb_close(fUsbInterface);
#endif

   // Response Calibration
   delete fResponseCalibration;

   // Time Calibration
   for (i = 0; i < fNumberOfTimeData; i++) {
      delete fTimeData[i];
   }
   delete[]fTimeData;
}

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

void DRSBoard::ConstructBoard()
{
   fDebug = 0;
   fWSRLoop = 1;
   fCtrlBits = 0;

   fExternalClockFrequency = 1000. / 30.;
   for (int i = 0; i < kNumberOfChips * kNumberOfChannelsV2; i++)
      fWaveTransferred[i] = false;
   strcpy(fCalibDirectory, ".");

   ReadSerialNumber();

   /* check for required firmware version */
   if (!HasCorrectFirmware())
      return;

   /* set correct reference clock */
   if (fBoardType == 5)
      fRefClock = 60;
   else
      fRefClock = 33;

   /* get mode from hardware */
   if (fDRSType == 4) {
      fDominoMode = (GetConfigReg() & BIT_CONFIG_DMODE) > 0;
   } else {
      fDominoMode = (GetCtrlReg() & BIT_DMODE) > 0;
   }
   fTriggerEnable1 = (GetConfigReg() & BIT_ENABLE_TRIGGER1) > 0;
   fTriggerEnable2 = (GetConfigReg() & BIT_ENABLE_TRIGGER2) > 0;
   fTriggerSource = ((GetConfigReg() & BIT_TR_SOURCE1) > 0) | (((GetConfigReg() & BIT_TR_SOURCE2) > 0) << 1);
   fReadoutMode = (GetConfigReg() & BIT_READOUT_MODE) > 0;
   fDominoActive = (GetConfigReg() & BIT_DACTIVE) > 0;
   ReadFrequency(0, &fFrequency);

   /* initialize number of channels */
   if (fDRSType == 4)
      fNumberOfReadoutChannels = 9;
   else
      fNumberOfReadoutChannels = 10;

   if (fBoardType == 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 (fBoardType == 2 || fBoardType == 3) {
      fDAC_COFS = 0;
      fDAC_DSA = 1;
      fDAC_DSB = 2;
      fDAC_TLEVEL = 3;
      fDAC_CLKOFS = 5;
      fDAC_ACALIB = 6;
      fDAC_ADCOFS = 7;
   } else if (fBoardType == 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;
   } else if (fBoardType == 5) {
      fDAC_ROFS_1 = 0;
      fDAC_CMOFS = 1;
      fDAC_CALN = 2;
      fDAC_CALP = 3;
      fDAC_BIAS = 4;
      fDAC_TLEVEL = 5;
      fDAC_ONOFS = 6;
   } else if (fBoardType == 6) {
      fDAC_ONOFS = 0;
      fDAC_CMOFS = 1;
      fDAC_CALN = 2;
      fDAC_CALP = 3;
      fDAC_ROFS_1 = 6;
      fDAC_BIAS = 7;
   }
   
   if (fDRSType == 4) {

      ReadCalibration();

   } else {
      // Response Calibration
      fResponseCalibration = new ResponseCalibration(this);

      // Time Calibration
      fTimeData = new DRSBoard::TimeData *[kNumberOfChips];
      fNumberOfTimeData = 0;
   }
}

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

void DRSBoard::ReadSerialNumber()
{
   unsigned char buffer[2];
   int number;

   // check magic number
   if (Read(T_STATUS, buffer, REG_MAGIC, 2) < 0) {
      printf("Cannot read from board\n");
      return;
   }

   number = (static_cast < int >(buffer[1]) << 8) +buffer[0];
   if (number != 0xC0DE) {
      printf("Invalid magic number: %04X\n", number);
      return;
   }
   
   // read board type
   Read(T_STATUS, buffer, REG_BOARD_TYPE, 2);
   fDRSType = buffer[0];
   fBoardType = buffer[1];

   // read firmware version
   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_BOARD, 2);
   number = (static_cast < int >(buffer[1]) << 8) +buffer[0];
   fBoardSerialNumber = number;

   // determine DRS type and board type for old boards from setial number
   if (fBoardType == 0) {
      // determine board version from serial number
      if (number >= 2000) {
         fBoardType = 5;
         fDRSType = 4;
      } else if (number >= 1000)
         fBoardType = 4;
      else if (number >= 100)
         fBoardType = 3;
      else if (number > 0)
         fBoardType = 2;
      else {
         fBoardType = 4;
         fDRSType = 3;
         fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS3;
      }
   }
   // retrieve firmware version
   if (fDRSType == 2)
      fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS2;
   if (fDRSType == 3)
      fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS3;
   if (fDRSType == 4)
      fRequiredFirmwareVersion = REQUIRED_FIRMWARE_VERSION_DRS4;

}

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

void DRSBoard::ReadCalibration(void)
{
   unsigned short buf[1024*9*2];
   int i, j;

   /* read offsets and gain from eeprom */
   ReadEEPROM(1, buf, 1024*32);
   for (i=0 ; i<8 ; i++)
      for (j=0 ; j<1024; j++) {
         fCellOffset[i][j] = buf[(i*1024+j)*2];
         fCellGain[i][j]   = buf[(i*1024+j)*2 + 1]/65535.0*0.3+0.8;
      }

   ReadEEPROM(2, buf, 1024*4);
   for (i=0 ; i<1 ; i++)
      for (j=0 ; j<1024; j++) {
         fCellOffset[i+8][j] = buf[(i*1024+j)*2];
         fCellGain[i+8][j]   = buf[(i*1024+j)*2 + 1]/65535.0*0.3+0.8;
      }

   if (abs(fCellOffset[0][0] - 32768) < 2000 &&
       abs(fCellOffset[7][0] - 32768) < 2000)
      fCellCalibrationValid = true;

   /* read timing calibration from eeprom */
   ReadEEPROM(0, buf, 1024*4);

   fCellT[0][0] = 0;
   for (i=1 ; i<1024; i++)
      fCellT[0][i] = fCellT[0][i-1] + buf[i*2+1]/10000.0;
   if (fCellT[0][1023] < 210 && fCellT[0][1023] > 195)
      fTimingCalibrationValid = true;
   if (!fTimingCalibrationValid)
      for (i=0 ; i<1024 ; i++)
        fCellT[0][i] = 1/fFrequency*i;
}

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

bool DRSBoard::HasCorrectFirmware()
{
   /* check for required firmware version */
   return (fFirmwareVersion >= fRequiredFirmwareVersion);
}

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

int DRSBoard::InitFPGA(void)
{

#ifdef HAVE_USB
   if (fTransport == TR_USB2) {
      unsigned char buffer[1];
      int i, status;

      /* blink Cy7C68013A LED and issue an FPGA reset */
      buffer[0] = 0;            // LED off
      musb_write(fUsbInterface, 1, buffer, 1, 100);
      Sleep(50);

      buffer[0] = 1;            // LED on
      musb_write(fUsbInterface, 1, buffer, 1, 100);

      /* wait until EEPROM page #0 has been read */
      for (i=0 ; i<100 ; i++) {
         Read(T_STATUS, &status, REG_STATUS, 4);
         if ((status & BIT_EEPROM_BUSY) == 0)
            break;
         Sleep(10);
      }
   }
#endif

   return 1;
}

int DRSBoard::Write(int type, unsigned int addr, void *data, int size)
{
   // Generic write function

   if (fTransport == TR_VME) {

#ifdef HAVE_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);
      }

      return size;
#endif                          // HAVE_VME

   } else if (fTransport == TR_USB) {
#ifdef HAVE_USB
      unsigned char buffer[64], ack;
      unsigned int base_addr;
      int i, j, n;

      if (type == T_CTRL)
         base_addr = USB_CTRL_OFFSET;
      else if (type == T_STATUS)
         base_addr = USB_STATUS_OFFSET;
      else if (type == T_RAM)
         base_addr = USB_RAM_OFFSET;
      else
         base_addr = 0;

      if (type != T_RAM) {

         /*---- register access ----*/

         if (size == 2) {
            /* word swapping: first 16 bit sit at upper address */
            if ((addr % 4) == 0)
               addr = addr + 2;
            else
               addr = addr - 2;
         }

         buffer[0] = USB_CMD_WRITE;
         buffer[1] = base_addr + addr;
         buffer[2] = size;

         for (i = 0; i < size; i++)
            buffer[3 + i] = *((unsigned char *) data + i);

         /* try 10 times */
         ack = 0;
         for (i = 0; i < 10; i++) {
            n = musb_write(fUsbInterface, 2, buffer, 3 + size, USB_TIMEOUT);
            if (n == 3 + size) {
               for (j = 0; j < 10; j++) {
                  /* wait for acknowledge */
                  n = musb_read(fUsbInterface, 1, &ack, 1, USB_TIMEOUT);
                  if (n == 1 && ack == 1)
                     break;

                  printf("Redo receive\n");
               }
            }

            if (ack == 1)
               return size;

            printf("Redo send\n");
         }
      } else {

         /*---- RAM access ----*/

         buffer[0] = USB_CMD_ADDR;
         buffer[1] = base_addr + addr;
         musb_write(fUsbInterface, 2, buffer, 2, USB_TIMEOUT);

         /* chop buffer into 60-byte packets */
         for (i = 0; i <= (size - 1) / 60; i++) {
            n = size - i * 60;
            if (n > 60)
               n = 60;
            buffer[0] = USB_CMD_WRITE12;
            buffer[1] = n;

            for (j = 0; j < n; j++)
               buffer[2 + j] = *((unsigned char *) data + j + i * 60);

            musb_write(fUsbInterface, 2, buffer, 2 + n, USB_TIMEOUT);

            for (j = 0; j < 10; j++) {
               /* wait for acknowledge */
               n = musb_read(fUsbInterface, 1, &ack, 1, USB_TIMEOUT);
               if (n == 1 && ack == 1)
                  break;

               printf("Redo receive acknowledge\n");
            }
         }

         return size;
      }
#endif                          // HAVE_USB
   } else if (fTransport == TR_USB2) {
#ifdef HAVE_USB
      unsigned int base_addr;
      int i;

      if (usb2_buffer == NULL)
         usb2_buffer = (unsigned char *) malloc(USB2_BUFFER_SIZE);
      assert(usb2_buffer);

      /* only accept even address and number of bytes */
      assert(addr % 2 == 0);
      assert(size % 2 == 0);

      /* check for maximum size */
      assert(size <= USB2_BUFFER_SIZE - 10);

      if (type == T_CTRL)
         base_addr = USB2_CTRL_OFFSET;
      else if (type == T_STATUS)
         base_addr = USB2_STATUS_OFFSET;
      else if (type == T_FIFO)
         base_addr = USB2_FIFO_OFFSET;
      else if (type == T_RAM)
         base_addr = USB2_RAM_OFFSET;
      else
         base_addr = 0;

      if (type != T_RAM && size == 2) {
         /* word swapping: first 16 bit sit at upper address */
         if ((addr % 4) == 0)
            addr = addr + 2;
         else
            addr = addr - 2;
      }

      addr += base_addr;

      usb2_buffer[0] = USB2_CMD_WRITE;
      usb2_buffer[1] = 0;

      usb2_buffer[2] = (addr >> 0) & 0xFF;
      usb2_buffer[3] = (addr >> 8) & 0xFF;
      usb2_buffer[4] = (addr >> 16) & 0xFF;
      usb2_buffer[5] = (addr >> 24) & 0xFF;

      usb2_buffer[6] = (size >> 0) & 0xFF;
      usb2_buffer[7] = (size >> 8) & 0xFF;
      usb2_buffer[8] = (size >> 16) & 0xFF;
      usb2_buffer[9] = (size >> 24) & 0xFF;

      for (i = 0; i < size; i++)
         usb2_buffer[10 + i] = *((unsigned char *) data + i);

      i = musb_write(fUsbInterface, 4, usb2_buffer, 10 + size, USB_TIMEOUT);
      if (i != 10 + size)
         printf("musb_write error: %d\n", i);

      return i;
#endif                          // HAVE_USB
   }

   return 0;
}

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

int DRSBoard::Read(int type, void *data, unsigned int addr, int size)
{
   // Generic read function

   if (fTransport == TR_VME) {

#ifdef HAVE_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);   // 2eVME 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                          // HAVE_VME
   } else if (fTransport == TR_USB) {
#ifdef HAVE_USB
      unsigned char buffer[64];
      unsigned int base_addr;
      int i, j, ret, n;

      if (type == T_CTRL)
         base_addr = USB_CTRL_OFFSET;
      else if (type == T_STATUS)
         base_addr = USB_STATUS_OFFSET;
      else if (type == T_RAM)
         base_addr = USB_RAM_OFFSET;
      else
         assert(0);             // FIFO not implemented

      if (type != T_RAM) {

         /*---- register access ----*/

         if (size == 2) {
            /* word swapping: first 16 bit sit at uppder address */
            if ((addr % 4) == 0)
               addr = addr + 2;
            else
               addr = addr - 2;
         }

         buffer[0] = USB_CMD_READ;
         buffer[1] = base_addr + addr;
         buffer[2] = size;

         musb_write(fUsbInterface, 2, buffer, 2 + size, USB_TIMEOUT);
         i = musb_read(fUsbInterface, 1, data, size, USB_TIMEOUT);

         if (i != size)
            return 0;

         return size;
      } else {

         /*---- RAM access ----*/

         /* in RAM mode, only the 2048-byte page can be selected */
         buffer[0] = USB_CMD_ADDR;
         buffer[1] = base_addr + (addr >> 11);
         musb_write(fUsbInterface, 2, buffer, 2, USB_TIMEOUT);

         /* receive data in 60-byte packets */
         for (i = 0; i <= (size - 1) / 60; i++) {
            n = size - i * 60;
            if (n > 60)
               n = 60;
            buffer[0] = USB_CMD_READ12;
            buffer[1] = n;
            musb_write(fUsbInterface, 2, buffer, 2, USB_TIMEOUT);

            ret = musb_read(fUsbInterface, 1, buffer, n, USB_TIMEOUT);

            if (ret != n) {
               /* try again */
               ret = musb_read(fUsbInterface, 1, buffer, n, USB_TIMEOUT);
               if (ret != n)
                  return 0;
            }

            for (j = 0; j < ret; j++)
               *((unsigned char *) data + j + i * 60) = buffer[j];
         }

         return size;
      }
#endif                          // HAVE_USB
   } else if (fTransport == TR_USB2) {
#ifdef HAVE_USB
      unsigned char buffer[10];
      unsigned int base_addr;
      int i;

      /* only accept even address and number of bytes */
      assert(addr % 2 == 0);
      assert(size % 2 == 0);

      /* check for maximum size */
      assert(size <= USB2_BUFFER_SIZE - 10);

      if (type == T_CTRL)
         base_addr = USB2_CTRL_OFFSET;
      else if (type == T_STATUS)
         base_addr = USB2_STATUS_OFFSET;
      else if (type == T_FIFO)
         base_addr = USB2_FIFO_OFFSET;
      else if (type == T_RAM)
         base_addr = USB2_RAM_OFFSET;
      else
         base_addr = 0;

      if (type != T_RAM && size == 2) {
         /* word swapping: first 16 bit sit at upper address */
         if ((addr % 4) == 0)
            addr = addr + 2;
         else
            addr = addr - 2;
      }

      addr += base_addr;

      buffer[0] = USB2_CMD_READ;
      buffer[1] = 0;

      buffer[2] = (addr >> 0) & 0xFF;
      buffer[3] = (addr >> 8) & 0xFF;
      buffer[4] = (addr >> 16) & 0xFF;
      buffer[5] = (addr >> 24) & 0xFF;

      buffer[6] = (size >> 0) & 0xFF;
      buffer[7] = (size >> 8) & 0xFF;
      buffer[8] = (size >> 16) & 0xFF;
      buffer[9] = (size >> 24) & 0xFF;

      i = musb_write(fUsbInterface, 4, buffer, 10, USB_TIMEOUT);
      if (i != 10)
         printf("musb_read error %d\n", i);

      i = musb_read(fUsbInterface, 8, data, size, USB_TIMEOUT);
      return i;
#endif                          // HAVE_USB
   }

   return 0;
}

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

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 (fDRSType == 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 (fDRSType == 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);
   } else if (fDRSType == 4) {
      // upper four bits of register must contain last channel to read out
      d = (firstChannel << 4) | lastChannel;
      Write(T_CTRL, REG_CHANNEL_SPAN, &d, 2);

      // set bit pattern for write shift register
      fChannelConfig = 0;
      switch (nConfigChannels) {
      case 1:
         fChannelConfig = 0x01;
         break;
      case 2:
         fChannelConfig = 0x11;
         break;
      case 4:
         fChannelConfig = 0x55;
         break;
      case 8:
         fChannelConfig = 0xFF;
         break;
      default:
         printf("Invalid channel configuration\n");
      }
      d = fChannelConfig | (fDominoMode << 8) | (1 << 9) | (fWSRLoop << 10) | (0xF8 << 8);

      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 */
   if (value < 0)
      value = 0;
   if (value > 2.5)
      value = 2.5;
   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 (fBoardType == 1)
      Read(T_STATUS, buffer, REG_RDAC3, 2);
   else if (fBoardType == 2 || fBoardType == 3 || fBoardType == 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()
{
   // Init FPGA on USB2 board
   InitFPGA();

   // Default values

   // 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;

   if (fBoardType == 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 (fBoardType == 2 || fBoardType == 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 (fBoardType == 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_ROFS_2, 1.05);        // differential input from Lecce splitter

      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
   
   } else if (fBoardType == 5) {
      // DRS4 USB Evaluation Board 1.1 + 2.0

      // set max. domino speed
      SetDAC(fDAC_DSA, 2.5);

      // set readout offset
      fROFS = 1.6;              // differential input range -0.5V ... +0.5V
      SetDAC(fDAC_ROFS_1, fROFS);

      // set common mode offset
      fCommonMode = 0.8;        // 0.8V +- 0.5V inside NMOS range
      SetDAC(fDAC_CMOFS, fCommonMode);

      // calibration voltage
      SetDAC(fDAC_CALP, fCommonMode);
      SetDAC(fDAC_CALN, fCommonMode);

      // OUT- offset
      SetDAC(fDAC_ONOFS, 1.25);

      SetDAC(fDAC_BIAS, 0.70);

   } else if (fBoardType == 6) {
      // DRS4 Mezzanine Board 1.0
      
      // set readout offset
      SetDAC(fDAC_ROFS_1, 1.57);        // differential input range -0.5V ... +0.5V

      // set common mode offset
      fCommonMode = 0.8;        // 0.8V +- 0.5V inside NMOS range
      SetDAC(fDAC_CMOFS, fCommonMode);

      // calibration voltage
      SetDAC(fDAC_CALN, fCommonMode);
      SetDAC(fDAC_CALP, fCommonMode);

      // OUT- offset
      SetDAC(fDAC_ONOFS, 1.25);

      SetDAC(fDAC_BIAS, 0.70);
   }

   /* set default number of channels per chip */
   if (fDRSType == 4)
      SetChannelConfig(0, fNumberOfReadoutChannels - 1, 8);
   else
      SetChannelConfig(0, fNumberOfReadoutChannels - 1, 12);

   // default settings
   fDominoMode = 1;
   fReadoutMode = 1;
   fTriggerEnable1 = 0;
   fTriggerEnable2 = 0;
   fTriggerSource = 0;
   fFrequency = 1;

   SetDominoMode(fDominoMode);
   SetReadoutMode(fReadoutMode);
   EnableTrigger(fTriggerEnable1, fTriggerEnable2);
   SetFrequency(fFrequency);

   /* disalbe 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 (fDRSType == 4) {
      unsigned short d;
      Read(T_CTRL, &d, REG_CONFIG, 2);
      fChannelConfig = d & 0xFF;

      d = fChannelConfig | (fDominoMode << 8) | (1 << 9) | (fWSRLoop << 10) | (0xF8 << 8);
      Write(T_CTRL, REG_CONFIG, &d, 2);
   } else {
      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
   //
   fDominoActive = 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 flag1, int flag2)
{
   // Enable external trigger
   fTriggerEnable1 = flag1;
   fTriggerEnable2 = flag2;
   if (flag1)
      fCtrlBits |= BIT_ENABLE_TRIGGER1;
   else
      fCtrlBits &= ~BIT_ENABLE_TRIGGER1;

   if (flag2)
      fCtrlBits |= BIT_ENABLE_TRIGGER2;
   else
      fCtrlBits &= ~BIT_ENABLE_TRIGGER2;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   return 1;
}

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

int DRSBoard::SetDelayedTrigger(int flag)
{
   // Select delayed trigger from trigger bus
   if (flag)
      fCtrlBits |= BIT_TRIGGER_DELAYED;
   else
      fCtrlBits &= ~BIT_TRIGGER_DELAYED;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   return 1;
}

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

int DRSBoard::SetTriggerLevel(double voltage, bool negative, int delay)
{
   short ticks;
   int delay0;

   fTcalLevel = negative;

   if (negative)
      fCtrlBits |= BIT_NEG_TRIGGER;
   else
      fCtrlBits &= ~BIT_NEG_TRIGGER;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   // Adjust trigger delay to middle of window
   delay0 = (int) ((kNumberOfBins/fFrequency) / 2);
   delay0 -= 33; // internal trigger delay is about 33 ns

   // convert delay in ns into ticks, ~4*580 ps per quad LUT
   ticks = (unsigned short) ((delay0 + delay) / 2.3 + 0.5);
   if (ticks > 511)
      ticks = 511;
   if (ticks < 0)
      ticks = 0;
   Write(T_CTRL, REG_TRG_DELAY, &ticks, 2);

   return SetDAC(fDAC_TLEVEL, voltage/2 + 0.8);
}

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

int DRSBoard::SetTriggerSource(int source)
{
   // Set trigger source 
   // 0=CH1, 1=CH2, 2=CH3, 3=CH4
   if (source & 1)
      fCtrlBits |= BIT_TR_SOURCE1;
   else
      fCtrlBits &= ~BIT_TR_SOURCE1;
   if (source & 2)
      fCtrlBits |= BIT_TR_SOURCE2;
   else
      fCtrlBits &= ~BIT_TR_SOURCE2;

   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::SetTranspMode(int flag)
{
   // Enable/disable transparent mode
   fTranspMode = flag;
   if (flag)
      fCtrlBits |= BIT_TRANSP_MODE;
   else
      fCtrlBits &= ~BIT_TRANSP_MODE;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   return 1;
}

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

int DRSBoard::SetStandbyMode(int flag)
{
   // Enable/disable standby mode
   fTranspMode = flag;
   if (flag)
      fCtrlBits |= BIT_STANDBY_MODE;
   else
      fCtrlBits &= ~BIT_STANDBY_MODE;

   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::IsPLLLocked()
{
   // Get running flag
   unsigned int status;

   Read(T_STATUS, &status, REG_STATUS, 4);
   if (GetBoardType() == 6)
      return ((status >> 1) & 0x0F) == 0x0F;
   return (status & BIT_PLL_LOCKED0) > 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)
{
   if (fDRSType == 4) {

      if (fBoardType == 6) {
         *f = fFrequency;
         return 1;
      }

      unsigned short ticks;

      Read(T_CTRL, &ticks, REG_FREQ_SET, 2);
      ticks += 2;

      /* convert rounded ticks back to frequency */
      if (ticks > 2)
         *f = 1.024 / ticks * fRefClock;
      else
         *f = 0;
   } else {
      // 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 (fDRSType == 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 (fDRSType == 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::ConfigureLMK(double demand, int flag)
{
   unsigned int data[] = { 0x80000100,   // RESET=1
                           0x0003FF00,   // CLKOUT0: EN=1, DIV=FF (=510) MUX=Divided
                           0x00010101,   // CLKOUT1: EN=1, DIV=1 (=2), MUX=Bypassed
                           0x00000102,   // CLKOUT2: DIV=1
                           0x00000103,   // CLKOUT3: DIV=1
                           0x00000104,   // CLKOUT4: DIV=1
                           0x00000105,   // CLKOUT5: DIV=1
                           0x00000106,   // CLKOUT6: DIV=1
                           0x00000107,   // CLKOUT7: DIV=1
                           0x0082000B,   // R11: DIV4=0
                           0x028288AD,   // R13: VCO settings
                           0x0800010E,   // R14: PLL settings
                           0xD400080F }; // R15: VCO Divider=5

   /* calculate dividing ratio */
   int vco_divider = 5;
   double vco = 30*8*vco_divider;
   int divider = (int) ((vco / vco_divider / (demand/2.048) / 2.0) + 0.5);
   static bool configured = false;

   data[1] = 0x00030000 | (divider << 8);

   fFrequency = vco/vco_divider/(divider*2)*2.048;

   for (int i=0 ; i<(int)(sizeof(data)/sizeof(unsigned int)) ; i++) {
      if (configured && i != 1)
         continue; // only update divider after initial configuration
      Write(T_CTRL, REG_LMK_LSB, &data[i], 2);
      Write(T_CTRL, REG_LMK_MSB, ((char *)&data[i])+2, 2);
      Sleep(10);
   }

   configured = true;

   return 1;
}

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

int DRSBoard::SetFrequency(double demand)
{
   // Set domino sampling frequency
   double freq, voltage, delta_voltage;
   unsigned short ticks;
   int i, index, timeout;
   int dominoModeSave = fDominoMode;
   int triggerEnableSave1 = fTriggerEnable1;
   int triggerEnableSave2 = fTriggerEnable2;

   if (fDRSType == 4) {
      /* allowed range is 100 MHz to 6 GHz */
      if (demand > 6 || demand < 0.1)
         return 0;

      if (fBoardType == 6) {
         return ConfigureLMK(demand, 0);
      }

      /* convert frequency in GHz into ticks counted by reference clock */
      if (demand == 0)
         ticks = 0;             // turn off frequency generation
      else
         ticks = static_cast < unsigned short >(1.024 / demand * fRefClock + 0.5);

      ticks -= 2;               // firmware counter need two additional clock cycles
      Write(T_CTRL, REG_FREQ_SET, &ticks, 2);
      ticks += 2;

      /* convert rounded ticks back to frequency */
      if (demand > 0)
         demand = 1.024 / ticks * fRefClock;
      fFrequency = demand;
   } else {                     // fDRSType == 4
      SetDominoMode(1);
      EnableTrigger(0, 0);
      EnableAcal(0, 0);

      fFrequency = demand;

      /* turn automatic adjustment off */
      fCtrlBits &= ~BIT_FREQ_AUTO_ADJ;

      /* disable external trigger */
      fCtrlBits &= ~BIT_ENABLE_TRIGGER1;
      fCtrlBits &= ~BIT_ENABLE_TRIGGER2;

      /* set start pulse length for future DRSBoard_domino_start() */
      if (fDRSType == 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();

      ticks = 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", GetBoardSerialNumber(),
                   index, demand);
            return 0;
         }
      }

      SetDominoMode(dominoModeSave);
      EnableTrigger(triggerEnableSave1, triggerEnableSave2);
   }

   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 (fBoardType == 1)
         Read(T_STATUS, &dac, REG_RDAC3, 2);
      else if (fBoardType == 2 || fBoardType == 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];

   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);
}

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

int DRSBoard::RAMTest(int flag)
{
#define N_WORDS  1024*9
#define N_DWORDS ((N_WORDS)/2)

   int i, j, n;
   unsigned int buffer[N_DWORDS], ret[N_DWORDS];
   unsigned short sbuffer[N_WORDS], sret[N_WORDS];
   time_t now;

   /* integrety test */
   printf("Buffer size: %d (%1.1lfk)\n", static_cast < int >(sizeof(ret)), sizeof(ret) / 1024.0);
   if (flag & 1) {
      if (fTransport == TR_VME) {
         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("RAM test: %d errors\n", n);
      } else if (fTransport == TR_USB2) {
         for (i = 0; i < N_WORDS; i++)
            sbuffer[i] = (rand() & 0xFFFF);     // random 16-bit values
         i = Write(T_RAM, 0, sbuffer, sizeof(sbuffer));

         memset(sret, 0, sizeof(sret));
         i = Read(T_RAM, sret, 0, sizeof(sret));
         if (i != sizeof(sret))
            printf("Read error: %d\n", i);

         for (i = n = 0; i < N_WORDS; i++) {
            if (sbuffer[i] != sret[i]) {
               printf("%03X: %X - %X\n", i, sbuffer[i], sret[i]);
               n++;
            }
         }

         printf("RAM test: %d errors\n", n);
      }
   }

   /* speed test */
   if (flag & 2) {
      /* read continously to determine speed */
      time(&now);
      while (now == time(NULL));
      time(&now);
      i = n = 0;
      if (fTransport == TR_VME) {
         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/sec, %1.2lf MB/sec, %d errors\n", static_cast < int >(i / sizeof(ret)),
                         i / 1024.0 / 1024.0, n);
               else
                  printf("%d read/sec, %1.2lf MB/sec\n", static_cast < int >(i / sizeof(ret)),
                         i / 1024.0 / 1024.0);
               time(&now);
               i = 0;
            }

            if (drs_kbhit())
               break;

         } while (1);

         while (drs_kbhit())
            getch();
      } else if (fTransport == TR_USB2) {
         do {
            memset(ret, 0, sizeof(sret));

            for (j = 0; j < 10; j++) {
               Read(T_RAM, sret, 0, sizeof(sret));
               i += sizeof(sret);
            }

            if (flag & 1) {
               for (j = 0; j < N_WORDS; j++)
                  if (sbuffer[j] != sret[j])
                     n++;
            }

            if (now != time(NULL)) {
               if (flag & 1)
                  printf("%d read/sec, %1.2lf MB/sec, %d errors\n", static_cast < int >(i / sizeof(ret)),
                         i / 1024.0 / 1024.0, n);
               else
                  printf("%d read/sec, %1.2lf MB/sec\n", static_cast < int >(i / sizeof(ret)),
                         i / 1024.0 / 1024.0);
               time(&now);
               i = 0;
            }

            if (drs_kbhit())
               break;

         } while (1);

         while (drs_kbhit())
            getch();
      }
   }
   return 0;
}

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

int DRSBoard::ChipTest()
{
   int i, j, t;
   double freq, old_freq, min, max, mean, rms;
   float  waveform[1024];

   Init();
   SetChannelConfig(0, 8, 8);
   SetDominoMode(1);
   SetReadoutMode(1);
   SetDominoActive(1);
   SetTranspMode(0);
   EnableTrigger(0, 0);
   EnableTcal(1, 0, 0);
   EnableAcal(1, 0);

   /* test 1 GHz */
   SetFrequency(1);
   StartDomino();
   Sleep(100);
   if (!(GetStatusReg() & BIT_PLL_LOCKED0)) {
      puts("PLL did not lock at 1 GHz");
      return 0;
   }

   /* test up to 6 GHz */
   for (freq = 5 ; freq < 6 ; freq += 0.1) {
      SetFrequency(freq);
      Sleep(10);
      if (!(GetStatusReg() & BIT_PLL_LOCKED0)) {
         printf("Max. frequency is %1.1lf GHz\n", old_freq);
         break;
      }
      ReadFrequency(0, &old_freq);
   }

   /* read and check at 0 calibration voltage */
   SetFrequency(5);
   Sleep(10);
   SoftTrigger();
   while (IsBusy());
   TransferWaves(0, 8);

   for (i=0 ; i<8 ; i++) {
      t = GetTriggerCell(0);
      GetWave(0, i, waveform, false, t, false);
      for (j=0 ; j<1024; j++)
         if (waveform[j] < -100 || waveform[j] > 100) {
            if (j<5) {
               /* skip this cells */
            } else {
               printf("Cell error on channel %d, cell %d: %1.1lf mV instead 0 mV\n", i, j, waveform[j]);
               return 0;
            }
         }
   }

   /* read and check at +0.5V calibration voltage */
   EnableAcal(1, 0.5);
   StartDomino();
   SoftTrigger();
   while (IsBusy());
   TransferWaves(0, 8);

   for (i=0 ; i<8 ; i++) {
      t = GetTriggerCell(0);
      GetWave(0, i, waveform, false, t, false);
      for (j=0 ; j<1024; j++)
         if (waveform[j] < 350) {
            if (j<5) {
               /* skip this cell */
            } else {
               printf("Cell error on channel %d, cell %d: %1.1lf mV instead 400 mV\n", i, j, waveform[j]);
               return 0;
            }
         }
   }

   /* read and check at -0.5V calibration voltage */
   EnableAcal(1, -0.5);
   StartDomino();
   Sleep(10);
   SoftTrigger();
   while (IsBusy());
   TransferWaves(0, 8);

   for (i=0 ; i<8 ; i++) {
      t = GetTriggerCell(0);
      GetWave(0, i, waveform, false, t, false);
      for (j=0 ; j<1024; j++)
         if (waveform[j] > -350) {
            if (j<5) {
               /* skip this cell */
            } else {
               printf("Cell error on channel %d, cell %d: %1.1lf mV instead -400mV\n", i, j, waveform[j]);
               return 0;
            }
         }
   }

   /* check clock channel */
   GetWave(0, 8, waveform, false, 0);
   min = max = mean = rms = 0;
   for (j=0 ; j<1024 ; j++) {
      if (waveform[j] > max)
         max = waveform[j];
      if (waveform[j] < min)
         min = waveform[j];
      mean += waveform[j];
   }
   mean /= 1024.0;
   for (j=0 ; j<1024 ; j++)
      rms += (waveform[j] - mean) * (waveform[j] - mean);
   rms = sqrt(rms/1024);

   if (max - min < 400) {
      printf("Error on clock channel amplitude: %1.1lf mV\n", max-min);
      return 0;
   }

   if (rms < 100 || rms > 200) {
      printf("Error on clock channel RMS: %1.1lf mV\n", rms);
      return 0;
   }

   return 1;
}

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

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

   // let DAC settle
   Sleep(100);
}

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

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)
{
   int offset;

   if (fTransport == TR_USB || fTransport == TR_USB2)
      offset = firstChannel * sizeof(short int) * kNumberOfBins;
   else
      offset = 0;               //in VME, always start from zero

   return TransferWaves(fWaveforms + offset, firstChannel, lastChannel);
}

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

int DRSBoard::TransferWaves(unsigned char *p, int firstChannel, int lastChannel)
{
   // Transfer all waveforms at once from VME or USB to location
   int i, n, offset, n_requested;

   if (lastChannel < 0 || lastChannel > kNumberOfChips * kNumberOfChannelsV2) {
      printf("Error: Invalid channel index %d\n", lastChannel);
      return 0;
   }

   if (firstChannel < 0 || firstChannel > kNumberOfChips * kNumberOfChannelsV2) {
      printf("Error: Invalid channel index %d\n", firstChannel);
      return 0;
   }

   if (fTransport == TR_VME) {
      /* in VME, always transfer all waveforms, since channels sit 'next' to each other */
      firstChannel = 0;
      lastChannel = kNumberOfChips * kNumberOfChannelsV2 - 1;
   }

   if (fTransport == TR_USB) {
      /* USB FPGA contains only 16 channels */
      if (lastChannel > 15)
         lastChannel = 15;
   }

   if (fTransport == TR_USB2) {
      /* Evaluation board contains only 9 channels */
      if (lastChannel > 8)
         lastChannel = 8;
   }

   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;

   // read trigger cells
   if (fDRSType == 4) {
      Read(T_STATUS, fTriggerCell, REG_STOP_CELL0, 2); // first chip for the moment...
   }
  
   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;

   /* remap channel */
   if (fBoardType == 1) {
      if (channel < 8)
         channel = 7 - channel;
      else
         channel = 16 - channel;
   } else {
      channel = channel;
   }

   // Read channel
   if (fTransport == TR_USB) {
      offset = kNumberOfBins * 2 * (chipIndex * 16 + channel);
      for (i = 0; i < kNumberOfBins; i++) {
         // 12-bit data
         waveform[i] = ((waveforms[i * 2 + 1 + offset] & 0x0f) << 8) + waveforms[i * 2 + offset];
      }
   } else if (fTransport == TR_USB2) {
      offset = kNumberOfBins * 2 * (chipIndex * 16 + channel);
      for (i = 0; i < kNumberOfBins; i++) {
         // 16-bit data
         waveform[i] = ((waveforms[i * 2 + 1 + offset] & 0xff) << 8) + waveforms[i * 2 + offset];
      }
   } else if (fTransport == TR_VME) {
      offset = (kNumberOfBins * 4) * channel;
      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);
      }
   } else {
      printf("Error: invalid transport %d\n", fTransport);
      return kInvalidTransport;
   }
   return kSuccess;
}

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

int DRSBoard::GetWave(unsigned int chipIndex, unsigned char channel, float *waveform)
{
   return GetWave(chipIndex, channel, waveform, true, fTriggerCell[chipIndex], true, 0, true);
}

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

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

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

int DRSBoard::GetWave(unsigned int chipIndex, unsigned char channel, float *waveform, bool responseCalib,
                      int triggerCell, bool adjustToClock, float threshold, bool offsetCalib)
{
   int ret, i;
   unsigned short waveS[kNumberOfBins];
   ret =
       GetWave(fWaveforms, chipIndex, channel, waveS, responseCalib, triggerCell, adjustToClock, threshold,
               offsetCalib);
   if (responseCalib)
      if (fDRSType == 4) {
         for (i = 0; i < kNumberOfBins; i++) {
            waveform[i] = static_cast < float >(waveS[i] / 65.535); // 16-bit corresponding to 1V
            waveform[i] -= 500; // mid range
         }
      } else {
         for (i = 0; i < kNumberOfBins; i++)
            waveform[i] = static_cast < float >(waveS[i] * GetPrecision());
      }
   else {
      for (i = 0; i < kNumberOfBins; i++) {
         if (fBoardType == 4) {
            waveform[i] = static_cast < float >(waveS[i] / 65.535); // 16-bit corresponding to 1V
            waveform[i] -= 570; // guess some offset
         } else if (fBoardType == 5 || fBoardType == 6) {
            waveform[i] = static_cast < float >(waveS[i] / 65.535); // 16-bit corresponding to 1V
            waveform[i] -= 500; // mid range
         } else
            waveform[i] = static_cast < float >(waveS[i]);
      }
   }
   return ret;
}

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

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

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

   return ret;
}

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

int DRSBoard::GetWave(unsigned char *waveforms, unsigned int chipIndex, unsigned char channel,
                      unsigned short *waveform, bool responseCalib, int triggerCell, bool adjustToClock,
                      float threshold, bool offsetCalib)
{
   if (!fWaveTransferred[chipIndex * kNumberOfChannelsV2 + 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, offsetCalib);
}

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

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::CalibrateWaveform(unsigned int chipIndex, unsigned char channel, unsigned short *adcWaveform,
                                unsigned short *waveform, bool responseCalib,
                                int triggerCell, bool adjustToClock, float threshold, bool offsetCalib)
{
   int j;
   double value;

   // calibrate waveform
   if (responseCalib) {
      if (GetDRSType() == 4) {
         for (j = 0 ; j < kNumberOfBins; j++) {
            value = adcWaveform[j] - (fCellOffset[channel][(j + triggerCell) % kNumberOfBins] - 32768);
            value = (value - 32768) / fCellGain[channel][(j + triggerCell) % kNumberOfBins] + 32768;
            if (value > 65535)
               value = 65535;
            if (value < 0)
               value = 0;

            if (adjustToClock)          
               waveform[(j + triggerCell) % kNumberOfBins] = (unsigned short) (value + 0.5);
            else
               waveform[j] = (unsigned short) (value + 0.5); 
         }
      } else {
         if (!fResponseCalibration->
             Calibrate(chipIndex, channel % 10, adcWaveform, waveform, triggerCell, threshold, offsetCalib))
            return kZeroSuppression;       // return immediately if below threshold
      }
   } else {
      for (j = 0; j < kNumberOfBins; j++) {
         if (adjustToClock) {
            // rotate waveform such that waveform[0] corresponds to bin #0 on the chip
            waveform[j] = adcWaveform[(kNumberOfBins-triggerCell+j) % kNumberOfBins];
         } else {
            waveform[j] = adcWaveform[j];
         }
      }
   }
   // fix bad cells for single turn mode
   if (GetDRSType() == 2) {
      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)
{
   if (fDRSType == 4) {
      unsigned short d;

      Read(T_STATUS, &d, REG_STOP_CELL0 + chipIndex * 4, 2);
      return d;
   } else
      return GetTriggerCell(fWaveforms, chipIndex);
}

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

int DRSBoard::GetTriggerCell(unsigned char *waveforms, unsigned int chipIndex)
{
   int j, triggerCell;
   bool calib;
   unsigned short baseLevel = 1000;
   unsigned short triggerChannel[1024];
   if (!fWaveTransferred[chipIndex * kNumberOfChannelsV2 + 8])
      return -1;

   GetADCWave(waveforms, chipIndex, 8, triggerChannel);
   calib = fResponseCalibration->SubtractADCOffset(chipIndex, 8, triggerChannel, triggerChannel, baseLevel);

   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 (fDRSType == 3) {
            if (triggerChannel[j] <= 2000 && triggerChannel[(j + 1) % kNumberOfBins] > 2000) {
               triggerCell = j;
               break;
            }
         } else {
            if (triggerChannel[j] >= 2000 && triggerChannel[(j + 1) % kNumberOfBins] < 2000) {
               triggerCell = j;
               break;
            }
         }
      }
   }
   if (triggerCell == -1) {
      return kInvalidTriggerSignal;
   }
   fTriggerCell[0] = triggerCell;
   return triggerCell;
}

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

int DRSBoard::GetStopCell(unsigned int chipIndex)
{
   unsigned short d;

   switch (chipIndex) {
   case 0:
      Read(T_STATUS, &d, REG_STOP_CELL0, 2);
      break;
   case 1:
      Read(T_STATUS, &d, REG_STOP_CELL1, 2);
      break;
   case 2:
      Read(T_STATUS, &d, REG_STOP_CELL2, 2);
      break;
   case 3:
      Read(T_STATUS, &d, REG_STOP_CELL3, 2);
      break;
   }

   return d;
}

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

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 BOARD_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 (fBoardType == 1)
         Read(T_STATUS, &data, REG_RDAC3, 2);
      else if (fBoardType == 2 || fBoardType == 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);
}

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

int DRSBoard::TestShift(int n)
{
   // Test shift register
   unsigned char buffer[3];

   memset(buffer, 0, sizeof(buffer));

#if 0
   buffer[0] = CMD_TESTSHIFT;
   buffer[1] = n;

   status = msend_usb(buffer, 2);
   if (status != 2)
      return status;

   status = mrecv_usb(buffer, sizeof(buffer));
   if (status != 1)
      return status;
#endif

   if (buffer[0] == 1)
      printf("Shift register %c works correctly\n", 'A' + n);
   else if (buffer[0] == 2)
      printf("SROUT%c does hot go high after reset\n", 'A' + n);
   else if (buffer[0] == 3)
      printf("SROUT%c does hot go low after 1024 clocks\n", 'A' + n);

   return 1;
}

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

unsigned int DRSBoard::GetCtrlReg()
{
   unsigned int status;

   Read(T_CTRL, &status, REG_CTRL, 4);
   return status;
}

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

unsigned short DRSBoard::GetConfigReg()
{
   unsigned short status;

   Read(T_CTRL, &status, REG_CONFIG, 2);
   return status;
}

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

unsigned int DRSBoard::GetStatusReg()
{
   unsigned int status;

   Read(T_STATUS, &status, REG_STATUS, 4);
   return status;
}

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

int DRSBoard::EnableTcal(int flag, int level, int source)
{
   fTcalMode = flag;
   fTcalLevel = level;
   fTcalSource = source;

   // Enable clock channel
   if (flag)
      fCtrlBits |= BIT_TCAL_EN;
   else
      fCtrlBits &= ~BIT_TCAL_EN;

   // Set output level, needed for gain calibration
   if (fDRSType == 4) {
      if (level)
         fCtrlBits |= BIT_NEG_TRIGGER;
      else
         fCtrlBits &= ~BIT_NEG_TRIGGER;
   }

   // Select clock source synchronous or asynchronous (2nd quartz)
   if (source) 
      fCtrlBits |= BIT_TCAL_SOURCE;
   else
      fCtrlBits &= ~BIT_TCAL_SOURCE;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   return 1;
}

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

int DRSBoard::SetRefclk(int source)
{
   // Select clock source to internal FPGA (0) or external P2 (1)
   if (source) 
      fCtrlBits |= BIT_REFCLK_SOURCE;
   else
      fCtrlBits &= ~BIT_REFCLK_SOURCE;

   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);

   return 1;
}

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

int DRSBoard::EnableAcal(int mode, double voltage)
{
   double t1, t2;

   fAcalMode = mode;
   fAcalVolt = voltage;

   if (mode == 0) {
      /* turn calibration off */
      SetCalibTiming(0, 0);
      if (fBoardType == 5 || fBoardType == 6) {
         /* turn voltages off (50 Ohm analog switch!) */
         SetDAC(fDAC_CALP, 0);
         SetDAC(fDAC_CALN, 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 (fDRSType == 2) {
      d = t_cal | (t_enable << 8);
      Write(T_CTRL, REG_CALIB_TIMING, &d, 2);
   }

   if (fDRSType == 3) {
      d = t_cal;
      Write(T_CTRL, REG_CALIB_TIMING, &d, 2);
   }

   return 1;
}

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

int DRSBoard::SetCalibVoltage(double value)
{
   // Set Calibration Voltage
   if (fBoardType == 5 || fBoardType == 6) {
      value = value * (1+fFrequency/65); // rough correction factor for input current
      SetDAC(fDAC_CALP, fCommonMode + value / 2);
      SetDAC(fDAC_CALN, fCommonMode - value / 2);
   } else
      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::SetBoardSerialNumber(unsigned short serialNumber)
{
   unsigned char buf[32768];

   /* merge serial number into eeprom page #0 */
   ReadEEPROM(0, buf, sizeof(buf));
   buf[0] = serialNumber & 0xFF;
   buf[1] = serialNumber >> 8;
   WriteEEPROM(0, buf, sizeof(buf));

   /* erase DPRAM */
   memset(buf, 0, sizeof(buf));
   Write(T_RAM, 0, buf, sizeof(buf));

   /* read back EEPROM */
   ReadEEPROM(0, buf, sizeof(buf));

   /* check if correctly set */
   if (((buf[1] << 8) | buf[0]) != serialNumber)
      return 0;

   fBoardSerialNumber = serialNumber;
   return 1;
}

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

int DRSBoard::ReadEEPROM(unsigned short page, void *buffer, int size)
{
   unsigned long status;
   int i;

   // write eeprom page number
   Write(T_CTRL, REG_EEPROM_PAGE, &page, 2);

   // execute eeprom read
   fCtrlBits |= BIT_EEPROM_READ_TRIG;
   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
   fCtrlBits &= ~BIT_EEPROM_READ_TRIG;

   // poll on eeprom_busy flag
   for (i=0 ; i<100 ; i++) {
      Read(T_STATUS, &status, REG_STATUS, 4);
      if ((status & BIT_EEPROM_BUSY) == 0)
         break;
      Sleep(10);
   }

   return Read(T_RAM, buffer, 0, size);
}

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

int DRSBoard::WriteEEPROM(unsigned short page, void *buffer, int size)
{
   unsigned long status;
   int i;

   // write eeprom page number
   Write(T_CTRL, REG_EEPROM_PAGE, &page, 2);

   // write eeprom page to RAM
   Write(T_RAM, 0, buffer, size);

   // execute eeprom write
   fCtrlBits |= BIT_EEPROM_WRITE_TRIG;
   Write(T_CTRL, REG_CTRL, &fCtrlBits, 4);
   fCtrlBits &= ~BIT_EEPROM_WRITE_TRIG;

   // poll on eeprom_busy flag
   for (i=0 ; i<100 ; i++) {
      Read(T_STATUS, &status, REG_STATUS, 4);
      if ((status & BIT_EEPROM_BUSY) == 0)
         break;
      Sleep(10);
   }

   return 1;
}

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

int DRSBoard::GetTime(unsigned int chipIndex, float *time, bool tcalibrated, bool rotated)
{
   int i;

   /* for DRS2, please use function below */
   if (fDRSType < 4)
      return GetTime(chipIndex, fFrequency, time, tcalibrated, rotated);

   if (!fTimingCalibrationValid || fFrequency < 4.5 || !tcalibrated) {
      for (i = 0; i < kNumberOfBins; i++)
         time[i] = static_cast < float >(i / fFrequency);
      return 1;
   }

   double t0 = fCellT[chipIndex][fTriggerCell[chipIndex]];

   for (i=0 ; i<kNumberOfBins ; i++) {
      if (rotated)
         time[i] = (float) (fCellT[chipIndex][(i+fTriggerCell[chipIndex]) % kNumberOfBins] - t0);
      else
         time[i] = (float) (fCellT[chipIndex][i] - t0);
      if (time[i] < 0)
         time[i] += (float) (kNumberOfBins / fFrequency);
   }
   return 1;
}

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

int DRSBoard::GetTime(unsigned int chipIndex, double freqGHz, float *time, bool tcalibrated, bool rotated)
{
   /* for DRS4, use function above */
   if (fDRSType == 4)
      return GetTime(chipIndex, time, tcalibrated, rotated);

   int i, irot;
   DRSBoard::TimeData * init;
   DRSBoard::TimeData::FrequencyData * freq;
   int frequencyMHz = (int)(freqGHz*1000);

   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 = (fTriggerCell[chipIndex] + i) % kNumberOfBins;
      if (fTriggerCell[chipIndex] + i < kNumberOfBins)
         time[i] = static_cast < float >((freq->fBin[irot] - freq->fBin[fTriggerCell[chipIndex]]) / fFrequency);
      else
      time[i] =
          static_cast <
          float
          >((freq->fBin[irot] - freq->fBin[fTriggerCell[chipIndex]] + 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, fBoardSerialNumber,
              fBoardSerialNumber, 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", GetBoardSerialNumber());
   }

   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;
}

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

int DRSBoard::AverageWaveforms(DRSCallback *pcb, int prog1, int prog2, float awf[kNumberOfCalibChannelsV4][1024], int n)
{
   int i, j, k, tc, prog, old_prog = 0;
   float wf[kNumberOfCalibChannelsV4][1024];

   if (pcb != NULL)
      pcb->Progress(prog1);

   for (j=0 ; j<kNumberOfCalibChannelsV4 ; j++)
      for (k=0 ; k<1024 ; k++)
         awf[j][k] = 0;

   for (i=0 ; i<n; i++) {
      for (j=0 ; j<10 ; j++)
         IsBusy();
      SoftTrigger();
      while (IsBusy());
      StartDomino();
      TransferWaves();
      tc = GetTriggerCell(0);

      for (j=0 ; j<kNumberOfCalibChannelsV4 ; j++) {
         GetWave(0, j, wf[j], false, tc, 1);

         for (k=0 ; k<1024 ; k++)
            awf[j][k] += wf[j][k];
      }

      prog = (int)(((double)i/n)*(prog2-prog1)+prog1);
      if (prog > old_prog) {
         old_prog = prog;
         if (pcb != NULL)
            pcb->Progress(prog);
      }
   }

   for (j=0 ; j<kNumberOfCalibChannelsV4 ; j++)
      for (k=0 ; k<1024 ; k++)
         awf[j][k] /= n;

   if (pcb != NULL)
      pcb->Progress(prog2);

   return 1;
}

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

int DRSBoard::CalibrateVolt(DRSCallback *pcb)
{
int    i, j;
float  wf1[kNumberOfCalibChannelsV4][1024], wf2[kNumberOfCalibChannelsV4][1024];
double f, g1, g2;
unsigned short buf[1024*kNumberOfCalibChannelsV4*2];
   
   f = fFrequency;

   Init();
   SetFrequency(f);
   SetDominoMode(1);
   SetDominoActive(1);
   SetReadoutMode(1);
   EnableTcal(0, 0, 0);
   StartDomino();

   /* measure offset */
   EnableAcal(1, 0);
   Sleep(100);
   AverageWaveforms(pcb, 0, 50, wf1, 100);

   /* convert offsets to 16-bit values */
   memset(fCellOffset, 0, sizeof(fCellOffset));
   for (i=0 ; i<8 ; i++)
      for (j=0 ; j<1024; j++) {
         fCellOffset[i][j] = (unsigned short) ((wf1[i][j]+500)/1000.0*65536 + 0.5);
      }

   /* measure gain at +400 mV */
   EnableAcal(1, 0.4);
   Sleep(100);
   AverageWaveforms(pcb, 50, 75, wf1, 50);

   /* measure gain at -400 mV */
   EnableAcal(1, -0.4);
   EnableTcal(0, 1, 0);
   Sleep(100);
   AverageWaveforms(pcb, 75, 100, wf2, 50);

   for (i=0 ; i<8 ; i++)
      for (j=0 ; j<1024; j++) {

         /* calculate gain */
         g1 = (wf1[i][j] - (fCellOffset[i][j] - 32768)/65536.0*1000) / 400.0;
         g2 = ((fCellOffset[i][j] - 32768)/65536.0*1000 - wf2[i][j]) / 400.0;

         fCellGain[i][j] = (g1 + g2) / 2;

         /* convert 0.8 .. 1.1 to 0..65535 */
         if (fCellGain[i][j] < 0.8 || fCellGain[i][j] > 1.1) {
            printf("Gain of %1.3lf for channel %d, cell %d out of range 0.8 ... 1.1\n",
               fCellGain[i][j], i, j);
         }
      }

   for (j=0 ; j<1024; j++) {

      /* calculate offset and gain for timing channel */
      fCellOffset[8][j] = (unsigned short) (((wf1[8][j]+wf2[8][j])/2+500)/1000.0*65536 + 0.5);
      fCellGain[8][j]   = (wf2[8][j] - wf1[8][j]) / 650.0;

      /* convert 0.8 .. 1.1 to 0..65535 */
      if (fCellGain[8][j] < 0.8 || fCellGain[8][j] > 1.1) {
         printf("Gain of %1.3lf for channel %d, cell %d out of range 0.8 ... 1.1\n",
            fCellGain[8][j], i, j);
      }
   }

   /* write calibration CH0-CH7 to EEPROM page 1 */
   ReadEEPROM(1, buf, 1024*32);
   for (i=0 ; i<8 ; i++)
      for (j=0 ; j<1024; j++) {
         buf[(i*1024+j)*2]   = fCellOffset[i][j];
         buf[(i*1024+j)*2+1] = (unsigned short) ((fCellGain[i][j] - 0.8) / 0.3 * 65535);
      }
   WriteEEPROM(1, buf, 1024*32);

   /* write calibration CH8 to EEPROM page 2 */
   ReadEEPROM(2, buf, 1024*4);
   for (j=0 ; j<1024; j++) {
      buf[j*2]   = fCellOffset[8][j];
      buf[j*2+1] = (unsigned short) ((fCellGain[8][j] - 0.8) / 0.3 * 65535);
   }
   WriteEEPROM(2, buf, 1024*4);

   fCellCalibrationValid = true;

   /* remove calibration voltage */
   EnableAcal(0, 0);
   EnableTcal(0, 0, 0);

   return 1;
}

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

int DRSBoard::CalibrateTiming(DRSCallback *pcb)
{
int    index, tc, i, i1, i2, j, nzx, zeroXing[1000], n_iter, clkon;
float  wf[1024];
double old[1024];
double damping, fSamp, fCalib, zeroLevel, t1, t2, dt, corr, inv_corr;
unsigned short buf[1024*2];
   
   fSamp   = fFrequency;
   fCalib  = 264;
   n_iter  = 5000;
   damping = 0.01;

   clkon = (GetCtrlReg() & BIT_TCAL_EN) > 0;

   Init();
   SetFrequency(fSamp);
   SetDominoMode(1);
   SetDominoActive(1);
   SetReadoutMode(1);
   EnableTcal(1, 0, 1);
   StartDomino();

   /* initialize time array */
   for (i=0 ; i<1024 ; i++)
      fCellT[0][i] = (float)1/fSamp*i; // [ns]

   for (index = 0 ; index < n_iter ; index++) {
      if (index % 10 == 0)
         pcb->Progress(100*index/n_iter);

      SoftTrigger();
      while (IsBusy());
      StartDomino();
      TransferWaves();
      tc = GetTriggerCell(0);

      GetWave(0, 8, wf, true, tc, true);

      /* calculate zero level */
      for (i=0,zeroLevel=0 ; i<1024 ; i++)
         zeroLevel += wf[i];
      zeroLevel /= 1024;

      /* correct for zero common mode */
      for (i=0 ; i<1024 ; i++)
         wf[i] -= (float)zeroLevel;

      /* find falling edge zero crossing with wrap-around */
      for (i=tc+3,nzx=0 ; i<tc+1023 && nzx < (int)(sizeof(zeroXing)/sizeof(int)) ; i++) {
         if (wf[(i+1) % 1024] < 0 && wf[i % 1024] > 0)
            zeroXing[nzx++] = i;
      }

      if (nzx < 40)
         return 0;

      for (i=0 ; i<nzx-1 ; i++) {
         i1 = zeroXing[i] % 1024;
         t1 = fCellT[0][i1] + (fCellT[0][i1+1] - fCellT[0][i1])*(1/(1-wf[i1+1]/wf[i1]));
         i2 = zeroXing[i+1] % 1024;
         t2 = fCellT[0][i2] + (fCellT[0][i2+1] - fCellT[0][i2])*(1/(1-wf[i2+1]/wf[i2]));

         /* wrap-around ? */
         if (t2 - t1 < 0)
            t2 += 1/fSamp*1024;

         /* calculate correction to nominal period in ns */
         corr = 1/fCalib*1000 - (t2 - t1);

         /* skip very large corrections (noise?) */
         if (fabs(corr) > 0.2)
            continue;

         /* apply damping factor */
         corr *= damping;

         /* calculate inverse correction */
         inv_corr = -corr;

         /* apply from (i1+1)+1 to i2 inclusive */
         i1 = zeroXing[i]+2;
         i2 = zeroXing[i+1];

         /* distribute correciton equally into bins inside the region ... */
         corr = corr / (i2-i1+1);

         /* ... and inverse correction into the outside bins */
         inv_corr = inv_corr / (1024 - (i2-i1+1));

         i1 = i1 % 1024;
         i2 = i2 % 1024;

         memcpy(old, fCellT, sizeof(old));

         for (j=0,t1=0 ; j<1024 ; j++) {
            if (j < 1023)
               dt = fCellT[0][j+1] - fCellT[0][j];
            else
               dt = 1/fSamp*1024 - fCellT[0][j];
            if ((i2 > i1 && (j >= i1 && j<= i2)) ||
                (i2 < i1 && (j >= i1 || j<= i2)) )
               dt += corr;
            else
               dt += inv_corr;

            fCellT[0][j] = t1;
            t1 += dt;
         }
      }
   }

   pcb->Progress(100);

   /* write timing calibration to EEPROM page 0 */
   ReadEEPROM(0, buf, 1024*4);
   for (i=0,t1=0 ; i<1024; i++) {
      t2 = fCellT[0][i] - t1;
      t2 = (unsigned short) (t2 * 10000 + 0.5);
      buf[i*2+1] = (unsigned short) t2;
      t1 += t2 / 10000.0;
   }
   WriteEEPROM(0, buf, 1024*4);

   fTimingCalibrationValid = true;

   /*
   {
      FILE *f;

      f = fopen("cellt.csv", "wt");
      for (i=0 ; i<1024 ; i++)
         fprintf(f, "%5.3lf\n", fCellT[0][i]);
      fclose(f);
   }
   */

   /* remove calibration voltage */
   EnableAcal(0, 0);
   EnableTcal(clkon, 0, 0);

   return 1;
}

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

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->GetDRSType() == 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->GetBoardSerialNumber());
   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->GetBoardSerialNumber(),
           fBoard->GetBoardSerialNumber(), 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 < kNumberOfCalibChannelsV3; 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->GetBoardSerialNumber());
   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->GetBoardSerialNumber(),
           fBoard->GetBoardSerialNumber(), 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 < kNumberOfCalibChannelsV4; 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->GetDRSType() == 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 < kNumberOfCalibChannelsV3; 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 < kNumberOfCalibChannelsV3; ii++) {
         fBoard->GetADCWave(chipNumber, ii, fWaveFormMode2[ii][j]);
      }
      CalibrationStart(voltage);
   }
   // Average Sample Points
   for (ii = 0; ii < kNumberOfCalibChannelsV3; 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++) {
      if (fBoard->Debug()) {
         printf("%02d:%02d\r", fNumberOfPoints - fCurrentPoint, fNumberOfSamples - i);
         fflush(stdout);
      }


      fBoard->SoftTrigger();
      while (fBoard->IsBusy());
      fBoard->StartDomino();
      Sleep(50);
      fBoard->TransferWaves();
      for (j = 0; j < kNumberOfCalibChannelsV4; j++) {
         fBoard->GetADCWave(chipNumber, j, fWaveFormMode3[j][i]);
      }
   }

   // Calculate averages
   for (i = 0; i < kNumberOfCalibChannelsV4; 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->GetDRSType() == 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 < kNumberOfCalibChannelsV3; 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 == kNumberOfCalibChannelsV3) {
      if (fShowStatistics) {
         for (i = 0; i < fNumberOfPoints; i++) {
            average = 0;
            averageError = 0;
            averageExt = 0;
            averageErrorExt = 0;
            for (j = 0; j < kNumberOfCalibChannelsV3 * 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 /= kNumberOfCalibChannelsV3 * kNumberOfBins;
            averageError =
                sqrt((averageError -
                      average * average / kNumberOfCalibChannelsV3 * kNumberOfBins) /
                     (kNumberOfCalibChannelsV3 * kNumberOfBins - 1));
            averageExt /= kNumberOfCalibChannelsV3 * kNumberOfBins;
            averageErrorExt =
                sqrt((averageErrorExt -
                      averageExt * averageExt / kNumberOfCalibChannelsV3 * kNumberOfBins) /
                     (kNumberOfCalibChannelsV3 * 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];
   static int error;

   CalibrationData *data = fCalibrationData[chipNumber];
   CalibrationData::CalibrationDataChannel * chn = data->fChannel[fCurrentFitChannel];

   if (fCurrentFitBin == 0 && fCurrentFitChannel == 0) {
      error = 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]);
   // exclude first two points (sometimes are on limit of FADC)
   DRSBoard::LinearRegression(fXXFit + 2, fYYFit + 2, fNumberOfPoints - 2, &par[1], &par[0]);

   chn->fOffset[fCurrentFitBin] = static_cast < unsigned short >(par[0] + 0.5);
   chn->fGain[fCurrentFitBin] = static_cast < unsigned short >(par[1] + 0.5);

   // Remember min/max of gain
   if (fCurrentFitBin == 0 && fCurrentFitChannel == 0)
      fGainMin = fGainMax = chn->fGain[0];
   if (chn->fGain[fCurrentFitBin] < fGainMin)
      fGainMin = chn->fGain[fCurrentFitBin];
   if (chn->fGain[fCurrentFitBin] > fGainMax)
      fGainMax = chn->fGain[fCurrentFitBin];

   // abort if outside normal region
   if (chn->fGain[fCurrentFitBin] / 4096.0 < 0.8 || chn->fGain[fCurrentFitBin] / 4096.0 > 1) {
      error++;

      if (error < 20)
         printf("Gain=%1.3lf for bin %d on channel %d on chip %d outside valid region\n",
                chn->fGain[fCurrentFitBin] / 4096.0, fCurrentFitBin, fCurrentFitChannel, chipNumber);
   }

   if (fCurrentFitChannel == 1 && fCurrentFitBin == 10) {
      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 == kNumberOfCalibChannelsV4) {

      if (fBoard->Debug()) {
         printf("Gain min=%1.3lf max=%1.3lf\n", fGainMin / 4096.0, fGainMax / 4096.0);
         fflush(stdout);
      }
      // allow up to three bad bins
      if (error > 3) {
         printf("Aborting calibration!\n");
         return true;
      }

      fFitted = true;
      fOffset = false;
      fCalibrationData[chipNumber]->fRead = true;
      fCalibrationData[chipNumber]->fHasOffsetCalibration = false;
      return true;
   }

   return false;
}

unsigned int millitime()
{
#ifdef _MSC_VER

   return (int) GetTickCount();

#else
   struct timeval tv;

   gettimeofday(&tv, NULL);

   return tv.tv_sec * 1000 + tv.tv_usec / 1000;
#endif
   return 0;
}

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

bool ResponseCalibration::OffsetCalibration(int chipNumber)
{
   if (!fFitted || fOffset)
      return true;
   if (fBoard->GetDRSType() == 3)
      return OffsetCalibrationV4(chipNumber);
   else
      return OffsetCalibrationV3(chipNumber);
}

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

bool ResponseCalibration::OffsetCalibrationV3(int chipNumber)
{
   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(0, 0.0);
   }
   // Loop Over Number Of Samples For Statistics
   t1 = millitime();
   fBoard->SoftTrigger();
   while (fBoard->IsBusy()) {
   }
   fBoard->TransferWaves();
   for (ii = 0; ii < kNumberOfCalibChannelsV3; ii++) {
      fBoard->GetADCWave(chipNumber, ii, fWaveFormOffsetADC[ii][fCurrentSample]);
      fBoard->CalibrateWaveform(chipNumber, ii, fWaveFormOffsetADC[ii][fCurrentSample],
                                fWaveFormOffset[ii][fCurrentSample], true, false, false, 0, true);
   }
   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 < kNumberOfCalibChannelsV3; 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;
}

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

bool ResponseCalibration::OffsetCalibrationV4(int chipNumber)
{
   int k, ii, j;
   float mean, error;
   CalibrationData *data = fCalibrationData[chipNumber];
   CalibrationData::CalibrationDataChannel * chn;

   /* switch DRS to input, hope that no real signal occurs */
   if (fCurrentSample == 0) {
      data->fHasOffsetCalibration = false;
      fBoard->SetCalibVoltage(0.0);
      fBoard->EnableAcal(0, 0.0);
      /* one dummy trigger for unknown reasons */
      fBoard->SoftTrigger();
      while (fBoard->IsBusy());
      fBoard->StartDomino();
      Sleep(50);
   }
   // Loop Over Number Of Samples For Statistics
   fBoard->SoftTrigger();
   while (fBoard->IsBusy());
   fBoard->TransferWaves();
   for (ii = 0; ii < kNumberOfCalibChannelsV4; ii++)
      fBoard->GetADCWave(chipNumber, ii, fWaveFormOffsetADC[ii][fCurrentSample]);

   fBoard->StartDomino();
   Sleep(50);
   fCurrentSample++;

   if (fBoard->Debug()) {
      printf("%02d\r", fNumberOfSamples - fCurrentSample);
      fflush(stdout);
   }

   if (fCurrentSample == fNumberOfSamples) {
      // Average Sample Points
      float *sample = new float[fNumberOfSamples];
      for (ii = 0; ii < kNumberOfCalibChannelsV3; ii++) {
         chn = data->fChannel[ii];
         for (k = 0; k < kNumberOfBins; k++) {
            for (j = 0; j < fNumberOfSamples; j++)
               sample[j] = static_cast < float >(fWaveFormOffsetADC[ii][j][k]);
            Average(1, sample, fNumberOfSamples, mean, error, 2);
            chn->fOffset[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 < kNumberOfCalibChannelsV3; ii++) {
      for (j = 0; j < kNumberOfBins; j++) {
         fResponseX[ii][j] = new float[fNumberOfPoints + fNumberOfPointsLowVolt];
      }
   }
   fResponseY = new float[fNumberOfPoints + fNumberOfPointsLowVolt];
   for (ii = 0; ii < kNumberOfCalibChannelsV3; ii++) {
      fWaveFormMode3[ii] = new unsigned short *[fNumberOfSamples];
      fWaveFormMode2[ii] = new unsigned short *[fNumberOfSamples];
      fWaveFormOffset[ii] = new unsigned 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 unsigned 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[kNumberOfCalibChannelsV3 * kNumberOfBins];
      fStatisticsApproxExt[i] = new float[kNumberOfCalibChannelsV3 * 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 < kNumberOfCalibChannelsV3; ii++) {
      for (j = 0; j < kNumberOfBins; j++) {
         delete fResponseX[ii][j];
      }
   }
   delete fResponseY;
   for (ii = 0; ii < kNumberOfCalibChannelsV3; 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,
                                    unsigned short *uWaveform, int triggerCell, float threshold, bool offsetCalib)
{
   int i;
   unsigned int NumberOfCalibChannels;
   int hasOffset;
   bool aboveThreshold;
   float wave, v;
   int j, irot;

   CalibrationData *data = fCalibrationData[chipIndex];
   CalibrationData::CalibrationDataChannel * chn;

   if (fBoard->GetDRSType() == 3)
      NumberOfCalibChannels = kNumberOfCalibChannelsV4;
   else
      NumberOfCalibChannels = kNumberOfCalibChannelsV3;

   if (channel >= NumberOfCalibChannels || data == NULL) {
      for (i = 0; i < kNumberOfBins; i++) {
         irot = i;
         if (triggerCell > -1)
            irot = (triggerCell + i) % kNumberOfBins;

         uWaveform[i] = adcWaveform[irot];
      }
      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

   short offset;

   // Calibrate
   for (i = 0; i < kNumberOfBins; i++) {
      if (fBoard->GetDRSType() != 3) {
         irot = i;
         if (triggerCell > -1)
            irot = (triggerCell + i) % kNumberOfBins;
         offset = offsetCalib ? chn->fOffset[irot] : 0;
         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] -= offset;
      } else {
         irot = i;
         if (triggerCell > -1)
            irot = (triggerCell + i) % kNumberOfBins;
#if 0                           /* not enabled yet for DRS3 */
         offset = offsetCalib ? chn->fOffset[irot] : 0;
#else
         offset = chn->fOffset[irot];
#endif
         v = static_cast < float >(adcWaveform[irot] - offset) / chn->fGain[irot];
         uWaveform[i] = static_cast < short >(v * 1000 / GetPrecision() + 0.5);
      }

      // 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;
   unsigned int NumberOfCalibChannels;
   CalibrationData *data = fCalibrationData[chipIndex];
   CalibrationData::CalibrationDataChannel * chn;

   if (fBoard->GetDRSType() == 3)
      NumberOfCalibChannels = kNumberOfCalibChannelsV4;
   else
      NumberOfCalibChannels = kNumberOfCalibChannelsV3;

   if (channel >= NumberOfCalibChannels || 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->GetDRSType() == 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->GetBoardSerialNumber(), fBoard->GetBoardSerialNumber(), 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->GetBoardSerialNumber());
      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->GetDRSType() != 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 < kNumberOfCalibChannelsV3; k++) {
      chn = data->fChannel[k];
      for (l = 0; l < kNumberOfBins; l++) {
         if (fBoard->GetDRSType() != 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->GetDRSType() == 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->GetDRSType() != 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->GetBoardSerialNumber(), fBoard->GetBoardSerialNumber(), 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->GetBoardSerialNumber());
      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 < kNumberOfCalibChannelsV4; 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->GetDRSType() == 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)) {
#ifdef _MSC_VER
      return mkdir(path);
#else
      return mkdir(path, 0711);
#endif                          // R__UNIX
   }
   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 < kNumberOfCalibChannelsV3; 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 < kNumberOfCalibChannelsV3; 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;
   }
}