/*===========================================================================*/
/*                                gainResults.C                              */
/* ------------------------------------------------------------------------- */
/* ISDCRoot script to extract line fit results to calibrate SPI gains.       */
/* ------------------------------------------------------------------------- */
/* Author  : Jurgen KNODLSEDER (CESR) (C) (all rights reserved)              */
/* Date    : 25-Sep-2003                                                     */
/* Version : 1.7.0                                                           */
/* ------------------------------------------------------------------------- */
/* History :                                                                 */
/* 1.0.0  11-Nov-2002  first version                                         */
/* 1.1.0  10-Mar-2003  read also OBT and ONTIME information                  */
/* 1.2.0  10-Mar-2003  read line characteristics from SPI.-LINE-SCT file     */
/* 1.3.0  17-Mar-2003  modify line list                                      */
/* 1.4.0  24-Mar-2003  modify line list (version 0003)                       */
/* 1.5.0  27-Mar-2003  modify line list (version 0004)                       */
/* 1.6.0  17-Sep-2003  modify line list (version 0005)                       */
/* 1.7.0  25-Sep-2003  don't shift 23.4 and 198.4 keV lines anymore          */
/*===========================================================================*/
#define ISDC_OK  0
#define DAL_SAVE 0
void gainResults(int argc = -2, char** argv = 0) {

    // Set path and file names
    char sctNameSE[] = "/users/jurgen/cesr_analysis/cal/spi_cal_lines_se_0005.fits";
    char sctNameME[] = "/users/jurgen/cesr_analysis/cal/spi_cal_lines_me_0005.fits";
    char srtNameSE[] = "spi_cal_se_results.fits";
    char srtNameME[] = "spi_cal_me_results.fits";
    char gainAscii[] = "std_calib.gain";

    // Script parameters
    int   tstdmp    = 0;
    float etol      = 1.0;

    // PHA 0 line parameters (negative energy = use this |value|)
    int   npha0     = 6;
    float pha0eng[] = { 23.438, 198.392, 309.88, 584.54, 882.51, 1764.4};

    // PHA 1 line parameters (negative energy = use this |value|)
    int   npha1     = 2;
    float pha1eng[] = {2223.27, 2754.03};

    // Declare variables
    int             status  = ISDC_OK;
    int             i;
    int             k;
    int             inx;
    int             dete;
    float           eng0[100];
    float           eng1[100];
    float           amp0[100][19];
    float           amp1[100][19];
    float           pha0[100][19];
    float           pha1[100][19];
    float           err0[100][19];
    float           err1[100][19];
    long            numSct;
    long            numSrt;
    long            numRows;
    char            name[256];
    char            srtDOLSE[256];
    char            sctDOLSE[256];
    char            srtDOLME[256];
    char            sctDOLME[256];
    float          *line_energy;
    unsigned long  *det_id;
    unsigned char  *e_range;
    unsigned short *obt_start;
    unsigned short *obt_end;
    double         *ontime;
    float          *line_par;
    float          *line_par_err;
    dal_element    *srt;
    dal_element    *sct;
    FILE           *fptr;
    
    // Initialise memory pointers
    line_energy  = NULL;
    det_id       = NULL;
    e_range      = NULL;
    obt_start    = NULL;
    obt_end      = NULL;
    ontime       = NULL;
    line_par     = NULL;
    line_par_err = NULL;

    // Main loop
    do {
    
      // Initialise result fields
      for (k = 0; k < 100; k++) {
        for (dete = 0; dete < 19; dete++) {
          pha0[k][dete] = -1.0;
          pha1[k][dete] = -1.0;
          amp0[k][dete] =  0.0;
          amp1[k][dete] =  0.0;
          err0[k][dete] =  0.0;
          err1[k][dete] =  0.0;
        }
      }

      // Build DOLs
      sprintf(srtDOLSE, "%s[SPI.-LINE-SRT]", srtNameSE);
      sprintf(sctDOLSE, "%s[SPI.-LINE-SCT]", sctNameSE);
      sprintf(srtDOLME, "%s[SPI.-LINE-SRT]", srtNameME);
      sprintf(sctDOLME, "%s[SPI.-LINE-SCT]", sctNameME);

      // Open SE HDUs
      status = DALobjectOpen(sctDOLSE, &sct, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to open [SPI.-LINE-SCT] SE table.\n", status);
        continue;
      }
      status = DALobjectOpen(srtDOLSE, &srt, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to open [SPI.-LINE-SRT] SE table.\n", status);
        continue;
      }

      // Compute buffer size
      numSct = 0;
      numSrt = 0;
      status = DALtableGetNumRows(sct, &numSct, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to get [SPI.-LINE-SCT] SE table size.\n", status);
        continue;
      }
      status = DALtableGetNumRows(srt, &numSrt, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to get [SPI.-LINE-SRT] SE table size.\n", status);
        continue;
      }
      if (numSct != numSrt) {
        printf("[SPI.-LINE-SCT] and [SPI.-LINE-SRT] SE tables have incompatible "
               "size (%d <=> %d). Stop.\n", numSct, numSrt);
        continue;
      }
      numRows = numSct;      
      if (numRows < 1) {
        printf("[SPI.-LINE-SCT] and [SPI.-LINE-SRT] SE tables are empty. Stop.\n");
        continue;
      }

      // Allocate memory
      line_energy  = new float[numRows];
      det_id       = new unsigned long[numRows];
      e_range      = new unsigned char[numRows];
      obt_start    = new unsigned short[numRows*4];
      obt_end      = new unsigned short[numRows*4];
      ontime       = new double[numRows];
      line_par     = new float[numRows*10];
      line_par_err = new float[numRows*10];
      if (line_energy  == NULL ||
          det_id       == NULL || 
          e_range      == NULL || 
          obt_start    == NULL ||
          obt_end      == NULL || 
          ontime       == NULL || 
          line_par     == NULL || 
          line_par_err == NULL) {
        printf("Unable to allocate memory.\n");
        continue;
      }

      // Read SCT columns
      status = DALtableGetCol(sct, "LINE_ENERGY", 0, NULL, NULL, line_energy, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to read [SPI.-LINE-SCT] SE table.\n", status);
        continue;
      }

      // Read SRT columns
      status = DALtableGetCol(srt, "DET_ID", 0, NULL, NULL, det_id, status);
      status = DALtableGetCol(srt, "E_RANGE", 0, NULL, NULL, e_range, status);
      status = DALtableGetCol(srt, "OBT_START", 0, NULL, NULL, obt_start, status);
      status = DALtableGetCol(srt, "OBT_END", 0, NULL, NULL, obt_end, status);
      status = DALtableGetCol(srt, "ONTIME", 0, NULL, NULL, ontime, status);
      status = DALtableGetCol(srt, "LINE_PAR", 0, NULL, NULL, line_par, status);
      status = DALtableGetCol(srt, "LINE_PAR_ERR", 0, NULL, NULL, line_par_err, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to read [SPI.-LINE-SRT] SE table.\n", status);
        continue;
      }

      // Close HDUs
      status = DALobjectClose(sct, DAL_SAVE, status);
      if (status != ISDC_OK) {
        printf("%d : Error while closing [SPI.-LINE-SCT] SE table.\n", status);
        continue;
      }
      status = DALobjectClose(srt, DAL_SAVE, status);
      if (status != ISDC_OK) {
        printf("%d : Error while closing [SPI.-LINE-SRT] SE table.\n", status);
        continue;
      }

      // Loop over lines
      for (i = 0; i < numRows; i++) {
      
        // Get line index
        inx  = i*10;
        dete = det_id[i];
      
        // PHA 0
        if (e_range[i] == 0) {
          for (k = 0; k < npha0; k++) {
            if (fabs(line_energy[i] - fabs(pha0eng[k])) < etol) {
              if (line_par[inx] > amp0[k][dete]) {
                eng0[k]       = (pha0eng[k] < 0.0) ? fabs(pha0eng[k]) : line_energy[i];
                amp0[k][dete] = line_par[inx];
                pha0[k][dete] = line_par[inx+1];
                err0[k][dete] = line_par_err[inx+1];
              }
            }
          }
        } // endif: PHA 0
	
        // PHA 1
        else {
          for (k = 0; k < npha1; k++) {
            if (fabs(line_energy[i] - fabs(pha1eng[k])) < etol) {
              if (line_par[inx] > amp1[k][dete]) {
                eng1[k]       = (pha1eng[k] < 0.0) ? fabs(pha1eng[k]) : line_energy[i];
                amp1[k][dete] = line_par[inx];
                pha1[k][dete] = line_par[inx+1];
                err1[k][dete] = line_par_err[inx+1];
              }
            }
          }
        } // endelse: PHA 1
	
      } // endfor: looped over lines

      // Delete arrays
      delete [] line_energy;
      delete [] det_id;
      delete [] e_range;
      delete [] obt_start;
      delete [] obt_end;
      delete [] ontime;
      delete [] line_par;
      delete [] line_par_err;

      // Open ME HDUs
      status = DALobjectOpen(sctDOLME, &sct, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to open [SPI.-LINE-SCT] ME table.\n", status);
        continue;
      }
      status = DALobjectOpen(srtDOLME, &srt, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to open [SPI.-LINE-SRT] ME table.\n", status);
        continue;
      }

      // Compute buffer size
      numSct = 0;
      numSrt = 0;
      status = DALtableGetNumRows(sct, &numSct, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to get [SPI.-LINE-SCT] ME table size.\n", status);
        continue;
      }
      status = DALtableGetNumRows(srt, &numSrt, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to get [SPI.-LINE-SRT] ME table size.\n", status);
        continue;
      }
      if (numSct != numSrt) {
        printf("[SPI.-LINE-SCT] and [SPI.-LINE-SRT] ME tables have incompatible "
               "size (%d <=> %d). Stop.\n", numSct, numSrt);
        continue;
      }
      numRows = numSct;      
      if (numRows < 1) {
        printf("[SPI.-LINE-SCT] and [SPI.-LINE-SRT] ME tables are empty. Stop.\n");
        continue;
      }

      // Allocate memory
      line_energy  = new float[numRows];
      det_id       = new unsigned long[numRows];
      e_range      = new unsigned char[numRows];
      obt_start    = new unsigned short[numRows*4];
      obt_end      = new unsigned short[numRows*4];
      ontime       = new double[numRows];
      line_par     = new float[numRows*10];
      line_par_err = new float[numRows*10];
      if (line_energy  == NULL ||
          det_id       == NULL || 
          e_range      == NULL || 
          obt_start    == NULL ||
          obt_end      == NULL || 
          ontime       == NULL || 
          line_par     == NULL || 
          line_par_err == NULL) {
        printf("Unable to allocate memory.\n");
        continue;
      }

      // Read SCT columns
      status = DALtableGetCol(sct, "LINE_ENERGY", 0, NULL, NULL, line_energy, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to read [SPI.-LINE-SCT] ME table.\n", status);
        continue;
      }

      // Read SRT columns
      status = DALtableGetCol(srt, "DET_ID", 0, NULL, NULL, det_id, status);
      status = DALtableGetCol(srt, "E_RANGE", 0, NULL, NULL, e_range, status);
      status = DALtableGetCol(srt, "OBT_START", 0, NULL, NULL, obt_start, status);
      status = DALtableGetCol(srt, "OBT_END", 0, NULL, NULL, obt_end, status);
      status = DALtableGetCol(srt, "ONTIME", 0, NULL, NULL, ontime, status);
      status = DALtableGetCol(srt, "LINE_PAR", 0, NULL, NULL, line_par, status);
      status = DALtableGetCol(srt, "LINE_PAR_ERR", 0, NULL, NULL, line_par_err, status);
      if (status != ISDC_OK) {
        printf("%d : Unable to read [SPI.-LINE-SRT] ME table.\n", status);
        continue;
      }

      // Close HDUs
      status = DALobjectClose(sct, DAL_SAVE, status);
      if (status != ISDC_OK) {
        printf("%d : Error while closing [SPI.-LINE-SCT] ME table.\n", status);
        continue;
      }
      status = DALobjectClose(srt, DAL_SAVE, status);
      if (status != ISDC_OK) {
        printf("%d : Error while closing [SPI.-LINE-SRT] ME table.\n", status);
        continue;
      }

      // Loop over lines
      for (i = 0; i < numRows; i++) {
      
        // Get line index
        inx  = i*10;
        dete = det_id[i];
      
        // PHA 0
        if (e_range[i] == 0) {
          for (k = 0; k < npha0; k++) {
            if (fabs(line_energy[i] - fabs(pha0eng[k])) < etol) {
              if (line_par[inx] > amp0[k][dete]) {
                eng0[k]       = (pha0eng[k] < 0.0) ? fabs(pha0eng[k]) : line_energy[i];
                amp0[k][dete] = line_par[inx];
                pha0[k][dete] = line_par[inx+1];
                err0[k][dete] = line_par_err[inx+1];
              }
            }
          }
        } // endif: PHA 0
	
        // PHA 1
        else {
          for (k = 0; k < npha1; k++) {
            if (fabs(line_energy[i] - fabs(pha1eng[k])) < etol) {
              if (line_par[inx] > amp1[k][dete]) {
                eng1[k]       = (pha1eng[k] < 0.0) ? fasb(pha1eng[k]) : line_energy[i];
                amp1[k][dete] = line_par[inx];
                pha1[k][dete] = line_par[inx+1];
                err1[k][dete] = line_par_err[inx+1];
              }
            }
          }
        } // endelse: PHA 1

      } // endfor: looped over lines

      // Write ASCII file
      fptr = fopen(gainAscii, "w");
      if (fptr == NULL) {
        printf("Error while opening ASCII file %s.\n", gainAscii);
        continue;
      }
      fprintf(fptr, "%d %d %d %d\n", obt_start[0], obt_start[1], obt_start[2], 
                                     obt_start[3]);
      fprintf(fptr, "%d %d %d %d\n", obt_end[0], obt_end[1], obt_end[2], 
                                     obt_end[3]);
      fprintf(fptr, "%.1f\n", ontime[0]);
      fprintf(fptr, "%d %d\n", npha0, npha1);
      if (npha0 > 0) {
        for (i = 0; i < npha0; i++)
          fprintf(fptr, "%9.3f  ", eng0[i]);
        fprintf(fptr, "\n");
        for (i = 0; i < 19; i++) {
          if (i < 10)
            fprintf(fptr, "%d  ", i);
          else
            fprintf(fptr, "%d ", i);
          for (k = 0; k < npha0; k++) {
            if (tstdmp)
              fprintf(fptr, "%9.3f %9.3f [%.4f]  ", pha0[k][i], err0[k][i], 
                      pha0eng[k]/pha0[k][i]);
            else
              fprintf(fptr, "%9.3f %9.3f  ", pha0[k][i], err0[k][i]);
          }
          fprintf(fptr, "\n");
        }
      }
      if (npha1 > 0) {
        for (i = 0; i < npha1; i++)
          fprintf(fptr, "%9.3f  ", eng1[i]);
        fprintf(fptr, "\n");
        for (i = 0; i < 19; i++) {
          if (i < 10)
            fprintf(fptr, "%d  ", i);
          else
            fprintf(fptr, "%d ", i);
          for (k = 0; k < npha1; k++) {
            if (tstdmp)
              fprintf(fptr, "%9.3f %9.3f [%.4f]  ", pha1[k][i], err1[k][i],
                      pha1eng[k]/pha1[k][i]);
            else
              fprintf(fptr, "%9.3f %9.3f  ", pha1[k][i], err1[k][i]);
          }
          fprintf(fptr, "\n");
        }
      }
      fclose (fptr);

      // Delete arrays
      delete [] line_energy;
      delete [] det_id;
      delete [] e_range;
      delete [] obt_start;
      delete [] obt_end;
      delete [] ontime;
      delete [] line_par;
      delete [] line_par_err;
       
    } while(0);

}