#include <RAT/LivetimeCuts.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)
  {
    DS::UniversalTime nullTemp(0,0,0);
    fNullTime = nullTemp;

    fPassFlag = true;
    fFoundRetriggers = false;
    fFoundNhitBurst = false;
    fFoundMuon = false;
    fFoundCLBurst = false;

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

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

    if (fProcessPass == 1) {
      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");

      DS::UniversalTime fStartTime;

      try {
        DBLinkPtr fRun = DB::Get()->GetLink("RUN");
        fRunType = fRun->GetI("runtype");
        fStart_times[0] = fRun->GetI("start_day");//Initialize start time to this value, so it wouldn't
        fStart_times[1] = 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"), fRun->GetD("stop_nsc"));
        DS::UniversalTime startTemp((UInt_t)fRun->GetI("start_day"), (UInt_t)fRun->GetI("start_sec"), fRun->GetD("start_nsc"));
        fValidEndTime = stopTemp;
        fStartTime = startTemp;
      }

      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.;
        fStartTime = fNullTime;
       }

      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 ProduceRunTableProc.cc line 213-241 for detail.
      fInvalidTime = ((fStart_times[0] == -1.) && (fStart_times[1] = -1.) && (fStart_times[2] = -1.));
      std::stringstream stream;
      stream<<"LIVETIME_"<<fRunNumber<<".ratdb";
      fLivetimeTableName = stream.str();

      UInt_t burstWindow = static_cast<unsigned int>(fLBurst->GetI("twin_nhitburst") * 1000);
      DS::UniversalTime nHitBurstWindowTemp(0,0, burstWindow);
      fNhitBurstWindow = nHitBurstWindowTemp;

      UInt_t retriggerWindow = static_cast<unsigned int>(fLRetrigger->GetI("twin_retrigger") * 1000);
      DS::UniversalTime retriggerWindowTemp(0,0, retriggerWindow);
      fRetriggerWindow = retriggerWindowTemp;

      UInt_t muonWindow = static_cast<unsigned int>(fLMuon->GetD("short_follower_time"));
      DS::UniversalTime muonWindowTemp(0, muonWindow, 0);
      fMuonWindow = muonWindowTemp;

      UInt_t CAENWindow = static_cast<unsigned int>(fLCAEN->GetD("twin_caenloss") * 1000);
      DS::UniversalTime CAENWindowTemp(0, 0, CAENWindow);
      fCAENWindow = CAENWindowTemp;

      UInt_t MMFWindow = static_cast<unsigned int>(fLMMF->GetD("time_window"));
      DS::UniversalTime MMFWindowTemp(0, 0, MMFWindow);
      fMissedMuonWindow = MMFWindowTemp;

      UInt_t PedWindow = static_cast<unsigned int>(fLPed->GetD("time_window"));
      DS::UniversalTime PedWindowTemp(0, 0, PedWindow);
      fPedWindow = PedWindowTemp;

      double AtmWindow = fLAtm->GetD("cut_time_window");
      DS::UniversalTime fAtmWindow(0, (UInt_t)(AtmWindow/1E9), (UInt_t)fmod(AtmWindow,1E9));



      //initialize muon parameters to null value.
      fMuonCut = fNullTime;
      fLastMuon = fNullTime;
      fMuonCutTime = fNullTime;
      fMissedMuonCut = fNullTime;
      fMissedMuonCutTime = fNullTime;

      fNhitBurstCutGTID[0] = -1;
      fNhitBurstCutGTID[1] = -1;
      fNhitBurstCutTime[2] = fNullTime;
      fCLBurstCutGTID[0] = -1;
      fCLBurstCutGTID[1] = -1;
      fCLBurstCutTime[2] = fNullTime;

      fPedestalCutTime[0] = fNullTime;
      fPedestalCutTime[1] = fNullTime;
      fPedestalCutTime[2] = fNullTime;

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


    try{//Reading Muon table.
        DBLinkPtr fMuonData = DB::Get()->GetLink("TPMUONFOLLOWER");
        fMuonReadout = fMuonData->GetIArray("gtid_muons");
        std::vector<int>::iterator fakeMuonSearch = find (fMuonReadout.begin(), fMuonReadout.end(), -1);
        if (fakeMuonSearch != fMuonReadout.end()) {
          DS::UniversalTime runStartTemp((UInt_t)fStart_times[0], (UInt_t)fStart_times[1], (UInt_t)fStart_times[2]);
          fCutStart.push_back(runStartTemp);
          fCutEnd.push_back(runStartTemp + fMuonWindow);
        } else {
          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);
                fLastMuon = lastMuonTemp + fMuonWindow;
                detail << "LivetimeCuts::BeginOfRun: Got last muon in previous run." << endl;
           }
         }
         catch(DBNotFoundError &e){
           detail << "LivetimeCuts::BeginOfRun: Could not find LAST_MUON ratdb file." << endl;
         }
        }
        fFoundMuon = true;
      }
      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 = true;
        }
      }
      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);
            DS::UniversalTime fLastAtmTime = fLastAtmTime + fAtmWindow;
            if ((fStartTime != fNullTime) && (fStartTime < fLastAtmTime)) {
              fCutStart.push_back(fStartTime);
              fCutEnd.push_back(fLastAtmTime);
              fMuonCutTime += (fLastAtmTime - fStartTime);
            }
            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]);
            DS::UniversalTime atmCutEnd = atmCutStart + fAtmWindow;
            fCutStart.push_back(atmCutStart);
            fCutEnd.push_back(atmCutEnd);
            fMuonCutTime += fAtmWindow;
          }
        }
      } catch(DBNotFoundError &e){
        warn << "LivetimeCuts::BeginOfRun: Could not find Atmospheric Neutrino ratdb file. Atmospheric Neutrino Cut will not be applied." << endl;
      }

    } else if (fProcessPass == 2) {
      fDCBitCollection.push_back(DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("retriggercut"));
      fDCBitCollection.push_back(DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("tpburstcut"));
      fDCBitCollection.push_back(DU::Utility::Get()->GetDataCleaningBits().GetBitIndex("firstevent"));
    }
  }

  void LivetimeCuts::SetI(const std::string& param, const int value)
  {
    DataCleaningProc::SetI(param,value);
    if (param == "pass") {
      fProcessPass= value;
    }
  }

  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;

    //Set the bit assuming the event passes all cut, if otherwise, will update
    //corresponding bit again later.
    if (fProcessPass == 2) {
      for (size_t iBit = 0; iBit < fDCBitCollection.size(); iBit++) {
        fBit = fDCBitCollection[iBit];
        UpdateMask(ev);
      }
    }

    //Do nothing with orphans or event before first valid events.
    if (ev.GetTrigType() == 0) return OKFALSE;

    if ((fProcessPass == 1) && (ev.GetGTID() == fValidFirstGTID)) {//if the event is the first valid event
      if (fLastMuon != fNullTime) {//if last muon is found
        DS::UniversalTime cutStartTemp((UInt_t)fStart_times[0], (UInt_t)fStart_times[1], (UInt_t)fStart_times[2]);
        if (fLastMuon > cutStartTemp) {//if last muon has muon follower in this run.
          fCutStart.push_back(cutStartTemp);
          fCutEnd.push_back(fLastMuon);
          fMuonCutTime += (fLastMuon - cutStartTemp);
        }
      }
    }

    //if the event is forced event, only update stop time but do nothing.
    if (ForcedEventsTagger(ev)) {
      fStop_times[0] = static_cast<double>(ev.GetUniversalTime().GetDays());
      fStop_times[1] = static_cast<double>(ev.GetUniversalTime().GetSeconds());
      fStop_times[2] = static_cast<double>(ev.GetUniversalTime().GetNanoSeconds());
      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 (fFoundRetriggers) {
      if(RetriggerIdentifier(ev)) {
        fPassFlag = false;
        if (fProcessPass == 2) {
          fBit = fDCBitCollection[0];
          UpdateMask(ev);
        }
      } else {
          if (ev.GetNhits() >= fNThreshold) {
            fRateCount[1]++;
          }
          if (!ev.DigitiserExists()) {
            fRateCount[2]++;
          }
      }
    }

    //Update fStop_times
    fStop_times[0] = static_cast<double>(ev.GetUniversalTime().GetDays());
    fStop_times[1] = static_cast<double>(ev.GetUniversalTime().GetSeconds());
    fStop_times[2] = static_cast<double>(ev.GetUniversalTime().GetNanoSeconds());

    //Call burst identifier for NhitBurst events and update mask if necessary.
    if (fFoundNhitBurst && (BurstIdentifier(ev, fNhitBurstWindow, fNhitBurstCutGTID, fNhitBurstCutTime, fFirstNhitBurstReadout, fLastNhitBurstReadout, fNhitBurstReadout))) {
      fPassFlag = false;
      if (fProcessPass == 2) {
        fBit = fDCBitCollection[1];
        UpdateMask(ev);
      }
    }

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

      //Call muon identifier for TPMuonFollower Cut, no UpdateMask
      if (fFoundMuon && (MuonIdentifier(ev, fMuonReadout, fMuonWindow, fMuonCutTime, fMuonCut))) {
        fPassFlag = false;
      }

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

    return fPassFlag ? OKFALSE : OKTRUE;
  }

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

    if (fProcessPass == 1) {

      DS::UniversalTime startTime((UInt_t)fStart_times[0], (UInt_t)fStart_times[1], (UInt_t)fStart_times[2]);
      DS::UniversalTime stopTime = fValidEndTime;
      RAT::DBTable *LIVETIME_table = new DBTable("LIVETIME");
      LIVETIME_table->SetRunRange(fRunNumber,fRunNumber);
      LIVETIME_table->SetPassNumber(-1);
      LIVETIME_table->SetS("index","");
      LIVETIME_table->SetI("version",1);
      LIVETIME_table->SetI("runtype",fRunType);

      if ((fCutStart.size() != fCutEnd.size()) || (fInvalidTime)) {
        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("pedestal",-1);
        LIVETIME_table->SetD("muon",-1);
        LIVETIME_table->SetD("overlap",-1.);
        LIVETIME_table->SetD("missedmuon",-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 {
        //Push in the last pair of pedestal cut.
        fCutStart.push_back(fPedestalCutTime[0]);
        fCutEnd.push_back(fPedestalCutTime[1]);

        //crop out the cut that exceed end of the run or prior to the start of run.
        for (size_t i = 0; i < fCutEnd.size(); i++) {
          if (fCutEnd[i] < fCutStart[i]) {
            fCutEnd[i] = fCutStart[i];
          }
          if ((fCutEnd[i] > stopTime) && (fCutStart[i] < stopTime)){
            if ((fCutEnd[i] - fCutStart[i]) == fMuonWindow) {
              fMuonCutTime -= (fCutEnd[i] - stopTime);//Crop out overcounting from muon follower cut.
            }
            fCutEnd[i] = stopTime;
          }
          if ((fCutStart[i] < startTime)){
            if (fCutEnd[i] > startTime) {
              fCutStart[i] = startTime;
            } else {
              fCutStart[i] = fCutEnd[i];
            }
          }
        }

        double rawCutTime = 0.0;
        double cutTime = 0.0;
        //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.
        for(size_t i = 0; i < fCutStart.size(); i++) {
          DS::UniversalTime currentStartTime = fCutStart[i];
          DS::UniversalTime currentEndTime = fCutEnd[i];
          rawCutTime += CalculateLivetime(currentStartTime, currentEndTime);
          //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++) {

            //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.
              fCutEnd[j] = fCutStart[j];
            } else if ((currentStartTime >= fCutStart[j]) && (currentEndTime <= fCutEnd[j])) {
              //if target cut fully covers current cut, crush the current cut:
              fCutEnd[i] = fCutStart[i];
            } 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] = currentStartTime;
            }
          }
          cutTime += CalculateLivetime(fCutStart[i], fCutEnd[i]);
        }

        double overlap = rawCutTime - cutTime;
        double rawlt = CalculateLivetime(startTime, stopTime);//raw livetime
        double retrigger = CalculateLivetime(fNullTime, fRetriggerCutTime[2]);
        double burst = CalculateLivetime(fNullTime, fNhitBurstCutTime[2]);
        double clburst = CalculateLivetime(fNullTime, fCLBurstCutTime[2]);
        double muon = CalculateLivetime(fNullTime, fMuonCutTime);
        double ped = CalculateLivetime(fNullTime, fPedestalCutTime[2]);
        double mmf = CalculateLivetime(fNullTime, fMissedMuonCutTime);
        double retrigger_sacrifice = CalculateSacrifice(0, rawlt);
        double burst_sacrifice = CalculateSacrifice(1, rawlt);
        double clburst_sacrifice = CalculateSacrifice(2, rawlt);
        double livetime = rawlt - cutTime;

        LIVETIME_table->SetD("raw_livetime",rawlt);
        LIVETIME_table->SetD("retrigger_cut",retrigger);
        LIVETIME_table->SetD("nhit_burst_cut",burst);
        LIVETIME_table->SetD("caen_loss_burst_cut",clburst);
        LIVETIME_table->SetD("overlap",overlap);
        LIVETIME_table->SetD("muon",muon);
        LIVETIME_table->SetD("pedestal",ped);
        LIVETIME_table->SetD("missedmuon",mmf);
        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();
    fMuonReadout.clear();
    fDCBitCollection.clear();
  }

  double LivetimeCuts::CalculateLivetime(DS::UniversalTime &start, DS::UniversalTime &end) {
    std::vector<double> TimeDiff(3);
    TimeDiff[0] = (double)end.GetDays() - (double)start.GetDays();
    TimeDiff[1] = (double)end.GetSeconds() - (double)start.GetSeconds();
    TimeDiff[2] = end.GetNanoSeconds() - start.GetNanoSeconds();

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

    return timeDifference;
  }

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

    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] = ev.GetUniversalTime();
      }

      burstSearch = find (fLastBurstReadout.begin(), fLastBurstReadout.end(),ev.GetGTID());
      if (burstSearch != fLastBurstReadout.end()){
        fBurstCutGTID[1] = ev.GetGTID();
        fBurstCutTime[1] = 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)) {
          if ((fBurstCutTime[0] - fBurstWindow) < (fBurstCutTime[1] + fBurstWindow)) {
            fCutStart.push_back(fBurstCutTime[0] - fBurstWindow);
            fCutEnd.push_back(fBurstCutTime[1] + fBurstWindow);
            fBurstCutTime[2] += (fCutEnd.back() - fCutStart.back());
          }
        }
        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] = ev.GetUniversalTime();
      if (fProcessPass == 2) {
        fPassFlag = false;
        fBit = fDCBitCollection[2];//Update mask for all the first event
        UpdateMask(ev);
      }
      return false;
    } else if (retriggerLastSearch != fLastRetriggerReadout.end()) {//if input event is last event
      fRetriggerCutGTID[1] = ev.GetGTID();
      fRetriggerCutTime[1] = ev.GetUniversalTime();
      if (fProcessPass == 1) {
        if (CorrespondenceCheck(fRetriggerCutGTID, fFirstRetriggerReadout, fLastRetriggerReadout)) {
          DS::UniversalTime perviousEV((UInt_t)fStop_times[0], (UInt_t)fStop_times[1], (UInt_t)fStop_times[2]);
          if (fRetriggerCutTime[0] < (perviousEV + fRetriggerWindow)) {
            fCutStart.push_back(fRetriggerCutTime[0]);
            fCutEnd.push_back(perviousEV + fRetriggerWindow);
            fRetriggerCutTime[2] += fCutEnd.back() - fCutStart.back();
          }
        }
      }
      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::MuonIdentifier(DS::EV& ev,
                                    std::vector<int>& muonReadout,
                                    DS::UniversalTime& muonWindow,
                                    DS::UniversalTime& muonCutTime,
                                    DS::UniversalTime& muonCut) {

    //Check if the input event is a muon follower of the last muon in previous run
    if (ev.GetUniversalTime() < (fLastMuon + muonWindow)) return true;

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

    if (muonSearch != muonReadout.end()) {//if the input event is muon
      fCutStart.push_back(currentTime);
      fCutEnd.push_back(currentTime + muonWindow);
      //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 != fNullTime) {
      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 = fNullTime;
        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 ext2_bit = DU::TrigBits::GetBitIndex("EXT2");
    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) || ((mtcd_trig_word & (1 << ext2_bit)) != 0)) {
      if (fProcessPass == 1) {
        DS::UniversalTime currentTime = ev.GetUniversalTime();
        if ((fPedestalCutTime[0] == fNullTime) && (fPedestalCutTime[1] == fNullTime)) {
          fPedestalCutTime[0] = currentTime - fPedWindow;
          fPedestalCutTime[1] = currentTime + fPedWindow;
          fPedestalCutTime[2] += fPedWindow + fPedWindow;
        } else{
          if ((currentTime - fPedWindow) < fPedestalCutTime[1]) {//Extend duration of current Pedestal Cut
            if (currentTime + fPedWindow > fPedestalCutTime[1]) {
              fPedestalCutTime[2] += currentTime + fPedWindow - fPedestalCutTime[1];
              fPedestalCutTime[1] = currentTime + fPedWindow;
            }
          } else{//Output previous Pedestal Cut, and start a new Pedestal Cut.
            fCutStart.push_back(fPedestalCutTime[0]);
            fCutEnd.push_back(fPedestalCutTime[1]);
            fPedestalCutTime[0] = currentTime - fPedWindow;
            fPedestalCutTime[1] = currentTime + fPedWindow;
            fPedestalCutTime[2] += fPedWindow + 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] = CalculateLivetime(fNullTime, fRetriggerWindow);
    cutWindow[1] = CalculateLivetime(fNullTime, fNhitBurstWindow);
    cutWindow[2] = CalculateLivetime(fNullTime, 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