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

/* ======================================================================== *\
!
! *
! * 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@uni-sw.gwdg.de>
!              Markus Gaug   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MHPedestalCam                                                            //
//                                                                         //
// Hold the Pedestal information for all pixels in the camera              //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MHPedestalCam.h"

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

#include "MParList.h"

#include "MHPedestalPix.h"

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

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

ClassImp(MHPedestalCam);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. Creates a MHPedestalPix object for each pixel
//
MHPedestalCam::MHPedestalCam(const char *name, const char *title) 
    : fExtractSlices(0.)
{
    fName  = name  ? name  : "MHPedestalCam";
    fTitle = title ? title : "";

    //
    // loop over all Pixels and create two histograms
    // one for the Low and one for the High gain
    // connect all the histogram with the container fHist
    //
    fArray = new TObjArray;
    fArray->SetOwner();

}

// --------------------------------------------------------------------------
//
// Delete the array conatining the pixel pedest information
//
MHPedestalCam::~MHPedestalCam()
{
  delete fArray;
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MHPedestalPix &MHPedestalCam::operator[](UInt_t i)
{
  return *(MHPedestalPix*)(fArray->At(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
const MHPedestalPix &MHPedestalCam::operator[](UInt_t i) const
{
  return *(MHPedestalPix*)(fArray->At(i));
}


// --------------------------------------------------------------------------
//
// Our own clone function is necessary since root 3.01/06 or Mars 0.4
// I don't know the reason
//
TObject *MHPedestalCam::Clone(const char *) const
{

  const Int_t n = fArray->GetSize();
  
  //
  // FIXME, this might be done faster and more elegant, by direct copy.
  //
  MHPedestalCam *cam = new MHPedestalCam;
  
  cam->fArray->Expand(n);
  
  for (int i=0; i<n; i++)
    {
      delete (*cam->fArray)[i];
      (*cam->fArray)[i] = (*fArray)[i]->Clone();
    }

  return cam;
}



// --------------------------------------------------------------------------
//
// To setup the object we get the number of pixels from a MGeomCam object
// in the Parameter list. 
// MHPedestalPix sets its parameters to 0. (other than default which is -1.)
//
Bool_t MHPedestalCam::SetupFill(const MParList *pList)
{

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

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

  fArray->Delete();
  return kTRUE;

}

// --------------------------------------------------------------------------
//
Bool_t MHPedestalCam::Fill(const MParContainer *par, const Stat_t w)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
  if (!signal)
    {
      gLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
      return kFALSE;
    }
  

  Float_t slices = signal->GetNumUsedHiGainFADCSlices();

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

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

  fExtractSlices = slices;

  const Int_t n = signal->GetSize();
  
  if (fArray->GetEntries()==0)
    {
      fArray->Expand(n);
      
      for (Int_t i=0; i<n; i++)
        {
          (*fArray)[i] = new MHPedestalPix;
          MPedestalPix  &pix  = (*fPedestals)[i];
          pix.InitUseHists();
          MHPedestalPix &hist = (*this)[i];
          hist.ChangeHistId(i);
          hist.InitBins();
        }
    }
  
  if (fArray->GetEntries() != n)
    {
      gLog << err << "ERROR - Size mismatch... abort." << endl;
      return kFALSE;
    }
  
  for (Int_t i=0; i<n; i++)
    {

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

      MHPedestalPix  &hist = (*this)[i];
      //
      // Don't fill signal per slice, we get completely screwed up 
      // with the sigma. Better fill like it is and renorm later
      //
      hist.FillHistAndArray(signal);
    }
  
  return kTRUE;
}

Bool_t MHPedestalCam::Finalize()
{
    for (Int_t i=0; i<fArray->GetSize(); i++)
    {

        MHPedestalPix &hist = (*this)[i];

        //
        // 1) Return if the charge distribution is already succesfully fitted
        //    or if the histogram is empty
        //
        if (hist.IsGausFitOK() || hist.IsEmpty())
          continue;

        //
        // 2) Fit the Hi Gain histograms with a Gaussian
        //
        hist.FitGaus();
        hist.Renorm(fExtractSlices);
        hist.CreateFourierSpectrum();

    }
    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 (fArray->GetSize() <= idx)
    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.*(ped-mean)/(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.*(pederr - meanerr)/(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;
}

void MHPedestalCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].DrawClone();
}