#include <RAT/FrontEndProc.hh>
#include <RAT/PMTPulse.hh>
#include <RAT/PMTWaveform.hh>
#include <RAT/DS/MCHit.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>
#include <RAT/BitManip.hh>
#include <RAT/DS/MCHit.hh>
#include <RAT/DS/MCPE.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>

#include <CLHEP/Random/RandGauss.h>
#include <G4ThreeVector.hh>

#include <cmath>
#include <algorithm>
#include <vector>

namespace RAT {

inline bool Cmp_PMTPulse_TimeAscending(const PMTPulse *a,
                                       const PMTPulse *b)
{
   double atime = a->GetPulseStartTime();
   double btime = b->GetPulseStartTime();
   return atime < btime;
}

FrontEndProc::FrontEndProc() : Processor("FrontEnd")
{

}

void
FrontEndProc::BeginOfRun(DS::Run&)
{
  fLdaq = DB::Get()->GetLink("DAQ");

  //sampling time in ns --- this is the size of a PMT time window
  fSamplingTimeDB = fLdaq->GetD("sampling_time");
  //short integration time in ns
  fShortIntTimeDB = fLdaq->GetD("shortint_time");
  //long integration time in ns
  fLongIntTimeDB = fLdaq->GetD("longint_time");
  //factor that attenuates low gain charges
  fLowGainFactorDB = fLdaq->GetD("lowgain_factor");
  //Pulse type: 0=square pulses, 1=real pulses
  fPulseTypeDB = fLdaq->GetI("pulse_type");
  //width of a PMT pulse in ns
  fPulseWidthDB = fLdaq->GetD("pulse_width");
  //offset of a PMT pulse in mV
  fPulseOffsetDB = fLdaq->GetD("pulse_offset");
  //mean of a PMT pulse in ns
  fPulseMeanDB = fLdaq->GetD("pulse_mean");
  //stepping time for discrimination
  fStepTimeDB = fLdaq->GetD("step_time");
  //Minimum pulse height to consider
  fPulseMinDB = fLdaq->GetD("pulse_min");
  //width of noise in adc counts
  fNoiseAmplDB = fLdaq->GetD("noise_amplitude");
  //time before discriminator fires that sampling gate opens
  fGDelayDB = fLdaq->GetD("gate_delay");

  detail << "FrontEndProc::BeginOfRun: DAQ constants loaded" << newline;
  detail << "  PMT Pulse type: " << (fPulseTypeDB==0 ? "square" : "realistic") << newline;
  detail << dformat("  PMT Pulse Width: ....................... %5.1f ns\n", fPulseWidthDB);
  detail << dformat("  PMT Pulse Offset: ...................... %5.1f ADC Counts\n", fPulseOffsetDB);
  detail << dformat("  PMT Pulse Mean: ........................ %5.1f\n", fPulseMeanDB);
  detail << dformat("  Min PMT Pulse Height: .................. %5.1f ns\n", fPulseMinDB);
  detail << dformat("  PMT Channel Short Integration Time: .... %6.2f ns\n", fShortIntTimeDB);
  detail << dformat("  PMT Channel Long Integration Time: ..... %6.2f ns\n", fLongIntTimeDB);
  detail << dformat("  PMT Channel Total Sample Time: ......... %6.2f ns\n", fSamplingTimeDB);
  detail << dformat("  PMT Channel Stepping Time: ............. %6.2f ns\n", fStepTimeDB);
  detail << dformat("  PMT Channel Low Gain Attenuation Factor: %6.2f\n", fLowGainFactorDB);
  detail << dformat("  Channel Gate Delay: .................... %5.1f ns\n", fGDelayDB);
  detail << dformat("  Hi Freq. Channel Noise: ................ %6.2f adc counts\n", fNoiseAmplDB);

  fFlasherCrossTalk.BeginOfRun();
}

  Processor::Result FrontEndProc::DSEvent(DS::Run& run, DS::Entry& ds)
{
  // Here we simulate the action of the front end electronics, whose job is to
  // simulate the charge spectrum and cluster hits on individual PMTs. We do
  // the clustering depending upon a sampling window, so that a given MCPMT can
  // have several samples corresponding to a single detected event. The output
  // is kept in the data structure in DS::MCHits and the samples in and array
  // within MCHits called MCSamples.

  if (!run.GetMCFlag()) {
    warn << "FrontEndProc::DSEvent : There is no MC branch for this Run. Most likely this is real data. Aborting..." << newline;
    warn << "FrontEndProc::DSEvent : Check the options in your macro, and your input file." << newline;
    return Processor::ABORT;
  }

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

  DS::MC& mc = ds.GetMC();

  //Outer loop is over all hit PMTs
  for (size_t imcpmt=0; imcpmt < mc.GetMCPMTCount(); imcpmt++) {
    DS::MCPMT &mcpmt = mc.GetMCPMT(imcpmt);

    //GDOG: Check whether we want to use Channel Hardwre Status (formerly DQXX) info
    //to determine tube status
    // If not, continue as you were
    // If so, check whether the tube is online and, if not, skip it
    int ipmt = mcpmt.GetID();
    bool recordhit = true;
    if( channelHardwareStatus.IsEnabled() )
      recordhit = channelHardwareStatus.IsDAQEnabled(ipmt);
    if(recordhit){

      //For each hit PMT, we need to build an array of PMT pulses based on the
      //photons that hit it, time sort them, and then sum them to see if we are over
      // the discriminator thresho.d


      PMTWaveform pmtwf;

      //Now loop over photons, and create a PMT pulse for each.
      double TimeNow;
      double PulseDuty=0.0;
      for (size_t ipe=0; ipe < mcpmt.GetMCPECount(); ipe++) {
          const DS::MCPE& mcphotoelectron = mcpmt.GetMCPE(ipe);
          TimeNow = mcphotoelectron.GetFrontEndTime();
          float charge =  mcphotoelectron.GetCharge();

          PMTPulse *pmtpulse;
          if (fPulseTypeDB==0){
            pmtpulse = new SquarePMTPulse; //square PMT pulses
          }
          else{
            pmtpulse = new RealPMTPulse; //real PMT pulses shape
          }
          pmtpulse->SetStepTime(fStepTimeDB);
          pmtpulse->SetPulseMin(fPulseMinDB);
          pmtpulse->SetPulseCharge(charge);
          pmtpulse->SetPulseWidth(fPulseWidthDB);
          pmtpulse->SetPulseOffset(fPulseOffsetDB);
          pmtpulse->SetPulseMean(fPulseMeanDB);
          pmtpulse->SetPulseStartTime(TimeNow); //also sets end time
          PulseDuty += pmtpulse->GetPulseEndTime()-pmtpulse->GetPulseStartTime();
          pmtwf.fPulse.push_back(pmtpulse);

      } //end loop over photons that have hit this PMT.

      //Sort pulses in time order
      std::sort(pmtwf.fPulse.begin(),pmtwf.fPulse.end(),Cmp_PMTPulse_TimeAscending);

      //Now we start to look for cases where the discriminator triggers
      //We create a new `sample' every time we do, which is the equivalent of
      //one PMT hit as far as electronics are concerned.  The sample width is the
      //length of SNO's integrator `RESET' pulse, which is also equivalent to the
      //long integration time
      while (pmtwf.fPulse.size()>0){

        int NextPulse=0;
        float qfuzz=0.0;

        TimeNow = pmtwf.fPulse[0]->GetPulseStartTime();
        double LastPulseTime =
          pmtwf.fPulse[pmtwf.fPulse.size()-1]->GetPulseEndTime();


        float NoiseAmpl = fNoiseAmplDB/sqrt(PulseDuty/fStepTimeDB);
        while( qfuzz < channelHardwareStatus.GetThresholdAboveNoise(ipmt) && TimeNow < LastPulseTime){
          float pmtsum = pmtwf.GetHeight(TimeNow); //gives height of pulses at TimeNow

          //Now add some hi freq noise to the sum, and check it against threshold
          qfuzz = pmtsum+NoiseAmpl*CLHEP::RandGauss::shoot();


          // move forward in time by step or if we're at baseline
          // move ahead to next pulse
          if (pmtsum==0.0){
            NextPulse = pmtwf.GetNext(TimeNow);
            if (NextPulse>0)
              {TimeNow = pmtwf.fPulse[NextPulse]->GetPulseStartTime();}
            else
              {TimeNow= LastPulseTime+1.;}
          }
          else{
            TimeNow += fStepTimeDB;
          }
        }

        if (TimeNow < LastPulseTime){

          //We have a hit above threshold
          DS::MCHit hit;

          hit.SetID(ipmt);
          //Now we begin integration, with information stored as an MC `sample'
          hit.SetTime(TimeNow);

          //We start integration gate earlier than discriminator firing
          double TimeStartGate = TimeNow - fGDelayDB;

          double sample_charge_qhs = 0.;
          double sample_charge_qhl = 0.;
          double sample_charge_qlx = 0.;

          //Set end times for different integration gates
          double TimeEndSIGate = TimeStartGate+fShortIntTimeDB;
          double TimeEndLIGate = TimeStartGate+fLongIntTimeDB;

          //Sampling gate may be longer than long integration gate-->enforces channel
          //deadtime.
          double TimeEndSampGate = TimeStartGate+fSamplingTimeDB;

          unsigned int ipulse=0;
          while (ipulse < pmtwf.fPulse.size() &&
                 pmtwf.fPulse[ipulse]->GetPulseStartTime()<TimeEndLIGate){

            double qhighshort=pmtwf.fPulse[ipulse]->Integrate(TimeStartGate,TimeEndSIGate);
            sample_charge_qhs+=qhighshort;
            sample_charge_qhl+=qhighshort+pmtwf.fPulse[ipulse]->Integrate(TimeEndSIGate,TimeEndLIGate);
            ipulse++;
          }

          // Now check for any pulses remaining before the end of the sampling window
          // But DO NOT add in the Q
          // we just want to skip them so they're not inc in the next sample
          while (ipulse < pmtwf.fPulse.size() &&
                 pmtwf.fPulse[ipulse]->GetPulseStartTime()<TimeEndSampGate ){ipulse++;}

          //now set all charges
          hit.SetQHS(sample_charge_qhs);
          hit.SetQHL(sample_charge_qhl);
          sample_charge_qlx = sample_charge_qhl/fLowGainFactorDB;
          hit.SetQLX(sample_charge_qlx);

          std::vector<PMTPulse*>::iterator it;
          for ( it=pmtwf.fPulse.begin() ; it < pmtwf.fPulse.begin()+ipulse; ++it )
            delete *it;

          pmtwf.fPulse.erase(pmtwf.fPulse.begin(),pmtwf.fPulse.begin()+ipulse);

          mc.GetUnbuiltMCHits().AddPMT( hit, pmtInfo.GetType( hit.GetID() ) );
        }
        else{ //we never got a pulse above discrminator threshold+
          std::vector<PMTPulse*>::iterator it;
          for ( it=pmtwf.fPulse.begin() ; it < pmtwf.fPulse.end(); ++it )
            delete *it;

          pmtwf.fPulse.erase(pmtwf.fPulse.begin(),pmtwf.fPulse.end());
        }
      } //end while loop over `pmtwf.fPulse.size()>0'

      //clean up just in case
      for (unsigned int i = 0; i<pmtwf.fPulse.size();i++){
        delete pmtwf.fPulse[i];
      }
    }// end test of tube status
  } //end loop over all hit PMTs

  // Generate flasher cross talk, adds hits to UnbuiltMCHits
  fFlasherCrossTalk.GenerateFlasherCrossTalk( run, ds);

  return OK;
}

} // namespace RAT