/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                               
// MCalibrationChargeBlindPix    
//
// Storage container of the fit results of the Blind Pixel signal 
// (from MHCalibrationChargeBlindPix). 
//
// The Flux is calculated in photons per mm^2 in the camera plane. 
//
// Currently, the following numbers are implemented:
// - gkBlindPixelArea: 100 mm^2
// - Average QE of Blind Pixel: 
//    gkBlindPixelQEGreen: 0.154
//    gkBlindPixelQEBlue : 0.226
//    gkBlindPixelQEUV   : 0.247
//    gkBlindPixelQECT1  : 0.247
// - Average QE Error of Blind Pixel: 
//    gkBlindPixelQEGreenErr: 0.015;
//    gkBlindPixelQEBlueErr : 0.02;
//    gkBlindPixelQEUVErr   : 0.02;
//    gkBlindPixelQECT1Err  : 0.02;
// - Attenuation factor Blind Pixel:
//    gkBlindPixelAttGreen :  1.97;
//    gkBlindPixelAttBlue  :  1.96;
//    gkBlindPixelAttUV    :  1.95;
//    gkBlindPixelAttCT1   :  1.95;
//
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeBlindPix.h"
#include "MCalibrationCam.h"

#include <TH1.h>

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

ClassImp(MCalibrationChargeBlindPix);

using namespace std;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelArea       = 100;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttGreen   = 1.97;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttBlue    = 1.96;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttUV      = 1.95;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttCT1     = 1.95;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedGreen    = 0.154;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedBlue     = 0.226;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedUV       = 0.247;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedCT1      = 0.247;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedGreenErr = 0.005;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedBlueErr  = 0.007;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedUVErr    = 0.01;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedCT1Err   = 0.01;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedGreen      = 0.192;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedBlue       = 0.27;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedUV         = 0.285;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedCT1        = 0.285;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedGreenErr   = 0.007;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedBlueErr    = 0.01;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedUVErr      = 0.012;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedCT1Err     = 0.012;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelCollectionEff      = 0.95; 
const Float_t MCalibrationChargeBlindPix::gkBlindPixelCollectionEffErr   = 0.02; 
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Calls:
// - Clear()
// - SetCoated()
//
MCalibrationChargeBlindPix::MCalibrationChargeBlindPix(const char *name, const char *title)
{

  fName  = name  ? name  : "MCalibrationChargeBlindPix";
  fTitle = title ? title : "Container of the fit results of the blind pixel";

  SetCoated();
  Clear();
}


// ------------------------------------------------------------------------
//
// Sets:
// - all flags to kFALSE
// - all variables to -1.
// - the fColor to MCalibrationCam::kNONE
//
// Calls: 
// - MCalibrationChargePix::Clear()
//
void MCalibrationChargeBlindPix::Clear(Option_t *o)
{

  fFluxInsidePlexiglass    = -1.;
  fFluxInsidePlexiglassVar = -1.;
  fLambda                  = -1.;
  fLambdaCheck             = -1.;
  fLambdaVar               = -1.;
  fMu0                     = -1.;
  fMu0Err                  = -1.;
  fMu1                     = -1.;
  fMu1Err                  = -1.;
  fSigma0                  = -1.;
  fSigma0Err               = -1.;
  fSigma1                  = -1.;
  fSigma1Err               = -1.;

  SetOscillating                   ( kFALSE );
  SetExcluded                      ( kFALSE );
  SetChargeFitValid                ( kFALSE );
  SetPedestalFitOK                 ( kFALSE );
  SetSinglePheFitOK                ( kFALSE );
  SetFluxInsidePlexiglassAvailable ( kFALSE );
  
  SetColor(MCalibrationCam::kNONE);

  MCalibrationChargePix::Clear();
}

void  MCalibrationChargeBlindPix::SetFluxInsidePlexiglassAvailable( const Bool_t b)
{
    b ? SETBIT(fFlags,kFluxInsidePlexiglassAvailable) : CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);
}

// --------------------------------------------------------------------------
//
// Set the Coated Bit from outside 
//
void  MCalibrationChargeBlindPix::SetCoated( const Bool_t b)
{
    b ? SETBIT(fFlags,kCoated) : CLRBIT(fFlags,kCoated);
}

// --------------------------------------------------------------------------
//
// Set the Oscillating Bit from outside 
//
void  MCalibrationChargeBlindPix::SetOscillating( const Bool_t b)
{
    b ? SETBIT(fFlags,kOscillating) : CLRBIT(fFlags,kOscillating);
}

// --------------------------------------------------------------------------
//
// Set the ChargeFitValid Bit from outside 
//
void  MCalibrationChargeBlindPix::SetChargeFitValid( const Bool_t b)
{
    b ? SETBIT(fFlags,kChargeFitValid) : CLRBIT(fFlags,kChargeFitValid);
}

// --------------------------------------------------------------------------
//
// Set the PedestalFitValid Bit from outside 
//
void  MCalibrationChargeBlindPix::SetPedestalFitOK( const Bool_t b)
{
    b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
}

// --------------------------------------------------------------------------
//
// Set the SinglePheFitValid Bit from outside 
//
void  MCalibrationChargeBlindPix::SetSinglePheFitOK( const Bool_t b)
{
    b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
}

// --------------------------------------------------------------------------
//
// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
// Return square root of fFluxInsidePlexiglassVar
// 
Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassErr() const
{
  if (fFluxInsidePlexiglassVar < 0.)
    return -1.;
  
  return TMath::Sqrt(fFluxInsidePlexiglassVar);
}

// --------------------------------------------------------------------------
//
// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
// Return -1 if fFluxInsidePlexiglass    is 0.
// Return fFluxInsidePlexiglassVar / fFluxInsidePlexiglass^2
// 
Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassRelVar() const
{
  if (fFluxInsidePlexiglassVar < 0.)
    return -1.;

  if (fFluxInsidePlexiglass == 0.)
    return -1.;
  
  return fFluxInsidePlexiglassVar / (fFluxInsidePlexiglass * fFluxInsidePlexiglass) ;
}

// --------------------------------------------------------------------------
//
// Return -1 if fLambdaVar is smaller than 0.
// Return square root of fLambdaVar
// 
Float_t MCalibrationChargeBlindPix::GetLambdaErr() const
{
  if (fLambdaVar < 0.)
    return -1.;
  
  return TMath::Sqrt(fLambdaVar);
}

// --------------------------------------------------------------------------
//
// Return -1 if fLambdaVar is smaller than 0.
// Return -1 if fLambda    is 0.
// Return fLambdaVar / (fLambda * fLambda )
// 
Float_t MCalibrationChargeBlindPix::GetLambdaRelVar() const
{
  if (fLambdaVar < 0.)
    return -1.;
  
  if (fLambda  == 0.)
    return -1.;
  
  return fLambdaVar / fLambda / fLambda ;
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEGreenErr is smaller than 0.
// Return -1 if gkBlindPixelQEGreen    is 0.
// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreen() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedGreen < 0.)
        return -1.;
      
      return gkBlindPixelQECoatedGreen;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedGreen < 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedGreen;
    }
  
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEBlueErr is smaller than 0.
// Return -1 if gkBlindPixelQEBlue    is 0.
// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlue() const
{
  if (IsCoated())
    {
      if (gkBlindPixelQECoatedBlue < 0.)
        return -1.;
      
      return gkBlindPixelQECoatedBlue;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedBlue < 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedBlue;
    }
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEUVErr is smaller than 0.
// Return -1 if gkBlindPixelQEUV    is 0.
// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2 
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUV() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedUV < 0.)
        return -1.;
      
      return gkBlindPixelQECoatedUV;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedUV < 0.)
        return -1.;
      
      
      return gkBlindPixelQEUnCoatedUV;
    }
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQECT1Err is smaller than 0.
// Return -1 if gkBlindPixelQECT1    is 0.
// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedCT1 < 0.)
        return -1.;
      
      return gkBlindPixelQECoatedCT1;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedCT1 < 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedCT1;
    }
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEGreenErr is smaller than 0.
// Return -1 if gkBlindPixelQEGreen    is 0.
// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreenRelVar() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedGreenErr < 0.)
        return -1.;
      
      if (gkBlindPixelQECoatedGreen  == 0.)
        return -1.;
      
      return gkBlindPixelQECoatedGreenErr * gkBlindPixelQECoatedGreenErr 
         / gkBlindPixelQECoatedGreen / gkBlindPixelQECoatedGreen ;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedGreenErr < 0.)
        return -1.;
      
      if (gkBlindPixelQEUnCoatedGreen  == 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedGreenErr * gkBlindPixelQEUnCoatedGreenErr 
         / gkBlindPixelQEUnCoatedGreen / gkBlindPixelQEUnCoatedGreen ;
    }
  
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEBlueErr is smaller than 0.
// Return -1 if gkBlindPixelQEBlue    is 0.
// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlueRelVar() const
{
  if (IsCoated())
    {
      if (gkBlindPixelQECoatedBlueErr < 0.)
        return -1.;
      
      if (gkBlindPixelQECoatedBlue  == 0.)
        return -1.;
      
      return gkBlindPixelQECoatedBlueErr * gkBlindPixelQECoatedBlueErr 
         / gkBlindPixelQECoatedBlue / gkBlindPixelQECoatedBlue ;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedBlueErr < 0.)
        return -1.;
      
      if (gkBlindPixelQEUnCoatedBlue  == 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedBlueErr * gkBlindPixelQEUnCoatedBlueErr 
         / gkBlindPixelQEUnCoatedBlue / gkBlindPixelQEUnCoatedBlue ;
    }
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQEUVErr is smaller than 0.
// Return -1 if gkBlindPixelQEUV    is 0.
// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2 
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUVRelVar() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedUVErr < 0.)
        return -1.;
      
      if (gkBlindPixelQECoatedUV  == 0.)
        return -1.;
      
      return gkBlindPixelQECoatedUVErr * gkBlindPixelQECoatedUVErr 
         / gkBlindPixelQECoatedUV / gkBlindPixelQECoatedUV ;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedUVErr < 0.)
        return -1.;
      
      if (gkBlindPixelQEUnCoatedUV  == 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedUVErr * gkBlindPixelQEUnCoatedUVErr 
         / gkBlindPixelQEUnCoatedUV / gkBlindPixelQEUnCoatedUV ;
    }
}

// --------------------------------------------------------------------------
//
// Return -1 if gkBlindPixelQECT1Err is smaller than 0.
// Return -1 if gkBlindPixelQECT1    is 0.
// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1RelVar() const
{

  if (IsCoated())
    {
      if (gkBlindPixelQECoatedCT1Err < 0.)
        return -1.;
      
      if (gkBlindPixelQECoatedCT1  == 0.)
        return -1.;
      
      return gkBlindPixelQECoatedCT1Err * gkBlindPixelQECoatedCT1Err 
         / gkBlindPixelQECoatedCT1 / gkBlindPixelQECoatedCT1 ;
    }
  else
    {
      if (gkBlindPixelQEUnCoatedCT1Err < 0.)
        return -1.;
      
      if (gkBlindPixelQEUnCoatedCT1  == 0.)
        return -1.;
      
      return gkBlindPixelQEUnCoatedCT1Err * gkBlindPixelQEUnCoatedCT1Err 
         / gkBlindPixelQEUnCoatedCT1 / gkBlindPixelQEUnCoatedCT1 ;
    }
}

// --------------------------------------------------------------------------
//
// Return gkBlindPixelCollectionEffErr^2 / (gkBlindPixelCollectionEff^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetBlindPixelCollectionEffRelVar() const
{

  return gkBlindPixelCollectionEffErr * gkBlindPixelCollectionEffErr 
    / gkBlindPixelCollectionEff / gkBlindPixelCollectionEff ;
}

// --------------------------------------------------------------------------
//
// Test bit kChargeFitValid
//
Bool_t MCalibrationChargeBlindPix::IsChargeFitValid()  const 
{
    return TESTBIT(fFlags,kChargeFitValid);
}

// --------------------------------------------------------------------------
//
// Test bit kCoated
//
Bool_t MCalibrationChargeBlindPix::IsCoated()  const 
{
    return TESTBIT(fFlags,kCoated);
}

// --------------------------------------------------------------------------
//
// Test bit kOscillating
//
Bool_t MCalibrationChargeBlindPix::IsOscillating()  const 
{
    return TESTBIT(fFlags,kOscillating);
}

// --------------------------------------------------------------------------
//
// Test bit kPedestalFitValid
//
Bool_t MCalibrationChargeBlindPix::IsPedestalFitOK()  const 
{
    return TESTBIT(fFlags,kPedestalFitOK);
}

// --------------------------------------------------------------------------
//
// Test bit kSinglePheFitValid
//
Bool_t MCalibrationChargeBlindPix::IsSinglePheFitOK()  const 
{
    return TESTBIT(fFlags,kSinglePheFitOK);
}

// --------------------------------------------------------------------------
//
// Test bit kFluxInsidePlexiglassAvailable
//
Bool_t  MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
}


// --------------------------------------------------------------------------
//
// Return kFALSE if IsChargeFitValid() is kFALSE
//
// Calculate fFluxInsidePlexiglass with the formula:
// - fFluxInsidePlexiglass    = fLambda 
//                            / GetBlindPixelCollectionEff()
//                            / GetBlindPixelQE() 
//                            * 10**gkBlindPixelAtt[color] 
//                            / gkBlindPixelArea 
// - fFluxInsidePlexiglassVar = sqrt( fLambdaVar / ( fLambda * fLambda ) 
//                             + GetBlindPixelQERelVar() 
//                             + GetBlindPixelCollectionEffRelVar()
//                               ) * fFluxInsidePlexiglass * * fFluxInsidePlexiglass
//
// If the fFluxInsidePlexiglass is smaller than 0., return kFALSE
// If the Variance is smaller than 0., return kFALSE
//
// SetFluxInsidePlexiglassAvailable() and return kTRUE
//
Bool_t MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
{

  if (IsChargeFitValid())
    return kFALSE;
  

  //
  // Start calculation of number of photons 
  // The blind pixel has exactly 100 mm^2 area (with negligible error), 
  //
  switch (fColor)
    {
    case MCalibrationCam::kGREEN:
      fFluxInsidePlexiglass    = fLambda / GetBlindPixelQEGreen() * TMath::Power(10,gkBlindPixelAttGreen);   
      // attenuation has negligible error
      fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEGreenRelVar();
      break;
    case MCalibrationCam::kBLUE:
      fFluxInsidePlexiglass    = fLambda / GetBlindPixelQEBlue() * TMath::Power(10,gkBlindPixelAttBlue);   
      // attenuation has negligible error
      fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEBlueRelVar();
      break;
    case MCalibrationCam::kUV:
      fFluxInsidePlexiglass    = fLambda / GetBlindPixelQEUV() * TMath::Power(10,gkBlindPixelAttUV);   
      // attenuation has negligible error
      fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEUVRelVar();
      break;
    case MCalibrationCam::kCT1:
    default:
      fFluxInsidePlexiglass    = fLambda / GetBlindPixelQECT1() * TMath::Power(10,gkBlindPixelAttCT1);   
      // attenuation has negligible error
      fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQECT1RelVar();
      break;
    }

  fFluxInsidePlexiglass    /= gkBlindPixelArea;
  fFluxInsidePlexiglass    /= gkBlindPixelCollectionEff;
  //
  // Finish calculation of errors -> convert from relative variance to absolute variance
  //
  fFluxInsidePlexiglassVar += GetBlindPixelCollectionEffRelVar();
  fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass;

  if (fFluxInsidePlexiglass < 0.)
      return kFALSE;

  if (fFluxInsidePlexiglassVar < 0.)
      return kFALSE;

  SetFluxInsidePlexiglassAvailable(kTRUE);  

  *fLog << inf << endl;
  *fLog << inf << GetDescriptor() 
        << ": Photon flux [ph/mm^2] inside Plexiglass: " 
        << Form("%5.3f%s%5.3f",fFluxInsidePlexiglass," +- ",GetFluxInsidePlexiglassErr()) << endl;

  return kTRUE;
}