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

/* ======================================================================== *\
!
! *
! * 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-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  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";

  fADC2Phot = 0.;
  fEvents = 0;

  fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500., -3., 3.);
  fHistRatio->GetXaxis()->SetTitle("log_{10}(fPassPhotCone / fSize) (in photons/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 << dbginf << "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)
{
  //
  // If it is no MC file display error and exit
  //
  if (!CheckRunType(pList))
    {
      *fLog << err << dbginf << "This is no MC file... aborting." << endl;
      return kFALSE;
    }

  fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam"));

  if ( !fCalCam )
    {
      *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MCalibrationCam") << "... aborting." << endl;
      return kFALSE;
    }

  fHillas = (MHillas*) pList->FindCreateObj(AddSerialNumber("MHillas"));
  if ( !fHillas)
    {
      *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MHillas") << "... aborting." << endl;
      return kFALSE;
    }

  fNew = (MNewImagePar*) pList->FindCreateObj(AddSerialNumber("MNewImagePar"));
  if ( !fNew)
    {
      *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MNewImagePar") << "... aborting." << endl;
      return kFALSE;
    }

  fMcEvt = (MMcEvt*) pList->FindCreateObj(AddSerialNumber("MMcEvt"));
  if ( !fMcEvt)
    {
      *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MMcEvt") << "... aborting." << endl;
      return kFALSE;
    }

  return kTRUE;

}

// --------------------------------------------------------------------------
//
// Check for the runtype.
// Search for MGeomCam and MMcFadcHeader.
//
Bool_t MMcCalibrationCalc::ReInit(MParList *pList)
{
  //
  // Now check the existence of all necessary containers.
  //

  fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
  if ( ! fGeom )
    {
      *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MGeomCam") << "... aborting." << endl;
      return kFALSE;
    }

  fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
  if (!fHeaderFadc)
    {
      *fLog << err << dbginf << 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 << endl << endl << dbginf << "You are trying to calibrate the data using a Camera file produced with added noise. 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();
  fEvents ++;

  fHistRatio->Fill(log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize()));

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the MCalibrationCam object
//
Int_t MMcCalibrationCalc::PostProcess()
{

  if (fEvents > 0)
    fADC2Phot /= fEvents;
  else
    {
      *fLog << err << dbginf << "No events were read! Aborting." << endl;
      return kFALSE;
    }

  //
  // For the calibration we no longer use the mean, but thepeak of the distribution:
  //

  Float_t summax = 0.;
  Int_t mode = 0;
  Int_t reach = 2;
  for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++)
    {
      Float_t sum = 0;
      for(Int_t k = ibin-reach; k <= ibin+reach; k++)
        sum += fHistRatio->GetBinContent(k);
      if (sum > summax)
        {
          summax = sum;
          mode = ibin;
        }
    }

  fADC2Phot = pow(10., fHistRatio->GetXaxis()->GetBinCenter(mode));

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

      Float_t factor =  fADC2Phot*calpix.GetMeanConversionBlindPixelMethod();

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

  return kTRUE;

}