source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityHiLoCam.cc@ 5895

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationIntensityHiLoCam                                               
//                                                               
// Storage container for intensity charge calibration results. 
// 
// Individual MCalibrationHiLoCam'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, MCalibrationHiLoCam,
//           MCalibrationHiLoPix, MCalibrationHiLoCalc, MCalibrationQECam
//           MHCalibrationHiLoPix, MHCalibrationHiLoCam              
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationIntensityHiLoCam.h"
#include "MCalibrationHiLoCam.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationHiLoPix.h"
#include "MCalibrationChargePix.h"

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

#include "MLogManip.h"

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

ClassImp(MCalibrationIntensityHiLoCam);

using namespace std;

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

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

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

TGraphErrors *MCalibrationIntensityHiLoCam::GetHiLoRatioVsCharge( const UInt_t pixid, const MCalibrationIntensityChargeCam &chargecam,
                                  const MCalibrationCam::PulserColor_t col)
{

  if (chargecam.GetSize() != GetSize())
    {
      *fLog << err << GetDescriptor() << ": Size mismatch between MCalibrationIntensityHiLoCam "
            << "and MCalibrationIntensityChargeCam. " << endl;
      return NULL;
    }
  
  Int_t size = CountNumEntries(col);
  
  if (size == 0)
    return NULL;

  if (size != chargecam.CountNumEntries(col))
    {
      *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityHiLoCam "
            << "and MCalibrationIntensityChargeCam. " << endl;
      return NULL;
    }
  
  const Int_t nvalid = chargecam.CountNumValidEntries(pixid,col);

  if (nvalid == 0)
    {
      *fLog << err << GetDescriptor() << ": Only un-calibrated events in pixel: " << pixid << endl;
      return NULL;
    }

  TArrayF ratio(nvalid);
  TArrayF ratioerr(nvalid);
  TArrayF sig(nvalid);
  TArrayF sigerr(nvalid);

  Int_t cnt = 0;
  
  for (Int_t i=0;i<GetSize();i++)
    {
      //
      // Get the calibration cam from the intensity cam
      //
      MCalibrationChargeCam *cam     = (MCalibrationChargeCam*)chargecam.GetCam(i);
      MCalibrationHiLoCam *hilocam = (MCalibrationHiLoCam*)GetCam(i);

      if (col != MCalibrationCam::kNONE)
        if (hilocam->GetPulserColor() != col)
          continue;
      //
      // Get the calibration pix from the calibration cam
      //
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)(*cam)[pixid];
      MCalibrationHiLoPix &relpix = (MCalibrationHiLoPix&)(*hilocam)[pixid];
      //
      // Don't use bad pixels
      //
      if (!pix.IsFFactorMethodValid())
        continue;
      //
      ratio[cnt]    = relpix.GetHiLoRatio();
      ratioerr[cnt] = relpix.GetHiLoRatioErr();
      //
      sig   [cnt] = pix.GetPheFFactorMethod();
      sigerr[cnt] = pix.GetPheFFactorMethodErr();
      cnt++;
    }

  TGraphErrors *gr = new TGraphErrors(nvalid,
                                     sig.GetArray(),ratio.GetArray(),
                                     sigerr.GetArray(),ratioerr.GetArray());
  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
  gr->GetXaxis()->SetTitle("<Photo-electrons> [1]");
  gr->GetYaxis()->SetTitle("HiLo Ratio [1]");      
  return gr;
}


TGraphErrors *MCalibrationIntensityHiLoCam::GetHiLoRatioVsChargePerArea( const Int_t aidx,const MCalibrationIntensityChargeCam &chargecam, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
{
  
  Int_t size = CountNumEntries(col);
  
  if (size == 0)
    return NULL;

  if (size != chargecam.CountNumEntries(col))
    {
      *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityHiLoCam "
            << "and MCalibrationIntensityChargeCam. " << endl;
      return NULL;
    }
  
  TArrayF ratio(size);
  TArrayF ratioerr(size);
  TArrayF sig(size);
  TArrayF sigerr(size);
  
  Int_t cnt = 0;

  for (Int_t i=0;i<GetSize();i++)
    {
      //
      // Get the calibration cam from the intensity cam
      //
      MCalibrationHiLoCam *hilocam = (MCalibrationHiLoCam*)GetCam(i);
      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);

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

      const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
      const Float_t phe          = apix.GetPheFFactorMethod();
      const Float_t pheerr       = apix.GetPheFFactorMethodErr();

      sig[cnt]       = phe;
      sigerr[cnt]    = pheerr;

      Double_t ratiool  = 0.;
      Double_t ratiool2 = 0.;
      Double_t var    = 0.;
      Int_t    num    = 0;
     //
      // Get the area calibration pix from the calibration cam
      //
      for (Int_t i=0; i<cam->GetSize(); i++)
        {
          const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[i];
          const MCalibrationHiLoPix &relpix = (MCalibrationHiLoPix&)(*hilocam)[i];
          //
          // Don't use bad pixels
          //
          if (!pix.IsFFactorMethodValid())
            continue;
          //
          //
          if (aidx != geom[i].GetAidx())
            continue;
          
          ratiool  += relpix.GetHiLoRatio();
          ratiool2 += relpix.GetHiLoRatio()*relpix.GetHiLoRatio();
          num++;
        }
      
      if (num > 1)
        {
          ratiool /= num;
          var  = (ratiool2 - ratiool*ratiool*num) / (num-1);

          ratio[cnt] = ratiool;
          if (var > 0.)
            ratioerr[cnt] = TMath::Sqrt(var);
          else
            ratioerr[cnt] = 0.;
        }
      else
        {
          ratio[cnt]    = -1.;
          ratioerr[cnt] = 0.;
        }
      cnt++;
    }
  
  TGraphErrors *gr = new TGraphErrors(size,
                                      sig.GetArray(),ratio.GetArray(),
                                      sigerr.GetArray(),ratioerr.GetArray());
  gr->SetTitle(Form("%s%3i","Area Index ",aidx));
  gr->GetXaxis()->SetTitle("<phes> [1]");
  gr->GetYaxis()->SetTitle("HiLo Ratio [1]");      
  return gr;
}