/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek 1/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz 4/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHSigmaTheta                                                            //
//                                                                          //
//  calculates - the 2D-histogram   sigmabar vs. Theta (Inner), and         //
//             - the 2D-histogram   sigmabar vs. Theta (Outer)              //
//             - the 3D-histogram   sigma, pixel no., Theta                 //
//             - the 3D-histogram   (sigma^2-sigmabar^2), pixel no., Theta  //
//             - the 2D-histogram   Theta    vs. Phi                        //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHSigmaTheta.h"

#include <TCanvas.h>

#include "MTime.h"
#include "MPointingPos.h"

#include "MBinning.h"
#include "MParList.h"

#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MBadPixelsCam.h"

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

#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"

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

ClassImp(MHSigmaTheta);

using namespace std;

// --------------------------------------------------------------------------
//
// Constructor. 
//
MHSigmaTheta::MHSigmaTheta(const char *name, const char *title)
{
    fName  = name  ? name  : "MHSigmaTheta";
    fTitle = title ? title : "2D histogram sigmabar vs. Theta";

    fThetaPhi = new TH2D;
    fThetaPhi->SetDirectory(NULL);
    fThetaPhi->UseCurrentStyle();
    fThetaPhi->SetName("2D-ThetaPhi");
    fThetaPhi->SetTitle("2D: \\Theta vs. \\Phi");
    fThetaPhi->SetXTitle("\\phi [\\circ]");
    fThetaPhi->SetYTitle("\\Theta [\\circ]");
    fThetaPhi->SetTitleOffset(1.2,"Y");

    fSigmaTheta = new TH2D;
    fSigmaTheta->SetDirectory(NULL);
    fSigmaTheta->UseCurrentStyle();
    fSigmaTheta->SetName("2D-ThetaSigmabar(Inner)");
    fSigmaTheta->SetTitle("2D: \\bar{\\sigma}, \\Theta");
    fSigmaTheta->SetXTitle("\\Theta [\\circ]");
    fSigmaTheta->SetYTitle("Sigmabar(Inner)");
    fSigmaTheta->SetTitleOffset(1.2,"Y");

    fSigmaThetaOuter = new TH2D;
    fSigmaThetaOuter->SetDirectory(NULL);
    fSigmaThetaOuter->UseCurrentStyle();
    fSigmaThetaOuter->SetName("2D-ThetaSigmabar(Outer)");
    fSigmaThetaOuter->SetTitle("2D: \\bar{\\sigma}, \\Theta");
    fSigmaThetaOuter->SetXTitle("\\Theta [\\circ]");
    fSigmaThetaOuter->SetYTitle("Sigmabar(Outer)");
    fSigmaThetaOuter->SetTitleOffset(1.2,"Y");

    fSigmaPixTheta = new TH3D;
    fSigmaPixTheta->SetDirectory(NULL);
    fSigmaPixTheta->UseCurrentStyle();
    fSigmaPixTheta->SetName("3D-ThetaPixSigma");
    fSigmaPixTheta->SetTitle("3D: \\Theta, Pixel Id, \\sigma");
    fSigmaPixTheta->SetXTitle("\\Theta [\\circ]");
    fSigmaPixTheta->SetYTitle("Pixel Id");
    fSigmaPixTheta->SetZTitle("Sigma");

    fDiffPixTheta = new TH3D;
    fDiffPixTheta->SetDirectory(NULL);
    fDiffPixTheta->UseCurrentStyle();
    fDiffPixTheta->SetName("3D-ThetaPixDiff");
    //fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}");
    fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, (Sigma2-Sigmabar2)/Area");
    fDiffPixTheta->SetXTitle("\\Theta [\\circ]");
    fDiffPixTheta->SetYTitle("Pixel Id");
    fDiffPixTheta->SetZTitle("(Sigma2 - Sigmabar2)/Area");

    // Define default binning
    fBinsPhi = new MBinning;
    fBinsPhi->SetEdges( 20, 0.0, 360.0);

    Double_t fThetaLo =  0.0;
    Double_t fThetaHi = 90.0;
    fBinsTheta = new MBinning;
    fBinsTheta->SetEdgesCos( 10,  fThetaLo, fThetaHi);   // theta
    //fBinsTheta->SetEdges( 1, fThetaLo, fThetaHi);   // theta

    fBinsSigma = new MBinning;
    fBinsSigma->SetEdges(   100,     0.0, 120.0);   // sigma 

    fBinsSigmabarIn  = new MBinning;
    fBinsSigmabarOut = new MBinning;;
    fBinsSigmabarIn->SetEdges( 100,     0.0,  25.0);   // sigmabar (inner)
    fBinsSigmabarOut->SetEdges(100,     0.0,  60.0);   // sigmabar (outer)

    fBinsPix = new MBinning;
    fBinsPix->SetEdges(578, -0.5, 577.5);

    fBinsDiff = new MBinning;
    fBinsDiff->SetEdges(   100,  -500.0, 1500.0);   // (sigma2-sigmabar2)/area

    SetBinning(fThetaPhi,        fBinsPhi,   fBinsTheta);
    SetBinning(fSigmaTheta,      fBinsTheta, fBinsSigmabarIn);
    SetBinning(fSigmaThetaOuter, fBinsTheta, fBinsSigmabarOut);
    SetBinning(fSigmaPixTheta,   fBinsTheta, fBinsPix, fBinsSigma);
    SetBinning(fDiffPixTheta,    fBinsTheta, fBinsPix, fBinsDiff);

    //--------------------------------------------
    fNamePedPhotCam = "MPedPhotCamFromData";
}


// --------------------------------------------------------------------------
//
// Destructor. 
//
MHSigmaTheta::~MHSigmaTheta()
{
  delete fThetaPhi;
  delete fSigmaTheta;
  delete fSigmaThetaOuter;
  delete fSigmaPixTheta;
  delete fDiffPixTheta;

  delete fBinsPhi;
  delete fBinsTheta;
  delete fBinsSigma;
  delete fBinsSigmabarIn;
  delete fBinsSigmabarOut;
  delete fBinsPix;
  delete fBinsDiff;
}

// --------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histogram
//
Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
{
    fCam = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

    fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
    if (!fPointPos)
    {
        *fLog << err << "MPointingPos not found... aborting." << endl;
        return kFALSE;
    }

    fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
    if (!fPed)
    {
        *fLog << err << AddSerialNumber(fNamePedPhotCam) 
              << "[MPedPhotCam] not found... aborting." << endl;
        return kFALSE;
    }
    fPed->InitSize(fCam->GetNumPixels()); 

    
    fBad = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
    if (!fBad)
    {  
       *fLog << err << "MBadPixelsCam not found... continue. " << endl; 
    }


    fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
    if (!fEvt)
    {
        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }


    //---------------------------------------------------------------
    *fLog << inf << "Name of MPedPhotCam container : " 
          << fNamePedPhotCam << endl;


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

    // Get Phi Binning
    MBinning* binsphi  = (MBinning*)plist->FindObject("BinningPhi", "MBinning");
    if (!binsphi)
    {
        *fLog << warn << "Object 'BinningPhi' [MBinning] not found... use default binning." << endl;
        binsphi = fBinsPhi;
    }

    // Get Theta Binning
    MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
    if (!binstheta)
    {
        *fLog << warn << "Object 'BinningTheta' [MBinning] not found... use default binning." << endl;
        binstheta = fBinsTheta;
    }

    // Get Sigma binning
    MBinning* binssig = (MBinning*)plist->FindObject("BinningSigma", "MBinning");
    if (!binssig)
    {
        *fLog << warn << "Object 'BinningSigma' [MBinning] not found... use default binning." << endl;
        binssig = fBinsSigma;
    }

    // Get SigmabarIn binning
    MBinning* binssigmain = (MBinning*)plist->FindObject("BinningSigmabarIn", "MBinning");
    if (!binssigmain)
    {
        *fLog << warn << "Object 'BinningSigmabarIn' [MBinning] not found... use default binning." << endl;
        binssigmain = fBinsSigmabarIn;
    }

    // Get SigmabarOut binning
    MBinning* binssigmaout = (MBinning*)plist->FindObject("BinningSigmabarOut", "MBinning");
    if (!binssigmaout)
    {
        *fLog << warn << "Object 'BinningSigmabarOut' [MBinning] not found... use default binning." << endl;
        binssigmaout = fBinsSigmabarOut;
    }

    // Get binning for (sigma^2-sigmabar^2)
    MBinning* binsdiff  = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
    if (!binsdiff)
    {
        *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... use default binning." << endl;
        binsdiff = fBinsDiff;
    }


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

    // Get binning for pixel number
    const UInt_t npix1 = fPed->GetSize()+1;
    //*fLog << "MHSigmaTheata::SetupFill(); npix1 = " << npix1 << endl;
    MBinning binspix("BinningPixel");
    binspix.SetEdges(npix1, -0.5, npix1-0.5);

    // Set binnings in histograms
    SetBinning(fThetaPhi,        binsphi,   binstheta);
    SetBinning(fSigmaTheta,      binstheta, binssigmain);
    SetBinning(fSigmaThetaOuter, binstheta, binssigmaout);
    SetBinning(fSigmaPixTheta,   binstheta, &binspix, binssig);
    SetBinning(fDiffPixTheta,    binstheta, &binspix, binsdiff);

    *fLog << "MHSigmaTheta::SetupFill(); binnings were set" << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Fill the histograms
//
//  ignore pixels if they are unused or blind
//
Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
{
    Double_t theta = fPointPos->GetZd();
    Double_t phi   = fPointPos->GetAz();

    fPed->ReCalc(*fCam, fBad);

    Double_t mysig      = (fPed->GetArea(0)).GetRms();
    Double_t mysigouter = (fPed->GetArea(1)).GetRms();

    //*fLog << "theta, mysig, mysigouter = " << theta << ",  " << mysig 
    //      << ",  " << mysigouter << endl;

    fThetaPhi->Fill(phi, theta);
    fSigmaTheta->Fill(theta, mysig);
    fSigmaThetaOuter->Fill(theta, mysigouter);

    MCerPhotPix *cerpix = 0;
    TIter Next(*fEvt);
    while ( (cerpix=(MCerPhotPix*)Next()) )
    {
        const Int_t id = cerpix->GetPixId();

        if (!cerpix->IsPixelUsed())
	{
          //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
          //      << id << endl;
          continue;
	}

        const MPedPhotPix &pix = (*fPed)[id];

        // ratio is the area of pixel 0 
        //          divided by the area of the current pixel
        const Double_t ratio = fCam->GetPixRatio(id);
        const Double_t sigma = pix.GetRms();

        fSigmaPixTheta->Fill(theta, (Double_t)id, sigma);

	Double_t diff;
        const Byte_t aidx = (*fCam)[id].GetAidx();
        if (aidx == 0)
	{
          // inner pixel
          diff = (sigma*sigma - mysig*mysig) * ratio;
	}
        else
	{
          // outer pixel
          diff = (sigma*sigma - mysigouter*mysigouter) * ratio;
	}

        fDiffPixTheta->Fill(theta, (Double_t)id, diff);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Update the projections and (if possible) set log scales before painting
//
void MHSigmaTheta::Paint(Option_t *opt)
{
    TVirtualPad *padsave = gPad;

    TH1D* h;

    padsave->cd(1);
    if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
    {
        ProjectionX(*h, *fSigmaTheta);
        if (h->GetEntries()!=0)
            gPad->SetLogy();
    }

    padsave->cd(4);
    if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
        ProjectionY(*h, *fSigmaTheta);

    gPad = padsave;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHSigmaTheta::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);
    AppendPad("");

    pad->Divide(3, 3);

    // draw the 2D histogram Sigmabar versus Theta
    TH1 *h;

    pad->cd(1);
    gPad->SetBorderMode(0);
    h = fSigmaTheta->ProjectionX("ProjX-Theta", -1, 9999, "E");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    h->SetTitle("Distribution of \\Theta");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("No.of events");
    h->SetTitleOffset(1.2,"Y");
    h->Draw(opt);
    h->SetBit(kCanDelete);

    pad->cd(2);
    gPad->SetBorderMode(0);
    h = fDiffPixTheta->Project3D("zx");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
    h->SetTitle("(Sigma2-Sigmabar2)/Area  vs. \\Theta (all pixels)");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
    h->SetTitleOffset(1.2,"Y");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(3);
    gPad->SetBorderMode(0);
    h = fSigmaPixTheta->Project3D("zx");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("Sigma");
    h->SetTitleOffset(1.2,"Y");
    h->Draw("box");
    h->SetBit(kCanDelete);

    //pad->cd(7);
    //gPad->SetBorderMode(0);
    //h = fSigmaTheta->ProjectionY("ProjY-sigma", -1, 9999, "E");
    //h->SetDirectory(NULL);
    //h->UseCurrentStyle();
    //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
    //h->SetXTitle("\\bar{\\sigma}_{ped}");
    //h->SetYTitle("No.of events");
    //h->SetTitleOffset(1.2,"Y");
    //h->Draw(opt);
    //h->SetBit(kCanDelete);

    pad->cd(5);
    gPad->SetBorderMode(0);
    h = fDiffPixTheta->Project3D("zy");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all  \\Theta)");
    h->SetTitle("(Sigma2-Sigmabar2)/Area  vs. pixel Id (all  \\Theta)");
    h->SetXTitle("Pixel Id");
    h->SetYTitle("(Sigma2 - sigmabar2) / Area");
    h->SetTitleOffset(1.2,"Y");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(6);
    gPad->SetBorderMode(0);
    h = fSigmaPixTheta->Project3D("zy");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    h->SetTitle("\\sigma_{ped} vs. pixel Id (all  \\Theta)");
    h->SetXTitle("Pixel Id");
    h->SetYTitle("Sigma");
    h->SetTitleOffset(1.2,"Y");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(4);
    gPad->SetBorderMode(0);
    fSigmaTheta->Draw();

    pad->cd(7);
    gPad->SetBorderMode(0);
    fSigmaThetaOuter->Draw();

    pad->cd(8);
    gPad->SetBorderMode(0);
    fThetaPhi->Draw();

    pad->cd(9);
    gPad->SetBorderMode(0);
    h = fThetaPhi->ProjectionX("ProjX-Phi", -1, 9999, "E");
    h->SetDirectory(NULL);
    h->UseCurrentStyle();
    h->SetTitle("Distribution of \\Phi");
    h->SetXTitle("\\phi [\\circ]");
    h->SetYTitle("No. of events");
    h->SetTitleOffset(1.2,"Y");
    h->Draw(opt);
    h->SetBit(kCanDelete);
}