source: trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc@ 4336

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MHCalibrationCam                                               
//
// Base class for camera calibration classes. Incorporates the TObjArray's:
// - fHiGainArray (for calibrated High Gains per pixel)
// - fLoGainArray (for calibrated Low Gains per pixel) 
// - fAverageHiGainAreas (for averaged High Gains events per camera area index)
// - fAverageLoGainAreas (for averaged High Gains events per camera area index)
// - fAverageHiGainSectors (for averaged High Gains events per camera sector )
// - fAverageLoGainSectors (for averaged High Gains events per camera sector )
// These TObjArray's are called by their default constructors, thus no objects 
// are created, until the derived class does so. 
//
// The corresponding operators: [],() and the operators GetAverageHiGainArea(), 
// GetAverageLoGainArea(), GetAverageHiGainSector() and GetAverageLoGainSector() 
// have to be cast to the corresponding class. It is assumed that all classes 
// dealing with calibration pixels derive from MHGausEvents.
//
/////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationCam.h"

#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TText.h>
#include <TPaveText.h>

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

#include "MCalibrationPix.h"
#include "MCalibrationCam.h"

#include "MHGausEvents.h"

#include "MBadPixelsPix.h"
#include "MBadPixelsCam.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MParList.h"

#include "MRawRunHeader.h"

ClassImp(MHCalibrationCam);

using namespace std;

const Int_t   MHCalibrationCam::fgAverageNbins    = 2000;
const Int_t   MHCalibrationCam::fgPulserFrequency = 500;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets:
// - all pointers to NULL
//
// Initializes and sets owner of:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
// Initializes:
// - fPulserFrequency to fgPulserFrequency
//
MHCalibrationCam::MHCalibrationCam(const char *name, const char *title)
    :  fBadPixels(NULL), fCam(NULL), fGeom(NULL), fRunHeader(NULL), fDebug(kFALSE)
{

    fHiGainArray = new TObjArray;
    fHiGainArray->SetOwner();
    
    fLoGainArray = new TObjArray;
    fLoGainArray->SetOwner();

    fAverageHiGainAreas = new TObjArray;
    fAverageHiGainAreas->SetOwner();

    fAverageLoGainAreas = new TObjArray;
    fAverageLoGainAreas->SetOwner();

    fAverageHiGainSectors = new TObjArray;
    fAverageHiGainSectors->SetOwner();

    fAverageLoGainSectors = new TObjArray;
    fAverageLoGainSectors->SetOwner();

    SetAverageNbins();
    SetPulserFrequency();

    fHiGainOverflow = 0;
    fLoGainOverflow = 0;
}

// --------------------------------------------------------------------------
//
// Deletes the TClonesArray of:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
MHCalibrationCam::~MHCalibrationCam()
{

  delete fHiGainArray;
  delete fLoGainArray;

  delete fAverageHiGainAreas;
  delete fAverageLoGainAreas;

  delete fAverageHiGainSectors;
  delete fAverageLoGainSectors;
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel (pixel number)
//
MHGausEvents &MHCalibrationCam::operator[](UInt_t i)
{
  return *static_cast<MHGausEvents*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel (pixel number)
//
const MHGausEvents &MHCalibrationCam::operator[](UInt_t i) const
{
  return *static_cast<MHGausEvents*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel (pixel number)
//
MHGausEvents  &MHCalibrationCam::operator()(UInt_t i)
{
  return *static_cast<MHGausEvents*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel (pixel number)
//
const MHGausEvents  &MHCalibrationCam::operator()(UInt_t i) const
{
  return *static_cast<MHGausEvents*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Returns the current size of the TObjArray fAverageHiGainAreas
// independently if the MHGausEvents is filled with values or not.
//
const Int_t MHCalibrationCam::GetAverageAreas() const
{
  return fAverageHiGainAreas->GetEntries();
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel Area (area number)
//
MHGausEvents  &MHCalibrationCam::GetAverageHiGainArea(UInt_t i)
{
  return *static_cast<MHGausEvents*>(fAverageHiGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel Area (area number)
//
const MHGausEvents  &MHCalibrationCam::GetAverageHiGainArea(UInt_t i) const
{
  return *static_cast<MHGausEvents *>(fAverageHiGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel Area (area number)
//
MHGausEvents  &MHCalibrationCam::GetAverageLoGainArea(UInt_t i)
{
  return *static_cast<MHGausEvents*>(fAverageLoGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel Area (area number)
//
const MHGausEvents  &MHCalibrationCam::GetAverageLoGainArea(UInt_t i) const
{
  return *static_cast<MHGausEvents*>(fAverageLoGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Returns the current size of the TObjArray fAverageHiGainSectors
// independently if the MHGausEvents is filled with values or not.
//
const Int_t MHCalibrationCam::GetAverageSectors() const
{
  return fAverageHiGainSectors->GetEntries();
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain Sector (sector number)
//
MHGausEvents  &MHCalibrationCam::GetAverageHiGainSector(UInt_t i)
{
  return *static_cast<MHGausEvents*>(fAverageHiGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain Sector (sector number)
//
const MHGausEvents  &MHCalibrationCam::GetAverageHiGainSector(UInt_t i) const
{
  return *static_cast<MHGausEvents*>(fAverageHiGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain Sector (sector number)
//
MHGausEvents  &MHCalibrationCam::GetAverageLoGainSector(UInt_t i)
{
  return *static_cast<MHGausEvents*>(fAverageLoGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain Sector (sector number)
//
const MHGausEvents  &MHCalibrationCam::GetAverageLoGainSector(UInt_t i) const
{
  return *static_cast<MHGausEvents*>(fAverageLoGainSectors->UncheckedAt(i));
}


const TArrayI &MHCalibrationCam::GetRunNumbers() const 
{
  return fRunNumbers;
}

// --------------------------------------------------------------------------
//
// Our own clone function is necessary since root 3.01/06 or Mars 0.4
// I don't know the reason. 
//
// Creates new MHCalibrationCam
// Deletes the TObjArray's and Clones them individually
// Copies the TArray's
// Copies the fPulserFrequency
//
TObject *MHCalibrationCam::Clone(const char *) const
{

  //  const Int_t nhi   = fHiGainArray->GetEntries();
  //  const Int_t nlo   = fLoGainArray->GetEntries();
  const Int_t navhi = fAverageHiGainAreas->GetEntries();
  const Int_t navlo = fAverageLoGainAreas->GetEntries();
  const Int_t nsehi = fAverageHiGainSectors->GetEntries();
  const Int_t nselo = fAverageLoGainSectors->GetEntries();
  
  //
  // FIXME, this might be done faster and more elegant, by direct copy.
  //
  MHCalibrationCam *cam = new MHCalibrationCam();

  //  cam->fHiGainArray->Expand(nhi);
  //  cam->fLoGainArray->Expand(nlo);
  cam->fAverageHiGainAreas->Expand(navhi);
  cam->fAverageLoGainAreas->Expand(navlo);
  cam->fAverageHiGainSectors->Expand(nsehi);
  cam->fAverageLoGainSectors->Expand(nselo);

  /*
  for (int i=0; i<nhi; i++)
    {
      delete (*cam->fHiGainArray)[i];
      (*cam->fHiGainArray)[i] = (*fHiGainArray)[i]->Clone();
    }
  for (int i=0; i<nlo; i++)
    {
      delete (*cam->fLoGainArray)[i];
      (*cam->fLoGainArray)[i] = (*fLoGainArray)[i]->Clone();
    }
  */
  
  for (int i=0; i<navhi; i++)
    {
      //      delete (*cam->fAverageHiGainAreas)[i];
      (*cam->fAverageHiGainAreas)[i] = (*fAverageHiGainAreas)[i]->Clone();
    }
  for (int i=0; i<navlo; i++)
    {
      //      delete (*cam->fAverageLoGainAreas)[i];
      (*cam->fAverageLoGainAreas)[i] = (*fAverageLoGainAreas)[i]->Clone();
    }
  for (int i=0; i<nsehi; i++)
    {
      //      delete (*cam->fAverageHiGainSectors)[i];
      (*cam->fAverageHiGainSectors)[i] = (*fAverageHiGainSectors)[i]->Clone();
    }
  for (int i=0; i<nselo; i++)
    {
      //      delete (*cam->fAverageLoGainSectors)[i];
      (*cam->fAverageLoGainSectors)[i] = (*fAverageLoGainSectors)[i]->Clone();
    }

  cam->fAverageAreaNum         = fAverageAreaNum;
  cam->fAverageAreaSat         = fAverageAreaSat;
  cam->fAverageAreaSigma       = fAverageAreaSigma;      
  cam->fAverageAreaSigmaVar    = fAverageAreaSigmaVar;   
  cam->fAverageAreaRelSigma    = fAverageAreaRelSigma;
  cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar;   
  cam->fAverageSectorNum       = fAverageSectorNum;      
  cam->fRunNumbers             = fRunNumbers;

  cam->fPulserFrequency        = fPulserFrequency;
  cam->fAverageNbins           = fAverageNbins;

  return cam;
}

// --------------------------------------------------------------------------
//
// Gets the pointers to:
// - MGeomCam
//
// Calls SetupHists(const MParList *pList)
//
// Calls Delete-Function of:
// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
//
Bool_t MHCalibrationCam::SetupFill(const MParList *pList)
{
  
  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeom)
  {
      *fLog << err << GetDescriptor() 
            << ": MGeomCam not found... aborting." << endl;
      return kFALSE;
  }

  fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
  if (!fRunHeader)
  {
    *fLog << warn << GetDescriptor() 
          << ": MRawRunHeader not found... will not store run numbers." << endl;
  }

  fHiGainArray->Delete();
  fLoGainArray->Delete();

  fAverageHiGainAreas->Delete();
  fAverageLoGainAreas->Delete();

  fAverageHiGainSectors->Delete();
  fAverageLoGainSectors->Delete();

  return SetupHists(pList);
}


Bool_t MHCalibrationCam::SetupHists(const MParList *pList)
{
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Gets or creates the pointers to:
// - MBadPixelsCam
//
// Searches pointer to:
// - MArrivalTimeCam
//
// Initializes, if empty to MArrivalTimeCam::GetSize() for:
// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
//
// Initializes, if empty to MGeomCam::GetNumAreas() for:
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
//
// Initializes, if empty to MGeomCam::GetNumSectors() for:
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
// 
// Initializes TArray's to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively
// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()):
// - MHCalibrationCam::fAverageAreaNum[area index]
// - MHCalibrationCam::fAverageSectorNum[area index]
//
// Calls InitializeHists() for every entry in:
// - MHCalibrationCam::fHiGainArray
// - MHCalibrationCam::fAverageHiGainAreas
// - MHCalibrationCam::fAverageHiGainSectors
//
// Sets Titles and Names for the Histograms 
// - MHCalibrationCam::fAverageHiGainAreas
// - MHCalibrationCam::fAverageHiGainSectors
// 
// Retrieves the run numbers from MRawRunHeader and stores them in fRunNumbers
//
Bool_t MHCalibrationCam::ReInit(MParList *pList)
{

  const Int_t npixels  = fGeom->GetNumPixels();
  const Int_t nsectors = fGeom->GetNumSectors();
  const Int_t nareas   = fGeom->GetNumAreas();

  fBadPixels = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam");
  if (!fBadPixels)
    {

      fBadPixels = (MBadPixelsCam*)pList->FindCreateObj(AddSerialNumber("MBadPixelsCam"));
      if (!fBadPixels)
        {
          gLog << err << "Cannot find nor create MBadPixelsCam ... abort." << endl;
          return kFALSE;
        }
      else 
        fBadPixels->InitSize(npixels);
    }

  //
  // The function TArrayF::Set() already sets all entries to 0.
  //
  fAverageAreaNum.        Set(nareas);
  fAverageAreaSat.        Set(nareas);           
  fAverageAreaSigma.      Set(nareas);      
  fAverageAreaSigmaVar.   Set(nareas);   
  fAverageAreaRelSigma.   Set(nareas);   
  fAverageAreaRelSigmaVar.Set(nareas);
  fAverageSectorNum.      Set(nsectors);
  fRunNumbers.            Set(fRunNumbers.GetSize()+1);

  for (Int_t aidx=0; aidx<nareas; aidx++)
      fAverageAreaNum[aidx] = 0;

  for (Int_t sector=0; sector<nsectors; sector++)
      fAverageSectorNum[sector] = 0;

  for (Int_t i=0; i<npixels; i++)
    {

      if ((*fBadPixels)[i].IsBad())
        continue;

      fAverageAreaNum  [(*fGeom)[i].GetAidx()  ]++;
      fAverageSectorNum[(*fGeom)[i].GetSector()]++;
    }

  if (fRunHeader)
    fRunNumbers[fRunNumbers.GetSize()-1] = fRunHeader->GetRunNumber();

  if (!ReInitHists(pList))
    return kFALSE;

  if (!fRunHeader)
    return kTRUE;
  
  for (Int_t i=0; i<fHiGainArray->GetEntries(); i++)
    {
      TH1F *h = (*this)[i].GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }

  for (Int_t i=0; i<fLoGainArray->GetEntries(); i++)
    {
      TH1F *h = (*this)(i).GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }
  
  for (Int_t j=0; j<nareas; j++)
    {
      TH1F *h = GetAverageHiGainArea(j).GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }
  
  for (Int_t j=0; j<nareas; j++)
    {
      TH1F *h = GetAverageLoGainArea(j).GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }
  
  for (Int_t j=0; j<nsectors; j++)
    {
      TH1F *h = GetAverageHiGainSector(j).GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }
  
  for (Int_t j=0; j<nsectors; j++)
    {
      TH1F *h = GetAverageLoGainSector(j).GetHGausHist();
      h->SetTitle( Form("%s%i%s", h->GetTitle(),fRunNumbers[fRunNumbers.GetSize()-1]," "));
    }

  return kTRUE;
}



Bool_t MHCalibrationCam::ReInitHists(MParList *pList)
{
  return kTRUE;
}



//--------------------------------------------------------------------------------
//
// Retrieves from MGeomCam:
// - number of pixels
// - number of pixel areas
// - number of sectors
//
// For all TObjArray's (including the averaged ones), the following steps are performed: 
//
// 1) Test size and return kFALSE if not matching
// 2) 
//
Bool_t MHCalibrationCam::Fill(const MParContainer *par, const Stat_t w)
{

  const Int_t npixels  = fGeom->GetNumPixels();
  const Int_t nareas   = fGeom->GetNumAreas();
  const Int_t nsectors = fGeom->GetNumSectors();
  
  if (fHiGainArray->GetEntries() != npixels)
    {
      gLog << err << "ERROR - Size mismatch... abort." << endl;
      return kFALSE;
    }
  
  if (fLoGainArray->GetEntries() != npixels)
    {
      gLog << err << "ERROR - Size mismatch... abort." << endl;
      return kFALSE;
    }
  
  if (fAverageHiGainAreas->GetEntries() != nareas)
    {
      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
      return kFALSE;
    }

  if (fAverageLoGainAreas->GetEntries() != nareas)
    {
      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
      return kFALSE;
    }

  if (fAverageHiGainSectors->GetEntries() != nsectors)
    {
      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
      return kFALSE;
    }

  if (fAverageLoGainSectors->GetEntries() != nsectors)
    {
      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
      return kFALSE;
    }

  return FillHists(par, w);
}

Bool_t MHCalibrationCam::FillHists(const MParContainer *par, const Stat_t w)
{
  *fLog << warn << GetDescriptor() << "FillHists not overloaded! Can't be used!" << endl;
  return kFALSE;
}

// --------------------------------------------------------------------------
//
// 1) FinalizeHists()
// 2) FinalizeBadPixels()
// 3) CalcAverageSigma()
//
Bool_t MHCalibrationCam::Finalize()
{
  if (!FinalizeHists())
    return kFALSE;

  FinalizeBadPixels();
  CalcAverageSigma();

  return kTRUE;
}

Bool_t MHCalibrationCam::FinalizeHists()
{
  return kTRUE;
}

void MHCalibrationCam::FinalizeBadPixels()
{
}


// -------------------------------------------------------------
//
// If MBadPixelsPix::IsBad():
// - calls MHGausEvents::SetExcluded()
//
// Calls:
// - MHGausEvents::InitBins()
// - MHGausEvents::ChangeHistId(i)
// - MHGausEvents::SetEventFrequency(fPulserFrequency)
// 
void MHCalibrationCam::InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i)
{

  if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
    hist.SetExcluded();
  
  hist.InitBins();
  hist.ChangeHistId(i);
  hist.SetEventFrequency(fPulserFrequency);

  TH1F *h = hist.GetHGausHist();
  h->SetTitle( Form("%s%s", h->GetTitle()," Runs: "));
}

void MHCalibrationCam::FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
                                       MBadPixelsPix::UncalibratedType_t fittyp, 
                                       MBadPixelsPix::UncalibratedType_t osctyp)
{
  
  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
    {
      
      MHGausEvents &hist = (*this)[i];
      
      if (hist.IsExcluded())
        continue;
      
      MCalibrationPix &pix    = calcam[i];
      MBadPixelsPix   &bad    = badcam[i];
      
      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
      
    }

  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
    {
      
      MHGausEvents     &hist = GetAverageHiGainArea(j);      
      MCalibrationPix  &pix  = calcam.GetAverageArea(j);
      MBadPixelsPix    &bad  = calcam.GetAverageBadArea(j);        
      
      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
  }
  

  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
    {
      
      MHGausEvents     &hist = GetAverageHiGainSector(j);      
      MCalibrationPix  &pix  = calcam.GetAverageSector(j);
      MBadPixelsPix    &bad  = calcam.GetAverageBadSector(j);        
      
      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
    }

}

void MHCalibrationCam::FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
                                            MBadPixelsPix::UncalibratedType_t fittyp, 
                                            MBadPixelsPix::UncalibratedType_t osctyp)
{
  
  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
    {
      
      MHGausEvents &hist = (*this)(i);
      
      if (hist.IsExcluded())
        continue;
      
      MCalibrationPix &pix    = calcam[i];
      MBadPixelsPix   &bad    = badcam[i];
      
      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
      
    }

  for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
    {
      
      MHGausEvents     &hist = GetAverageLoGainArea(j);      
      MCalibrationPix  &pix  = calcam.GetAverageArea(j);
      MBadPixelsPix    &bad  = calcam.GetAverageBadArea(j);        
      
      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
  }
  

  for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
    {
      
      MHGausEvents     &hist = GetAverageLoGainSector(j);      
      MCalibrationPix  &pix  = calcam.GetAverageSector(j);
      MBadPixelsPix    &bad  = calcam.GetAverageBadSector(j);        
      
      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
    }
}

//------------------------------------------------------------
//
// For all averaged areas, the fitted sigma is multiplied with the square root of 
// the number involved pixels
//
void MHCalibrationCam::CalcAverageSigma()
{
  
  for (UInt_t j=0; j<fGeom->GetNumAreas(); j++)
    {
  
      MCalibrationPix &pix    = fCam->GetAverageArea(j);

      const Float_t numsqr    = TMath::Sqrt((Float_t)fAverageAreaNum[j]);
      fAverageAreaSigma[j]    = pix.GetSigma    () * numsqr;
      fAverageAreaSigmaVar[j] = pix.GetSigmaErr () * pix.GetSigmaErr() * numsqr;

      pix.SetSigma   (fAverageAreaSigma[j]);
      pix.SetSigmaVar(fAverageAreaSigmaVar[j]);

      fAverageAreaRelSigma   [j]  = fAverageAreaSigma[j]    / pix.GetMean();
      fAverageAreaRelSigmaVar[j]  = fAverageAreaSigmaVar[j] / (fAverageAreaSigma[j]*fAverageAreaSigma[j]);
      fAverageAreaRelSigmaVar[j] += pix.GetMeanRelVar();
      fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j];
    }
}

// ---------------------------------------------------------------------------
//
// Returns if the histogram is empty and sets the following flag:
// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
//
// Fits the histograms with a Gaussian, in case of failure 
// calls MHGausEvents::RepeatFit(), in case of repeated failure 
// calls MHGausEvents::BypassFit() and sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
// 
// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 
// In case no, sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
//
// Retrieves the results and store them in MCalibrationPix
//
void MHCalibrationCam::FitHiGainHists(MHGausEvents &hist, 
                                      MCalibrationPix &pix, 
                                      MBadPixelsPix &bad, 
                                      MBadPixelsPix::UncalibratedType_t fittyp,
                                      MBadPixelsPix::UncalibratedType_t osctyp)
{


  if (hist.IsEmpty())
    return;

  //
  // 2) Fit the Hi Gain histograms with a Gaussian
  //
  if (!hist.FitGaus())
    //
    // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
    //
    if (!hist.RepeatFit())
      {
        hist.BypassFit();
        bad.SetUncalibrated( fittyp );
      }
  
  hist.Renorm();
  //
  // 4) Check for oscillations
  // 
  hist.CreateFourierSpectrum();
  

  if (!hist.IsFourierSpectrumOK())
    bad.SetUncalibrated( osctyp );
  
  //
  // 5) Retrieve the results and store them in this class
  //
  pix.SetHiGainMean       ( hist.GetMean()      );
  pix.SetHiGainMeanVar    ( hist.GetMeanErr() * hist.GetMeanErr()   );
  pix.SetHiGainSigma      ( hist.GetSigma()     );
  pix.SetHiGainSigmaVar   ( hist.GetSigmaErr()* hist.GetSigmaErr()  );
  pix.SetHiGainProb       ( hist.GetProb()      );
  pix.SetHiGainNumBlackout( hist.GetBlackout()  );
  pix.SetHiGainNumPickup  ( hist.GetPickup()    );
  
  if (IsDebug())
    {
      *fLog << dbginf << GetDescriptor() << ": ID " << hist.GetPixId() 
            << " HiGainSaturation: " << pix.IsHiGainSaturation() 
            << " HiGainMean: " << hist.GetMean()
            << " HiGainMeanErr: " << hist.GetMeanErr()
            << " HiGainMeanSigma: " << hist.GetSigma()
            << " HiGainMeanSigmaErr: " << hist.GetSigmaErr()
            << " HiGainMeanProb: " << hist.GetProb()
            << " HiGainNumBlackout: " << hist.GetBlackout()
            << " HiGainNumPickup  : " << hist.GetPickup  ()
            << endl;
    }

}


// ---------------------------------------------------------------------------
//
// Returns if the histogram is empty and sets the following flag:
// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
//
// Fits the histograms with a Gaussian, in case of failure 
// calls MHGausEvents::RepeatFit(), in case of repeated failure 
// calls MHGausEvents::BypassFit() and sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
// 
// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 
// In case no, sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
//
// Retrieves the results and store them in MCalibrationPix
//
void MHCalibrationCam::FitLoGainHists(MHGausEvents &hist, 
                                      MCalibrationPix &pix, 
                                      MBadPixelsPix &bad, 
                                      MBadPixelsPix::UncalibratedType_t fittyp,
                                      MBadPixelsPix::UncalibratedType_t osctyp)
{

  if (hist.IsEmpty())
      return;
  

  //
  // 2) Fit the Hi Gain histograms with a Gaussian
  //
  if (!hist.FitGaus())
    //
    // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
    //
    if (!hist.RepeatFit())
      {
        hist.BypassFit();
        bad.SetUncalibrated( fittyp );
      }
  
  //
  // 4) Check for oscillations
  // 
  hist.CreateFourierSpectrum();
  
  if (!hist.IsFourierSpectrumOK())
    bad.SetUncalibrated( osctyp );
  
  //
  // 5) Retrieve the results and store them in this class
  //
  pix.SetLoGainMean       ( hist.GetMean()      );
  pix.SetLoGainMeanVar    ( hist.GetMeanErr()  * hist.GetMeanErr()   );
  pix.SetLoGainSigma      ( hist.GetSigma()     );
  pix.SetLoGainSigmaVar   ( hist.GetSigmaErr() * hist.GetSigmaErr()  );
  pix.SetLoGainProb       ( hist.GetProb()      );
  pix.SetLoGainNumBlackout( hist.GetBlackout()  );
  pix.SetLoGainNumPickup  ( hist.GetPickup()    );
  
  if (IsDebug())
    {
      *fLog << dbginf << GetDescriptor() << "ID: " << hist.GetPixId() 
            << " HiGainSaturation: " << pix.IsHiGainSaturation() 
            << " LoGainMean: " << hist.GetMean()
            << " LoGainMeanErr: " << hist.GetMeanErr()
            << " LoGainMeanSigma: " << hist.GetSigma()
            << " LoGainMeanSigmaErr: " << hist.GetSigmaErr()
            << " LoGainMeanProb: " << hist.GetProb()
            << " LoGainNumBlackout: " << hist.GetBlackout()
            << " LoGainNumPickup  : " << hist.GetPickup  ()
            << endl;
    }

}



// --------------------------------------------------------------------------
//
// Dummy, needed by MCamEvent
//
Bool_t MHCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// What MHCamera needs in order to draw an individual pixel in the camera
//
void MHCalibrationCam::DrawPixelContent(Int_t idx) const
{
}

// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
//
void MHCalibrationCam::Draw(const Option_t *opt)
{

  const Int_t nareas = fAverageHiGainAreas->GetEntries();
  if (nareas == 0)
    return;

  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);  
  pad->SetBorderMode(0);

  pad->Divide(2,nareas);

  for (Int_t i=0; i<nareas;i++) 
    {

      pad->cd(2*(i+1)-1);
      GetAverageHiGainArea(i).Draw(opt);

      if (!fAverageAreaSat[i])
        DrawAverageSigma(fAverageAreaSat[i], i,
                         fAverageAreaSigma[i],    fAverageAreaSigmaVar[i],
                         fAverageAreaRelSigma[i], fAverageAreaRelSigmaVar[i]);

      pad->cd(2*(i+1));
      GetAverageLoGainArea(i).Draw(opt);
      
      if (fAverageAreaSat[i])
        DrawAverageSigma(fAverageAreaSat[i], i,
                         fAverageAreaSigma[i], fAverageAreaSigmaVar[i],
                         fAverageAreaRelSigma[i], fAverageAreaRelSigmaVar[i]);
    }

}

// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays a TPaveText with the re-normalized sigmas of the average area
//
void MHCalibrationCam::DrawAverageSigma(Bool_t sat, Bool_t inner,
                                              Float_t sigma, Float_t sigmavar,
                                              Float_t relsigma, Float_t relsigmavar) const 
{
  
  if (sigma != 0 && sigmavar >= 0 && relsigmavar >= 0.)
    {
      
      TPad *newpad = new TPad("newpad","transparent",0,0,1,1);
      newpad->SetFillStyle(4000);
      newpad->Draw();
      newpad->cd();
      
      TPaveText *text = new TPaveText(sat? 0.1 : 0.35,0.7,sat ? 0.4 : 0.7,1.0);
      text->SetTextSize(0.07);
      const TString line1 = Form("%s%s%s",inner ? "Outer" : "Inner",
                                 " Pixels ", sat ? "Low Gain" : "High Gain");
      TText *txt1 = text->AddText(line1.Data());
      const TString line2 = Form("#sigma per pix: %2.2f #pm %2.2f",sigma,TMath::Sqrt(sigmavar));
      TText *txt2 = text->AddText(line2.Data());
      const TString line3 = Form("Rel. #sigma per pix: %2.2f #pm %2.2f",relsigma,TMath::Sqrt(relsigmavar));
      TText *txt3 = text->AddText(line3.Data());
      text->Draw("");
      
      text->SetBit(kCanDelete);
      txt1->SetBit(kCanDelete);
      txt2->SetBit(kCanDelete);
      txt3->SetBit(kCanDelete);
      newpad->SetBit(kCanDelete);
    }
}