/* ======================================================================== *\
!
! *
! * 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:
// - fArea: 100 mm^2
// - Average QE of Blind Pixel: 
//    fQEGreen: 0.154
//    fQEBlue : 0.226
//    fQEUV   : 0.247
//    fQECT1  : 0.247
// - Average QE Error of Blind Pixel: 
//    fQEGreenErr: 0.015;
//    fQEBlueErr : 0.02;
//    fQEUVErr   : 0.02;
//    fQECT1Err  : 0.02;
// - Attenuation factor Blind Pixel:
//    fAttGreen :  1.97;
//    fAttBlue  :  1.96;
//    fAttUV    :  1.95;
//    fAttCT1   :  1.95;
//
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeBlindPix.h"

#include <TH1.h>

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

ClassImp(MCalibrationChargeBlindPix);

using namespace std;
const Float_t MCalibrationChargeBlindPix::fgArea       = 100;
const Float_t MCalibrationChargeBlindPix::fgAttGreen   = 1.97;
const Float_t MCalibrationChargeBlindPix::fgAttBlue    = 1.96;
const Float_t MCalibrationChargeBlindPix::fgAttUV      = 1.95;
const Float_t MCalibrationChargeBlindPix::fgAttCT1     = 1.95;
const Float_t MCalibrationChargeBlindPix::fgAttErr     = 0.01;
const Float_t MCalibrationChargeBlindPix::fgQEGreen    = 0.154;
const Float_t MCalibrationChargeBlindPix::fgQEBlue     = 0.226;
const Float_t MCalibrationChargeBlindPix::fgQEUV       = 0.247;
const Float_t MCalibrationChargeBlindPix::fgQECT1      = 0.247;
const Float_t MCalibrationChargeBlindPix::fgQEErrGreen = 0.005;
const Float_t MCalibrationChargeBlindPix::fgQEErrBlue  = 0.007;
const Float_t MCalibrationChargeBlindPix::fgQEErrUV    = 0.01;
const Float_t MCalibrationChargeBlindPix::fgQEErrCT1   = 0.01;
const Float_t MCalibrationChargeBlindPix::fgCollEffGreen = 0.99; 
const Float_t MCalibrationChargeBlindPix::fgCollEffBlue  = 0.93; 
const Float_t MCalibrationChargeBlindPix::fgCollEffUV    = 0.90; 
const Float_t MCalibrationChargeBlindPix::fgCollEffCT1   = 0.90; 
const Float_t MCalibrationChargeBlindPix::fgCollEffErr   = 0.05; 
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Calls:
// - Clear()
//
// For backward-compatibility reasons, quantum eff., coll. eff. and att. 
// are intialized from the static members. This should, however, be 
// overwritten by a class deriving from MCalibrationChargeBlindCam. 
//
MCalibrationChargeBlindPix::MCalibrationChargeBlindPix(const char *name, const char *title)
{

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

  Clear();

  fArea = fgArea;

  fAtt       .Set( MCalibrationCam::gkNumPulserColors );       
  fAttErr    .Set( MCalibrationCam::gkNumPulserColors );    
  fQE        .Set( MCalibrationCam::gkNumPulserColors );        
  fQEErr     .Set( MCalibrationCam::gkNumPulserColors );     
  fCollEff   .Set( MCalibrationCam::gkNumPulserColors );     
  fCollEffErr.Set( MCalibrationCam::gkNumPulserColors );

  SetAtt        ( fgAttGreen,     MCalibrationCam::kGREEN );
  SetAtt        ( fgAttBlue,      MCalibrationCam::kBLUE  );
  SetAtt        ( fgAttUV   ,     MCalibrationCam::kUV    );
  SetAtt        ( fgAttCT1  ,     MCalibrationCam::kCT1   );

  SetAttErr     ( fgAttErr  ,     MCalibrationCam::kGREEN );
  SetAttErr     ( fgAttErr ,      MCalibrationCam::kBLUE  );
  SetAttErr     ( fgAttErr  ,     MCalibrationCam::kUV    );
  SetAttErr     ( fgAttErr  ,     MCalibrationCam::kCT1   );
  
  SetQE         ( fgQEGreen,      MCalibrationCam::kGREEN );
  SetQE         ( fgQEBlue,       MCalibrationCam::kBLUE  );
  SetQE         ( fgQEUV   ,      MCalibrationCam::kUV    );
  SetQE         ( fgQECT1  ,      MCalibrationCam::kCT1   );

  SetQEErr      ( fgQEErrGreen,   MCalibrationCam::kGREEN );
  SetQEErr      ( fgQEErrBlue,    MCalibrationCam::kBLUE  );
  SetQEErr      ( fgQEErrUV   ,   MCalibrationCam::kUV    );
  SetQEErr      ( fgQEErrCT1  ,   MCalibrationCam::kCT1   );

  SetCollEff    ( fgCollEffGreen, MCalibrationCam::kGREEN );
  SetCollEff    ( fgCollEffBlue,  MCalibrationCam::kBLUE  );
  SetCollEff    ( fgCollEffUV   , MCalibrationCam::kUV    );
  SetCollEff    ( fgCollEffCT1  , MCalibrationCam::kCT1   );

  SetCollEffErr ( fgCollEffErr,   MCalibrationCam::kGREEN );
  SetCollEffErr ( fgCollEffErr,   MCalibrationCam::kBLUE  );
  SetCollEffErr ( fgCollEffErr,   MCalibrationCam::kUV    );
  SetCollEffErr ( fgCollEffErr,   MCalibrationCam::kCT1   );

}


// ------------------------------------------------------------------------
//
// 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);

  MCalibrationPix::Clear();
}

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


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

// -----------------------------------------------------
//
// copy 'constructor'
//
void MCalibrationChargeBlindPix::Copy(TObject& object) const
{

  MCalibrationChargeBlindPix &pix = (MCalibrationChargeBlindPix&)object;

  //
  // Copy the data members
  //
  pix.fArea       = fArea;                     
  pix.fAtt        = fAtt;                      
  pix.fAttErr     = fAttErr;                   
  pix.fQE         = fQE;                       
  pix.fQEErr      = fQEErr;                    
  pix.fCollEff    = fCollEff;                  
  pix.fCollEffErr = fCollEffErr;               
 
  pix.fLambda                  = fLambda;                   
  pix.fLambdaCheck             = fLambdaCheck;              
  pix.fLambdaCheckErr          = fLambdaCheckErr;           
  pix.fLambdaVar               = fLambdaVar;                
  pix.fFluxInsidePlexiglass    = fFluxInsidePlexiglass;     
  pix.fFluxInsidePlexiglassVar = fFluxInsidePlexiglassVar;  
  pix.fMu0       = fMu0;                      
  pix.fMu0Err    = fMu0Err;                   
  pix.fMu1       = fMu1;                      
  pix.fMu1Err    = fMu1Err;                   
  pix.fSigma0    = fSigma0;                   
  pix.fSigma0Err = fSigma0Err;                
  pix.fSigma1    = fSigma1;                   
  pix.fSigma1Err = fSigma1Err;                
 
  pix.fColor     = fColor;
  pix.fPixId     = fPixId;
  pix.fFlags     = fFlags;
  
}


// --------------------------------------------------------------------------
//
// 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
// 
const 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
// 
const 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
// 
const 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 )
// 
const Float_t  MCalibrationChargeBlindPix::GetLambdaRelVar() const
{
  if (fLambdaVar < 0.)
    return -1.;
  
  if (fLambda  == 0.)
    return -1.;
  
  return fLambdaVar / fLambda / fLambda ;
}

// --------------------------------------------------------------------------
//
// Return TMath::Power(10,fAtt[fColor])
//
const Float_t MCalibrationChargeBlindPix::GetAtt() const
{
  return TMath::Power(10,fAtt[fColor]);   
}

// --------------------------------------------------------------------------
//
// Return -1 if fAttErr[fColor] is smaller than 0.
// Error of TMath::Power(10,fAtt[fColor]) = TMath::Power(10,fAtt[fColor])*ln(10.)*fAttErr[fColor]
// Return fAttErr^2 / (fAtt^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetAttRelVar() const
{

  const Float_t err = fAttErr[fColor];

  if (err < 0.)
    return -1.;
  
  return err*err*2.3;
}

// --------------------------------------------------------------------------
//
// Return fQE[fColor]
//
const Float_t MCalibrationChargeBlindPix::GetQE() const
{
  return fQE[fColor];
}

// --------------------------------------------------------------------------
//
// Return -1 if fQEErr[fColor] is smaller than 0.
// Return fQEErr^2 / (fQE^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetQERelVar() const
{

  if (fQEErr[fColor] < 0.)
    return -1.;
  
  return fQEErr[fColor]* fQEErr[fColor] / GetQE() / GetQE();
}

// --------------------------------------------------------------------------
//
// Return fCollEff[fColor]
//
const Float_t MCalibrationChargeBlindPix::GetCollEff() const
{
  return fCollEff[fColor];
}

// --------------------------------------------------------------------------
//
// Return -1 if fCollEffErr[fColor] is smaller than 0.
// Return fCollEffErr^2 / (fCollEff^2 )
// 
const Float_t MCalibrationChargeBlindPix::GetCollEffRelVar() const
{

  if (fCollEffErr[fColor] < 0.)
    return -1.;
  
  return fCollEffErr[fColor]* fCollEffErr[fColor] / GetCollEff() / GetCollEff();
}

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

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

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

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

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


// --------------------------------------------------------------------------
//
// Return kFALSE if IsChargeFitValid() is kFALSE
//
// Calculate fFluxInsidePlexiglass with the formula:
// - fFluxInsidePlexiglass    = fLambda 
//                            / GetCollEff()
//                            / GetQE() 
//                            * GetAtt()
//                            / fArea 
// - fFluxInsidePlexiglassVar = sqrt( fLambdaVar / ( fLambda * fLambda ) 
//                             + GetQERelVar() 
//                             + GetCollEffRelVar()
//                             + GetAttRelVar()
//                               ) * 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), 
  //
  fFluxInsidePlexiglass    = fLambda / GetQE() * GetAtt() / GetCollEff() / fArea;   

  if (fFluxInsidePlexiglass < 0.)
      return kFALSE;

  fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetQERelVar() + GetAttRelVar() + GetCollEffRelVar();

  //
  // Finish calculation of errors -> convert from relative variance to absolute variance
  //
  fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass;

  if (fFluxInsidePlexiglassVar < 0.)
      return kFALSE;

  SetFluxInsidePlexiglassAvailable(kTRUE);  

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

  return kTRUE;
}

void MCalibrationChargeBlindPix::Print(Option_t *opt) const
{
  
  *fLog << all << GetDescriptor() 
        << Form("%s%3i","BlindPixel: ",GetPixId())
        << Form("%s%4.2f%s%4.2f","  Lambda: ",GetLambda(),"+-",GetLambdaErr())
        << Form("%s%4.2f%s%4.2f","  Mu0: ",GetMu0(),"+-",GetMu0Err())
        << Form("%s%4.2f%s%4.2f","  Mu1: ",GetMu1(),"+-",GetMu1Err()) 
        << Form("%s%4.2f%s%4.2f","  Sigma0: ",GetSigma0(),"+-",GetSigma0Err())
        << Form("%s%4.2f%s%4.2f","  Sigma1: ",GetSigma1(),"+-",GetSigma1Err())
        << endl;
  *fLog << all
        << " Pedestal Fit OK? :" << IsPedestalFitOK() 
        << Form("%s%4.2f%s%4.2f","  Lambda (Check): " ,GetLambdaCheck(),"+-",GetLambdaCheckErr()) 
        << endl;
  *fLog << all
        << " Flux available? :" << IsFluxInsidePlexiglassAvailable() 
        << Form("%s%4.2f%s%4.2f","  Flux: " ,GetFluxInsidePlexiglass(),"+-",GetFluxInsidePlexiglassErr())
        << endl;
}









