Index: trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.cc	(revision 4996)
+++ trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.cc	(revision 4996)
@@ -0,0 +1,581 @@
+/* ======================================================================== *\
+!
+! *
+! * 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 MHCalibrationPix-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
+// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
+// In case this does not make the fit valid, the histogram means and RMS's are
+// taken directly (see MHCalibrationPix::BypassFit()).
+//
+// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
+// from the mean are counted as Pickup events (stored in MHCalibrationPix::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"
+
+#include "TH1.h"
+
+ClassImp(MHPedestalCam);
+
+using namespace std;
+// --------------------------------------------------------------------------
+//
+// Default constructor. 
+//
+// Initializes:
+// - fExtractHiGainSlices to 0.
+// - fExtractLoGainSlices to 0.
+// - the event frequency to 1200 Hz.
+// - the fRenorm flag to kTRUE
+//
+MHPedestalCam::MHPedestalCam(const char *name, const char *title) 
+    : fExtractHiGainSlices(0.), fExtractLoGainSlices(0.)
+{
+
+  fName  = name  ? name  : "MHPedestalCam";
+  fTitle = title ? title : "";
+  
+  SetPulserFrequency(1200);
+  SetRenorm();
+
+}
+
+
+
+// --------------------------------------------------------------------------
+//
+// 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
+// 
+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 ");
+
+        GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
+
+        InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices);
+
+      }
+  }
+
+  if (fAverageLoGainAreas->GetEntries()==0)
+  {
+    fAverageLoGainAreas->Expand(nareas);
+    
+    for (Int_t j=0; j<nareas; j++)
+      {
+        (*fAverageLoGainAreas)[j] = 
+          new MHPedestalPix("AverageLoGainArea",
+                                      "Pedestals average Area idx ");
+
+        GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
+
+        InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices);
+
+      }
+  }
+
+  if (fAverageHiGainSectors->GetEntries()==0)
+  {
+      fAverageHiGainSectors->Expand(nsectors);
+
+      for (Int_t j=0; j<nsectors; j++)
+      {
+	  (*fAverageHiGainSectors)[j] = 
+            new MHPedestalPix("AverageHiGainSector",
+                                        "Pedestals average sector ");
+
+          GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
+
+          InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices);
+
+      }
+  }
+
+  if (fAverageLoGainSectors->GetEntries()==0)
+  {
+      fAverageLoGainSectors->Expand(nsectors);
+
+      for (Int_t j=0; j<nsectors; j++)
+      {
+	  (*fAverageLoGainSectors)[j] = 
+            new MHPedestalPix("AverageLoGainSector",
+                                        "Pedestals average sector ");
+
+          GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
+
+          InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices);
+          
+      }
+  }
+
+  return kTRUE;
+}
+
+
+// -------------------------------------------------------------
+//
+// If MBadPixelsPix::IsBad():
+// - calls MHCalibrationPix::SetExcluded()
+//
+// Calls:
+// - MHGausEvents::InitBins()
+// - MHCalibrationPix::ChangeHistId(i)
+// - MHCalibrationPix::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);
+
+  if (fRenorm)
+    hist.SetNSlices(nslices);
+
+  hist.SetProbLimit(0.);
+
+  TH1F *h = hist.GetHGausHist();
+  h->SetTitle( Form("%s%s", h->GetTitle()," Runs: "));
+}
+// -------------------------------------------------------------------------------
+//
+// 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++)
+  {
+      MHCalibrationPix &histhi = (*this)[i];
+      MHCalibrationPix &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++)
+    {
+      MHCalibrationPix &histhi = GetAverageHiGainArea(j);
+      histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
+
+      MHCalibrationPix &histlo = GetAverageLoGainArea(j);
+      histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]);
+    }
+  
+  for (UInt_t j=0; j<nsectors; j++)
+    {
+      MHCalibrationPix &histhi = GetAverageHiGainSector(j);
+      histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
+
+      MHCalibrationPix &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();
+}
Index: trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.h
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.h	(revision 4996)
+++ trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.h	(revision 4996)
@@ -0,0 +1,42 @@
+#ifndef MARS_MHPedestalCam
+#define MARS_MHPedestalCam
+
+#ifndef MARS_MHCalibrationCam
+#include "MHCalibrationCam.h"
+#endif
+
+class MGeomCam;
+class MPedestalCam;
+class MHPedestalPix;
+class MHPedestalCam : public MHCalibrationCam
+{
+
+private:
+
+  void   InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices);
+  
+  Bool_t ReInitHists(MParList *pList);
+  Bool_t FillHists(const MParContainer *par, const Stat_t w=1);
+  Bool_t FinalizeHists();
+  
+  Float_t fExtractHiGainSlices;           // Number of FADC slices used for high gain signal extraction
+  Float_t fExtractLoGainSlices;           // Number of FADC slices used for low  gain signal extraction
+
+  Bool_t  fRenorm;
+
+  MPedestalCam *fPedestals;               //! Pedestal Cam filled by MPedCalcPedRun
+  
+public:
+  MHPedestalCam(const char *name=NULL, const char *title=NULL);
+  ~MHPedestalCam() {}
+  
+  Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
+  void DrawPixelContent(Int_t idx) const;
+  
+  void SetRenorm( const Bool_t b=kTRUE )  {  fRenorm = b; }
+
+  ClassDef(MHPedestalCam, 1)	// Histogram class for Charge Camera Pedestals 
+};
+
+#endif
+