/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Thomas Bretz    08/2004 <mailto:tbretz@astro.uni-wuerzburg.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!
//
//
//   This class can calibrate data and/or pedestals. To switch off calibration of data
//   set the Calibration Mode to kSkip. To switch on pedestal calibration call either
//     SetPedestalFlag(MCalibrateData::kRun)    (calibration is done once in ReInit)
//     SetPedestalFlag(MCalibrateData::kEvent)  (calibration is done for each event)
//
//   By calling SetNamePedADCContainer() and/or SetNamePedPhotContainer() you
//   can control the name of the MPedestalCam and/or MPedPhotCam container which is used.
//
//   Assume you want to calibrate "MPedestalCam" once and "MPedestalFromLoGain" [MPedestalCam]
//   event-by-event, so:
//     MCalibrateData cal1;
//     cal1.SetCalibrationMode(MCalibrateData::kSkip);
//     cal1.SetPedestalFlag(MCalibrateData::kRun);
//     MCalibrateData cal2;
//     cal2.SetCalibrationMode(MCalibrateData::kSkip);
//     cal2.SetNamePedADCContainer("MPedestalFromLoGain");
//     cal2.SetNamePedPhotContainer("MPedPhotFromLoGain")
//     cal2.SetPedestalFlag(MCalibrateData::kEvent);
//
//
//   Input Containers:
//    [MPedestalCam]
//    [MExtractedSignalCam]
//    [MCalibrationChargeCam]
//    [MCalibrationQECam]
//    MBadPixelsCam
//
//   Output Containers:
//    [MPedPhotCam]
//    [MCerPhotEvt]
//
// See also: MJCalibration, MJPedestal, MJExtractSignal, MJExtractCalibTest
// 
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrateData.h"

#include <fstream>

#include <TEnv.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;

const TString MCalibrateData::fgNamePedADCContainer  = "MPedestalCam";
const TString MCalibrateData::fgNamePedPhotContainer = "MPedPhotCam";
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets all pointers to NULL
// 
// Initializes:
// - fCalibrationMode to kDefault
// - fPedestalFlag to kRun
// - fNamePedADCRunContainer    to "MPedestalCam"  
// - fNamePedPhotRunContainer   to "MPedPhotCam"    
//
MCalibrateData::MCalibrateData(CalibrationMode_t calmode,const char *name, const char *title) 
    : fGeomCam(NULL),   fPedestal(NULL),
      fBadPixels(NULL), fCalibrations(NULL), fQEs(NULL), fSignals(NULL), 
      fPedPhot(NULL),   fCerPhotEvt(NULL), fPedestalFlag(kNo)
{

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

  SetNamePedADCContainer();
  SetNamePedPhotContainer();
}

// --------------------------------------------------------------------------
//
// 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

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

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

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

    fCalibrations = 0;
    fQEs          = 0;
    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;
        }
    }

    fPedestal = 0;
    fPedPhot  = 0;
    if (fPedestalFlag)
    {
        fPedestal = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedADCContainer), "MPedestalCam");
        if (!fPedestal)
        {
            *fLog << err << AddSerialNumber(fNamePedADCContainer) << " [MPedestalCam] not found ... aborting" << endl;
            return kFALSE;
        }

        fPedPhot = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam", AddSerialNumber(fNamePedPhotContainer));
        if (!fPedPhot)
            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;
    }

    // Sizes might have changed
    if (fPedestalFlag && (Int_t)fPedestal->GetSize() != fSignals->GetSize())
    {
        *fLog << err << "Size mismatch of MPedestalCam and MCalibrationCam... abort." << endl;
        return kFALSE;
    }

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

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

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

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

    if (TestPedestalFlag(kRun))
        Calibrate(kFALSE, kTRUE);

    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) const
{
    //
    // 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.;

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

        hiloconv   = pix.GetConversionHiLo   ();
        hiloconverr= pix.GetConversionHiLoErr();

        if ((*fBadPixels)[pixidx].IsUnsuitable())
            return kFALSE;

        calibConv    = pix.GetMeanConvFADC2Phe();
        calibConvVar = pix.GetMeanConvFADC2PheVar();
        calibFFactor = pix.GetMeanFFactorFADC2Phot();

        MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];

        switch(fCalibrationMode)
        {
        case kFlatCharge:
            {
                MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
                calibConv    = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx));
                calibConvVar = (avpix.GetMeanRelVar() + 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())
                return kFALSE;
            calibQE     = qe.GetQECascadesBlindPixel   ( zenith );
            calibQEVar  = qe.GetQECascadesBlindPixelVar( zenith );
            break;

        case kPinDiode:
            if (!qe.IsPINDiodeMethodValid())
                return kFALSE;
            calibQE     = qe.GetQECascadesPINDiode   ( zenith );
            calibQEVar  = qe.GetQECascadesPINDiodeVar( zenith );
            break;

        case kFfactor:
            if (!pix.IsFFactorMethodValid())
                return kFALSE;
            calibQE     = qe.GetQECascadesFFactor   ( zenith );
            calibQEVar  = qe.GetQECascadesFFactorVar( zenith );
            break;

        case kCombined:
            if (!qe.IsCombinedMethodValid())
                return kFALSE;
            calibQE     = qe.GetQECascadesCombined   ( zenith );
            calibQEVar  = qe.GetQECascadesCombinedVar( zenith );
            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.)
    {
        // Now doing:
        // calibConvVar  = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
        // calibConvVar *= (calibConv*calibConv);
        calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE);
    }

    return kTRUE;
}

Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
{
    if (!data && !pedestal)
        return kTRUE;

    const UInt_t  npix       = fSignals->GetSize();
    const Float_t slices     = fSignals->GetNumUsedHiGainFADCSlices();
    const Float_t sqrtslices = TMath::Sqrt(slices);

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

    UInt_t skip = 0;
    for (UInt_t pixidx=0; pixidx<npix; pixidx++)
    {
        if (!GetConversionFactor(pixidx, hiloconv, hiloconverr,
                                 calibConv, calibConvErr, calibFFactor))
        {
            skip++;
            continue;
        }

        if (data)
        {
            const MExtractedSignalPix &sig = (*fSignals)[pixidx];

            Float_t signal = 0;
            Float_t signalErr = 0.;

            if (sig.IsLoGainUsed())
            {
                signal    = sig.GetExtractedSignalLoGain()*hiloconv;
                signalErr = signal*hiloconverr;
            }
            else
            {
                if (sig.GetExtractedSignalHiGain() <= 9999.)
                    signal = sig.GetExtractedSignalHiGain();
            }

            const Float_t nphot    = signal*calibConv;
            const Float_t 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();
        }
        if (pedestal)
        {
            const MPedestalPix &ped = (*fPedestal)[pixidx];

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

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

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

    if (skip>npix*0.9)
    {
        *fLog << warn << "WARNING - GetConversionFactor has skipped more than 90% of the pixels... skip." << endl;
        return kCONTINUE;
    }

    if (pedestal)
        fPedPhot->SetReadyToSave();

    if (data)
    {
        fCerPhotEvt->FixSize();
        fCerPhotEvt->SetReadyToSave();
    }
    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;
     }
     */

    return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent));
}

// --------------------------------------------------------------------------
//
// 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 (TestPedestalFlag(kEvent))
        out << "   " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kEvent)" << endl;
    if (TestPedestalFlag(kRun))
        out << "   " << GetUniqueName() << ".EnablePedestalType(MCalibrateData::kRun)" << endl;

    if (fCalibrationMode != kDefault)
    {
        out << "   " << GetUniqueName() << ".SetCalibrationMode(MCalibrateData::";
        switch (fCalibrationMode)
        {
        case kSkip:       out << "kSkip";               break;
        case kNone:       out << "kNone";               break;
        case kFlatCharge: out << "kFlatCharge";         break;
        case kBlindPixel: out << "kBlindPixel";         break;
        case kFfactor:    out << "kFfactor";            break;
        case kPinDiode:   out << "kPinDiode";           break;
        case kCombined:   out << "kCombined";           break;
        case kDummy:      out << "kDummy";              break;
        default:          out << (int)fCalibrationMode; break;
        }
        out << ")" << endl;
    }

    if (fNamePedADCContainer != fgNamePedADCContainer)
    {
        out << "   " << GetUniqueName() << ".SetNamePedADCContainer(";
        out << fNamePedADCContainer.Data() << ");" << endl;
    }

    if (fNamePedPhotContainer  != fgNamePedPhotContainer)
    {
        out << "   " << GetUniqueName() << ".SetNamePedPhotContainer(";
        out << fNamePedPhotContainer.Data() << ");" << endl;
    }
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MJPedestal.MCalibrateDate.PedestalFlag: no,run,event
//   MJPedestal.MCalibrateDate.CalibrationMode: skip,none,flatcharge,blindpixel,ffactor,pindiode,combined,dummy,default
//
Int_t MCalibrateData::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "PedestalFlag", print))
    {
        rc = kTRUE;
        TString s = GetEnvValue(env, prefix, "PedestalFlag", "");
        s.ToLower();
        if (s.BeginsWith("no"))
            SetPedestalFlag(kNo);
        if (s.BeginsWith("run"))
            SetPedestalFlag(kRun);
        if (s.BeginsWith("event"))
            SetPedestalFlag(kEvent);
    }

    if (IsEnvDefined(env, prefix, "CalibrationMode", print))
    {
        rc = kTRUE;
        TString s = GetEnvValue(env, prefix, "CalibrationMode", "");
        s.ToLower();
        if (s.BeginsWith("skip"))
            SetCalibrationMode(kSkip);
        if (s.BeginsWith("none"))
            SetCalibrationMode(kNone);
        if (s.BeginsWith("flatcharge"))
            SetCalibrationMode(kFlatCharge);
        if (s.BeginsWith("blindpixel"))
            SetCalibrationMode(kBlindPixel);
        if (s.BeginsWith("ffactor"))
            SetCalibrationMode(kFfactor);
        if (s.BeginsWith("pindiode"))
            SetCalibrationMode(kPinDiode);
        if (s.BeginsWith("combined"))
            SetCalibrationMode(kCombined);
        if (s.BeginsWith("dummy"))
            SetCalibrationMode(kDummy);
        if (s.BeginsWith("default"))
            SetCalibrationMode();
    }

    return rc;
}