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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHFlux                                                                 //
//                                                                          //
//  calculates absolute photon fluxes                                       //
//             from the distributions of the estimated energy               //
//                      for the different bins in some variable 'Var'       //
//                      (Var = Theta or time)                               //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHFlux.h"

#include <TStyle.h>

#include <TF1.h>
#include <TH2.h>
#include <TProfile.h>


#include <TCanvas.h>

#include "MTime.h"

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

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

ClassImp(MHFlux);


// --------------------------------------------------------------------------
//
// Default Constructor. It sets the variable name (Theta or time)
//                      and the units for the variable
// 
MHFlux::MHFlux(const TH2D &h2d,  const Bool_t Draw,
	       const TString varname, const TString unit)
    : fHOrig(), fHUnfold(), fHFlux()
{
    if (varname.IsNull() || unit.IsNull())
        *fLog << warn << dbginf << "varname or unit not defined" << endl;

    fVarname = varname;
    fUnit    = unit;

    TString strg(varname);
    strg += unit;

    // char txt[100];

    // original distribution of E-est for different bins 
    //                       of the variable (Theta or time)
    // sprintf(txt, "gammas vs. E-est and %s",varname);

    TString strg1 = "no.of gammas vs. E-est and ";
    strg1 += varname;
    
    ((TH2D&)h2d).Copy(fHOrig);

    fHOrig.SetName("E-est");
    fHOrig.SetTitle(strg1);

    fHOrig.SetDirectory(NULL);
    fHOrig.SetXTitle("E-est [GeV]");
    fHOrig.SetYTitle(strg);
    fHOrig.Sumw2();

    SetBinning((TH2*)&fHOrig, (TH2*)&h2d);


    // unfolded distribution of E-unfold for different bins 
    //                       of the variable (Theta or time)
    // sprintf(txt, "gammas vs. E-unfold and %s",varname);
    TString strg2 = "no.of gammas vs. E-unfold and ";
    strg2 += varname;

    fHUnfold.SetName("E-unfolded");
    fHUnfold.SetTitle(strg2);

    fHUnfold.SetDirectory(NULL);
    fHUnfold.SetXTitle("E-unfold [GeV]");
    fHUnfold.SetYTitle(strg);
    fHUnfold.Sumw2();
    
    SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);


    // absolute photon flux vs. E-unfold
    //          for different bins of the variable (Theta or time)
    //
    // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
    TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and ";
    strg3 += varname;

    fHFlux.SetName("photon flux");
    fHFlux.SetTitle(strg3);

    fHFlux.SetDirectory(NULL);
    fHFlux.SetXTitle("E-unfold [GeV]");
    fHFlux.SetYTitle(strg);
    fHFlux.Sumw2();

    SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);


    // copy fHOrig into fHUnfold in case no unfolding is done
    const Int_t nEunf   = fHUnfold.GetNbinsX();
    const Int_t nVar    = fHUnfold.GetNbinsY();
    for (int m=1; m<=nEunf; m++)
    {
      for (int n=1; n<=nVar; n++)
      {
        Double_t cont  = fHOrig.GetBinContent(m,n);
        Double_t dcont = fHOrig.GetBinError(m,n);
        fHUnfold.SetBinContent(m,n,cont);
        fHUnfold.SetBinError(m,n,dcont);
      }
    }
  //..............................................
  // draw the No.of photons vs. E-est 
  // for the individual bins of the variable Var

  if (Draw == kTRUE)
  {
    const Int_t nVar = fHOrig.GetNbinsY();

    for (int n=1; n<=nVar; n++)
    {
      TString strg0("Orig-");
      strg0 += fVarname;
      TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");

      strg0  = fVarname;
      strg0 += "-bin ";
      strg0 += n;

      TString strg1("No.of photons vs. E-est for ");
      strg1 += strg0;

      new TCanvas(strg0, strg1);
      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
      // gROOT->SetSelectedPad(NULL);

      gPad->SetLogx();

      h.SetName(strg0);
      h.SetTitle(strg1);
      h.SetXTitle("E-est [GeV]");
      h.SetYTitle("No.of photons");
      h.DrawCopy();

      // c.Modified();
      // c.Update();
    }
  }
  //........................
}

// -------------------------------------------------------------------------
//
// Dummy Fill (has to be included because in base class MH Fill is set to 0 
//             (abstract member function));
// without the dummy Fill one gets the error message :
//
// Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope 
// FILE:macros/flux.C LINE:465
// Possible candidates are...
// filename       line:size busy function type and name  (in MHFlux)
// filename       line:size busy function type and name  (in MH)
// filename       line:size busy function type and name  (in MParContainer)
// filename       line:size busy function type and name  (in TObject)
//
Bool_t MHFlux::Fill(const MParContainer *par)
{
   return kTRUE;
}


// -------------------------------------------------------------------------
//
// Unfold the distribution in E-est
//
void MHFlux::Unfold(const Bool_t Draw)
{
  //..............................................
  // draw the No.of photons vs. E-unfold 
  // for the individual bins of the variable Var

  if (Draw == kTRUE)
  {
    const Int_t nVar = fHUnfold.GetNbinsY();

    for (int n=1; n<=nVar; n++)
    {
      TString strg0("Unfold-");
      strg0 += fVarname;
      TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
      strg0  = fVarname;
      strg0 += "-bin ";
      strg0 += n;

      TString strg1("No.of photons vs. E-unfold for ");
      strg1 += strg0;

      new TCanvas(strg0, strg1);

      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
      // gROOT->SetSelectedPad(NULL);

      gPad->SetLogx();

      h.SetName(strg0);
      h.SetTitle(strg1);
      h.SetXTitle("E-unfold [GeV]");
      h.SetYTitle("No.of photons");
      h.DrawCopy();

      // c.Modified();
      // c.Update();
    }
  }
  //........................
}


// -------------------------------------------------------------------------
//
// Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by
//                       the width of the energy interval     (deltaE)
//                       the effective ontime                 (*teff)
//                       and the effective collection area    (*aeff)
//
void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
                      const TH2D *aeff, const Bool_t Draw)
{
  // Note that fHUnfold  has bins in Eunf and Var 
  //           *teff     has bins in Var  (the same bins in Var as fHUnfold)
  //           *thetabar has bins in Var  (the same bins in Var as fHUnfold)
  //           *aeff     has bins in Etru and Theta
  //                     (where in general the binning in Etru is different
  //                      from the binning in Eunf)
  // The variable Var may be 'time' or 'Theta'

  // Draw = kTRUE means the differential flux vs E-unf should be drawn
  //              for the individual bins of the variable Var

    const TAxis &axex = *((TH2*)aeff)->GetXaxis();
    const TAxis &axey = *((TH2*)aeff)->GetYaxis();

  //....................................
  // define dummy histogram *aeff
  ((TH1*)aeff)->Sumw2();
  MBinning binsetru("BinningEtru");
  binsetru.SetEdgesLog(10, 10, 1e3);

  MBinning binsthetatru("BinningThetatru");
  binsthetatru.SetEdges(7, -2.5, 32.5);
  //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
  SetBinning((TH2*)aeff, &binsetru, &binsthetatru);

  const Int_t netru  = aeff->GetNbinsX();
  const Int_t ntheta = aeff->GetNbinsY();

  for (int j=1; j<=netru; j++)
  {
      for (int k=1; k<=ntheta; k++)
      {
          Double_t cont = 10000.0;
          ((TH1*)aeff)->SetBinContent(j, k, cont);

          Double_t dcont = 100.0;
          ((TH1*)aeff)->SetBinError(j, k, dcont);
      }
  }
  // *fLog << "Dummy aeff : netru =" << netru << ",  ntheta = " << ntheta << endl;
  //....................................

  // number of Eunf and Var bins   (histograms : fHUnfold, fHFlux)
  const Int_t nEunf    = fHFlux.GetNbinsX();
  const Int_t nVar     = fHFlux.GetNbinsY();

  // number of Etru and Theta bins (histogram *aeff of collection area)

  const Int_t nEtru    = aeff->GetNbinsX();
  const Int_t nTheta   = aeff->GetNbinsY();

  //*fLog << "nEunf =" << nEunf << ",  nVar = "   << nVar   << endl;
  //*fLog << "nEtru =" << nEtru << ",  nTheta = " << nTheta << endl;

  //...................................................................
  // calculate effective collection area 
  //    for the Eunf and Var bins of the histogram fHUnfold
  //    from the histogram *aeff, which has bins in Etru and Theta
  // the result is the histogram fHAeff
  //

  TH2D fHAeff;
  fHAeff.Sumw2();
  SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);

  Double_t *aeffbar  = new Double_t[nEtru];
  Double_t *daeffbar = new Double_t[nEtru];

  Double_t aeffEunfVar;
  Double_t daeffEunfVar;

  //------   start n loop   ------
  for (int n=1; n<=nVar; n++)
  {
    Double_t Thetabar = thetabar->GetBinContent(n);
    Double_t cosThetabar = cos(Thetabar);

    // determine Theta bins (k1, k2, k3) for interpolation in Theta
    // k0 denotes the Theta bin from whicvh the error is copied
    Int_t k0=0;
    Int_t k1=0;
    Int_t k2=0;
    Int_t k3=0;

    for (int k=3; k<=nTheta; k++)
    {
        Double_t Thetalow = axey.GetBinLowEdge(k);
        if (Thetabar < Thetalow)
        {
            k1 = k-2;
            k2 = k-1;
            k3 = k;
            k0 = k2;
            break;
        }
    }

    if (k3 == 0)
    {
        k1 = nTheta-2;
        k2 = nTheta-1;
        k3 = nTheta;
        k0 = k2;
    }

    if (Thetabar <  axey.GetBinLowEdge(2))
      k0 = 1;
    else
        if (Thetabar >  axey.GetBinLowEdge(nTheta))
            k0 = nTheta;

    Double_t Thetamin = axey.GetBinLowEdge(1);
    Double_t Thetamax = axey.GetBinLowEdge(nTheta+1);
    if (Thetabar < Thetamin  ||  Thetabar > Thetamax)
    {
      *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
            << ",  Thetamin =" << Thetamin
            << ",  Thetamax =" << Thetamax << endl;
    } 

    //*fLog << "Var bin "   << n  << ":"  <<  endl;
    //*fLog << "Thetabar= " << Thetabar 
    //      << ",  k0= "    << k0
    //      << ",  k1= "    << k1
    //      << ",  k2= "    << k2
    //      << ",  k3= "    << k3         <<  endl;
 

    // calculate effective collection area at Theta = Thetabar
    // by quadratic interpolation in cos(Theta);
    // do this for each bin of Etru 
    //

    for (int j=1; j<=nEtru; j++)
    {
        double c0 = 0;
        double c1 = 0;
        double c2 = 0;

        const double t1 = cos( axey.GetBinCenter (k1) );
        const double t2 = cos( axey.GetBinCenter (k2) );
        const double t3 = cos( axey.GetBinCenter (k3) );

        const double a1 = aeff->GetBinContent(j, k1);
        const double a2 = aeff->GetBinContent(j, k2);
        const double a3 = aeff->GetBinContent(j, k3);

        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
        aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
        daeffbar[j] = aeff->GetBinError(j,k0);

        //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar
        //      << ",  abar= "  << aeffbar[j]
        //      << ",  dabar= " << daeffbar[j] << endl;
    }

    //---   start m loop ---
    // calculate effective collection area at (E = Ebar, Theta = Thetabar)
    // by quadratic interpolation in log10(Etru)
    // do this for each bin of Eunf
    //
    for (int m=1; m<=nEunf; m++)
    {
        Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )+
                                     log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
        Double_t Ebar = pow(10.0, log10Ebar);

        // determine Etru bins (j1, j2, j3) for interpolation in E
        // j0 denotes the Etru bin from which the error is copied
        Int_t j0=0;
        Int_t j1=0;
        Int_t j2=0;
        Int_t j3=0;

        for (int j=3; j<=nEtru; j++)
        {
            Double_t Elow = axex.GetBinLowEdge(j);
            if (Ebar < Elow)
            {
                j1 = j-2;
                j2 = j-1;
                j3 = j;
                j0 = j2;
                break;
            }
        }

        if (j3 == 0)
        {
            j1 = nEtru-2;
            j2 = nEtru-1;
            j3 = nEtru;
            j0 = j2;
        }

        if (Ebar <  axex.GetBinLowEdge(2))
            j0 = 1;
        else
            if (Ebar >  axex.GetBinLowEdge(nEtru))
                j0 = nEtru;

        Double_t Etrumin = axex.GetBinLowEdge(1);
        Double_t Etrumax = axex.GetBinLowEdge(nEtru+1);
        if (Ebar < Etrumin  ||  Ebar > Etrumax)
        {
            *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
                << ",  Etrumin =" << Etrumin
                << ",  Etrumax =" << Etrumax << endl;
        }

        //*fLog << "Var bin "   << n  << ":"  <<  endl;
        //*fLog << "Ebar= " << Ebar
        //      << ",  j1= "    << j1
        //      << ",  j2= "    << j2
        //      << ",  j3= "    << j3         <<  endl;


        double c0=0.0;
        double c1=0.0;
        double c2=0.0;

        const double t1 = 0.5 * ( log10( axex.GetBinLowEdge (j1)  )+
                                  log10( axex.GetBinLowEdge (j1+1)) );
        const double t2 = 0.5 * ( log10( axex.GetBinLowEdge (j2)  )+
                                  log10( axex.GetBinLowEdge (j2+1)) );
        const double t3 = 0.5 * ( log10( axex.GetBinLowEdge (j3)  )+
                                  log10( axex.GetBinLowEdge (j3+1)) );

        const double a1 = aeffbar[j1];
        const double a2 = aeffbar[j2];
        const double a3 = aeffbar[j3];

        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
        aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
        daeffEunfVar = daeffbar[j0];

        //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar
        //      << ",  aeffEunfVar = "  << aeffEunfVar
        //      << ",  daeffEunfVar = " << daeffEunfVar << endl;

        fHAeff.SetBinContent(m,n,aeffEunfVar);
        fHAeff.SetBinError(m,n,daeffEunfVar);
    }
    //---   end m loop ---
  }
  //------   end n loop   ------
  delete aeffbar;

  //...................................................................
  // now calculate the absolute gamma flux
  //
  for (int m=1; m<=nEunf; m++)
  {
      Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);

      for (int n=1; n<=nVar; n++)
      {
          Double_t Ngam   = fHUnfold.GetBinContent(m,n);
          Double_t dNgam  = fHUnfold.GetBinError(m,n);

          Double_t Aeff   = fHAeff.GetBinContent(m,n);
          Double_t dAeff  = fHAeff.GetBinError(m,n);

          Double_t Effon  = teff->GetBinContent(n);
          Double_t dEffon = teff->GetBinError(n);

          Double_t Cont, dCont;
          if (Ngam>0 && DeltaE>0 && Effon>0 && Aeff>0)
          {
              Cont  = Ngam / (DeltaE * Effon * Aeff);
              dCont = Cont * sqrt( dNgam *dNgam  / (Ngam*Ngam) +
                                   dEffon*dEffon / (Effon*Effon) +
                                   dAeff *dAeff  / (Aeff*Aeff)  );
          }
          else
          {
              Cont  = 1.e-20;
              dCont = 1.e-20;
          }

          fHFlux.SetBinContent(m,n,Cont);
          fHFlux.SetBinError(m,n,dCont);

          //*fLog << "Eunf bin "    << m      << ",  Var bin " << n
          //      << ":  Ngam = "   << Ngam   << ",  Flux = "
          //      << Cont  << ", dFlux = " << dCont << endl;
          //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
          //      << ",  Aeff = "   << Aeff   << endl;
      }
  }

  //..............................................
  // draw the differential photon flux vs. E-unf 
  // for the individual bins of the variable Var

  if (Draw == kTRUE)
  {
      for (int n=1; n<=nVar; n++)
      {
          TString strg0("Flux-");
          strg0 += fVarname;

          TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");

          TString strg1("Photon flux vs. E-unfold for ");
          TString strg2 = fVarname;

          strg2 += "-bin ";
          strg2 += n;

          TString strg3 = strg1 + strg2;
          new TCanvas(strg2, strg3);
          // TCanvas &c = *MakeDefCanvas(txt, txt);
          // gROOT->SetSelectedPad(NULL);

          gPad->SetLogx();

          h.SetName(strg2);
          h.SetTitle(strg3);
          h.SetXTitle("E-unfold [GeV]            ");
          h.SetYTitle("photons / (s m2 GeV)");
          h.DrawCopy();

          // c.Modified();
          // c.Update();
      }
  }
  //........................
}

// -------------------------------------------------------------------------
//
// Draw copies of the histograms
//
TObject *MHFlux::DrawClone(Option_t *opt) const
{
    TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots");
    c.Divide(2, 2);

    gROOT->SetSelectedPad(NULL);

    c.cd(1);
    ((TH2*)&fHOrig)->DrawCopy("");
    gPad->SetLogx();

    c.cd(2);
    ((TH2*)&fHUnfold)->DrawCopy("");
    gPad->SetLogx();

    c.cd(3);
    ((TH2*)&fHFlux)->DrawCopy("");
    gPad->SetLogx();

    c.Modified();
    c.Update();

    return &c;
}

// -------------------------------------------------------------------------
//
// Draw the histograms
//
void MHFlux::Draw(Option_t *opt)
{
    if (!gPad)
        MakeDefCanvas("flux", "orig-unfold-flux plots");

    gPad->Divide(2,2);

    gPad->cd(1);
    fHOrig.Draw(opt);

    gPad->cd(2);
    fHUnfold.Draw(opt);

    gPad->cd(3);
    fHFlux.Draw(opt);

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

// -------------------------------------------------------------------------
//
// Quadratic interpolation
//
// *** calculate the parameters of a parabula 
//                      y = a + b*x + c*x**2 = F(x)
//     such that       yi = F(xi)       for (i=1,3)
//
Bool_t MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
                     Double_t y1, Double_t y2, Double_t y3,
                     Double_t *a, Double_t *b, Double_t *c)
{
    const double det =
        + x2*x3*x3 + x1*x2*x2 + x3*x1*x1
        - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;

    if (det == 0.0)
    {
        *a = 0;
        *b = 0;
        *c = 0;
        return kFALSE;
    }

    const double det1 = 1.0/det;

    const double ai11 = x2*x3*x3 - x3*x2*x2;
    const double ai12 = x3*x1*x1 - x1*x3*x3;
    const double ai13 = x1*x2*x2 - x2*x1*x1;

    const double ai21 = x2*x2 - x3*x3;
    const double ai22 = x3*x3 - x1*x1;
    const double ai23 = x1*x1 - x2*x2;

    const double ai31 = x3 - x2;
    const double ai32 = x1 - x3;
    const double ai33 = x2 - x1;

    *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
    *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1;
    *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;

    //yt1 = *a + *b * x1 + *c * x1*x1;
    //yt2 = *a + *b * x2 + *c * x2*x2;
    //yt3 = *a + *b * x3 + *c * x3*x3;

    //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl;
    //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl;
    //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2
    //	  << ",  yt3 = " << yt3 << endl;
    //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c
    //      << ",  *errflag = " << *errflag << endl;

    return kTRUE;
}
