source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc@ 5611

/* ======================================================================== *\
!
! *
! * 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-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationIntensityChargeCam                                               
//                                                               
// Storage container for intensity charge calibration results. 
// 
// Individual MCalibrationChargeCam's can be retrieved with: 
// - GetCam() yielding the current cam.
// - GetCam("name") yielding the current camera with name "name".
// - GetCam(i) yielding the i-th camera.
//
// See also: MCalibrationIntensityCam, MCalibrationChargeCam,
//           MCalibrationChargePix, MCalibrationChargeCalc, MCalibrationQECam
//           MCalibrationBlindCam, MCalibrationChargePINDiode
//           MHCalibrationChargePix, MHCalibrationChargeCam              
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationIntensityChargeCam.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"

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

#include <TOrdCollection.h>
#include <TGraphErrors.h>
#include <TH2F.h>
#include <TF1.h>

ClassImp(MCalibrationIntensityChargeCam);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationIntensityChargeCam::MCalibrationIntensityChargeCam(const char *name, const char *title)
{

  fName  = name  ? name  : "MCalibrationIntensityChargeCam";
  fTitle = title ? title : "Results of the Intensity Calibration";
  
  InitSize(1);
}

// -------------------------------------------------------------------
//
// Add MCalibrationChargeCam's in the ranges from - to. 
//
void MCalibrationIntensityChargeCam::Add(const UInt_t from, const UInt_t to)
{
  for (UInt_t i=from; i<to; i++)
    fCams->AddAt(new MCalibrationChargeCam,i);
}

// -------------------------------------------------------------------
//
// Returns a TGraphErrors with the number of photo-electrons vs. 
// the extracted signal of pixel "pixid". 
//
TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
{
  
  const Int_t size = GetSize();
  
  TArrayF phe(size);
  TArrayF pheerr(size);
  TArrayF sig(size);
  TArrayF sigerr(size);
  
  for (Int_t i=0;i<size;i++)
    {
      //
      // Get the calibration cam from the intensity cam
      //
      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);

      if (col != MCalibrationCam::kNONE)
        if (cam->GetPulserColor() != col)
          continue;
      //
      // Get the calibration pix from the calibration cam
      //
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
      //
      // Don't use bad pixels
      //
      if (!pix.IsFFactorMethodValid())
        continue;
      //
      phe[i]    = pix.GetPheFFactorMethod();
      pheerr[i] = pix.GetPheFFactorMethodErr();
      //
      // For the calculation of Q, we have to use the 
      // converted value!
      //
      sig   [i] = pix.GetConvertedMean();
      sigerr[i] = pix.GetConvertedMeanErr();
    }
  
  TGraphErrors *gr = new TGraphErrors(size,
                                     sig.GetArray(),phe.GetArray(),
                                     sigerr.GetArray(),pheerr.GetArray());
  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
  gr->GetXaxis()->SetTitle("Q [FADC counts]");
  gr->GetYaxis()->SetTitle("photo-electrons [1]");      
  return gr;
}

// -------------------------------------------------------------------
//
// Returns a TGraphErrors with the number of photo-electrons vs. 
// the extracted signal over all pixels with area index "aidx". 
//
// The points represent the means of the pixels values, while the error bars
// the sigma of the pixels values.  
//
TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
{
  
  const Int_t size = GetSize();
  
  TArrayF phe(size);
  TArrayF pheerr(size);
  TArrayF sig(size);
  TArrayF sigerr(size);
  
  for (Int_t i=0;i<size;i++)
    {
      //
      // Get the calibration cam from the intensity cam
      //
      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);

      if (col != MCalibrationCam::kNONE)
        if (cam->GetPulserColor() != col)
          continue;

      //
      // Get the area calibration pix from the calibration cam
      //
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(cam->GetAverageArea(aidx));

      phe[i]    = pix.GetPheFFactorMethod();
      pheerr[i] = pix.GetPheFFactorMethodErr();
      //
      // For the calculation of Q, we have to use the 
      // converted value!
      //
      sig   [i] = pix.GetConvertedMean();
      sigerr[i] = pix.GetConvertedMeanErr();
    }
  
  TGraphErrors *gr = new TGraphErrors(size,
                                      sig.GetArray(),phe.GetArray(),
                                      sigerr.GetArray(),pheerr.GetArray());
  gr->SetTitle(Form("%s%3i","Area Index ",aidx));
  gr->GetXaxis()->SetTitle("Q [FADC counts]");
  gr->GetYaxis()->SetTitle("photo-electrons [1]");      
  return gr;
}

// -------------------------------------------------------------------
//
// Returns a TGraphErrors with the 'Razmik plot' of pixel "pixid". 
// The Razmik plot shows the value of 'R' vs. 1/Q where:
//
//        sigma^2         F^2
//  R =   -------    =  ------
//         <Q>^2        <m_pe>
//
// and 1/Q is the inverse (mean) extracted signal
//
TGraphErrors *MCalibrationIntensityChargeCam::GetRazmikPlot( const UInt_t pixid )
{
  
  const Int_t size = GetSize();
  
  TArrayF r(size);
  TArrayF rerr(size);
  TArrayF oneoverq(size);
  TArrayF oneoverqerr(size);
  
  for (Int_t i=0;i<size;i++)
    {
      //
      // Get the calibration cam from the intensity cam
      //
      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
      //
      // Get the calibration pix from the calibration cam
      //
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
      //
      // Don't use bad pixels
      //
      if (!pix.IsFFactorMethodValid())
        continue;
      //
      // For the calculation of R, use the un-converted values, like 
      // in the calibration, since: 
      //                 C^2*sigma^2     sigma^2
      //  R(lowgain) =   -----------  =  ------  = R
      //                 C^2*<Q>^2       <Q>^2
      //
      const Float_t mean = pix.GetMean();
      const Float_t meanerr = pix.GetMeanErr();
      const Float_t rsigma = pix.GetRSigma();
      const Float_t rsigmaerr = pix.GetRSigmaErr();
      r[i]    = rsigma*rsigma/mean/mean;
      const Float_t rrelvar = 4.*rsigmaerr*rsigmaerr/rsigma/rsigma + 4.*meanerr*meanerr/mean/mean;
      rerr[i] = rrelvar * r[i] * r[i];
      rerr[i] = rerr[i] <= 0 ? 0. : TMath::Sqrt(rerr[i]);
      //
      // For the calculation of 1/Q, we have to use the 
      // converted value!
      //
      const Float_t q  = pix.GetConvertedMean();
      const Float_t qe = pix.GetConvertedMeanErr();
      oneoverq   [i] = 1./q;
      oneoverqerr[i] = qe / q / q;
    }
  
  TGraphErrors *gr = new TGraphErrors(size,
                                     oneoverq.GetArray(),r.GetArray(),
                                     oneoverqerr.GetArray(),rerr.GetArray());
  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
  gr->GetXaxis()->SetTitle("1/Q [FADC counts^{-1}]");
  gr->GetYaxis()->SetTitle("\sigma_{red}^{2}/Q^{2}");      
  return gr;
}

// -------------------------------------------------------------------
//
// Returns a 2-dimensional histogram with the fit results of the 
// 'Razmik plot' for each pixel of area index "aidx" (see GetRazmikPlot())
// 
// The results of the polynomial fit of grade 1 are: 
//
// x-axis: Offset (Parameter 0 of the polynomial)
// y-axis: Slope  (Parameter 1 of the polynomial)
//
// The offset is a measure of how well-known the supposed additional contributions 
// to the value "reduced sigma" are. Because a photo-multiplier is a linear instrument, 
// the excess fluctuations are linear w.r.t. the signal amplitude and can be expressed by 
// the proportionality constant F (the "F-Factor"). 
// Adding noise from outside (e.g. night sky background) modifies the recorded noise, but 
// not the mean extracted signal, due to the AC-coupling. Thus, noise contributions from outside
// (e.g. calculating the pedestal RMS)have to be subtracted from the recorded signal fluctuations 
// in order to retrieve the linearity relation: 
//
// sigma(signal)^2 / mean(signal)^2  = sigma^2 / <Q>^2 = F^2 / <n_phe>               (1) 
//
// Any systematic offset in the sigma(signal) will produce an offset in the "Razmik plot"), 
// characterized by the Offset of the polynomial fit. Thus, in an ideal case, all pixels have their
// "offset" centered very closely around zero.
//
// The "slope" is the proportionality constant F^2, multiplied with the conversion factor 
// phe's to mean signal (because the "Razmik plot" plots the left side of eq. (1) w.r.t. 
// 1/<Q> instead of 1/<n_phe>. However, the mean number of photo-electrons <n_phe> can be 
// expressed by <Q> with the relation:
//
//  <n_phe> = c_phe * <Q>                                                            (2)
//
// Thus: 
//
// 1/<n_phe> = 1/c_phe * 1/<Q>                                                       (3) 
// 
// and:
//
// Slope = F^2 / c_phe
//
// In the ideal case of having equal photo-multipliers and a perfectly flat-fielded camera,
// the "slope" -values should thus all be closely centered around F^2/c_phe. 
// 
TH2F *MCalibrationIntensityChargeCam::GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom)
{
  
  TH2F *hist = new TH2F("hist","R vs. Inverse Charges - Fit results",45,-0.02,0.02,45,0.,30.);
  hist->SetXTitle("Offset [FADC counts^{-1}]");
  hist->SetYTitle("F^{2} / <n_phe>/<Q> [FADC count / phe]");  
  hist->SetFillColor(kRed+aidx);

  MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();

  for (Int_t npix=0;npix<cam->GetSize();npix++)
    {

      if (geom[npix].GetAidx() == aidx)
        {
          TGraph *gr = GetRazmikPlot(npix);
          gr->Fit("pol1","Q");
          hist->Fill(gr->GetFunction("pol1")->GetParameter(0),gr->GetFunction("pol1")->GetParameter(1));
        }
    }
  return hist;
}