source: trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc@ 3211

/* ======================================================================== *\
!
! *
! * 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): Abelardo Moralejo, 12/2003 <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MMcCalibrationCalc
//
//  Input Containers:
//   MRawRunHeader
//   MMcFadcHeader
//   MHillas
//   MNewImagePar
//   MMcEvt
//
//  Output Containers:
//   MCalibrationCam
//
/////////////////////////////////////////////////////////////////////////////
#include "MMcCalibrationCalc.h"

#include <TH1.h>

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

#include "MParList.h"

#include "MCalibrationPix.h"
#include "MCalibrationCam.h"

#include "MGeomCam.h"
#include "MRawRunHeader.h"

#include "MHillas.h"
#include "MNewImagePar.h"

#include "MMcEvt.hxx"
#include "MMcFadcHeader.hxx"

ClassImp(MMcCalibrationCalc);

using namespace std;

MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title)
{
    fName  = name  ? name  : "MMcCalibrationCalc";
    fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container";

    fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.);
    fHistRatio->SetXTitle("log_{10}(fPassPhotCone / fSize) [phot/ADC count]");
}

// --------------------------------------------------------------------------
//
// Check for the run type. Return kTRUE if it is a MC run or if there
// is no MC run header (old camera files) kFALSE in case of a different
// run type
//
Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const
{
    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!run)
    {
        *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
        return kTRUE;
    }

    return  run->GetRunType() == kRTMonteCarlo;
}

// --------------------------------------------------------------------------
//
// Make sure, that there is an MCalibrationCam Object in the Parameter List.
//
Int_t MMcCalibrationCalc::PreProcess(MParList *pList)
{
    fHistRatio->Reset();
    fADC2Phot = 0;

    fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam"));
    if (!fCalCam)
    {
        *fLog << err << AddSerialNumber("MCalibrationCam") << "not found... aborting." << endl;
        return kFALSE;
    }

    fHillas = (MHillas*) pList->FindObject(AddSerialNumber("MHillas"));
    if ( !fHillas)
    {
        *fLog << err << AddSerialNumber("MHillas") << "not found... aborting." << endl;
        return kFALSE;
    }

    fNew = (MNewImagePar*)pList->FindObject(AddSerialNumber("MNewImagePar"));
    if (!fNew)
    {
        *fLog << err << AddSerialNumber("MNewImagePar") << "not found... aborting." << endl;
        return kFALSE;
    }

    fMcEvt = (MMcEvt*) pList->FindObject(AddSerialNumber("MMcEvt"));
    if (!fMcEvt)
    {
        *fLog << err << AddSerialNumber("MMcEvt") << "not found... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Check for the runtype.
// Search for MGeomCam and MMcFadcHeader.
//
Bool_t MMcCalibrationCalc::ReInit(MParList *pList)
{
  //
  // If it is no MC file display error and exit
  //
  if (!CheckRunType(pList))
  {
      *fLog << err << "MMcCalibrationCalc can only be used with MC files... aborting." << endl;
      return kFALSE;
  }

  //
  // Now check the existence of all necessary containers.
  //
  fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
  if (!fGeom)
  {
      *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
      return kFALSE;
  }

  fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
  if (!fHeaderFadc)
  {
      *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
      return kFALSE;
  }

  for (UInt_t ipix = 0; ipix < fGeom->GetNumPixels(); ipix++)
  {
      if (fHeaderFadc->GetPedestalRmsHigh(ipix) > 0 ||
          fHeaderFadc->GetPedestalRmsLow(ipix)  > 0 )
      {
          *fLog << err << "Trying to calibrate the data using a Camera file produced with added noise." << endl;
          *fLog << "Please use a noiseless file for calibration... aborting." << endl << endl;
          return kFALSE;
      }
  }

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Obtain average ratio of photons in camera to image Size.
//
Int_t MMcCalibrationCalc::Process()
{
    //
    // Exclude events with some saturated pixel
    //
    if (fNew->GetNumSaturatedPixels()>0)
        return kTRUE;

    //
    // Exclude events with low Size (larger fluctuations)
    // FIXME? The present cut (1000 "inner-pixel-counts") is somehow
    // arbitrary. Might it be optimized?
    //
    if (fHillas->GetSize()<1000)
        return kTRUE;

    fADC2Phot += fMcEvt->GetPassPhotCone()/fHillas->GetSize();

    fHistRatio->Fill(TMath::Log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize()));

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the MCalibrationCam object
//
Int_t MMcCalibrationCalc::PostProcess()
{
    const Stat_t n = fHistRatio->GetEntries();
    if (n<1)
    {
        *fLog << err << "No events read... aborting." << endl;
        return kFALSE;
    }

    fADC2Phot /= n;

    //
    // For the calibration we no longer use the mean,
    // but the peak of the distribution:
    //
    const Int_t reach = 2;

    Stat_t summax = 0;
    Int_t  mode   = 0;

    // FIXME: Is this necessary? We could use GetMaximumBin instead..
    for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++)
    {
        const Stat_t sum = fHistRatio->Integral(ibin-reach, ibin+reach);

        if (sum <= summax)
            continue;

        summax = sum;
        mode = ibin;
    }

    fADC2Phot = TMath::Power(10, fHistRatio->GetBinCenter(mode));

    const Int_t num = fCalCam->GetSize();
    for (int i=0; i<num; i++)
    {
        MCalibrationPix &calpix = (*fCalCam)[i];

        const Float_t factor = fADC2Phot*calpix.GetMeanConversionBlindPixelMethod();

        calpix.SetConversionBlindPixelMethod(factor, 0., 0.);
    }

    return kTRUE;
}