/* ======================================================================== *\
!
! *
! * 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  07/2004 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHTelAxisFromStars
//
// This class contains histograms for MTelAxisFromStars
//
/////////////////////////////////////////////////////////////////////////////
#include "MHTelAxisFromStars.h"

#include <math.h>

#include <TH1.h>
#include <TH2.h>
#include <TPad.h>
#include <TStyle.h>
#include <TCanvas.h>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MBinning.h"

#include "MStarCam.h"
#include "MStarPos.h"
#include "MSkyCamTrans.h"
#include "MSrcPosCam.h"


ClassImp(MHTelAxisFromStars);

using namespace std;

// --------------------------------------------------------------------------
//
// Setup the histograms 
//
MHTelAxisFromStars::MHTelAxisFromStars(const char *name, const char *title)
    : fMm2Deg(1), fUseMmScale(kFALSE)
{
    fName  = name  ? name  : "MHTelAxisFromStars";
    fTitle = title ? title : "Plots for MTelAxisFromStars";

    fNStars     = new TH1D("NStars", "No. of stars", 10, -0.5, 9.5);
    fNStars->SetDirectory(NULL);
    fNStars->SetXTitle("Numder of stars");    
    fNStars->SetYTitle("Counts");

    fNdoF     = new TH1D("NdoF", "No. of degrees of freedom", 20, -0.5, 19.5);
    fNdoF->SetDirectory(NULL);
    fNdoF->SetXTitle("Numder of degrees of freedom");    
    fNdoF->SetYTitle("Counts");

    fLog10Chi2 = new TH1D("log10Chi2", "log10Chi2", 60, -10, 5);
    fLog10Chi2->SetDirectory(NULL);
    fLog10Chi2->SetXTitle("log10(Chi2)");
    fLog10Chi2->SetYTitle("Counts");

    fChi2Prob = new TH1D("Chi2-Prob", "Chi2 probability",     40,  0.0,  1.0);
    fChi2Prob->SetDirectory(NULL);
    fChi2Prob->SetXTitle("Chi2 probability");
    fChi2Prob->SetYTitle("Counts");


    fNumIter  = new TH1D("NumIter", "Number of iterations",   50, -0.5, 49.5);
    fNumIter->SetDirectory(NULL);
    fNumIter->SetXTitle("Number of iterations");
    fNumIter->SetYTitle("Counts");


    fLambda   = new TH1D("Lambda", "Scale factor lambda",     80, 0.90, 1.10);
    fLambda->SetDirectory(NULL);
    fLambda->SetXTitle("Scale factor lambda");
    fLambda->SetYTitle("Counts");


    fAlfa     = new TH1D("Alfa", "Rotation angle alfa",      100, -2.5,  2.5);
    fAlfa->SetDirectory(NULL);
    fAlfa->SetXTitle("Rotation angle alfa [\\circ]");
    fAlfa->SetYTitle("Counts");


    fShift = new TH2D("Shift", "Sky-Cam transformnation : (x,y) shift", 
                      72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
    fShift->SetDirectory(NULL);
    fShift->SetZTitle("Counts");
    fShift->SetXTitle("x-shift [\\circ]");
    fShift->SetYTitle("y-shift [\\circ]");


    fEstPos1 = new TH2D("EstPos1", "Estimated position1", 
                      72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
    fEstPos1->SetDirectory(NULL);
    fEstPos1->SetZTitle("Counts");
    fEstPos1->SetXTitle("Estimated position1    x [\\circ]");
    fEstPos1->SetYTitle("Estimated position1    y [\\circ]");

    fEstPos2 = new TH2D("EstPos2", "Estimated position2", 
                      72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
    fEstPos2->SetDirectory(NULL);
    fEstPos2->SetZTitle("Counts");
    fEstPos2->SetXTitle("Estimated position2    x [\\circ]");
    fEstPos2->SetYTitle("Estimated position2    y [\\circ]");

    fEstPos3 = new TH2D("EstPos3", "Estimated position3", 
                      72, -534.0707, 534.0707, 72, -534.0707, 534.0707);
    fEstPos3->SetDirectory(NULL);
    fEstPos3->SetZTitle("Counts");
    fEstPos3->SetXTitle("Estimated position3    x [\\circ]");
    fEstPos3->SetYTitle("Estimated position3    y [\\circ]");


    // default input type : results from Gauss fit
    fInputType = 1;
}

// --------------------------------------------------------------------------
//
// Set the type of the input
//
//     type = 0        calculated star positions (by averaging)
//     type = 1        fitted star positions (by Gauss fit)
//
void MHTelAxisFromStars::SetInputType(Int_t type)
{
  *fLog << all << "MHTelAxisFromStars::SetInputType; input type is set equal to : " 
        << type ;

  if (type == 0)
    *fLog << " (calculated star positions)" << endl;
  else
    *fLog << " (fitted star positions)" << endl;


  fInputType = type;
}


// --------------------------------------------------------------------------
//
// Delete the histograms
//
MHTelAxisFromStars::~MHTelAxisFromStars()
{
  //*fLog << "MHTelAxisFromStars::Destructor" << endl;

    delete fNStars;
    delete fNdoF;
    delete fLog10Chi2;
    delete fChi2Prob;
    delete fNumIter;
    delete fLambda;
    delete fAlfa;

    delete fShift;
    delete fEstPos1;
    delete fEstPos2;
    delete fEstPos3;
}

// --------------------------------------------------------------------------
//
// Setup the Binning for the histograms. 
//
// Find the pointers to the containers.
//
// Use this function if you want to set the conversion factor which
// is used to convert the mm-scale in the camera plane into the deg-scale.
// The conversion factor is part of the camera geometry. Please create a 
// corresponding MGeomCam container.
//
Bool_t MHTelAxisFromStars::SetupFill(const MParList *plist)
{

   fStarCam = (MStarCam*)plist->FindObject("MStarCam", "MStarCam");
   if (!fStarCam)
   {
       *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MStarCam' not found... aborting." << endl;
       return kFALSE;
   }

   
   fSourceCam = (MStarCam*)plist->FindObject("MSourceCam", "MStarCam");
   if (!fSourceCam)
   {
       *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MSourceCam' not found... " << endl;
   }
   


    fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
    if (!fSrcPos)
    {
        *fLog << err << "MHTelAxisFromStars::SetupFill;  MSrcPosCam not found...  aborting" << endl;
        return kFALSE;
    }

    fSkyCamTrans = (MSkyCamTrans*)plist->FindObject(AddSerialNumber("MSkyCamTrans"));
    if (!fSkyCamTrans)
    {
        *fLog << err << "MHTelAxisFromStars::SetupFill;  MSkyCamTrans not found...  aborting" << endl;
        return kFALSE;
    }


   //------------------------------------------------------------------
    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!geom)
    {
        *fLog << warn << GetDescriptor() 
          << ": No Camera Geometry available. Using mm-scale for histograms." 
          << endl;
        SetMmScale(kTRUE);
    }
    else
    {
        fMm2Deg = geom->GetConvMm2Deg();
        SetMmScale(kFALSE);
    }

    ApplyBinning(*plist, "NStars",   fNStars);
    ApplyBinning(*plist, "NdoF",     fNdoF);
    ApplyBinning(*plist, "Log10Chi2",fLog10Chi2);
    ApplyBinning(*plist, "Chi2Prob", fChi2Prob);
    ApplyBinning(*plist, "NumIter",  fNumIter);
    ApplyBinning(*plist, "Lambda",   fLambda);
    ApplyBinning(*plist, "Alfa",     fAlfa);

    const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
    if (!bins)
    {
        float r = geom ? geom->GetMaxRadius() : 600;
        r *= 0.9;
        if (!fUseMmScale)
            r *= fMm2Deg;

        MBinning b;
        b.SetEdges(61, -r, r);
        SetBinning(fShift,   &b, &b);
        SetBinning(fEstPos1, &b, &b);
        SetBinning(fEstPos2, &b, &b);
        SetBinning(fEstPos3, &b, &b);
    }
    else
    {
        SetBinning(fShift,   bins, bins);
        SetBinning(fEstPos1, bins, bins);
        SetBinning(fEstPos2, bins, bins);
        SetBinning(fEstPos3, bins, bins);
    }

    //-------------------------------------------    
    // copy names from MStarCam to the histograms
    MStarPos* starpos;
    Int_t istar = 0;
    TIter Next(fSourceCam->GetList());

    while ((starpos=(MStarPos*)Next())) {
      fStarnames[istar] =  starpos->GetName();
      //*fLog << "istar, star name = " << istar << ",  " 
      //      << fStarnames[istar] << endl;
      istar++;
      if (istar >= fNstarnames) break;
    }

    if (fSourceCam)
    {
      MStarPos *starSource = 0;
      TIter nextSource(fSourceCam->GetList());
      while ( (starSource = (MStarPos*)nextSource()) )
      {
         if     ( fNstarnames > 0  &&  starSource->GetName() == fStarnames[0] )
           fEstPos1->SetName(starSource->GetName());

         else if( fNstarnames > 1  &&  starSource->GetName() == fStarnames[1] )
           fEstPos2->SetName(starSource->GetName());

         else if( fNstarnames > 2  &&  starSource->GetName() == fStarnames[2] )
           fEstPos3->SetName(starSource->GetName());
      }
    }
    

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Use this function to setup your own conversion factor between degrees
// and millimeters. The conversion factor should be the one calculated in
// MGeomCam. Use this function with Caution: You could create wrong values
// by setting up your own scale factor.
//
void MHTelAxisFromStars::SetMm2Deg(Float_t mmdeg)
{
    if (mmdeg<0)
    {
        *fLog << warn << dbginf 
              << "Warning - Conversion factor < 0 - nonsense. Ignored." 
              << endl;
        return;
    }

    if (fMm2Deg>=0)
        *fLog << warn << dbginf 
              << "Warning - Conversion factor already set. Overwriting" 
              << endl;

    fMm2Deg = mmdeg;
}

// --------------------------------------------------------------------------
//
// With this function you can convert the histogram ('on the fly') between
// degrees and millimeters.
//
void MHTelAxisFromStars::SetMmScale(Bool_t mmscale)
{
    if (fUseMmScale == mmscale)
        return;

    if (fMm2Deg<0)
    {
        *fLog << warn << dbginf 
         << "Warning - Sorry, no conversion factor for conversion available." 
         << endl;
        return;
    }

    const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg;
    MH::ScaleAxis(fShift,   scale, scale);
    MH::ScaleAxis(fEstPos1, scale, scale);
    MH::ScaleAxis(fEstPos2, scale, scale);
    MH::ScaleAxis(fEstPos3, scale, scale);

    if (mmscale)
    {
        fShift->SetXTitle("x-shift [mm]");
        fShift->SetYTitle("y-shift [mm]");

        fEstPos1->SetXTitle("Estimated position1 x [mm]");
        fEstPos1->SetYTitle("Estimated position1 y [mm]");

        fEstPos2->SetXTitle("Estimated position2 x [mm]");
        fEstPos2->SetYTitle("Estimated position2 y [mm]");

        fEstPos3->SetXTitle("Estimated position3 x [mm]");
        fEstPos3->SetYTitle("Estimated position3 y [mm]");
    }
    else
    {
        fShift->SetXTitle("x-shift [\\circ]");
        fShift->SetYTitle("y-shift [\\circ]");

        fEstPos1->SetXTitle("Estimated position1 x [\\circ]");
        fEstPos1->SetYTitle("Estimated position1 y [\\circ]");

        fEstPos2->SetXTitle("Estimated position2 x [\\circ]");
        fEstPos2->SetYTitle("Estimated position2 y [\\circ]");

        fEstPos3->SetXTitle("Estimated position3 x [\\circ]");
        fEstPos3->SetYTitle("Estimated position3 y [\\circ]");
    }

    fUseMmScale = mmscale;
}

// --------------------------------------------------------------------------
//
// Fill the histograms 
// 
//
Bool_t MHTelAxisFromStars::Fill(const MParContainer *par, const Stat_t w)
{
    Int_t Ndof = fSkyCamTrans->GetNdof();
    if (Ndof < 0)
      return kTRUE;

    const Double_t scale = fUseMmScale ? 1 : fMm2Deg;


    fNStars   ->Fill(fSkyCamTrans->GetNStars(),               w);
    fNdoF     ->Fill(fSkyCamTrans->GetNdof(),                 w);
    if (fSkyCamTrans->GetChiSquare() > 0.0)
      fLog10Chi2->Fill( log10(fSkyCamTrans->GetChiSquare() ), w);
    fChi2Prob ->Fill(fSkyCamTrans->GetChiSquareProb(),        w);
    fNumIter  ->Fill(fSkyCamTrans->GetNumIter(),              w);
    fLambda   ->Fill(fSkyCamTrans->GetLambda(),               w);
    fAlfa     ->Fill(fSkyCamTrans->GetAlfa(),                 w);


    Double_t lowx;
    Double_t higx;
    Double_t lowy;
    Double_t higy;
    Double_t x;
    Double_t y;
    Int_t    nbix;
    Int_t    nbiy;

    nbix = fShift->GetXaxis()->GetNbins();
    lowx = fShift->GetXaxis()->GetBinLowEdge(1);
    higx = fShift->GetXaxis()->GetBinLowEdge(nbix+1);

    nbiy = fShift->GetYaxis()->GetNbins();
    lowy = fShift->GetYaxis()->GetBinLowEdge(1);
    higy = fShift->GetYaxis()->GetBinLowEdge(nbiy+1);

    x = scale * (fSkyCamTrans->GetShiftD())[0];
    y = scale * (fSkyCamTrans->GetShiftD())[1];
    if (x>lowx  && x<higx  && y>lowy  && y<higy)
    {
      fShift    ->Fill(x, y, w);
    }

    if (fSourceCam)
    {
      MStarPos *starSource = 0;
      TIter nextSource(fSourceCam->GetList());

      while( (starSource = (MStarPos*)nextSource()) )
      {
        if (fInputType == 1)
        {
          x = scale * starSource->GetMeanXFit(); 
          y = scale * starSource->GetMeanYFit(); 
        }
        else 
        {
          x = scale * starSource->GetMeanXCalc(); 
          y = scale * starSource->GetMeanYCalc(); 
        }

	if (x>lowx  && x<higx  && y>lowy  && y<higy)
        {
         if     ( fNstarnames > 0  &&  starSource->GetName() == fStarnames[0] )
           fEstPos1->Fill(x, y, w);

         else if( fNstarnames > 1  &&  starSource->GetName() == fStarnames[1] )
           fEstPos2->Fill(x, y, w);

         else if( fNstarnames > 2  &&  starSource->GetName() == fStarnames[2] )
           fEstPos3->Fill(x, y, w);
	}
      }
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Finalize :
//
//    it is called in the postprocessing
//    reset pointers in order to allow cloning of the object
//
Bool_t MHTelAxisFromStars::Finalize()
{
  //*fLog << "MHTelAxisFromStars::Finalize;  fSourceCam = " 
  //      << fSourceCam  << endl; 

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the histograms.
// Be careful: The histograms belongs to this object and won't get deleted
// together with the canvas.
//

void MHTelAxisFromStars::Draw(Option_t *opt)
{
  TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
  pad->SetBorderMode(0);

  AppendPad("");
  
  //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 1200, 900);
  //gROOT->SetSelectedPad(NULL);

    pad->Divide(4,3);

    pad->cd(1);
    gPad->SetBorderMode(0);
    fNStars->Draw(opt); 

    pad->cd(2);
    gPad->SetBorderMode(0);
    fNdoF->Draw(opt); 

    pad->cd(3);
    gPad->SetBorderMode(0);
    fLog10Chi2->Draw(opt); 

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

    pad->cd(5);
    gPad->SetBorderMode(0);
    fNumIter->Draw(opt);

    pad->cd(6);
    gPad->SetBorderMode(0);
    fLambda->Draw(opt);

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

    pad->cd(8);
    gPad->SetBorderMode(0);
    //SetColors();
    //fShift->Draw("colz");
    fShift->Draw("");


    //-----------------------------------------------
    // plot the expected positions of some sources
    //*fLog << "fSourcsCam = " << fSourceCam << endl;

    if (fSourceCam)
    {
      *fLog << "MHTelAxisFromSrars::Draw; plotting" << endl;

      pad->cd(9);
      gPad->SetBorderMode(0);
      //SetColors();
      fEstPos1->Draw("");

      pad->cd(10);
      gPad->SetBorderMode(0);
      //SetColors();
      fEstPos2->Draw("");

      pad->cd(11);
      gPad->SetBorderMode(0);
      //SetColors();
      fEstPos3->Draw("");
    }

    pad->Modified();
    pad->Update();
}

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

TH1 *MHTelAxisFromStars::GetHistByName(const TString name)
{
    if (name.Contains("NStars", TString::kIgnoreCase))
        return fNStars;
    if (name.Contains("NdoF", TString::kIgnoreCase))
        return fNdoF;
    if (name.Contains("Log10Chi2", TString::kIgnoreCase))
        return fLog10Chi2;
    if (name.Contains("Chi2Prob", TString::kIgnoreCase))
        return fChi2Prob;
    if (name.Contains("NumIter", TString::kIgnoreCase))
        return fNumIter;
    if (name.Contains("Lambda", TString::kIgnoreCase))
        return fLambda;
    if (name.Contains("Alfa", TString::kIgnoreCase))
        return fAlfa;

    if (name.Contains("Shift", TString::kIgnoreCase))
        return fShift;
    if (name.Contains("EstPos1", TString::kIgnoreCase))
        return fEstPos1;
    if (name.Contains("EstPos2", TString::kIgnoreCase))
        return fEstPos2;
    if (name.Contains("EstPos3", TString::kIgnoreCase))
        return fEstPos3;

    return NULL;
}

// --------------------------------------------------------------------------
//
// Setup an inversed deep blue sea palette.
//

void MHTelAxisFromStars::SetColors() const
{
    gStyle->SetPalette(51, NULL);
    Int_t c[50];
    for (int i=0; i<50; i++)
        c[49-i] = gStyle->GetColorPalette(i);
    gStyle->SetPalette(50, c);
}


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

void MHTelAxisFromStars::Paint(Option_t *opt)
{
    SetColors();
    MH::Paint();
}

//==========================================================================