// ---------------------------------------------------------
//
// fiberQA.cc - main program to read data in from an evio
//              event file produced by running the script
//              fiberQA_sequencer.sh to automatically
//              measure the response of 15 fibers to
//              the fast laser diode pulser.
// 
// authors: richard.t.jones at uconn.edu,
// version: april 26, 2014
// ---------------------------------------------------------

#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();
    
    ROOTfile->cd();
    ROOTfile->Write();
    ROOTfile->Close();
    delete ROOTfile;
    cout<<endl<<"\nClose 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();

  // histograms for individual fibers in high gain mode

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

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

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

  // histograms for individual channels when all lit

  for (int id=1; id < 9; id++) {
     char name[10];
     sprintf(name,"hall%d",id);
     char title[256];
     sprintf(title,"channel %d from run %s",id,fname);
     if (id < 6) {
        hall[id] = new TH1D(name,title,500,-5.,95.);
     }
     else {
        hall[id] = new TH1D(name,title,500,-5.,1995.);
     }
  }

}

void fillHistograms()
{
   double pintegral[9];
   for (int channel=1; channel < 9; channel++) {
      int ped_window[] = {1, 12};
      double mean[2] = {0,0};
      double rms[2] = {0,0};
      int count = 0;
      for (int i=ped_window[0]; i < ped_window[1]; ++i) {
         double h0 = fadc_raw[channel][i-1];
         double h1 = fadc_raw[channel][window_size-i];
         mean[0] += h0;
         mean[1] += h1;
         rms[0] += h0*h0;
         rms[1] += h1*h1;
         ++count;
      }
      mean[0] /= count;
      mean[1] /= count;
      rms[0] /= count;
      rms[1] /= count;
      rms[0] = sqrt(rms[0] - mean[0]*mean[0] + 1e-10);
      rms[1] = sqrt(rms[1] - mean[1]*mean[1] + 1e-10);
      if (rms[0] > 0.03 * mean[0] || rms[1] > 0.03 * mean[1]) {
         pintegral[channel] = -99;
         continue;
      }
      double pedestal = mean[0];

      int pulse_window[2];
      if (channel > 5) {
         pulse_window[0] = 107 +28;
         pulse_window[1] = 118 +28;
      }
      else {
         pulse_window[0] = 88 +28;
         pulse_window[1] = 99 +28;
      }

      double pulse_sum = 0;
      for (int i=pulse_window[0]; i < pulse_window[1]; ++i) {
         pulse_sum += fadc_raw[channel][i];
      }
      pulse_sum = pulse_sum/(pulse_window[1] - pulse_window[0]) - pedestal;
      double adc_thresh = -8;
      if (pulse_sum < adc_thresh) {
         pintegral[channel] = -99;
         continue;
      }
      pintegral[channel] = pulse_sum;
   }

   double row_select_thresh = 3;
   double row_select_maximum = 300;
   int rows_hit = 0;
   for (int row=1; row < 6; ++row) {
      if (pintegral[row] > row_select_thresh &&
          pintegral[6] > row_select_thresh &&
          pintegral[6] < pintegral[row]*6.5+row_select_maximum)
      {
         hhigh[row]->Fill(pintegral[6]);
         hhigh[row+5]->Fill(pintegral[7]);
         hhigh[row+10]->Fill(pintegral[8]);
         hlow[row]->Fill(pintegral[row]);
         rows_hit++;
      }
      else if (pintegral[6] > pintegral[row]*6.5+row_select_maximum*1.2) {
         hall[row]->Fill(pintegral[row]);
         if (row == 1) {
            hall[6]->Fill(pintegral[6]);
            hall[7]->Fill(pintegral[7]);
            hall[8]->Fill(pintegral[8]);
         }
      }
   }
   if (rows_hit > 1) {
      cout << "weird rows hit " << rows_hit << endl;
      cout << "pintegral[6] is " << pintegral[6] << endl;
   }
   else {
      totalRows += rows_hit;
   }
}