/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Markus Gaug  09/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationCalc
//
//   Task to calculate the calibration conversion factors from the FADC
//   time slices. The integrated time slices have to be delivered by an 
//   MExtractedSignalCam. The pedestals by an MPedestalCam.
//
//   The output container MCalibrationCam holds one entry of type MCalibrationPix 
//   for every pixel. It is filled in the following way:
//
//   ProProcess: Search for MPedestalCam, MExtractedSignalCam
//               Initialize MCalibrationCam
//               Initialize pulser light wavelength
//               
//   ReInit:     MCalibrationCam::InitSize(NumPixels) is called which allocates
//               memory in a TClonesArray of type MCalibrationPix
//               Initialize number of used FADC slices
//               Optionally exclude pixels from calibration               
//
//   Process:    Every MCalibrationPix holds a histogram class,
//               MHCalibrationPixel which itself hold histograms of type:
//               HCharge(npix) (distribution of summed FADC time slice
//               entries)
//               HTime(npix) (distribution of position of maximum)
//               HChargevsN(npix) (distribution of charges vs. event number.
//
//  PostProcess:  All histograms HCharge(npix) are fitted to a Gaussian
//                All histograms HTime(npix) are fitted to a Gaussian
//                The histogram HBlindPixelCharge (blind pixel) is fitted to
//                a single PhE fit
//
//                The histograms of the PIN Diode are fitted to Gaussians
//
//                Fits can be excluded via the commands:
//                MalibrationCam::SkipBlindPixelFits()  (skip all blind
//                pixel fits)
//
//                Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
//                can be skipped with the command:
//                MalibrationCam::SkipHiLoGainCalibration()
//
//  Input Containers:
//   MRawEvtData
//   MPedestalCam
//   MExtractedSignalCam
//
//  Output Containers:
//   MCalibrationCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationCalc.h"

// FXIME: Usage of fstream is a preliminary workaround!
#include <fstream>

// FXIME: This has to be removed!!!! (YES, WHEN WE HAVE ACCESS TO THE DATABASE!!!!!)
#include "MCalibrationConfig.h"

#include <TSystem.h>
#include <TH1.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

#include "MGeomCam.h"
#include "MRawRunHeader.h"
#include "MRawEvtPixelIter.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include "MCalibrationCam.h"
#include "MCalibrationPix.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

#include "MCalibrationBlindPix.h"

ClassImp(MCalibrationCalc);

using namespace std;

const UInt_t MCalibrationCalc::fgBlindPixelId    = 559;
const UInt_t MCalibrationCalc::fgPINDiodeId      = 9999;
const Byte_t MCalibrationCalc::fgSaturationLimit = 254;
const Int_t  MCalibrationCalc::fgBlindPixelFirst = 3;
const Int_t  MCalibrationCalc::fgBlindPixelLast  = 14;
const Int_t  MCalibrationCalc::fgBlindPixelSinglePheCut = 400;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
    : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
      fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL)
{

    fName  = name  ? name  : "MCalibrationCalc";
    fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";

    AddToBranchList("MRawEvtData.fHiGainPixId");
    AddToBranchList("MRawEvtData.fLoGainPixId");
    AddToBranchList("MRawEvtData.fHiGainFadcSamples");
    AddToBranchList("MRawEvtData.fLoGainFadcSamples");

    Clear();
    SetBlindPixelRange();
}

void MCalibrationCalc::Clear(const Option_t *o)
{
  
    SETBIT(fFlags, kUseBlindPixelFit);
    SETBIT(fFlags, kUseQualityChecks);
    SETBIT(fFlags, kHiLoGainCalibration);

    CLRBIT(fFlags, kHiGainOverFlow);
    CLRBIT(fFlags, kLoGainOverFlow);

    fBlindPixelFirst   = 0;
    fBlindPixelLast    = 0;

    fNumBlindPixelSinglePhe = 0;
    fNumBlindPixelPedestal  = 0;  

    fNumHiGainSamples  = 0;
    fNumLoGainSamples  = 0;
    fConversionHiLo    = 0;
    fNumExcludedPixels = 0;

    fColor = kECT1;
}

void MCalibrationCalc::SetBlindPixelRange(Int_t first, Int_t last)
{
  
  fBlindPixelFirst = first;
  fBlindPixelLast  = last;

}


MCalibrationBlindPix *MCalibrationCalc::GetBlindPixel() const
{
    return fCalibrations->GetBlindPixel();
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MHCalibrationBlindPixel
//  - MCalibrationCam
//
// The following output containers are only searched, but not created
//
//  - MTime
//
Int_t MCalibrationCalc::PreProcess(MParList *pList)
{

    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRawEvt)
    {
      *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
    }

    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!runheader)
      *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    else
      if (runheader->GetRunType() == kRTMonteCarlo)
        {
          return kTRUE;
        }
    
    fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
    if (!fCalibrations)
      {
        *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;        
        return kFALSE;
      }

    fEvtTime      = (MTime*)pList->FindObject("MTime");

    switch (fColor)
      {
      case kEBlue:
        fCalibrations->SetColor(MCalibrationCam::kECBlue);
	break;        
      case kEGreen:
        fCalibrations->SetColor(MCalibrationCam::kECGreen);      
	break;
      case kEUV:
        fCalibrations->SetColor(MCalibrationCam::kECUV);            
	break;
      case kECT1:
        fCalibrations->SetColor(MCalibrationCam::kECCT1);            
	break;
      default:
        fCalibrations->SetColor(MCalibrationCam::kECCT1); 
      }

    fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPedestals)
      {
        *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
        return kFALSE;
      }


    fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
    if (!fSignals)
      {
        *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
        return kFALSE;
      }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//
Bool_t MCalibrationCalc::ReInit(MParList *pList )
{
 
    fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRunHeader)
    {
      *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
      return kFALSE;
    }


    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
      *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
      return kFALSE;
    }


    fCalibrations->SetGeomCam(cam);

    fNumHiGainSamples =  fSignals->GetNumUsedHiGainFADCSlices();
    fNumLoGainSamples =  fSignals->GetNumUsedLoGainFADCSlices();
    fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);

    UInt_t npixels = cam->GetNumPixels();

    for (UInt_t i=0; i<npixels; i++)
      {
        
        MCalibrationPix &pix = (*fCalibrations)[i];
        pix.DefinePixId(i);

        pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
                                    fSignals->GetLastUsedSliceHiGain());
        pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
                                    fSignals->GetLastUsedSliceLoGain());
        
        if (!TESTBIT(fFlags,kUseQualityChecks))
          pix.SetExcludeQualityCheck();

        // Exclude the blind pixel and the PIN Diode from normal pixel calibration:
        if (i == fgBlindPixelId)
          pix.SetExcluded();

        if (i == fgPINDiodeId)
          pix.SetExcluded();

     }
    
    //
    // Look for file to exclude pixels from analysis
    //
    if (!fExcludedPixelsFile.IsNull())
      {
        
        fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
        
        //
        // Initialize reading the file
        //
        ifstream in(fExcludedPixelsFile.Data(),ios::in);

        if (in)
          {
            *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
            //
            // Read the file and count the number of entries
            //
            UInt_t pixel = 0;
            
            while (++fNumExcludedPixels)
              {
                
                in >> pixel;

                if (!in.good())
                  break;
                //
                // Check for out of range
                //
                if (pixel > npixels)
                  {
                    *fLog << warn << "WARNING: To be excluded pixel: " << pixel 
                          << " is out of range " << endl;
                    continue;
                  }
                //
                // Exclude pixel
                //
                MCalibrationPix &pix = (*fCalibrations)[pixel];
                pix.SetExcluded();
                
                *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
                
              }
            
            if (--fNumExcludedPixels == 0)
              *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data() 
                    << "'" << " is empty " << endl;
            else
              fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
            
          }
        else
          *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
      }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MCerPhotEvt container.
//
Int_t MCalibrationCalc::Process()
{

  //
  // Initialize pointers to blind pixel, PIN Diode and individual pixels
  //
  MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
  
  MRawEvtPixelIter pixel(fRawEvt);
  
  //
  // Create a loop to fill the calibration histograms
  // Search for: a signal in MExtractedSignalCam 
  // Fill histograms with:
  //             charge
  //             charge vs. event nr.
  //
  while (pixel.Next())
    {
      
      const UInt_t pixid = pixel.GetPixelId();
      
      MCalibrationPix &pix = (*fCalibrations)[pixid];
      
      MExtractedSignalPix &sig =  (*fSignals)[pixid];
      
      const Float_t sumhi  = sig.GetExtractedSignalHiGain();
      const Float_t sumlo  = sig.GetExtractedSignalLoGain();

      Float_t abstime = 0.;
#if 0
      Float_t reltime = 0.;

      if (TESTBIT(fFlags,kUseTimes))
        {

          //
          // Have a look in MArrivalTime, 
          // otherwise search the position of maximum bin 
          // in MRawEvtData
          //
          if (fArrivalTime)
            {
              abstime = (*fArrivalTime)[pixid];
              reltime = abstime - (*fArrivalTime)[1];
            }
          
          else
            {
              if (pixid == 1)
                referencetime = (Float_t)pixel.GetIdxMaxHiGainSample();
              if (sig.IsLoGainUsed())
                {
                  abstime = (Float_t)pixel.GetIdxMaxLoGainSample();
                  //            reltime = abstime - referencetime;
                }
              else
                {
                  abstime = (Float_t)pixel.GetIdxMaxHiGainSample();
                  //                  reltime = abstime - referencetime;
                }
              //            }
              //        }  /* if Use Times */

#endif
      switch(pixid)
        {
          
        case fgBlindPixelId:

          if (TESTBIT(fFlags,kUseBlindPixelFit))
            {
          
              Byte_t *ptr = pixel.GetHiGainSamples();

              Int_t blindpixelsumhi = 0;
              Int_t blindpixelsumlo = 0;
              //
              // We need a dedicated signal extractor for the blind pixel
              //
              Int_t diff  = 0;
              Int_t last  = fBlindPixelLast;
              Int_t first = fBlindPixelFirst;

              if (last > 15)
                {
                  diff = last - 15;
                  last = 15;
                }

              Byte_t *start   = ptr + first - 1;
              Byte_t *end     = ptr + last  - 1;

              ptr = start;

              Int_t sum = 0;
              
              while (ptr<=end)
                {
                  sum += *ptr;
                  ptr++;
                }

              if (diff > 0)
                {
                  ptr = pixel.GetLoGainSamples();
                  end = ptr + diff - 1;
                  
                  while (ptr<=end)
                    {
                      sum += *ptr;
                      ptr++;
                    }
                }
              
              blindpixelsumhi = sum;

              ptr = pixel.GetLoGainSamples();

              start   = ptr + fBlindPixelFirst - 1;
              end     = ptr + fBlindPixelLast;

              ptr = start;

              sum = 0;
              
              while (++ptr<end)
                sum += *ptr;

              blindpixelsumlo = sum;

              //              if (!CalcSignalBlindPixel(hiptr, blindpixelsumhi))
              //                return kFALSE;

              if (!blindpixel.FillCharge(blindpixelsumhi)) 
                *fLog << warn << 
                  "Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsumhi << endl;
              
              //              Byte_t *loptr = pixel.GetLoGainSamples();
              //              CalcSignalBlindPixel(loptr, blindpixelsumlo);
              
              blindpixel.FillGraphs(blindpixelsumhi,blindpixelsumlo);
          
              TH1I *hist;

              if (blindpixelsumhi > fBlindPixelSinglePheCut)
                {
                  hist = (blindpixel.GetHist())->GetHSinglePheFADCSlices();
                  fNumBlindPixelSinglePhe++;
                }
              else
                {
                  hist = (blindpixel.GetHist())->GetHPedestalFADCSlices();
                  fNumBlindPixelPedestal++;
                }
              
              ptr = pixel.GetHiGainSamples();
              for (Int_t i=1;i<16;i++)
                hist->Fill(i,*ptr++);

              ptr = pixel.GetLoGainSamples();
              for (Int_t i=16;i<31;i++)
                hist->Fill(i,*ptr++);

            } /* if use blind pixel */
          
          break;
        default:
              
          if (pix.IsExcluded())
            continue;

          pix.FillGraphs(sumhi,sumlo);

          if (sig.IsLoGainUsed())
            {
              
              if (!pix.FillChargeLoGain(sumlo))
                *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid 
                      << " signal = " << sumlo << endl;

              if (!pix.FillAbsTimeLoGain(abstime)) 
                *fLog << warn << "Could not fill Lo Gain Abs. Time of pixel: " 
                      << pixid << " time = " << abstime << endl;
              /*
                if (!pix.FillRelTimeLoGain(reltime)) 
                *fLog << warn << "Could not fill Lo Gain Rel. Time of pixel: " 
                    << pixid << " time = " << reltime << endl;
              */
            } /* if (sig.IsLoGainUsed()) */
          else
            {
              if (!pix.FillChargeHiGain(sumhi))
                *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid 
                      << " signal = " << sumhi << endl;
              
              if (!pix.FillAbsTimeHiGain(abstime))
                *fLog << warn << "Could not fill Hi Gain Abs. Time of pixel: " 
                      << pixid << " time = " << abstime << endl;
              /*
              if (!pix.FillRelTimeHiGain(reltime))
                    *fLog << warn << "Could not fill Hi Gain Rel. Time of pixel: " 
                          << pixid << " time = " << reltime << endl;
              */
            } /* else (sig.IsLoGainUsed()) */
          break;
          
        } /* switch(pixid) */
      
    } /* while (pixel.Next()) */

  return kTRUE;
}

Int_t MCalibrationCalc::PostProcess()
{

  *fLog << inf << endl;

  *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;

  //
  // Cut edges to make fits and viewing of the hists easier  
  //
  fCalibrations->CutEdges();

  // 
  // Fit the blind pixel
  //
  if (TESTBIT(fFlags,kUseBlindPixelFit))
    {
      //
      // Get pointer to blind pixel
      //
      MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
      
      *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;

      //
      // retrieve mean and sigma of the blind pixel pedestal, 
      // so that we can use it for the fit
      //
      Float_t pedestal;
      Float_t pederr;
      Float_t pedsigma;
      Float_t pedsigmaerr;
      
      const ULong_t nentries    = fPedestals->GetTotalEntries();
      const Int_t   nslices     = fBlindPixelLast-fBlindPixelFirst+1;
      const Float_t sqrslice    = TMath::Sqrt((Float_t)nslices);

      MPedestalPix    &pedpix  = (*fPedestals)[fgBlindPixelId];
      pedestal    = pedpix.GetPedestal()*nslices;
      pederr      = pedpix.GetPedestalRms()*nslices/nentries;
      pedsigma    = pedpix.GetPedestalRms()*sqrslice;
      pedsigmaerr = pederr/2.;
      //
      // retrieve the histogram containers
      //
      MHCalibrationBlindPixel *hist = blindpixel.GetHist();
      
      hist->SetMeanPedestal(pedestal);
      hist->SetMeanPedestalErr(pederr);
      hist->SetSigmaPedestal(pedsigma);
      hist->SetSigmaPedestalErr(pedsigmaerr);
  
      if (!blindpixel.FitCharge())
        {
          *fLog << warn << "Could not fit the blind pixel! " << endl;
          *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
          fCalibrations->SetBlindPixelMethodValid(kFALSE);
        }
      else
        fCalibrations->SetBlindPixelMethodValid(kTRUE);
      
      if (blindpixel.CheckOscillations())
        fCalibrations->SetBlindPixelMethodValid(kFALSE);

      TH1I *sphehist = hist->GetHSinglePheFADCSlices();
      TH1I *pedhist  = hist->GetHPedestalFADCSlices();

      if (fNumBlindPixelSinglePhe > 1)
        sphehist->Scale(1./fNumBlindPixelSinglePhe);
      if (fNumBlindPixelPedestal > 1)
        pedhist->Scale(1./fNumBlindPixelPedestal);
      
      blindpixel.DrawClone();
    }
  else 
    *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;

  *fLog << err << "total: " << GetNumExecutions() << " sphe: " << fNumBlindPixelSinglePhe << " ped: " << fNumBlindPixelPedestal << endl;

  
  *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;

  //
  // loop over the pedestal events and check if we have calibration
  //
  for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
    {

      MCalibrationPix &pix = (*fCalibrations)[pixid];

      //
      // Check if the pixel has been excluded from the fits
      //
      if (pix.IsExcluded())
        continue;

      //
      // get the pedestals
      //
      const Float_t ped    = (*fPedestals)[pixid].GetPedestal();
      const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms();

      //
      // set them in the calibration camera
      //
      pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples);

      //
      // perform the Gauss fits to the charges
      //
      pix.FitCharge();

      //
      // check also for oscillations
      // 
      pix.CheckOscillations();
      
    }

  if (TESTBIT(fFlags,kUseBlindPixelFit) && fCalibrations->IsBlindPixelMethodValid())
    {
      if (!fCalibrations->CalcFluxInsidePlexiglass())
        {
          *fLog << err 
                << "Could not calculate the number of photons from the blind pixel " << endl;
          *fLog << err 
                << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
        fCalibrations->SetBlindPixelMethodValid(kFALSE);          
        }
    }
  else
    *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;


  if (fCalibrations->IsPINDiodeMethodValid())
    {
      if (!fCalibrations->CalcFluxOutsidePlexiglass())
        {
          *fLog << err 
                << "Could not calculate the number of photons from the PIN Diode " << endl;
          fCalibrations->SetPINDiodeMethodValid(kFALSE);
        }
    }
  else
    *fLog << inf << GetDescriptor() << ": Skipping PIN Diode Calibration! " << endl;

  fCalibrations->SetReadyToSave();
  
  return kTRUE;
}


Bool_t MCalibrationCalc::CalcSignalBlindPixel(Byte_t *ptr, Float_t &signal) const
{

  Byte_t *newptr  = ptr;
  Byte_t *start   = ptr + fBlindPixelFirst-1;
  Byte_t *end     = ptr + fBlindPixelLast;

  Byte_t sum = 0;
  Int_t  sat = 0;
  
  newptr = start;

  while (ptr<end)
    {
      sum += *ptr;
      if (*ptr++ >= fgSaturationLimit)
        sat++;
    }

  if (sat)
    {
      *fLog << err << "HI Gain Saturation occurred in the blind pixel! " 
            << " Do not know yet how to treat this ... aborting " << endl;
      *fLog << err << "If you need absolutely any other kind of calibration, " 
            << " use SkipBlindPixelFit() " << endl;
      return kFALSE;
    }

  signal = (Float_t)sum;

  return kTRUE;
}
