// ---------------------------------------------------------
//
// fiberQA.cc - main program to read data in from an evio
//              event file generated by the UConn TAGM dark
//              box CODA 2.6 based data acquisition system.
//              Runs are generated by executing the script
//              fiberQA_sequencer.sh to automatically bias
//              individual rows one-by-one in a 6-column
//              cluster of 30 fibers that are continously
//              being illuminated by a fast laser diode
//              pulser.
// 
// authors: richard.t.jones at uconn.edu,
// version: april 26, 2014
//
// notes:
// 1) The input is a single evio file that is assumed to
//    end with the extension ".evio" which is named as the
//    only required command line argument.
// 2) The output is a root file created in the same
//    directory as the input file, and with the same name
//    with the extension ".evio" changed to ".root". This
//    file contains the following histograms that are filled
//    using pulse data extracted from the evio data file.
//    The titles listed below indicate what is histogrammed.
//
//     name      title
//     ---------------------------------------------------
//     hnorm     laser pulse normalization factor
//     hhighN    high gain pulse height for fiber N
//     hlowN     low gain pulse height for fiber N
//     hallN     individual channel pulse height for
//               fiber N, only valid for N=1..5.
// ---------------------------------------------------------

#define TRIGGERED_BY_PULSER 1

#include <pulse_analyzer.h>
#include <netinet/in.h>
#include <evioFileChannel.hxx>
#include <evioUtil.hxx>
#include <unistd.h>
#include <vector>
#include <math.h>

// root stuff
#include <TFile.h>
#include <TH1D.h>
#include <TROOT.h>

#include "fiberQA.h"

using namespace std;
using namespace evio;

int evtCount = 0;
int totalRows = 0;
int verbose_print = 0;

// prototypes
void analyzeEvent(evioDOMTree &eventTree);
void analyzeBank(evioDOMNodeP bankPtr);
void initHistograms(const char* filename);
void fillHistograms();

// Define bit masks
// F1 bit masks
#define F1_RAW_DATA          0xFF000000

int main(int argc, char **argv) {

  if (argc == 1) {
     printf("Usage: fiberQA <path_to_evio_file>\n");
     exit(1);
  }

  char filename[128];
  strcpy(filename, argv[1]);

  try {
    char fname[128],fname1[128];
    sprintf(fname1,"%s",filename);
    strtok(fname1,".");
    sprintf(fname,"%s.root",fname1);
    initHistograms(fname);
    
    cout << "open data file " << filename << endl;
    evioFileChannel chan(filename);
    chan.open();

    while (chan.read()) {
      memset(fadc_raw,0,sizeof(fadc_raw));
      evtCount++;

#if STOP_AFTER_NEVENTS
      if (evtCount > STOP_AFTER_NEVENTS) {
         break;
      }
#endif

      if (evtCount%10000 == 0) {
         cout << "fiberQA read " << evtCount << " events" << endl;
      }
 
      evioDOMTree eventTree(chan);
      analyzeEvent(eventTree);
    }
    
    chan.close();

    cout << endl
         << "normalization pulse height mean = " << hnorm->GetMean()
         << ", rms = " << hnorm->GetRMS()
         << endl;

    ROOTfile->cd();
    ROOTfile->Write();
    ROOTfile->Close();
    delete ROOTfile;
    cout << endl << "Close ROOT file" << endl;

  } catch (evioException *e) {
    cerr << endl << e->toString() << endl << endl;
    exit(EXIT_FAILURE);
  }

  cout << endl << endl 
       << " ***fiberQA done after " << totalRows << " events analyzed***"
       << endl << endl;
  exit(EXIT_SUCCESS);
}

void analyzeEvent(evioDOMTree &eventTree) {

  evioDOMNodeListP bankList = eventTree.getNodeList();
  for_each (bankList->begin(), bankList->end(), analyzeBank);

}

void analyzeBank(evioDOMNodeP bankPtr) {

  int SlotID  =  -1;
  int nboard = -1;
  int trigger_numb = -1;

  switch (bankPtr->tag) {
    
  case 0:          // ignore
    break;

  case 3:          // header bank
    {
      const vector<int32_t> *vec = bankPtr->getVector<int32_t>();
      if (vec == NULL) {
         cerr << "?unable to get header bank vector" << endl;
         return;
      }
      if (verbose_print) {
         cout << dec << "found header bank:  run number " << (*vec)[0] << endl
                     << "                    event number " << (*vec)[1]
              << endl << endl;
      }
      break;
    }
    
  case 6:        // F1TDC data block
    {
      const vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec == NULL) {
         cerr << "?unable to get tdc bank vector" << endl;
         return;
      }

      if (verbose_print) {
        cout << " -------------------" << endl;
        cout << " F1TDC Vector Size = " << vec->size() << endl;
        cout << " -------------------" << endl;      
      }

      for (unsigned i = 0; i < vec->size(); i++){
        if (verbose_print) {
           cout << hex << (*vec)[i] << dec << endl;
        }
        if ((*vec)[i] == 0xf1daffff) continue;

        unsigned int  data_type = ((*vec)[i] & 0x800000) >> 23;
        if (verbose_print) {
          cout << " Data type = " << data_type << endl;
        }
        if (data_type == 0) { // Event header
          unsigned int event_number = ((*vec)[i] & 0x3f0000) >> 16;
          unsigned int trigger_time = ((*vec)[i] & 0xff80) >> 7;
          unsigned int trigger_fifo_overflow  =  ((*vec)[i] & 0x400000) >> 22;
          unsigned int f1_chip = ((*vec)[i] & 0x38) >> 3;
          unsigned int f1_chan = ((*vec)[i] & 0x7);  
          if (verbose_print) {
             cout << "Event Header:    Event number = " << event_number
                  << "  Trigger time = " << trigger_time
                  << " Chip = " << f1_chip << "    Channel = " << f1_chan
                  <<  " Trigger FIFO overflow = " << trigger_fifo_overflow
                  << endl;
           }
        }
        else if (data_type == 1) { // Data
           if (verbose_print) {
              cout << endl;
           }

           unsigned int slot_id  = ((*vec)[i] & 0xF8000000) >> 27;
           unsigned int chip_res_locked   =  ((*vec)[i] & 0x4000000) >> 26;
           unsigned int fifo_overflows    =  ((*vec)[i] & 0x3000000) >> 24;        
        
           unsigned int f1_chip = ((*vec)[i] & 0x380000) >> 19;
           unsigned int f1_chan = ((*vec)[i] & 0x70000) >> 16;
           unsigned int f1_time = ((*vec)[i] & 0xffff);

           if (verbose_print) {
              cout << "Data word:   Slot ID  = " << slot_id
                   <<   "  Chip = " << f1_chip <<   "  Channel = " << f1_chan
                   <<  "  Resolution Locked = " << chip_res_locked
                   <<  "  FIFO overflows =  " << fifo_overflows
                   <<  "  cTime = " << f1_time << endl;
           }

           if (chip_res_locked  == 0) {
              cout << " FATAL:  Chip Resolution is NOT locked  "
                   <<  "F1 CHIP =  " << f1_chip
                   << "  Channel =  " << f1_chan << endl;
             continue;
           }

           if (fifo_overflows  != 0) {
              cout << " FATAL:  FIFO overflows  " 
                   <<  "F1 CHIP =  " << f1_chip
                   << "  Channel =  " <<  f1_chan
                   << "  FIFO overflow " << fifo_overflows << endl;
             continue;
           }
           
           if ( (f1_chip < 0) || (f1_chip > 7) ) {
              cout << "F1TDC: Wrong Chip  Number " << f1_chip << endl;
              continue;
           }
        }
      }
      break;
    }

  case 5:    // FADC data block
    {
      const vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec == NULL) {
         cerr << "?unable to get fadc bank vector" << endl;
         return;
      }

      unsigned int data_new_type  =   0;
      unsigned int data_type      =   0;
      int channel = -1 ;
      unsigned int time_stamp     =   0;

      unsigned int pulse_number   =  -1;
      unsigned int first_sample   =  -1;

      if (verbose_print) {
         cout << " Vector Size = " << vec->size() << endl;
      }

      float adc[2] = {0,0};

      for (unsigned i = 0; i < vec->size(); i++) {
        
        // Bank header
        // bit  31       =  1
        // bits 27- 30   =  0
        // bits 22 - 26  =  SlotID 
        // bits 11 - 21  =  Number of events in block;  bits 0 - 10 event block number

        if ( (*vec)[i] >> 27 == 0x10 ) {        /* Bank header 1000 0 */
           SlotID =  ((*vec)[i] & 0x7c00000) >> 22;
           nboard++;
           if (verbose_print) {
              cout << " FADC found in slot:  " << SlotID << endl;
           }
        }
                
        // Event header
        // bit  31       =  1
        // bits 27- 30   =  0x1001
        // bits 0 - 26   =  Trigger number 

        else if ( (*vec)[i] >> 27 == 0x12 ) {       /* Event header 1001 0 */
           trigger_numb =  ((*vec)[i] & 0x7FFFFFF);
           if (verbose_print) {
              cout << " FADC found in slot:  " << SlotID
                   << " Trigger time " << trigger_numb << endl;
           }
        }

        else if ( (*vec)[i] & 0x80000000 ) {       /* data type defining word */
           data_type = ((*vec)[i] & 0x78000000) >> 27;
           // data_type == 4  - window raw data
           // data_type == 6  - pulse raw data
           data_new_type = 1;
        }

        else {
           data_new_type = 0;
        };

        if (data_type == 4 ) {         // window raw data
           if (data_new_type == 1) {
              unsigned int fadc_width = ((*vec)[i] & 0xFFF);
              channel = ((*vec)[i] & 0x7800000) >> 23;
              if (verbose_print) {
                 cout << " Data Type = " << data_type
                      << " Fadc_width  = " << fadc_width
                      << " Channel = " << channel << endl;
              }
              time_stamp = 1;
           }
           else {
              unsigned int valid_1 = 1;
              unsigned int valid_2 = 1;
              unsigned int adc_1 = ((*vec)[i] & 0x1FFF0000) >> 16;
              if ((*vec)[i] & 0x20000000) {
                 valid_1 = 0;
              }
              unsigned int adc_2 = ((*vec)[i] & 0x1FFF);
              if ( (*vec)[i] & 0x2000 ) {
                 valid_2 = 0;
              }
              adc[0] = (adc_1 < 4000)? float(adc_1) : 0;
              adc[1] = (adc_2 < 4000)? float(adc_2) : 0;

              if (verbose_print) {
                 cout << " Time stamp = " << time_stamp 
                      << " Nboard = " << nboard 
                      << " SlotID = " << SlotID 
                      << " Channel = " << channel
                      << "  ADC1 = " <<  adc[0]
                      <<  " " << adc[1]
                      <<  endl;
              }

              if ((channel >= 0) && (channel < channels)) {
                 fadc_raw[channel][time_stamp*2-2] = adc[0];
                 fadc_raw[channel][time_stamp*2-1] = adc[1];
              }
              else {
                 //cout << " Wrong channel number  " << channel
                 //     << endl; exit(EXIT_FAILURE);
              }
              time_stamp++;
           }
        }

        else if (data_type == 6 ) { // raw pulse data
           if (data_new_type == 1) {
              channel       = ((*vec)[i] & 0x7800000) >> 23;
              pulse_number  = ((*vec)[i] & 0x600000) >> 21;
              first_sample  = ((*vec)[i] & 0x3FF);

              if (verbose_print) {
                 cout << endl;
                 cout << " Data Type = " << data_type
                      << " Channel  = " << channel
                      << "  Pulse number = " <<  pulse_number
                      << "   First sample =  " << first_sample << endl;
                 cout << endl;
              }
            
              time_stamp = 1;
           }
           else {
              unsigned int valid_1 = -1;
              unsigned int valid_2 = -1;
              unsigned int adc_1 = ((*vec)[i] & 0x1FFF0000) >> 16;            
              valid_1 = ((*vec)[i] & 0x20000000 ) >> 29;
              unsigned int adc_2 = ((*vec)[i] & 0x1FFF);
              valid_2 = ( (*vec)[i] & 0x2000 ) >> 13;
              float adc[2];
              adc[0] = (adc_1 < 4000)? float(adc_1) : 0;
              adc[1] = (adc_2 < 4000)? float(adc_2) : 0;
           }
           time_stamp++;
        }

        if ( data_type == 7 ) { // Pulse Integral data
           channel    = ((*vec)[i] & 0x7800000) >> 23;
           unsigned int pulse_numb = ((*vec)[i] & 0x200000) >> 21;
           unsigned int pulse_int  = ((*vec)[i] & 0x3FFFF);

           if (verbose_print) {
              cout << " Channel = " << channel
                   << " Pulse number = " << pulse_numb
                   << " Pulse int = " << pulse_int << endl;
           }
        }
      }
      fillHistograms();
      break;
    }

  default:
    {
      break;
    }
  }

}

void initHistograms(const char* fname)
{
  
  ROOTfile = new TFile(fname,"RECREATE","new");
  cout << "Opened ROOT file " << fname << " ..." << endl;
  ROOTfile->cd();

  // extract the basename from fname

  char fname_copy[256];
  strncpy(fname_copy, fname, 255);
  char *base_fname = strtok(fname_copy,"/");
  char *next_token;
  while ((next_token = strtok(0,"/"))) {
     base_fname = next_token;
  }

  // histogram for normalization pulse

  hnorm = new TH1D("hnorm", "laser pulse normalization", 8000, -5, 1995.);

  // histograms for individual fibers in high gain mode

  for (int id=1; id <= 30; id++) {
     char name[10];
     sprintf(name,"hhigh%d",id);
     char title[256];
     sprintf(title,"fiber %d high gain from run %s",id,base_fname);
     hhigh[id] = new TH1D(name,title,8000,-5.,1995.);
  }

  // histograms for individual fibers in low gain mode (column 1 only)

  for (int id=1; id <= 5; id++) {
     char name[10];
     sprintf(name,"hlow%d",id);
     char title[256];
     sprintf(title,"fiber %d low gain from run %s",id,base_fname);
     hlow[id] = new TH1D(name,title,400,-5.,95.);
  }

  // histograms for individual channels when all lit

  for (int id=1; id <= 11; id++) {
     char name[10];
     sprintf(name,"hall%d",id);
     char title[256];
     sprintf(title,"channel %d from run %s",id,base_fname);
     if (id <= 5) {
        hall[id] = new TH1D(name,title,1000,-5.,245.);
     }
     else {
        hall[id] = new TH1D(name,title,8000,-5.,1995.);
     }
  }

}

void fillHistograms()
{
   double pheight[channels];
   pulse_analyzer pana(fadc_raw[0], channels, window_size);
   for (int channel=0; channel < 12; channel++) {
#if TRIGGERED_BY_PULSER
      pheight[channel] = pana.pheight_trig_pulser(channel);
#else
      pheight[channel] = pana.pheight_trig_summer(channel);
#endif
   }
   double pnorm = pana.pulse_normalization();
   hnorm->Fill(pheight[0]);

   double row_select_thresh = 10;
   int rows_hit = 0;
   int row_hit = 0;
   for (int row=1; row <= 5; ++row) {
      if (pheight[row] < row_select_thresh)
         continue;
      row_hit = row;
      ++rows_hit;
   }

   if (rows_hit == 1) {
      hhigh[row_hit]->Fill(pheight[6]/pnorm);
      hhigh[row_hit+5]->Fill(pheight[7]/pnorm);
      hhigh[row_hit+10]->Fill(pheight[8]/pnorm);
      hhigh[row_hit+15]->Fill(pheight[9]/pnorm);
      hhigh[row_hit+20]->Fill(pheight[10]/pnorm);
      hhigh[row_hit+25]->Fill(pheight[11]/pnorm);
      hlow[row_hit]->Fill(pheight[row_hit]/pnorm);
   }
   else if (rows_hit == 5) {
      for (int row=1; row <= 11; ++row) {
         hall[row]->Fill(pheight[row]);
      }
   }
   ++totalRows;
}