/* ======================================================================== *\
!
! *
! * 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   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MCalibrationChargePix                                                   //
//                                                                         //
// Storage container to hold informations about the calibration values     //
// values of one Pixel (PMT).                                              //
//                                                                         //
// 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.
//   with the Munich definition of the F-Factor, thus:
//   F = Sigma(Out)/Mean(Out) * Mean(In)/Sigma(In)
//   Mean F-Factor  = 1.15
//   Error F-Factor = 0.02
//
// - Average QE: (email David Paneque, 14.2.04):
//
//  The conversion factor that comes purely from QE folded to a Cherenkov
//  spectrum has to be multiplied by:
//  * Plexiglass window -->> 0.96 X 0.96
//  * PMT photoelectron collection efficiency -->> 0.9
//  * Light guides efficiency -->> 0.94
//
//  Concerning the light guides effiency estimation... Daniel Ferenc 
//  is preparing some work (simulations) to estimate it. Yet so far, he has 
//  been busy with other stuff, and this work is still UNfinished.
//
//  The estimation I did comes from:
//  1) Reflectivity of light guide walls is 85 % (aluminum)
//  2) At ZERO degree light incidence, 37% of the light hits such walls 
//    (0.15X37%= 5.6% of light lost)
//  3) When increasing the light incidence angle, more and more light hits 
//     the walls.
//
//  However, the loses due to larger amount of photons hitting the walls is more 
//  or less counteracted by the fact that more and more photon trajectories cross 
//  the PMT photocathode twice, increasing the effective sensitivity of the PMT.
//
//  Jurgen Gebauer did some quick measurements about this issue. I attach a 
//  plot. You can see that the angular dependence is (more or less) in agreement 
//  with a CosTheta function (below 20-25 degrees),
//  which is the variation of teh entrance window cross section. So, in 
//  first approximation, no loses when increasing light incidence angle; 
//  and therefore, the factor 0.94.
//
//  So, summarizing... I would propose the following conversion factors 
//  (while working with CT1 cal box) in order to get the final number of photons 
//  from the detected measured size in ADC counts.
// 
//  Nph = ADC * FmethodConversionFactor / ConvPhe-PhFactor
// 
//  FmethodConversionFactor ; measured for individual pmts
// 
//  ConvPhe-PhFactor = 0.98 * 0.23 * 0.90 * 0.94 * 0.96 * 0.96 = 0.18
// 
//  I would not apply any smearing of this factor (which we have in nature), 
//  since we might be applying it to PMTs in the totally wrong direction.
// 
//
//  Error of all variables are calculated by error-propagation. Note that internally, 
//  all error variables contain Variances in order to save the CPU-intensive square rooting 
// 
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargePix.h"

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

#include "MBadPixelsPix.h"

ClassImp(MCalibrationChargePix);

using namespace std;

const Float_t MCalibrationChargePix::gkElectronicPedRms         = 1.5;
const Float_t MCalibrationChargePix::gkElectronicPedRmsErr      = 0.3;
const Float_t MCalibrationChargePix::gkFFactor                  = 1.15;
const Float_t MCalibrationChargePix::gkFFactorErr               = 0.02;

const Float_t MCalibrationChargePix::gkConversionHiLo           = 10.;
const Float_t MCalibrationChargePix::gkConversionHiLoErr        = 2.5;

const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit    = 5.;
// --------------------------------------------------------------------------
//
// Default Constructor: 
//
MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
    : fPixId(-1),
      fFlags(0)
{

  fName  = name  ? name  : "MCalibrationChargePix";
  fTitle = title ? title : "Container of the fit results of MHCalibrationChargePixs ";

  Clear();

  //
  // At the moment, we don't have a database, yet, 
  // so we get it from the configuration file
  //
  SetConversionHiLo();
  SetConversionHiLoErr();

  SetPheFFactorMethodLimit();
  
}

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

  SetHiGainSaturation       ( kFALSE );
  SetExcluded               ( kFALSE );
  SetBlindPixelMethodValid  ( kFALSE );
  SetFFactorMethodValid     ( kFALSE );
  SetPINDiodeMethodValid    ( kFALSE );
  SetCombinedMethodValid    ( kFALSE );

  fHiGainMeanCharge                 =  -1.;
  fHiGainMeanChargeVar              =  -1.;
  fHiGainSigmaCharge                =  -1.;
  fHiGainSigmaChargeVar             =  -1.;
  fHiGainChargeProb                 =  -1.;

  fLoGainMeanCharge                 =  -1.;
  fLoGainMeanChargeVar              =  -1.;
  fLoGainSigmaCharge                =  -1.;
  fLoGainSigmaChargeVar             =  -1.;
  fLoGainChargeProb                 =  -1.;

  fRSigmaCharge                     =  -1.;
  fRSigmaChargeVar                  =  -1.;
  
  fHiGainNumPickup                  =  -1;
  fLoGainNumPickup                  =  -1;

  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fPedVar                           =  -1.;

  fLoGainPedRms                     =  -1.;
  fLoGainPedRmsVar                  =  -1.;

  fAbsTimeMean                      =  -1.;
  fAbsTimeRms                       =  -1.;

  fPheFFactorMethod                 =  -1.;
  fPheFFactorMethodVar              =  -1.;

  fMeanConversionFFactorMethod      =  -1.;
  fMeanConversionBlindPixelMethod   =  -1.;
  fMeanConversionPINDiodeMethod     =  -1.;
  fMeanConversionCombinedMethod     =  -1.;

  fConversionFFactorMethodVar       =  -1.;
  fConversionBlindPixelMethodVar    =  -1.;
  fConversionPINDiodeMethodVar      =  -1.;
  fConversionCombinedMethodVar      =  -1.;

  fSigmaConversionFFactorMethod     =  -1.;
  fSigmaConversionBlindPixelMethod  =  -1.;
  fSigmaConversionPINDiodeMethod    =  -1.;
  fSigmaConversionCombinedMethod    =  -1.;

  fTotalFFactorFFactorMethod        =  -1.;
  fTotalFFactorBlindPixelMethod     =  -1.;
  fTotalFFactorPINDiodeMethod       =  -1.;
  fTotalFFactorCombinedMethod       =  -1.;

}


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

  fPed       = ped;    
  fPedRms    = pedrms;
  fPedVar    = pederr*pederr;
}
  

void  MCalibrationChargePix::SetMeanCharge( const Float_t f ) 
{
    if (IsHiGainSaturation())
	fLoGainMeanCharge = f;
    else
	fHiGainMeanCharge = f;
}

void  MCalibrationChargePix::SetMeanChargeErr( const Float_t f ) 
{
    if (IsHiGainSaturation())
	fLoGainMeanChargeVar = f*f;
    else
	fHiGainMeanChargeVar = f*f;

}

void  MCalibrationChargePix::SetSigmaCharge( const Float_t f )  
{
    if (IsHiGainSaturation())
	fLoGainSigmaCharge = f;
    else
	fHiGainSigmaCharge = f;
}


void  MCalibrationChargePix::SetSigmaChargeErr( const Float_t f )   
{
    if (IsHiGainSaturation())
	fLoGainSigmaChargeVar = f*f;
    else
	fHiGainSigmaChargeVar = f*f;

}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionFFactorMethod  = c;
  fConversionFFactorMethodVar   = err*err;
  fSigmaConversionFFactorMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationChargePix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionCombinedMethod  = c;
  fConversionCombinedMethodVar   = err*err;
  fSigmaConversionCombinedMethod = sig;
}


// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationChargePix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionBlindPixelMethod  = c;
  fConversionBlindPixelMethodVar   = err*err;
  fSigmaConversionBlindPixelMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionPINDiodeMethod  = c ;
  fConversionPINDiodeMethodVar   = err*err;
  fSigmaConversionPINDiodeMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the Hi Gain Saturation Bit from outside
//
void MCalibrationChargePix::SetHiGainSaturation(Bool_t b)
{
    b ?  SETBIT(fFlags, kHiGainSaturation) : CLRBIT(fFlags, kHiGainSaturation); 
}


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

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

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

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

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )  
{ 
  b ?  SETBIT(fFlags, kCombinedMethodValid) : CLRBIT(fFlags, kCombinedMethodValid); 
}


Float_t  MCalibrationChargePix::GetPedRms()  const
{
  return IsHiGainSaturation() ? fLoGainPedRms : fPedRms;
}

Float_t  MCalibrationChargePix::GetPedRmsErr()  const
{
  return IsHiGainSaturation() ? TMath::Sqrt(fLoGainPedRmsVar) : TMath::Sqrt(fPedVar)/2.;
}

Float_t  MCalibrationChargePix::GetPedErr()  const
{
  return  TMath::Sqrt(fPedVar);
}

Float_t  MCalibrationChargePix::GetMeanCharge()   const 
{
    return IsHiGainSaturation() ? GetLoGainMeanCharge() : GetHiGainMeanCharge() ;
}

Float_t  MCalibrationChargePix::GetMeanChargeErr()   const 
{
    return IsHiGainSaturation() ? GetLoGainMeanChargeErr() : GetHiGainMeanChargeErr() ;
}

Float_t  MCalibrationChargePix::GetChargeProb()   const 
{
    return IsHiGainSaturation() ? fLoGainChargeProb : fHiGainChargeProb ;
}

Float_t  MCalibrationChargePix::GetSigmaCharge()   const 
{
    return IsHiGainSaturation() ? GetLoGainSigmaCharge() : GetHiGainSigmaCharge() ;
}

Float_t  MCalibrationChargePix::GetSigmaChargeErr()   const 
{
    return IsHiGainSaturation() ? GetLoGainSigmaChargeErr() : GetHiGainSigmaChargeErr() ;
}

Float_t MCalibrationChargePix::GetHiGainMeanChargeErr()  const
{
  return TMath::Sqrt(fHiGainMeanChargeVar);
}

Float_t MCalibrationChargePix::GetLoGainMeanCharge()  const 
{
  return fLoGainMeanCharge * fConversionHiLo;
}

Float_t MCalibrationChargePix::GetLoGainMeanChargeErr()  const
{

  const Float_t chargeRelVar     =  fLoGainMeanChargeVar
                                 /( fLoGainMeanCharge * fLoGainMeanCharge );

  const Float_t conversionRelVar =  fConversionHiLoVar
                                 /( fConversionHiLo   * fConversionHiLo   );

  return TMath::Sqrt(chargeRelVar+conversionRelVar) * GetLoGainMeanCharge();
}

Float_t MCalibrationChargePix::GetLoGainSigmaCharge()  const 
{
  return fLoGainSigmaCharge * fConversionHiLo;
}

Float_t MCalibrationChargePix::GetLoGainSigmaChargeErr()  const
{

  const Float_t sigmaRelVar     =  fLoGainSigmaChargeVar
                                /( fLoGainSigmaCharge * fLoGainSigmaCharge );

  const Float_t conversionRelVar =  fConversionHiLoVar
                                 /( fConversionHiLo   * fConversionHiLo   );

  return TMath::Sqrt(sigmaRelVar+conversionRelVar) * GetLoGainSigmaCharge();
}

Float_t MCalibrationChargePix::GetHiGainSigmaChargeErr()  const
{
  return TMath::Sqrt(fHiGainSigmaChargeVar);
}

Float_t MCalibrationChargePix::GetRSigmaCharge()  const
{
  return IsHiGainSaturation() ? fRSigmaCharge*fConversionHiLo : fRSigmaCharge ;
} 

Float_t MCalibrationChargePix::GetRSigmaChargeErr()  const
{
 if (IsHiGainSaturation())
    {
      const Float_t rsigmaRelVar     =  fRSigmaChargeVar
                                    /( fRSigmaCharge * fRSigmaCharge );
      const Float_t conversionRelVar =  fConversionHiLoVar
                                     /( fConversionHiLo   * fConversionHiLo   );
      return TMath::Sqrt(rsigmaRelVar+conversionRelVar) * GetRSigmaCharge();
    }
 else
   return TMath::Sqrt(fRSigmaChargeVar);

}

Float_t MCalibrationChargePix::GetConversionHiLoErr()  const
{
  if (fConversionHiLoVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionHiLoVar);
}

Float_t MCalibrationChargePix::GetPheFFactorMethodErr()  const
{
  if (fPheFFactorMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fPheFFactorMethodVar);
}

Float_t MCalibrationChargePix::GetConversionCombinedMethodErr()  const
{
  if (fConversionCombinedMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionCombinedMethodVar);
}

Float_t MCalibrationChargePix::GetConversionPINDiodeMethodErr()  const
{
  if (fConversionPINDiodeMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionPINDiodeMethodVar);
}

Float_t MCalibrationChargePix::GetConversionBlindPixelMethodErr()  const
{
  if (fConversionBlindPixelMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionBlindPixelMethodVar);
}

Float_t MCalibrationChargePix::GetConversionFFactorMethodErr()  const
{
  if (fConversionFFactorMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionFFactorMethodVar);
}

Float_t MCalibrationChargePix::GetTotalFFactorCombinedMethodErr()  const
{
  if (fTotalFFactorCombinedMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fTotalFFactorCombinedMethodVar);
}

Float_t MCalibrationChargePix::GetTotalFFactorPINDiodeMethodErr()  const
{
  if (fTotalFFactorPINDiodeMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fTotalFFactorPINDiodeMethodVar);
}

Float_t MCalibrationChargePix::GetTotalFFactorBlindPixelMethodErr()  const
{
  if (fTotalFFactorBlindPixelMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fTotalFFactorBlindPixelMethodVar);
}

Float_t MCalibrationChargePix::GetTotalFFactorFFactorMethodErr()  const
{
  if (fTotalFFactorFFactorMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fTotalFFactorFFactorMethodVar);
}

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

Bool_t MCalibrationChargePix::IsHiGainSaturation()    const
{ 
   return TESTBIT(fFlags,kHiGainSaturation);  
}

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

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

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

Bool_t MCalibrationChargePix::IsCombinedMethodValid()  const
{ 
  return TESTBIT(fFlags, kCombinedMethodValid);    
}

//
// 
//
Bool_t MCalibrationChargePix::CalcReducedSigma()
{

  const Float_t sigmacharge    = IsHiGainSaturation() ? fLoGainSigmaCharge    : fHiGainSigmaCharge   ;
  const Float_t sigmachargevar = IsHiGainSaturation() ? fLoGainSigmaChargeVar : fHiGainSigmaChargeVar;

  const Float_t sigmaSquare     =      sigmacharge     * sigmacharge;
  const Float_t sigmaSquareVar  = 4.*  sigmachargevar  * sigmaSquare;

  Float_t pedRmsSquare ;         
  Float_t pedRmsSquareVar;
  
  if (IsHiGainSaturation())
    {
      pedRmsSquare     =      fLoGainPedRms    * fLoGainPedRms;
      pedRmsSquareVar  =  4.* fLoGainPedRmsVar * pedRmsSquare;
    }
  else
  {
      pedRmsSquare       =  fPedRms * fPedRms;                                          
      pedRmsSquareVar    =  fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
  }
  //
  // Calculate the reduced sigmas
  //
  const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
  if (rsigmachargesquare <= 0.)
    {
      *fLog << warn 
            << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " 
            << fPixId << endl;
      return kFALSE;
    }
  

  fRSigmaCharge    = TMath::Sqrt(rsigmachargesquare);
  fRSigmaChargeVar = 0.25 * (sigmaSquareVar + pedRmsSquareVar) / rsigmachargesquare;

  return kTRUE;
}

//
// Calculate the number of photo-electrons after the F-Factor method
// Calculate the errors of the F-Factor method
//
Bool_t MCalibrationChargePix::CalcFFactorMethod()
{

  if (fRSigmaCharge < 0.)
    {
      SetFFactorMethodValid(kFALSE);
      return kFALSE;
    }
  
  const Float_t charge    = IsHiGainSaturation() ? fLoGainMeanCharge    : fHiGainMeanCharge   ;
  const Float_t chargevar = IsHiGainSaturation() ? fLoGainMeanChargeVar : fHiGainMeanChargeVar;

  //
  // Square all variables in order to avoid applications of square root
  //
  // First the relative error squares
  //
  const Float_t chargeSquare        =     charge * charge;
  const Float_t chargeSquareRelVar  = 4.* chargevar/ chargeSquare;

  const Float_t ffactorsquare       =    gkFFactor    * gkFFactor;
  const Float_t ffactorsquareRelVar = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;

  const Float_t rsigmaSquare        =     fRSigmaCharge    * fRSigmaCharge;
  const Float_t rsigmaSquareRelVar  = 4.* fRSigmaChargeVar / rsigmaSquare;

  //
  // Calculate the number of phe's from the F-Factor method
  // (independent on Hi Gain or Lo Gain)
  //
  fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmaSquare;

  if (fPheFFactorMethod < fPheFFactorMethodLimit)
    {
      SetFFactorMethodValid(kFALSE);
      return kFALSE;
    }
  
  //
  // Calculate the Error of Nphe
  //
  fPheFFactorMethodVar =  (ffactorsquareRelVar + chargeSquareRelVar + rsigmaSquareRelVar) 
                         * fPheFFactorMethod * fPheFFactorMethod;

  SetFFactorMethodValid(kTRUE);
  return kTRUE;
}


void MCalibrationChargePix::CalcLoGainPedestal(Float_t logainsamples)
{

  fElectronicPedRms       = gkElectronicPedRms    * TMath::Sqrt(logainsamples);
  fElectronicPedRmsVar    = gkElectronicPedRmsErr * gkElectronicPedRmsErr * logainsamples;
  
  Float_t pedRmsSquare      = fPedRms * fPedRms;
  Float_t pedRmsSquareVar   = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
  
  //
  // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it 
  // from the HI GAIN (all calculation per slice up to now):  
  //
  // We extract the pure NSB contribution:
  //
  const Float_t elecRmsSquare    =    fElectronicPedRms    * fElectronicPedRms;
  const Float_t elecRmsSquareVar = 4.*fElectronicPedRmsVar * elecRmsSquare;
  
  Float_t nsbSquare             =  pedRmsSquare    - elecRmsSquare;
  Float_t nsbSquareRelVar       = (pedRmsSquareVar + elecRmsSquareVar)
                                 / (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 convertedNsbSquare      =    nsbSquare       / conversionSquare;
  const Float_t convertedNsbSquareVar   =    nsbSquareRelVar
                                	    * convertedNsbSquare * convertedNsbSquare;
    
  pedRmsSquare     = convertedNsbSquare    + elecRmsSquare;
  pedRmsSquareVar  = convertedNsbSquareVar + elecRmsSquareVar;
  
  fLoGainPedRms    = TMath::Sqrt(pedRmsSquare);
  fLoGainPedRmsVar = 0.25 * pedRmsSquareVar /  pedRmsSquare;
  
}
 
 
 
