source: trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc@ 3027

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

/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationCam                                               
//                                                               
// Hold the whole Calibration results of the camera:
//                                                               
// 1) MCalibrationCam initializes a TClonesArray whose elements are 
//    pointers to MCalibrationPix Containers
// 2) It initializes a pointer to an MCalibrationBlindPix container
// 3) It initializes a pointer to an MCalibrationPINDiode container
//
// 4)  
// 
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationCam.h"

#include <TH2.h>
#include <TCanvas.h>
#include <TClonesArray.h>

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

#include "MGeomCam.h"

#include "MCalibrationPix.h"
#include "MCalibrationConfig.h"
#include "MCalibrationBlindPix.h"
#include "MCalibrationPINDiode.h"

#include "MHCalibrationPixel.h"

ClassImp(MCalibrationCam);

using namespace std;

const Int_t   MCalibrationCam::gkBlindPixelId   =  559;
const Int_t   MCalibrationCam::gkPINDiodeId     = 9999;
const Float_t MCalibrationCam::gkTimeSliceWidth = 3.3;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
// Later, a call to MCalibrationCam::InitSize(Int_t size) has to be performed
//
// Creates an MCalibrationBlindPix container 
// Creates an MCalibrationPINDiode container
//
MCalibrationCam::MCalibrationCam(const char *name, const char *title)
    : fOffsets(NULL),
      fSlopes(NULL),
      fOffvsSlope(NULL)
{
    fName  = name  ? name  : "MCalibrationCam";
    fTitle = title ? title : "Storage container for the Calibration Information in the camera";

    fPixels     = new TClonesArray("MCalibrationPix",1);
    fBlindPixel = new MCalibrationBlindPix();
    fPINDiode   = new MCalibrationPINDiode();

    Clear();
}

// --------------------------------------------------------------------------
//
// Delete the TClonesArray of MCalibrationPix containers
// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
//
// Delete the histograms if they exist
//
MCalibrationCam::~MCalibrationCam()
{

  //
  // delete fPixels should delete all Objects stored inside
  // 
  delete fPixels;
  delete fBlindPixel;
  delete fPINDiode;

  if (fOffsets)
    delete fOffsets;
  if (fSlopes)
    delete fSlopes;
  if (fOffvsSlope)
    delete fOffvsSlope;

}

// -------------------------------------------------------------------
//
// This function simply allocates memory via the ROOT command:
// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*),
//                                      fSize * sizeof(TObject*));
// newSize corresponds to size in our case
// fSize is the old size (in most cases: 1)
//
void MCalibrationCam::InitSize(const UInt_t i)
{
  
  //
  // check if we have already initialized to size
  //
  if (CheckBounds(i))
    return;
  
  fPixels->ExpandCreate(i);

}

// --------------------------------------------------------------------------
//
// This function returns the current size of the TClonesArray 
// independently if the MCalibrationPix is filled with values or not.
//
// It is the size of the array fPixels.
//
Int_t MCalibrationCam::GetSize() const
{
  return fPixels->GetEntriesFast();
}

// --------------------------------------------------------------------------
//
// Check if position i is inside the current bounds of the TClonesArray
//
Bool_t MCalibrationCam::CheckBounds(Int_t i) const 
{
  return i < GetSize();
} 


// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MCalibrationPix &MCalibrationCam::operator[](Int_t i)
{
  return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MCalibrationPix &MCalibrationCam::operator[](Int_t i) const
{
  return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
}


// --------------------------------------
//
void MCalibrationCam::Clear(Option_t *o)
{

  fPixels->ForEach(TObject, Clear)();
  fBlindPixel->Clear();
  fPINDiode->Clear();

  fMeanPhotInsidePlexiglass          = -1.;
  fMeanPhotErrInsidePlexiglass       = -1.;
  fMeanPhotOutsidePlexiglass         = -1.;
  fMeanPhotErrOutsidePlexiglass      = -1.;

  fNumExcludedPixels                 = 0;

  CLRBIT(fFlags,kBlindPixelMethodValid);
  CLRBIT(fFlags,kPINDiodeMethodValid);
  CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
  CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);  

  return;
}

void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b)
{

  if (b) 
    SETBIT(fFlags, kBlindPixelMethodValid); 
  else    
    CLRBIT(fFlags, kBlindPixelMethodValid); 
  
}

void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b)
{

  if (b) 
    SETBIT(fFlags, kPINDiodeMethodValid); 
  else    
    CLRBIT(fFlags, kPINDiodeMethodValid); 
  
  
}

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

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


Bool_t  MCalibrationCam::IsNumPhotInsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
}

Bool_t  MCalibrationCam::IsNumPhotOutsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
}



// --------------------------------------------------------------------------
//
// Print first the well fitted pixels 
// and then the ones which are not FitValid
//
void MCalibrationCam::Print(Option_t *o) const
{

  *fLog << all << GetDescriptor() << ":" << endl;
  int id = 0;
  
  *fLog << all << "Succesfully calibrated pixels:" << endl;
  *fLog << all << endl;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      
      if (pix->IsChargeFitValid() && !pix->IsExcluded()) 
        {

          Float_t rsigma = pix->GetRSigmaSquare();
          if (rsigma > 0.)
            rsigma = TMath::Sqrt(rsigma);

          *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " 
                << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " 
                << pix->GetSigmaCharge() << " Reduced Sigma: " << rsigma 
                << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl;
          id++;
        }
    }
  
  *fLog << all << id << " succesful pixels :-))" << endl;
  id = 0;
  
  *fLog << all << endl;
  *fLog << all << "Pixels with errors:" << endl;
  *fLog << all << endl;
  
  TIter Next2(fPixels);
    while ((pix=(MCalibrationPix*)Next2()))
      {
        
        if (!pix->IsChargeFitValid() && !pix->IsExcluded())
          {

            Float_t rsigma = pix->GetRSigmaSquare();
            if (rsigma > 0.)
              rsigma = TMath::Sqrt(rsigma);
            
            *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " 
                  << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " 
                  << pix->GetSigmaCharge() << " Reduced Sigma: " << rsigma << endl;
            id++;
          }
      }
    *fLog << all << id << " pixels with errors :-((" << endl;
    
  *fLog << all << endl;
  *fLog << all << "Excluded pixels:" << endl;
  *fLog << all << endl;
  
  TIter Next3(fPixels);
    while ((pix=(MCalibrationPix*)Next3()))
      if (pix->IsExcluded())
        *fLog << all << pix->GetPixId() << endl;

  *fLog << all << fNumExcludedPixels << " excluded pixels " << endl;
}

// --------------------------------------------------------------------------
//
// Return true if pixel is inside bounds of the TClonesArray fPixels
//
Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const 
{
  if (!CheckBounds(idx))
    return kFALSE;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Return true if pixel has already been fitted once (independent of the result)
//
Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const 
{

  if (!CheckBounds(idx))
    return kFALSE;

  return (*this)[idx].IsFitted();
}

// --------------------------------------------------------------------------
//
// Sets the user ranges of all histograms such that 
// empty bins at the edges are not used. Additionally, it rebins the 
// histograms such that in total, 50 bins are used.
//
void MCalibrationCam::CutEdges()
{

  fBlindPixel->GetHist()->CutAllEdges();
  fPINDiode->GetHist()->CutAllEdges();

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      pix->GetHist()->CutAllEdges();
    }

  return;
}
  

// The types are as follows:
// 
// 0: Fitted Charge
// 1: Error of fitted Charge
// 2: Sigma of fitted Charge
// 3: Error of Sigma of fitted Charge
// 4: Returned probability of Gauss fit to Charge distribution
// 5: Mean arrival time
// 6: Sigma of the arrival time
// 7: Chi-square of the Gauss fit to the arrival times
// 8: Pedestal
// 9: Pedestal RMS
// 10: Reduced Sigma Square
// 11: Number of Photo-electrons after the F-Factor method
// 12: Error on the Number of Photo-electrons after the F-Factor method
// 13: Mean conversion factor after the F-Factor method
// 14: Error on the conversion factor after the F-Factor method
// 15: Number of Photons after the Blind Pixel method
// 16: Mean conversion factor after the Blind Pixel method
//
Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{

  if (idx > GetSize())
    return kFALSE;

  if ( (!(*this)[idx].IsChargeFitValid()) || (*this)[idx].IsExcluded())
    return kFALSE;
  
  if (idx == gkBlindPixelId)
    return kFALSE;

  if (idx == gkPINDiodeId)
    return kFALSE;

  switch (type)
    {
    case 0:
      val = (*this)[idx].GetCharge();
      break;
    case 1:
      val = (*this)[idx].GetErrCharge();
      break;
    case 2:
      val = (*this)[idx].GetSigmaCharge();
      break;
    case 3:
      val = (*this)[idx].GetErrSigmaCharge();
      break;
    case 4:
      val = (*this)[idx].GetChargeProb();
      break;
    case 5:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetTime() * gkTimeSliceWidth;
      break;
    case 6:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetSigmaTime() * gkTimeSliceWidth;
      break;
    case 7:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetTimeProb();
      break;
    case 8:
      val = (*this)[idx].GetPed();
      break;
    case 9:
      val = (*this)[idx].GetPedRms();
      break;
    case 10:
      if ((*this)[idx].GetRSigmaSquare() > 0.)
        val = TMath::Sqrt((*this)[idx].GetRSigmaSquare());
      else
        val = -1.;
      break;
    case 11:
      val = (*this)[idx].GetPheFFactorMethod();
      break;
    case 12:
      val = (*this)[idx].GetPheFFactorMethodError();
      break;
    case 13:
      val = (*this)[idx].GetMeanConversionFFactorMethod();
      break;
    case 14:
      val = (*this)[idx].GetErrorConversionFFactorMethod();
      break;
    case 15:
      if (idx < 397)
        val = (double)fMeanPhotInsidePlexiglass;
      else
        val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
      break;
    case 16:
      if (idx < 397)
        val = (*this)[idx].GetMeanConversionBlindPixelMethod();
      else 
        val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea;
      break;
    case 17:
      if ( (*this)[idx].GetRSigmaSquare() > 0. && (*this)[idx].GetCharge() > 0. )
        val = TMath::Sqrt((*this)[idx].GetRSigmaSquare()) / (*this)[idx].GetCharge();
      else
        val = -1.;
      break;
    case 18:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeMean();
      break;
    case 19:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeMeanErr();
      break;
    case 20:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeRms();
      break;
    case 21:
      if (!(*this)[idx].IsTimeFitValid())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.);
      break;
    default:
      return kFALSE;
    }
  return val!=-1.;
}

// --------------------------------------------------------------------------
//
// What MHCamera needs in order to draw an individual pixel in the camera
//
void MCalibrationCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].Draw();
}


// --------------------------------------------------------------------------
//
//
//
Bool_t MCalibrationCam::CalcNumPhotInsidePlexiglass()
{

  if (!fBlindPixel->IsFitOK())
    return kFALSE;
  
  const Float_t mean = fBlindPixel->GetLambda();
  const Float_t merr = fBlindPixel->GetErrLambda();
  
  switch (fColor)
    {
    case kECGreen:
      fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEGreen)     // real photons
                            *TMath::Power(10,gkCalibrationBlindPixelAttGreen) // correct for absorption 
                            * gkCalibrationInnerPixelArea;                    // correct for area

      
      break;
    case kECBlue:
      fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEBlue )
                            *TMath::Power(10,gkCalibrationBlindPixelAttBlue)
                            * gkCalibrationInnerPixelArea;
      break;
    case kECUV:
      fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEUV )
                            *TMath::Power(10,gkCalibrationBlindPixelAttUV)
                            * gkCalibrationInnerPixelArea;
      break;
    case kECCT1:
    default:
      fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQECT1 )
                            *TMath::Power(10,gkCalibrationBlindPixelAttCT1)
                            * gkCalibrationInnerPixelArea;
      break;
    }

  SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable);

  *fLog << inf << endl;
  *fLog << inf << "Mean number of Photons for an Inner Pixel (inside Plexiglass): " 
        << fMeanPhotInsidePlexiglass << endl;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.))
        {

          Float_t conversion      = fMeanPhotInsidePlexiglass/pix->GetCharge();
          Float_t conversionerr   = 0.;
          Float_t conversionsigma = 0.;
          pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);

          if (conversionerr/conversion < 0.1) 
            pix->SetBlindPixelMethodValid();
        }
    }
  return kTRUE;
}


Bool_t MCalibrationCam::CalcNumPhotOutsidePlexiglass()
{

  if (!fPINDiode->IsChargeFitValid())
    return kFALSE;
  
  const Float_t mean = fPINDiode->GetCharge();
  const Float_t merr = fPINDiode->GetErrCharge();
  
  switch (fColor)
    {
    case kECGreen:
      fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEGreen)   // real photons
                            * gkCalibrationInnerPixelvsPINDiodeArea;        // correct for area
      break;
    case kECBlue:
      fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEBlue )
                            * gkCalibrationInnerPixelvsPINDiodeArea;
      break;
    case kECUV:
      fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEUV )
                            * gkCalibrationInnerPixelvsPINDiodeArea;
      break;
    case kECCT1:
    default:
      fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQECT1 )
                            * gkCalibrationInnerPixelvsPINDiodeArea;
      break;
    }

  SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);  

  *fLog << inf << endl;
  *fLog << inf << mean << " Mean number of Photons for an Inner Pixel (outside Plexiglass): " 
        << fMeanPhotOutsidePlexiglass << endl;
  *fLog << inf << endl;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      
      if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.))
        pix->SetConversionPINDiodeMethod(fMeanPhotOutsidePlexiglass/pix->GetCharge(), 0., 0.);
    }
  return kTRUE;
}



Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{
  
  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  if (!IsNumPhotInsidePlexiglassAvailable())
    if (!CalcNumPhotInsidePlexiglass())
      return kFALSE;

  mean  = (*this)[ipx].GetMeanConversionBlindPixelMethod();
  err   = (*this)[ipx].GetErrorConversionBlindPixelMethod();
  sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod();

  return kTRUE;
}


Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{
  
  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod();

  if (conv < 0.)
    return kFALSE;

  mean  = conv;
  err   = (*this)[ipx].GetErrorConversionFFactorMethod();
  sigma = (*this)[ipx].GetSigmaConversionFFactorMethod();

  return kTRUE;
}


//-----------------------------------------------------------------------------------
//
// Calculates the conversion factor between the integral of FADCs slices 
// (as defined in the signal extractor MExtractSignal.cc)
// and the number of photons reaching the plexiglass for one Inner Pixel 
//
// FIXME: The PINDiode is still not working and so is the code 
//
Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  if (!IsNumPhotOutsidePlexiglassAvailable())
    if (!CalcNumPhotOutsidePlexiglass())
      return kFALSE;

  mean  = (*this)[ipx].GetMeanConversionPINDiodeMethod();
  err   = (*this)[ipx].GetErrorConversionPINDiodeMethod();
  sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod();

  return kFALSE;

}

//-----------------------------------------------------------------------------------
//
// Calculates the best combination of the three used methods possible 
// between the integral of FADCs slices 
// (as defined in the signal extractor MExtractSignal.cc)
// and the number of photons reaching one Inner Pixel. 
// The procedure is not yet defined.
//
// FIXME: The PINDiode is still not working and so is the code 
//
Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  return kFALSE;

}


void MCalibrationCam::DrawHiLoFits()
{

  if (!fOffsets)
    fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.);
  if (!fSlopes)
    fSlopes  = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.);
  if (!fOffvsSlope)
    fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.);

  TIter Next(fPixels);
  MCalibrationPix *pix;
  MHCalibrationPixel *hist;
  while ((pix=(MCalibrationPix*)Next()))
    {
      hist = pix->GetHist();
      hist->FitHiGainvsLoGain();
      fOffsets->Fill(hist->GetOffset(),1.);
      fSlopes->Fill(hist->GetSlope(),1.);
      fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.);
    }

   TCanvas *c1 = new TCanvas();

   c1->Divide(1,3);
   c1->cd(1);
   fOffsets->Draw();
   gPad->Modified();
   gPad->Update();

   c1->cd(2);
  fSlopes->Draw();
  gPad->Modified();
  gPad->Update();

  c1->cd(3);
  fOffvsSlope->Draw("col1");
  gPad->Modified();
  gPad->Update();
}