/* ======================================================================== *\
!
! *
! * 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   11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MCalibrationPix                                                         //
//                                                                         //
// This is the storage container to hold informations about the pedestal   //
// (offset) value of one Pixel (PMT).                                      //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationPix.h"
#include "MCalibrationConfig.h"

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

ClassImp(MCalibrationPix);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor: 
//
// The following values are initialized to meaningful values:
//
// - The Electronic Rms to 1.5 per FADC slice
// - The uncertainty about the Electronic RMS to 0.3 per slice
// - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17. 
//   We use here the Square of the Munich definition, thus:
//   Mean F-Factor  = 1.15*1.15 = 1.32
//   Error F-Factor = 2.*0.02   = 0.04
//
MCalibrationPix::MCalibrationPix(const char *name, const char *title)
    : fPixId(-1),
      fElectronicPedRms(1.5),
      fErrElectronicPedRms(0.3),
      fFactor(1.32),
      fFactorError(0.04),
      fChargeLimit(3.),
      fChargeErrLimit(0.),
      fChargeRelErrLimit(1.),
      fFlags(0)
{

  fName  = name  ? name  : "MCalibrationPixel";
  fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";

  //
  // At the moment, we don't have a database, yet, 
  // so we get it from the configuration file
  //
  fConversionHiLo      = gkConversionHiLo;
  fConversionHiLoError = gkConversionHiLoError;  

  fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");

  if (!fHist)
    *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;

  Clear();
}

MCalibrationPix::~MCalibrationPix() 
{
  delete fHist;
}


// ------------------------------------------------------------------------
//
// Invalidate values
//
void MCalibrationPix::Clear(Option_t *o)
{

  fHist->Reset();

  CLRBIT(fFlags, kHiGainSaturation);
  CLRBIT(fFlags, kExcluded);
  CLRBIT(fFlags, kFitValid);
  CLRBIT(fFlags, kFitted);
  CLRBIT(fFlags, kBlindPixelMethodValid);
  CLRBIT(fFlags, kFFactorMethodValid);
  CLRBIT(fFlags, kPINDiodeMethodValid);

  fCharge                           =  -1.;
  fErrCharge                        =  -1.;
  fSigmaCharge                      =  -1.;
  fErrSigmaCharge                   =  -1.;
  fRSigmaSquare                     =  -1.;
  fChargeProb                       =  -1.;
  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fErrPedRms                        =   0.;
  fTime                             =  -1.;
  fSigmaTime                        =  -1.;
  fTimeChiSquare                    =  -1.;
  fPheFFactorMethod                 =  -1.;
  fPheFFactorMethodError            =  -1.;
  fConversionFFactorMethod          =  -1.;
  fConversionBlindPixelMethod       =  -1.;
  fConversionPINDiodeMethod         =  -1.;
  fConversionErrorFFactorMethod     =  -1.;
  fConversionErrorBlindPixelMethod  =  -1.;
  fConversionErrorPINDiodeMethod    =  -1.;
  fConversionSigmaFFactorMethod     =  -1.;
  fConversionSigmaBlindPixelMethod  =  -1.;
  fConversionSigmaPINDiodeMethod    =  -1.;

}


void MCalibrationPix::DefinePixId(Int_t i)
{
  
  fPixId = i;
  fHist->ChangeHistId(i);
  
}


// --------------------------------------------------------------------------
//
// Set the pedestals from outside
//
void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms)
{

  fPed    = ped;    
  fPedRms = pedrms;
  
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
{
  fConversionFFactorMethod = c;
  fConversionErrorFFactorMethod = err;
  fConversionSigmaFFactorMethod = sig;
}


// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
{
  fConversionBlindPixelMethod = c;
  fConversionErrorBlindPixelMethod = err;
  fConversionSigmaBlindPixelMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
{
  fConversionPINDiodeMethod = c ;
  fConversionErrorPINDiodeMethod = err;
  fConversionSigmaPINDiodeMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the Hi Gain Saturation Bit from outside (only for MC)
//
void MCalibrationPix::SetHiGainSaturation(Bool_t b)
{

  if (b) 
    { 
      SETBIT(fFlags, kHiGainSaturation); 
      fHist->SetUseLoGain(1); 
    }
  else    
    {
      CLRBIT(fFlags, kHiGainSaturation); 
      fHist->SetUseLoGain(0); 
    }
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetExcluded(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded); 
}


// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetExcludeQualityCheck(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck); 
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetFitValid(Bool_t b )    
{ 
  b ?  SETBIT(fFlags, kFitValid) : CLRBIT(fFlags, kFitValid); 
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetFitted(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted); 
}
    
// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetFFactorMethodValid(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b )  
{ 
  b ?  SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); 
}


Bool_t MCalibrationPix::IsExcluded() const
 { 
   return TESTBIT(fFlags,kExcluded);  
 }

Bool_t MCalibrationPix::IsFitValid() const 
{
  return TESTBIT(fFlags, kFitValid);  
}

Bool_t MCalibrationPix::IsFitted() const 
{ 
  return TESTBIT(fFlags, kFitted);    
}

Bool_t MCalibrationPix::IsBlindPixelMethodValid() const 
{ 
  return TESTBIT(fFlags, kBlindPixelMethodValid); 
}

Bool_t MCalibrationPix::IsFFactorMethodValid()  const 
{ 
  return TESTBIT(fFlags, kFFactorMethodValid);     
}

Bool_t MCalibrationPix::IsPINDiodeMethodValid() const 
{ 
  return TESTBIT(fFlags, kPINDiodeMethodValid);    
}


// --------------------------------------------------------------------------
//
// 1) Return if the charge distribution is already succesfully fitted  
//    or if the histogram is empty
// 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid 
//    possible remaining cosmics to spoil the fit.
// 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
// 4) Fit the histograms with a Gaussian
// 5) In case of failure print out the fit results
// 6) Retrieve the results and store them in this class
// 7) Calculate the number of photo-electrons after the F-Factor method
// 8) Calculate the errors of the F-Factor method
//
// The fits are declared valid (fFitValid = kTRUE), if:
//
// 1) Pixel has a fitted charge greater than 5*PedRMS
// 2) Pixel has a fit error greater than 0.
// 3) Pixel has a fit Probability greater than 0.0001 
// 4) Pixel has a charge sigma bigger than its Pedestal RMS
// 5) If FitTimes is used, 
//    the mean arrival time is at least 1.0 slices from the used edge slices 
// (this stage is only performed in the times fit)
//
// If the histogram is empty, all values are set to -1.
//
// The conversion factor after the F-Factor method is declared valid, if:
//
// 1) fFitValid is kTRUE
// 2) Conversion Factor is bigger than 0.
// 3) The error of the conversion factor is smaller than 10%
//
Bool_t MCalibrationPix::FitCharge() 
{

  //
  // 1) Return if the charge distribution is already succesfully fitted  
  //    or if the histogram is empty
  //
  if (fHist->IsFitOK() || fHist->IsEmpty())
    return kTRUE;

  //
  // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid 
  //    possible remaining cosmics to spoil the fit.
  //
  //  if (fPed && fPedRms)
  //    fHist->SetLowerFitRange(1.5*fPedRms);
  //  else
  //    *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;

  //
  // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
  //
  if (fHist->UseLoGain())
    SetHiGainSaturation();
      
  //
  // 4) Fit the Lo Gain histograms with a Gaussian
  //
  if(fHist->FitCharge())
    {
      SETBIT(fFlags,kFitted);
    }
  else
    {
      *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl;
      //          
      // 5) In case of failure print out the fit results
      //
      //	  fHist->PrintChargeFitResult();
      CLRBIT(fFlags,kFitted);
    }

  //
  // 6) Retrieve the results and store them in this class
  //
  fCharge         = fHist->GetChargeMean();
  fErrCharge      = fHist->GetChargeMeanErr(); 
  fSigmaCharge    = fHist->GetChargeSigma();
  fErrSigmaCharge = fHist->GetChargeSigmaErr();
  fChargeProb     = fHist->GetChargeProb();

  if (CheckChargeFitValidity())
    SETBIT(fFlags,kFitValid);
  else
    {
      CLRBIT(fFlags,kFitValid);
      return kFALSE;
    }
  
  //
  // 7) Calculate the number of photo-electrons after the F-Factor method
  // 8) Calculate the errors of the F-Factor method
  //
  if ((fPed > 0.)  && (fPedRms > 0.))
    {
      
      //
      // Square all variables in order to avoid applications of square root
      //
      // First the relative error squares
      //
      const Float_t chargeSquare              =       fCharge*   fCharge;
      const Float_t chargeSquareRelErrSquare  = 4.*fErrCharge*fErrCharge / chargeSquare;

      const Float_t fFactorRelErrSquare       = fFactorError * fFactorError / (fFactor * fFactor);
      //
      // Now the absolute error squares
      //
      const Float_t sigmaSquare               =       fSigmaCharge*   fSigmaCharge;
      const Float_t sigmaSquareErrSquare      = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare;

      const Float_t elecRmsSquare             =       fElectronicPedRms*   fElectronicPedRms;
      const Float_t elecRmsSquareErrSquare    = 4.*fErrElectronicPedRms*fErrElectronicPedRms * elecRmsSquare;

      Float_t pedRmsSquare                    =       fPedRms*   fPedRms;
      Float_t pedRmsSquareErrSquare           = 4.*fErrPedRms*fErrPedRms * pedRmsSquare;

      if (TESTBIT(fFlags,kHiGainSaturation))
	{
          
          //
          // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it 
          // from the Hi Gain: 
          //
          // We extract the pure NSB contribution:
          //
          Float_t nsbSquare             =  pedRmsSquare          - elecRmsSquare;
          Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
                                          / (nsbSquare * nsbSquare) ;

	  if (nsbSquare < 0.)
            nsbSquare = 0.;
          
          //
          // Now, we divide the NSB by the conversion factor and 
          // add it quadratically to the electronic noise
          //
          const Float_t conversionSquare             =    fConversionHiLo     *fConversionHiLo;
          const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare;

          //
          // Calculate the new "Pedestal RMS"
          //
          const Float_t convertedNsbSquare          =  nsbSquare             / conversionSquare;
          const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
                                                      * convertedNsbSquare * convertedNsbSquare;

          pedRmsSquare           = convertedNsbSquare           + elecRmsSquare;
          pedRmsSquareErrSquare  = convertedNsbSquareErrSquare  + elecRmsSquareErrSquare;
          
	}    /* if (kHiGainSaturation) */
      
      //
      // Calculate the reduced sigmas
      //
      fRSigmaSquare = sigmaSquare - pedRmsSquare;
      if (fRSigmaSquare <= 0.)
        {
          *fLog << warn 
                << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel " 
                << fPixId << endl;
          if (TESTBIT(fFlags,kHiGainSaturation))
            ApplyLoGainConversion();
          return kFALSE;
        }

      const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
                                               / (fRSigmaSquare * fRSigmaSquare) ;
  
      //
      // Calculate the number of phe's from the F-Factor method
      // (independent on Hi Gain or Lo Gain)
      //
      fPheFFactorMethod = fFactor * chargeSquare / fRSigmaSquare;

      const Float_t pheFFactorRelErrSquare =  fFactorRelErrSquare
                                            + chargeSquareRelErrSquare
                                            + rSigmaSquareRelErrSquare ;
          
      fPheFFactorMethodError        =  TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;

      //
      // Calculate the conversion factors 
      //
      if (TESTBIT(fFlags,kHiGainSaturation))
        ApplyLoGainConversion();
      
      const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge);
      
      fConversionFFactorMethod        =  fPheFFactorMethod / fCharge ;
      fConversionErrorFFactorMethod   =  ( pheFFactorRelErrSquare + chargeRelErrSquare )
                                         * fConversionFFactorMethod * fConversionFFactorMethod;
      
      if ( IsFitValid()                     &&
           (fConversionFFactorMethod > 0.) &&
           (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) )
        SETBIT(fFlags,kFFactorMethodValid);
      else 
        CLRBIT(fFlags,kFFactorMethodValid);
      
    } /*   if ((fPed > 0.)  && (fPedRms > 0.)) */

  return kTRUE;
  
}

//
// The check return kTRUE if:
//
// 1) Pixel has a fitted charge greater than 5*PedRMS
// 2) Pixel has a fit error greater than 0.
// 3) Pixel has a fitted charge greater its charge error
// 4) Pixel has a fit Probability greater than 0.0001 
// 5) Pixel has a charge sigma bigger than its Pedestal RMS
// 
Bool_t MCalibrationPix::CheckChargeFitValidity()
{

  if (TESTBIT(fFlags,kExcludeQualityCheck))
    return kTRUE;

  Float_t equivpedestal = GetPedRms();

  if (TESTBIT(fFlags,kHiGainSaturation))
    equivpedestal /= fConversionHiLo;
      
  if (fCharge < fChargeLimit*equivpedestal)
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << fChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
      return kFALSE;
    }
  
  if (fErrCharge < fChargeErrLimit) 
    {
      *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
            << fChargeErrLimit << " in Pixel " << fPixId << endl;
      return kFALSE;
    }
      
  if (fCharge < fChargeRelErrLimit*fErrCharge) 
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << fChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
      return kFALSE;
    }
      
  if (!fHist->IsFitOK()) 
    {
      *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel " << fPixId << endl;
      return kFALSE;
    }

  if (fSigmaCharge < equivpedestal)
    {
      *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " << fPixId << endl;
      return kFALSE;
    }
  return kTRUE;
}

//
// The check returns kTRUE if:
//
// The mean arrival time is at least 1.0 slices from the used edge slices 
//
Bool_t MCalibrationPix::CheckTimeFitValidity()
{

  if (TESTBIT(fFlags,kExcludeQualityCheck))
    return kTRUE;

  Float_t lowerrange;
  Float_t upperrange;

  if (TESTBIT(fFlags,kHiGainSaturation))
    {
      lowerrange = (Float_t)fHist->GetTimeLowerFitRangeLoGain()+1.;
      upperrange = (Float_t)fHist->GetTimeUpperFitRangeLoGain()+1.;
    }
  else
    {
      lowerrange = (Float_t)fHist->GetTimeLowerFitRangeHiGain()+1.;
      upperrange = (Float_t)fHist->GetTimeUpperFitRangeHiGain()+1.;
    }


  if (fTime < lowerrange)
    {
      *fLog << warn 
            << "WARNING: Mean Fitted Time inside or smaller than first used FADC slice in Pixel " 
            << fPixId << " time: " << fTime << " Range: " << lowerrange << endl;
      return kFALSE;
    }

  if (fTime > upperrange)
    {
      *fLog << warn 
            << "WARNING: Mean Fitted Time inside or greater than last used FADC slice in Pixel " 
            << fPixId << " time: " << fTime << " Range: " << upperrange << endl;
      return kFALSE;
    }

  return kTRUE;
}


//
// The check returns kTRUE if:
//
// 
//
Bool_t MCalibrationPix::CheckOscillations()
{


  return kTRUE;
}



void MCalibrationPix::ApplyLoGainConversion()
{
  
  const Float_t chargeRelErrSquare     =   fErrCharge*fErrCharge
                                         /(   fCharge *  fCharge);
  const Float_t sigmaRelErrSquare      =   fErrSigmaCharge*fErrSigmaCharge
                                         /(   fSigmaCharge *  fSigmaCharge);
  const Float_t conversionRelErrSquare =   fConversionHiLoError*fConversionHiLoError 
                                         /(fConversionHiLo    * fConversionHiLo);
  
  fCharge         *= fConversionHiLo;
  fErrCharge       = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
  
  fSigmaCharge    *= fConversionHiLo;
  fErrSigmaCharge =  TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
  
}



// --------------------------------------------------------------------------
//
// 1) Fit the arrival times
// 2) Retrieve the results
// 3) Note that because of the low number of bins, the NDf is sometimes 0, so 
//    Root does not give a reasonable Probability, the Chisquare is more significant
//
// This fit has to be done AFTER the Charges fit, 
// otherwise only the Hi Gain will be fitted, even if there are no entries
//
// 
Bool_t MCalibrationPix::FitTime() 
{

  //
  // Fit the Low Gain
  //
  if (TESTBIT(fFlags,kHiGainSaturation))
    {
      if(!fHist->FitTimeLoGain())
	{
	  *fLog << warn << "WARNING: Could not fit Lo Gain times of pixel " << fPixId << endl;
          //	  fHist->PrintTimeFitResult();
	  return kFALSE;
	}
    }

  //
  // Fit the High Gain
  //
  else
    {
      if(!fHist->FitTimeHiGain())
	{
	  *fLog << warn << "WARNING: Could not fit Hi Gain times of pixel " << fPixId << endl;
          //	  fHist->PrintTimeFitResult();
	  return kFALSE;
	}
    }
    
  fTime          = fHist->GetTimeMean();
  fSigmaTime     = fHist->GetTimeSigma();
  fTimeChiSquare = fHist->GetTimeChiSquare();
  fTimeProb      = fHist->GetTimeProb();

  if (CheckTimeFitValidity())
    SETBIT(fFlags,kFitValid);
  else
    CLRBIT(fFlags,kFitValid);

  return kTRUE;
}

