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

/* ======================================================================== *\
   !
   ! *
   ! * 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 "MParList.h"

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

#include "MCalibrationPix.h"
#include "MCalibrationCam.h"
#include "MGeomCam.h"
#include "MRawRunHeader.h"
#include "MMcFadcHeader.hxx"
#include "MHillas.h"
#include "MNewImagePar.h"
#include "MMcEvt.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("HistRatio", "Ratio fPassPhotCone/fSize", 1000., 0., 10.);
}

// --------------------------------------------------------------------------
//
// 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 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 MGeomCam... aborting." << endl;
      return kFALSE;
    }

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

  //
  // FIXME: Check that the relevant parameters in the FADC header are the
  // same for all the analyzed data!
  //

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Obtain average ratio of photons in camera to image Size.
//
Int_t MMcCalibrationCalc::Process()
{

  if ( fNew->GetNumSaturatedPixels() > 0 )
    return kTRUE;

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

  fHistRatio->Fill(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;
    }

  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;

}