/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHAlpha
//
// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
// check the result when it is filles in the MStatusDisplay
// For more information see MHFalseSource::FitSignificance
//
// For convinience (fit) the output significance is stored in a
// container in the parlisrt
//
// PRELIMINARY!
//
//////////////////////////////////////////////////////////////////////////////
#include "MHAlpha.h"

#include <TF1.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TPaveStats.h>
#include <TStopwatch.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MHillasSrc.h"
#include "MEnergyEst.h"
#include "MTime.h"
#include "MObservatory.h"
#include "MPointingPos.h"
#include "MAstroSky2Local.h"
#include "MStatusDisplay.h"
#include "MParameters.h"
#include "MHMatrix.h"

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

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

ClassImp(MHAlpha);
ClassImp(MAlphaFitter);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHAlpha::MHAlpha(const char *name, const char *title)
    : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
    fPointPos(0), fTimeEffOn(0), fTime(0),
    fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE),
    //fEnergyMin(-1), fEnergyMax(-1), fSizeMin(-1), fSizeMax(-1),
    fMatrix(0)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHAlpha";
    fTitle = title ? title : "Alpha plot";

    fHAlpha.SetName("Alpha");
    fHAlpha.SetTitle("Alpha");
    fHAlpha.SetXTitle("\\Theta [deg]");
    fHAlpha.SetYTitle("E_{est} [GeV]");
    fHAlpha.SetZTitle("|\\alpha| [\\circ]");
    fHAlpha.SetDirectory(NULL);
    fHAlpha.UseCurrentStyle();

    // Main histogram
    fHAlphaTime.SetName("Alpha");
    fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
    fHAlphaTime.SetYTitle("Counts");
    fHAlphaTime.UseCurrentStyle();
    fHAlphaTime.SetDirectory(NULL);


    //fHEnergy.SetName("Energy");
    //fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
    //fHEnergy.SetXTitle("E_{est} [GeV]");
    fHEnergy.SetYTitle("N_{exc}");
    fHEnergy.SetDirectory(NULL);
    fHEnergy.UseCurrentStyle();

    fHTheta.SetName("Theta");
    fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
    fHTheta.SetXTitle("\\Theta [\\circ]");
    fHTheta.SetYTitle("N_{exc}");
    fHTheta.SetDirectory(NULL);
    fHTheta.UseCurrentStyle();

    // effective on time versus time
    fHTime.SetName("Time");
    fHTime.SetTitle(" N_{exc} vs. Time ");
    fHTime.SetXTitle("Time");
    fHTime.SetYTitle("N_{exc} [s]");
    fHTime.UseCurrentStyle();
    fHTime.SetDirectory(NULL);
    fHTime.GetYaxis()->SetTitleOffset(1.2);
    fHTime.GetXaxis()->SetLabelSize(0.033);
    fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    fHTime.GetXaxis()->SetTimeDisplay(1);
    fHTime.Sumw2();

    MBinning binsa, binse, binst;
    binsa.SetEdges(18, 0, 90);
    binse.SetEdgesLog(25, 10, 100000);
    binst.SetEdgesCos(50, 0, 60);
    binse.Apply(fHEnergy);
    binst.Apply(fHTheta);
    binsa.Apply(fHAlphaTime);

    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
}

void MHAlpha::FitEnergySpec(Bool_t paint)
{
    /*
    TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
    f.SetParameter(0,  -2.1);
    f.SetParameter(1,  1400); // 50
    f.SetParameter(2,  fHEnergy.Integral()); // 50
    f.SetLineWidth(2);

    fHEnergy.Fit(&f,  "0QI", "", 100, 2400); // Integral? 90, 2800

    const Double_t chi2 = f.GetChisquare();
    const Int_t    NDF  = f.GetNDF();

    // was fit successful ?
    const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    f.SetLineColor(ok ? kGreen : kRed);

    if (paint)
    {
        f.Paint("same");

        if (ok)
        {
        TString s;
        s += Form("\\alpha=%.1f ",   f.GetParameter(0));
        s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
        s += Form("N=%.1e",          f.GetParameter(2));

        TLatex text(0.25, 0.94, s.Data());
        text.SetBit(TLatex::kTextNDC);
        text.SetTextSize(0.04);
        text.Paint();
        }
    }*/
}

void MHAlpha::FitEnergyBins(Bool_t paint)
{
    // Do not store the 'final' result in fFit
    MAlphaFitter fit(fFit);

    const Int_t n = fHAlpha.GetNbinsY();

    fHEnergy.SetEntries(0);

    for (int i=1; i<=n; i++)
    {
        if (fit.FitEnergy(fHAlpha, fOffData, i))
        {
            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
            fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
        }
    }
}

void MHAlpha::FitThetaBins(Bool_t paint)
{
    // Do not store the 'final' result in fFit
    MAlphaFitter fit(fFit);

    const Int_t n = fHAlpha.GetNbinsX();

    fHTheta.SetEntries(0);

    for (int i=1; i<=n; i++)
    {
        if (fit.FitTheta(fHAlpha, fOffData, i))
        {
            fHTheta.SetBinContent(i, fit.GetEventsExcess());
            fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
        }
    }
}

Bool_t MHAlpha::SetupFill(const MParList *pl)
{
    fHAlpha.Reset();
    fHEnergy.Reset();
    fHTheta.Reset();
    fHTime.Reset();

    if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
    {
        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha");
        if (!hoff)
            *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
        else
        {
            *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
            SetOffData(*hoff);
        }
    }

    fHillas = 0;
    /*
    if (fSizeMin>=0 || fSizeMax>=0)
    {
        fHillas = (MHillas*)pl->FindObject("MHillas");
        if (!fHillas)
        {
            *fLog << warn << "Size cut set, but MHillas not found... abort." << endl;
            return kFALSE;
        }
    }
    */
    fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
    if (!fEnergy)
    { /*
        if (fEnergyMin>=0 || fEnergyMax>=0)
        {
            *fLog << warn << "Energy cut set, but MEnergyEst not found... abort." << endl;
            return kFALSE;
        } */

        *fLog << warn << "MEnergyEst not found... " << flush;

        if (!fHillas)
            fHillas = (MHillas*)pl->FindObject("MHillas");
        if (fHillas)
            *fLog << "using SIZE instead!" << endl;
        else
            *fLog << "ignored." << endl;

        fHEnergy.SetName("Size");
        fHEnergy.SetTitle(" N_{exc} vs. Size ");
        fHEnergy.SetXTitle("Size [\\gamma]");
    }
    else
    {
        fHEnergy.SetName("Energy");
        fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
        fHEnergy.SetXTitle("E_{est} [GeV]");
    }

    fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
    if (!fPointPos)
        *fLog << warn << "MPointingPos not found... ignored." << endl;

    fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
    if (!fTimeEffOn)
        *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
    else
        *fTimeEffOn = MTime(); // FIXME: How to do it different?

    fTime = (MTime*)pl->FindObject("MTime");
    if (!fTime)
        *fLog << warn << "MTime not found... ignored." << endl;

    fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MinimizationValue");
    if (!fResult)
        return kFALSE;

    //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
    //if (!fExcess)
    //    return kFALSE;

    fLastTime = MTime();

    MBinning binst, binse, binsa;
    binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
    binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
    binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');

    if (!fOffData)
    {
        if (!fPointPos)
            binst.SetEdges(1, 0, 60);
        else
            binst.SetEdges(*pl, "BinningTheta");

        if (!fEnergy && !fHillas)
            binse.SetEdges(1, 10, 100000);
        else
            binse.SetEdges(*pl, "BinningEnergyEst");

        binsa.SetEdges(*pl, "BinningAlpha");
    }

    binse.Apply(fHEnergy);
    binst.Apply(fHTheta);
    binsa.Apply(fHAlphaTime);
    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);

    MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
    if (!fit)
        *fLog << warn << "MAlphaFitter not found... ignored." << endl;
    else
        fit->Copy(fFit);

    fFit.Print();

    return kTRUE;
}

void MHAlpha::InitAlphaTime(const MTime &t)
{
    //
    // If this is the first call we have to initialize the time-histogram
    //
    MBinning bins;
    bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
    bins.Apply(fHTime);

    fLastTime=t;
}

void MHAlpha::UpdateAlphaTime(Bool_t final)
{
    if (!fTimeEffOn)
        return;

    const Int_t steps = 10;

    static int rebin = steps;

    if (!final)
    {
        if (fLastTime==MTime() && *fTimeEffOn!=MTime())
        {
            InitAlphaTime(*fTimeEffOn);
            return;
        }

        if (fLastTime!=*fTimeEffOn)
        {
            fLastTime=*fTimeEffOn;
            rebin--;
        }

        if (rebin>0)
            return;
    }

    // Check new 'last time'
    MTime *time = final ? fTime : fTimeEffOn;

    if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
    {
        *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl;
        *fLog << "than upper edge of histogram... skipped." << endl;
        *fLog << "This should not happen. Maybe you started you eventloop with" << endl;
        *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
        rebin++;
        return;
    }

    // Fit time histogram
    MAlphaFitter fit(fFit);

    TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
    const Bool_t rc = fit.ScaleAndFit(fHAlphaTime, h);
    if (h)
        delete h;

    if (!rc)
        return;

    // Reset Histogram
    fHAlphaTime.Reset();
    fHAlphaTime.SetEntries(0);

    // Get number of bins
    const Int_t n = fHTime.GetNbinsX();

    //
    // Prepare Histogram
    //
    if (final)
        time->Plus1ns();

    // Enhance binning
    MBinning bins;
    bins.SetEdges(fHTime, 'x');
    bins.AddEdge(time->GetAxisTime());
    bins.Apply(fHTime);

    //
    // Fill histogram
    //
    fHTime.SetBinContent(n+1, fit.GetEventsExcess());
    fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);

    *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;

    rebin = steps;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
{
    Double_t alpha, energy, theta;
    Double_t size=-1;

    if (fMatrix)
    {
        alpha  = (*fMatrix)[fMap[0]];
        energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
        size   = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
        //<0 ? 1000 : (*fMatrix)[fMap[1]];
        theta  = 0;

        if (energy<0)
            energy=size;
        if (size<0)
            size=energy;

        if (energy<0 && size<0)
            energy = size = 1000;
    }
    else
    {
        const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
        if (!par)
        {
            *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
            return kFALSE;
        }

        alpha  = hil->GetAlpha();
        if (fHillas)
            size = fHillas->GetSize();
        energy = fEnergy   ? fEnergy->GetEnergy() : (fHillas?fHillas->GetSize():1000);
        theta  = fPointPos ? fPointPos->GetZd()   : 0;
    }

    // enhance histogram if necessary
    UpdateAlphaTime();

    // Fill histograms
    fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);

    // Check cuts
    /*
    if ( (fEnergyMin>=0 && energy<fEnergyMin) ||
         (fEnergyMax>=0 && energy>fEnergyMax) ||
         (fSizeMin>=0   && size  <fSizeMin)   ||
         (fSizeMax>=0   && size  >fSizeMin) )
         return kTRUE;
         */

    if (!fSkipHistTime)
        fHAlphaTime.Fill(TMath::Abs(alpha), w);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Paint the integral and the error on top of the histogram
//
void MHAlpha::PaintText(Double_t val, Double_t error) const
{
    TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
    text.SetBit(TLatex::kTextNDC);
    text.SetTextSize(0.04);
    text.Paint();
}

// --------------------------------------------------------------------------
//
// Update the projections
//
void MHAlpha::Paint(Option_t *opt)
{
    // Note: Any object cannot be added to a pad twice!
    //       The object is searched by gROOT->FindObject only in
    //       gPad itself!
    TVirtualPad *padsave = gPad;

    TH1D *h0=0;

    TString o(opt);
    if (o==(TString)"proj")
    {
        TPaveStats *st=0;
        for (int x=0; x<4; x++)
        {
            TVirtualPad *p=padsave->GetPad(x+1);
            if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
                continue;

            if (st->GetOptStat()==11)
                continue;

            const Double_t y1 = st->GetY1NDC();
            const Double_t y2 = st->GetY2NDC();
            const Double_t x1 = st->GetX1NDC();
            const Double_t x2 = st->GetX2NDC();

            st->SetY1NDC((y2-y1)/3+y1);
            st->SetX1NDC((x2-x1)/3+x1);
            st->SetOptStat(11);
        }

        padsave->cd(1);
        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
        if (fOffData)
        {
            TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
            const Double_t alpha = fFit.Scale(*hoff, *hon);

            hon->SetMaximum();
            hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);

            // BE CARFEULL: This is a really weird workaround!
            hoff->SetMaximum(alpha);

            if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
            {
                h0->Reset();
                h0->Add(hoff, hon, -1);
                const Float_t min = h0->GetMinimum()*1.05;
                hon->SetMinimum(min<0 ? min : 0);
            }
        }
        FitEnergyBins();
        FitThetaBins();
    }

    if (o==(TString)"alpha")
        if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
        {
            // Check whether we are dealing with on-off analysis
            TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff");

            // BE CARFEULL: This is a really weird workaround!
            const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();

            // Do not store the 'final' result in fFit
            MAlphaFitter fit(fFit);
            fit.Fit(*h0, hoff, scale, kTRUE);
            fit.PaintResult();
        }

    if (o==(TString)"time")
        PaintText(fHTime.Integral(), 0);

    if (o==(TString)"theta")
    {
        TH1 *h = (TH1*)gPad->FindObject("Alpha_x");
        if (h)
        {
            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x");
            h2->SetDirectory(0);
            h2->Scale(fHTheta.Integral()/h2->Integral());
            h->Reset();
            h->Add(h2);
            delete h2;
        }
        PaintText(fHTheta.Integral(), 0);
    }

    if (o==(TString)"energy")
    {
        TH1 *h = (TH1*)gPad->FindObject("Alpha_y");
        if (h)
        {
            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y");
            h2->SetDirectory(0);
            h2->Scale(fHEnergy.Integral()/h2->Integral());
            h->Reset();
            h->Add(h2);
            delete h2;
        }

        if (fHEnergy.GetMaximum()>1)
        {
            gPad->SetLogx();
            gPad->SetLogy();
        }
        FitEnergySpec(kTRUE);
    }

    gPad = padsave;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHAlpha::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);

    // Do the projection before painting the histograms into
    // the individual pads
    AppendPad("proj");

    pad->SetBorderMode(0);
    pad->Divide(2,2);

    TH1D *h=0;

    pad->cd(1);
    gPad->SetBorderMode(0);

    h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
    h->SetBit(TH1::kNoTitle);
    h->SetStats(kTRUE);
    h->SetXTitle("\\alpha [\\circ]");
    h->SetYTitle("Counts");
    h->SetDirectory(NULL);
    h->SetMarkerStyle(kFullDotMedium);
    h->SetBit(kCanDelete);
    h->Draw("");

    if (fOffData)
    {
        h->SetMarkerColor(kGreen);

        h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
        h->SetBit(TH1::kNoTitle);
        h->SetXTitle("\\alpha [\\circ]");
        h->SetYTitle("Counts");
        h->SetDirectory(NULL);
        h->SetMarkerStyle(kFullDotMedium);
        h->SetBit(kCanDelete);
        h->SetMarkerColor(kRed);
        h->Draw("same");

        h = (TH1D*)h->Clone("ProjAlphaOnOff");
        h->SetBit(TH1::kNoTitle);
        h->SetXTitle("\\alpha [\\circ]");
        h->SetYTitle("Counts");
        h->SetDirectory(NULL);
        h->SetMarkerStyle(kFullDotMedium);
        h->SetBit(kCanDelete);
        h->SetMarkerColor(kBlue);
        h->Draw("same");

        TLine lin;
        lin.SetLineStyle(kDashed);
        lin.DrawLine(0, 0, 90, 0);
    }

    // After the Alpha-projection has been drawn. Fit the histogram
    // and paint the result into this pad
    AppendPad("alpha");

    if (fHEnergy.GetNbinsX()>1)
    {
        pad->cd(2);
        gPad->SetBorderMode(0);
        fHEnergy.Draw();

        AppendPad("energy");

        h = (TH1D*)fHAlpha.Project3D("y");
        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
        h->SetXTitle("E [GeV]");
        h->SetYTitle("Counts");
        h->SetDirectory(NULL);
        h->SetMarkerStyle(kFullDotMedium);
        h->SetBit(kCanDelete);
        h->SetMarkerColor(kCyan);
        h->SetLineColor(kCyan);
        h->Draw("Psame");
    }
    else
        delete pad->GetPad(2);

    if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
    {
        pad->cd(3);
        gPad->SetBorderMode(0);
        fHTime.Draw();
        AppendPad("time");
    }
    else
        delete pad->GetPad(3);

    if (fHTheta.GetNbinsX()>1)
    {
        pad->cd(4);
        gPad->SetBorderMode(0);
        fHTheta.Draw();

        AppendPad("theta");

        h = (TH1D*)fHAlpha.Project3D("x");
        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
        h->SetXTitle("\\theta [\\circ]");
        h->SetYTitle("Counts");
        h->SetDirectory(NULL);
        h->SetMarkerStyle(kFullDotMedium);
        h->SetBit(kCanDelete);
        h->SetMarkerColor(kCyan);
        h->SetLineColor(kCyan);
        h->Draw("Psame");
    }
    else
        delete pad->GetPad(4);
}

void MHAlpha::DrawAll()
{
    // FIXME: Do in Paint if special option given!
    TCanvas *c = new TCanvas;
    Int_t n = fHAlpha.GetNbinsY();
    Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
    c->Divide(nc, nc, 0, 0);

    // Do not store the 'final' result in fFit
    MAlphaFitter fit(fFit);

    for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
    {
        c->cd(i);

        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
        hon->SetBit(TH1::kNoTitle);
        hon->SetStats(kTRUE);
        hon->SetXTitle("\\alpha [\\circ]");
        hon->SetYTitle("Counts");
        hon->SetDirectory(NULL);
        hon->SetMarkerStyle(kFullDotMedium);
        hon->SetBit(kCanDelete);
        hon->Draw("");

        TH1D *hof = 0;
        Double_t alpha = 1;

        if (fOffData)
        {
            hon->SetMarkerColor(kGreen);

            hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
            hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
            hof->SetXTitle("\\alpha [\\circ]");
            hof->SetYTitle("Counts");
            hof->SetDirectory(NULL);
            hof->SetMarkerStyle(kFullDotMedium);
            hof->SetBit(kCanDelete);
            hof->SetMarkerColor(kRed);
            hof->Draw("same");

            alpha = fit.Scale(*hof, *hon);

            hon->SetMaximum();
            hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);

            TH1D *diff = new TH1D(*hon);
            diff->Add(hof, -1);
            diff->SetBit(TH1::kNoTitle);
            diff->SetXTitle("\\alpha [\\circ]");
            diff->SetYTitle("Counts");
            diff->SetDirectory(NULL);
            diff->SetMarkerStyle(kFullDotMedium);
            diff->SetBit(kCanDelete);
            diff->SetMarkerColor(kBlue);
            diff->Draw("same");

            TLine lin;
            lin.SetLineStyle(kDashed);
            lin.DrawLine(0, 0, 90, 0);

            const Float_t min = diff->GetMinimum()*1.05;
            hon->SetMinimum(min<0 ? min : 0);
        }

        if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
            *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
        /*
        if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
        {
            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
            fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
        }*/
    }

}


Bool_t MHAlpha::Finalize()
{
    //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
    //h->SetDirectory(0);
    //Bool_t rc = fFit.Fit(*h);
    //delete h;
    if (!fFit.FitAlpha(fHAlpha, fOffData))
    {
        *fLog << warn << "Histogram empty." << endl;
        return kTRUE;
    }

    // Store the final result in fFit
    fFit.Print("result");

    fResult->SetVal(-fFit.GetSignificance());

    if (!fSkipHistEnergy)
    {
        *fLog << inf << "Processing energy bins..." << endl;
        FitEnergyBins();
    }
    if (!fSkipHistTheta)
    {
        *fLog << inf << "Processing theta bins..." << endl;
        FitThetaBins();
    }
    if (!fSkipHistTime)
    {
        *fLog << inf << "Processing time bins..." << endl;
        UpdateAlphaTime(kTRUE);
        MH::RemoveFirstBin(fHTime);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of
// MMcEvt. This function adds all necessary columns to the
// given matrix. Afterward you should fill the matrix with the corresponding
// data (eg from a file by using MHMatrix::Fill). If you now loop
// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
// will take the values from the matrix instead of the containers.
//
// It takes fSkipHist* into account!
//
void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
{
    if (fMatrix)
        return;

    fMatrix = mat;

    fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
    fMap[1] = -1;
    fMap[2] = -1;
    fMap[3] = -1;
    fMap[4] = -1;

    if (!fSkipHistEnergy)
        if (type==0)
        {
            fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
            fMap[2] = -1;
        }
        else
        {
            fMap[1] = -1;
            fMap[2] = fMatrix->AddColumn("MHillas.fSize");
        }

    if (!fSkipHistTheta)
        fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");

   // if (!fSkipHistTime)
   //     fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
}

void MHAlpha::StopMapping()
{
    fMatrix = NULL; 
}