/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Hendrik Bartko  08/2004 <mailto:hbartko@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

///////////////////////////////////////////////////////////////////////////////////
//
//   MCalibrateData
//
//   This task takes the integrated charge from MExtractedSignal and applies
//   the calibration constants from MCalibrationCam to convert the summed FADC 
//   slices into photons. The number of photons obtained is stored in MCerPhotEvt. 
//   Optionally, the calibration of pedestals from an MPedestalCam container into 
//   an MPedPhotCam container can be chosen with the member functions 
//   SetPedestalType(). Default is 'kRun', i.e. calibration of pedestals from a 
//   dedicated pedestal run.
//   In case, the chosen pedestal type is kRun or kEvt, in ReInit() the MPedPhotCam 
//   container is filled using the information from MPedestalCam, MExtractedSignalCam, 
//   MCalibrationChargeCam and MCalibrationQECam
//
//   Selection of different calibration methods is allowed through the 
//   SetCalibrationMode() member function (default: kFfactor)
//
//   The calibration modes which exclude non-valid pixels are the following: 
//
//   kFfactor:    calibrates using the F-Factor method
//   kBlindpixel: calibrates using the BlindPixel method 
//   kBlindpixel: calibrates using the BlindPixel method 
//   kFlatCharge: perform a charge flat-flatfielding. Outer pixels are area-corrected.
//   kDummy:      calibrates with fixed conversion factors of 1 and errors of 0.
//
//   The calibration modes which include all pixels regardless of their validity is:
//
//   kNone:       calibrates with fixed conversion factors of 1 and errors of 0.
//
//   Use the kDummy and kNone methods ONLY FOR DEBUGGING!   
//
//   Input Containers:
//    MPedestalCam
//    MExtractedSignalCam
//    MCalibrationChargeCam
//    MCalibrationQECam
//
//   Output Containers:
//    MPedPhotCam
//    MCerPhotEvt
//
// See also: MJCalibration, MJPedestal, MJExtractSignal, MJExtractCalibTest
// 
//////////////////////////////////////////////////////////////////////////////
#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) 
    : fGeomCam(NULL), fPedestal(NULL), fBadPixels(NULL), fCalibrations(NULL), fSignals(NULL), 
      fPedPhot(NULL), fCerPhotEvt(NULL), fCalibrationMode(calmode), fFlags(kRun), fNamePedADCRunContainer("MPedestalCam"), fNamePedADCEventContainer("MPedestalCamFromData"), fNamePedPhotRunContainer("MPedPhotCam"), fNamePedPhotEventContainer("MPedPhotCamFromData")
{
    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)
{

  // input containers


  if (TestFlag(kRun))
  {
    fPedestal = (MPedestalCam*)pList->FindObject(fNamePedADCRunContainer, AddSerialNumber("MPedestalCam"));
    if (!fPedestal)
    {  
      *fLog << err << "container 'MPedestalCam'" 
	    << " not found ... aborting" << endl;
      return kFALSE;
    }
  }


  if (TestFlag(kEvent))
  {
    fPedestalFromData = (MPedestalCam*)pList->FindObject(fNamePedADCEventContainer, AddSerialNumber("MPedestalCam"));
    if (!fPedestalFromData)
    {
      *fLog << err << "container 'MPedestalCamFromData'" 
	    << " 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;
        }
    }

    // output containers

    if (TestFlag(kRun))
    {
      fPedPhot = (MPedPhotCam*)pList->FindCreateObj(AddSerialNumber("MPedPhotCam"), fNamePedPhotRunContainer);
      if (!fPedPhot)
      {
	*fLog << err << "container 'MPedPhotCam'" 
	      << " not found ... aborting" << endl;
        return kFALSE;
      }
    }


    if (TestFlag(kEvent))
    {  
      fPedPhotFromData = (MPedPhotCam*)pList->FindCreateObj(AddSerialNumber("MPedPhotCam"), fNamePedPhotEventContainer);
      if (!fPedPhotFromData)
      {
	*fLog << err << "container 'MPedPhotCamFromData'" 
	      << " not found ... aborting" << endl;
        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;
    }


  if (TestFlag(kRun))
  {
    if (!CalibratePedestal(fPedestal,fPedPhot)) return kFALSE;
  }

  return kTRUE;
}



// --------------------------------------------------------------------------
//
// Calibrate the Pedestal values
// 
//
Bool_t MCalibrateData::CalibratePedestal(MPedestalCam *pedADCCam, MPedPhotCam *pedPhotCam)
{
  //---------------------------------------------
  // fill MPedPhot container using the informations from
  // MPedestalCam, MExtractedSignalCam and MCalibrationCam

  fNumUsedHiGainFADCSlices = fSignals->GetNumUsedHiGainFADCSlices();

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

 

  for (Int_t pixid=0; pixid<pedADCCam->GetSize(); pixid++)
  {
    const MPedestalPix    &ped = (*pedADCCam)[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 calibConvVar;
    Float_t calibFFactor;

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

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

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

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

  }

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

  pedPhotCam->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 &calibConvVar, 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.;
  calibConvVar  = 0.;
  calibFFactor = 0.;
  Float_t calibQE       = 1.;
  Float_t calibQEVar    = 0.;

  Float_t avMean        = 1.;
  Float_t avMeanRelVar  = 0.;



  if (fCalibrationMode == kFlatCharge)
  {
    MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
    avMean       =  avpix.GetMean();
    avMeanRelVar =  avpix.GetMeanRelVar();
  }


  if(fCalibrationMode!=kNone)
  {
	  
    MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];

    hiloconv   = pix.GetConversionHiLo   ();
    hiloconverr= pix.GetConversionHiLoErr();
	  
    if (fBadPixels)
    {
      MBadPixelsPix         &bad = (*fBadPixels)[pixidx];
      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
	return kFALSE;
    }
            
    calibConv      = pix.GetMeanConvFADC2Phe();
    calibConvVar   = pix.GetMeanConvFADC2PheVar();
    calibFFactor   = pix.GetMeanFFactorFADC2Phot();

    MCalibrationQEPix &qe  = (MCalibrationQEPix&) (*fQEs)[pixidx];
          
    switch(fCalibrationMode)
    {
    case kFlatCharge:
      calibConv        = avMean / pix.GetMean() / fGeomCam->GetPixRatio(pixidx) ;
      calibConvVar     = (avMeanRelVar + pix.GetMeanRelVar()) * calibConv * calibConv;
      if (pix.IsFFactorMethodValid())
      {
	const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
	if (convmin1 > 0)
	  calibFFactor *= TMath::Sqrt(convmin1);
	else
	  calibFFactor = -1.;
      }
      break;
    case kBlindPixel:
      if (qe.IsBlindPixelMethodValid())
      {
	calibQE      = qe.GetQECascadesBlindPixel   ( zenith );
	calibQEVar   = qe.GetQECascadesBlindPixelVar( zenith );
      }
      else
	return kFALSE;
      break;
    case kPinDiode:
      if (qe.IsPINDiodeMethodValid())
      {
	calibQE      = qe.GetQECascadesPINDiode   ( zenith );
	calibQEVar   = qe.GetQECascadesPINDiodeVar( zenith );
      }
      else
	return kFALSE;
      break;
    case kFfactor:
      if (pix.IsFFactorMethodValid())
      {
	calibQE      = qe.GetQECascadesFFactor   ( zenith );
	calibQEVar   = qe.GetQECascadesFFactorVar( zenith );
      }
      else
	return kFALSE;
      break;
    case kCombined:
      if (qe.IsCombinedMethodValid())
      {
	calibQE      = qe.GetQECascadesCombined   ( zenith );
	calibQEVar   = qe.GetQECascadesCombinedVar( zenith );
      }
      else
	return kFALSE;
      break;
    case kDummy:
      hiloconv     = 1.;
      hiloconverr  = 0.;
      calibQE      = 1.;
      calibQEVar   = 0.;
      break;
              
    } /* switch calibration mode */
  } /* if(fCalibrationMode!=kNone) */
  else
  {
    hiloconv       = 1.;
    hiloconverr    = 0.;
    calibConv      = 1./fGeomCam->GetPixRatio(pixidx);
    calibConvVar   = 0.;
    calibFFactor   = 0.;
    calibQE        = 1.;
    calibQEVar     = 0.;
  }     
  
  calibConv    /= calibQE;
  calibConvVar /= 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;
  Float_t calibFFactor;

  for (UInt_t pixidx=0; pixidx<npix; pixidx++)
    {
      if ( !GetConversionFactor(pixidx, hiloconv, hiloconverr,
	    calibrationConversionFactor, calibrationConversionFactorErr, calibFFactor) )
        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 = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
      //   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();
  


  if(TestFlag(kEvent)) 
  {
    if (!CalibratePedestal(fPedestalFromData,fPedPhotFromData)) return kFALSE;
  }

  /*
  //---------------------------------------------
  // fill MPedPhot(FromData) 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<fPedestalFromData->GetSize(); pixid++)
  {
    const MPedestalPix    &ped = (*fPedestalFromData)[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;

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

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

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

  fPedPhotFromData->SetReadyToSave();

  }
     */

  return kTRUE;
}
