// PMTNoiseProc.cc
// Contact person: Ed Callaghan <ejc3@berkeley.edu>
// See PMTNoiseProc.hh for more details.
//--------------------------------------------------------//

#include <RAT/PMTNoiseProc.hh>

using std::vector;
using namespace CLHEP;
using namespace std;

namespace RAT
{
  PMTNoiseProc::PMTNoiseProc() : Processor( "noisecalibProc" ) { }

  void PMTNoiseProc::SetI( const std::string& param, const int value  )
  {
    if (param == "simulation") {
      if (value == 1 || value == 0) {
        fSimulationFlag = value;
      } else {
        throw ParamInvalid(param,"Must be either 0 or 1.");
      }
    } else if (param == "average_min") {
      if (value >= 0)
        fAverageMinBound = value;
      else
        throw ParamInvalid(param,"Must be greater than or equal to 0.");
    } else if (param == "average_max") {
      if (value >= 0)
        fAverageMaxBound = value;
      else
        throw ParamInvalid(param,"Must be greater than or equal to 0.");
    } else {
      throw ParamUnknown(param);
    }
  }

  void PMTNoiseProc::SetS( const std::string& param, const std::string& value  )
  {
    if (param == "file") {
      fFilename = value;
      info << "PMTNoiseProc: Using output filename " << fFilename.c_str() << "\n";
    }
    else if (param == "include"){
      includeTriggers.push_back(value);
    }
    else if (param == "exclude"){
      excludeTriggers.push_back(value);
    } else {
      throw ParamUnknown(param);
    }
  }

  void PMTNoiseProc::BeginOfRun( DS::Run& )
  {
    const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
    fPMTCount = pmtInfo.GetCount();
    fHitCounter.resize(fPMTCount);

    fPGTCounter = 0; // Number of events with ony PulseGT (PGT) and EXTASY in the trigger word
    fAverageMinBound = 50; // Default cut 50Hz - 50kHz
    fAverageMaxBound = 50000;

    // Resize charge distribution counters
    fQHSArray.resize(fPMTCount);
    fQHLArray.resize(fPMTCount);

    for ( int k=0; k<fPMTCount; k++ ) {
      fQHSArray[k].clear();
      fQHLArray[k].clear();
    }

    fQHSMin.assign(fPMTCount,9999.0);
    fQHLMin.assign(fPMTCount,9999.0);
    fQHSMax.assign(fPMTCount,0.0);
    fQHLMax.assign(fPMTCount,0.0);

    fPGTMask = RAT::DU::TrigBits::GetMask("PulseGT");
    if (fSimulationFlag)
        fPGTMask |= RAT::DU::TrigBits::GetMask("EXTASY");

    DBLinkPtr db = DB::Get()->GetLink("PMTCAL");
    fBlacklistTriggers = db->GetS("noise_blacklist");
    for (size_t i = 0 ; i < excludeTriggers.size() ; i++)
        fBlacklistTriggers += "," + excludeTriggers[i];

    fBlacklist = 0; // default blacklist is all physics, prescale, pedestal
    fBlacklist |= RAT::DU::TrigBits::GetMask(fBlacklistTriggers);
    info << "PMTNoiseProc: Excluding triggers " << fBlacklistTriggers << " from noise estimate\n";
    for (size_t i = 0 ; i < includeTriggers.size() ; i++){
        fBlacklist |= !(RAT::DU::TrigBits::GetMask(includeTriggers[i]));
        info << "PMTNoiseProc: Allowing " << includeTriggers[i] << " triggers...\n";

    }
  }

  void PMTNoiseProc::EndOfRun( DS::Run& run )
  {
    const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
    const RAT::DU::ChanHWStatus& chanHWStatus = RAT::DU::Utility::Get()->GetChanHWStatus();

    // Create charge distribution histograms
    TH1F *qhsPeakDist = new TH1F("qhsPeakDist","QHS Peak Bin Value Distribution",300,-50,250);
    TH1F *qhsThreshDist = new TH1F("qhsThreshDist","QHS Threshold Value Distribution",300,-50,250);
    TH1F *qhsHHPDist = new TH1F("qhsHHPDist","QHS High Half Point Distribution",300,-50,250);
    TH1F *qhlPeakDist = new TH1F("qhlPeakDist","QHL Peak Bin Value Distribution",300,-50,250);
    TH1F *qhlThreshDist = new TH1F("qhlThreshDist","QHL Threshold Value Distribution",300,-50,250);
    TH1F *qhlHHPDist = new TH1F("qhlHHPDist","QHL High Half Point Distribution",300,-50,250);

    vector<int> qhsPeakVec;
    vector<int> qhsThreshVec;
    vector<int> qhsHHPVec;
    vector<int> qhlPeakVec;
    vector<int> qhlThreshVec;
    vector<int> qhlHHPVec;
    qhsPeakVec.resize(fPMTCount);
    qhsThreshVec.resize(fPMTCount);
    qhsHHPVec.resize(fPMTCount);
    qhlPeakVec.resize(fPMTCount);
    qhlThreshVec.resize(fPMTCount);
    qhlHHPVec.resize(fPMTCount);

    double peakQHS = 0.0;
    double threshQHS = 0.0;
    double hhpQHS = 0.0;
    double peakQHL = 0.0;
    double threshQHL = 0.0;
    double hhpQHL = 0.0;

    int rawOnlinePMTs = 0; // total number of online normal/hqe pmts

    double eventLength = DB::Get()->GetLink("DAQ")->GetD("triggergate"); // The length of an event in ns
    double averageNoise = 0.0; // Estimate for the average noise rate (not OWLs, no bounds on rate)

    double averageNHitRaw = 0.0; // Calculate average nhit for run (not OWLs, no bounds on rate)
    double averageNHitNormal = 0.0; // Calculate average nhit for run (bounded average, not OWL)
    double averageNHitTriggerEnabled = 0.0; // Calculate average nhit for normal/hqe pmts with n20 and n100 enabled

    vector<double> averageHitCounter; // PGT hit counters for every PMT averaged over the time
    averageHitCounter.resize(fPMTCount);

    for ( int i=0; i<fPMTCount; i++ ) {
      if ( pmtInfo.GetType(i) == DU::PMTInfo::NORMAL || pmtInfo.GetType(i) == DU::PMTInfo::OWL
           || pmtInfo.GetType(i) == DU::PMTInfo::HQE ) {
        if (fPGTCounter != 0)
            averageHitCounter[i] = static_cast<double>(fHitCounter[i])*second / (static_cast<double>(fPGTCounter)*eventLength);
        else
            averageHitCounter[i] = 0;

        if ( chanHWStatus.IsTubeOnline(i) && pmtInfo.GetType(i) != DU::PMTInfo::OWL ) { // Only consider normal/hqe online PMTs for average rates
          if (fPGTCounter != 0)
            averageNHitRaw += static_cast<double>(fHitCounter[i])/(static_cast<double>(fPGTCounter));
          else
            averageNHitRaw = 0;
          averageNoise += averageHitCounter[i]; // Divide by rawOnlinePMTs after sum
          rawOnlinePMTs++;

          if ( averageHitCounter[i] <= fAverageMaxBound && averageHitCounter[i] >= fAverageMinBound ){
            if (fPGTCounter != 0)
              averageNHitNormal += static_cast<double>(fHitCounter[i])/(static_cast<double>(fPGTCounter));
            else
              averageNHitNormal = 0;
          }

          if ( chanHWStatus.IsNHit20Enabled(i) && chanHWStatus.IsNHit100Enabled(i) ){
            if (fPGTCounter != 0)
              averageNHitTriggerEnabled += static_cast<double>(fHitCounter[i])/(static_cast<double>(fPGTCounter));
            else
              averageNHitTriggerEnabled = 0;
          }
        }
      }
    }

    for ( int i=0; i<fPMTCount; i++ ) {
      // Obtain peak bin value, threshold value, and high half point for QHS
      int resultQHS = GainFitPoints(i,0,peakQHS,threshQHS,hhpQHS);
      if ( resultQHS == 0 ) {
        qhsPeakDist->Fill(static_cast<int>(peakQHS),1);
        qhsThreshDist->Fill(static_cast<int>(threshQHS),1);
        qhsHHPDist->Fill(static_cast<int>(hhpQHS),1);

        qhsPeakVec[i] = static_cast<int>(peakQHS);
        qhsThreshVec[i] = static_cast<int>(threshQHS);
        qhsHHPVec[i] = static_cast<int>(hhpQHS);
      } else {
        qhsPeakDist->Fill(-50,1);
        qhsThreshDist->Fill(-50,1);
        qhsHHPDist->Fill(-50,1);

        qhsPeakVec[i] = -50;
        qhsThreshVec[i] = -50;
        qhsHHPVec[i] = -50;
      }
      // Obtain peak bin value, threshold value, and high half point for QHL
      int resultQHL = GainFitPoints(i,1,peakQHL,threshQHL,hhpQHL);
      if ( resultQHL == 0 ) {
        qhlPeakDist->Fill(static_cast<int>(peakQHL),1);
        qhlThreshDist->Fill(static_cast<int>(threshQHL),1);
        qhlHHPDist->Fill(static_cast<int>(hhpQHL),1);

        qhlPeakVec[i] = static_cast<int>(peakQHL);
        qhlThreshVec[i] = static_cast<int>(threshQHL);
        qhlHHPVec[i] = static_cast<int>(hhpQHL);
      } else {
        qhlPeakDist->Fill(-50,1);
        qhlThreshDist->Fill(-50,1);
        qhlHHPDist->Fill(-50,1);

        qhlPeakVec[i] =-50;
        qhlThreshVec[i] = -50;
        qhlHHPVec[i] = -50;
      }
    }

    std::string file_name = "QHS_QHL_" + fFilename + ".root";
    fFile = TFile::Open(file_name.c_str(),"RECREATE");
    fFile->cd();

    qhsPeakDist->GetXaxis()->SetTitle("QHS Peak (cap)");
    qhsThreshDist->GetXaxis()->SetTitle("QHS Threshold (cap)");
    qhsHHPDist->GetXaxis()->SetTitle("QHS HHP (cap)");
    qhlPeakDist->GetXaxis()->SetTitle("QHL Peak (cap)");
    qhlThreshDist->GetXaxis()->SetTitle("QHL Threshold (cap)");
    qhlHHPDist->GetXaxis()->SetTitle("QHL HHP (cap)");

    qhsPeakDist->GetYaxis()->SetTitle("Counts");
    qhsThreshDist->GetYaxis()->SetTitle("Counts");
    qhsHHPDist->GetYaxis()->SetTitle("Counts");
    qhlPeakDist->GetYaxis()->SetTitle("Counts");
    qhlThreshDist->GetYaxis()->SetTitle("Counts");
    qhlHHPDist->GetYaxis()->SetTitle("Counts");

    qhsPeakDist->Write();
    qhsThreshDist->Write();
    qhsHHPDist->Write();
    qhlPeakDist->Write();
    qhlThreshDist->Write();
    qhlHHPDist->Write();

    fFile->Close();

    if( rawOnlinePMTs > 0 ) {
      averageNoise = averageNoise/rawOnlinePMTs;
    } else {
      averageNoise = 0;
    }

    RAT::DBTable TWtable("NOISE_RUN_LEVEL");
    TWtable.SetI("version",3);
    TWtable.SetIndex("");
    TWtable.SetRunRange(run.GetRunID(), run.GetRunID());
    TWtable.SetPassNumber(-1);
    TWtable.SetI("PGT_n",fPGTCounter);
    TWtable.SetD("average_noiserate",averageNoise);
    TWtable.SetIArray("perpmt_pgt", fHitCounter);
    TWtable.SetDArray("perpmt_noiserate",averageHitCounter);
    TWtable.SaveAs("noise_rate_" + fFilename + ".ratdb");

    RAT::DBTable QHtable("NOISE_RUN_LEVEL_MONITORING");
    QHtable.SetI("version",2);
    QHtable.SetIndex("");
    QHtable.SetRunRange(run.GetRunID(), run.GetRunID());
    QHtable.SetPassNumber(-1);
    QHtable.SetIArray("qhs_peak",qhsPeakVec);
    QHtable.SetIArray("qhs_thresh",qhsThreshVec);
    QHtable.SetIArray("qhs_hhp",qhsHHPVec);
    QHtable.SetIArray("qhl_peak",qhlPeakVec);
    QHtable.SetIArray("qhl_thresh",qhlThreshVec);
    QHtable.SetIArray("qhl_hhp",qhlHHPVec);
    QHtable.SetI("online_pmts",rawOnlinePMTs);
    QHtable.SetI("PGT_n",fPGTCounter);
    QHtable.SetD("average_nhit_raw",averageNHitRaw);
    QHtable.SetD("average_nhit_normal",averageNHitNormal);
    QHtable.SetD("average_nhit_trigenabled",averageNHitTriggerEnabled);
    QHtable.SaveAs("noise_rate_monitoring_" + fFilename + ".ratdb");

    delete qhsPeakDist;
    delete qhsThreshDist;
    delete qhsHHPDist;
    delete qhlPeakDist;
    delete qhlThreshDist;
    delete qhlHHPDist;
  }

  int PMTNoiseProc::GainFitPoints(int PMTNum, int chargeFlag, double &peakBinVal,
                             double &threshVal, double &hhpVal)
  {
    // Set up the gain calculation as it was done in SNO. This can be
    // modified later when/if we think up something cooler.

    // Make sure these are 0
    peakBinVal = 0.0;
    threshVal = 0.0;
    hhpVal = 0.0;

    int peakBin = 0;
    int spanBins = 0;
    int window = 0;
    int vectorSize = 0;
    int qMax = 0;
    int qMin = 1e9;
    int peakWidth = 0;
    double maxBinSum = 0.0;
    double binSum = 0.0;
    double peakBinHeightVal = 0.0;
    double lowHPPeakAvg = 0.0;
    vector<double> localChargeVector;

    // First fill the ChargeVector with the QHS or QHL values
    // Combine all sources to take advantage of maximum statistics
    if ( chargeFlag == 0 ) {
      localChargeVector.insert(localChargeVector.end(),
                               fQHSArray[PMTNum].begin(),
                               fQHSArray[PMTNum].end());
      // std::cout<<"Size of vector for QHS dist: "<<fQHSArray[PMTNum].size()<<std::endl;
      if ( fQHSMax[PMTNum] > qMax ) qMax=fQHSMax[PMTNum];
      if ( fQHSMin[PMTNum] < qMin ) qMin=fQHSMin[PMTNum];
      // std::cout<<"QHS: qMax and qMin after pass: "<<qMax<<","<<qMin<<std::endl;
    } else if ( chargeFlag == 1 ) {
      localChargeVector.insert(localChargeVector.end(),
                               fQHLArray[PMTNum].begin(),
                               fQHLArray[PMTNum].end());
      // std::cout<<"Size of vector for QHL dist: "<<fQHLArray[PMTNum].size()<<std::endl;
      if ( fQHLMax[PMTNum] > qMax ) qMax=fQHLMax[PMTNum];
      if ( fQHLMin[PMTNum] < qMin ) qMin=fQHLMin[PMTNum];
      // std::cout<<"QHL: qMax and qMin after pass: "<<qMax<<","<<qMin<<std::endl;
    } else {
      return -1;
    }

    vectorSize = localChargeVector.size();
    if ( (qMax-qMin) < 2 ) {
      return -1;
    }
    // Set up the temporary charge histogram, for now with 1 bin per ADC bin
    TH1D *tempChargeHist = new TH1D("tempChargeHist","Charge Histogram",
                                     (qMax-qMin),qMin,qMax);
    for ( int j=0; j<vectorSize; j++ ) {
        tempChargeHist->Fill(localChargeVector.at(j));
    }

    // First identify the "span" as in J.Cameron's thesis: the 100 count wide
    // sliding window that contains the most hits
    for ( int j=0; j<tempChargeHist->GetNbinsX()-100; j++ ) {
      binSum = 0.0;
      for ( int k=0; k<100; k++ ) {
        binSum+=tempChargeHist->GetBinContent(j+k);
      }
      if ( binSum > maxBinSum ) {
        maxBinSum = binSum;
        peakBin = j;
      }
    }

    // Slide windows of increasing sizes until we found the smallest one that
    // contains 65% of the data.
    // Start with a width of 20 bins
    for ( int width=20 ;width<100 ;width++ ) {
      for ( int j=peakBin; j<peakBin+100-width; j++ ) {
        binSum = 0.0;
        for ( int k=0; k<width; k++ ) {
          binSum+=tempChargeHist->GetBinContent(j+k);
        }
        if ( binSum > 0.65*maxBinSum ) {
          peakWidth = width-1;
          width = 100;
          peakBin = j;
          break;
        }
      }
    }
    spanBins = peakWidth;
    int peakbinwidth = static_cast<int>(spanBins*0.4);
    // 0.4: as defined in J. Cameron's thesis
    // 1) Find peak
    // Peak is found by using sliding window technique,
    // the width of the window is peakbinwidth
    peakBin = 0;
    maxBinSum = 0.0;
    int maxbinnumber = tempChargeHist->GetNbinsX()-peakbinwidth;
    for ( int j=0; j<maxbinnumber; j++ ) {
      binSum = 0.0;
      for ( int k=0; k < peakbinwidth; k++ ) {
        binSum+=tempChargeHist->GetBinContent(j+k);
      }
      if (binSum > maxBinSum ) {
        maxBinSum = binSum;
        peakBin = j;
      }
    }
    int peakbinplace = peakBin+static_cast<int>(peakbinwidth*0.5);
    // Get the peak position
    peakBinVal = tempChargeHist->GetBinCenter(peakBin+(peakbinwidth*0.5));
    peakBinHeightVal = maxBinSum/static_cast<double>(peakbinwidth);

    // 2) Find HHP
    for ( int k=peakBin; k<tempChargeHist->GetNbinsX()-peakbinwidth; k++ ) {
      // Find when the sum over the peakbinwidth bins drops
      // To 1/2 of the sum over peakbinwidth bins of the max window
      binSum = 0.0;
      for ( int m=0; m<peakbinwidth; m++) {
        binSum+=tempChargeHist->GetBinContent(k+m);
      }
      if ( binSum < (0.5*maxBinSum) ) {
        // The high half point
        hhpVal = tempChargeHist->GetBinCenter(k+(peakbinwidth*0.5));
        break;
      }
    }

    // 3) Use the SNO method to find the threshold
    int new_avg_sum = 0;
    if ( spanBins < 40 ) {
      window = 1;
    } else if ( spanBins < 80 ) {
      window = 2;
    } else {
      window = 3;
    }
    for ( int k=static_cast<int>(peakbinplace); k>window; k-- ) {
      binSum = 0.0;
      int temp_low_edge = 0;
      temp_low_edge = k-(static_cast<int>(window/2.0));
      for ( int m=0; m<window; m++) {
        binSum+=tempChargeHist->GetBinContent(static_cast<int>(temp_low_edge+m));
      }
      if ( (binSum/static_cast<double>(window)) < (0.5*peakBinHeightVal) ) {
        lowHPPeakAvg = 0.0;
        for ( int n=k; n<peakbinplace; n++ ) {
          lowHPPeakAvg+=tempChargeHist->GetBinContent(n);
        }
        lowHPPeakAvg/=static_cast<double>(peakbinplace-k);
        new_avg_sum = 0;
        for ( int o=k; o>window; o-- ) {
          new_avg_sum = 0;
          temp_low_edge = o-(static_cast<int>(window/2.0));
          for ( int p=0; p<window; p++ ) {
            new_avg_sum+=static_cast<int>
            (tempChargeHist->GetBinContent(temp_low_edge+p));
          }
          if ( (new_avg_sum/static_cast<double>(window) ) < 0.5*lowHPPeakAvg ) {
            threshVal = o;
            o = window;
            k = window;
            break;
          }
        }
      }
    }
    threshVal = tempChargeHist->GetBinCenter(threshVal); // The threshold
    // Clean up
    delete tempChargeHist;
    return 0;
  }

  Processor::Result PMTNoiseProc::DSEvent(DS::Run& run, DS::Entry& ds)
  {
    int nevC = ds.GetEVCount();
    if ( nevC > 0 ) {
      RAT::DS::EV& ev = ds.GetEV(0);

      // count event if a PGT, and not on blacklist
      unsigned int trig = ev.GetTrigType();
      if ((trig & fPGTMask) && !(trig & fBlacklist)){
        int cpmtc = ev.GetCalPMTs().GetCount();
        for ( int ipmt=0; ipmt<cpmtc; ipmt++ ) {
          // find the PMTs that have fired and increment the corresponding fHitCounter by 1.
          const DS::PMTCal& pmtCal = ev.GetCalPMTs().GetPMT(ipmt);
          if( !pmtCal.GetCrossTalkFlag() ) {  // if not a cross-talk hit
            int cpmtid = pmtCal.GetID();
            double qhs = pmtCal.GetQHS();
            double qhl = pmtCal.GetQHL();
            // Keep track of min and max QHS and QHL, will come in handy later!
            if (qhs < fQHSMin[cpmtid]) fQHSMin[cpmtid]=qhs;
            if (qhs > fQHSMax[cpmtid]) fQHSMax[cpmtid]=qhs;
            if (qhl < fQHLMin[cpmtid]) fQHLMin[cpmtid]=qhl;
            if (qhl > fQHLMax[cpmtid]) fQHLMax[cpmtid]=qhl;
            // Fill the QHS and QHL vectors for this PMT
            fQHSArray[cpmtid].push_back(qhs);
            fQHLArray[cpmtid].push_back(qhl);
            fHitCounter[cpmtid]+=1;
          }
        }
        //increase the counter to keep track of how many events
        fPGTCounter+=1;
      }
    }
    return OK;
  }
} //namespace RAT