/* ======================================================================== *\
!
! *
! * 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                                                    //
//                                                                         //
// This is the storage container to hold informations about the calibration//
// blind pixel                                                             //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeBlindPix.h"

#include <TH1.h>

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

ClassImp(MCalibrationChargeBlindPix);

using namespace std;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelArea = 100;
// Average QE of Blind Pixel (three colours)
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEGreen = 0.154;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEBlue  = 0.226;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUV    = 0.247;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECT1   = 0.247;
// Average QE Error of Blind Pixel (three colours)
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEGreenErr = 0.015;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEBlueErr  = 0.02;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUVErr    = 0.02;
const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECT1Err   = 0.02;
// Attenuation factor Blind Pixel (three colours)
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::fgLambdaCheckLimit   = 0.2;
const Float_t MCalibrationChargeBlindPix::fgLambdaErrLimit     = 0.2;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MCalibrationChargeBlindPix::MCalibrationChargeBlindPix(const char *name, const char *title)
{

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

  SetLambdaCheckLimit();
  SetLambdaErrLimit();

  Clear();
}


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

  fLambda      = -1.;
  fLambdaCheck = -1.;
  fMu0         = -1.;
  fMu1         = -1.;
  fSigma0      = -1.;
  fSigma1      = -1.;
  fLambdaErr   = -1.;
  fMu0Err      = -1.;
  fMu1Err      = -1.;
  fSigma0Err   = -1.;
  fSigma1Err   = -1.;

  fMeanFluxInsidePlexiglass          = -1.;
  fMeanFluxErrInsidePlexiglass       = -1.;

  SetOscillating                   ( kFALSE );
  SetExcluded                      ( kFALSE );
  SetChargeFitValid                ( kFALSE );
  SetPedestalFitOK                 ( kFALSE );
  SetSinglePheFitOK                ( kFALSE );
  SetFluxInsidePlexiglassAvailable ( kFALSE );

}

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

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

void  MCalibrationChargeBlindPix::SetChargeFitValid( const Bool_t b)
{
    b ? SETBIT(fFlags,kChargeFitValid) : CLRBIT(fFlags,kChargeFitValid);
}
void  MCalibrationChargeBlindPix::SetPedestalFitOK( const Bool_t b)
{
    b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
}

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

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

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

Bool_t MCalibrationChargeBlindPix::IsOscillating()  const 
{
    return TESTBIT(fFlags,kOscillating);
}

Bool_t MCalibrationChargeBlindPix::IsChargeFitValid()  const 
{
    return TESTBIT(fFlags,kChargeFitValid);
}

Bool_t MCalibrationChargeBlindPix::IsPedestalFitOK()  const 
{
    return TESTBIT(fFlags,kPedestalFitOK);
}

Bool_t MCalibrationChargeBlindPix::IsSinglePheFitOK()  const 
{
    return TESTBIT(fFlags,kSinglePheFitOK);
}

Bool_t  MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
}


//
// The check return kTRUE if:
//
// 1) fLambda and fLambdaCheck are separated relatively by fLambdaCheckLimit
// 2) BlindPixel has an fLambdaErr smaller than  fLambdaErrLimit
// 
Bool_t MCalibrationChargeBlindPix::CheckChargeFitValidity()
{

  if (2.*(fLambdaCheck-fLambda)/(fLambdaCheck+fLambda) < fLambdaCheckLimit)
    {
      *fLog << warn << "WARNING: Lambda and Lambda-Check differ by more than "
            << fLambdaCheckLimit << " in the Blind Pixel " << endl;
      return kFALSE;
    }
  
  if (fLambdaErr < fLambdaErrLimit) 
    {
      *fLog << warn << "WARNING: Error of Fitted Lambda is greater than "
            << fLambdaErrLimit << " in Blind Pixel " << endl;
      return kFALSE;
    }
      
  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), 
  //
  fMeanFluxInsidePlexiglass    = fLambda*gkBlindPixelArea;

  // Start calculation of number of photons relative Variance 
  fMeanFluxErrInsidePlexiglass  = fLambdaErr*fLambdaErr/fLambda/fLambda;
  
  switch (fColor)
    {
    case kGREEN:
      fMeanFluxInsidePlexiglass    /= gkBlindPixelQEGreen;   
      fMeanFluxErrInsidePlexiglass += gkBlindPixelQEGreenErr*gkBlindPixelQEGreenErr
                                    / gkBlindPixelQEGreen  /   gkBlindPixelQEGreen;   

      fMeanFluxInsidePlexiglass  *= TMath::Power(10,gkBlindPixelAttGreen); // correct for absorption 
      // attenuation has negligible error
      break;
    case kBLUE:
      fMeanFluxInsidePlexiglass    /= gkBlindPixelQEBlue;   
      fMeanFluxErrInsidePlexiglass += gkBlindPixelQEBlueErr*gkBlindPixelQEBlueErr
                                    / gkBlindPixelQEBlue  /   gkBlindPixelQEBlue;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttBlue); // correct for absorption 
      // attenuation has negligible error
      break;
    case kUV:
      fMeanFluxInsidePlexiglass    /= gkBlindPixelQEUV;   
      fMeanFluxErrInsidePlexiglass += gkBlindPixelQEUVErr*gkBlindPixelQEUVErr
                                    / gkBlindPixelQEUV  /   gkBlindPixelQEUV;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttUV); // correct for absorption 
      // attenuation has negligible error
      break;
    case kCT1:
    default:
      fMeanFluxInsidePlexiglass    /= gkBlindPixelQECT1;   
      fMeanFluxErrInsidePlexiglass += gkBlindPixelQECT1Err*gkBlindPixelQECT1Err
                                    / gkBlindPixelQECT1  /   gkBlindPixelQECT1;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkBlindPixelAttCT1); // correct for absorption 
      // attenuation has negligible error
      break;
    }

  *fLog << inf << endl;
  *fLog << inf << " Photon flux [ph/mm^2] inside Plexiglass: " 
        << fMeanFluxInsidePlexiglass << endl;

  if (fMeanFluxInsidePlexiglass < 0.)
      return kFALSE;

  if (fMeanFluxErrInsidePlexiglass < 0.)
      return kFALSE;

  SetFluxInsidePlexiglassAvailable(kTRUE);  

  // Finish calculation of errors -> convert from relative variance to absolute error
  fMeanFluxErrInsidePlexiglass = TMath::Sqrt(fMeanFluxErrInsidePlexiglass);
  fMeanFluxErrInsidePlexiglass *= fMeanFluxInsidePlexiglass;

  *fLog << inf << " Error on photon flux [ph/mm^2] inside Plexiglass: " 
        << fMeanFluxErrInsidePlexiglass << endl;
  *fLog << inf << endl;

  return kTRUE;
}