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

/* ======================================================================== *\
!
! *
! * 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 "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.)
{
    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;
}

// --------------------------------------------------------------------------
//
// 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].GetRCharge() > 0. && (*this)[idx].GetErrRCharge() > 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;

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

        *fLog << id << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() << " Charge: " 
              << pix->GetCharge() << " Reduced Charge: " << pix->GetRCharge() << " +- " 
              << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigma() << endl;

        id++;
    }
}


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

  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].GetTimeProb();
        break;
    case 8:
       val = (*this)[idx].GetPed();
        break;
    case 9:
       val = (*this)[idx].GetPedRms();
        break;
    case 10:
      val = (*this)[idx].GetRCharge();
      break;
    case 11:
      val = (*this)[idx].GetRSigma();
      break;
    case 12:
      val = (*this)[idx].GetPheFFactorMethod();
      break;
    case 13:
      val = (*this)[idx].GetConversionFFactorMethod();
      break;
    case 14:
      val = (double)fMeanNrPhotInnerPix;
      break;
    case 15:
      if ((fMeanNrPhotInnerPix > 0. ) && ((*this)[idx].GetRCharge() > 100.))
        val = fMeanNrPhotInnerPix / (*this)[idx].GetRCharge();
      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->IsValid())
    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;
    }

  fMeanNrPhotAvailable = kTRUE;
  return kTRUE;
}



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

  if (!fMeanNrPhotAvailable)
    if (!CalcNrPhotInnerPixel())
      return kFALSE;

  mean = fMeanNrPhotInnerPix / (*this)[ipx].GetRCharge();

  //
  // Not yet ready , sorry 
  //
  err  = 100000000000.;

  return kTRUE;
}