/* ======================================================================== *\
!
! *
! * 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 kEvent, 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 <fstream>

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

const TString MCalibrateData::fgNamePedADCRunContainer     = "MPedestalCam"        ;   
const TString MCalibrateData::fgNamePedADCEventContainer   = "MPedestalCamFromData"; 
const TString MCalibrateData::fgNamePedPhotRunContainer    = "MPedPhotCam"         ;  
const TString MCalibrateData::fgNamePedPhotEventContainer  = "MPedPhotCamFromData" ;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets all pointers to NULL
// 
// Initializes:
// - fCalibrationMode to kDefault
// - fPedestalFlag to kRun
// - fNamePedADCRunContainer    to "MPedestalCam"  
// - fNamePedADCEventContainer  to "MPedestalCamFromData" 
// - fNamePedPhotRunContainer   to "MPedPhotCam"    
// - fNamePedPhotEventContainer to "MPedPhotCamFromData"
//
MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title) 
    : fGeomCam(NULL), fPedestal(NULL), fPedestalFromData(NULL), 
      fBadPixels(NULL), fCalibrations(NULL), fQEs(NULL), fSignals(NULL), 
      fPedPhot(NULL), fPedPhotFromData(NULL), fCerPhotEvt(NULL)
{

  fName  = name  ? name  : "MCalibrateData";
  fTitle = title ? title : "Task to calculate the number of photons in one event";
  
  SetCalibrationMode();
  SetPedestalType();

  SetNamePedADCRunContainer   ();
  SetNamePedADCEventContainer ();
  SetNamePedPhotRunContainer  ();
  SetNamePedPhotEventContainer();

}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//
//  - MGeomCam
//  - MPedestalCam
//  - MCalibrationChargeCam
//  - MCalibrationQECam
//  - MExtractedSignalCam
//  - MBadPixelsCam
//
// 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;
}

// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//
//  - MGeomCam
//
// 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)
{

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

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

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

  if(fCalibrationMode == kCombined && !fQEs->IsCombinedMethodValid())
    { 
      *fLog << warn << GetDescriptor() 
            << ": Combined calibration method not valid, switching to default F-factor method" << endl;
      fCalibrationMode = kFfactor;
    }

  //
  // output information or warnings:
  //
  switch(fCalibrationMode)
    {
    case kBlindPixel:
      break;
    case kFfactor:
      break;
    case kPinDiode:
      *fLog << err << GetDescriptor() 
                    << ": PIN Diode Calibration mode not yet available " << endl;
      return kFALSE;
      break;
    case kCombined:
      *fLog << err << GetDescriptor() 
                    << ": Combined Calibration mode not yet available " << endl;
      return kFALSE;
      break;
    case kFlatCharge:
      *fLog << warn << GetDescriptor() 
            << ": WARNING: Flat-fielding charges - only for muon calibration! " << endl;
      break;
    case kDummy:
      *fLog << warn << GetDescriptor() 
            << ": WARNING: Dummy calibration, no calibration applied!!" << endl;
      break;
    case kNone:
      *fLog << warn << GetDescriptor() 
            << ": WARNING: No calibration applied!!" << endl;
      break;
    default:
      *fLog << warn << GetDescriptor() 
            << ": WARNING: Calibration mode value ("
            << fCalibrationMode << ") not known" << endl;
      return kFALSE;
    }

  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
  //
  const Float_t slices = 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()    * slices;
    Float_t pedrms = ped.GetPedestalRms() * TMath::Sqrt(slices);

    //
    // 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())
            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.;
          break;
          
        } /* switch calibration mode */
    } /* if(fCalibrationMode!=kNone) */
  else
    {
      calibConv  = 1./fGeomCam->GetPixRatio(pixidx);
    }     
  
  calibConv    /= calibQE;

  if (calibConv != 0. && calibQE != 0.)
    {
      calibConvVar  = calibConvVar/calibConv/calibConv + calibQEVar/calibQE/calibQE;
      calibConvVar *= calibConv*calibConv;
    }
  
  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));
      //
      // The following part is the commented first version of the error calculation
      // Contact Markus Gaug for questions (or wait for the next documentation update...)
      //
      /*
        nphotErr = signal    > 0 ? signalErr*signalErr / (signal * signal)  : 0.
                 + calibConv > 0 ? calibConvVar  / (calibConv * calibConv ) : 0.;
        nphotErr  = TMath::Sqrt(nphotErr) * nphot;
      */

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

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MCalibrateData::StreamPrimitive(ofstream &out) const
{

  out << "   " << ClassName() << " " << GetUniqueName() << "(\"";
  out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;

  if (TestFlag(kEvent))
    out << "   " << GetUniqueName() << ".EnablePedestalType(kEvent)" << endl;
  if (TestFlag(kRun))
    out << "   " << GetUniqueName() << ".EnablePedestalType(kRun)" << endl;

  if (fCalibrationMode != kDefault)
    out << "   " << GetUniqueName() << ".SetCalibrationMode(fCalibrationMode)" << endl;

  if (fNamePedADCRunContainer != fgNamePedADCRunContainer)
    {
      out << "   " << GetUniqueName() << ".SetNamePedADCRunContainer(";
      out << fNamePedADCRunContainer.Data() << ");" << endl;
    }
  
  if (fNamePedADCEventContainer != fgNamePedADCEventContainer)
    {
      out << "   " << GetUniqueName() << ".SetNamePedADCEventContainer(";
      out << fNamePedADCEventContainer.Data() << ");" << endl;
    }

  
  if (fNamePedPhotRunContainer  != fgNamePedPhotRunContainer)
    {
      out << "   " << GetUniqueName() << ".SetNamePedPhotRunContainer(";
      out << fNamePedPhotRunContainer.Data() << ");" << endl;
    }


  if (fNamePedPhotEventContainer != fNamePedPhotEventContainer)
    {
      out << "   " << GetUniqueName() << ".SetNamePedPhotEventContainer(";
      out <<  fNamePedPhotEventContainer.Data() << ");" << endl;
    }

}

  
