// PCACal.cc
// Contact person: Freija Descamps <fbdescamps@lbl.gov>
// See PCACal.hh for more details
//———————————————————————//

#include <G4UImessenger.hh>

#include <TFile.h>
#include <TMath.h>

#include <RAT/PCACal.hh>
#include <RAT/PCABits.hh>
#include <RAT/DB.hh>
#include <RAT/DS/DataNotFound.hh>
#include <RAT/Log.hh>

#include <iostream>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <Randomize.hh>

using namespace std;
namespace RAT {

  // Time is set to DS::INVALID if calibration fails
  const double kFailValue = DS::INVALID;
  const double kFailCheck = DS::INVALID;
  const int kPCABitOffset = 32;

  PCACal::PCACal()
  {
    // Make sure these are all zero
    fQHS = 0.;
    fQHL = 0.;
    fQLX = 0.;
    fTime = 0.;
    // Set all Booleans
    fTimeWalkResultOK=false;
    fGainFitresultOK=false;
    fDoTimeWalk = false;
    fDoGainFit = false;
  }

  PCACal::~PCACal()
  {
  }

  void PCACal::BeginOfRun( DS::Run& run )
  {
    PCABits pcabits;
    pcabits.BeginOfRun(run);
    // Set up the necessary links to the database
    fCalbank = DB::Get()->GetLink("PMTCAL");
    // The requested PCA calibration
    fPCAtype = fCalbank->GetI("pca_type");
    // The requested tolerance level
    fTolLevel = fCalbank->GetI("pca_validation");
    // These are the TimeWalk PCA calibration constants
    fTimeWalkBank = DB::Get()->GetLink("PCA_TW");
    //Check if we are calibrating MC
    fIsMC = run.GetMCFlag();
    if(fIsMC) info<<" PCACal::BeginOfRun: Analyzing MC. PCA calibration won't be applied. "<<newline;

    // Find out which type of calibration is required
    if(fPCAtype%10 == 1) fDoTimeWalk = true;
    if((static_cast<int>((fPCAtype%100)/10) == 1))
      {
        fDoGainFit = true; // No npe is extracted at this time
        PCACal::DontDie(GAINFIT,0,0);
      }

    if(fDoTimeWalk)
      {
        // Find out how many interpolation points there are
        fPCANpointsTimeWalk=fTimeWalkBank->GetI("tw_npoints");
        // Read in the array with the timewalk interpolation points
        fTimeWalkPoints=fTimeWalkBank->GetIArray("twinter");
        // Read in the array with the timewalk status words
        fTimeWalkStatus=fTimeWalkBank->GetIArray("PCATW_status");
        // Find out if these are SNO calibration constants
        fIsSNO=fTimeWalkBank->GetI("is_sno");
        if( fIsSNO )
          {
            // Warn the user that we are using SNO constants
            PCACal::DontDie(SNOCONST,0,0);
          }
        // In order to save time during the calibration, extract all the
        // interpolation points at this point and store them! Trust me,
        // this makes everything much faster.
        // Loop through all channels and extract the interpolation points and the
        // fit parameters.
        for( unsigned i=0; i<fTimeWalkStatus.size(); i++ )
          {
            int locWord=0;
            int packResult=0;
            double charge=0.;
            double time=0.;
            double gradient=0.;
            double intercept=0.;
            for( int j=0; j<fPCANpointsTimeWalk; j++)
              {
                locWord=i*(fPCANpointsTimeWalk+1)+j;
                packResult=pcabits.PCAPacker(fPack.UNPACK,
                                              fPack.QTINTERPOLATION,
                                              fTimeWalkPoints[locWord],
                                              charge,
                                              time);
                if(packResult == 0)
                  {
                    fTimeWalkCharge[i].push_back(charge);
                    fTimeWalkTime[i].push_back(time);
                  }
              }

            locWord=i*(fPCANpointsTimeWalk+1)+fPCANpointsTimeWalk;
            packResult=pcabits.PCAPacker(fPack.UNPACK,
                                          fPack.GRADIENT_INTERCEPT,
                                          fTimeWalkPoints[locWord],
                                          intercept,
                                          gradient);
            if(packResult == 0)
              {
                fTimeWalkIntercept[i]=intercept;
                fTimeWalkGradient[i]=gradient;
              }
          }
      }
  }

  DS::PMTCal PCACal::CalibratePMT( const DS::PMTCal& pmt )
  {

    fStatus = pmt.GetStatus(); // Add to the status, not replace
    CalibrateHit( pmt.GetQHS(), pmt.GetQHL(), pmt.GetQLX(), pmt.GetTime(), pmt.GetID() );

    // Undo PCA calibration if MC. In this way we do not
    // perform the PCA calibration (which would be incorrect,
    // provided the PCA Uncal is not in place yet) but the
    // PCAStatus word gets correctly set for MC.
    if(fIsMC) fTime = pmt.GetTime();

    // Create new PMTCal and fill it
    DS::PMTCal newHit;
    newHit.SetID( pmt.GetID() );
    newHit.SetQHS( GetQHS() );
    newHit.SetQHL( GetQHL() );
    newHit.SetQLX( GetQLX() );
    newHit.SetTime( GetTime() );
    newHit.SetCellID( pmt.GetCellID() );
    newHit.SetChanFlags( pmt.GetChanFlags() );
    newHit.SetCrossTalkFlag( pmt.GetCrossTalkFlag() );
    newHit.SetStatus( GetStatus() );
    return newHit;
  }

  void PCACal::CalibrateHit(float qhs, float qhl, float qlx,
                            double time, int lcn)
  {

    fQHS = qhs;
    fQHL = qhl;
    fQLX = qlx;
    fTime = time;
    if(fPCAtype == 0)
      {
        // We really should never get here... Just in case though...
        PCACal::DontDie(NOPCA,fPCAtype,0);
        return;
      }
    if( fDoTimeWalk )
      {
        // Check if this is a valid QHS/T point
        if( fQHS == DS::INVALID && fTime==DS::INVALID )
          {
            // Invalid QHS/T point, probably due to failed ECA
            // calibration. For the moment, just leave as-is.
            // Mark the PCA as failed though in the status word
            // Bit 0 marks PCA failure.
            fStatus.Set( kPCABitOffset + 0);
            return;
          }
        else
          {
            // Perform the PCA time walk correction to the time
            PCACal::TimeWalkCal(lcn);
          }
      }

    if( fDoGainFit )
      {
        // Nothing here so far, charges are unaffected by PCA
        // This can change to calculating the QHS/HHP.. which
        // is a measure of npe.
      }

    // After calibration is done, check how well we did
    PCACal::ValidateCalibrate(lcn);

    // Call the tolerance level function only if needed.
    if( fTolLevel != 0 || fIsSNO )
      {

        // The tolerance level defines when a hit will be
        // set to (T,QHS) (DS::INVALID,QHS) on top of the failed
        // time-walk calibrations defined in PCACAL::TimeWalkCal(lcn)
        PCACal::CalStatus();
      }
    return;
  }

  // TimeWalk: timewalk correction, inspired by the SNO TimeWalk code
  void PCACal::TimeWalkCal(const int lcn)
  {
    // Set up local variables
    double thisq=0.0;
    double thist=0.0;
    double nextq=0.0;
    double nextt=0.0;
    double qstep=0.0;
    double twalk=0.0;

    // Check the first interpolation point for this channel
    // FYI: For now, the assumption is that, if the first point is undefined,
    // that all interpolation points will be. Assumptions are never good...
    // in the near future, we need to check if there are any channels that
    // have the first interpolation point undefined but not the subsequent ones.
    if( fTimeWalkCharge[lcn].empty() || fTimeWalkTime[lcn].empty() )
      {
        // No interpolation points available -> Do not correct by delays
        // For this case PCA time = ECA time

        // Flag this as a failed hit due to no TimeWalk calibration data
        fStatus.Set(kPCABitOffset + 2);
        PCACal::DontDie(FIRSTINT,lcn,0);
        return; // there is no calibration data available
      }
    // Set thisq to the first interpolation point
    thisq=fTimeWalkCharge[lcn][0];
    thist=fTimeWalkTime[lcn][0];
    // If the first interpolation point is defined, go ahead and
    // retrieve it!
    if( fTimeWalkCharge[lcn].size()<2 || fTimeWalkTime[lcn].size()<2 )
      {

        // No interpolation points available -> Do not correct by delays
        // For this case PCA time = ECA time

        // Flag this as a failed hit: only 1 calibration data-point available
        fStatus.Set(kPCABitOffset + 3);
        PCACal::DontDie(SECONDINT,lcn,0);
        // One single interpolation point is not enough
        // to perform a time calibration
        return;
      }
    nextq=fTimeWalkCharge[lcn][1];
    nextt=fTimeWalkTime[lcn][1];

    // We now have two valid interpolation points, there is much rejoicing
    // Define qstep as the difference between these two points
    qstep=nextq-thisq; // Distance in QHS between the two interpolation points

    /* Calculate TW effect based on PCA constants. There are 3 cases:
     * 1. QHS too low (<thisq-qstep): this is probably noise and we provide
     * our best guess, which is probably wrong but will account on cable delays properly.
     * 2. QHS low (<thisq): extrapolate back from the two first points
     * 3. QHS within reasonable values: interpolate if QHS is within range or
     *extrapolate using the TW gradient if QHS>>
     */
    if( fQHS<thisq-qstep )
      {
        // fQHS is too far below the first interpolation point,
        // it is probably noise so we don't actually now what is
        // the Q-T relation due to different pulse shape nature.
        // For this case provide our best guess given by the TW at -qstep
        twalk=thist-(nextt-thist);
        // Flag this as a failed hit. Responsibility to remove this is left
        // to the analyzer
        fStatus.Set(kPCABitOffset + 4);
        PCACal::DontDie(LOWQHS,lcn,fQHS);
      }
    // Now if the T/QHS point is situated below thisq, go
    // ahead and extract the TimeWalk correction for this point.
    // Also allow the case where we are sitting on top of
    // an interpolation point
    else if(fQHS<=thisq)
      {
        // Extrapolate backwards
        // This is the actual time correction:
        twalk=thist+(nextt-thist)*(fQHS-thisq)/qstep;
      }

    // If the T/QHS point is above the first interpolation point
    // we need to enter interpolation mode... also known as the
    // PCA loop from hell!
    else
      {

        for( int i=2;i<fPCANpointsTimeWalk;i++ )
          {
            // We start by using linear interpolation to calculate
            // twalk and then figure out next if this is the twalk
            // to be used.
            twalk=thist+(nextt-thist)*(fQHS-thisq)/(nextq-thisq); // 'this twalk'
            if(fQHS<=nextq)
              {
                // Allow also the case on interpolation point

                break; // Keep this twalk

                // else: get the next interpolation point and try again
              }
            // Only get the next interpolation point if we have not
            // reached the end of the available interpolation points.
            // Check if these next points are valid. It is confusing, but if the
            // interpolation is valid, the next if statement will be true and we'll
            // go back to the start of this loop, where the twalk is calculated and
            // we check if the T/QHS point is below the qnext point.
            if(i != fPCANpointsTimeWalk-1 &&
               fTimeWalkCharge[lcn].size() >= static_cast<unsigned int>(i+1) &&
               fTimeWalkTime[lcn].size() >= static_cast<unsigned int>(i+1)){
              // we have not yet exhausted the available interpolation points
              thisq=nextq;
              thist=nextt;
              nextq=fTimeWalkCharge[lcn][i+1];
              nextt=fTimeWalkTime[lcn][i+1];
            }
            // This loop is only entered if no more valid interpolation points are available
            else
              {
              // The next point is invalid. Therefore, in the current scheme
              // we will assume that there are no more valid interpolation points
                if( fQHS<thisq+qstep/2.0 )
                  {
                    // Close enough to last point so use the extrapolation (this twalk)
                    break; // twalk found
                  }
                else
                  {
                    // use the gradient and intercept if available
                    if( fTimeWalkIntercept.count(lcn)>0 && fTimeWalkGradient.count(lcn)>0 )
                      {
                        double gradient=fTimeWalkGradient[lcn];
                        double intercept=fTimeWalkIntercept[lcn];

                        // Now check the sanity of the intercept and gradient
                        if( gradient>0 || intercept==0 )
                          {
                            // It is pathological behaviour that the time
                            // correction would go up with charge,
                            // flag and use twalk=thist;
                            fStatus.Set(kPCABitOffset + 6);
                            twalk=thist; // Overwrite with the last time point
                            break; // twalk found
                          }
                        if( gradient<-0.10 )
                          {
                            // FIXME when real data comes in,
                            // set this limit correctly.
                            // It is pathological behaviour that the time
                            // correction would go down so quickly with charge,
                            // flag and use twalk=thist
                            fStatus.Set(kPCABitOffset + 7);
                            twalk=thist;
                            break; // twalk found
                          }
                        // The correction in time to be applied to data-point:
                        twalk = intercept+gradient*fQHS;
                        break; // twalk found
                      }
                    else
                      {
                        // Else: last resort! There is no fit to the last
                        // interpolation points saved. Use the last stored time,
                        // this seems sloppy. Being in the high charge tail,
                        // however, this should not have a big effect (<0.5ns).
                        // Still, would be nice to make a better decision here.
                        // FIXME when real data comes in, no more 'should'
                        // The correction in time to be applied to data-point:
                        twalk=thist;
                        // Flag this as a failed high-charge hit
                        fStatus.Set(kPCABitOffset + 5);
                      }
                    continue;
                  }
              }
            qstep=nextq-thisq;
          }
      } // else fQHS > thisq

    // use the extracted twalk to correct the time
    // (twalk contains the electronic delay and the time walk
    // from PCA calibration data). It does not contain the transit
    // time for this event.
    fTime=fTime-twalk;
    return;
  }

  void PCACal::ValidateTimeWalk(int lcn){
    // First, if we are using old SNO calibration constants,
    // set the flag and exit, no further validation is done.
    if( fIsSNO )
      {
        fStatus.Set(kPCABitOffset + 31);
        // The SNO status word is as follows (see thesis J. Cameron):
        // (100 x tube_type) + (10 x tube_status) + fit_quality
        // The tube type is 1 for a normal tube and > 1 for all other tube types.
        // Tube status:
        //      0 passed
        //      1 offline
        //      2 zero occupancy
        //      3 low occupancy
        //      4 low ECA
        // Fit Quality
        //      0 not fitted
        //      1 passed
        //      2 no points
        //      3 failed cut
        //      4 failed RMS
        //
        // So '101' is a perfect score in SNO-lingo
        // Anything else should be indicated in the statusword.
        if( fTimeWalkStatus[lcn] != 101 )
          {
            fStatus.Set(kPCABitOffset + 30);
          }
      }
    else
      // Only set the status word bits if we are using SNO+
      // calibration constants. If one uses SNO calibration
      // constants, only one bit will be set: 31.
      {
        // Check the PCA TimeWalk status words and flag accordingly
        // 1. Check RMS
        // RMS way too high
        if( BitManip::TestBit(fTimeWalkStatus[lcn],7) )
          {
            fStatus.Set(kPCABitOffset + 8);
            return;
          }
        // RMS too high
        if( BitManip::TestBit(fTimeWalkStatus[lcn],8) )
          {
            fStatus.Set(kPCABitOffset + 9);
          }
        // 2. High Q tail not fitted
        if( BitManip::TestBit(fTimeWalkStatus[lcn],9) )
          {
            fStatus.Set(kPCABitOffset + 5);
          }
        // 3. Check TStep
        if( BitManip::TestBit(fTimeWalkStatus[lcn],10) )
          {
            fStatus.Set(kPCABitOffset + 10);
          }
        // 4. Check FLate
        if( BitManip::TestBit(fTimeWalkStatus[lcn],11) )
          {
            fStatus.Set(kPCABitOffset + 11);
          }
        // 5. Check FOut
        if( BitManip::TestBit(fTimeWalkStatus[lcn],12) )
          {
            fStatus.Set(kPCABitOffset + 12);
          }
        // 6. Set general TimeWalk test bit
        if( fStatus.GetBits(kPCABitOffset + 8,kPCABitOffset + 12) != 0){
          fStatus.Set(kPCABitOffset + 13);
        }
      }
    return;
  }

  void PCACal::ValidateCalibrate(int lcn){
    // Check the PCA GainFit status words and flag
    if(fDoTimeWalk)
      {
        PCACal::ValidateTimeWalk(lcn);
      }
    // Set general FAIL bit
    if( fStatus.GetBits(kPCABitOffset + 2,kPCABitOffset + 4) != 0){
      fStatus.Set(kPCABitOffset + 0);
    }
    // Set general WARN bit
    if( fStatus.GetBits(kPCABitOffset + 5,kPCABitOffset + 7) != 0){
      fStatus.Set(kPCABitOffset + 1);
    }
    return;
  }


  void PCACal::CalStatus(){
    // This function looks at the tolerance level that was set by the user in
    // PMTCalibration.ratdb and sets the Q/T values to DS::INVALID accordingly.
    if( fDoTimeWalk )
      {
        if( fStatus.Get( kPCABitOffset + 30 ) )
          {
            // SNO calibration failed
            fTime = kFailValue;
          }

        if( fTolLevel==1 )
          {
            if( fStatus.Get( kPCABitOffset + 1 ) )
              {
                // WARN level
                fTime = kFailValue;
              }
          }
        if( fTolLevel==2 )
          {
            if( fStatus.Get( kPCABitOffset + 13 ) )
              {
                // TEST level
                fTime = kFailValue;
              }
          }
        if( fTolLevel==3 )
          {
            if( (fStatus.Get( kPCABitOffset + 13 )) || (fStatus.Get( kPCABitOffset+ 1 )) )
              {
                // TEST or WARN level
                fTime = kFailValue;
                return;
              }
          }
      }
    return;
  }

  // Handle non-fatal errors
  void PCACal::DontDie(LiveFlagErr flagErr, int info1, int info2)
  {
    switch(flagErr)
      {
      case NOPCA:
        // PCA selected but no PCA calibs requested
        warn <<newline;
        warn << YELLOW << "PCA calibration SETBACK:" << CLR << newline;
        warn << "PCACal::DontDie"<<newline;
        warn << "You asked me to do PCA, but have not selected any PCA calibrations..."<<newline;
        warn << "The PCA flag is set to "<<info1<<newline;
        warn <<newline;
        break;
      case GAINFIT:
        // DoGainFit selected, but this is not yet implemented
        detail <<newline;
        detail << YELLOW << "PCA calibration SETBACK:" << CLR << newline;
        detail << "PCACal::DontDie"<<newline;
        detail << "You asked me to do PCA Gain Fit..."<<newline;
        detail << "QHS, QHL and QLX will just be copied without changes. "<<newline;
        detail <<newline;
        break;
      case FIRSTINT:
        // The first interpolation point is not set
        debug <<newline;
        debug << "PCACal::DontDie"<<newline;
        debug << "PCACal - PMT "<<info1<<" : The first interpolation point is unavailable."<<newline;
        debug <<newline;
        break;
      case SECONDINT:
        // Only 1 interpolation point available
        debug <<newline;
        debug << "PCACal::DontDie"<<newline;
        debug <<"PCACal - PMT "<<info1<<" : There is only one available interpolation point."<<newline;
        debug <<newline;
        break;
      case LOWQHS:
        // fQHS is too low, it is probably noise
        debug <<newline;
        debug <<"PCACal - PMT "<<info1<<" : fQHS is too low, it is probably noise. QHS is "<<info2<<newline;
        debug <<newline;
        break;
      case SNOCONST:
        // PCA selected but with SNO constants
        detail <<newline;
        detail << YELLOW << "PCA calibration WARNING:" << CLR << newline;
        detail << "PCACal::DontDie"<<newline;
        detail << "Using SNO PCA constants."<<newline;
        detail <<newline;
        break;
      } //switch(flagErr)
  }

} // namespace RAT