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

/* ======================================================================== *\
!
! *
! * 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:
//  (containers mus exist already, they are filled with new values).
//   MCalibrationChargeCam
//   MCalibrationQECam  
//
/////////////////////////////////////////////////////////////////////////////
#include "MMcCalibrationCalc.h"

#include <TH1.h>

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

#include "MParList.h"

#include "MCalibrationChargePix.h"
#include "MCalibrationChargeCam.h"

#include "MCalibrationQEPix.h"
#include "MCalibrationQECam.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 MCalibrationChargeCam and MCalibrationQECam containers";

    fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPhotElfromShower/fSize)", 1500, -3., 3.);
    fHistRatio->SetXTitle("log_{10}(fPhotElfromShower / fSize) [photel/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->IsMonteCarloRun();
}

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

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

    fQECam = (MCalibrationQECam*) pList->FindObject(AddSerialNumber("MCalibrationQECam"));
    if (!fQECam)
    {
        *fLog << err << AddSerialNumber("MCalibrationQECam") << "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?
    //   

    const Float_t size = fHillas->GetSize(); 
    // Size will at this point be in ADC counts (still uncalibrated)

    if (size < 1000)
        return kTRUE;

    fPhot2PhotEl += (Float_t) fMcEvt->GetPhotElfromShower() /
      (Float_t) fMcEvt->GetPassPhotCone();

    fHistRatio->Fill(TMath::Log10(fMcEvt->GetPhotElfromShower()/size));

    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;
    }

    fPhot2PhotEl /= n;   // Average quantum efficiency

    //
    // Find peak of log10(photel/Size) histogram:
    //
    const Int_t reach = 2;
    Stat_t summax = 0;
    Int_t  mode   = 0;
    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;
    }

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

    const Int_t num = fCalCam->GetSize();

    for (int i=0; i<num; i++)
    {
        MCalibrationChargePix &calpix = (MCalibrationChargePix&)(*fCalCam)[i];
        MCalibrationQEPix     &qepix  = (MCalibrationQEPix&)    (*fQECam) [i];

        qepix.SetAverageQE(fPhot2PhotEl);

        qepix.SetAvNormFFactor(1.);
        // This factor should convert the default average QE for different colors to 
        // average QE for a spectrum like that of Cherenkov light (see the documentration 
        // of MCalibrationQEPix). 
        // Here we obtain average QE using already a Cherenkov spectrum so AvNormFFacto 
        // must be 1.

        
        Float_t factor = fADC2PhotEl;

        //
        // We take into account the (possibly) different gain of outer pixels:
        //
        if (fGeom->GetPixRatio(i) < 1.)
          factor *= fHeaderFadc->GetAmplitud()/fHeaderFadc->GetAmplitudOuter();

        calpix.SetMeanConvFADC2Phe(factor);
        calpix.SetMeanConvFADC2PheVar(0.);

        calpix.SetMeanFFactorFADC2Phot(0.);
    }

    return kTRUE;
}