/* ======================================================================== *\
!
! *
! * 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): Markus Gaug  02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationChargeCalc
//
//   Task to calculate the calibration conversion factors from the FADC
//   time slices. The integrated time slices have to be delivered by an 
//   MExtractedSignalCam. The pedestals by an MPedestalCam.
//
//   The output container MCalibrationCam holds one entry of type MCalibrationChargePix 
//   for every pixel. It is filled in the following way:
//
//   ProProcess: Initialize MCalibrationCam
//               Initialize pulser light wavelength
//               
//   ReInit:     MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates
//               memory in a TClonesArray of type MCalibrationChargePix)
//               Initializes pointer to MBadPixelsCam
//
//   Process:    Nothing done by this class, histograms are filled by 
//               MHCalibrationChargeCam
//
//   PstProcess:  Fit results from MHCalibrationChargeCam are retrieved
//                and used for the calculation of the reduced sigma, 
//                the F-Factor method, the blind pixel method (photon flux 
//                                     inside plexiglass) and
//                                     the PINDiode method (photon flux
//                                     outside plexiglass)
//
//                Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
//                can be skipped with the command:
//                MalibrationCam::SkipHiLoGainCalibration()
//
//  Input Containers:
//   MRawEvtData
//   MPedestalCam
//   MBadPixelsCam
//
//  Output Containers:
//   MCalibrationCam
//   MBadPixelsCam
//
//
//  Preliminary description of the calibration in photons (email from 12/02/04)
//
//  Why calibrating in photons:
//  ===========================
//  
//  At the Barcelona meeting in 2002, we decided to calibrate the camera in
//  photons. This for the following reasons:
//  
//  * The physical quantity arriving at the camera are photons. This is
//    the direct physical information from the air shower. The photons
//    have a flux and a spectrum.
//  
//  * The photon fluxes depend mostly on the shower energy (with
//    corrections deriving from the observation conditions), while the photon
//    spectra depend mostly on the observation conditions: zenith angle,
//    quality of the air, also the impact parameter of the shower.
//  
//  * The photomultiplier, in turn, has different response properties
//    (quantum efficiencies) for photons of different colour. (Moreover,
//    different pixels have slightly different quantum efficiencies).
//    The resulting number of photo-electrons is then amplified (linearly)
//    with respect to the photo-electron flux.
//  
//  * In the ideal case, one would like to disentagle the effects
//    of the observation conditions from the primary particle energy (which
//    one likes to measure). To do so, one needs:
//  
//    1) A reliable calibration relating the FADC counts to the photo-electron
//       flux -> This is accomplished with the F-Factor method.
//  
//    2) A reliable calibration of the wavelength-dependent quantum efficiency
//       -> This is accomplished with the combination of the three methods,
//          together with QE-measurements performed by David in order to do
//          the interpolation.
//  
//    3) A reliable calibration of the observation conditions. This means:
//       - Tracing the atmospheric conditions   -> LIDAR
//       - Tracing the observation zenith angle -> Drive System
//   4) Some knowlegde about the impact parameter:
//       - This is the only part which cannot be accomplished well with a
//         single telescope. We would thus need to convolute the spectrum
//         over the distribution of impact parameters.
//  
//  
//  How an ideal calibration would look like:
//  =========================================
//  
//  We know from the combined PIN-Diode and Blind-Pixel Method the response of
//  each pixel to well-measured light fluxes in three representative
//  wavelengths (green, blue, UV). We also know the response to these light
//  fluxes in photo-electrons. Thus, we can derive:
//  
//  - conversion factors to photo-electrons
//  - conversion factors to photons in three wavelengths.
//  
//  Together with David's measurements and some MC-simulation, we should be
//  able to derive tables for typical Cherenkov-photon spectra - convoluted
//  with the impact parameters and depending on the athmospheric conditions
//  and the zenith angle (the "outer parameters").
//  
//  From these tables we can create "calibration tables" containing some
//  effective quantum efficiency depending on these outer parameters and which
//  are different for each pixel.
//  
//  In an ideal MCalibrate, one would thus have to convert first the FADC
//  slices to Photo-electrons and then, depending on the outer parameters,
//  look up the effective quantum efficiency and get the mean number of
//  photons which is then used for the further analysis.
//  
//  How the (first) MAGIC calibration should look like:
//  ===================================================
//  
//  For the moment, we have only one reliable calibration method, although
//  with very large systematic errors. This is the F-Factor method. Knowing
//  that the light is uniform over the whole camera (which I would not at all
//  guarantee in the case of the CT1 pulser), one could in principle already
//  perform a relative calibration of the quantum efficiencies in the UV.
//  However, the spread in QE at UV is about 10-15% (according to the plot
//  that Abelardo sent around last time. The spread in photo-electrons is 15%
//  for the inner pixels, but much larger (40%) for the outer ones.
//  
//  I'm not sure if we can already say that we have measured the relative
//  difference in quantum efficiency for the inner pixels and produce a first
//  QE-table for each pixel. To so, I would rather check in other wavelengths
//  (which we can do in about one-two weeks when the optical transmission of
//  the calibration trigger is installed).
//  
//  Thus, for the moment being, I would join Thomas proposal to calibrate in
//  photo-electrons and apply one stupid average quantum efficiency for all
//  pixels. This keeping in mind that we will have much preciser information
//  in about one to two weeks.
//  
//  
//  What MCalibrate should calculate and what should be stored:
//  ===========================================================
//  
//  It is clear that in the end, MCerPhotEvt will store photons.
//  MCalibrationCam stores the conversionfactors to photo-electrons and also
//  some tables of how to apply the conversion to photons, given the outer
//  parameters. This is not yet implemented and not even discussed.
//  
//  To start, I would suggest that we define the "average quantum efficiency"
//  (maybe something like 25+-3%) and apply them equally to all
//  photo-electrons. Later, this average factor can be easily replaced by a
//  pixel-dependent factor and later by a (pixel-dependent) table.
//  
//  
//  
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeCalc.h"

#include <TSystem.h>
#include <TH1.h>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MRawRunHeader.h"
#include "MRawEvtPixelIter.h"

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

#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"
#include "MCalibrationChargePINDiode.h"
#include "MCalibrationChargeBlindPix.h"

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

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


ClassImp(MCalibrationChargeCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
    : fPedestals(NULL), fCam(NULL), 
      fRawEvt(NULL), fRunHeader(NULL), fGeom(NULL), 
      fBadPixels(NULL), fEvtTime(NULL),
      fSignals(NULL), fPINDiode(NULL), fBlindPixel(NULL)
{

    fName  = name  ? name  : "MCalibrationChargeCalc";
    fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";

    AddToBranchList("MRawEvtData.fHiGainPixId");
    AddToBranchList("MRawEvtData.fLoGainPixId");
    AddToBranchList("MRawEvtData.fHiGainFadcSamples");
    AddToBranchList("MRawEvtData.fLoGainFadcSamples");

    Clear();
}

void MCalibrationChargeCalc::Clear(const Option_t *o)
{
  
    SETBIT(fFlags, kUseQualityChecks);
    SETBIT(fFlags, kHiLoGainCalibration);

    fNumHiGainSamples  = 0;
    fNumLoGainSamples  = 0;
    fConversionHiLo    = 0;

}


// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MCalibrationCam
//
// The following output containers are only searched, but not created
//
//  - MTime
//
Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
{

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

    fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
    if (!fCam)
        return kFALSE;

    fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
    if (!fPINDiode)
        return kFALSE;

    fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
    if (!fBlindPixel)
        return kFALSE;

    fCam->SetPINDiode(fPINDiode);
    fCam->SetBlindPixel(fBlindPixel);

    fEvtTime = (MTime*)pList->FindObject("MTime");

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

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

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//
Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
{

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

  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeom)
    {
      *fLog << err << "No MGeomCam found... aborting." << endl;
      return kFALSE;
    }
  
  fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
  if (!fBadPixels)
    {
      *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
      return kFALSE;
    }
  
  fCam->SetGeomCam(fGeom);

  fNumHiGainSamples  =  fSignals->GetNumUsedHiGainFADCSlices();
  fNumLoGainSamples  =  fSignals->GetNumUsedLoGainFADCSlices();
  fSqrtHiGainSamples =  TMath::Sqrt((Float_t)fNumHiGainSamples);
  
  UInt_t npixels     = fGeom->GetNumPixels();
  
  for (UInt_t i=0; i<npixels; i++)
    {
      MCalibrationChargePix &pix = (*fCam)[i];
      MBadPixelsPix         &bad = (*fBadPixels)[i];

      pix.DefinePixId(i);

      if (bad.IsBad())
        {
          pix.SetExcluded();
          continue;
        }
      
      pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
                                  fSignals->GetLastUsedSliceHiGain());
      pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
                                  fSignals->GetLastUsedSliceLoGain());
    }

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MCerPhotEvt container.
//
Int_t MCalibrationChargeCalc::Process()
{
  return kTRUE;
}

Int_t MCalibrationChargeCalc::PostProcess()
{

  //
  // loop over the pedestal events and check if we have calibration
  //
  Int_t nvalid   = 0;
  Float_t avinnerped  = 0;
  Float_t avinnerprms = 0;
  Float_t avinnernum  = 0;
  Float_t avouterped  = 0;
  Float_t avouterprms = 0;
  Float_t avouternum  = 0;

  for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
    {

      MCalibrationChargePix &pix = (*fCam)[pixid];
      MBadPixelsPix         &bad = (*fBadPixels)[pixid];

      //
      // Check if the pixel has been excluded from the fits
      //
      if (pix.IsExcluded())
        continue;

      //
      // get the pedestals
      //
      const Float_t ped    = (*fPedestals)[pixid].GetPedestal();
      const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms();
      const Float_t num    = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());

      if (fGeom->GetPixRatio(pixid) == 1.)
      {
	  avinnerped  += ped;
	  avinnerprms += prms;
	  avinnernum++;
      }
      else
      {
	  avouterped  += ped;
	  avouterprms += prms;
	  avouternum++;
      }
      //
      // set them in the calibration camera
      //
      if (pix.IsHiGainSaturation())
      {
	  pix.SetPedestal(ped  * fNumLoGainSamples,
			  prms * TMath::Sqrt((Float_t)fNumLoGainSamples),
			  prms * fNumLoGainSamples / num);
	  pix.SetNumLoGainSamples((Float_t)fNumLoGainSamples);
          pix.ApplyLoGainConversion();
      }
      else
      {
	  pix.SetPedestal(ped  * fNumHiGainSamples,
			  prms * TMath::Sqrt((Float_t)fNumHiGainSamples),
			  prms * fNumHiGainSamples / num);
      }

      pix.CheckChargeValidity(&bad);
      pix.CheckTimeValidity(&bad);
      
      if (!bad.IsCalibrationOK())
        continue;

      nvalid++;

      if (!pix.CalcReducedSigma())
	  continue;

      pix.CalcFFactorMethod();

    }



  //
  // The Michele check ...
  //
  if (nvalid == 0)
  {
      *fLog << err << GetDescriptor() << ": Dear Michele! All pixels have non-valid calibration. " 
	    << "Did you forget to fill the histograms (filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
    return kFALSE;
  }

  MCalibrationChargePix *avinnerpix = fCam->GetAverageInnerPix();
  MCalibrationChargePix *avouterpix = fCam->GetAverageOuterPix();
  //
  // set the pedestans in the calibration camera
  //
  if (avinnerpix->IsHiGainSaturation())
  {
      avinnerpix->SetPedestal(avinnerped/avinnernum  * fNumLoGainSamples,
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples),
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumLoGainSamples/avinnernum));
      avinnerpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
      avinnerpix->ApplyLoGainConversion();
  }
  else
  {
      avinnerpix->SetPedestal(avinnerped/avinnernum  * fNumHiGainSamples,
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples),
			 avinnerprms/avinnernum * TMath::Sqrt((Float_t)fNumHiGainSamples/avinnernum));
  }

  if (avouterpix->IsHiGainSaturation())
  {
      avouterpix->SetPedestal(avouterped/avouternum  * fNumLoGainSamples,
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples),
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumLoGainSamples/avouternum));
      avouterpix->SetNumLoGainSamples((Float_t)fNumLoGainSamples);
      avouterpix->ApplyLoGainConversion();
  }
  else
  {
      avouterpix->SetPedestal(avouterped/avouternum  * fNumHiGainSamples,
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples),
			      avouterprms/avouternum * TMath::Sqrt((Float_t)fNumHiGainSamples/avouternum));
  }
  
  if (!avinnerpix->CheckChargeValidity() || !avinnerpix->CheckTimeValidity())
      if (!avinnerpix->CalcReducedSigma())
	  avinnerpix->CalcFFactorMethod();

  if (!avouterpix->CheckChargeValidity() || !avouterpix->CheckTimeValidity())
      if (!avouterpix->CalcReducedSigma())
	  avouterpix->CalcFFactorMethod();


  if (!fBlindPixel->CheckChargeFitValidity())
  {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
      fCam->SetBlindPixelMethodValid(kFALSE);
  }
  else
  { 
      if (!fBlindPixel->CalcFluxInsidePlexiglass())
      {
	  *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
	  fCam->SetBlindPixelMethodValid(kFALSE);
      }
      else
      {
	  fCam->SetBlindPixelMethodValid(kTRUE);
	  fCam->ApplyBlindPixelCalibration();
      }
  }

  if (!fPINDiode->CheckChargeFitValidity() || !fPINDiode->CheckTimeFitValidity())
  {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, charge fit not valid " << endl;
      fCam->SetPINDiodeMethodValid(kFALSE);
  }
  else
  { 
      if (!fPINDiode->CalcFluxOutsidePlexiglass())
      {
	  *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, will skip PIN Diode Calibration " << endl;
	  fCam->SetPINDiodeMethodValid(kFALSE);
      }
      else
      {
	  fCam->SetPINDiodeMethodValid(kTRUE);
	  fCam->ApplyPINDiodeCalibration();
      }
  }
  fCam->SetReadyToSave();
  
return kTRUE;
}
