// ECAProc.cc
// Contact person: Gabriel Orebi Gann <gorebigann@lbl.gov> - Author and contact
// Contact person: Javier Caravaca <jcaravaca@berkeley.edu> - Contact
// See ECAProc.hh for more details
//———————————————————————//

#include <CLHEP/Units/SystemOfUnits.h>
#include <G4UImessenger.hh>
#include <Randomize.hh>


#include <RAT/ECAProc.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/ChanHWStatus.hh>

#include <TObject.h>

#include <sstream>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;

namespace RAT {
const int ncrate=19;
ECAProc::ECAProc() : Processor("ECAProc")
{

  // Init. variables
  fNumSkip = 0;
  fHitSkip = 0;
  fTinj = 0.;
  fArrayID = -1;
  fNSet = 0;
  //kate
  fTimeSkipped = 0;
  fLastPatternID = 0;

  fNcrate = 19;
  fNcard = 16;
  fNchannel = 32;
  fNcell = 16;
  fNTotCell = fNcrate * fNcard *fNchannel * fNcell;
  fTotChannel = fNcrate * fNcard *fNchannel;
  fNChanNotTested = 0;

  //Javi retro-compatibility
  isSNOdata = false;

  fNTACFitPars = 4;
  fNMaxEV = 10;
  fNDSEvents = 0;
  fBits = new BitManip();

  //Read ECA_Generation file and set flags
  fECAbank = DB::Get()->GetLink("ECAgen");
  fECAflag = fECAbank->GetI("eca_type");
  fECAmerge = fECAbank->GetI("merge_eca");
  fRunMerge = fECAbank->GetI("merge_run");
  fECAdiff = fECAbank->GetI("eca_diff");
  fECAlogverb = fECAbank->GetI("eca_logverbosity");
  fDAQns_offset = fECAbank->GetD("daq_to_ns1");
  fDAQns_max =    fECAbank->GetD("daq_to_ns2");
  fDAQns_scale1 = fECAbank->GetD("daq_to_ns3");
  fDAQns_scale2 = fECAbank->GetD("daq_to_ns4");
  fDoDroop = fECAbank->GetI("eca_dodroop");
  fSuperDebug = (int) fECAbank->GetI("extra_diagnostics");

  // Name output files according to calibration type and run #
  // Check whether those output files already exist and, if so, increment file name
  // Hang on to most recent existing name, for accessing previous values
  fECArootname = "Null.root";
  fECAdebugname = "Null_debug.root";

  //Do not increment pass numbers anymore. This will be handled by the DB
  fPassNo = 0;
  fRunRangeMax = RAT::DB::DBINFINITY;
  fECAIndex = "";

  //Set # TSLP points
  if(fECAbank->GetI("isSNOdata")){
    warn << RED << " Analyzing SNO data "<< CLR << newline;
    isSNOdata = true;
    fNSlopePoints = 31;
  }
  else{
    fNSlopePoints = fECAbank->GetI("NTSlopePoints");
  }
  info << "Setting number of TSlope points to " << fNSlopePoints << newline;

  fNdim = 0;
  fNdim = static_cast<int>(fNSlopePoints / 32) +1;

  // GetECA pattern information
  fECAPatt = DB::Get()->GetLink("ECApatt");
  fNSet = fECAPatt->GetI("NumSets");
  fNPatt = fECAPatt->GetIArray("NumPatt");
  // Load patterns into memory
  ECAProc::LoadPatterns();

  if(fSuperDebug == 1){ //extra histograms for ped or tslopes
    warn << "You are running with extra diagnostics turned on!" << newline;
  }
  if(fSuperDebug == 2){ //plots of tslopes & fits, and a ratdb file w/ the errors on the tslope points
    warn << "You are running with the option to plot all tslopes & fits!" << newline;
    warn << "You can only do this for a few crates." << newline;
    warn << "If you are processing TSlopes from the entire detector, THE PROCESSOR WILL RUN OUT OF MEMORY AND CRASH!" << newline;
  }

  // Initialise arrays for results to 0
  fOverallMask = 0;  // status of ECA run
  // Set bit 30 to 1 to flag job failure
  fOverallMask = fBits->SetBit(fOverallMask, 30);

  // Only init arrays relevant to this calibration run type (save space)
  // AND only init arrays of the previous pass of this run IF we want to merge

  fTestChannel.resize(fTotChannel);

  fCellStatus.resize(fNTotCell);
  for(int n = 0; n<fNTotCell;++n){
    fCellStatus[n]=0;
    if(n<fTotChannel) fTestChannel[n]=0;
  }

  if(fECAflag==1){
    fNevCell.resize(fNTotCell);
    fCellQHS.resize(fNTotCell);
    fCellQHL.resize(fNTotCell);
    fCellQLX.resize(fNTotCell);
    fCellTAC.resize(fNTotCell);
    fCellSigQHS.resize(fNTotCell);
    fCellSigQHL.resize(fNTotCell);
    fCellSigQLX.resize(fNTotCell);
    fCellSigTAC.resize(fNTotCell);
    fCellWidthQHS.resize(fNTotCell);
    fCellWidthQHL.resize(fNTotCell);
    fCellWidthQLX.resize(fNTotCell);
    fCellWidthTAC.resize(fNTotCell);
    fEVQHS.resize(fNTotCell);
    fEVQHL.resize(fNTotCell);
    fEVQLX.resize(fNTotCell);
    fEVTAC.resize(fNTotCell);

    fCellNo.resize(fTotChannel);
    for(int n = 0; n<fTotChannel;++n){
      fCellNo[n].resize(fNcell);
      for(int nc=0; nc<fNcell;++nc){
        fCellNo[n][nc]=0;
      }
    }

    for(int n = 0; n<fNTotCell;++n){
      fNevCell[n] = 0;
      fCellQHS[n]=0.0;
      fCellQHL[n]=0.0;
      fCellQLX[n]=0.0;
      fCellTAC[n]=0.0;
      fCellSigQHS[n]=0.0;
      fCellSigQHL[n]=0.0;
      fCellSigQLX[n]=0.0;
      fCellSigTAC[n]=0.0;
      fCellWidthQHS[n]=0.0;
      fCellWidthQHL[n]=0.0;
      fCellWidthQLX[n]=0.0;
      fCellWidthTAC[n]=0.0;

      fEVQHS[n].resize(fNMaxEV);
      fEVQHL[n].resize(fNMaxEV);
      fEVQLX[n].resize(fNMaxEV);
      fEVTAC[n].resize(fNMaxEV);

      for(int ne = 0; ne<fNMaxEV;++ne){
        fEVQHS[n][ne] = 0;
        fEVQHL[n][ne] = 0;
        fEVQLX[n][ne] = 0;
        fEVTAC[n][ne] = 0;
      }
    }
  }

  if(fECAflag==2){
    fNevTACSlp.resize(fNTotCell);
    fCellBad.resize(fNTotCell);
    fCellSusp.resize(fNTotCell);
    fCellCurl.resize(fNTotCell);
    fPrevBad.resize(fNTotCell);
    for(int n=0;n<fNTotCell;++n){
      fCellBad[n].resize(fNdim);
      fCellSusp[n].resize(fNdim);
      fCellCurl[n].resize(fNdim);
      fPrevBad[n].resize(fNdim);
    }
    fCellTACPoints.resize(fNTotCell);
    fCellTACrms.resize(fNTotCell);
    fCellTACwidth.resize(fNTotCell);
    fCellTACFits.resize(fNTotCell);
    fCellTACFitErrs.resize(fNTotCell); //kate extra output
    fPrevTACPoints.resize(fNTotCell);
    fPrevTACFits.resize(fNTotCell);
    //fLastTAC.resize(fNTotCell);
    fEVTACPoints.resize(fNTotCell);
    fEVTACX.resize(fNSlopePoints);
    for(int n = 0; n<fNTotCell;++n){
      for(int nd=0;nd<fNdim;++nd){
        fCellBad[n][nd]=0;
        fCellSusp[n][nd]=0;
        fCellCurl[n][nd]=0;
      }
      //fLastTAC[n]=0;
      fNevTACSlp[n].resize(fNSlopePoints);
      fCellTACPoints[n].resize(fNSlopePoints);
      fCellTACrms[n].resize(fNSlopePoints);
      fCellTACwidth[n].resize(fNSlopePoints);
      fCellTACFits[n].resize(fNTACFitPars);
      fCellTACFitErrs[n].resize(fNTACFitPars); //kate
      fPrevTACFits[n].resize(fNTACFitPars);
      fEVTACPoints[n].resize(fNSlopePoints);
      fPrevTACPoints[n].resize(fNSlopePoints);
      for(int t = 0; t<fNSlopePoints;++t){
        fNevTACSlp[n][t] = 0;
        fCellTACPoints[n][t] = 0.0;
        fCellTACrms[n][t] = 0.0;
        fCellTACwidth[n][t] = 0.0;
        fEVTACPoints[n][t].resize(fNMaxEV);
        if(n==0)fEVTACX[t] = 0.;
        for(int ne = 0; ne<fNMaxEV;++ne){
          fEVTACPoints[n][t][ne] = 0;
        }
        if(t<fNTACFitPars){
          fCellTACFits[n][t] = 0;
          fCellTACFitErrs[n][t] = 0;//kate
        }
      }
    }
  }

  //kate
  //initialize the "debug" files that contains extra information to aid in figuring out
  //ECA issues in commissioning
  if(fSuperDebug == 1){
    //this is for all 'extra' histograms/plots that aren't tslopes
    //because the tslope files are huge and take forever to open

    //these are the 'debug histos' -- extra info
    fTrigWordHistSkippedEV = new TH1F("fTrigWordHistSkippedEV", "fTrigWordHistSkippedEV", 26, 0., 26.);
    fTrigWordHist = new TH1F("fTrigWordHist", "fTrigWordHist", 26, 0., 26.);

    //this is for tslopes only
    if(fECAflag == 2)
    {
      fResetTimes = new TH1F("fResetTimes", "fResetTimes", 50, -10., 510.);

      //these are to plot which triggers are active at which times
      fb0_NHit100Lo = new TH1F("fb0_NHit100Lo", "fb0_NHit100Lo", 51, 0., 510.);
      fb1_NHit100Med = new TH1F("fb1_NHit100Med", "fb1_NHit100Med", 51, 0., 510.);
      fb2_NHit100Hi = new TH1F("fb2_NHit100Hi", "fb2_NHit100Hi", 51, 0., 510.);
      fb3_NHit20 = new TH1F("fb3_NHit20", "fb3_NHit20", 51, 0., 510.);
      fb4_NHit20LB = new TH1F("fb4_NHit20LB", "fb4_NHit20LB", 51, 0., 510.);
      fb5_ESumLo = new TH1F("fb5_ESumLo", "fb5_ESumLo", 51, 0., 510.);
      fb6_ESumHi = new TH1F("fb6_ESumHi", "fb6_ESumHi", 51, 0., 510.);
      fb7_OWLN = new TH1F("fb7_OWLN", "fb7_OWLN", 51, 0., 510.);
      fb8_OWLELo = new TH1F("fb8_OWLELo", "fb8_OWLELo", 51, 0., 510.);
      fb9_OWLEHi = new TH1F("fb9_OWLEHi", "fb9_OWLEHi", 51, 0., 510.);
      fb10_PulseGT = new TH1F("fb10_PulseGT", "fb10_PulseGT", 51, 0., 510.);
      fb11_Prescale = new TH1F("fb11_Prescale", "fb11_Prescale", 51, 0., 510.);
      fb12_Pedestal = new TH1F("fb12_Pedestal", "fb12_Pedestal", 51, 0., 510.);
      fb13_Pong = new TH1F("fb13_Pong", "fb13_Pong", 51, 0., 510.);
      fb14_Sync = new TH1F("fb14_Sync", "fb14_Sync", 51, 0., 510.);
      fb15_EXTASY = new TH1F("fb15_EXTASY", "fb15_EXTASY", 51, 0., 510.);
      fb16_Ext2 = new TH1F("fb16_Ext2", "fb16_Ext2", 51, 0., 510.);
      fb17_Ext3 = new TH1F("fb17_Ext3", "fb17_Ext3", 51, 0., 510.);
      fb18_Ext4 = new TH1F("fb18_Ext4", "fb18_Ext4", 51, 0., 510.);
      fb19_Ext5 = new TH1F("fb19_Ext5", "fb19_Ext5", 51, 0., 510.);
      fb20_Ext6 = new TH1F("fb20_Ext6", "fb20_Ext6", 51, 0., 510.);
      fb21_Ext7 = new TH1F("fb21_Ext7", "fb21_Ext7", 51, 0., 510.);
      fb22_Ext8PulseAsync = new TH1F("fb22_Ext8PulseAsync", "fb22_Ext8PulseAsync", 51, 0., 510.);
      fb23_SpecialRaw = new TH1F("fb23_SpecialRaw", "fb23_SpecialRaw", 51, 0., 510.);
      fb24_NCD = new TH1F("fb24_NCD", "fb24_NCD", 51, 0., 510.);
      fb25_SoftGT = new TH1F("fb25_SoftGT", "fb25_SoftGT", 51, 0., 510.);
      fb26_MissedTrig = new TH1F("fb26_MissedTrig", "fb26_MissedTrig", 51, 0., 510.);

    }//tslope-specific
  } //fSuperDebug==1

}


ECAProc::~ECAProc()
{
  //Nothing yet
}

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

  std::stringstream eca_pattern_rate;
  eca_pattern_rate << fPattSet <<"_"<<fECArate;
  fECAIndex = eca_pattern_rate.str();

  // IF we are going to want to merge with prev results
  // OR we want to do `diff' tests
  // We need the previous file, so identify it
  // Iff fRunMerge=0, that means use a previous pass of this run
  // If the counter (countpass) >0, then a file already existed for this run, good
  // Otherwise, load the highest pass of run fRunMerge

  // Flag whether we succeed
  bool LoadPrevOK = false;
  int GetPrevious = 0;
  if(fECAmerge==1 || fECAdiff==1) GetPrevious=1;

  if(GetPrevious==1 && fRunMerge==0){  // Load previous results
    DBLinkPtr currentbank;
    if(fECAflag==1){
      currentbank = DB::Get()->GetLink("ECA_PDST",fECAIndex);
      try{
        currentbank->GetS("index");
      }
      catch(DBNotFoundError){
        ECAProc::Die("Previous ECA constants with index "+fECAIndex+" do not exist");
      }
      fPrevStatus = currentbank->GetIArray("pdst_status");
      fPrevQHS = currentbank->GetDArray("pdst_qhs");
      fPrevQHL = currentbank->GetDArray("pdst_qhl");
      fPrevQLX = currentbank->GetDArray("pdst_qlx");
      fPrevTAC = currentbank->GetDArray("pdst_tac");
      if(fECAmerge==1){
        fPrevSigQHS = currentbank->GetDArray("pdst_qhswidth");
        fPrevSigQHL = currentbank->GetDArray("pdst_qhlwidth");
        fPrevSigQLX = currentbank->GetDArray("pdst_qlxwidth");
        fPrevSigTAC = currentbank->GetDArray("pdst_tacwidth");
      }
      fECAPrevDQID = currentbank->GetIArray("pdst_dqid");
      LoadPrevOK = true;
    }
    if(fECAflag==2){
      currentbank = DB::Get()->GetLink("ECA_TSLP",fECAIndex);
      try{
        currentbank->GetS("index");
      }
      catch(DBNotFoundError){
        ECAProc::Die("Previous ECA constants with index "+fECAIndex+" do not exist");
      }
      // CAN ONLY compare to / merge with previous run that has same number TSlope points
      fPrevNSlopePoints = currentbank->GetI("tslp_ntimes");
      if(fNSlopePoints != fPrevNSlopePoints)ECAProc::Die("", 12, 0, 0);

      fPrevBad.resize(fNTotCell);
      for(int n=0;n<fNTotCell;++n){
        fPrevBad[n].resize(fNdim);
      }

      std::vector<double> holdpoints;
      std::vector<double> holdfits;
      std::vector<int> holdbad;
      fPrevEVTACX = currentbank->GetDArray("tslp_times");
      fPrevStatus = currentbank->GetIArray("tslp_status");
      holdbad = currentbank->GetIArray("tslp_bad");
      fECAPrevDQID = currentbank->GetIArray("tslp_dqid");
      holdpoints = currentbank->GetDArray("tslp_points");
      holdfits = currentbank->GetDArray("tslp_fitpars");
      int c1=0;
      int c2=0;
      int c3 = 0;
      for(int n = 0; n<fNTotCell;++n){
        for(int t = 0; t<fNSlopePoints;++t){
          fPrevTACPoints[n][t] = holdpoints[c1];
          ++c1;
        }
        for(int t = 0; t<fNTACFitPars; ++t){
          fPrevTACFits[n][t] = holdfits[c2];
          ++c2;
        }
        for(int nd=0;nd<fNdim;++nd){
          fPrevBad[n][nd] = holdbad[c3];
          ++c3;
        }
      }
      LoadPrevOK = true;
    }
  } // Loaded previous results

  if(GetPrevious==1 && !LoadPrevOK){
    // We wanted the previous results (for merging or diffing)
    // but we don't have any
    ECAProc::Die("", 1, 0,0);
  }

  // Calculate medians
  if(fECAflag==1){
    for(int nn = 0; nn<fNTotCell; ++nn){
      int nch = fBits->GetLCN(nn);

      if(fNevCell[nn]==0){
        fCellQHS[nn] = -9999.0;
        fCellQHL[nn] = -9999.0;
        fCellQLX[nn] = -9999.0;
        fCellTAC[nn] = -9999.0;
        fCellSigQHS[nn] = -9999.0;
        fCellSigQHL[nn] = -9999.0;
        fCellSigQLX[nn] = -9999.0;
        fCellSigTAC[nn] = -9999.0;
        fCellWidthQHS[nn] = -9999.0;
        fCellWidthQHL[nn] = -9999.0;
        fCellWidthQLX[nn] = -9999.0;
        fCellWidthTAC[nn] = -9999.0;
      }

      if(fTestChannel[nch]==0){
        if(nn%16==0){ //Print only for first cell
          fNChanNotTested++;
        }
        continue;
      }

      if(fNevCell[nn]>0){
        ECAProc::MultiCalc(fNevCell[nn], fEVQHS[nn]);
        fCellQHS[nn] = fMedian;
        fCellWidthQHS[nn] = fWidth2;
        fCellSigQHS[nn] = fSig;
        ECAProc::MultiCalc(fNevCell[nn], fEVQHL[nn]);
        fCellQHL[nn] = fMedian;
        fCellWidthQHL[nn] = fWidth2;
        fCellSigQHL[nn] = fSig;
        ECAProc::MultiCalc(fNevCell[nn], fEVQLX[nn]);
        fCellQLX[nn] = fMedian;
        fCellWidthQLX[nn] = fWidth2;
        fCellSigQLX[nn] = fSig;
        ECAProc::MultiCalc(fNevCell[nn], fEVTAC[nn]);
        fCellTAC[nn] = fMedian;
        fCellWidthTAC[nn] = fWidth2;
        fCellSigTAC[nn] = fSig;

      }
    } // !cells

    //run reassignment for invalid fSig values
    ECAProc::ReassignSigs();

    // How did we do for each cell?
    // NB this is done per cell, but the loop is inside the subroutine
    //    because some tests require checks on multiple cells
    ECAProc::PedDiagnostic();

    info <<fNChanNotTested<<" channels not tested " << newline;

  } // ! Pedestal runs (fECAflag==1)

  // Calculate medians at each point on TAC slope, and do fits
  if(fECAflag==2){
    for(int nn = 0; nn<fNTotCell; ++nn){
      int nch = fBits->GetLCN(nn);

      for(int np = 0; np < fNSlopePoints; ++np){
        if(fNevTACSlp[nn][np] == 0){
          fCellTACPoints[nn][np] = -9999.0;
          fCellTACrms[nn][np] = -9999.0;
          fCellTACwidth[nn][np] = -9999.0;
          continue;
        }
      }

      if(fTestChannel[nch]==0){
        if(nn%16==0){
          fNChanNotTested++;
        }
        continue;
      }

      for(int np = 0; np < fNSlopePoints; ++np){
        ECAProc::MultiCalc(fNevTACSlp[nn][np], fEVTACPoints[nn][np]);
        fCellTACPoints[nn][np] = fMedian;
        fCellTACrms[nn][np] = fWidth;
        fCellTACwidth[nn][np] = fWidth2;
      }

      // Filter bad points out of the fits
      fCellBad[nn] = ECAProc::TSlopeFilter(nn);

      // Perform fits to the TAC points
      ECAProc::TSlopeFit(nn);

    } // end loop cells

    // Cell diagnostics
    ECAProc::TSlopeDiagnostic();

    /*// Identify last good point (i.e. not bad or susp) in curl region
    for(int nn = 0; nn<fNTotCell; ++nn){
      fLastTAC[nn] = ECAProc::IdentifyLast(nn);
    } // !cells*/


    // Combine fCellBad and fCellSusp (NOT fCellCurl) words into one (fCellBad)
    for(int nn = 0; nn<fNTotCell; ++nn){
      fCellBad[nn] = ECAProc::TSlopeBadCombined(nn);
    }

    info <<fNChanNotTested<<" channels not tested " << newline;

  }  // Tslopes runs (fECAflag==2)

  // Set bit 30 to 0 to flag job completion
  fOverallMask = fBits->ClearBit(fOverallMask, 30);

  ECAProc::CreateECATable();
  ECAProc::CreateECAMoni();
  ECAProc::CloseECARoot();

  //kate
  if(fSuperDebug ==1)
  {
    ECAProc::WriteECADebugHistosRoot();
    ECAProc::CloseECADebug();
  }
  if(fSuperDebug == 2 )
  {
    ECAProc::FillECADebug();
    ECAProc::WriteTSlopesRoot();
  }

  // Write summary to log file
  ECAProc::Summarise();
  // Close log file
  fECAlogfile.close();


}



void ECAProc::BeginOfRun( DS::Run& run )
{

  fRunNo = run.GetRunID();
  fRunMask = static_cast< UInt_t >( DB::Get()->GetLink("RUN")->GetI("runtype") );

  //Init variables
  fGTDelC = 0; // will be read from ECA header info
  fGTDelF = 0; // will be read from ECA header info
  fPattern = 0; // will be read from ECA header info
  fCalType = 0; // will be read from ECA header info
  fPattSet = 0; // will be read from ECA header info
  fECArate = 0; // will be read from Trigger header info

  //Set filenames
  if(fECAflag==1){ // Doing Pedestal run
    fECAmoniname = dformat("ECA_MONI_%i_%i.pmtdb", fRunNo, fPassNo);
    fECAdbname = dformat("PDST_%i_%i.ratdb", fRunNo, fPassNo);
    fECArootname = dformat("PDSThists_%i_%i.root", fRunNo, fPassNo);
  }
  else if(fECAflag==2){ // Doing Tslope run
    fECAmoniname = dformat("ECA_MONI_%i_%i.pmtdb", fRunNo, fPassNo);
    fECAdbname = dformat("TSLP_%i_%i.ratdb", fRunNo, fPassNo);
    fECArootname = dformat("TSLPhists_%i_%i.root", fRunNo, fPassNo);
  }
  else{
    fECAdbname = "ECA_INVALID.ratdb";
  }

  info<<"ECAProc::ECAProc: Writing ECA results to: "<<fECAlogname<<", "<<fECArootname<<", "<<fECAdbname<<", "<<fECAmoniname<<newline;

  fECAlogname = dformat("ECA_log_%i_%i.txt", fRunNo, fPassNo);
  ECAProc::InitLogOut(fECAlogname);
  ECAProc::InitECARoot(fECArootname);
  if((fSuperDebug==1) || (fSuperDebug==2)){
    ECAProc::InitECADebugOut(fECAdebugdbname);
    ECAProc::WriteECADebugHeader();
  }

  // Initialise output files
  if(fSuperDebug == 1){
    fECAdebughistosname = dformat("ECA_debughistos_%i_%i.root", fRunNo, fPassNo);
    fdebughistosROOTfile = new TFile(fECAdebughistosname.c_str(), "RECREATE");
  }
  //plot tslopes for every cell
  if(fSuperDebug == 2)
  {
    fTSlopesRootName = dformat("ECA_debug_%i_%i.root", fRunNo, fPassNo);
    fTSlopesRootFile = new TFile(fTSlopesRootName.c_str(), "RECREATE");
    //this is where tslope point errors are printed
    fECAdebugdbname = dformat("DBUG_%i_%i.ratdb", fRunNo, fPassNo);
  }


  int GetPrevious = 0;
  if(fECAmerge==1 || fECAdiff==1) GetPrevious=1;

  // Get Run Number
  info << "ECAProc::ECAProc: processing run " << fRunNo << newline;

  // Get DQID info i.e. board IDs and put into a new array for the ECA
  // (We only need 5 of the 6 IDs)
  fDQXXbank = DB::Get()->GetLink("PMT_DQXX");
  warn << RED << " DQXX file: " << fDQXXbank->GetIArray("run_range").at(0) << CLR << newline;
  fDQID = fDQXXbank->GetIArray("dqid");
  if(static_cast<int>(fDQID.size()) != fNcrate*fNcard*6)ECAProc::Die("",3,fDQID.size(),fNcrate*fNcard*6);
  //  G4cout<<dqid.size()<<" DQID entries"<<G4endl;
  int keepval[6] = {1,0,1,1,1,1};
  int counter = 0;
  //int counter2 = 0;
  fECADQID.resize(fNcrate*fNcard*5);
  if(GetPrevious==1)fECAPrevDQID.resize(fNcrate*fNcard*5);
  for(int ncr=0;ncr<fNcrate;++ncr){
    for(int nca=0;nca<fNcard;++nca){
      for(int i=0;i<6;++i){
        if(keepval[i]==1){
          int offset = 96*ncr + 6*nca + i;
          //if(offset!=counter2)G4cout<<"BADBADBADBAD"<<G4endl;
          fECADQID[counter] = fDQID[offset];
          if(GetPrevious==1)fECAPrevDQID[counter] = 0;
          ++counter;
        }
        //++counter2;
      }
    }
  }
  fDQID.resize(0);

}

// ID last good point on TAC slope
int ECAProc::IdentifyLast(int cccc)
{
  std::vector<int> bad = fCellBad[cccc];
  std::vector<int> susp = fCellSusp[cccc];
  std::vector<int> curl = fCellCurl[cccc];
  int value = 0;

  // Find last good point
  int i = fNSlopePoints - 1;
  while(i>-1 && (fBits->TestBit(bad,i) || fBits->TestBit(susp,i))){
    i--;
  }

  // No good points?
  if(i==-1){
    value = 999;
    return value;
  }

  // If we're still in the curl region, this is our guy
  if(fBits->TestBit(curl,i)){
    value = i;
    return value;
  }

  // Last good point was outside curl region, so this cell has no curl
  // Set value to (*hmm, ponder*) something useful... has to be > any possible Npoints
  value = 998;
  return value;

}

// Combine points flagged as either bad or suspicious into one word
std::vector<int> ECAProc::TSlopeBadCombined(int cccc)
{
  std::vector<int> word;
  word.resize(fNdim);
  for(int id=0;id<fNdim;++id){
    word[id] = 0;
  }
  std::vector<int> bad = fCellBad[cccc];
  std::vector<int> susp = fCellSusp[cccc];
  //int curl = fCellCurl[cccc];
  for(int i=0;i<fNSlopePoints;++i){
    if(fBits->TestBit(bad,i) || fBits->TestBit(susp,i)){// || fBits->TestBit(curl,i)){
      word = fBits->SetBit(word,i);
    }
  }
  return word;
}

void ECAProc::PedDiagnostic()
{

  std::vector<double> qhs_valid = fECAbank->GetDArray("qhs_valid");
  std::vector<double> qhl_valid = fECAbank->GetDArray("qhl_valid");
  std::vector<double> qlx_valid = fECAbank->GetDArray("qlx_valid");
  std::vector<double> tac_valid = fECAbank->GetDArray("tac_valid");
  int MaxN0ev  = fECAbank->GetI("PedMaxN0ev");
  int MaxN160ev = fECAbank->GetI("PedMaxN160ev");
  int MaxNboardID = fECAbank->GetI("PedMaxNboardID");
  std::vector<int> MaxNFlag  = fECAbank->GetIArray("PedMaxNFlag");
  std::vector<double> MinAvePed = fECAbank->GetDArray("PedMinAvePed"); //67,69,71,73
  std::vector<double> MaxAvePed = fECAbank->GetDArray("PedMaxAvePed"); //68,70,72,74
  std::vector<double> MaxAveWidth = fECAbank->GetDArray("PedMaxAveWidth"); //75-78
  std::vector<double> MaxAveDiff = fECAbank->GetDArray("PedMaxAveDiff");   //79-82
  int MaxNChan = fECAbank->GetI("PedMaxNChan");   //83
  int MaxNfSigChan = fECAbank->GetI("PedMaxNfSigChan");   //83
  const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();

  float testvals[4][4];  // [quantity][test]
  // Negate value for lower median, so all the tests can be `<'
  float test[4] = {-1.0, 1.0, 1.0, 1.0};
  for(int i=0;i<4;++i){
    testvals[0][i] = test[i]*qhs_valid[i];
    testvals[1][i] = test[i]*qhl_valid[i];
    testvals[2][i] = test[i]*qlx_valid[i];
    testvals[3][i] = test[i]*tac_valid[i];
  }

  int NBadCh = 0;
  int NBadChCa = 0;
  int NBadSeq = 0;
  int NBadMB = 0;
  int NBadDB = 0;
  int NEmptyCh = 0;
  int NPartCh = 0;
  int NBadPedCh[4][5];  //[quantity][test]
  for(int i=0;i<4;++i){
    for(int j=0;j<5;++j){
      NBadPedCh[i][j] = 0;
    }
  }
  int NfSigUsedCh = 0; //kate

  float runningsum[5][5]; // [quantity][test] eg [qhs][median]
  int nrunsum[5][5];
  for(int i=0;i<5;++i){
    for(int j=0;j<5;++j){
      runningsum[i][j] = 0.0;
      nrunsum[i][j] = 0;
    }
  }

  // Per crate sums/averages
  float average[ncrate][4][4];  // [quant][ped,width,diff,absdiff]
  int nave[ncrate][4][4];
  int NCrBadPed[ncrate][5][4];
  int *NCrSeq = new int[fNcrate];
  int *NCrMB = new int[fNcrate];
  int *NCrDB = new int[fNcrate];
  int *NCrNev = new int[fNcrate];
  int *NCr0ev = new int[fNcrate];
  int *NCrTot = new int[fNcrate];
  for(int ic=0;ic<fNcrate;++ic){
    NCrSeq[ic] = 0;
    NCrMB[ic] = 0;
    NCrDB[ic] = 0;
    NCrNev[ic] = 0;
    NCr0ev[ic] = 0;
    NCrTot[ic] = 0;
    for(int j=0;j<5;++j){
      for(int k=0;k<4;++k){
        NCrBadPed[ic][j][k] = 0;
        if(j<4)average[ic][j][k] = 0.0;
        if(j<4)nave[ic][j][k] = 0;
      }
    }
  }

  //flag cells using default fSig values
  for(int i=0; i < (int) fCellsToReassign.size(); i++){
    int icell = fCellsToReassign[i];
    if(fTestChannel[icell]==0){continue;}
    fCellStatus[icell] = fBits->SetBit(fCellStatus[icell],28); //cell flag
    int cccc0 = fBits->GetCCCC(icell,0);
    fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],29); //channel flag
  }

  for(int nch = 0; nch<fTotChannel; ++nch){
    int cccc0 = fBits->GetCCCC(nch,0);
    // We can get the crate/card/channel for cell 0
    // NB ichan is the channel number on the card i.e. [0,31]
    //    nch is the overall channel number i.e. [0,fTotChannel]
    int icrate = fBits->GetCratewCell(cccc0);
    int icard = fBits->GetCardwCell(cccc0);
    int ichan = fBits->GetChannelwCell(cccc0);
    //G4cout<<"CCC = "<<icrate<<" "<<icard<<" "<<ichan<<G4endl;

    bool dodiff = false;
    if(fECAdiff==1) dodiff = true;

    int NBadCe = 0;
    int NEmptyCe = 0;
    int NPartCe = 0;
    //int NTotBadPed[20], NTotBadPedCh[


    for(int nce = 0; nce<fNcell;++nce){
      int cccc = fBits->GetCCCC(nch,nce);

      // fECADQID[i]:       i = 80*iCrate + 5*iCard + iElement
      //   iElement = 0 : MB ID for card
      //              1 : DC ID for DC0
      //              2 : DC ID for DC1
      //              3 : DC ID for DC2
      //              4 : DC ID for DC3

      if(fECAdiff==1){
        // 1.a. Test MB ID
        int offset = 80*icrate + 5*icard;
        if(fECADQID[offset]!=fECAPrevDQID[offset]){
          dodiff = false;  // Don't do diff tests because the MB has changed
          fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],16);
        }

        // 1.b. Test DB ID
        int daughterID = (ichan/8) + 1;
        offset = offset + daughterID;
        if(fECADQID[offset]!=fECAPrevDQID[offset]){
          dodiff = false;  // Don't do diff tests because the DB has changed
          fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],17);
        }

        // 2. Test if channel got good data
        if( (fPrevStatus[cccc] & (1<<13) ) ||  // Cell with less than 10 events
            (fPrevStatus[cccc] & (1<<14) ) ||  // Cell with 0 events
            (fPrevStatus[cccc] & (1<<25) )  )  // Channel with no data
          dodiff = false; //Do not perform diff test
      }

      // 3. Test if channel is getting data
      bool seq = channelHardwareStatus.IsECAChannelActive(nch);
      if(seq){

        // 4. Test # events
        if(nce==0)fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],27);
        // Fill hist fNCeEv
        fHists.fNCeEv->Fill(fNevCell[cccc]);

        if(fNevCell[cccc]<10 && fNevCell[cccc]!=0){
          fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],13);
          fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],26);
        }
        if(fNevCell[cccc]>0){
          fCellStatus[cccc0] = fBits->ClearBit(fCellStatus[cccc0],27);

          // 5. Pedestal / width / difference tests
          float currvals[4][4]; // [quantity][test]
          // Passes test if currvals[][] < testvals[][] (hence negate first value)
          currvals[0][0] = -1.0*fCellQHS[cccc];
          currvals[1][0] = -1.0*fCellQHL[cccc];
          currvals[2][0] = -1.0*fCellQLX[cccc];
          currvals[3][0] = -1.0*fCellTAC[cccc];
          currvals[0][1] = fCellQHS[cccc];
          currvals[1][1] = fCellQHL[cccc];
          currvals[2][1] = fCellQLX[cccc];
          currvals[3][1] = fCellTAC[cccc];
          currvals[0][2] = fCellWidthQHS[cccc];
          currvals[1][2] = fCellWidthQHL[cccc];
          currvals[2][2] = fCellWidthQLX[cccc];
          currvals[3][2] = fCellWidthTAC[cccc];
          currvals[0][3] = 0.0;
          currvals[1][3] = 0.0;
          currvals[2][3] = 0.0;
          currvals[3][3] = 0.0;
          if(dodiff){
            currvals[0][3] = (fCellQHS[cccc] - fPrevQHS[cccc]);
            currvals[1][3] = (fCellQHL[cccc] - fPrevQHL[cccc]);
            currvals[2][3] = (fCellQLX[cccc] - fPrevQLX[cccc]);
            currvals[3][3] = (fCellTAC[cccc] - fPrevTAC[cccc]);
          }

          for(int nq=0; nq<4; ++nq){
            // Fill hists
            fHists.fPedestal[nq]->Fill(currvals[nq][1]);
            fHists.fPWidth[nq]->Fill(currvals[nq][2]);
            fHists.fDiff[nq]->Fill(currvals[nq][3]);
            // Droop plots
            if(fDoDroop==1){
              for(int nev=0;nev<fNevCell[cccc];++nev){
                if(nev>9)continue;
                if(nq==0){
                  fHists.fDroop[nq][nev]->Fill(currvals[nq][1]-fEVQHS[cccc][nev]);
                }
                if(nq==1){
                  fHists.fDroop[nq][nev]->Fill(currvals[nq][1]-fEVQHL[cccc][nev]);
                }
                if(nq==2){
                  fHists.fDroop[nq][nev]->Fill(currvals[nq][1]-fEVQLX[cccc][nev]);
                }
                if(nq==3){
                  fHists.fDroop[nq][nev]->Fill(currvals[nq][1]-fEVTAC[cccc][nev]);
                }
              }
              // Test channel
              // leave out crate 17, card 15 -> This is the FECD
              if(!(icrate==17&&icard==15)){
                for(int nev=0;nev<fNevCell[cccc];++nev){
                  if(nev>9)continue;
                  if(nq==0){
                    fHists.fDroopValid[nq][nev]->Fill(currvals[nq][1]-fEVQHS[cccc][nev]);
                  }
                  if(nq==1){
                    fHists.fDroopValid[nq][nev]->Fill(currvals[nq][1]-fEVQHL[cccc][nev]);
                  }
                  if(nq==2){
                    fHists.fDroopValid[nq][nev]->Fill(currvals[nq][1]-fEVQLX[cccc][nev]);
                  }
                  if(nq==3){
                    fHists.fDroopValid[nq][nev]->Fill(currvals[nq][1]-fEVTAC[cccc][nev]);
                  }
                }
                if(nq==3){  // Do cell ordering for TAC only
                  for(int ncell=0;ncell<fNcell;++ncell){
                    if(fCellNo[nch][ncell]==nce){
                      fHists.fDroopByCell[ncell]->Fill(currvals[nq][1]-fEVTAC[cccc][0]);
                    }
                  }
                }
              }
            } //If DoDroop

            // Increment crate-averages
            for(int nt=0;nt<3;++nt){
              average[icrate][nq][nt] += currvals[nq][nt+1];
              nave[icrate][nq][nt] += 1;
            }
            average[icrate][nq][3] += fabs(currvals[nq][3]);
            nave[icrate][nq][3] += 1;

            // Do Tests and increment run-averages
            for(int nt=0; nt<4; ++nt){ // PedMin, PedMax, Width, Diff
              //Do not include bad cells in the run average
              if(nt==3){
                int bitid = 1 + 3*nq;
                if(!fBits->TestBit(fCellStatus[cccc],bitid) &&
                   !fBits->TestBit(fPrevStatus[cccc],bitid)){
                  runningsum[nq][nt] += test[nt]*fabs(currvals[nq][nt]);
                  nrunsum[nq][nt] += 1;
                  runningsum[nq][4] += currvals[nq][3];
                }
              }
              else{
                runningsum[nq][nt] += test[nt]*fabs(currvals[nq][nt]);
                nrunsum[nq][nt] += 1;
              }
              if(currvals[nq][nt]>testvals[nq][nt]){
                // Failed the test!  Set appropriate bits.
                // Quant nq is bits nq*3-2 to nq*3 (tests 0 & 1 are the first bit)
                int ntest = nt;
                if(ntest>0)ntest -= 1;
                // Quant nq, test `ntest', is bit nq*3-2+ntest
                int bitid = (nq+1)*3 -2 +ntest;
                fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],bitid);
                // Also set bits in cell0 mask for the channel
                // Quant nq is bit 18+nq
                // Test ntest is bit 22+ntest
                bitid = 18+nq;
                fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],bitid);
                bitid = 22+ntest;
                fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],bitid);
              }
            }
          }

          //Javi
          fHists.fQHSvsQHL->Fill(currvals[0][1],currvals[1][1]);

        } // ! fNevCell >0
        else{ //fNevCell[cccc]==0
          fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],14);
          // Medians should already be set to -9999.0
          if(fCellQHS[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellQHL[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellQLX[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellTAC[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellWidthQHS[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellWidthQHL[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellWidthQLX[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }
          if(fCellWidthTAC[cccc] != -9999.0){
            ECAProc::Die("",2,0,0);
          }

        } // ! test # events

      } // getting data
      else{ // not getting data
        fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],25);
      } // ! test of ECAChannel active



    } // !cells for this channel

    // Check first cell of this channel to determine `channel status' bits
    int channel0 = fBits->GetBits(fCellStatus[cccc0], 16, 12);
    // If any were bad, set the overall channel status bit for cell0
    //Also need to include bit 29 in this check
    if(( channel0 != 0 ) || ( fBits->TestBit(fCellStatus[cccc0],29) ))fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],15);

    for(int nce = 0; nce<fNcell;++nce){
      int cccc = fBits->GetCCCC(nch,nce);
      int cellstat = fBits->GetBits(fCellStatus[cccc], 1, 14);
      // If any tests were bad, set the overall cell status bit
      if(( cellstat != 0 ) || ( fBits->TestBit(fCellStatus[cccc],28) )){
        fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],0);
        ++NBadCe;

        if(fBits->TestBit(fCellStatus[cccc],13))++NPartCe;
        if(fBits->TestBit(fCellStatus[cccc],14))++NEmptyCe;
        // Counters for number of cells failing each test for each quant?
      }
      // Copy channel status from cell0
      if(nce>0){
        for(int nbit=15;nbit<28;++nbit){
          if(fBits->TestBit(fCellStatus[cccc0],nbit)){
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],nbit);
          }
        }
        //also need to get bit 29
        if(fBits->TestBit(fCellStatus[cccc0],29)){
          fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],29);
        }
      }
    }

    // Increment some channel counters
    int test01 = fBits->GetBits(fCellStatus[cccc0], 18, 7);
    int test02 = fBits->GetBits(fCellStatus[cccc0], 26, 2);
    if(test01!=0 || test02!=0)++NBadChCa;
    if(fBits->TestBit(fCellStatus[cccc0],15)){
      ++NBadCh;
      ++NCrTot[icrate];
      if(fBits->TestBit(fCellStatus[cccc0],26)){++NPartCh;  ++NCrNev[icrate];}
      if(fBits->TestBit(fCellStatus[cccc0],27)){++NEmptyCh; ++NCr0ev[icrate];}
      if(fBits->TestBit(fCellStatus[cccc0],25)){++NBadSeq;  ++NCrSeq[icrate];}
      if(fBits->TestBit(fCellStatus[cccc0],16)){++NBadMB;   ++NCrMB[icrate];}
      if(fBits->TestBit(fCellStatus[cccc0],17)){++NBadDB;   ++NCrDB[icrate];}
      if(fBits->TestBit(fCellStatus[cccc0],29)) ++NfSigUsedCh;
      bool qhsORt = false;
      for(int nq=0; nq<4; ++nq){
        if(fBits->TestBit(fCellStatus[cccc0],18+nq)){
          ++NBadPedCh[nq][3];
          bool donechped = false;
          for(int ntest=0;ntest<3;++ntest){
            bool donechtest = false;
            bool qhsORt2 = false;
            for(int nce=0; nce<fNcell;++nce){
              int cccc = fBits->GetCCCC(nch,nce);
              int bitid = (nq+1)*3 -2 +ntest;
              if(!donechtest && fBits->TestBit(fCellStatus[cccc],bitid)){
                if(!donechped){
                  donechped = true;
                  ++NBadPedCh[nq][4];
                  ++NCrBadPed[icrate][nq][0];  // This channel failed any test for quant nq
                  if((nq==0||nq==3)&&!qhsORt){  // Not yet failed for QHS+TAC combo
                    qhsORt = true;
                    ++NCrBadPed[icrate][4][0];   // Failed either QHS/TAC on any test
                  }
                }
                ++NBadPedCh[nq][ntest];
                ++NCrBadPed[icrate][nq][ntest+1];  // This channel failed this test for quant nq
                donechtest = true;
                if((nq==0||nq==3)&&!qhsORt2){
                  qhsORt2 = true;
                  ++NCrBadPed[icrate][4][ntest+1]; // Failed QHS/TAC on this test
                }
              }
            }
          }
        }
        if(NBadPedCh[nq][3]!=NBadPedCh[nq][4]){
          ECAProc::Die("",4,NBadPedCh[nq][3],NBadPedCh[nq][4]);
        }
      }
    }

  }   // !channels

  // Loop over crates and fill hists
  for(int ic=0; ic<fNcrate; ++ic){
    fHists.fNBadNev->Fill(NCrNev[ic]);
    fHists.fNBad0ev->Fill(NCr0ev[ic]);
    fHists.fNBadSeq->Fill(NCrSeq[ic]);
    fHists.fNBadMB->Fill(NCrMB[ic]);
    fHists.fNBadDB->Fill(NCrDB[ic]);
    fHists.fNBad->Fill(NCrTot[ic]);

    for(int ip=0; ip<4; ++ip){
      for(int it=0; it<4; ++it){
        if((nave[ic][ip][it])>0.){average[ic][ip][it] = average[ic][ip][it] / nave[ic][ip][it];}
      }
      fHists.fPedAve[ip]->Fill(average[ic][ip][0]);
      fHists.fWidthAve[ip]->Fill(average[ic][ip][1]);
      fHists.fDiffAve[ip]->Fill(average[ic][ip][2]);
      fHists.fAbsDiffAve[ip]->Fill(average[ic][ip][3]);
    }
    for(int ip=0; ip<5; ++ip){
      for(int it=0; it<4; ++it){
        fHists.fNBadPed[ip][it]->Fill(NCrBadPed[ic][ip][it]);
      }
    }
  }

  // NOW WE DO THE FULL RUN STATUS CHECKS

  // Calculate run averages
  // (NB nt=0,1 should have the same averages: both ave medians)
  for(int nq=0; nq<4; ++nq){
    for(int nt=0; nt<4; ++nt){
      if(nrunsum[nq][nt]>0){
        runningsum[nq][nt] = runningsum[nq][nt]/nrunsum[nq][nt];
        if(nt==3)runningsum[nq][4] = runningsum[nq][4]/nrunsum[nq][nt];
      }
      else{
        runningsum[nq][nt] = 0.0;
        if(nt==3)runningsum[nq][4] = 0.0;
      }

      // Test medians (bits 7-10)
      if(nt==1){
        if((runningsum[nq][nt]<MinAvePed[nq])||(runningsum[nq][nt]>MaxAvePed[nq])){
          fOverallMask = fBits->SetBit(fOverallMask, 7+nq);
        }
      }

      // Test widths (bits 11-14)
      if(nt==2){
        if(runningsum[nq][nt]>MaxAveWidth[nq]){
          fOverallMask = fBits->SetBit(fOverallMask, 11+nq);
        }
      }

      // Test abs diffs (bits 15-18)
      if(nt==3){
        if(runningsum[nq][nt]>MaxAveDiff[nq]){
          fOverallMask = fBits->SetBit(fOverallMask, 15+nq);
        }
      }

    }
  }

  // Tests on overall run stats:

  // 0. Test UC bit
  if(fBits->TestBit(fRunMask, 21)){
    fOverallMask = fBits->SetBit(fOverallMask, 26);
  }

  // 1. Number of empty channels
  if(NEmptyCh > MaxN0ev){
    fOverallMask = fBits->SetBit(fOverallMask, 0);
  }

  // 2. Number of partial channels (i.e. <160 events)
  if(NPartCh > MaxN160ev){
    fOverallMask = fBits->SetBit(fOverallMask, 1);
  }

  // 3. Number of board ID failures (DB+MB)
  if((NBadDB+NBadMB) > MaxNboardID){
    fOverallMask = fBits->SetBit(fOverallMask, 2);
  }

  // 4. Number of flagged channels in total
  if(NBadCh > MaxNChan){
    fOverallMask = fBits->SetBit(fOverallMask, 19);
  }

  // 5. Number of channels failing the tests for QHS/QHL/QLX/TAC
  for(int nq=0;nq<4;++nq){
    if(NBadPedCh[nq][3] > MaxNFlag[nq]){
      fOverallMask = fBits->SetBit(fOverallMask, 3+nq);
    }
  }

  // 6. Number of channels using fSig defaults (i.e. Nev<4) is greater than max
  if(NfSigUsedCh > MaxNfSigChan){
    fOverallMask = fBits->SetBit(fOverallMask, 20);
  }

  delete[] NCrSeq;
  delete[] NCrMB;
  delete[] NCrDB;
  delete[] NCrNev;
  delete[] NCr0ev;
  delete[] NCrTot;

}

void ECAProc::TSlopeDiagnostic()
{
  std::vector<double> cubic_min = fECAbank->GetDArray("cubic_min");
  std::vector<double> cubic_max = fECAbank->GetDArray("cubic_max");
  std::vector<double> cubic_diff = fECAbank->GetDArray("cubic_diff");
  int cubic_minNfit = fECAbank->GetI("cubic_minNfit");
  float width_TStep = fECAbank->GetD("width_TStep");
  int width_maxN    = fECAbank->GetI("width_maxN");
  int Nev_TStep    = fECAbank->GetI("Nev_TStep");
  int Nev_maxN     = fECAbank->GetI("Nev_maxN");
  int susp_maxN     = fECAbank->GetI("susp_maxN");
  int NRailedPoints  = fECAbank->GetI("NRailedPoints");
  //float maxchisq    = fECAbank->GetD("MaxChisq");
  //float resid_TStep = fECAbank->GetD("Resid_TStep");
  int TslpMaxNboardID  = fECAbank->GetI("TslpMaxNboardID");
  int TslpMaxN0ev  = fECAbank->GetI("TslpMaxN0ev");
  int TslpMaxN0evTStep  = fECAbank->GetI("TslpMaxN0evTStep");
  int TslpMaxNWidth  = fECAbank->GetI("TslpMaxNWidth");
  int TslpMaxNFit  = fECAbank->GetI("TslpMaxNFit");
  int TslpMaxNSusp  = fECAbank->GetI("TslpMaxNSusp");
  int TslpMaxNTAC  = fECAbank->GetI("TslpMaxNTAC");
  std::vector<int> TslpMaxNCubic  = fECAbank->GetIArray("TslpMaxNCubic");
  std::vector<double> TslpMinAveCub  = fECAbank->GetDArray("TslpMinAveCub");
  std::vector<double> TslpMaxAveCub  = fECAbank->GetDArray("TslpMaxAveCub");
  std::vector<double> TslpMaxAveCubDiff  = fECAbank->GetDArray("TslpMaxAveCubDiff");
  int TslpMaxNChan  = fECAbank->GetI("TslpMaxNChan");
  int TslpMaxNRailedChan = fECAbank->GetI("TslpMaxNRailedChan");
  int TslpMaxNSuspPostResetHitsChan = fECAbank->GetI("TslpMaxNSuspPostResetHitsChan");
  int TslpMaxNFirstPointBadChan = fECAbank->GetI("TslpMaxNFirstPointBadChan");

  const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();

  int NBadCh = 0;
  int NEmptyCh = 0;
  int NPartCh = 0;
  int NFratiCh = 0;
  int NBadUsedCh = 0;
  int NBadWidthCh = 0;
  int NSuspCh = 0;
  int NBadFitCh[4][3]; // test / parameter
  int NBadSeq = 0;
  int NBadMB = 0;
  int NBadDB = 0;
  int NRailedCh = 0; //kate
  int NFirstPointBadCh = 0; //kate
  int NSuspPostResetHitsCh = 0; //kate

  float runningsum[3][3];
  int nrunsum[2][3];

  for(int nfp=0;nfp<3;++nfp){
    for(int nt=0;nt<2;++nt){
      NBadFitCh[nt][nfp] = 0;
      runningsum[nt][nfp] = 0.0;
      nrunsum[nt][nfp] = 0;
    }
    NBadFitCh[2][nfp] = 0;
    NBadFitCh[3][nfp] = 0;
    runningsum[2][nfp] = 0.0;
  }

  // Per crate values
  float Crrunningsum[ncrate][3][3];
  int Crnrunsum[ncrate][3][3];
  int NCrBadFitCh[ncrate][2][3];

  int *NCrBadUsedCh = new int[fNcrate];
  int *NCrBadWidthCh = new int[fNcrate];
  int *NCrPartCh = new int[fNcrate];
  int *NCrEmptyCh = new int[fNcrate];
  int *NCrBadSeq = new int[fNcrate];
  int *NCrBadMB = new int[fNcrate];
  int *NCrBadDB = new int[fNcrate];
  int *NCrSuspCh = new int[fNcrate];
  int *NCrBadCh = new int[fNcrate];

  for(int ic=0;ic<fNcrate;++ic){
    NCrBadUsedCh[ic] = 0;
    NCrBadWidthCh[ic] = 0;
    NCrPartCh[ic] = 0;
    NCrEmptyCh[ic] = 0;
    NCrBadSeq[ic] = 0;
    NCrBadMB[ic] = 0;
    NCrBadDB[ic] = 0;
    NCrSuspCh[ic] = 0;
    NCrBadCh[ic] = 0;
    for(int nfp=0;nfp<3;++nfp){
      for(int nt=0;nt<3;++nt){
        if(nt<2)NCrBadFitCh[ic][nt][nfp] = 0;
        Crrunningsum[ic][nt][nfp] = 0.0;
        Crnrunsum[ic][nt][nfp] = 0;
      }
    }
  }

  //float tresid; // per cell, per point XXXXXXX ?? Frati thing

  for(int nch = 0; nch<fTotChannel; ++nch){
    int cccc0 = fBits->GetCCCC(nch,0);
    // We can get the crate/card/channel for cell 0
    // NB ichan is the channel number on the card i.e. [0,31]
    //    nch is the overall channel number i.e. [0,fTotChannel]
    int icrate = fBits->GetCratewCell(cccc0);
    int icard = fBits->GetCardwCell(cccc0);
    int ichan = fBits->GetChannelwCell(cccc0);

    bool dodiff = false;
    if(fECAdiff==1){
      dodiff = true;
      if(!ECAProc::MatchPrevTinj("diff"))dodiff = false;
    }

    int NBadCe = 0;
    int NEmptyCe = 0;
    int NPartCe = 0;
    //int NFrati = 0;
    int NBadUsed = 0;
    int NBadWidth = 0;
    int NSusp = 0;
    int NBadFit[3][3];
    for(int nfp=0;nfp<3;++nfp){
      for(int nt=0;nt<3;++nt){
        NBadFit[nt][nfp] = 0;
      }
    }

    for(int nce = 0; nce<fNcell;++nce){
      int cccc = fBits->GetCCCC(nch,nce);

      //trying to get robust criteria to skip post-reset points
      //ideally, we would just read in GT valid for each cell and
      //skip points with time after GT valid, but this is not available yet
      int first9999index = -1;
      for(int np=0; np<fNSlopePoints; ++np){
         if(first9999index == -1){
            if(fCellTACPoints[cccc][np] == -9999.0){
               first9999index = np;
            }
         }
         else break;
      }
      int second9999index = -1;
      for(int np=first9999index+1; np<fNSlopePoints; ++np){
         if(second9999index == -1){
            if(fCellTACPoints[cccc][np] == -9999.0){
               second9999index = np;
            }
         }
         else break;
      }

      //this is for the 'first point is bad' flag
      bool firstpointis9999 = 0;
      if(first9999index == 0) firstpointis9999 = 1;

      //this is for the 'post-reset hits' not consistent with noise flag
      int numpostresethits = 0;
      if(!firstpointis9999){
        for(int np=first9999index; np<fNSlopePoints; ++np){
          if(first9999index == -1) break;
          if(fCellTACPoints[cccc][np] != -9999.0){
            numpostresethits++;
          }
        }
      }
      else{
        for(int np=second9999index; np<fNSlopePoints; ++np){
          if(second9999index == -1) break;
          if(fCellTACPoints[cccc][np] != -9999.0){
            numpostresethits++;
          }
        }
      }

      if(fECAdiff==1){
        // 1.a. Test MB ID
        int offset = 80*icrate + 5*icard;
        if(fECADQID[offset]!=fECAPrevDQID[offset]){
          dodiff = false;  // Don't do diff tests because the MB has changed
          fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],24);
        }

        // 1.b. Test DB ID
        int daughterID = (ichan/8) + 1;
        offset = offset + daughterID;
        if(fECADQID[offset]!=fECAPrevDQID[offset]){
          dodiff = false;  // Don't do diff tests because the DB has changed
          fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],25);
        }

        // 2. Test if channel got good data
        if( (fPrevStatus[cccc] & (1<<1) ) ||   // Cell with 0 events
            (fPrevStatus[cccc] & (1<<26) ) ||  // Cell no recieving data
            (fPrevStatus[cccc] & (1<<11) )  )  // Cell with too few events
          dodiff = false; //Do not perform diff test

      }

      // 3. Test if channel is getting data
      bool seq = channelHardwareStatus.IsECAChannelActive(nch);
      if(seq){

        // 4. Test # events
        if(nce==0)fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],14);
        fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],1);
        int ncetoofew = 0;
        for(int np=0; np<fNSlopePoints; ++np){
          // Fill hist fNEvPerTStep
          fHists.fNEvPerTStep->Fill(fNevTACSlp[cccc][np]);

          if(fNevTACSlp[cccc][np]>0){
            fCellStatus[cccc] = fBits->ClearBit(fCellStatus[cccc],1);
          }
          //want to check only if we're not in the post-reset region
          //this is loose criteria but good enough for now
          //ideally / eventually read in GT valid here instead
          if( np < first9999index && !firstpointis9999 ) {
            if(fNevTACSlp[cccc][np]<Nev_TStep){
              ++ncetoofew;
              fCellSusp[cccc] = fBits->SetBit(fCellSusp[cccc],np);
            }
          }
          if( np < second9999index && firstpointis9999 ) {
            if(fNevTACSlp[cccc][np]<Nev_TStep){
              ++ncetoofew;
              fCellSusp[cccc] = fBits->SetBit(fCellSusp[cccc],np);
            }
          }
        }
        //Fill hist fNBadStepNev
        fHists.fNBadStepNev->Fill(ncetoofew);

        // If bit 1 in the cell word is not set i.e. we did see some events:
        if(!(fBits->TestBit(fCellStatus[cccc],1))){
          fCellStatus[cccc0] = fBits->ClearBit(fCellStatus[cccc0],14);
          if(ncetoofew>Nev_maxN){
            // Too many points had too few events
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],11);
            fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],23);
          }

          // 5. Test cubic fit: Npoints, parameters in range and diff from prev
          int nused = ECAProc::TSlopeNGood(cccc);
          // Fill hist fNUsed
          fHists.fNUsed->Fill(nused);
          if(nused<cubic_minNfit){
            // Not enough points used for the fit
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],2);
            fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],15);
          }

          if(nused>0){   // Do tests on fit params
            for(int nfp=0; nfp<3; ++nfp){
              runningsum[0][nfp] += fCellTACFits[cccc][nfp+1]; // only test linear+
              nrunsum[0][nfp] += 1;
              Crrunningsum[icrate][0][nfp] += fCellTACFits[cccc][nfp+1]; // only test linear+
              Crnrunsum[icrate][0][nfp] += 1;
              if((fCellTACFits[cccc][nfp+1]<cubic_min[nfp])
                 || (fCellTACFits[cccc][nfp+1]>cubic_max[nfp])){
                int bitid = 8 - nfp*2;
                fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],bitid);
                bitid = 19 - nfp;
                fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],bitid);
                fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],20);
              }
            }
            if(fECAdiff==1){  // Can we do difference tests?
              // Only do these tests if the previous run doesn't have bits 1/26 set
              bool testprev1 = fBits->TestBit(fPrevStatus[cccc],1);
              bool testprev2 = fBits->TestBit(fPrevStatus[cccc],26);
              if(testprev1 || testprev2) dodiff = false;
            }
            if(dodiff){  // Do difference tests
              for(int nfp=0; nfp<3; ++nfp){
                runningsum[1][nfp] += (fCellTACFits[cccc][nfp+1] -
                                       fPrevTACFits[cccc][nfp+1]);
                runningsum[2][nfp] += fabs(fCellTACFits[cccc][nfp+1] -
                                       fPrevTACFits[cccc][nfp+1]);
                nrunsum[1][nfp] += 1;
                Crrunningsum[icrate][1][nfp] += (fCellTACFits[cccc][nfp+1] -
                                       fPrevTACFits[cccc][nfp+1]);
                Crrunningsum[icrate][2][nfp] += fabs(fCellTACFits[cccc][nfp+1] -
                                       fPrevTACFits[cccc][nfp+1]);
                Crnrunsum[icrate][1][nfp] += 1;
                Crnrunsum[icrate][2][nfp] += 1;
                if(fabs(fCellTACFits[cccc][nfp+1] - fPrevTACFits[cccc][nfp+1])
                    > cubic_diff[nfp]){
                  int bitid = 9 - nfp*2;
                  fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],bitid);
                  bitid = 19 - nfp;
                  fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],bitid);
                  fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],21);
                }
                // Fill hists fCubFitPar[nfp][1]
                fHists.fCubFitPar[nfp][1]->Fill(fCellTACFits[cccc][nfp+1] - fPrevTACFits[cccc][nfp+1]);
              }
              fHists.fCubIntDiff->Fill(fCellTACFits[cccc][0] - fPrevTACFits[cccc][0]);
            }
          } //! nused>0

          // 6. Test widths of each point
          int ncebadwidth = 0;

          for(int np=0; np<fNSlopePoints; ++np){
          //again only check widths / suspicious criteria if
          // we're not in the post-reset region
          //ideally / eventually read in GT valid here instead
            if( np < first9999index && !firstpointis9999 ) {
            // Fill hist fWidth
              fHists.fWidth->Fill(fCellTACwidth[cccc][np]);
              fHists.fWidthCrates[icrate]->Fill(fCellTACwidth[cccc][np]);
              fHists.fWidthPoints[np]->Fill(fCellTACwidth[cccc][np]);
              if(fCellTACwidth[cccc][np]>width_TStep){
                ++ncebadwidth;
                fCellSusp[cccc] = fBits->SetBit(fCellSusp[cccc],np);
               }
            }

            if( np < second9999index && firstpointis9999 ) {
            // Fill hist fWidth
              fHists.fWidth->Fill(fCellTACwidth[cccc][np]);
              fHists.fWidthCrates[icrate]->Fill(fCellTACwidth[cccc][np]);
              fHists.fWidthPoints[np]->Fill(fCellTACwidth[cccc][np]);
              if(fCellTACwidth[cccc][np]>width_TStep){
                ++ncebadwidth;
                fCellSusp[cccc] = fBits->SetBit(fCellSusp[cccc],np);
              }
            }
          }
          // Fill hist fNBadCeWidth
          fHists.fNBadCeWidth->Fill(ncebadwidth);

          if(ncebadwidth>width_maxN){
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],10);
            fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],22);
          }

          // 7. Frati TAC jump test??? XXXXXXXXXXXXXXXXXXXX
          /*frati_flag=.false.
          if(Frati_ChiSq.ge.rcons(ldceca+start_ldceca+15))
     +                frati_flag=.true.
          if(Frati_ChiSq.lt.rcons(ldceca+start_ldceca+15)) then
            call hfill(2602,Frati_ChiSq,0.,1.)
            n_residual(ice)=0
            do K=1,NS
              if((qtgen(K).gt.59.).and.(qtgen(K).lt.350.).and.
     +                         (med_qt(K,ice,ich,ica,4).gt.0)) then
                do K2=2,4
                  fit_coeffs(K2)=Frati_C_Fit(K2)
                enddo
                fit_coeffs(1)=Frati_C_Fit(1)
     +                          -1.*med_qt(K,ice,ich,ica,4)
                t_residual(K,ice)=qtgen(K)
                call poly_solve(fit_coeffs,3,100,0.001,
     +                          t_residual(K,ice),err_dummy,iretcode)
                if(iretcode.eq.1) then
                  t_residual(K,ice)=-9999.
                else
                  t_residual(K,ice)=t_residual(K,ice)-qtgen(K)
                endif
                if(abs(t_residual(K,ice)).ge.
     +             rcons(ldceca+start_ldceca+16))
     +                               n_residual(ice)=n_residual(ice)+1
                call hfill(2601,t_residual(K,ice),0.,1.)
              endif
            enddo
            call hfill(2603,n_residual(ice)*1.,0.,1.)
            if(n_residual(ice).ge.icons(ldceca+start_ldceca+17))
     +                                             frati_flag=.true.
          endif
          if(frati_flag) then
            maskT(ice)=ibset(maskT(ice),12)
            maskT(1)=ibset(maskT(1),start_ch_mask+14)
          endif*/
          // 8. Test # suspicious points (susp = bad width or too few ev)
          int ncesusp = 0;
          for(int np=0; np<fNSlopePoints; ++np){
            if(fBits->TestBit(fCellSusp[cccc],np))++ncesusp;
          }

          // Fill hist fNSuspStep
          fHists.fNSuspStep->Fill(ncesusp);

          //Fill debug hist: first reset time
          //indeally would just read in GT valid from ORCA, but this gives the general idea
          if(fSuperDebug == 1){
            if(firstpointis9999) fResetTimes->Fill(fEVTACX[second9999index]);
            else fResetTimes->Fill(fEVTACX[first9999index]);
          }

          if(ncesusp>susp_maxN){
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],3);
            fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],16);
          }

          // 9. Test if any railed points
          int ncerailed = 0;
          for(int np=0; np<fNSlopePoints; ++np){
            if(fCellTACPoints[cccc][np] >= 4095){
              ++ncerailed;
            }
          }
          if(ncerailed >= NRailedPoints){
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],12); //bit 12: cell has tslope = railed
            fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],27); //bit 27: chan has >=1 railed tslope
          }

          //10. test if the first point is bad (i.e got no events)
          if(firstpointis9999)
            {
              fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],28); //bit 28: cell has bad first point on tslope
              fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],29);//bit 29: chan has >=1 cell w/ bad first point tslope
            }

          //11. test if there is any behavior in post-reset points that is not consisent with noise
          // criteria 1: if there are many hits, set the flag

          fHists.fPResetNPoints->Fill(numpostresethits);

          if(numpostresethits > 2)
            {
              fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],30); //bit 30: cell has susp post-reset
              fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],31); //bit 31: chan has >=1 cell w/ susp post-reset
            }
          //criteria 2: if the hits are higher TAC value than expected, set the flag

          //Do not use first bad point if any
          int used9999index = -1.;
          if(!firstpointis9999)
            used9999index = first9999index;
          else
            used9999index = second9999index;

          for(int np=used9999index +1; np<fNSlopePoints; ++np){
            float thispoint = fCellTACPoints[cccc][np];
            if( thispoint != -9999.0) {
              float timepastreset = fEVTACX[np] - fEVTACX[used9999index];
              // find the surrounding times of timepastreset from earlier in the TSlope
              // e.g if thispoint occured 60ns after GTValid, find the index/indicies that correspond
              // to 60ns from the beginning of the TSlope
              int ilow = -1;
              int ihigh = fNSlopePoints-1;
              for(int i=0;i<fNSlopePoints;++i){
                if((timepastreset >= fEVTACX[i]) && (i>ilow)){
                  ilow = i;
                }
                if((timepastreset <= fEVTACX[i]) && (i<ihigh)){
                  ihigh = i;
                }
              }
              //what ADC value it should have (i.e. y-interpolation for x-value = timepastreset
              float comparisonTACValue = -1.;
              if(ilow == ihigh){ //TAC has an ADC value at this time, no need to interpolate
                comparisonTACValue = fCellTACPoints[cccc][ilow];
              }
              else if(ilow < ihigh){
                comparisonTACValue = fCellTACPoints[cccc][ilow] + (fCellTACPoints[cccc][ihigh] - fCellTACPoints[cccc][ilow]) * ( (timepastreset - fEVTACX[ilow]) / (fEVTACX[ihigh] - fEVTACX[ilow]) );
              }
              else continue; //something has gone wrong skip this point
              float upperlim = comparisonTACValue + width_TStep*3;

              //Javi
              if(numpostresethits < 2) {
                fHists.fPResetDiffComp->Fill(thispoint - comparisonTACValue);
                // info <<" set nevents post reset"<< newline;
                fHists.fPResetNev->Fill(fNevTACSlp[cccc][np]);
              }

              if(thispoint > upperlim){
                fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],30); //bit 30: cell has susp post-reset
                fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],31); //bit 31: chan has >=1 cell w/ susp post-reset

                //Javi
                if(numpostresethits < 2) {
                  fHists.fPResetDiffCompFailed->Fill(thispoint - comparisonTACValue);
                  // info<<" set nevents post reset failed"<< newline;
                  fHists.fPResetNevFailed->Fill(fNevTACSlp[cccc][np]);
                  // for(int ievt=0;ievt<fNMaxEV;ievt++)
                  // fHists.FPResetEVTAC[cccc]->Fill(fEVTACPoints[nn][np][ievt]);
                }
              }
            }
          }
        }  // ! Events observed on this cell
      } // ! seq enabled
      else{ // not getting data
        fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],26);
      } // ! test of ECAChannel active
    } // !cells for this channel

    // Check first cell of this channel to determine `channel status' bits
    int channel0 = fBits->GetBits(fCellStatus[cccc0], 14, 14);
    // If any were bad, set the overall channel status bit for cell0
    if( (channel0 != 0) || (fBits->TestBit(fCellStatus[cccc0],29)) || (fBits->TestBit(fCellStatus[cccc0],31)) ) fCellStatus[cccc0] = fBits->SetBit(fCellStatus[cccc0],13);

    // Increment some per-cell counters, and copy cell0 channel status to the others
    for(int nce = 0; nce<fNcell;++nce){
      int cccc = fBits->GetCCCC(nch,nce);
      int cellstat = fBits->GetBits(fCellStatus[cccc], 1, 12);
      // If any tests were bad, set the overall cell status bit
      // don't forget the new cell-level flags at bit 28 and 30!
      if( (cellstat != 0) || (fBits->TestBit(fCellStatus[cccc],28)) || (fBits->TestBit(fCellStatus[cccc],30)) ){
        fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],0);
        ++NBadCe;

        if(fBits->TestBit(fCellStatus[cccc],1))++NEmptyCe;
        if(fBits->TestBit(fCellStatus[cccc],2))++NBadUsed;
        if(fBits->TestBit(fCellStatus[cccc],3))++NSusp;
        if(fBits->TestBit(fCellStatus[cccc],10))++NBadWidth;
        if(fBits->TestBit(fCellStatus[cccc],11))++NPartCe;
        for(int nfp=0; nfp<3; ++nfp){  // 3 cubic fit pars
          for(int nt=0;nt<2;++nt){     // 2 tests (range & diff)
            int bitid = 8 - 2*nfp + nt;
            if(fBits->TestBit(fCellStatus[cccc],bitid))++NBadFit[nt][nfp];
          }
          int bitid = 8 - 2*nfp;
          int testchunk = fBits->GetBits(fCellStatus[cccc], bitid, 2);
          if(testchunk!=0)++NBadFit[2][nfp];
        }
      }
      // Copy channel status from cell0
      if(nce>0){
        for(int nbit=13;nbit<28;++nbit){
          if(fBits->TestBit(fCellStatus[cccc0],nbit)){
            fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],nbit);
          }
        }
        //don't forget new flags
        if(fBits->TestBit(fCellStatus[cccc0],29)){
          fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],29);
        }
        if(fBits->TestBit(fCellStatus[cccc0],31)){
          fCellStatus[cccc] = fBits->SetBit(fCellStatus[cccc],31);
        }
      }
    }

    // Increment some channel counters
    if(fBits->TestBit(fCellStatus[cccc0],13)){
      ++NBadCh;
      ++NCrBadCh[icrate];
      if(fBits->TestBit(fCellStatus[cccc0],14)){
        ++NEmptyCh;
        ++NCrEmptyCh[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],15)){
        ++NBadUsedCh;
        ++NCrBadUsedCh[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],16)){
        ++NSuspCh;
        ++NCrSuspCh[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],22)){
        ++NBadWidthCh;
        ++NCrBadWidthCh[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],23)){
        ++NPartCh;
        ++NCrPartCh[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],24)){
        ++NBadMB;
        ++NCrBadMB[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],25)){
        ++NBadDB;
        ++NCrBadDB[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],26)){
        ++NBadSeq;
        ++NCrBadSeq[icrate];
      }
      if(fBits->TestBit(fCellStatus[cccc0],27)){ //kate
        ++NRailedCh;
      }
      if(fBits->TestBit(fCellStatus[cccc0],29)){ //kate
        ++NFirstPointBadCh;
      }
      if(fBits->TestBit(fCellStatus[cccc0],31)){ //kate
        ++NSuspPostResetHitsCh;
      }
      for(int nfp=0;nfp<3;++nfp){
        if(fBits->TestBit(fCellStatus[cccc0],19-nfp)){
          ++NBadFitCh[3][nfp];
          bool badfit = false;  // ?
          for(int nt=0; nt<2; ++nt){
            bool badtest = false;  // chan has not yet failed this test for this par
            for(int nce=0; nce<fNcell;++nce){
              int cccc = fBits->GetCCCC(nch,nce);
              int bitid = 8 - 2*nfp + nt;
              bool testbit = fBits->TestBit(fCellStatus[cccc],bitid);
              if((!badtest) && testbit){
                if(!badfit){
                  badfit = true;
                  ++NBadFitCh[2][nfp];
                }
                ++NBadFitCh[nt][nfp];
                ++NCrBadFitCh[icrate][nt][nfp];
                badtest = true;
              }
            }
          }
        }
        if(NBadFitCh[3][nfp]!=NBadFitCh[2][nfp]){
          ECAProc::Die("",6,NBadFitCh[3][nfp],NBadFitCh[2][nfp]);
        }
      }

      /*
        //Frati test not currently used, but keep incase we want to implement it
        if(fBits->TestBit(fCellStatus[cccc0],27)){
        ++NFratiCh;
      }
      */
    }

  }   // !channels

  // Loop over crates and fill hists
  for(int ic=0; ic<fNcrate; ++ic){
    fHists.fNBadChNfit->Fill(NCrBadUsedCh[ic]);
    fHists.fNBadChWidth->Fill(NCrBadWidthCh[ic]);
    fHists.fNBadChNev->Fill(NCrPartCh[ic]);
    fHists.fNBadCh0ev->Fill(NCrEmptyCh[ic]);
    fHists.fNBadChSeq->Fill(NCrBadSeq[ic]);
    fHists.fNBadChMB->Fill(NCrBadMB[ic]);
    fHists.fNBadChDB->Fill(NCrBadDB[ic]);
    fHists.fNBadChSusp->Fill(NCrSuspCh[ic]);
    fHists.fNBadCh->Fill(NCrBadCh[ic]);
    for(int ip=0; ip<3; ++ip){  // swap order: runningsum[test][fit par]
      for(int it=0; it<3; ++it){          //but  hist[par][test]
        float avg = 0.;
        if((Crnrunsum[ic][it][ip])>0.){avg = Crrunningsum[ic][it][ip] / Crnrunsum[ic][it][ip];}
        fHists.fCubFitParAve[ip][it]->Fill(avg);
        if(it<2)fHists.fNBadChCub[ip][it]->Fill(NCrBadFitCh[ic][it][ip]);
      }
    }
  }

  // NOW WE DO THE FULL RUN STATUS CHECKS

  // Calculate run averages
  for(int nfp=0; nfp<3; ++nfp){
    for(int nt=0; nt<2; ++nt){  // actual / diff from prev / abs diff from prev
      if(nrunsum[nt][nfp]>0){   //
        runningsum[nt][nfp] = runningsum[nt][nfp]/nrunsum[nt][nfp];
        if(nt==1)runningsum[2][nfp] = runningsum[2][nfp]/nrunsum[nt][nfp];
      }
      else{
        runningsum[nt][nfp] = 0.0;
        if(nt==1)runningsum[2][nfp] = 0.0;
      }
    }

    // Test absolute values of fit pars
    if((runningsum[0][nfp]<TslpMinAveCub[nfp]) ||
      (runningsum[0][nfp]>TslpMaxAveCub[nfp])){
      fOverallMask = fBits->SetBit(fOverallMask, 12-nfp);
    }

    // Test fit parameter abs diffs
    if(runningsum[2][nfp]>TslpMaxAveCubDiff[nfp]){
      fOverallMask = fBits->SetBit(fOverallMask, 15-nfp);
    }

  }

  // Tests on overall run stats:

  // 0. Test UC bit
  if(fBits->TestBit(fRunMask, 21)){
    fOverallMask = fBits->SetBit(fOverallMask, 26);
  }

  // 1. Number of board ID failures (DB+MB)
  if((NBadDB+NBadMB) > TslpMaxNboardID){
    fOverallMask = fBits->SetBit(fOverallMask, 0);
  }

  // 2. Number of empty channels
  if(NEmptyCh > TslpMaxN0ev){
    fOverallMask = fBits->SetBit(fOverallMask, 1);
  }

  // 3. Number of partial channels (i.e. too few events)
  if(NPartCh > TslpMaxN0evTStep){
    fOverallMask = fBits->SetBit(fOverallMask, 2);
  }

  // 4. Number of bad width channels
  if(NBadWidthCh > TslpMaxNWidth){
    fOverallMask = fBits->SetBit(fOverallMask, 3);
  }

  // 5. Number of channels failing Nev in fit
  if(NBadUsedCh > TslpMaxNFit){
    fOverallMask = fBits->SetBit(fOverallMask, 4);
  }

  // 6. Number of suspicious channels
  if(NSuspCh > TslpMaxNSusp){
    fOverallMask = fBits->SetBit(fOverallMask, 5);
  }

  // 7. Number of channels failing Frati TAC jump test
  if(NFratiCh > TslpMaxNTAC){
    fOverallMask = fBits->SetBit(fOverallMask, 6);
  }

  // 8. Total flagged channels
  if(NBadCh > TslpMaxNChan){
    fOverallMask = fBits->SetBit(fOverallMask, 16);
  }

  // 9. Number of channels failing cubic fit par tests
  for(int nfp=0;nfp<3;++nfp){
    if(NBadFitCh[3][nfp] > TslpMaxNCubic[nfp]){
      fOverallMask = fBits->SetBit(fOverallMask, 9-nfp);
    }
  }

  // 10. Number of railed channels //kate
  if(NRailedCh > TslpMaxNRailedChan){
    fOverallMask = fBits->SetBit(fOverallMask, 17);
  }

  // 11. Number of railed channels //kate
  if(NFirstPointBadCh > TslpMaxNFirstPointBadChan){
    fOverallMask = fBits->SetBit(fOverallMask, 18);
  }

  // 12. Number of railed channels //kate
  if( NSuspPostResetHitsCh> TslpMaxNSuspPostResetHitsChan){
    fOverallMask = fBits->SetBit(fOverallMask, 19);
  }


  delete[] NCrBadUsedCh;
  delete[] NCrBadWidthCh;
  delete[] NCrPartCh;
  delete[] NCrEmptyCh;
  delete[] NCrBadSeq;
  delete[] NCrBadMB;
  delete[] NCrBadDB;
  delete[] NCrSuspCh;
  delete[] NCrBadCh;
}



std::vector<int> ECAProc::TSlopeFilter(int cccc)
{
  // Diagnostic pars
  int MinNoEvPP = fECAbank->GetI("cubic_minN_TStep");

  int ntot = fNSlopePoints;
  std::vector<int> bad = fCellBad[cccc];
  std::vector<int> curl = fCellCurl[cccc];

  for(int ni = 0; ni<ntot; ++ni){
    // If the median at this point is -9999, exclude from this check (no events)
    if(fCellTACPoints[cccc][ni] == -9999.0){
      bad = fBits->SetBit(bad, ni);
    }
    // Check whether the points are monotonically increasing
    //increase[ni] = true;
    // SNO+ method: keep all points that monot. inc. UNTIL one turns over
    else{
      if(ni>0){
        if((fCellTACPoints[cccc][ni]<fCellTACPoints[cccc][ni-1]) || (fBits->TestBit(curl,ni-1))){
          //increase[ni] = false;
          curl = fBits->SetBit(curl, ni);

          //SNO method: tag 4 points before the turn over as curl region as well
          if( (isSNOdata) && !fBits->TestBit(curl,ni-1) ){
            curl = fBits->SetBit(curl, ni-1);
            curl = fBits->SetBit(curl, ni-2);
            curl = fBits->SetBit(curl, ni-3);
            curl = fBits->SetBit(curl, ni-4);
          }
        }
      }
    }

    //SNO method: ALWAYS tag the 4 last points
    if(isSNOdata && ni==ntot-1){
      curl = fBits->SetBit(curl, ni);
      curl = fBits->SetBit(curl, ni-1);
      curl = fBits->SetBit(curl, ni-2);
      curl = fBits->SetBit(curl, ni-3);
    }

    // If we have <min # events for this point, flag as bad
    if(fNevTACSlp[cccc][ni] < MinNoEvPP){
      bad = fBits->SetBit(bad, ni);
    }
  }

  fCellBad[cccc] = bad;
  fCellCurl[cccc] = curl;

  return bad;

}

int ECAProc::TSlopeNGood(int cccc)
{
  // THIS ROUTINE IS ONLY USEFUL IF TSLOPEFILTER HAS BEEN CALLED FIRST
  // To set flags for the bad points in fCellBad[cccc]
  int ntot = fNSlopePoints;
  int nbad = 0;
  bool badbit = true;
  // Count up bad points!!
  for(int nb=0; nb<ntot;++nb){
    badbit = (fBits->TestBit(fCellBad[cccc],nb) || fBits->TestBit(fCellCurl[cccc],nb));
    if(badbit) ++nbad;
  }
  int ngood = ntot - nbad;

  return ngood;
}

void ECAProc::TSlopeFit(int cccc)
{
  int ntot = fNSlopePoints;
  int ngood = ECAProc::TSlopeNGood(cccc);
  float *yval = new float[ngood];
  float *xval = new float[ngood];
  float *yerr = new float[ngood];
  int bit = 0;

  //use these to set range on fit plot
  float firstgood = 0;
  float lastgood = 0;

  // Only fill arrays with the good points
  for(int ia=0; ia<ntot;++ia){
    bool badbit = fBits->TestBit(fCellBad[cccc],ia) || fBits->TestBit(fCellCurl[cccc],ia);
    if(!badbit){
      yval[bit] = fCellTACPoints[cccc][ia];
      yerr[bit] = fCellTACrms[cccc][ia];
      xval[bit] = fEVTACX[ia];
      ++bit;
    }
  }

  firstgood = xval[0];
  lastgood = xval[ngood-1];

  // Linear fit
  float *param = new float[4];
  float *parerr = new float[4];
  for(int np = 0; np<4; ++np){
    param[np] = 0.0;
    parerr[np] = 0.0;
  }
  float valid = 0.0;
  ECAProc::LinearFit(ngood, xval, yval, param, valid);
  float chisq = 0.0;
  //save these parameters for the tslope plots
  float lin0 = param[0];
  float lin1 = param[1];


  // Cubic fit, iff we have >3 good points on the slope
  if(ngood > 3){
    ECAProc::CubicFit(ngood, xval, yval, yerr, param, parerr, chisq);
  }
  else{
    for(int i=0;i<ngood;++i){
      chisq += pow(((param[1]*xval[i] + param[0] - yval[i])/yerr[i]),2);
    }
  }

  if(fSuperDebug==1){
    std::cout << "For cell" << cccc << " CubicFit Chi2 = " << chisq << std::endl;
  }

  // Copy fit pars into final array and fill hists fCubFitPar[np][0], fCubInt
  for(int np = 0; np<4; ++np){
    fCellTACFits[cccc][np] = param[np];
    fCellTACFitErrs[cccc][np] = parerr[np]; //kate
    if(np<3)fHists.fCubFitPar[np][0]->Fill(param[np+1]);
    if(np==0)fHists.fCubInt->Fill(param[np]);
  }

  //kate
  if(fSuperDebug == 2)
  {
    //Get the TSlopes
    float *yvald = new float[ntot];
    float *xvald = new float[ntot];
    float *yerrd = new float[ntot];
    float *xerrd = new float[ntot];
    int pt = 0;

    for(int ia=0; ia<ntot;++ia){ //put all points on tlsope plot in the output
      yvald[pt] = fCellTACPoints[cccc][ia];
      yerrd[pt] = fCellTACrms[cccc][ia];
      xvald[pt] = fEVTACX[ia];
      xerrd[pt] = 0;
      ++pt;
    }
    TGraphErrors* tslope_cccc = new TGraphErrors(fNSlopePoints, xvald, yvald, xerrd, yerrd);
    std::string info = dformat("cell%i", cccc);
    tslope_cccc->GetYaxis()->SetTitle("ADC Counts");
    std::string name = info + std::string(" time (ns)");
    tslope_cccc->GetXaxis()->SetTitle(name.c_str());
    tslope_cccc->SetName(info.c_str());

    //LINEAR FIT
    TF1* linearfit_cccc = new TF1(info.c_str(), "[0] + [1]*x", firstgood, lastgood);
    linearfit_cccc->SetParameter(0, lin0);
    linearfit_cccc->SetParameter(1, lin1);

    ///CUBIC FIT
    TF1 *cubicfit_cccc = new TF1(info.c_str(),"[0] + [1]*x + [2]*x*x + [3]*x*x*x",firstgood,lastgood);
    if(ngood >3){
      cubicfit_cccc->SetParameter(0, param[0]);
      cubicfit_cccc->SetParameter(1, param[1]);
      cubicfit_cccc->SetParameter(2, param[2]);
      cubicfit_cccc->SetParameter(3, param[3]);
      cubicfit_cccc->SetParError(0, parerr[0]);
      cubicfit_cccc->SetParError(1, parerr[1]);
      cubicfit_cccc->SetParError(2, parerr[2]);
      cubicfit_cccc->SetParError(3, parerr[3]);
    }
    else{
      cubicfit_cccc->SetParameter(0, 0.0);
      cubicfit_cccc->SetParameter(1, 0.0);
      cubicfit_cccc->SetParameter(2, 0.0);
      cubicfit_cccc->SetParameter(3, 0.0);
      cubicfit_cccc->SetParError(0, 0.0);
      cubicfit_cccc->SetParError(1, 0.0);
      cubicfit_cccc->SetParError(2, 0.0);
      cubicfit_cccc->SetParError(3, 0.0);
    }
    //can add option here to save only tslopes that have out of range parameters
    fTSlopes.push_back(tslope_cccc);
    fLinearFits.push_back(linearfit_cccc);
    fCubicFits.push_back(cubicfit_cccc);
  }//end if debug->save cubic fits

  delete[] param;
  delete[] parerr;
}
void ECAProc::WriteECADebugHistosRoot()  // added by kate
{

  // Write graphs and fits to file
  fdebughistosROOTfile->cd();

  //write trigger histograms
  fTrigWordHistSkippedEV->Write();
  fTrigWordHist->Write();

  //write triggers vs time (tslopes only)
  if(fECAflag == 2)
  {
    fResetTimes->Write();
    fb0_NHit100Lo->Write();
    fb1_NHit100Med->Write();
    fb2_NHit100Hi->Write();
    fb3_NHit20->Write();
    fb4_NHit20LB->Write();
    fb5_ESumLo->Write();
    fb6_ESumHi->Write();
    fb7_OWLN->Write();
    fb8_OWLELo->Write();
    fb9_OWLEHi->Write();
    fb10_PulseGT->Write();
    fb11_Prescale->Write();
    fb12_Pedestal->Write();
    fb13_Pong->Write();
    fb14_Sync->Write();
    fb15_EXTASY->Write();
    fb16_Ext2->Write();
    fb17_Ext3->Write();
    fb18_Ext4->Write();
    fb19_Ext5->Write();
    fb20_Ext6->Write();
    fb21_Ext7->Write();
    fb22_Ext8PulseAsync->Write();
    fb23_SpecialRaw->Write();
    fb24_NCD->Write();
    fb25_SoftGT->Write();
    fb26_MissedTrig->Write();
  }
  //close file
  fdebughistosROOTfile->Close();
}

void ECAProc::WriteTSlopesRoot()  // added by kate
{
  // Write tslope plotss and fits to file
  fTSlopesRootFile->cd();

  if(fECAflag==2){
    int size = fTSlopes.size();
    for(int i=0; i < size; i++)
    {
      fTSlopes[i]->SetMarkerSize(2);
      fTSlopes[i]->SetLineWidth(0);
      fTSlopes[i]->Write();
      fLinearFits[i]->Write("",TObject::kSingleKey);
      fCubicFits[i]->Write("",TObject::kSingleKey);
    }

  } //tslope plots extra diagnostics/debug ooutput

  // Close file
  fTSlopesRootFile->Close();
}



void ECAProc::LinearFit(int n, float* &x, float* &y, float* &p, float valid)
{
  float xsum = 0.0;
  float x2sum = 0.0;
  float ysum = 0.0;
  float xysum = 0.0;
  for(int i = 0; i < n; ++i){
    xsum += x[i];
    ysum += y[i];
    x2sum += pow(x[i],2);
    xysum += (x[i]*y[i]);
  }
  valid = n*x2sum - pow(xsum,2);
  if(valid==0.0){
    p[0] = 0.0;
    p[1] = 0.0;
    return;
  }
  p[0] = (x2sum*ysum - xsum*xysum)/valid;
  p[1] = (n*xysum - xsum*ysum) / valid;
  return;
}

void ECAProc::CubicFit(int n, float* &x, float* &y, float* &ye, float* &p, float* &pe, float& chisq)
{
  float *xe = new float[n];
  for(int i=0;i<n;++i){
    xe[i] = 0.0;
  }
  TGraphErrors *points = new TGraphErrors(n,x,y,xe,ye);
  TF1 *fit = new TF1("fit","[0] + [1]*x + [2]*x*x + [3]*x*x*x",0.0,410.0);
  fit->SetParameters(p[0],p[1],p[2],p[3]);
  points->Fit(fit, "MERQ");
  chisq = fit->GetChisquare();
  for(int i=0;i<4;++i){
    p[i] = fit->GetParameter(i);
    pe[i] = fit->GetParError(i);
  }
  delete points;
}

Processor::Result ECAProc::DSEvent(DS::Run& run, DS::Entry& ds)
{
  if(fECAflag == 0){
    warn<<"ECAProc::DSEvent: No ECA process selected, exiting now\n";
    return OK;
  }
  //  info<<"ECAProc::DSEvent: " <<ds.GetEVCount()<<" events\n";

  // Check whether this event has an EPED branch and, if so, get info
  if(ds.HeaderInfoExists()){
    const DS::HeaderInfo& header = ds.GetHeaderInfo();
    //Retrieve the pulser rate of the first trigger headers only if
    //ECA header exists
    if( header.TrigHeaderExists() && header.ECAHeaderExists() && fECArate == 0){
      const DS::TrigHeader& trigheaders = header.GetTrigrHeader();
      //Make up trigger rate for SNO data
      if(isSNOdata)
        fECArate = 10;
      else
        fECArate = trigheaders.GetPulserRate();
    }

    if(header.ECAHeaderExists()){
      const DS::ECAHeader& eped = header.GetECAHeader();
      if( eped.GetFlag()<3 ){ //values with flag>2 are set incorrectly so keep with the same headers than before
        fGTDelC = eped.GetGTDelayCoarse();
        fGTDelF = eped.GetGTDelayFine();
        fPedWidth = eped.GetChargePedestalWidth();
        fPattern = eped.GetPatternID();
        fCalType = eped.GetCalType();
      }
      // info << "Pattern " <<fLastPatternID<<" "<<fPattern<<newline;

      //Javi: add SNO compatibility
      if((!isSNOdata && (unsigned long)(fLastPatternID+1) == fPattern)
         || (isSNOdata && (unsigned long)fLastPatternID != fPattern)){
        if(fECAlogverb>0) info << "New step: " << fPattern << newline;
        if(fECAlogverb>0) info << "fArrayID = " << fArrayID << newline;
        fTinj = 0.;
        fArrayID = -1;
      }
      fLastPatternID = fPattern;

      // CalType is in the format XY where X = ped/Tslope and Y = pattern set identification

      //kate - skip events w/ outof range ECAType
      int ECATypeFromfCalType = (fCalType / 10) % 10;
      if((static_cast<int>(ECATypeFromfCalType) < 1) | (static_cast<int>(ECATypeFromfCalType) > 3)){
        ECAProc::DontDie("",8,ECATypeFromfCalType,fECAflag);
        return OK; //skip to next event
      }
      else if(static_cast<int>((fCalType/10))!=fECAflag) ECAProc::Die("",8,fCalType,fECAflag);

      fPattSet = fCalType % 10;

      /* Calculate injected T for this point */
      // SNO+: Internal delay = 10ns, Coarse delay -> ns, fine delay -> ps
      float time = 10. + static_cast<float>(fGTDelC) + 1e-3*static_cast<float>(fGTDelF);

      // SNO: Internal delay = 18.35ns (wrong and corrected in SNO+), Coarse delay -> ns, fine delay -> ps
      if(isSNOdata) time = 18.35 + 10.*(255. - static_cast<float>(fGTDelC)) + (25./256.0)*static_cast<float>(fGTDelF);

      // Only need to test time if doing a TSlope run
      if(fECAflag==2){
        if(time>fTinj){ // next (or first) Tinj point
          ++fArrayID;
          fTinj = time;
          fEVTACX[fArrayID] = fTinj;
        }
        if(time<fTinj)ECAProc::Die("",9,time,fTinj);
      } // tslope run
    } //exist eped info
    else{
      ECAProc::DontDie("",9,0,1);
      return OK;
    }
  } //exist header info
  else
  {
    ECAProc::DontDie("",9,1,0);
    return OK;
  }


  // GET channel pattern information
  // kate edit to catch segfaults (GDOG: use ratdb values instead of hard coded)
  bool fail_set = 0;
  bool fail_patt_num = 0;
  if(fPattSet < 1 || (int)fPattSet > fNSet+1)
  {
    fail_set = true;
  }
  else if(fPattSet-1 >= fNPatt.size() || (int)fPattern > fNPatt[fPattSet-1])
  {
    fail_patt_num = true;
  }

  if(fail_set || fail_patt_num)
  {
    ECAProc::Die("", 11, fPattSet, fPattern, 0);
  }
  else fCurrentPattern = fPatterns[fPattSet-1][fPattern];


  // Loop over events (presumably there should only ever be 1 in here?)
  for (size_t iev=0; iev<ds.GetEVCount(); ++iev){
    DS::EV& ev = ds.GetEV(iev);
    size_t nchan = ev.GetUncalPMTs().GetAllCount();

    // Get Pedestal bit and check this is a true pedestal
    int trigType = ev.GetTrigType();

    //plot triggers as function of time
    if(fSuperDebug == 1)
    {
      if(fECAflag == 2)
      {
        if(fBits->TestBit(trigType,0)) fb0_NHit100Lo->Fill(fTinj);
        if(fBits->TestBit(trigType,1)) fb1_NHit100Med->Fill(fTinj);
        if(fBits->TestBit(trigType,2)) fb2_NHit100Hi->Fill(fTinj);
        if(fBits->TestBit(trigType,3)) fb3_NHit20->Fill(fTinj);
        if(fBits->TestBit(trigType,4)) fb4_NHit20LB->Fill(fTinj);
        if(fBits->TestBit(trigType,5)) fb5_ESumLo->Fill(fTinj);
        if(fBits->TestBit(trigType,6)) fb6_ESumHi->Fill(fTinj);
        if(fBits->TestBit(trigType,7)) fb7_OWLN->Fill(fTinj);
        if(fBits->TestBit(trigType,8)) fb8_OWLELo->Fill(fTinj);
        if(fBits->TestBit(trigType,9)) fb9_OWLEHi->Fill(fTinj);
        if(fBits->TestBit(trigType,10)) fb10_PulseGT->Fill(fTinj);
        if(fBits->TestBit(trigType,11)) fb11_Prescale->Fill(fTinj);
        if(fBits->TestBit(trigType,12)) fb12_Pedestal->Fill(fTinj);
        if(fBits->TestBit(trigType,13)) fb13_Pong->Fill(fTinj);
        if(fBits->TestBit(trigType,14)) fb14_Sync->Fill(fTinj);
        if(fBits->TestBit(trigType,15)) fb15_EXTASY->Fill(fTinj);
        if(fBits->TestBit(trigType,16)) fb16_Ext2->Fill(fTinj);
        if(fBits->TestBit(trigType,17)) fb17_Ext3->Fill(fTinj);
        if(fBits->TestBit(trigType,18)) fb18_Ext4->Fill(fTinj);
        if(fBits->TestBit(trigType,19)) fb19_Ext5->Fill(fTinj);
        if(fBits->TestBit(trigType,20)) fb20_Ext6->Fill(fTinj);
        if(fBits->TestBit(trigType,21)) fb21_Ext7->Fill(fTinj);
        if(fBits->TestBit(trigType,22)) fb22_Ext8PulseAsync->Fill(fTinj);
        if(fBits->TestBit(trigType,23)) fb23_SpecialRaw->Fill(fTinj);
        if(fBits->TestBit(trigType,24)) fb24_NCD->Fill(fTinj);
        if(fBits->TestBit(trigType,25)) fb25_SoftGT->Fill(fTinj);
        if(fBits->TestBit(trigType,26)) fb26_MissedTrig->Fill(fTinj);

      }//end if tslope run
    }//end fSuperDebug for histograming triggers at times

    // want both the pedestal and the PGT bit set
    if(!((fBits->TestBit(trigType,10))&&(fBits->TestBit(trigType,12)))){
      //info << "Event skipped" << newline;
      ++fNumSkip;
      // fill some trigger word histograms to check contamination
      if(fSuperDebug == 1){
        if(fBits->TestBit(trigType,0)) fTrigWordHistSkippedEV->Fill(0);
        if(fBits->TestBit(trigType,1)) fTrigWordHistSkippedEV->Fill(1);
        if(fBits->TestBit(trigType,2)) fTrigWordHistSkippedEV->Fill(2);
        if(fBits->TestBit(trigType,3)) fTrigWordHistSkippedEV->Fill(3);
        if(fBits->TestBit(trigType,4)) fTrigWordHistSkippedEV->Fill(4);
        if(fBits->TestBit(trigType,5)) fTrigWordHistSkippedEV->Fill(5);
        if(fBits->TestBit(trigType,6)) fTrigWordHistSkippedEV->Fill(6);
        if(fBits->TestBit(trigType,7)) fTrigWordHistSkippedEV->Fill(7);
        if(fBits->TestBit(trigType,8)) fTrigWordHistSkippedEV->Fill(8);
        if(fBits->TestBit(trigType,9)) fTrigWordHistSkippedEV->Fill(9);
        if(fBits->TestBit(trigType,10)) fTrigWordHistSkippedEV->Fill(10);
        if(fBits->TestBit(trigType,11)) fTrigWordHistSkippedEV->Fill(11);
        if(fBits->TestBit(trigType,12)) fTrigWordHistSkippedEV->Fill(12);
        if(fBits->TestBit(trigType,13)) fTrigWordHistSkippedEV->Fill(13);
        if(fBits->TestBit(trigType,14)) fTrigWordHistSkippedEV->Fill(14);
        if(fBits->TestBit(trigType,15)) fTrigWordHistSkippedEV->Fill(15);
        if(fBits->TestBit(trigType,16)) fTrigWordHistSkippedEV->Fill(16);
        if(fBits->TestBit(trigType,17)) fTrigWordHistSkippedEV->Fill(17);
        if(fBits->TestBit(trigType,18)) fTrigWordHistSkippedEV->Fill(18);
        if(fBits->TestBit(trigType,19)) fTrigWordHistSkippedEV->Fill(19);
        if(fBits->TestBit(trigType,20)) fTrigWordHistSkippedEV->Fill(20);
        if(fBits->TestBit(trigType,21)) fTrigWordHistSkippedEV->Fill(21);
        if(fBits->TestBit(trigType,22)) fTrigWordHistSkippedEV->Fill(22);
        if(fBits->TestBit(trigType,23)) fTrigWordHistSkippedEV->Fill(23);
        if(fBits->TestBit(trigType,24)) fTrigWordHistSkippedEV->Fill(24);
        if(fBits->TestBit(trigType,25)) fTrigWordHistSkippedEV->Fill(25);
        if(fBits->TestBit(trigType,26)) fTrigWordHistSkippedEV->Fill(26);
      }// if debug
      continue;
    }

    //fill hist for NOT skipped events
    if(fSuperDebug == 1){
      if(fBits->TestBit(trigType,0)) fTrigWordHist->Fill(0);
      if(fBits->TestBit(trigType,1)) fTrigWordHist->Fill(1);
      if(fBits->TestBit(trigType,2)) fTrigWordHist->Fill(2);
      if(fBits->TestBit(trigType,3)) fTrigWordHist->Fill(3);
      if(fBits->TestBit(trigType,4)) fTrigWordHist->Fill(4);
      if(fBits->TestBit(trigType,5)) fTrigWordHist->Fill(5);
      if(fBits->TestBit(trigType,6)) fTrigWordHist->Fill(6);
      if(fBits->TestBit(trigType,7)) fTrigWordHist->Fill(7);
      if(fBits->TestBit(trigType,8)) fTrigWordHist->Fill(8);
      if(fBits->TestBit(trigType,9)) fTrigWordHist->Fill(9);
      if(fBits->TestBit(trigType,10)) fTrigWordHist->Fill(10);
      if(fBits->TestBit(trigType,11)) fTrigWordHist->Fill(11);
      if(fBits->TestBit(trigType,12)) fTrigWordHist->Fill(12);
      if(fBits->TestBit(trigType,13)) fTrigWordHist->Fill(13);
      if(fBits->TestBit(trigType,14)) fTrigWordHist->Fill(14);
      if(fBits->TestBit(trigType,15)) fTrigWordHist->Fill(15);
      if(fBits->TestBit(trigType,16)) fTrigWordHist->Fill(16);
      if(fBits->TestBit(trigType,17)) fTrigWordHist->Fill(17);
      if(fBits->TestBit(trigType,18)) fTrigWordHist->Fill(18);
      if(fBits->TestBit(trigType,19)) fTrigWordHist->Fill(19);
      if(fBits->TestBit(trigType,20)) fTrigWordHist->Fill(20);
      if(fBits->TestBit(trigType,21)) fTrigWordHist->Fill(21);
      if(fBits->TestBit(trigType,22)) fTrigWordHist->Fill(22);
      if(fBits->TestBit(trigType,23)) fTrigWordHist->Fill(23);
      if(fBits->TestBit(trigType,24)) fTrigWordHist->Fill(24);
      if(fBits->TestBit(trigType,25)) fTrigWordHist->Fill(25);
      if(fBits->TestBit(trigType,26)) fTrigWordHist->Fill(26);
    }// if debug

    //info << "Event seen!" << newline;

    // Loop over each hit
    // On ALL channels, not just regular, inward-facing tubes
    for(size_t inh = 0; inh < nchan; ++inh){
      DS::PMTUncal& pmt = ev.GetUncalPMTs().GetAllPMT(inh);

      int lcn = pmt.GetID();
      int cell = pmt.GetCellID();
      int cccc = pmt.GetCrateCardChannelCell();

      // Check whether this channel is being ped-pinged (if not, it's coinc. noise)
      bool goodhit = false;
      for(UInt_t igh = 0; igh < fCurrentPattern.size(); ++igh){
        if(lcn == fCurrentPattern[igh]) goodhit = true;
      }
      if(!goodhit){
        if(fECAlogverb>0) ecalog << " Hit skipped in Cr " << pmt.GetCrate() << " Card " << pmt.GetCard() << " Channel " <<pmt.GetChannel() << newline;
        ++fHitSkip;
        continue;
      }

      // Mark this channel as having been tested
      fTestChannel[lcn] = 1;

      // Get uncalibrated charges and time
      unsigned short qhs = pmt.GetQHS();
      unsigned short qhl = pmt.GetQHL();
      unsigned short qlx = pmt.GetQLX();
      unsigned short tac = pmt.GetTime();

      if(fECAflag==1){
        // check that we don't already have >10 events in this cell
        // if so, skip it
        if(fNevCell[cccc]==fNMaxEV){continue;}

        // For first 16 events on this channel, store cell#
        int nchanev = 0;
        for(int ic=0; ic<fNcell; ++ic){
          int tempcccc = fBits->GetCCCC(lcn,ic);
          nchanev += fNevCell[tempcccc];
        }
        if(nchanev<16){
          fCellNo[lcn][nchanev] = cell;
        }

        fEVQHS[cccc][fNevCell[cccc]] = qhs;
        fEVQHL[cccc][fNevCell[cccc]] = qhl;
        fEVQLX[cccc][fNevCell[cccc]] = qlx;
        fEVTAC[cccc][fNevCell[cccc]] = tac;
        // Add 1 to the event count for this cell
        fNevCell[cccc] += 1;
      }
      if(fECAflag==2){
        if(fArrayID<0)ECAProc::Die("",10,0,0);
        // Check # events at this point, and fill arrays
        // First event at this point, fill arrays:
        if(fNevTACSlp[cccc][fArrayID]==fNMaxEV){continue;} // Already got 10 points
        fEVTACPoints[cccc][fArrayID][fNevTACSlp[cccc][fArrayID]] = tac;
        // Add 1 to the event count for this point of this cell
        fNevTACSlp[cccc][fArrayID] += 1;
      }
    } // ! loop over channels for this event
  } // ! event loop
  ++fNDSEvents;
  return OK;
}

void ECAProc::InitLogOut(std::string _logfilename)
{
  int count = ecalog.device_count();
  for (int i = 0; i < count; i++)ecalog.remove(0); // Always remove head device
  fECAlogname = _logfilename;
  fECAlogfile.open(fECAlogname);
  ecalog.add(fECAlogfile);
}

void ECAProc::InitECADebugOut(std::string _debugfilename)
{
  int count = ecadebugrecord.device_count();
  for (int i = 0; i < count; i++)ecadebugrecord.remove(0); // Always remove head device
  fECAdebugdbname = _debugfilename;
  fECAdebugdbfile.open(fECAdebugdbname);
  ecadebugrecord.add(fECAdebugdbfile);
}

void ECAProc::InitECARoot(std::string _rootfilename)
{
  fECArootname = _rootfilename;
  TString hold = fECArootname;
  // Open root file
  fECArootfile = new TFile(hold, "RECREATE");
  // Initialise required histograms
  fHists.Beautify();
  fHists.InitialiseHists(fECAflag, fDoDroop);
}

void ECAProc::WriteECADebugHeader()
{
  if(fECAflag==0){
    ecadebugrecord << "You don't have an ECA calibration type selected"<<newline;
  }
  if(fECAflag==1){
    ecadebugrecord<<"You're running pedestal, so there is no extra info to put in a ratdb file" <<newline;
  }
  if(fECAflag==2){
    ecadebugrecord<<"////////////////////////////////////////////////////////////////////" <<newline;
    ecadebugrecord<<"// Last svn revision: $Id:$" <<newline;
    ecadebugrecord<<"////////////////////////////////////////////////////////////////////" <<newline;
    ecadebugrecord<<"///" <<newline;
    ecadebugrecord<<"/// \\brief   TSLP Debug File output by the ECA for run "<< fRunNo <<newline;
    ecadebugrecord<<"///" <<newline;
    ecadebugrecord<<"///" <<newline;
    ecadebugrecord<<"////////////////////////////////////////////////////////////////////" <<newline;
    ecadebugrecord<<"{" <<newline;
    ecadebugrecord<<"type: \"ECA_DBUG\"," <<newline;
    ecadebugrecord<<"run_range : [ "<< fRunNo << ", 1000000]," <<newline;
    ecadebugrecord<<"" <<newline;
    ecadebugrecord<<"" <<newline;
    ecadebugrecord<<"" <<newline;
  }
}

void ECAProc::Summarise()
{
  // Init status words
  fECABits.InitialiseStatusWords(fECAflag);
  // Print out non-zero bits of overall status mask
  // Or every bit if user has requested extra verbosity
  if(fOverallMask!=0){
    ecalog<<"ECA run-status word != 0:"<<newline;
    for(int n=0;n<32;++n){
      int bit = fBits->GetBits(fOverallMask, n, 1);
      if(fECAlogverb>0||bit!=0){
        ecalog<<"Bit "<<n<<" = "<<bit<<" : "<<fECABits.fOverallMask[n]<<newline;
      }
    }
  }
  else{
    ecalog<<"ECA run-status word = 0 i.e. SUCCESS: a clean run"<<newline;
  }

  // Print out number of events skipped because not pedestal events
  ecalog<<fNumSkip<< " non-pedestal events skipped."<<newline;
  ecalog<<fHitSkip<< " non-pedestal hits skipped." << newline;
  //kate
  ecalog<<fTimeSkipped<< " events skipped because time out of order " <<newline;
  // Print out number of channels no tested
  ecalog<<fNChanNotTested<< " channels not tested! " <<newline;


}

void ECAProc::CreateECATable()
{

  if(fECAmerge==1){  // We want to merge the results
    ECAProc::MergeECAResults();
  }

  RAT::DBTable *ECA_table;
  if(fECAflag==0){
    ECA_table = new DBTable("ECA_INVALID");
  }
  if(fECAflag==1){
    ECA_table = new DBTable("ECA_PDST",fECAIndex);
    ECA_table->SetRunRange(fRunNo,fRunRangeMax);
    //pass = -2 means that this table will be pushed to RATDB in truncated mode
    ECA_table->SetPassNumber(-2);
    ECA_table->SetI("version",3);
    ECA_table->SetI("run_status",fOverallMask);
    // PDST: we want pedestals for 3Q+T, and status word, per CELL i.e. 16*channels
    ECA_table->SetIArray("pdst_status",fCellStatus);
    ECA_table->SetDArray("pdst_qhs",fCellQHS);
    ECA_table->SetDArray("pdst_qhl",fCellQHL);
    ECA_table->SetDArray("pdst_qlx",fCellQLX);
    ECA_table->SetDArray("pdst_tac",fCellTAC);
    ECA_table->SetDArray("pdst_qhswidth",fCellWidthQHS);
    ECA_table->SetDArray("pdst_qhlwidth",fCellWidthQHL);
    ECA_table->SetDArray("pdst_qlxwidth",fCellWidthQLX);
    ECA_table->SetDArray("pdst_tacwidth",fCellWidthTAC);
    ECA_table->SetIArray("pdst_dqid",fECADQID);
  }
  if(fECAflag==2){
    ECA_table = new DBTable("ECA_TSLP",fECAIndex);
    ECA_table->SetRunRange(fRunNo,fRunRangeMax);
    //pass = -2 means that this table will be pushed to RATDB in truncated mode
    ECA_table->SetPassNumber(-2);
    ECA_table->SetI("version",3);
    ECA_table->SetI("run_status",fOverallMask);
    // TSLP: 3 status words per cell, plus array of up to N points + 4 fit pars for slope
    ECA_table->SetIArray("tslp_status",fCellStatus);
    ECA_table->SetDArray("tslp_times",fEVTACX);
    ECA_table->SetI("tslp_ntimes",fNSlopePoints);
    std::vector<int> temp_int;
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNdim;++t){
        temp_int.push_back(fCellBad[n][t]);
      }
    }
    ECA_table->SetIArray("tslp_bad",temp_int);
    temp_int.clear();
    std::vector<double> temp;
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNSlopePoints;++t){
        temp.push_back(fCellTACPoints[n][t]);
      }
    }
    ECA_table->SetDArray("tslp_points",temp);
    temp.clear();
    //round doubles to reduce RATDB table sizes
    double precision [] = {0.1, 1e-3, 1e-6, 1e-9};
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNTACFitPars;++t){
        temp.push_back( floor (fCellTACFits[n][t]/precision[t] + .5)*precision[t] );
      }
    }
    ECA_table->SetDArray("tslp_fitpars",temp);
    temp.clear();
    ECA_table->SetIArray("tslp_dqid",fECADQID);

  } //endif tslps

  ECA_table->SaveAs(fECAdbname);

}


void ECAProc::CloseECADebug()
{
  ecadebugrecord<<"" <<newline;
  ecadebugrecord<<"" <<newline;
  ecadebugrecord<<"}" <<newline;
  fECAdebugdbfile.close();
}



void ECAProc::FillECADebug()
{
  if(fECAflag==0){
    ecadebugrecord<<"So this file is empty empty empty"<<newline;
    return;
  }
  if(fECAflag==1){
    ecadebugrecord<< "" << newline;
  }

  if(fECAflag==2){
    ecadebugrecord<<"\"tslp_parerrs\": [";
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNTACFitPars;++t){
        fprintf(fECAdebugdbfile, "%4.3f, ",fCellTACFitErrs[n][t]);
      }
    }
    ecadebugrecord<<"]," <<newline;
  }
}

void ECAProc::CloseECARoot()
{
  // Write hists to file
  // Close file
  fECArootfile->cd();
  fHists.SetAxes(fECAflag, fDoDroop);
  fHists.WriteHists(fECAflag, fECArootfile, fDoDroop);
  fECArootfile->Close();
}

void ECAProc::CreateECAMoni()
{

  RAT::DBTable *ECA_monitoring;
  if(fECAflag==0){
    ECA_monitoring = new DBTable("ECA_INVALID");
  }
  if(fECAflag==1){
    ECA_monitoring = new DBTable("ECA_MONI",fECAIndex);
    ECA_monitoring->SetRunRange(fRunNo,fRunRangeMax);
    //pass = -2 means that this table will be pushed to RATDB in truncated mode
    ECA_monitoring->SetPassNumber(-2);
    ECA_monitoring->SetI("version",1);
    ECA_monitoring->SetI("run_status",fOverallMask);
    ECA_monitoring->SetIArray("chan_tested",fTestChannel);

    // PDST: we want pedestals for 3Q+T, and status word, per CELL i.e. 16*channels
    ECA_monitoring->SetIArray("pdst_nevents",fNevCell);
    std::vector<int> temp_int;
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNMaxEV;++t){
        temp_int.push_back(fEVQHS[n][t]);
      }
    }
    ECA_monitoring->SetIArray("pdst_evt_qhs",temp_int);
    temp_int.clear();
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNMaxEV;++t){
        temp_int.push_back(fEVQHL[n][t]);
      }
    }
    ECA_monitoring->SetIArray("pdst_evt_qhl",temp_int);
    temp_int.clear();
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNMaxEV;++t){
        temp_int.push_back(fEVQLX[n][t]);
      }
    }
    ECA_monitoring->SetIArray("pdst_evt_qlx",temp_int);
    temp_int.clear();
    for(int n=0;n<fNTotCell;++n){
      for(int t = 0; t<fNMaxEV;++t){
        temp_int.push_back(fEVTAC[n][t]);
      }
    }
    ECA_monitoring->SetIArray("pdst_evt_tac",temp_int);
  } else if(fECAflag == 2){

    ECA_monitoring = new DBTable("ECA_MONI",fECAIndex);
    //Empty. No useful info to store at the moment.

  }

  ECA_monitoring->SaveAs(fECAmoniname);

}

// Merge the results we have from this run with previous results
// Identify cells that don't have results as those with value = 0
// Check that if one value is 0, they all are
// New: We now have a flag for channels we tested in this calibration
// i.e. fTestChannel[lcn] so check that too
// Scenarios:
// a) All values == 0 && fTestChannel == 0: empty, use previous
// b) No values == 0 && fTestChannel == 1: properly full, do nothing
// c) (Exclusive || of the above) || Some values == 0: problem...
void ECAProc::MergeECAResults()
{
  if(fECAflag==0){
    return;
  }
  if(fECAflag==1){
    for(int nch = 0; nch<fTotChannel;++nch){
      int testch = 0;
      for(int nce = 0; nce<fNcell;++nce){
        int n = fBits->GetCCCC(nch,nce);
        int test = 0;  // no values empty
        /*if(fCellStatus[n]==0){    //CANNOT TEST ON STATUS WORD SINCE IT COULD BE 0 ANYWAY
          fCellStatus[n]=fPrevStatus[n]; // So set it with the tests on other values instead
          ++test;
          }*/
        if(fCellQHS[n]==0.0){
          fCellQHS[n]=fPrevQHS[n];
          fCellSigQHS[n]=fPrevSigQHS[n];
          fCellStatus[n]=fPrevStatus[n];
          ++test;
          ++testch;
        }
        if(fCellQHL[n]==0.0){
          fCellQHL[n]=fPrevQHL[n];
          fCellSigQHL[n]=fPrevSigQHL[n];
          fCellStatus[n]=fPrevStatus[n];
          ++test;
          ++testch;
        }
        if(fCellQLX[n]==0.0){
          fCellQLX[n]=fPrevQLX[n];
          fCellSigQLX[n]=fPrevSigQLX[n];
          fCellStatus[n]=fPrevStatus[n];
          ++test;
          ++testch;
        }
        if(fCellTAC[n]==0.0){
          fCellTAC[n]=fPrevTAC[n];
          fCellSigTAC[n]=fPrevSigTAC[n];
          fCellStatus[n]=fPrevStatus[n];
          ++test;
          ++testch;
        }
        if(!(test==0 || test==4)){
          // Some values are 0, but not all == error??
          ECAProc::DontDie("", 1, n, test);
          // Set bit 31 in status mask, for a 'non-smooth' merge
          fOverallMask = fBits->SetBit(fOverallMask, 31);
          // Use values from previous ECA run for all quantities
          fCellStatus[n]=fPrevStatus[n];
          fCellQHS[n]=fPrevQHS[n];
          fCellQHL[n]=fPrevQHL[n];
          fCellQLX[n]=fPrevQLX[n];
          fCellTAC[n]=fPrevTAC[n];
          fCellSigQHS[n]=fPrevSigQHS[n];
          fCellSigQHL[n]=fPrevSigQHL[n];
          fCellSigQLX[n]=fPrevSigQLX[n];
          fCellSigTAC[n]=fPrevSigTAC[n];
        }
        if(test==4 && fTestChannel[nch]==1){
          // No filled values BUT we did test this channel!
          ECAProc::Die("", 5, nch, fTestChannel[nch]);
        }
        if(test==0 && fTestChannel[nch]==0){
          // All values filled, BUT we did NOT test this channel!
          ECAProc::Die("", 5, nch, fTestChannel[nch]);
        }
      }// !cells
      if(testch>0){
        //Any cells in this channel had some empty values
        //Replace DQID values with previous
        for(int nd=0;nd<5;++nd){
          int icrate = fBits->GetCratewCell(fBits->GetCCCC(nch,0));
          int icard  = fBits->GetCardwCell(fBits->GetCCCC(nch,0));
          int offset = 80*icrate + 5*icard + nd;
          fECADQID[offset] = fECAPrevDQID[offset];
        }
        if(testch!=4*fNcell){
          // only some cells were empty for this channel...
          ECAProc::DontDie("",5,nch,testch);
          // Set bit 31 in status mask, for a 'non-smooth' merge
          fOverallMask = fBits->SetBit(fOverallMask, 31);
        }
      } // !this channel had some empty cells
    } // !channels
  }
  if(fECAflag==2){
    bool domerge = ECAProc::MatchPrevTinj("merge");
    if(domerge){
      for(int nch = 0; nch<fTotChannel;++nch){
        int testch = 0;
        for(int nce = 0; nce<fNcell;++nce){
          int n = fBits->GetCCCC(nch,nce);

          int test = 0;  // no values empty
          // Can't test status words, since good words are 0 anyway
          // Test first point in TAC slope, and the first 2 fit pars
          if(fCellTACPoints[n][0]==0.0){
            ++test;
            ++testch;
          }
          if(fCellTACFits[n][0]==0.0){
            ++test;
            ++testch;
          }
          if(fCellTACFits[n][1]==0.0){
            ++test;
            ++testch;
          }
          if(!(test==0 || test==3)){
            // Some values are 0, but not all == error??
            ECAProc::DontDie("", 2, n, test);
            // Set bit 31 in status mask, for a 'non-smooth' merge
            fOverallMask = fBits->SetBit(fOverallMask, 31);
          }
          if(test==3 && fTestChannel[nch]==1){
            // No filled values BUT we did test this channel!
            ECAProc::Die("", 5, nch, fTestChannel[nch]);
          }
          if(test==0 && fTestChannel[nch]==0){
            // All values filled, BUT we did NOT test this channel!
            ECAProc::Die("", 5, nch, fTestChannel[nch]);
          }
          if(test>0){
            // Some values are zero, so use values from previous run for everything
            fCellStatus[n]=fPrevStatus[n];
            fCellBad[n]=fPrevBad[n];
            //fLastTAC[n] = fPrevLast[n];
            for(int t = 0; t<fNSlopePoints;++t){
              fCellTACPoints[n][t] = fPrevTACPoints[n][t];
            }
            for(int t = 0; t<fNTACFitPars;++t){
              fCellTACFits[n][t] = fPrevTACFits[n][t];
            }
          }
        } //!cells

        if(testch>0){
          //Any cells in this channel had some empty values
          //Replace DQID values with previous
          for(int nd=0;nd<5;++nd){
            int icrate = fBits->GetCratewCell(fBits->GetCCCC(nch,0));
            int icard  = fBits->GetCardwCell(fBits->GetCCCC(nch,0));
            int offset = 80*icrate + 5*icard + nd;
            fECADQID[offset] = fECAPrevDQID[offset];
          }
          if(testch!=3*fNcell){
            // only some cells were empty for this channel...
            ECAProc::DontDie("",6,nch,testch);
            // Set bit 31 in status mask, for a 'non-smooth' merge
            fOverallMask = fBits->SetBit(fOverallMask, 31);
          }
        } // !this channel had some empty cells

      } //!channels
    } // domerge check
  }
}



// Return the median and widths of the n values of vector arg
void ECAProc::MultiCalc(int n, std::vector<int> arg)
{
  fMedian = 0.0;
  fWidth  = 0.0;
  fWidth2 = 0.0;
  fSig = 0.0;
  if(n==0){return;}

  int *array = &arg[0]; //Javi: root wants an array

  if(n%2==1){  // Odd number, return middle value
    fMedian = TMath::KOrdStat(n, array, int(n/2));
  }
  else{        // Even number, return mean of middle 2 values
    fMedian = 0.5*(TMath::KOrdStat(n, array, int(n/2 - 1)) + TMath::KOrdStat(n, array, int(n/2)));
  }

// Return the width of the n values of vector arg
// Span around median containing 50% of points
  fWidth = TMath::KOrdStat(n, array, int(0.75*n)) - TMath::KOrdStat(n, array, int(0.25*n));
  if(fWidth == 0.0){
    fWidth = (1.0/1.95)*(TMath::KOrdStat(n, array, int(0.9*n)) - TMath::KOrdStat(n, array,int(0.1*n)));
  }
  if(fWidth == 0.0){
    fWidth = 1.0;
  }
  fWidth = int(10*fWidth*1.48/2.0) / 10.0;
  if(fWidth < 0.5) fWidth = 0.5;

// Return the 2nd width measure of the n values of vector arg
// Span around median containing 20% of points
  fWidth2 = TMath::KOrdStat(n, array, int(0.7*n)) - TMath::KOrdStat(n, array, int(0.4*n));

  // Estimate of 1 sigma width:
  // 40% of a gaussian's data lies inside 0.52*sigma of the range
  // plus for some reason estimating sigma from percentiles induces a small bias...(?)

  //fSig is invalid in this calculation of n < 4, so set those values to -1 and fix them later.
  //fSig is only used in Ped runs
  if( n < 4 ) fSig = -1.;
  else{
  fSig = TMath::KOrdStat(n, array, int(0.7*n)-1) - TMath::KOrdStat(n, array, int(0.3*n)-1);
  fSig = fSig * 1.0055 / (2.0*0.52);
  }

  //  delete[] array;

}

void ECAProc::ReassignSigs()
{

  for(int i=0; i<fNTotCell; i++){
    if(fCellSigQHS[i] == -1.) fCellsToReassign.push_back(i);
  }


  std::vector<float> qhslimits;
  qhslimits = fECAbank->GetFArrayFromD("qhs_width_ranges");
  std::vector<float> newqhssigs;
  newqhssigs = fECAbank->GetFArrayFromD("qhs_replacement_sigs");

  std::vector<float> qhllimits;
  qhllimits = fECAbank->GetFArrayFromD("qhl_width_ranges");
  std::vector<float> newqhlsigs;
  newqhlsigs = fECAbank->GetFArrayFromD("qhl_replacement_sigs");

  std::vector<float> qlxlimits;
  qlxlimits = fECAbank->GetFArrayFromD("qlx_width_ranges");
  std::vector<float> newqlxsigs;
  newqlxsigs = fECAbank->GetFArrayFromD("qlx_replacement_sigs");

  std::vector<float> taclimits;
  taclimits = fECAbank->GetFArrayFromD("tac_width_ranges");
  std::vector<float> newtacsigs;
  newtacsigs = fECAbank->GetFArrayFromD("tac_replacement_sigs");

  for(int i=0; i < (int) fCellsToReassign.size(); i++){
    int icell = fCellsToReassign[i];

    float test = fCellWidthQHS[icell];
    if((test >= 0.) && (test < qhslimits[0])) fCellSigQHS[icell] = newqhssigs[0];
    else if((test >= qhslimits[0]) && (test < qhslimits[1])) fCellSigQHS[icell] = newqhssigs[1];
    else if((test >= qhslimits[1]) && (test < qhslimits[2])) fCellSigQHS[icell] = newqhssigs[2];
    else if(test >= qhslimits[2]) fCellSigQHS[icell] = newqhssigs[3];
    else{
      std::cout << fCellWidthQHS[icell] << "  ahhhhh! " << std::endl;
    }

    test = fCellWidthQHL[icell];
    if((test >= 0.) && (test < qhllimits[0])) fCellSigQHL[icell] = newqhlsigs[0];
    else if((test >= qhllimits[0]) && (test < qhllimits[1])) fCellSigQHL[icell] = newqhlsigs[1];
    else if((test >= qhllimits[1]) && (test < qhllimits[2])) fCellSigQHL[icell] = newqhlsigs[2];
    else if(test >= qhllimits[2]) fCellSigQHL[icell] = newqhlsigs[0];
    else{
      std::cout << fCellWidthQHL[icell] << "  ahhhhh! " << std::endl;
    }

    test = fCellWidthQLX[icell];
    if((test >= 0.) && (test < qlxlimits[0])) fCellSigQLX[icell] = newqlxsigs[0];
    else if((test >= qlxlimits[0]) && (test < qlxlimits[1])) fCellSigQLX[icell] = newqlxsigs[1];
    else if((test >= qlxlimits[1]) && (test < qlxlimits[2])) fCellSigQLX[icell] = newqlxsigs[2];
    else if(test >= qlxlimits[2]) fCellSigQLX[icell] = newqlxsigs[0];
    else{
      std::cout << fCellWidthQLX[icell] << "  ahhhhh! " << std::endl;
    }

    test = fCellWidthTAC[icell];
    if((test >= 0.) && (test < taclimits[0])) fCellSigTAC[icell] = newtacsigs[0];
    else if((test >= taclimits[0]) && (test < taclimits[1])) fCellSigTAC[icell] = newtacsigs[1];
    else if(test >= taclimits[1]) fCellSigTAC[icell] = newtacsigs[0];
    else{
      std::cout << fCellWidthTAC[icell] << "  ahhhhh! " << std::endl;
    }
  }
}




bool ECAProc::FileExists(std::string filename)
{
  struct stat stFileInfo;
  bool open;
  int intStat;

  // Attempt to get the file attributes
  intStat = stat(filename.c_str(),&stFileInfo);
  if(intStat == 0) {
    open = true;
  } else {
    open = false;
  }

  return(open);
}

bool ECAProc::MatchPrevTinj(std::string msg)
{
  bool dotest = true;

  if(msg == "diff")msg = "test differences from";
  if(msg == "merge")msg = "merge with";

  if(fNSlopePoints != fPrevNSlopePoints){
    ECAProc::DontDie(msg,7,fPrevNSlopePoints,fNSlopePoints);
    dotest = false;
    return dotest;
  }

  for(int i=0;i<fNSlopePoints;++i){
    if(fEVTACX[i] != fPrevEVTACX[i]){
      ECAProc::DontDie(msg,7,fPrevEVTACX[i],fEVTACX[i]);
      dotest = false;
      return dotest;
    }
  }

  return dotest;
}


  //Load Patterns
void ECAProc::LoadPatterns(){

  fPatterns.resize(fNSet);
  for(int ni=0;ni<fNSet;++ni){
    fPatterns[ni].resize(fNPatt[ni]);
  }

  fPatterns[0][0] = fECAPatt->GetIArray("bb1");
  fPatterns[0][1] = fECAPatt->GetIArray("bb2");
  fPatterns[0][2] = fECAPatt->GetIArray("bb3");
  fPatterns[0][3] = fECAPatt->GetIArray("bb4");
  fPatterns[0][4] = fECAPatt->GetIArray("bb5");
  fPatterns[0][5] = fECAPatt->GetIArray("bb6");
  fPatterns[0][6] = fECAPatt->GetIArray("bb7");
  fPatterns[0][7] = fECAPatt->GetIArray("bb8");
  fPatterns[1][0] = fECAPatt->GetIArray("s1");
  fPatterns[1][1] = fECAPatt->GetIArray("s2");
  fPatterns[1][2] = fECAPatt->GetIArray("s3");
  fPatterns[1][3] = fECAPatt->GetIArray("s4");
  fPatterns[1][4] = fECAPatt->GetIArray("s5");
  fPatterns[1][5] = fECAPatt->GetIArray("s6");
  fPatterns[1][6] = fECAPatt->GetIArray("s7");
  fPatterns[1][7] = fECAPatt->GetIArray("s8");
  fPatterns[1][8] = fECAPatt->GetIArray("s9");
  fPatterns[1][9] = fECAPatt->GetIArray("s10");
  fPatterns[1][10] = fECAPatt->GetIArray("s11");
  fPatterns[1][11] = fECAPatt->GetIArray("s12");
  fPatterns[1][12] = fECAPatt->GetIArray("s13");
  fPatterns[1][13] = fECAPatt->GetIArray("s14");
  fPatterns[1][14] = fECAPatt->GetIArray("s15");
  fPatterns[1][15] = fECAPatt->GetIArray("s16");
  fPatterns[1][16] = fECAPatt->GetIArray("s17");
  fPatterns[1][17] = fECAPatt->GetIArray("s18");
  fPatterns[1][18] = fECAPatt->GetIArray("s19");
  fPatterns[1][19] = fECAPatt->GetIArray("s20");
  fPatterns[1][20] = fECAPatt->GetIArray("s21");
  fPatterns[1][21] = fECAPatt->GetIArray("s22");
  fPatterns[1][22] = fECAPatt->GetIArray("s23");
  fPatterns[1][23] = fECAPatt->GetIArray("s24");
  fPatterns[1][24] = fECAPatt->GetIArray("s25");
  fPatterns[1][25] = fECAPatt->GetIArray("s26");
  fPatterns[1][26] = fECAPatt->GetIArray("s27");
  fPatterns[1][27] = fECAPatt->GetIArray("s28");
  fPatterns[1][28] = fECAPatt->GetIArray("s29");
  fPatterns[1][29] = fECAPatt->GetIArray("s30");
  fPatterns[1][30] = fECAPatt->GetIArray("s31");
  fPatterns[1][31] = fECAPatt->GetIArray("s32");
  fPatterns[2][0] = fECAPatt->GetIArray("d1");
  fPatterns[2][1] = fECAPatt->GetIArray("d2");
  fPatterns[2][2] = fECAPatt->GetIArray("d3");
  fPatterns[2][3] = fECAPatt->GetIArray("d4");
  fPatterns[2][4] = fECAPatt->GetIArray("d5");
  fPatterns[2][5] = fECAPatt->GetIArray("d6");
  fPatterns[2][6] = fECAPatt->GetIArray("d7");
  fPatterns[2][7] = fECAPatt->GetIArray("d8");
  fPatterns[2][8] = fECAPatt->GetIArray("d9");
  fPatterns[2][9] = fECAPatt->GetIArray("d10");
  fPatterns[2][10] = fECAPatt->GetIArray("d11");
  fPatterns[2][11] = fECAPatt->GetIArray("d12");
  fPatterns[2][12] = fECAPatt->GetIArray("d13");
  fPatterns[2][13] = fECAPatt->GetIArray("d14");
  fPatterns[2][14] = fECAPatt->GetIArray("d15");
  fPatterns[2][15] = fECAPatt->GetIArray("d16");
  fPatterns[2][16] = fECAPatt->GetIArray("d17");
  fPatterns[2][17] = fECAPatt->GetIArray("d18");
  fPatterns[2][18] = fECAPatt->GetIArray("d19");
  fPatterns[3][0] = fECAPatt->GetIArray("p1");
  fPatterns[3][1] = fECAPatt->GetIArray("p2");
  fPatterns[3][2] = fECAPatt->GetIArray("p3");
  fPatterns[3][3] = fECAPatt->GetIArray("p4");
  fPatterns[3][4] = fECAPatt->GetIArray("p5");
  fPatterns[3][5] = fECAPatt->GetIArray("p6");
  fPatterns[3][6] = fECAPatt->GetIArray("p7");
  fPatterns[3][7] = fECAPatt->GetIArray("p8");
  fPatterns[3][8] = fECAPatt->GetIArray("p9");
  //kate adding channel by channel patterns
  fPatterns[4][0] = fECAPatt->GetIArray("ch0");
  fPatterns[4][1] = fECAPatt->GetIArray("ch1");
  fPatterns[4][2] = fECAPatt->GetIArray("ch2");
  fPatterns[4][3] = fECAPatt->GetIArray("ch3");
  fPatterns[4][4] = fECAPatt->GetIArray("ch4");
  fPatterns[4][5] = fECAPatt->GetIArray("ch5");
  fPatterns[4][6] = fECAPatt->GetIArray("ch6");
  fPatterns[4][7] = fECAPatt->GetIArray("ch7");
  fPatterns[4][8] = fECAPatt->GetIArray("ch8");
  fPatterns[4][9] = fECAPatt->GetIArray("ch9");
  fPatterns[4][10] = fECAPatt->GetIArray("ch10");
  fPatterns[4][11] = fECAPatt->GetIArray("ch11");
  fPatterns[4][12] = fECAPatt->GetIArray("ch12");
  fPatterns[4][13] = fECAPatt->GetIArray("ch13");
  fPatterns[4][14] = fECAPatt->GetIArray("ch14");
  fPatterns[4][15] = fECAPatt->GetIArray("ch15");
  fPatterns[4][16] = fECAPatt->GetIArray("ch16");
  fPatterns[4][17] = fECAPatt->GetIArray("ch17");
  fPatterns[4][18] = fECAPatt->GetIArray("ch18");
  fPatterns[4][19] = fECAPatt->GetIArray("ch19");
  fPatterns[4][20] = fECAPatt->GetIArray("ch20");
  fPatterns[4][21] = fECAPatt->GetIArray("ch21");
  fPatterns[4][22] = fECAPatt->GetIArray("ch22");
  fPatterns[4][23] = fECAPatt->GetIArray("ch23");
  fPatterns[4][24] = fECAPatt->GetIArray("ch24");
  fPatterns[4][25] = fECAPatt->GetIArray("ch25");
  fPatterns[4][26] = fECAPatt->GetIArray("ch26");
  fPatterns[4][27] = fECAPatt->GetIArray("ch27");
  fPatterns[4][28] = fECAPatt->GetIArray("ch28");
  fPatterns[4][29] = fECAPatt->GetIArray("ch29");
  fPatterns[4][30] = fECAPatt->GetIArray("ch30");
  fPatterns[4][31] = fECAPatt->GetIArray("ch31");
  //kate - adding half-crate masks
  fPatterns[5][0] = fECAPatt->GetIArray("hc1");
  fPatterns[5][1] = fECAPatt->GetIArray("hc2");
  fPatterns[5][2] = fECAPatt->GetIArray("hc3");
  fPatterns[5][3] = fECAPatt->GetIArray("hc4");
  fPatterns[5][4] = fECAPatt->GetIArray("hc5");
  fPatterns[5][5] = fECAPatt->GetIArray("hc6");
  fPatterns[5][6] = fECAPatt->GetIArray("hc7");
  fPatterns[5][7] = fECAPatt->GetIArray("hc8");
  fPatterns[5][8] = fECAPatt->GetIArray("hc9");
  fPatterns[5][9] = fECAPatt->GetIArray("hc10");
  fPatterns[5][10] = fECAPatt->GetIArray("hc11");
  fPatterns[5][11] = fECAPatt->GetIArray("hc12");
  fPatterns[5][12] = fECAPatt->GetIArray("hc13");
  fPatterns[5][13] = fECAPatt->GetIArray("hc14");
  fPatterns[5][14] = fECAPatt->GetIArray("hc15");
  fPatterns[5][15] = fECAPatt->GetIArray("hc16");
  fPatterns[5][16] = fECAPatt->GetIArray("hc17");
  fPatterns[5][17] = fECAPatt->GetIArray("hc18");
  fPatterns[5][18] = fECAPatt->GetIArray("hc19");
  fPatterns[5][19] = fECAPatt->GetIArray("hc20");
  fPatterns[5][20] = fECAPatt->GetIArray("hc21");
  fPatterns[5][21] = fECAPatt->GetIArray("hc22");
  fPatterns[5][22] = fECAPatt->GetIArray("hc23");
  fPatterns[5][23] = fECAPatt->GetIArray("hc24");
  fPatterns[5][24] = fECAPatt->GetIArray("hc25");
  fPatterns[5][25] = fECAPatt->GetIArray("hc26");
  fPatterns[5][26] = fECAPatt->GetIArray("hc27");
  fPatterns[5][27] = fECAPatt->GetIArray("hc28");
  fPatterns[5][28] = fECAPatt->GetIArray("hc29");
  fPatterns[5][29] = fECAPatt->GetIArray("hc30");
  fPatterns[5][30] = fECAPatt->GetIArray("hc31");
  fPatterns[5][31] = fECAPatt->GetIArray("hc32");
  fPatterns[5][32] = fECAPatt->GetIArray("hc33");
  fPatterns[5][33] = fECAPatt->GetIArray("hc34");
  fPatterns[5][34] = fECAPatt->GetIArray("hc35");
  fPatterns[5][35] = fECAPatt->GetIArray("hc36");
  fPatterns[5][36] = fECAPatt->GetIArray("hc37");
  fPatterns[5][37] = fECAPatt->GetIArray("hc38");

}





// Handle errors
void ECAProc::DontDie(std::string message, int flag_err, float info1, float info2)
{

  if (fECAlogverb == 0) return;

  if(flag_err < 8) warn << BYELLOW << " ECA SETBACK: " << CLR << newline;
  ecalog << "----- ECA SETBACK: -----" <<newline;

  if(flag_err==1){  // Problem in the merging of old and new files: PDST
    int lcn = fBits->GetLCN(static_cast<int>(info1));
    int icell = fBits->GetCell(static_cast<int>(info1));
    ecalog << "  ECA for LCN "<<lcn<<", cell "<<icell;
    ecalog <<"  only has non-zero values for "<<newline;
    ecalog <<"  "<<  info2 <<" of the 4 quantities (QHS, QHL, QLX, TAC)"<<newline;
    ecalog <<"  => Using values from previous ECA run"<< newline;
    ecalog <<"  and setting bit 31 in the ECA run-status word"<< newline;
    ecalog << newline;
    warn << "  ECA for LCN "<<lcn<<", cell "<<icell;
    warn <<"  only has non-zero values for "<<newline;
    warn <<"  "<< info2 <<" of the 4 quantities (QHS, QHL, QLX, TAC)"<<newline;
    warn <<"  => Using values from previous ECA run"<< newline;
    warn <<"  and setting bit 31 in the ECA run-status word"<< newline;
    warn << newline;
  }

  if(flag_err==2){  // Problem in the merging of old and new files: TSLP
    int lcn = fBits->GetLCN(static_cast<int>(info1));
    int icell = fBits->GetCell(static_cast<int>(info1));
    ecalog << "  ECA for LCN "<<lcn<<", cell "<<icell;
    ecalog <<"  only has non-zero values for "<<newline;
    ecalog <<"  "<<  info2 <<" of the 3 quantities checked"<<newline;
    ecalog <<"  (first point on TAC slope + 0th and 1st order fit pars)"<<newline;
    ecalog <<"  => Using values from previous ECA run"<< newline;
    ecalog <<"  and setting bit 31 in the ECA run-status word"<< newline;
    ecalog << newline;
    warn << "  ECA for LCN "<<lcn<<", cell "<<icell;
    warn <<"  only has non-zero values for "<<newline;
    warn <<"  "<<  info2 <<" of the 3 quantities checked"<<newline;
    warn <<"  (first point on TAC slope + 0th and 1st order fit pars)"<<newline;
    warn <<"  => Using values from previous ECA run"<< newline;
    warn <<"  and setting bit 31 in the ECA run-status word"<< newline;
    warn << newline;
  }


  if(flag_err==4){  // Number of TAC points for this cell != expectation, skip it
    ecalog << "  Number of TAC points for cell "<<info1<<" is not what we expect"<<newline;
    ecalog << "  Npoints = "<<info2<<" where we expect "<<fNSlopePoints<<newline;
//    ecalog << "  Skipping this cell.";
    ecalog << newline;
    warn << "  Number of TAC points for cell "<<info1<<" is not what we expect"<<newline;
    warn << "  Npoints = "<<info2<<" where we expect "<<fNSlopePoints<<newline;
//    warn << "  Skipping this cell.";
    warn << newline;
  }

  if(flag_err==5){  // Problem in the merging of old and new files: PDST
    ecalog << "  ECA for LCN "<<info1;
    ecalog <<"  only has non-zero values for "<<info2<<newline;
    ecalog <<"  of the 4 quantities (QHS, QHL, QLX, TAC) over "<<fNcell<<" cells"<<newline;
    ecalog <<"  => Using values from previous ECA run for board IDs"<< newline;
    ecalog <<"  and setting bit 31 in the ECA run-status word"<< newline;
    ecalog << newline;
    warn << "  ECA for LCN "<<info1;
    warn <<"  only has non-zero values for "<<info2<<newline;
    warn <<"  of the 4 quantities (QHS, QHL, QLX, TAC) over "<<fNcell<<" cells"<<newline;
    warn <<"  => Using values from previous ECA run for board IDs"<< newline;
    warn <<"  and setting bit 31 in the ECA run-status word"<< newline;
    warn << newline;
  }

  if(flag_err==6){  // Problem in the merging of old and new files: TSLP
    ecalog << "  ECA for LCN "<<info1;
    ecalog <<"  only has non-zero values for "<<info2<<newline;
    ecalog <<"  of the 3 quantities checked over "<<fNcell<<" cells"<<newline;
    ecalog <<"  (first point on TAC slope + 0th and 1st order fit pars)"<<newline;
    ecalog <<"  => Using values from previous ECA run for board IDs"<< newline;
    ecalog <<"  and setting bit 31 in the ECA run-status word"<< newline;
    ecalog << newline;
    warn << "  ECA for LCN "<<info1;
    warn <<"  only has non-zero values for "<<info2<<newline;
    warn <<"  of the 3 quantities checked over "<<fNcell<<" cells"<<newline;
    warn <<"  (first point on TAC slope + 0th and 1st order fit pars)"<<newline;
    warn <<"  => Using values from previous ECA run for board IDs"<< newline;
    warn <<"  and setting bit 31 in the ECA run-status word"<< newline;
    warn << newline;
  }
  if(flag_err==7){  // EVTacX doesn't match previous
    ecalog << "  You have asked me to "<<message<<" previous results"<<newline;
    ecalog << "  But the previous injected T values do not match current ones:"<<newline;
    ecalog << "  Previous/current number of steps / Tinj value = "<<info1<<" / "<<info2<<newline;
    ecalog << "  So I cannot do the requested task."<<newline;
    ecalog << "  Continuing with ECAProc."<<newline;
    warn << "  You have asked me to "<<message<<" previous results"<<newline;
    warn << "  But the previous injected T values do not match current ones:"<<newline;
    warn << "  Previous/current number of steps / Tinj value = "<<info1<<" / "<<info2<<newline;
    warn << "  So I cannot do the requested task."<<newline;
    warn << "  Continuing with ECAProc."<<newline;
  }

  if(flag_err==8){ //ecatype is set to 0
    ecalog << "  CalType set to non-ECA type:" << info1 <<newline;
    ecalog << "  Skipping this event..."<<newline;
    ecalog<<newline;
  }

  if(flag_err==9){ //no header info
    if(info1) ecalog << "  This event has no HeaderInfo" <<newline;
    if(info2) ecalog << "  This event has no EPEDInfo" <<newline;
    ecalog << "  Skipping this event..."<<newline;
    ecalog<<newline;
  }
  if(flag_err==10){ //Tinj out of order
    ecalog << "   This event has time out of order, skipping..." << newline;
  }
}

// Handle FATAL errors
void ECAProc::Die(std::string message, int flag_err, float info1, float info2, int return_code)
{
  warn << RED << " ECA PROBLEM: " << CLR << newline;
  ecalog << "----- ECA PROBLEM: -----" <<newline;
  if(flag_err==0){
    ecalog<< message << newline;
    warn<< message << newline;
  }

  if(flag_err==1){  // We couldn't find a previous run to merge our results with.  Poo.
    ecalog << "  You have asked me to load results from the previous run"<<newline;
    ecalog << "  but I cannot find a file to load!"<<newline;
    ecalog << newline;
    warn << "  You have asked me to load results from the previous run"<<newline;
    warn << "  but I cannot find a file to load!"<<newline;
    warn << newline;
  }

  if(flag_err==2){  // No events in cell, but median not set to -9999.0 == CODE BUG
    ecalog << "  There are no events in this cell"<<newline;
    ecalog << "  and yet the median is not set to -9999.0!"<<newline;
    ecalog << "  It must be a bug in the code!!"<<newline;
    ecalog<<newline;
    warn << "  There are no events in this cell"<<newline;
    warn << "  and yet the median is not set to -9999.0!"<<newline;
    warn << "  It must be a bug in the code!!"<<newline;
    warn <<newline;
  }

  if(flag_err==3){  // Length of DQID array is not what we expect it to be
    ecalog << "  We have the wrong number of elements in DQID:"<<newline;
    ecalog << "  "<<info1<<" entries and we expect "<<info2<<newline;
    ecalog << "  Dying miserably now!!"<<newline;
    ecalog<<newline;
    warn << "  We have the wrong number of elements in DQID:"<<newline;
    warn << "  "<<info1<<" entries and we expect "<<info2<<newline;
    warn << "  Dying miserably now!!"<<newline;
    warn<<newline;
  }

  if(flag_err==4){  // MisCalc # channels with bad pedestals
    ecalog << "  You are calculating the # channels with bad peds wrongly!"<<newline;
    ecalog << "  I get "<<info1<<" and snoman-method gets "<<info2<<newline;
    ecalog << "  Dying miserably now!!"<<newline;
    ecalog<<newline;
    warn << "  You are calculating the # channels with bad peds wrongly!"<<newline;
    warn << "  I get "<<info1<<" and snoman-method gets "<<info2<<newline;
    warn << "  Dying miserably now!!"<<newline;
    warn<<newline;
  }

  if(flag_err==5){  // Error: Discrepancy between channels we tested, and values calc
    ecalog << "  Error for channel "<<info1<<newline;
    if(info2==0){
      ecalog << "  This channel was NOT tested and yet you have non-zero values!"<<newline;
    }
    if(info2==1){
      ecalog << "  This channel WAS tested and yet you have some zero values!"<<newline;
    }
    ecalog << "  Dying miserably now!!"<<newline;
    ecalog<<newline;
    warn << "  Error for channel "<<info1<<newline;
    if(info2==0){
      warn << "  This channel was NOT tested and yet you have non-zero values!"<<newline;
    }
    if(info2==1){
      warn << "  This channel WAS tested and yet you have some zero values!"<<newline;
    }
    warn << "  Dying miserably now!!"<<newline;
    warn<<newline;
  }

  if(flag_err==6){  // MisCalc # channels with bad tslope fits
    ecalog << "  You are calculating the # channels with bad fits wrongly!"<<newline;
    ecalog << "  I get "<<info1<<" and snoman-method gets "<<info2<<newline;
    ecalog << "  Dying miserably now!!"<<newline;
    ecalog<<newline;
    warn << "  You are calculating the # channels with bad fits wrongly!"<<newline;
    warn << "  I get "<<info1<<" and snoman-method gets "<<info2<<newline;
    warn << "  Dying miserably now!!"<<newline;
    warn<<newline;
  }

  if(flag_err==7){  // MisCalc whether channel is active
    ecalog << "  You are calculating whether channel "<<info1<<" is enabled INCORRECTLY"<<newline;
    ecalog << "  Calcs from DQCH/CR and DQXX disagree.  Check."<<newline;
    ecalog << "  Dying miserably now!!"<<newline;
    ecalog<<newline;
    warn << "  You are calculating whether channel "<<info1<<" is enabled INCORRECTLY"<<newline;
    warn << "  Calcs from DQCH/CR and DQXX disagree.  Check."<<newline;
    warn << "  Dying miserably now!!"<<newline;
    warn<<newline;
  }

  if(flag_err==8){  // Disagreement in calib type
    ecalog << "  Discrepancy in calibration type:"<<newline;
    ecalog << "  EPED Info thinks this is type "<<info1<<newline;
    ecalog << "  But you requested type "<<info2<<newline;
    ecalog<<newline;
    warn << "  Discrepancy in calibration type:"<<newline;
    warn << "  EPED Info thinks this is type "<<info1<<newline;
    warn << "  But you requested type "<<info2<<newline;
    warn<<newline;
  }

  if(flag_err==9){  // Tinj points out of order
    ecalog << "  Injected times are out of order:"<<newline;
    ecalog << "  Current EPED gives Tinj = "<<info1<<" but the prev Tinj was "<<info2<<newline;
    ecalog << "  This is a problem."<<info2<<newline;
    ecalog<<newline;
    warn << "  Injected times are out of order:"<<newline;
    warn << "  Current EPED gives Tinj = "<<info1<<" but the prev Tinj was "<<info2<<newline;
    warn << "  This is a problem."<<info2<<newline;
    warn<<newline;
  }

  if(flag_err==10){  // No EPED to tell us Tslopes
    ecalog << "  I have got to the first ECA event without an EPED to tell me the injected T"<<newline;
    ecalog << "  This is a problem."<<info2<<newline;
    ecalog<<newline;
    warn << "  I have got to the first ECA event without an EPED to tell me the injected T"<<newline;
    warn << "  This is a problem."<<info2<<newline;
    warn<<newline;
  }

  if(flag_err==11){  // Pattern Set or Pattern Number are out of range
    ecalog << "  PatternSet or Pattern Number was out of range "<<newline;
    ecalog << "  You called PattSet " << info1 << " and Pattern Number " << info2 << newline;
    ecalog << "  Which means you tried to index fPatterns[" << info1-1 << "][" << info2 << "]" << newline;
    ecalog<<newline;
    warn << "  PatternSet or Pattern Number was out of range "<<newline;
    warn << "  You called PattSet " << info1 << " and Pattern Number " << info2 << newline;
    warn << "  Which means you tried to index fPatterns[" << info1-1 << "][" << info2-1 << "]" << newline;
    warn<<newline;
  }

  if(flag_err==12){  // Previous run has diff number TSlope point, can't merge / diff results.
    ecalog << "  You have asked me to load results from the previous TSlope run"<<newline;
    ecalog << "  but this run has a different number of TSlope points."<<newline;
    ecalog << "  Please re-run without the requested merge/diff."<<newline;
    ecalog << newline;
    warn << "  You have asked me to load results from the previous TSlope run"<<newline;
    warn << "  but this run has a different number of TSlope points."<<newline;
    warn << "  Please re-run without the requested merge/diff."<<newline;
    warn << newline;
  }

  exit(return_code);

}
} // namespace RAT