source: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc@ 2662

/* ======================================================================== *\
!
! *
! * 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 Calibration information for all pixels in the camera           //
//                                                                         //
// The Classes MCalibrationPix's are stored in a TClonesArray              //
// The Class MCalibrationBlindPix and the MCalibrationPINDiode             // 
//           are stored in separate pointers                               //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationCam.h"
#include "MCalibrationPix.h"
#include "MCalibrationBlindPix.h"
#include "MCalibrationConfig.h"

#include <TClonesArray.h>

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

#include "TCanvas.h"

#include "MGeomCam.h"

ClassImp(MCalibrationCam);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
// Creates a TClonesArray of MHCalibrationPix objects, initialized to 0 entries
// Creates an MCalibrationBlindPix and an MCalibrationPINDiode
//
MCalibrationCam::MCalibrationCam(const char *name, const char *title)
    : fMeanNrPhotAvailable(kFALSE),
      fMeanNrPhotInnerPix(-1.),
      fMeanNrPhotInnerPixErr(-1.),
      fOffsets(NULL),
      fSlopes(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();
}

// --------------------------------------------------------------------------
//
// Delete the TClonesArray or MCalibrationPix's
// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
//
MCalibrationCam::~MCalibrationCam()
{

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

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

}

void MCalibrationCam::DrawHiLoFits()
{

  fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-300.,300.);
  fSlopes  = new TH1D("mm","Slopess of the HiGain LoGain Fit",100,-3.,3.);

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {

      fOffsets->Fill(pix->GetHist()->GetOffset());
      fSlopes->Fill(pix->GetHist()->GetSlope());

    }

  TCanvas *c1 = new TCanvas();
  c1->Divide(1,2);
  c1->cd(1);
  fOffsets->Draw();
  gPad->Modified();
  gPad->Update();

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


// --------------------------------------------------------------------------
//
// This function return the size of the FILLED MCalibrationCam
// It is NOT the size of the array fPixels !!!
// Note that with every call to AddPixels, GetSize() might change
//
Int_t MCalibrationCam::GetSize() const
{

  //
  // Here we get the number of "filled" fPixels!!
  //
  return fPixels->GetEntriesFast();
}

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

// --------------------------------------------------------------------------
//
// Return a pointer to the pixel with the requested idx. 
// NULL if it doesn't exist. 
//
MCalibrationPix *MCalibrationCam::GetCalibrationPix(Int_t idx) const
{
   if (idx<0)
     return NULL;

  if (!CheckBounds(idx))
    return NULL;
  
  return (MCalibrationPix*)fPixels->At(idx);
}



// --------------------------------------------------------------------------
//
// Check if position i is inside bounds
//
Bool_t MCalibrationCam::CheckBounds(Int_t i) const 
{
    return i < fPixels->GetEntriesFast();
} 

Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const 
{
  if (!CheckBounds(idx))
    return kFALSE;

  return kTRUE;
}

Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const 
{
  return ((*this)[idx].GetCharge() > 0. && (*this)[idx].GetErrCharge() > 0.);
}



void MCalibrationCam::Clear(Option_t *o)
{
    fPixels->ForEach(TObject, Clear)();
}

//
// Perform the fits of the charges 
// If i=-1, then all pixels will be fitted
// Otherwise only the one with index i
// 
// The number of succesful fits is returned
//
UShort_t MCalibrationCam::FitCharge(Int_t i)
{

  UShort_t nsuccess = 0;

  // invalidate i if it exceeds the number of entries in fPixels
  if (i > fPixels->GetEntriesFast())      
    {
      *fLog << warn << "Tried to fit pixel out of allowed range " << endl;
      return 0;
    }
  
  if (i == -1)    // loop over all events
    {

      TIter Next(fPixels);
      MCalibrationPix *pix;
      while ((pix=(MCalibrationPix*)Next()))
        {
          if (pix->FitCharge())
            nsuccess++;
        }
    }
  else                  // fit only the pixel with index i
    {
      if((*this)[i].FitCharge())
        nsuccess++;
    }
  
  return nsuccess;
  
}

//
// Perform the fits of the charges of all pixels 
// plus the blind pixel and the PIN Diode
// 
// The number of succesful fits is returned
//
UShort_t MCalibrationCam::FitAllCharge()
{

  // FIXME: Once the blind pixel is fully working, 
  //        there must be some penalty in case the
  //        fit does not succeed
  //
  UShort_t nsuccess = 0;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      if (pix->FitCharge())
        nsuccess++;
    }
  
  if (fBlindPixel->FitCharge())
        nsuccess++;

  if (fPINDiode->FitCharge())
        nsuccess++;

  return nsuccess;

}



//
// Perform the fits of the arrival times
// If i=-1, then all pixels will be fitted
// Otherwise only the one with index i
// 
// The number of succesful fits is returned
//
UShort_t MCalibrationCam::FitTime(Int_t i)
{

  UShort_t nsuccess = 0;

  // invalidate i if it exceeds the number of entries in fPixels
  if (i > fPixels->GetEntriesFast())      
    {
      *fLog << warn << "Tried to fit pixel out of allowed range " << endl;
      return 0;
    }
  
  if (i == -1)                // loop over all events
    {

      TIter Next(fPixels);
      MCalibrationPix *pix;
      while ((pix=(MCalibrationPix*)Next()))
        {
         if (pix->FitTime())
            nsuccess++;
        }
    }
  else                     // fit only the pixel with index i
    {
      if((*this)[i].FitTime())
        nsuccess++;
    }
  
  return nsuccess;

}


//
// Perform the fits of the times of all pixels 
// plus the blind pixel and the PIN Diode
// 
// The number of succesful fits is returned
//
UShort_t MCalibrationCam::FitAllTime()
{

  // FIXME: Once the blind pixel is fully working, 
  //        there must be some penalty in case the
  //        fit does not succeed
  //
  UShort_t nsuccess = 0;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      if (pix->FitTime())
        nsuccess++;
    }
  
  if (fBlindPixel->FitTime())
        nsuccess++;

  if (fPINDiode->FitTime())
       nsuccess++;

  return nsuccess;

}



void MCalibrationCam::CutEdges()
{

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

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

  return;
}

// ---------------------------------------------------------------------
//
// 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: 0)
//
void MCalibrationCam::InitSize(Int_t size)
{
  
  //
  // check if we have already initialized to size
  //
  if (CheckBounds(size))
    return;
  
  //
  // It is important to use Expand and NOT ExpandCreate, because 
  // we want to keep all pixel not used with a NULL pointer.
  //
  fPixels->ExpandCreate(size);

}
  
void MCalibrationCam::Print(Option_t *o) const
{
    *fLog << all << GetDescriptor() << ":" << endl;
    int id = 0;

    *fLog << "Succesfully calibrated pixels:" << endl;
    *fLog << endl;

    TIter Next(fPixels);
    MCalibrationPix *pix;
    while ((pix=(MCalibrationPix*)Next()))
    {

      if (pix->GetCharge() >= 0.)
        {
          *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() 
                << " Reduced Charge: " << pix->GetCharge() << " +- " 
                << pix->GetSigmaCharge() << " Reduced Sigma: " << TMath::Sqrt(pix->GetRSigmaSquare()) << endl;
          id++;
        }
    }

    *fLog << id << " succesful pixels :-))" << endl;
    id = 0;

    *fLog << endl;
    *fLog << "Pixels with errors:" << endl;
    *fLog << endl;

    TIter Next2(fPixels);
    while ((pix=(MCalibrationPix*)Next2()))
    {

      if (pix->GetCharge() == -1.)
        {
          *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() 
                << " Reduced Charge: " << pix->GetCharge() << " +- " 
                << pix->GetSigmaCharge() << " Reduced Sigma: " << TMath::Sqrt(pix->GetRSigmaSquare()) << endl;
          id++;
        }
    }
    *fLog << id << " pixels with errors :-((" << endl;

}


Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{

  if (idx > GetSize())
    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:
      val = (*this)[idx].GetTime();
      break;
    case 6:
      val = (*this)[idx].GetSigmaTime();
      break;
    case 7:
      val = (*this)[idx].GetTimeChiSquare();
      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 15:
      val = ((*this)[idx].GetSigmaCharge()/(*this)[idx].GetCharge())*
            ((*this)[idx].GetSigmaCharge()/(*this)[idx].GetCharge());
      break;
    case 11:
      val = (*this)[idx].GetPheFFactorMethod();
      break;
    case 12:
      val = (*this)[idx].GetConversionFFactorMethod();
      break;
    case 13:
      if (idx < 397)
        val = (double)fMeanNrPhotInnerPix;
      else
        val = (double)fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea;
      break;
    case 14:
      if ((fMeanNrPhotInnerPix > 0. ) && ((*this)[idx].GetCharge() != -1.))
        {
          if (idx < 397)
            val = fMeanNrPhotInnerPix / (*this)[idx].GetCharge();
          else 
            val = fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea / (*this)[idx].GetCharge();
        } 
      else
        {
          val = -1.;
        }
      break;
    default:
      return kFALSE;
    }
  return val>=0;
}

void MCalibrationCam::DrawPixelContent(Int_t idx) const
{

  (*this)[idx].Draw();

}

Bool_t MCalibrationCam::CalcNrPhotInnerPixel()
{

  if (!fBlindPixel->IsFitOK())
    return kFALSE;
  
  const Float_t mean = fBlindPixel->GetLambda();
  //  const Float_t merr = fBlindPixel->GetErrLambda();
  
  switch (fColor)
    {
    case kECGreen:
      fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEGreen) // real photons
                            *TMath::Power(10,gkCalibrationBlindPixelAttGreen) // correct for absorption 
                            / gkCalibrationBlindPixelArea;          // correct for area
      break;
    case kECBlue:
      fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEBlue )
                            *TMath::Power(10,gkCalibrationBlindPixelAttBlue)
                            / gkCalibrationBlindPixelArea;
      break;
    case kECUV:
      fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEUV )
                            *TMath::Power(10,gkCalibrationBlindPixelAttUV)
                            / gkCalibrationBlindPixelArea;
      break;
    case kECCT1:
    default:
      fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQECT1 )
                            *TMath::Power(10,gkCalibrationBlindPixelAttCT1)
                            / gkCalibrationBlindPixelArea;
      break;
    }

  fMeanNrPhotAvailable = kTRUE;

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      
      if((pix->GetCharge() > 0.) && (fMeanNrPhotInnerPix > 0.))
        pix->SetConversionBlindPixelMethod(fMeanNrPhotInnerPix/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 (!fMeanNrPhotAvailable)
    if (!CalcNrPhotInnerPixel())
      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].GetConversionFFactorMethod();

  if (conv < 0.)
    return kFALSE;

  mean = conv;

  //
  // Not yet ready , sorry 
  //
  err  = -1.;
  sigma = -1.;

  return kTRUE;
}


Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  return kTRUE;

}

Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  return kTRUE;

}