source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityRelTimeCam.cc@ 5710

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationIntensityRelTimeCam                                               
//                                                               
// Storage container for intensity charge calibration results. 
// 
// Individual MCalibrationRelTimeCam'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, MCalibrationRelTimeCam,
//           MCalibrationRelTimePix, MCalibrationRelTimeCalc, MCalibrationQECam
//           MHCalibrationRelTimePix, MHCalibrationRelTimeCam              
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationIntensityRelTimeCam.h"
#include "MCalibrationRelTimeCam.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationRelTimePix.h"
#include "MCalibrationChargePix.h"

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

#include "MLogManip.h"

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

ClassImp(MCalibrationIntensityRelTimeCam);

using namespace std;

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

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

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

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

  if (chargecam.GetSize() != GetSize())
    {
      *fLog << err << GetDescriptor() << ": Size mismatch between MCalibrationIntensityRelTimeCam "
            << "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 MCalibrationIntensityRelTimeCam "
            << "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 res(nvalid);
  TArrayF reserr(nvalid);
  TArrayF sig(nvalid);
  TArrayF sigerr(nvalid);

  const Float_t sqrt2   = 1.414;
  const Float_t fadc2ns = 3.333;
  
  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);
      MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);

      if (col != MCalibrationCam::kNONE)
        if (relcam->GetPulserColor() != col)
          continue;
      //
      // Get the calibration pix from the calibration cam
      //
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)(*cam)[pixid];
      MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[pixid];
      //
      // Don't use bad pixels
      //
      if (!pix.IsFFactorMethodValid())
        continue;
      //
      res[cnt]    = relpix.GetTimePrecision() / sqrt2 * fadc2ns;
      reserr[cnt] = relpix.GetTimePrecisionErr() / sqrt2 * fadc2ns;
      //
      sig   [cnt] = pix.GetPheFFactorMethod();
      sigerr[cnt] = pix.GetPheFFactorMethodErr();
      cnt++;
    }

  TGraphErrors *gr = new TGraphErrors(nvalid,
                                     sig.GetArray(),res.GetArray(),
                                     sigerr.GetArray(),reserr.GetArray());
  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
  gr->GetXaxis()->SetTitle("<Photo-electrons> [1]");
  gr->GetYaxis()->SetTitle("Time Resolution [ns]");      
  return gr;
}


TGraphErrors *MCalibrationIntensityRelTimeCam::GetTimeResVsChargePerArea( 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 MCalibrationIntensityRelTimeCam "
            << "and MCalibrationIntensityChargeCam. " << endl;
      return NULL;
    }
  
  const Float_t sqrt2   = 1.414;
  const Float_t fadc2ns = 3.333;
  
  TArrayF res(size);
  TArrayF reserr(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
      //
      MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);

      if (col != MCalibrationCam::kNONE)
        if (relcam->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 resol  = 0.;
      Double_t resol2 = 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 MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[i];
          //
          // Don't use bad pixels
          //
          if (!pix.IsFFactorMethodValid())
            continue;
          //
          //
          if (aidx != geom[i].GetAidx())
            continue;
          
          resol  += relpix.GetTimePrecision();
          resol2 += relpix.GetTimePrecision()*relpix.GetTimePrecision();
          num++;
        }
      
      if (num > 1)
        {
          resol /= num;
          var  = (resol2 - resol*resol*num) / (num-1);

          res[cnt] = resol * fadc2ns / sqrt2;
          if (var > 0.)
            reserr[cnt] = TMath::Sqrt(var)/ sqrt2 * fadc2ns;
          else
            reserr[cnt] = 0.;
        }
      else
        {
          res[cnt]    = -1.;
          reserr[cnt] = 0.;
        }
      cnt++;
    }
  
  TGraphErrors *gr = new TGraphErrors(size,
                                      sig.GetArray(),res.GetArray(),
                                      sigerr.GetArray(),reserr.GetArray());
  gr->SetTitle(Form("%s%3i","Area Index ",aidx));
  gr->GetXaxis()->SetTitle("<phes> [1]");
  gr->GetYaxis()->SetTitle("Time Resolution [ns]");      
  return gr;
}