/* ======================================================================== *\
!
! *
! * 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. 
//
MCalibrationPix::MCalibrationPix(const char *name, const char *title)
    : fPixId(-1),
      fCharge(-1.),
      fErrCharge(-1.),
      fSigmaCharge(-1.),
      fErrSigmaCharge(-1.),
      fRSigmaSquare(-1.),
      fChargeProb(-1.),
      fPed(-1.),
      fPedRms(-1.),
      fErrPedRms(0.),
      fElectronicPedRms(1.5),
      fErrElectronicPedRms(0.3),
      fTime(-1.),
      fSigmaTime(-1.),
      fTimeChiSquare(-1.),
      fFactor(1.3),
      fPheFFactorMethod(-1.),
      fConversionFFactorMethod(-1.),
      fConversionBlindPixelMethod(-1.),
      fConversionPINDiodeMethod(-1.),
      fConversionErrorFFactorMethod(-1.),
      fConversionErrorBlindPixelMethod(-1.),
      fConversionErrorPINDiodeMethod(-1.),
      fConversionSigmaFFactorMethod(-1.),
      fConversionSigmaBlindPixelMethod(-1.),
      fConversionSigmaPINDiodeMethod(-1.),
      fHiGainSaturation(kFALSE),
      fFitValid(kFALSE),        
      fFitted(kFALSE), 
      fBlindPixelMethodValid(kFALSE),
      fFFactorMethodValid(kFALSE),
      fPINDiodeMethodValid(kFALSE)
{

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

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

}

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


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


// ------------------------------------------------------------------------
//
// Invalidate values
//
void MCalibrationPix::Clear(Option_t *o)
{
  fHist->Reset();
}


// --------------------------------------------------------------------------
//
// 1) Return if the charge distribution is already succesfully fitted  
// 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 fit is declared as valid, if:
//
// 1) Pixel has a fitted charge greater than 3*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.5 slices from the edge
// 
// The F-Factor method is declared as valid, if:
//
// 1) Pixel is FitValid
// 2) Conversion Factor is bigger than 0.
// 3) The error of the conversion factor is smaller than 10%
//
Bool_t MCalibrationPix::FitCharge() 
{

  if (fHist->IsFitOK())
    return kTRUE;

  if (fPed && fPedRms)
    fHist->SetLowerFitRange(1.5*fPedRms);
  else
    *fLog << warn << "Cannot set lower fit range: Pedestals not available" << endl;

  if (fHist->UseLoGain())
    {

      SetHiGainSaturation();

      if(!fHist->FitChargeLoGain())
	{
	  *fLog << warn << "Could not fit Lo Gain charges of pixel " << fPixId << endl;
	  fHist->PrintChargeFitResult();
	}
    }
  else 
    {
      if(!fHist->FitChargeHiGain())
	{
	  *fLog << warn << "Could not fit Hi Gain charges of pixel " << fPixId << endl;
	  fHist->PrintChargeFitResult();
	}
    }
  

  fCharge         = fHist->GetChargeMean();
  fErrCharge      = fHist->GetChargeMeanErr(); 
  fSigmaCharge    = fHist->GetChargeSigma();
  fErrSigmaCharge = fHist->GetChargeSigmaErr();
  fChargeProb     = fHist->GetChargeProb();

  if (fCharge <= 0.)
    {
      *fLog << warn << "Cannot apply calibration: Mean Fitted Charges are smaller than 0 in pixel " 
            << fPixId << endl;
      return kFALSE;
    }

  if (fErrCharge > 0.) 
    fFitted = kTRUE;

  if ( (fCharge > 3.*GetPedRms()) &&
       (fErrCharge  > 0.)        &&
       (fHist->IsFitOK())        &&
       (fSigmaCharge > fPedRms)  &&
       (fTime > fHist->GetTimeLowerFitRange()+1.) &&
       (fTime < fHist->GetTimeUpperFitRange()-1.) )
    fFitValid =  kTRUE;

  //
  // Apply the F-Factor Method
  //
  if ((fPed > 0.)  && (fPedRms > 0.))
    {
      
      //
      // Square all variables in order to avoid applications of square root
      //
      const Float_t chargesquare              =       fCharge*fCharge;
      const Float_t chargessquarerelerrsquare = 4.*fErrCharge*fErrCharge/chargesquare;

      const Float_t sigmasquare         =       fSigmaCharge*fSigmaCharge;
      const Float_t sigmasquareerr      = 2.*fErrSigmaCharge*fSigmaCharge;

      const Float_t pedrmssquare        =       fPedRms*fPedRms;
      const Float_t pedrmssquareerr     = 2.*fErrPedRms*fPedRms;

      const Float_t elecrmssquare       =       fElectronicPedRms*fElectronicPedRms;
      const Float_t elecrmssquareerr    = 2.*fErrElectronicPedRms*fElectronicPedRms;

      const Float_t conversionsquare    =    gkConversionHiLo     *gkConversionHiLo;
      const Float_t conversionsquareerr = 2.*gkConversionHiLoError*gkConversionHiLo;

      const Float_t ffactorrelerrsquare = fFactorError * fFactorError / fFactor / fFactor;

      Float_t rsigmasquarerelerrsquare  = 0.;
      Float_t   pheffactorrelerrsquare  = 0.;

      if (fHiGainSaturation)
	{

          //
          // 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 nsbsquareerrsquare   = pedrmssquareerr  * pedrmssquareerr +
                                         elecrmssquareerr * elecrmssquareerr;

	  if (nsbsquare < 0.)
            nsbsquare = 0.;
          
          //
          // Now, we divide the NSB by the conversion factor and 
          // add it quadratically to the electronic noise
          //
          Float_t logainrmssquare          = nsbsquare/conversionsquare + elecrmssquare;
          Float_t logainrmssquareerrsquare = nsbsquareerrsquare/conversionsquare/conversionsquare 
                                           + elecrmssquareerr * elecrmssquareerr 
                                           + conversionsquareerr*conversionsquareerr
                                             / (conversionsquare*conversionsquare)
                                             / (conversionsquare*conversionsquare);
          //
          // Calculate the reduced sigma with the new "Pedestal RMS"
          //
          fRSigmaSquare            = sigmasquare - logainrmssquare;
          rsigmasquarerelerrsquare = (sigmasquareerr * sigmasquareerr + logainrmssquareerrsquare)
                                     / (fRSigmaSquare * fRSigmaSquare);

          if (fRSigmaSquare > 0.)
            {
              fPheFFactorMethod      = fFactor * chargesquare / fRSigmaSquare;
              pheffactorrelerrsquare =        ffactorrelerrsquare
                                      + chargessquarerelerrsquare
                                      +  rsigmasquarerelerrsquare ;
            }
          
	}    /* if (fHiGainSaturation) */
      else   
        {
          fRSigmaSquare            = sigmasquare - pedrmssquare;
          rsigmasquarerelerrsquare = (  sigmasquareerr *  sigmasquareerr 
                                     + pedrmssquareerr * pedrmssquareerr) 
                                   / (fRSigmaSquare * fRSigmaSquare);

          fPheFFactorMethod      = fFactor * chargesquare / fRSigmaSquare;
          //
          // We first calculate the squared relative error in order to save the 
          // TMath::Sqrt
          //
          pheffactorrelerrsquare = (       ffactorrelerrsquare
                                   + chargessquarerelerrsquare
                                   +  rsigmasquarerelerrsquare );
          

        }   /* if (fHiGainSaturation) */
      

      if (fPheFFactorMethod <= 0.)
        {
          *fLog << warn << "Cannot apply F-Factor calibration: Number of PhE smaller than 0 in pixel " 
                << fPixId << endl;
          return kFALSE;
        }
      
      fConversionFFactorMethod      =  fPheFFactorMethod /  fCharge ;
      fConversionErrorFFactorMethod =  fConversionFFactorMethod  *
                                       ( pheffactorrelerrsquare +
                                       (fErrCharge * fErrCharge / chargesquare ) );

      fPheFFactorMethodError        =  TMath::Sqrt(pheffactorrelerrsquare) * fPheFFactorMethod;

      if ( IsFitValid()                    &&
          (fConversionFFactorMethod > 0.) &&
           (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) )
        fFFactorMethodValid = kTRUE;
      

    } /*   if ((fPed > 0.)  && (fPedRms > 0.)) */

  return kTRUE;
  
}


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

  fPed    = ped;    
  fPedRms = pedrms;
  
}

// --------------------------------------------------------------------------
//
// 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 (fHiGainSaturation)
    {
      if(!fHist->FitTimeLoGain())
	{
	  *fLog << warn << "Could not fit Lo Gain times of pixel " << fPixId << endl;
	  fHist->PrintTimeFitResult();
	  return kFALSE;
	}
    }

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

  return kTRUE;
}

