#include <RAT/LivetimeCuts.hh>
#include <RAT/DataCleaningUtility.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/Processor.hh>
#include <RAT/DBJsonLoader.hh>
#include <RAT/DB.hh>
#include <RAT/BitManip.hh>
#include <RAT/DS/UniversalTime.hh>

#include <RAT/DBTable.hh>
#include <RAT/DBTextLoader.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/DataCleaningBits.hh>
#include <RAT/Log.hh>
#include <RAT/json.hh>

#include <cstddef>
#include <string>
#include <math.h>

using namespace std;

namespace RAT {

  void LivetimeCuts::BeginOfRun(DS::Run& run)
  {
    fPassFlag = true;
    fFoundRetriggers = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("retriggercut");
    fFoundNhitBurst = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("tpburstcut");
    fFoundMuon = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("tpmuonfollowercut-short");
    fFoundCLBurst = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("missingcaendata");
    fFoundMC = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("missedcountburstcut");
    fFoundAtmospheric = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("atmospheric");
    fFoundPedestal = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("pedcut");
    fFoundMissedMuon = DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("missedmuonfollower");



    fRetriggerCutGTID[0] = -1;
    fRetriggerCutGTID[1] = -1;
    fRetriggerCutTime[2] = 0.0;//initialize accumulated value to 0.

    //Initialization customized burst type
    std::string burstType = "";
    if (fBurstType == "ambe") {
      burstType = "_" + fBurstType;
    } else {
      if (fBurstType != "") {
        warn << "LivetimeCuts::BeginOfRun: burst type " << fBurstType <<" not recognized, default burst parameters will be used." << endl;
      }
    }

    //Read Analysis Mask and compare with livetime mask
    if (fLivetimeMask == -1) {
      fLivetimeMask = GetDataCleaningWord("analysis_mask");
    }

    std::vector<int> temp;
    if (!(ReadFile(run, "retriggercut", fFirstRetriggerReadout, fLastRetriggerReadout, temp))) {
      fFoundRetriggers = -1;
    }
    if (!(ReadFile(run, "burstcut", fFirstNhitBurstReadout, fLastNhitBurstReadout, fNhitBurstReadout))) {
      fFoundNhitBurst = -1;
    }
    if (!(ReadFile(run, "caenlosscut", fFirstCLBurstReadout, fLastCLBurstReadout, fCLBurstReadout))) {
      fFoundCLBurst = -1;
    }
    temp.clear();

    DBLinkPtr fLBurst = DB::Get()->GetLink("DATACLEANING","tpburstcut");
    DBLinkPtr fLRetrigger = DB::Get()->GetLink("DATACLEANING","retriggercut");
    DBLinkPtr fLMuon = DB::Get()->GetLink("DATACLEANING","tpmuonfollowercut");
    DBLinkPtr fLCAEN = DB::Get()->GetLink("DATACLEANING","caenlosscut");
    DBLinkPtr fLMMF = DB::Get()->GetLink("DATACLEANING","missedmuonfollower");
    DBLinkPtr fLPed = DB::Get()->GetLink("DATACLEANING","pedcut");
    DBLinkPtr fLAtm = DB::Get()->GetLink("DATACLEANING","atmospheric");

    try {
      DBLinkPtr fRun = DB::Get()->GetLink("RUN");
      fRunType = fRun->GetI("runtype");
      fStart_times[0] = double(fRun->GetI("start_day"));//Initialize start time to this value, so it wouldn't
      fStart_times[1] = double(fRun->GetI("start_sec"));//read 0 if the frist valid GTID is not found.
      fStart_times[2] = fRun->GetD("start_nsc");
      DS::UniversalTime stopTemp((UInt_t)fRun->GetI("stop_day"), (UInt_t)fRun->GetI("stop_sec"), (UInt_t)fRun->GetD("stop_nsc"));
      DS::UniversalTime startTemp((UInt_t)fRun->GetI("start_day"), (UInt_t)fRun->GetI("start_sec"), (UInt_t)fRun->GetD("start_nsc"));
      fLength = 9999999999.0;
      fLength = ExtractTimeIndex(stopTemp);
    }

    catch(DBNotFoundError &e){
      detail << "LivetimeCuts::BeginOfRun: Could not find RUN table, no Livetime Cut will be applied." << endl;
      fStart_times[0] = -1.;//Initialize start time to this value, so it wouldn't
      fStart_times[1] = -1.;//read 0 if the frist valid GTID is not found.
      fStart_times[2] = -1.;
      fLength = -1.;
     }

    fValidFirstGTID = run.GetValidGTID();
    fRunNumber = run.GetRunID();
    fRunType = run.GetRunType();
    fNThreshold = fLBurst->GetI("nhit_threshold");
    //If this statement is true, it means the RUN table tell us that the current run's run time is
    //invalid (refer to /src/daq/ProduceRunTableProc.cc line 213-241 for more detail):
    //
    // case 1:
    //   warn<<" ProduceRunTableProc: Start or stop times are zero. Setting invalid times."<<newline;
    //   break;
    // case 2:
    //   warn<<" ProduceRunTableProc: Start or stop times are invalid. Setting invalid times."<<newline;
    //   break;
    // case 3:
    //   warn<<" ProduceRunTableProc: Run length is zero. Setting invalid times."<<newline;
    //   break;
    // case 4:
    //   warn<<" ProduceRunTableProc: Run length largely differs from DetectorDB. Setting invalid times."<<newline;
    //   break;

    std::stringstream stream;
    stream<<"LIVETIME_"<<fRunNumber<<".ratdb";
    fLivetimeTableName = stream.str();

    fNhitBurstWindow = static_cast<double>(fLBurst->GetI("twin_nhitburst" + burstType)) * 1e-6;
    fRetriggerWindow = static_cast<double>(fLRetrigger->GetI("twin_retrigger")) * 1e-6;
    fMuonWindow = static_cast<double>(fLMuon->GetD("short_follower_time"));
    fCAENWindow = static_cast<double>(fLCAEN->GetD("twin_caenloss")) * 1e-6;
    fMissedMuonWindow = static_cast<double>(fLMMF->GetD("time_window")) * 1e-9;
    fPedWindow = static_cast<double>(fLPed->GetD("time_window")) * 1e-9;
    fAtmWindow = (fLAtm->GetD("cut_time_window")) * 1e-9;



    //Initializing all Follower Type Cut parameter
    fMuonCutTime = 0.0;//Muon
    fAtmCutTime = 0.0;
    fMissedMuonCut = 0.0;//MissedMuon
    fMissedMuonCutTime = 0.0;
    fMCCut = 0.0;//Missed Count Burst
    fMCCutTime = 0.0;

    //Initializing all Livetime Type Cut Parameter
    fNhitBurstCutGTID[0] = -1;//Nhit Burst Cut
    fNhitBurstCutGTID[1] = -1;
    fNhitBurstCutTime[2] = 0.0;
    fCLBurstCutGTID[0] = -1;//CAEN Loss Livetime Cut
    fCLBurstCutGTID[1] = -1;
    fCLBurstCutTime[2] = 0.0;

    //Initializing Pedestal Cut
    fPedestalCutTime[0] = 0.0;
    fPedestalCutTime[1] = 0.0;
    fPedestalCutTime[2] = 0.0;

    fRateCount[0] = 0;
    fRateCount[1] = 0;
    fRateCount[2] = 0;

    fStop_times = 0.0;


  try{//Reading Muon table.
      DBLinkPtr fMuonData = DB::Get()->GetLink("TPMUONFOLLOWER");
      std::vector<int> fMuonDay = fMuonData->GetIArray("day_muons");
      std::vector<int> fMuonSec = fMuonData->GetIArray("sec_muons");
      std::vector<int> fMuonNsc = fMuonData->GetIArray("ns_muons");
      if ((fMuonSec.size() == fMuonDay.size()) && (fMuonSec.size() == fMuonNsc.size())) {
        for (size_t muonIndex = 0; muonIndex < fMuonDay.size(); muonIndex++){
          DS::UniversalTime muonCutStart((UInt_t)fMuonDay[muonIndex], (UInt_t)fMuonSec[muonIndex], (UInt_t)fMuonNsc[muonIndex]);
          double muonStartDouble = ExtractTimeIndex(muonCutStart);
          double muonEndDouble = ExtractTimeIndex(muonCutStart, fMuonWindow);
          fMuonCutTime += PushToCuttingRecord(muonStartDouble, muonEndDouble, fFoundMuon, true);
        }
      }
      std::vector<int> fMuonReadout = fMuonData->GetIArray("gtid_muons");
      std::vector<int>::iterator fakeMuonSearch = find (fMuonReadout.begin(), fMuonReadout.end(), -1);
      if (fakeMuonSearch == fMuonReadout.end()) {
        try{//Reading the last muon table
           DBLinkPtr fLastMuonData = DB::Get()->GetLink("LAST_MUON");
           if ((fLastMuonData->GetI("day_last_muon") != -1) &&
              (fLastMuonData->GetI("sec_last_muon") != -1) &&
              (fLastMuonData->GetI("ns_last_muon") != -1) ) {
                UInt_t cutLastMuonDay = static_cast<unsigned int>( fLastMuonData->GetI("day_last_muon") );
                UInt_t cutLastMuonSec = static_cast<unsigned int>( fLastMuonData->GetI("sec_last_muon") );
                UInt_t cutLastMuonNs  = static_cast<unsigned int>( fLastMuonData->GetI("ns_last_muon") );
                DS::UniversalTime lastMuonTemp(cutLastMuonDay, cutLastMuonSec, cutLastMuonNs);
                double fLastMuonBegin = ExtractTimeIndex(lastMuonTemp, fMuonWindow);
                double fLastMuonEnd = ExtractTimeIndex(lastMuonTemp, fMuonWindow);
                fMuonCutTime += PushToCuttingRecord(fLastMuonBegin, fLastMuonEnd, fFoundMuon, true);
                detail << "LivetimeCuts::BeginOfRun: Got last muon in previous run." << endl;
           }
         }
         catch(DBNotFoundError &e){
           detail << "LivetimeCuts::BeginOfRun: Could not find LAST_MUON ratdb file." << endl;
         }
      }
    }
    catch(DBNotFoundError &e){
      warn << "LivetimeCuts::BeginOfRun: Could not find TPMUONFOLLOWER ratdb file. Muon Cut will not be applied." << endl;
    }

    try{//Reading Missed Muon table.
      DBLinkPtr fMissedMuonData = DB::Get()->GetLink("MISSEDMUONFOLLOWER");
      fMissedMuonReadout = fMissedMuonData->GetIArray("gtid_missedmuons");
      if (fMissedMuonReadout.size() == 0) {
        fFoundMissedMuon = -1;
      }
    }
    catch(DBNotFoundError &e){
      warn << "LivetimeCuts::BeginOfRun: Could not find MISSEDMUONFOLLOWER ratdb file. Missed Muon Cut will not be applied." << endl;
    }


    try{//Reading the last atmospheric neutrino table
     DBLinkPtr fLastAtm = DB::Get()->GetLink("LAST_ATMOSPHERIC");
     if ((fLastAtm->GetI("last_day_atmospheric") != -1) &&
        (fLastAtm->GetI("last_sec_atmospheric") != -1) &&
        (fLastAtm->GetI("last_nsec_atmospheric") != -1)) {
          UInt_t cutLastAtmDay = static_cast<unsigned int>( fLastAtm->GetI("day_last_atmospheric") );
          UInt_t cutLastAtmSec = static_cast<unsigned int>( fLastAtm->GetI("sec_last_atmospheric") );
          UInt_t cutLastAtmNs  = static_cast<unsigned int>( fLastAtm->GetI("ns_last_atmospheric") );
          DS::UniversalTime lastAtmTemp(cutLastAtmDay, cutLastAtmSec, cutLastAtmNs);
          double fLastAtmTimeBegin = ExtractTimeIndex(lastAtmTemp);
          double fLastAtmTimeEnd = ExtractTimeIndex(lastAtmTemp, fAtmWindow);
          fAtmCutTime += PushToCuttingRecord(fLastAtmTimeBegin, fLastAtmTimeEnd, fFoundAtmospheric, true);
          detail << "LivetimeCuts::BeginOfRun: Got Last Atmospheric Neutrino in previous run." << endl;
     }
   }
   catch(DBNotFoundError &e){
     warn << "LivetimeCuts::BeginOfRun: Could not find Last Atmospheric Neutrino ratdb file." << endl;
   }

   try{//Reading Atmospheric Neutrino table.
      DBLinkPtr fAtmNuData = DB::Get()->GetLink("ATMOSPHERICS");
      std::vector<int> fAtmNuDay = fAtmNuData->GetIArray("days_atmospheric");
      std::vector<int> fAtmNuSec = fAtmNuData->GetIArray("secs_atmospheric");
      std::vector<int> fAtmNuNsc = fAtmNuData->GetIArray("nsecs_atmospheric");
      if ((fAtmNuSec.size() == fAtmNuDay.size()) && (fAtmNuSec.size() == fAtmNuNsc.size())) {
        for (size_t atmIndex = 0; atmIndex < fAtmNuDay.size(); atmIndex++){
          DS::UniversalTime atmCutStart((UInt_t)fAtmNuDay[atmIndex], (UInt_t)fAtmNuSec[atmIndex], (UInt_t)fAtmNuNsc[atmIndex]);
          double atmStartDouble = ExtractTimeIndex(atmCutStart);
          double atmEndDouble = ExtractTimeIndex(atmCutStart, fAtmWindow);
          fAtmCutTime += PushToCuttingRecord(atmStartDouble, atmEndDouble, fFoundAtmospheric, true);
        }
      }
    } catch(DBNotFoundError &e){
      warn << "LivetimeCuts::BeginOfRun: Could not find Atmospheric Neutrino ratdb file. Atmospheric Neutrino Cut will not be applied." << endl;
    }

    //Missed Count Burst
    DBLinkPtr fMCBurst = DB::Get()->GetLink("MISSED_COUNT_BURST");
    fMCTimeCut = (fMCBurst->GetD("time_cut")) * 1e-9;;
    try{
      DBLinkPtr fMissedCountBurst = DB::Get()->GetLink("MISSED_BURST");
      fMCGTIDs = fMissedCountBurst->GetIArray("missed_counts_burst_gtid");
    }
    catch(DBNotFoundError &e){
      warn << "MissedCountsBurst::BeginOfRun: Cound not find burst ratdb file." << newline;
      fFoundMC = -1;
    }
  }

  double LivetimeCuts::ExtractTimeIndex(DS::UniversalTime inputTime, double offset) {
    std::vector<double> TimeDiff(3);
    TimeDiff[0] = (double)inputTime.GetDays() - (double)fStart_times[0];
    TimeDiff[1] = (double)inputTime.GetSeconds() - (double)fStart_times[1];
    TimeDiff[2] = (double)inputTime.GetNanoSeconds() - (double)fStart_times[2];

    double timeDifference = (double)((TimeDiff[0])*86400) + (double)(TimeDiff[1]) + (double)((TimeDiff[2]) * 1.0E-9);
    timeDifference += offset;

    if ((timeDifference < 0) || (timeDifference > fLength)) {
      return -1.0;
    }

    return timeDifference;
  }

  double LivetimeCuts::PushToCuttingRecord(double begin, double end, int dcbit, bool softPush) {
    if (!((fLivetimeMask >> dcbit) & 1)) {
      return 0.0;
    }
    if (softPush) {
      if ((begin == -1.0) && (end == -1.0)) {
        return 0.0;
      }
      begin = std::max(begin, 0.0);
      if (end == -1.0) {
        end = fLength;
      }
    }
    bool validBegin = (begin >= 0.0) && (begin < fLength);
    bool validEnd = (end >= 0.0) && (end > begin);

    if ((validBegin) && (validEnd)) {
      fCutStart.push_back(begin);
      fCutEnd.push_back(end);
      return double(end - begin);
    }
    return 0.0;
  }


  bool LivetimeCuts::ReadFile(DS::Run& run,
                              std::string cutName,
                              std::vector<int>& firstReadout,
                              std::vector<int>& lastReadout,
                              std::vector<int>& burstReadout) {
    bool foundDBTable = false;
    bool foundJSONFile = false;
    bool foundCut = false;
    try{
      DBLinkPtr fData = DB::Get()->GetLink("LIVETIME_CUT",cutName);
      if (fData->GetI("cutapplied") == 1) {
        firstReadout = fData->GetIArray("firstburst");
        lastReadout = fData->GetIArray("lastburst");
        if (cutName != "retriggercut") burstReadout = fData->GetIArray("burst");
        foundCut = true;
      } else {
        std::string readDetail = "LivetimeCuts::BeginOfRun: No livetime is cut by " + cutName;
        detail << readDetail << endl;
      }
      foundDBTable = true;
    }

    catch(DBNotFoundError &e){
      std::string readWarn = "LivetimeCuts::BeginOfRun: Could not find " + cutName + " info, Looking for .json file";
      warn << readWarn << endl;
    }

    if (!foundDBTable) {
      std::ostringstream fname;
      fname << "livetimecut_" << cutName << "_" << run.GetRunID() << ".json";
      try {
        std::vector<DBTable *> contents = DBJsonLoader::parse(fname.str());
        json::Value firstEventTemp = contents[0]->GetJSON("firstburst");
        json::Value lastEventTemp = contents[0]->GetJSON("lastburst");
        json::Value burstEventTemp = contents[0]->GetJSON("burst");
        foundJSONFile = true;

        if ((int)firstEventTemp.getArraySize() != 0 ) {
          foundCut = true;
          for (int i=0; i<(int)burstEventTemp.getArraySize(); i++) {
            int gtid = burstEventTemp[i].cast<int>();
            burstReadout.push_back(gtid);
          }

          for (int j = 0; j < (int)firstEventTemp.getArraySize(); j++) {
            int firstgtid = firstEventTemp[j].cast<int>();
            int lastgtid = lastEventTemp[j].cast<int>();
            firstReadout.push_back(firstgtid);
            lastReadout.push_back(lastgtid);
          }
        } else {
          std::string readJSONDetail = "LivetimeCuts::BeginOfRun: No livetime is cut by " + cutName;
          detail << readJSONDetail << endl;
        }
      }

      catch(FileError &e){
        std::string readJSONWarn = "LivetimeCuts::BeginOfRun: Could not find " + cutName + ".json file";
        warn << readJSONWarn << endl;
      }
    }

    if ((!foundDBTable) && (!foundJSONFile)) {
      std::string readNoInfo = "LivetimeCuts::BeginOfRun: Could not find " + cutName + " info, " + cutName + " will not be applied..";
      warn << readNoInfo << endl;
    }
    return foundCut;
  }

  Processor::Result LivetimeCuts::DSEvent(DS::Run&, DS::Entry& ds)
  {
    bool pass = true; // If any triggered event fails, they all fail
    for( size_t iEV = 0; iEV < ds.GetEVCount(); iEV++ ) {
      if( Event(ds, ds.GetEV(iEV)) != OKTRUE )
        pass = false;
    }
    return pass ? OKTRUE : OKFALSE;
  }

  Processor::Result LivetimeCuts::Event(DS::Entry&, DS::EV& ev)
  {

    fPassFlag = true;

    if (fValidFirstGTID != -1) {
      if (ev.GetGTID() == fValidFirstGTID) {
        fStart_times[0] = (double)ev.GetUniversalTime().GetDays();
        fStart_times[1] = (double)ev.GetUniversalTime().GetSeconds();
        fStart_times[2] = (double)ev.GetUniversalTime().GetNanoSeconds();
      }
      fValidFirstGTID = -1;
    }
    //Do nothing with orphans or event before first valid events.
    if (ev.GetTrigType() == 0) return OKFALSE;

    //if the event is forced event, only update stop time but do nothing.
    if (ForcedEventsTagger(ev)) {
      fStop_times = ExtractTimeIndex(ev.GetUniversalTime());
      return OKFALSE;
    }

    //Call retrigger identifier. If the event is retrigger event, update mask. Call this before
    //the acutal fStop_times update because RetriggerIdentifier needs previous event's time.
    if(RetriggerIdentifier(ev)) {
      fPassFlag = false;
    } else {
        if (ev.GetNhits() >= fNThreshold) {
          fRateCount[1]++;
        }
        if (!ev.DigitiserExists()) {
          fRateCount[2]++;
        }
    }

    //Update fStop_times
    fStop_times = ExtractTimeIndex(ev.GetUniversalTime());

    //Call burst identifier for NhitBurst events and update mask if necessary.
    if (BurstIdentifier(ev, fNhitBurstWindow, fNhitBurstCutGTID, fNhitBurstCutTime, fFirstNhitBurstReadout, fLastNhitBurstReadout, fNhitBurstReadout, fFoundNhitBurst)) {
      fPassFlag = false;
    }

    //Call burst identifier for CAEN Loss Burst events, no UpdateMask.
    if (BurstIdentifier(ev, fCAENWindow, fCLBurstCutGTID, fCLBurstCutTime, fFirstCLBurstReadout, fLastCLBurstReadout, fCLBurstReadout, fFoundCLBurst)) {
      fPassFlag = false;
    }

    //Call the follower cut identifier for missed count follower cut
    if (FollowerIdentifier(ev, fMCGTIDs, fMCTimeCut, fMCCutTime, fMCCut, fFoundMC)) {
      fPassFlag = false;
    }

    //Call missed muon identifier for MissedMuonFollower Cut, no UpdateMask
    if (FollowerIdentifier(ev, fMissedMuonReadout, fMissedMuonWindow, fMissedMuonCutTime, fMissedMuonCut, fFoundMissedMuon)) {
      fPassFlag = false;
    }

    return fPassFlag ? OKFALSE : OKTRUE;
  }

void LivetimeCuts::EndOfRun(DS::Run& run) {

    fPedestalCutTime[2] += PushToCuttingRecord(fPedestalCutTime[0], fPedestalCutTime[1], fFoundPedestal, true);

    RAT::DBTable *LIVETIME_table = new DBTable("LIVETIME");
    LIVETIME_table->SetRunRange(fRunNumber,fRunNumber);
    LIVETIME_table->SetPassNumber(-1);
    LIVETIME_table->SetS("index","");
    LIVETIME_table->SetI("version",2);
    LIVETIME_table->SetI("runtype",fRunType);

    if ((fCutStart.size() != fCutEnd.size()) || (fLength == -1.)) {
      warn << "LivetimeCuts::EndOfRun: Livetime Calculation Failed!" << "\n";
      warn << "LivetimeCuts::EndOfRun: Possible Cause 1: Cut Start Array and Cut End Array size doesn't match." << "\n";
      warn << "LivetimeCuts::EndOfRun: Possible Cause 2: ProduceRunTableProc indicates that the given run has an invalid run length." << "\n";

      LIVETIME_table->SetD("raw_livetime",-1);
      LIVETIME_table->SetD("retrigger_cut",-1);
      LIVETIME_table->SetD("nhit_burst_cut",-1);
      LIVETIME_table->SetD("caen_loss_burst_cut",-1);
      LIVETIME_table->SetD("overlap",-1);
      LIVETIME_table->SetD("muon",-1);
      LIVETIME_table->SetD("pedestal",-1);
      LIVETIME_table->SetD("missedmuon",-1);
      LIVETIME_table->SetD("atmospheric",-1);
      LIVETIME_table->SetD("missedcountburst",-1);
      LIVETIME_table->SetD("livetime",-1);
      LIVETIME_table->SetD("retrigger_cut_sacrifice",-1);
      LIVETIME_table->SetD("nhit_burst_cut_sacrifice",-1);
      LIVETIME_table->SetD("caen_loss_burst_cut_sacrifice",-1);

    } else {
      double rawCutTime = 0.0;
      double cutTime = 0.0;

      for(size_t i = 0; i < fCutStart.size(); i++) {
        rawCutTime += fCutEnd[i] - fCutStart[i];
      }

      //The Overlap Reservior make sure that we do not overcount if more
      //than 2 livetime cuts happened during the same period. This algorithm will
      //make sure that no two cut will occupy the same period of time.
      //==========================================================================
      //Define a few parameters:
      //   begin  :begin of the index in the entire cutting pool(e.g 0)
      //   end    :end of the index in the entire cutting pool(e.g fCutStart.size(),
      //           if there are 10 cuts in the cutting pool, then end will be 9)
      //   i      :The index of the cut we are working on.(e.g. 6 if we are working on the 7th cut in the cutting pool)
      //This algorithm contains a double-iteration structure. First, a 'current' cut is selected by iterating from [begin]
      //to [end]. For each selected current cut at index [i], the second level iteration is performed from [i] to [end].
      //This second-level iteration gives us a 'target' cut. 'target' cut and 'current' cut are compared to eliminate any
      //overlap, the 4 possible overlapping situation is discussed and illustrated in the following code.
      for(size_t i = 0; i < fCutStart.size(); i++) {
        double currentStartTime = fCutStart[i];
        double currentEndTime = fCutEnd[i];

        //ignoring zero cut
        if (currentStartTime == currentEndTime) continue;

          //index j belongs to the target cut we want to compare
        for(size_t j = (i + 1); j < fCutStart.size(); j++) {

          //Ignore all zero-cut to increase processing speed
          if (fCutStart[j] == fCutEnd[j]) continue;

          //No overlap
          if ((currentStartTime >= fCutEnd[j]) || (currentEndTime <= fCutStart[j])) continue;

          if ((currentStartTime <= fCutStart[j]) && (currentEndTime >= fCutEnd[j])) {
            //if current cut fully covers target cut, collapse the target cut into a 0 livetime cut.
            //Before: target:       =====
            //       current:   ==============
            //After:  target:
            //       current:   ==============
            fCutEnd[j] = fCutStart[j];
          } else if ((currentStartTime >= fCutStart[j]) && (currentEndTime <= fCutEnd[j])) {
            //if target cut fully covers current cut, collapse the current cut:
            //Before: target:   ==============
            //       current:       ======
            //After:  target:   ==============
            //       current:
            fCutEnd[i] = fCutStart[i];
            break;
          } else if (currentEndTime >= fCutEnd[j]) {
            //if target cut and current cut partially overlap:
            //Before: target:   ==========
            //       current:         ========
            //After:  target:   ======
            //       current:         ========
            fCutEnd[j] = currentStartTime;
          } else {
            //if target cut and current cut partially overlap:
            //Before: target:   ==========
            //       current: ======
            //After:  target:       ======
            //       current: ======
            fCutStart[j] = currentEndTime;
          }
        }
        cutTime += (fCutEnd[i] - fCutStart[i]);
      }

      double overlap = rawCutTime - cutTime;
      double retrigger_sacrifice = CalculateSacrifice(0, fLength);
      double burst_sacrifice = CalculateSacrifice(1, fLength);
      double clburst_sacrifice = CalculateSacrifice(2, fLength);
      double livetime = fLength - cutTime;

      LIVETIME_table->SetD("raw_livetime",fLength);
      LIVETIME_table->SetD("retrigger_cut",fRetriggerCutTime[2]);
      LIVETIME_table->SetD("nhit_burst_cut",fNhitBurstCutTime[2]);
      LIVETIME_table->SetD("caen_loss_burst_cut",fCLBurstCutTime[2]);
      LIVETIME_table->SetD("overlap",overlap);
      LIVETIME_table->SetD("muon",fMuonCutTime);
      LIVETIME_table->SetD("pedestal",fPedestalCutTime[2]);
      LIVETIME_table->SetD("missedmuon",fMissedMuonCutTime);
      LIVETIME_table->SetD("atmospheric",fAtmCutTime);
      LIVETIME_table->SetD("missedcountburst",fMCCutTime);
      LIVETIME_table->SetD("livetime",livetime);
      LIVETIME_table->SetD("retrigger_cut_sacrifice", retrigger_sacrifice);
      LIVETIME_table->SetD("nhit_burst_cut_sacrifice",burst_sacrifice);
      LIVETIME_table->SetD("caen_loss_burst_cut_sacrifice",clburst_sacrifice);
    }
    LIVETIME_table->SaveAs(fLivetimeTableName);

    // Clear all vectors
    fCutStart.clear();
    fCutEnd.clear();
    fFirstNhitBurstReadout.clear();
    fLastNhitBurstReadout.clear();
    fNhitBurstReadout.clear();
    fFirstCLBurstReadout.clear();
    fLastCLBurstReadout.clear();
    fCLBurstReadout.clear();
    fFirstRetriggerReadout.clear();
    fLastRetriggerReadout.clear();
    fDCBitCollection.clear();
  }

  bool LivetimeCuts::BurstIdentifier(DS::EV& ev,
                                     double fBurstWindow,
                                     int fBurstCutGTID[2],
                                     double fBurstCutTime[3],
                                     std::vector<int>& fFirstBurstReadout,
                                     std::vector<int>& fLastBurstReadout,
                                     std::vector<int>& fBurstReadout,
                                     int fBurstBit) {

    std::vector<int>::iterator burstSearch;
    burstSearch = find (fBurstReadout.begin(), fBurstReadout.end(), ev.GetGTID());
    if (burstSearch == fBurstReadout.end()) return false;

    if (fProcessPass == 1) {
      burstSearch = find (fFirstBurstReadout.begin(), fFirstBurstReadout.end(), ev.GetGTID());
      if (burstSearch != fFirstBurstReadout.end()) {
        fBurstCutGTID[0] = ev.GetGTID();
        fBurstCutTime[0] = ExtractTimeIndex(ev.GetUniversalTime());
      }

      burstSearch = find (fLastBurstReadout.begin(), fLastBurstReadout.end(),ev.GetGTID());
      if (burstSearch != fLastBurstReadout.end()){
        fBurstCutGTID[1] = ev.GetGTID();
        fBurstCutTime[1] = ExtractTimeIndex(ev.GetUniversalTime());
      }

      //Burst Cut allows the Last Burst event to happen prior to First Burst Event, so wait
      //until both events are set to output livetime cut.
      if ((fBurstCutGTID[0] != -1) && (fBurstCutGTID[1] != -1)) {
        //The first_last event pair has to pass correspondence check. if not, abandon this pair.
        if (CorrespondenceCheck(fBurstCutGTID, fFirstBurstReadout, fLastBurstReadout)) {
          fBurstCutTime[2] += PushToCuttingRecord((fBurstCutTime[0] - fBurstWindow),(fBurstCutTime[1] + fBurstWindow), fBurstBit);
        }
        fBurstCutGTID[0] = -1;
        fBurstCutGTID[1] = -1;
      }
    }

    return true;
  }

  bool LivetimeCuts::RetriggerIdentifier(DS::EV& ev) {

    std::vector<int>::iterator retriggerFirstSearch = find (fFirstRetriggerReadout.begin(), fFirstRetriggerReadout.end(), ev.GetGTID());
    std::vector<int>::iterator retriggerLastSearch = find (fLastRetriggerReadout.begin(), fLastRetriggerReadout.end(), ev.GetGTID());

    if (retriggerFirstSearch != fFirstRetriggerReadout.end()) {//If input event is first event
      fRetriggerCutGTID[0] = ev.GetGTID();
      fRetriggerCutTime[0] = ExtractTimeIndex(ev.GetUniversalTime());
      return false;
    } else if (retriggerLastSearch != fLastRetriggerReadout.end()) {//if input event is last event
      fRetriggerCutGTID[1] = ev.GetGTID();
      fRetriggerCutTime[1] = ExtractTimeIndex(ev.GetUniversalTime());
      if (CorrespondenceCheck(fRetriggerCutGTID, fFirstRetriggerReadout, fLastRetriggerReadout)) {
        fRetriggerCutTime[2] += PushToCuttingRecord(fRetriggerCutTime[0],(fStop_times + fRetriggerWindow), fFoundRetriggers);
      }
      fRetriggerCutGTID[0] = -1;
      fRetriggerCutGTID[1] = -1;
      return true;
    } else {
      //if first retrigger GTID is set, it means input event is a retrigger event, otherwise it's not a retrigger event.
      return (fRetriggerCutGTID[0] != -1);
    }
  }

  //Define two type of event that is relevant to this identifier:
  //Muon: the event itself is a muon
  //Muon follower: the event is not(necessarily) a muon, but it happened within fMuonWindow of a muon
  bool LivetimeCuts::FollowerIdentifier(DS::EV& ev,
                                    std::vector<int>& muonReadout,
                                    double muonWindow,
                                    double& muonCutTime,
                                    double& muonCut,
                                    int followerBit) {

    std::vector<int>::iterator muonSearch;
    muonSearch = find (muonReadout.begin(), muonReadout.end(), ev.GetGTID());
    double currentTime = ExtractTimeIndex(ev.GetUniversalTime());

    if (muonSearch != muonReadout.end()) {//if the input event is muon
      if (PushToCuttingRecord(currentTime, currentTime + muonWindow, followerBit) != 0.0) {
      //This deal with the situation that a event is both a muon and a muon follower, to crop out
      //the overlap between two muon livetime cut in accumulated muon cut time.
        muonCutTime += muonCut > currentTime? (currentTime + muonWindow - muonCut) : muonWindow;
        muonCut = currentTime + muonWindow;
      }
      return false;
    } else if (muonCut != 0.0) {
      if (currentTime <= muonCut) {//if the input event is a muon follower
        return true;
      } else {//First event out of muon follower time window, reset muonCut to null.
        muonCut = 0.0;
        return false;
      }
    } else {//The event is neither muon nor muon follower.
      return false;
    }
  }

  bool LivetimeCuts::ForcedEventsTagger(DS::EV& ev) {

    Int_t mtcd_trig_word = ev.GetTrigType();
    int pulsegt_bit = DU::TrigBits::GetBitIndex("PulseGT");
    int pedestal_bit = DU::TrigBits::GetBitIndex("Pedestal");
    int extasy_bit = DU::TrigBits::GetBitIndex("EXTASY");

    if (((mtcd_trig_word & (1 << pulsegt_bit)) != 0) || ((mtcd_trig_word & (1 << extasy_bit)) != 0)) return true;
    if ((mtcd_trig_word & (1 << pedestal_bit)) != 0) {
      double currentTime = ExtractTimeIndex(ev.GetUniversalTime());
      if (currentTime == -1.0) return true;
      if ((fPedestalCutTime[0] == 0.0) && (fPedestalCutTime[1] == 0.0)) {
        fPedestalCutTime[0] = currentTime - fPedWindow;
        fPedestalCutTime[1] = currentTime + fPedWindow;
      } else{
        if (((currentTime - fPedWindow) < fPedestalCutTime[1]) && (currentTime + fPedWindow > fPedestalCutTime[1])) {//Extend duration of current Pedestal Cut
          fPedestalCutTime[1] = currentTime + fPedWindow;
        } else{//Output previous Pedestal Cut, and start a new Pedestal Cut.
          fPedestalCutTime[2] +=  PushToCuttingRecord(fPedestalCutTime[0], fPedestalCutTime[1], fFoundPedestal,true);
          fPedestalCutTime[0] = currentTime - fPedWindow;
          fPedestalCutTime[1] = currentTime + fPedWindow;
        }
      }
      return true;
    }
    fRateCount[0]++;
    return false;
  }

  bool LivetimeCuts::CorrespondenceCheck(int eventSet[2], std::vector<int>& firstArray, std::vector<int>& secondArray) {

    if ((eventSet[0] == -1) || (eventSet[1] == -1)) return false;

    std::vector<int>::iterator firstArraySearch = find (firstArray.begin(), firstArray.end(), eventSet[0]);
    std::vector<int>::iterator secondArraySearch = find (secondArray.begin(), secondArray.end(), eventSet[1]);

    int firstCheck = (int)(firstArray.end() - firstArraySearch);//The position of input first_event in the first event readout array
    int lastCheck = (int)(secondArray.end() - secondArraySearch);//The position of input last_event in the last event readout array

    if ((firstCheck == 0) || (lastCheck == 0)) {//event not found
      warn << "TPBurstCut::CorrespondenceCheck: Start of Burst or End of Burst Event not found!" << "\n";
      return false;
    } else if (firstCheck != lastCheck) {//event position is not identical, so not correspondent.
       warn << "TPBurstCut::CorrespondenceCheck: Start of Burst and End of Burst Events not Correspondent!" << "\n";
      return false;
    } else {//pass correspondence check.
      firstArray.erase(firstArraySearch);//For optimization, removing good pairs out of the readout
      secondArray.erase(secondArraySearch);//to reduce the search time for further call.
      return true;
    }
  }

  double LivetimeCuts::CalculateSacrifice(int cut, double runLength) {
    double rate = 0.0;
    double cutWindow[3];
    double fractionLT = 0.0;

    cutWindow[0] = fRetriggerWindow;
    cutWindow[1] = fNhitBurstWindow;
    cutWindow[2] = fCAENWindow;

    rate = fRateCount[cut]/runLength;

    if ((cut == 0) || (cut == 2)) {
      fractionLT = 3.0 * pow(rate * cutWindow[cut], 2) / 2.0;
    } else if (cut == 1) {
      fractionLT = 11.0 * pow(rate * cutWindow[cut], 6) / 720.0;
    }

    return (fractionLT * runLength);
  }

} // namespace RAT