source: trunk/MagicSoft/Mars/manalysis/MHPedestalCam.cc@ 4200

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz  12/2000 <mailto:tbretz@atsro.uni-wuerzburh.de>
!   Author(s): Markus Gaug   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHPedestalCam
//
// Fills the extracted pedestals of MExtractedSignalCam into
// the MHGausEvents-classes MHPedestalPix for every:
//
// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
//   or MHCalibrationCam::fHiGainArray, respectively, depending if
//   MArrivalTimePix::IsLoGainUsed() is set.
//
// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
//   MHCalibrationCam::fAverageHiGainAreas
//
// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
//   and MHCalibrationCam::fAverageHiGainSectors
//
// Every pedestal is directly taken from MExtractedSignalCam, filled by the
// appropriate extractor.
//
// The pedestals are filled into a histogram and an array, in order to perform
// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
// event-by-event basis and written into the corresponding average pixels.
//
// The histograms are fitted to a Gaussian, mean and sigma with its errors
// and the fit probability are extracted. If none of these values are NaN's and
// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
// the fit is declared valid.
// Otherwise, the fit is repeated within ranges of the previous mean
// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit())
// In case this does not make the fit valid, the histogram means and RMS's are
// taken directly (see MHGausEvents::BypassFit()).
//
// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas
// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup)
//
// The number of summed FADC slices is taken to re-normalize
// the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()):
// - Mean Pedestal        / slice = Mean Pedestal        / Number slices
// - Mean Pedestal Error  / slice = Mean Pedestal Error  / Number slices
// - Sigma Pedestal       / slice = Sigma Pedestal       / Sqrt (Number slices)
// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices)
//
// The class also fills arrays with the signal vs. event number, creates a fourier
// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
// projected fourier components follow an exponential distribution.
//
// This same procedure is performed for the average pixels.
//
// The following results are written into MPedestalCam:
//
// - MCalibrationPix::SetMean()
// - MCalibrationPix::SetMeanErr()
// - MCalibrationPix::SetSigma()
// - MCalibrationPix::SetSigmaErr()
// - MCalibrationPix::SetProb()
// - MCalibrationPix::SetNumPickup()
//
// For all averaged areas, the fitted sigma is multiplied with the square root of
// the number involved pixels in order to be able to compare it to the average of
// sigmas in the camera.
//
/////////////////////////////////////////////////////////////////////////////
#include "MHPedestalCam.h"
#include "MHPedestalPix.h"

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

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

#include "MParList.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

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

#include "MCalibrationPedCam.h"

ClassImp(MHPedestalCam);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Initializes:
// - fExtractHiGainSlices to 0.
// - fExtractLoGainSlices to 0.
//
MHPedestalCam::MHPedestalCam(const char *name, const char *title) 
    : fExtractHiGainSlices(0.), fExtractLoGainSlices(0.)
{
    fName  = name  ? name  : "MHPedestalCam";
    fTitle = title ? title : "";

}



// --------------------------------------------------------------------------
//
// Searches pointer to:
// - MPedestalCam
// - MExtractedSignalCam
//
// Searches or creates:
// - MCalibrationPedCam
//
// Retrieves from MExtractedSignalCam:
// - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
// - number of used Low  Gain FADC slices (MExtractedSignalCam::GetNumUsedLoGainFADCSlices())
//
// Initializes, if empty to MGeomCam::GetNumPixels():
// - 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
// 
// Calls MHCalibrationCam::InitPedHists() for every entry in:
// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
//
// Sets Titles and Names for the Histograms 
// - MHCalibrationCam::fAverageHiGainAreas
// - MHCalibrationCam::fAverageHiGainSectors
// 
// Sets number of bins to MHCalibrationCam::fAverageNbins for:
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
//
Bool_t MHPedestalCam::ReInitHists(MParList *pList)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
  if (!signal)
    {
      gLog << err << "Cannot find MExtractedSignalCam... abort." << endl;
      return kFALSE;
    }
  

  fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");

  if (!fPedestals)
    {
      gLog << err << "Cannot find MPedestalCam ... abort." << endl;
      return kFALSE;
    }

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

  fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam");
  if (!fCam)
  {
      fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationPedCam"));
      if (!fCam)
          return kFALSE;

      fCam->Init(*fGeom);
  }
  
  Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
  Float_t sliceslo = signal->GetNumUsedLoGainFADCSlices();

  if (sliceshi == 0.)
    {
      gLog << err << "Number of used signal slices in MExtractedSignalCam is zero  ... abort." 
           << endl;
      return kFALSE;
    }

  if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
    {
      gLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam  ... abort." 
           << endl;
      return kFALSE;
    }

  if (fExtractLoGainSlices != 0. && sliceslo != fExtractLoGainSlices )
    {
      gLog << err << "Number of used Low Gain signal slices changed in MExtractedSignalCam  ... abort." 
           << endl;
      return kFALSE;
    }

  fExtractHiGainSlices = sliceshi;
  fExtractLoGainSlices = sliceslo;

  if (fHiGainArray->GetEntries()==0)
  {
      fHiGainArray->Expand(npixels);
      for (Int_t i=0; i<npixels; i++)
      {
          (*fHiGainArray)[i] = new MHPedestalPix;
          InitPedHists((MHPedestalPix&)(*this)[i],i,fExtractHiGainSlices);
          
          if ((*fBadPixels)[i].IsBad())
            {
              (*this)[i].SetExcluded();
            }
          (*fPedestals)[i].InitUseHists();
      }

  }

  if (fLoGainArray->GetEntries()==0)
  {
      fLoGainArray->Expand(npixels);
      for (Int_t i=0; i<npixels; i++)
      {
          (*fLoGainArray)[i] = new MHPedestalPix;
          InitPedHists((MHPedestalPix&)(*this)(i),i,fExtractLoGainSlices);

          if ((*fBadPixels)[i].IsBad())
            (*this)(i).SetExcluded();
      }
  }

  if (fAverageHiGainAreas->GetEntries()==0)
  {
    fAverageHiGainAreas->Expand(nareas);
    
    for (Int_t j=0; j<nareas; j++)
      {
        (*fAverageHiGainAreas)[j] = 
          new MHPedestalPix("AverageHiGainArea",
                                      "Average Pedestals area idx ");

        InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices);

        GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
        GetAverageHiGainArea(j).SetNbins(fAverageNbins);
      }
  }

  if (fAverageLoGainAreas->GetEntries()==0)
  {
    fAverageLoGainAreas->Expand(nareas);
    
    for (Int_t j=0; j<nareas; j++)
      {
        (*fAverageLoGainAreas)[j] = 
          new MHPedestalPix("AverageLoGainArea",
                                      "Pedestals average Area idx ");

        InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices);

        GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
        GetAverageLoGainArea(j).SetNbins(fAverageNbins);
      }
  }

  if (fAverageHiGainSectors->GetEntries()==0)
  {
      fAverageHiGainSectors->Expand(nsectors);

      for (Int_t j=0; j<nsectors; j++)
      {
          (*fAverageHiGainSectors)[j] = 
            new MHPedestalPix("AverageHiGainSector",
                                        "Pedestals average sector ");

          InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices);

          GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
          GetAverageHiGainSector(j).SetNbins(fAverageNbins);
      }
  }

  if (fAverageLoGainSectors->GetEntries()==0)
  {
      fAverageLoGainSectors->Expand(nsectors);

      for (Int_t j=0; j<nsectors; j++)
      {
          (*fAverageLoGainSectors)[j] = 
            new MHPedestalPix("AverageLoGainSector",
                                        "Pedestals average sector ");

          InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices);
          
          GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
          GetAverageLoGainSector(j).SetNbins(fAverageNbins);
      }
  }

  return kTRUE;
}


// -------------------------------------------------------------
//
// If MBadPixelsPix::IsBad():
// - calls MHGausEvents::SetExcluded()
//
// Calls:
// - MHGausEvents::InitBins()
// - MHGausEvents::ChangeHistId(i)
// - MHGausEvents::SetEventFrequency(fPulserFrequency)
// - MHPedestalPix::SetNSlices(nslices)
//
void MHPedestalCam::InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices)
{
  hist.InitBins();
  hist.ChangeHistId(i);
  hist.SetEventFrequency(fPulserFrequency);
  hist.SetNSlices(nslices);
  hist.SetProbLimit(0.);
}



// -------------------------------------------------------------------------------
//
// Retrieves pointer to MExtractedSignalCam:
//
// Retrieves from MGeomCam:
// - number of pixels
// - number of pixel areas
// - number of sectors
//
// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively
// with the signals MExtractedSignalCam::GetExtractedSignalHiGain() and 
// MExtractedSignalCam::GetExtractedSignalLoGain(), respectively.
//
Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
  if (!signal)
    {
      gLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
      return kFALSE;
    }
  
  const UInt_t npixels  = fGeom->GetNumPixels();
  const UInt_t nareas   = fGeom->GetNumAreas();
  const UInt_t nsectors = fGeom->GetNumSectors();

  TArrayF sumareahi(nareas);
  TArrayF sumarealo(nareas);
  TArrayF sumsectorhi(nsectors);
  TArrayD sumsectorlo(nsectors);
  TArrayI numareahi(nareas);
  TArrayI numarealo(nareas);
  TArrayI numsectorhi(nsectors);
  TArrayI numsectorlo(nsectors);

  for (UInt_t i=0; i<npixels; i++)
  {
      MHGausEvents &histhi = (*this)[i];
      MHGausEvents &histlo = (*this)(i);

      if (histhi.IsExcluded())
          continue;

      const MExtractedSignalPix &pix = (*signal)[i];

      const Float_t pedhi = pix.GetExtractedSignalHiGain();
      const Float_t pedlo = pix.GetExtractedSignalLoGain();

      const Int_t aidx   = (*fGeom)[i].GetAidx();
      const Int_t sector = (*fGeom)[i].GetSector();

      histhi.FillHistAndArray(pedhi) ;
      sumareahi  [aidx]   += pedhi;
      numareahi  [aidx]   ++;
      sumsectorhi[sector] += pedhi;
      numsectorhi[sector] ++;

      histlo.FillHistAndArray(pedlo);
      sumarealo  [aidx]   += pedlo;
      numarealo  [aidx]   ++;
      sumsectorlo[sector] += pedlo;
      numsectorlo[sector] ++;

    }
  
  for (UInt_t j=0; j<nareas; j++)
    {
      MHGausEvents &histhi = GetAverageHiGainArea(j);
      histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);

      MHGausEvents &histlo = GetAverageLoGainArea(j);
      histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]);
    }
  
  for (UInt_t j=0; j<nsectors; j++)
    {
      MHGausEvents &histhi = GetAverageHiGainSector(j);
      histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);

      MHGausEvents &histlo = GetAverageLoGainSector(j);
      histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]);
    }
  
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calls:
// - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
// - MHCalibrationCam::FitLoGainArrays() with Bad Pixels flags 0
// 
Bool_t MHPedestalCam::FinalizeHists()
{

  FitHiGainArrays((*fCam),*fBadPixels, 
                  MBadPixelsPix::kHiGainNotFitted,
                  MBadPixelsPix::kHiGainOscillating);
  FitLoGainArrays((*fCam),*fBadPixels,
                  MBadPixelsPix::kLoGainNotFitted,
                  MBadPixelsPix::kLoGainOscillating);


  return kTRUE;
  

}

// ------------------------------------------------------------------
//
// The types are as follows:
// 
// Fitted values:
// ============== 
//
// 0: Fitted Charge
// 1: Error of fitted Charge
// 2: Sigma of fitted Charge
// 3: Error of Sigma of fitted Charge
//
//
// Useful variables derived from the fit results:
// =============================================
//
// 4: Returned probability of Gauss fit to Charge distribution
// 5: Relative differenc of calculated pedestal (per slice) and fitted (per slice)
// 6: Error of the Relative differenc of calculated pedestal (per slice) and fitted (per slice)
// 7: Relative difference of the error of the mean pedestal (per slice) - calculated and fitted
// 8: Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
// 9: Error of Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
// 10: Relative difference of the error of the pedestal RMS (per slice) - calculated and fitted
//
// Localized defects:
// ==================
//
// 11: Gaus fit not OK
// 12: Fourier spectrum not OK
//
Bool_t MHPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{

  if (fHiGainArray->GetSize() <= idx)
    return kFALSE;

  if ((*this)[idx].IsExcluded())
    return kFALSE;

  const Float_t ped      = (*fPedestals)[idx].GetPedestal();
  const Float_t rms      = (*fPedestals)[idx].GetPedestalRms();

  const Float_t entsqr    =  TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());

  const Float_t pederr   = rms/entsqr;
  const Float_t rmserr   = rms/entsqr/2.;

  const Float_t mean     = (*this)[idx].GetMean();
  const Float_t meanerr  = (*this)[idx].GetMeanErr();
  const Float_t sigma    = (*this)[idx].GetSigma() ;
  const Float_t sigmaerr = (*this)[idx].GetSigmaErr();

  switch (type)
    {
    case 0:
      val = mean;
      break;
    case 1:
      val = meanerr;
      break;
    case 2:
      val = sigma;
      break;
    case 3:
      val = sigmaerr;
      break;
    case 4:
      val = (*this)[idx].GetProb();
      break;
    case 5:
      val = 2.*(mean-ped)/(ped+mean);
      break;
    case 6:
      val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean))
            *2./(ped+mean)/(ped+mean);
      break;
    case 7:
      val = 2.*(meanerr-pederr)/(pederr + meanerr);
      break;
    case 8:
      val = 2.*(sigma-rms)/(sigma+rms);
      break;
    case 9:
      val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma))
            *2./(rms+sigma)/(rms+sigma);
      break;
    case 10:
      val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr);
      break;
    case 11:
      if (!(*this)[idx].IsGausFitOK())
        val = 1.;
      break;
    case 12:
      if (!(*this)[idx].IsFourierSpectrumOK())
        val = 1.;
      break;
    default:
      return kFALSE;
    }
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calls MHGausEvents::DrawClone() for pixel idx
//
void MHPedestalCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].DrawClone();
}