/* ======================================================================== *\
!
! *
! * 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 char *varname, const char *unit)
 : fHOrig(), fHUnfold(), fHFlux()
{
    if (varname == NULL  ||  unit == NULL)
    {
      *fLog << "MHFlux : 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");

      char txt0[100];
      strg0  = fVarname;
      strg0 += "-bin %d";
      sprintf(txt0, strg0, n);

      TString strg1("No.of photons vs. E-est for ");
      strg1 += txt0;
      new TCanvas(txt0,strg1);
      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
      // gROOT->SetSelectedPad(NULL);

      gPad->SetLogx();

      h.SetName(txt0);
      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);
// 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");

      char txt0[100];
      strg0  = fVarname;
      strg0 += "-bin %d";
      sprintf(txt0, strg0, n);

      TString strg1("No.of photons vs. E-unfold for ");
      strg1 += txt0;
      new TCanvas(txt0,strg1);

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

      gPad->SetLogx();

      h.SetName(txt0);
      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

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

  Int_t    errflag;
  Double_t c0, c1, c2;  
  Double_t t1, t2, t3;
  Double_t a1, a2, a3;

  Double_t *aeffbar;
  aeffbar = new Double_t[nEtru];
  Double_t *daeffbar;
  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, k1=0, k2=0, k3=0;
    for (int k=3; k<=nTheta; k++)
    {
      Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->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 <  ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
      k0 = 1;
    else if (Thetabar >  ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta))
      k0 = nTheta;

    Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1);   
    Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->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++)
    {
      c0 = 0.0;
      c1 = 0.0;
      c2 = 0.0;

      t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) );
      t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) );
      t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) );

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

      Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
      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, j1=0, j2=0, j3=0;

      for (int j=3; j<=nEtru; j++)
      {
        Double_t Elow = ((TH1*)aeff)->GetXaxis()->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 <  ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2))
        j0 = 1;
      else if (Ebar >  ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru))
        j0 = nEtru;

      Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1);   
      Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->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;


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

      t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1)  )
                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) );

      t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2)  )
                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) );

      t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3)  )
                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) );


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

      Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
      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.0  &&  DeltaE > 0.0  &&  Effon > 0.0  &&  Aeff > 0.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");

      char txt[100];
      TString strg1("Photon flux vs. E-unfold for ");
      TString strg2 = fVarname;
      strg2 += "-bin %d";
      sprintf(txt, strg2, n);
      TString strg3 = strg1 + txt;
  
      new TCanvas(txt, strg3);
      // TCanvas &c = *MakeDefCanvas(txt, txt);
      // gROOT->SetSelectedPad(NULL);

      gPad->SetLogx();

      h.SetName(txt);
      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)
//
void 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, Int_t *errflag)
{
  double  ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
  double  det,det1;
  //double  yt1,yt2,yt3;

  det =   x2*x3*x3 + x1*x2*x2 + x3*x1*x1
        - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;

  if (det != 0.0)
  {
    *errflag = 0;
    det1 = 1.0/det;

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

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

    ai31 = x3 - x2;
    ai32 = x1 - x3;
    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;
  }

  *errflag = 1;
  *a = 0.0;
  *b = 0.0;
  *c = 0.0;
  return;
}
