// -*- C++/l -*-
//
// Description: 
//
// Author:      Nick Hastings <hastings@triumf.ca>
// Created:     Tue May 11 15:06:13 JST 2010
//

#include <sstream>
#include <cmath>

#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TF1.h"

#include "ISlowControlDatabase.hxx" // Get KeepDarkNoise ODB setting


#include "IFGDChanNumber.hxx"
#include "IFGDDataQuality.hxx"

#include "IFgdRawDataManager.hxx"

#include "IFGDDarkNoiseAndPed.hxx"
#include "IOnlinePulseProcessor.hxx"
#include "IEdgePulseFinder.hxx"
#include "IAdc2PePredicted.hxx"

class IGausNorm
{
  public:
    IGausNorm(void)
    {
    }
    Double_t operator()(Double_t* x, Double_t* par )
    {
  // Area par[0]
  // Mean par[1]
  // Sigma par[2]
  const Double_t a = (x[0]-par[1])/par[2];
  return fOneOnRoot2pi*par[0]/par[2]*exp( -0.5*a*a );
    }
    static Double_t area( Double_t height, Double_t sigma )
    {
      return height*sigma/fOneOnRoot2pi;
    }
  private:
    static const Double_t fOneOnRoot2pi = 0.3989422804014326779; // 1/sqrt( 2 * pi )
};
IDNChan::IDNChan(unsigned int miniCrate, unsigned int feb, unsigned int febCh,
     unsigned int bins, float xmin, float xmax )
  : fString("gaus"),
    fChanNum( miniCrate, feb, febCh ),
    fFitStatus(-999),
    fPrediction(0.0)
{
  const std::string n = fChanNum.Name();
  fH = TH1D( n.c_str(), n.c_str(), bins, xmin, xmax );
  fH.GetXaxis()->SetTitle("ADC");
  fF = TF1( (n+fString).c_str(), IGausNorm(), 10.0, 110.0, 3);
  fF.SetParNames("Area", "Mean", "Sigma");
}
void IDNChan::Reset(void)
{
  fH.Reset("M");

  const unsigned  int n = 3;
  for ( unsigned int i = 0; i < n; i++ ) {
    fF.SetParError(i, 0.0 );
  }
  fF.SetParameter(0, 100.0);
  fF.SetParameter(1, 40.0);
  fF.SetParameter(2, 6.0);
  fF.SetNDF( -1*static_cast<int>(n) );
  fF.SetChisquare(0.0);
}

void IDNChan::Fit(void)
{
  if ( fH.GetEntries() < 100.0 )
    return;

  COMETLog("IDNChan::Fit(): fitting" << fH.GetName() );

  const unsigned int firstbin = 18;
  const unsigned int lastbin  = 120;
  unsigned int maxbin = 0;
  float maxval = -9999.9;

  for ( unsigned int i = firstbin; i <= lastbin; i++ ) {
    if ( maxval < fH.GetBinContent(i) ) {
      maxbin = i;
      maxval = fH.GetBinContent(i);
    }
  }
  if ( maxval < 1.0 )
    return;

  float x = fH.GetBinCenter(maxbin);

  COMETVerbose("x of max bin = " << x );
  const float dx = 12.0;

  fF.SetParameter(0, IGausNorm::area( maxval, fF.GetParameter(2) ) );
  fF.SetParameter(1, x);

  int fitstatus = 0;
  for ( unsigned int fitnum = 0; fitnum<=2 && fitstatus==0; fitnum++ ) {
    for ( unsigned int i = 0 ; i < 3 ; i++ ) {
      COMETVerbose( fF.GetParName(i) << " = " << fF.GetParameter(i) );
    }
    COMETVerbose( "xmin = " <<  x-dx );
    COMETVerbose( "xmax = " << x+dx );
    fF.SetRange( x-dx, x+dx );
    TFitResultPtr result_ptr = fH.Fit( &fF, "RGOFFL", "", x-dx, x+dx);
    fitstatus = int(result_ptr); // int() is actually an overloaded
               // operator of the TFitResultsPtr class
    x = fF.GetParameter(1);
  }

  fFitStatus = fitstatus;
  if ( fFitStatus == 0 )
    return;

  COMETLog( "Warning, " << fH.GetName() << " fit failed with status: " << fFitStatus );
  for ( unsigned int i = 0 ; i < 3 ; i++ ) {
    COMETLog( fF.GetParName(i) << " = " << fF.GetParameter(i) );
  }
}

std::pair<double,double> IDNChan::GetFitMean(void) const
{
  return std::pair<double,double>(fF.GetParameter(1),
          fF.GetParError(1));
}
std::pair<double,double> IDNChan::GetFitEntries(void) const
{
  return std::pair<double,double>(fabs(fF.GetParameter(0)),
          fF.GetParError(0));
}
std::pair<double,double> IDNChan::GetMean(void) const
{
  double entries = fH.GetEntries();
  if ( entries < 1 )
    return std::pair<double,double>(0.0,10000.0);

  return std::pair<double,double>(fH.GetMean(),
          fH.GetRMS()/sqrt(entries) );
}
int IDNChan::GetFitStatus(void) const
{
  return fFitStatus;
}

///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
IAdcToPECalibSOff::IAdcToPECalibSOff()
{
  COMETVerbose("IAdcToPECalibSOff::InitChans(void)");
  for ( unsigned i = 0; i < IFGDChanNumber::NMiniCrate(); i++ ) {
    for ( unsigned j = 0; j < IFGDChanNumber::NFeb(); j++ ) {
      for ( unsigned k = 0; k < IFGDChanNumber::NFebCh(); k++ ) {
        fDNChans.push_back(new IDNChan(i,j,k));
      }
    }
  }
}

IAdcToPECalibSOff::~IAdcToPECalibSOff()
{
  for ( std::vector<IDNChan*>::iterator it = fDNChans.begin();
        it != fDNChans.end(); it++ ) {
    IDNChan* c = *it;
    delete c;
    c = 0;
  }
  fDNChans.clear();
}

void IAdcToPECalibSOff::Fit(void)
{
  COMETVerbose("fDNChans.size() = " << fDNChans.size() );
  for ( std::vector<IDNChan*>::iterator it = fDNChans.begin();
        it != fDNChans.end(); it++ ) {
    IDNChan* ch = *it;
    if ( !ch ) {
      COMETError("IDNChan* is a null pointer!");
      continue;
    }
    ch->Fit();
  }
}

int IAdcToPECalibSOff::processEvent(IFgdRawEvent* event)
{
  // Only analyse Dark Noise / Pedestal events or FGDLED trigger events
  // const UInt_t tb = event->GetTriggerBits();
  // if ( ! (tb & COMET::Trigger::kPedestal) &&
        //      ! (tb & COMET::Trigger::kFGDLED) )
  // 	return 0;

  // COMETVerbose("IAdcToPECalibSOff::processEvent(): Got pedestal trigger");

  // Some variables used for excluding pulses soon after large pulses on
  // the same waveform... fugly
  unsigned int index_prev = 999999;
  const int recover_time = 200;
  int t_prev = - recover_time;
  const float pulse_height_large = 500.0;
  float pulse_height_prev = 0.0;
    
  for( unsigned int i=0; i < event->getNumberOfPulses(); i++ ){
    IFgdPulse* p = event->getPulse(i);
    IFGDChanNumber cn( p->GetMinicrate(),
           p->GetFEB(),
           p->GetFebCh() );
    unsigned int index = cn.Index();
    
    if ( index != index_prev )
      pulse_height_prev = 0.0;
    
    if ( index >= fDNChans.size() ) {
      COMETError( "IAdcToPECalibSOff index out of range");
      continue;
    }

    // Don't use pulses near earlier large pulses in same waveform
    // Algorithm releys on pulses from the same channel being
    // adjacent in the pulse list.
    if ( pulse_height_prev > pulse_height_large && 
         (p->getStartBin()-t_prev)<recover_time ) {
      continue;
    }
    //COMETVerbose("Got Pulse");
    IDNChan* c = fDNChans[index];
    const float pulse_height =  p->getPulseHeightLo();

    // Only consider pulses after bin 20, since before this is the
    // region used to define the baseline and as such pulses here
    // are rubbish (NB. bin "20" is my guess)
    if ( p->getStartBin() > 20 ) {
      c->GetHist().Fill( pulse_height );
    }

    // Update the fugly variables
    t_prev = p->getStartBin();
    pulse_height_prev = pulse_height;
    index_prev = index;
  }
  return 0;
}
void IAdcToPECalibSOff::Reset(void)
{
  for ( std::vector<IDNChan*>::iterator it = fDNChans.begin();
        it != fDNChans.end(); it++ ) {
    (**it).Reset();
  }
}

///////////////////////////////////////////////////////////////
IFGDDarkNoiseAndPed::IFGDDarkNoiseAndPed( std::string fileName )
  : fManager ( new IFgdRawDataManager    ()    ),
    fTaskOPP ( new IOnlinePulseProcessor (4,4) ),
          fTaskEPF ( new IEdgePulseFinder      (4,4) ),
     fTaskAdcToPECalib ( new IAdcToPECalibSOff() ),
    fAdcToPEPredicted ( new IAdc2PePredicted() ),
    fFileName( fileName ),
    fRun(-1),
    fSubRun(-1),
    fTfirst(-1),
    fTlast(-1),
    fNChans    ( IFGDChanNumber::NIndex() ),
    fMean      ( new Double_t[fNChans] ),
    fDMean     ( new Double_t[fNChans] ),
    fFitMean   ( new Double_t[fNChans] ),
    fFitDMean  ( new Double_t[fNChans] ),
    fFitEntries( new Double_t[fNChans] ),
    fPredict   ( new Double_t[fNChans] ),
    fChi2ndf   ( new Double_t[fNChans] ),
    fFitStatus ( new Int_t[fNChans] ),
    fEntries   ( new Double_t[fNChans] ),
    fMiniCrate ( new Int_t[fNChans] ),
    fFeb       ( new Int_t[fNChans] ),
    fFebCh     ( new Int_t[fNChans] ),
    fIndex     ( new Int_t[fNChans] )
{
  std::vector<IDNChan*>& chans = fTaskAdcToPECalib->GetDNChans();

  if ( chans.size() != fNChans ) {
    COMETError( "IFGDDarkNoiseAndPed: number of channel missmatch.");
  }
}

IFGDDarkNoiseAndPed::~IFGDDarkNoiseAndPed()
{
  delete fManager;
  delete fTaskOPP;
  delete fTaskEPF;
  delete fTaskAdcToPECalib;
  delete fAdcToPEPredicted;

  delete [] fMean      ;
  delete [] fDMean     ;
  delete [] fFitMean   ;
  delete [] fFitDMean  ;
  delete [] fFitEntries;
  delete [] fPredict   ;
  delete [] fChi2ndf   ;
  delete [] fFitStatus ;
  delete [] fEntries   ;
  delete [] fMiniCrate ;
  delete [] fFeb       ;
  delete [] fFebCh     ;
  delete [] fIndex     ;
}

int IFGDDarkNoiseAndPed::Process( COMET::ICOMETEvent& event )
{
  // Only look at  pedestal (dark noise) triggers or FGD LED triggers
  // if (!event.CheckTrigger( COMET::Trigger::kPedestal ) &&
        //     !event.CheckTrigger( COMET::Trigger::kFGDLED ) ) 
  // 	return 0;

  // COMETVerbose("IFGDDarkNoiseAndPed::Process(): Got ped trigger");
  

  const COMET::ICOMETContext& c = event.GetContext();
  fAdcToPEPredicted->LoadEventInfo( c );
  const int run    = c.GetRun();
  const int subRun = c.GetSubRun();
  const int time   = c.GetTimeStamp();

  // Start of subrun that is not the first subrun
  if ( fRun != -1 && (run != fRun || subRun != fSubRun) ) {
    // Fill the tree using saved results from previous event
    EndPeriod();
  }
  // Start of any run
  if ( fRun == -1 || run != fRun || subRun != fSubRun ) {
    fTfirst = time;
    StartPeriod();
  }
  fRun = run;
  fSubRun = subRun;
  fTlast  = time;
  fNTriggers++;
  
  // Analyse	
  COMET::IHandle<COMET::IDigitContainer> fgdDigits = event.GetDigits("fgd");
  if ( !fgdDigits ) {
    COMETVerbose( "IFGDDarkNoiseAndPed: !event.GetDigits(\"fgd\")" );
    return 0;
  }
  
  if ( fManager->processEvent( fgdDigits ) != 0 ) {
    COMETVerbose( "IFGDDarkNoiseAndPed: fManager->processEvent(fgdDigits)!=0" );
  }

  std::vector<IDNChan*>& dnChans = fTaskAdcToPECalib->GetDNChans();
  for ( std::vector<IDNChan*>::iterator it = dnChans.begin();
        it != dnChans.end() ; it++ ) {
    IDNChan& ch = **it;
    const IFGDChanNumber& cn = ch.GetChanNum();
    ch.SetPrediction(fAdcToPEPredicted->
        CalculateSingleAvalancheHeight
        ( cn.MiniCrate(),
          cn.Feb(),
          cn.FebCh() )
        );
  }
  return 0;
}
int IFGDDarkNoiseAndPed::Initialize( void )
{
  COMETVerbose("IFGDDarkNoiseAndPed::Initialize( void )");
  fManager->addTask(fTaskOPP);
  fManager->addTask(fTaskEPF);
  fManager->addTask(fTaskAdcToPECalib);
  fManager->init();

  fFile = new TFile( fFileName.c_str(), "RECREATE" );
  fTree = new TTree( "dn", "FGD Dark Noise" );
  fTree->Branch( "run",    &fRun,    "run/I" );
  fTree->Branch( "subrun", &fSubRun, "subrun/I" );
  fTree->Branch( "t1", &fTfirst, "t1/I" );
  fTree->Branch( "t2", &fTlast,  "t2/I" );
  fTree->Branch( "nch", &fNChans, "nch/I");
  fTree->Branch( "ntrig", &fNTriggers, "ntrig/I");
  fTree->Branch( "keepdarknoise", &fKeepDarkNoise, "keepdarknoise/I");
  
  
  fTree->Branch( "mean",      fMean, "mean[nch]/D");
  fTree->Branch( "dmean",     fDMean, "dmean[nch]/D");
  fTree->Branch( "fitmean",   fFitMean, "fitmean[nch]/D");
  fTree->Branch( "fitdmean",  fFitDMean, "fitdmean[nch]/D");
  fTree->Branch( "fitentries",   fFitEntries, "fitentries[nch]/D");
  fTree->Branch( "predict",   fPredict, "predict[nch]/D");
  fTree->Branch( "chi2ndf",   fChi2ndf, "chi2ndf[nch]/D");
  fTree->Branch( "fitstatus", fFitStatus, "fitstatus[nch]/I");
  fTree->Branch( "entries",   fEntries, "entries[nch]/D");
  
  fTree->Branch( "minicrate", fMiniCrate, "minicrate[nch]/I");
  fTree->Branch( "feb",       fFeb,   "feb[nch]/I");
  fTree->Branch( "febch",     fFebCh, "febch[nch]/I");
  fTree->Branch( "index",     fIndex, "index[nch]/I");
  return 0;
}
int IFGDDarkNoiseAndPed::Finalize(void)
{
  COMETVerbose("IFGDDarkNoiseAndPed::Finalize( void )");
  fManager->end();
  EndPeriod();

  COMETVerbose( "Writing fTree to fFile");
  fFile->WriteTObject(fTree);
  fFile->Close();
  return 0;
}
void IFGDDarkNoiseAndPed::SetFileName( const std::string& fileName )
{
  COMETVerbose("IFGDDarkNoiseAndPed::SetFileName(" << fileName << ")" );
  fFileName = fileName;
}
bool IFGDDarkNoiseAndPed::SetOption( std::string option, std::string value )
{
  if ( option == "fgd-dn-root-file" ) {
    SetFileName( value );
    return true;
  }
  return false;
}
void IFGDDarkNoiseAndPed::StartPeriod(void)
{
  COMETLog("IFGDDarkNoiseAndPed::StartPeriod(void) " );
  fNTriggers = 0;
  for ( unsigned int i = 0; i < fNChans; i++ ) {
    fMean[i] = -9999.0;
    fDMean[i] = -9999.0;
    fFitMean[i] = -9999.0;
    fFitDMean[i] = -9999.0;
    fFitEntries[i] = -9999.0;
    fPredict[i] = -9999.0;
    fChi2ndf[i] = -9999.0;
    fFitStatus[i] = -9999;
    fEntries[i] = -9999.0;
    fMiniCrate[i] = -9999;
    fFeb[i]       = -9999;
    fFebCh[i]     = -9999;
    fIndex[i]     = -9999;
  }
  fTaskAdcToPECalib->Reset();
}
void IFGDDarkNoiseAndPed::EndPeriod(void)
{
  COMETLog("IFGDDarkNoiseAndPed::EndPeriod(void) " << fRun
     << " " << fSubRun << " " << fTfirst << " " << fTlast );
  fTaskAdcToPECalib->Fit();

  const std::vector<IDNChan*>& dnchans = fTaskAdcToPECalib->GetDNChans();
  const unsigned int n = dnchans.size();

  COMETVerbose( "n = " << n );
  COMETVerbose( "fNChans = " << fNChans );
  if (  n != fNChans ) {
    COMETError( "Wrong number of chans: " << n << "!="
          << fNChans );
    return;
  }

  // Query the ODB to find out what the fraction of the dark noise
  // triggers were kept
  ISlowControlMySQLRow* rowCMBSettings =
    ISlowControlDatabase::ODB().GetTable( fTlast, "FGD_CMB_Settings" );
  
  if ( rowCMBSettings ) {
    fKeepDarkNoise = rowCMBSettings->GetValueInteger("KeepDarkNoise");
  }
  else {
    COMETError( "Error: Unable to get KeepDarkNoise value from ODB" );
    fKeepDarkNoise = -1;
  }
  
  for ( unsigned int i = 0; i < n; i++ )
  {
    const IDNChan& ch = *(dnchans[i]);
    const unsigned int index = ch.GetChanNum().Index();

    if ( i != index ) {
      COMETError( "Index missmach : " << i << "!=" << index );
      return;
    }
    const TH1D& h = ch.GetHist();
    const TF1&  f = ch.GetFunc();
    std::pair<double,double> mdm     = ch.GetMean();
    std::pair<double,double> fit_mdm = ch.GetFitMean();
    std::pair<double,double> fit_entries = ch.GetFitEntries();
      
    fMean[i] = mdm.first;
    fDMean[i] = mdm.second;
    fFitMean[i] = fit_mdm.first;
    fFitDMean[i] = fit_mdm.second;
    fPredict[i] = ch.GetPrediction();
    fFitEntries[i] = fit_entries.first;
    const int ndf = f.GetNDF();
    fChi2ndf[i] = (ndf!=0) ? f.GetChisquare()/static_cast<float>(ndf)
      : 999999.9;
    fFitStatus[i] = ch.GetFitStatus();
    fEntries[i]   = h.GetEntries();
    const IFGDChanNumber& cn = ch.GetChanNum();
    fMiniCrate[i] = cn.MiniCrate();
    fFeb[i]       = cn.Feb();
    fFebCh[i]     = cn.FebCh();
    fIndex[i]     = i;
  }
  fTree->Fill();

  // Make a directory for the histograms
  std::ostringstream oss;
  oss << std::setw(8) << std::setfill('0') << fRun << "_"
      << std::setw(4) << std::setfill('0') << fSubRun;
  const std::string dname = oss.str();
  TDirectory* d = fFile->mkdir( dname.c_str() );
  // Write the histograms to the newly created directory
  COMETLog("IFGDDarkNoiseAndPed: writing histograms to " << dname );	
  for ( std::vector<IDNChan*>::const_iterator it = dnchans.begin();
        it != dnchans.end(); it++ ) {
    const IDNChan& c  = **it;
    d->WriteTObject( &c.GetHist() );
  }
}

// Local Variables:
// c-basic-offset: 8
// c-file-offsets: ((access-label . -2)(inclass . 4)(innamespace . 2))
// End: