source: trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc@ 4241

/* ======================================================================== *\
!
! *
! * 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): Javier Lopez    12/2003 <mailto:jlopez@ifae.es>
!   Author(s): Javier Rico     01/2004 <mailto:jrico@ifae.es>
!   Author(s): Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Markus Gaug     04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrateData
//
//   This task takes the integrated charge from MExtractedSignal and apply
//   the calibration constants from MCalibraitionCam to the charge. Then
//   stores number of photons obtained in MCerPhotEvt. Selection of different
//   calibration methods is allowed through SetCalibrationMode member function
//
//   in ReInit the MPedPhotCam container is filled using the information from
//   MPedestalCam, MExtractedSignalCam and MCalibrationCam
//
//   Input Containers:
//    MPedestalCam
//    MExtractedSingalCam
//    MCalibrationChargeCam
//    MCalibrationQECam
//
//   Output Containers:
//    MPedPhotCam
//    MCerPhotEvt
//
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrateData.h"

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

#include "MParList.h"
#include "MH.h"

#include "MGeomCam.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

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

#include "MCalibrationQECam.h"
#include "MCalibrationQEPix.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

#include "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

#include "MCerPhotEvt.h"

ClassImp(MCalibrateData);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title) 
    : fCam(NULL), fPedestal(NULL), fBadPixels(NULL), fCalibrations(NULL), fSignals(NULL), 
      fPedPhot(NULL), fCerPhotEvt(NULL), fCalibrationMode(calmode)
{
    fName  = name  ? name  : "MCalibrateData";
    fTitle = title ? title : "Task to calculate the number of photons in one event";
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MGeomCam
//  - MPedestalCam
//  - MCalibrationChargeCam
//  - MCalibrationQECam
//  - MExtractedSignalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MPedPhotCam
//  - MCerPhotEvt
//
Int_t MCalibrateData::PreProcess(MParList *pList)
{
    fPedestal = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
    if (!fPedestal)
    {
      *fLog << err << AddSerialNumber("MPedestalCam") << " not found ... aborting" << endl;
        return kFALSE;
    }


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

    fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));

    if (!fBadPixels)
    {
      *fLog << err << AddSerialNumber("MBadPixelsCam") << " not found ... aborting" << endl;
        return kFALSE;
    }

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

    fPedPhot = (MPedPhotCam*)pList->FindCreateObj(AddSerialNumber("MPedPhotCam"));
    if (!fPedPhot)
        return kFALSE;

    fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj(AddSerialNumber("MCerPhotEvt"));
    if (!fCerPhotEvt)
        return kFALSE;
    
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Check for validity of the selected calibration method, switch to a 
// different one in case of need
//
// fill the MPedPhotCam container using the information from MPedestalCam,
// MExtractedSignalCam and MCalibrationCam
//
//
Bool_t MCalibrateData::ReInit(MParList *pList)
{

  if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
    {
      *fLog << warn << GetDescriptor() << "Warning: Blind pixel calibration method not valid, switching to F-factor method" << endl;
      fCalibrationMode = kFfactor;
    }

  if(fCalibrationMode == kPinDiode && !fQEs->IsPINDiodeMethodValid())
    { 
      *fLog << warn << GetDescriptor() << "Warning: PIN diode calibration method not valid, switching to F-factor method" << endl;
      fCalibrationMode = kFfactor;
    }

  //---------------------------------------------
  // fill MPedPhot container using the informations from
  // MPedestalCam, MExtractedSignalCam and MCalibrationCam

  fNumUsedHiGainFADCSlices = fSignals->GetNumUsedHiGainFADCSlices();

  // is pixid equal to pixidx ?
  if (    (Int_t)(fPedestal->GetSize()) != fSignals->GetSize())
  {
    *fLog << err << "MCalibrateData::ReInit(); sizes of MPedestalCam and MCalibrationCam are different" 
          << endl;
  } 

  *fLog << all << "MCalibrateData::ReInit(); fill MPedPhotCam container"
        << endl;
  *fLog << all << "     fNumUsedHiGainADCSlices = " 
        <<  fNumUsedHiGainFADCSlices << endl;
  *fLog << all << "     pixid, calibrationConversionFactor, ped, pedRMS, pedphot, pedphotRMS :"
        << endl;
  for (Int_t pixid=0; pixid<fPedestal->GetSize(); pixid++)
  {
    const MPedestalPix    &ped = (*fPedestal)[pixid];

    // pedestals/(used FADC slices)   in [ADC] counts
    Float_t pedes  = ped.GetPedestal()    * fNumUsedHiGainFADCSlices;
    Float_t pedrms = ped.GetPedestalRms() * sqrt(fNumUsedHiGainFADCSlices);

    //----------------------------------
    // get phe/ADC conversion factor

    Float_t hiloconv;
    Float_t hiloconverr;
    Float_t calibConv;
    Float_t calibFFactor;

    if ( !GetConversionFactor(pixid, hiloconv, hiloconverr,
                              calibConv, calibFFactor ))
      continue;

    //---------------------------------- 

    // pedestals/(used FADC slices)   in [number of photons] 
    Float_t pedphot    = pedes  * calibConv;
    Float_t pedphotrms = pedrms * calibConv;

    (*fPedPhot)[pixid].Set(pedphot, pedphotrms);

    *fLog << all << pixid << ",  " << calibConv << ",  "
          << ped.GetPedestal() << ",  " << ped.GetPedestalRms() << ",  " 
          << pedphot << ",  " << pedphotrms << endl;
  }

  //---------------------------------------------

  fPedPhot->SetReadyToSave();

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Get conversion factor and its error from MCalibrationCam
// 
//
Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx,
 Float_t &hiloconv, Float_t &hiloconverr,
 Float_t &calibConv, Float_t &calibFFactor)
{

  //
  // For the moment, we use only a dummy zenith for the calibration:
  //
  const Float_t zenith = -1.;

  hiloconv     = 1.;
  hiloconverr  = 0.;
  calibConv    = 1.;
  calibFFactor = 0.;

  if(fCalibrationMode!=kNone)
  {
    MCalibrationChargePix &pix   = (MCalibrationChargePix&)(*fCalibrations)[pixidx];       
    MCalibrationQEPix     &qepix = (MCalibrationQEPix&)    (*fQEs)         [pixidx];       
    MBadPixelsPix         &bad   = (*fBadPixels)[pixidx];
    
    if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
      return kFALSE;
    
    hiloconv   = pix.GetConversionHiLo();
    hiloconverr= pix.GetConversionHiLoErr();

    const Float_t pheConv = pix.GetMeanConvFADC2Phe();
    calibFFactor          = pix.GetMeanFFactorFADC2Phot();
    Float_t calibQE       = 1.;
    
    switch(fCalibrationMode)
    {
      case kBlindPixel:
        calibQE      = qepix.GetQECascadesBlindPixel   ( zenith );
        //        calibQEVar   = qepix.GetQECascadesBlindPixelVar( zenith );
      break;
      case kFfactor:
        calibQE      = qepix.GetQECascadesFFactor   ( zenith );
        //        calibQEVar   = qepix.GetQECascadesFFactorVar( zenith );
      break;
      default:
      *fLog << warn << "MCalibrateData::GetConversionFactor; Warning: Calibration mode value ("<<fCalibrationMode<<") not known" << endl;
              break;
    }
    calibConv = pheConv / calibQE;

  }

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Apply the calibration factors to the extracted signal according to the 
// selected calibration method
//
Int_t MCalibrateData::Process()
{
  /*
    if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
    {
        // FIXME: MExtractedSignal must be of variable size -
        //        like MCerPhotEvt - because we must be able
        //        to reduce size by zero supression
        //        For the moment this check could be done in ReInit...
        *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
        return kFALSE;
    }
  */

  UInt_t npix = fSignals->GetSize();

  Float_t hiloconv;
  Float_t hiloconverr;
  Float_t calibrationConversionFactor;
  Float_t calibrationConversionFactorErr;
  
  for (UInt_t pixidx=0; pixidx<npix; pixidx++)
    {
      if ( !GetConversionFactor(pixidx, hiloconv, hiloconverr,
            calibrationConversionFactor, calibrationConversionFactorErr) )
        continue;
      
      MExtractedSignalPix &sig =  (*fSignals)[pixidx];
      
      Float_t signal;
      Float_t signalErr = 0.;
      Float_t nphot,nphotErr;
            
      if (sig.IsLoGainUsed())
        {
          signal    = sig.GetExtractedSignalLoGain()*hiloconv;
          signalErr = signal*hiloconverr;
        }
      else
        {
          if (sig.GetExtractedSignalHiGain() > 9999.)
            {
              signal = 0.;
              signalErr = 0.;
            }
          else
            signal = sig.GetExtractedSignalHiGain();
        }
      
      nphot    = signal*calibrationConversionFactor;
      nphotErr = signal*calibrationConversionFactorErr
        *signal*calibrationConversionFactorErr
        +signalErr*calibrationConversionFactor
        *signalErr*calibrationConversionFactor;
      
      nphotErr  = TMath::Sqrt(nphotErr);
      
      MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);

      if (sig.GetNumHiGainSaturated() > 0)
        cpix->SetPixelHGSaturated();

      if (sig.GetNumLoGainSaturated() > 0)
        cpix->SetPixelSaturated();
    }
  
  fCerPhotEvt->FixSize();
  fCerPhotEvt->SetReadyToSave();
  
  return kTRUE;
}