source: trunk/MagicSoft/Mars/mhist/MH.cc@ 1835

/* ======================================================================== *\
!
! *
! * 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  07/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MH                                                                      //
//                                                                          //
//  This is a base tasks for mars histograms. It defines a common interface //
//  for filling the histograms with events (MH::Fill) which is used by a    //
//  common 'filler' And a SetupFill member function which may be used       //
//  by MFillH. The idea is:                                                 //
//   1) If your Histogram can become filled by one single container         //
//      (like MHHillas) you overload MH::Fill and it gets called with       //
//      a pointer to the container with which it should be filled.          //
//                                                                          //
//   2) You histogram needs several containers to get filled. Than you      //
//      have to overload MH::SetupFill and get the necessary objects from   //
//      the parameter list. Use this objects in Fill to fill your           //
//      histogram.                                                          //
//                                                                          //
//  If you want to create your own histogram class the new class must be    //
//  derived from MH (instead of the base MParContainer) and you must        //
//  the fill function of MH. This is the function which is called to fill   //
//  the histogram(s) by the data of a corresponding parameter container.    //
//                                                                          //
//  Remark: the static member function (eg MakeDefCanvas) can be called     //
//          from everywhere using: MH::MakeDefCanvas(...)                   //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MH.h"

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TGaxis.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TPaveStats.h>

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

#include "MParList.h"
#include "MParContainer.h"

#include "MBinning.h"

ClassImp(MH);

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title only. Typically you won't
// need to change this.
//
MH::MH(const char *name, const char *title)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MH";
    fTitle = title ? title : "Base class for Mars histograms";
}

// --------------------------------------------------------------------------
//
// If you want to use the automatic filling of your derived class you
// must overload this function. If it is not overloaded you cannot use
// FillH with this class. The argument is a pointer to a container
// in your paremeter list which is specified in the MFillH constructor.
// If you are not going to use it you should at least add
//   Bool_t MH::Fill(const MParContainer *) { return kTRUE; }
// to your class definition.
//
Bool_t MH::Fill(const MParContainer *par)
{
    *fLog << warn << GetDescriptor() << ": Fill not overloaded! Can't be used!" << endl;
    return kFALSE;
}

// --------------------------------------------------------------------------
//
// This virtual function is ment as a generalized interface to retrieve
// a pointer to a root histogram from the MH-derived class.
//
TH1 *MH::GetHistByName(const TString name)
{
    *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl;
    return NULL;
}

// --------------------------------------------------------------------------
//
// This is a function which should replace the creation of default
// canvases like root does. Because this is inconvinient in some aspects.
// need to change this.
// You can specify a name for the default canvas and a title. Also
// width and height can be given.
// MakeDefCanvas looks for a canvas with the given name. If now name is
// given the DefCanvasName of root is used. If no such canvas is existing
// it is created and returned. If such a canvas already exists a new canvas
// with a name plus anumber is created (the number is calculated by the
// number of all existing canvases plus one)
//
TCanvas *MH::MakeDefCanvas(TString name, const char *title,
                           const UInt_t w, const UInt_t h)
{
    const TList *list = (TList*)gROOT->GetListOfCanvases();

    if (name.IsNull())
        name = gROOT->GetDefCanvasName();

    if (list->FindObject(name))
        name += Form(" <%d>", list->GetSize()+1);

    return new TCanvas(name, title, w, h);
}

// --------------------------------------------------------------------------
//
// This function works like MakeDefCanvas(name, title, w, h) but name
// and title are retrieved from the given TObject.
//
TCanvas *MH::MakeDefCanvas(const TObject *obj,
                           const UInt_t w, const UInt_t h)
{
    return MakeDefCanvas(obj->GetName(), obj->GetTitle(), w, h);
}

// --------------------------------------------------------------------------
//
// Applies a given binning to a 1D-histogram
//
void MH::SetBinning(TH1 *h, const MBinning *binsx)
{
    //
    // Another strange behaviour: TAxis::Set deletes the axis title!
    //
    TAxis &x = *h->GetXaxis();

#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    TString xtitle = x.GetTitle();
#endif

    //
    // This is a necessary workaround if one wants to set
    // non-equidistant bins after the initialization
    // TH1D::fNcells must be set correctly.
    //
    h->SetBins(binsx->GetNumBins(), 0, 1);

    //
    // Set the binning of the current histogram to the binning
    // in one of the two given histograms
    //
    x.Set(binsx->GetNumBins(), binsx->GetEdges());
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    x.SetTitle(xtitle);
#endif
}

// --------------------------------------------------------------------------
//
// Applies given binnings to the two axis of a 2D-histogram
//
void MH::SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy)
{
    TAxis &x = *h->GetXaxis();
    TAxis &y = *h->GetYaxis();

    //
    // Another strange behaviour: TAxis::Set deletes the axis title!
    //
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    TString xtitle = x.GetTitle();
    TString ytitle = y.GetTitle();
#endif

    //
    // This is a necessary workaround if one wants to set
    // non-equidistant bins after the initialization
    // TH1D::fNcells must be set correctly.
    //
    h->SetBins(binsx->GetNumBins(), 0, 1,
               binsy->GetNumBins(), 0, 1);

    //
    // Set the binning of the current histogram to the binning
    // in one of the two given histograms
    //
    x.Set(binsx->GetNumBins(), binsx->GetEdges());
    y.Set(binsy->GetNumBins(), binsy->GetEdges());
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    x.SetTitle(xtitle);
    y.SetTitle(ytitle);
#endif
}

// --------------------------------------------------------------------------
//
// Applies given binnings to the three axis of a 3D-histogram
//
void MH::SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
{
    //
    // Another strange behaviour: TAxis::Set deletes the axis title!
    //
    TAxis &x = *h->GetXaxis();
    TAxis &y = *h->GetYaxis();
    TAxis &z = *h->GetZaxis();

#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    TString xtitle = x.GetTitle();
    TString ytitle = y.GetTitle();
    TString ztitle = z.GetTitle();
#endif

    //
    // This is a necessary workaround if one wants to set
    // non-equidistant bins after the initialization
    // TH1D::fNcells must be set correctly.
    //
    h->SetBins(binsx->GetNumBins(), 0, 1,
               binsy->GetNumBins(), 0, 1,
               binsz->GetNumBins(), 0, 1);

    //
    // Set the binning of the current histogram to the binning
    // in one of the two given histograms
    //
    x.Set(binsx->GetNumBins(), binsx->GetEdges());
    y.Set(binsy->GetNumBins(), binsy->GetEdges());
    z.Set(binsz->GetNumBins(), binsz->GetEdges());
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    x.SetTitle(xtitle);
    y.SetTitle(ytitle);
    z.SetTitle(ztitle);
#endif
}

// --------------------------------------------------------------------------
//
// Applies given binning (the n+1 edges)  to the axis of a 1D-histogram
//
void MH::SetBinning(TH1 *h, const TArrayD &binsx)
{
    MBinning bx;
    bx.SetEdges(binsx);
    SetBinning(h, &bx);
}

// --------------------------------------------------------------------------
//
// Applies given binning (the n+1 edges) to the two axis of a
// 2D-histogram
//
void MH::SetBinning(TH2 *h, const TArrayD &binsx, const TArrayD &binsy)
{
    MBinning bx;
    MBinning by;
    bx.SetEdges(binsx);
    by.SetEdges(binsy);
    SetBinning(h, &bx, &by);
}

// --------------------------------------------------------------------------
//
// Applies given binning (the n+1 edges) to the three axis of a
// 3D-histogram
//
void MH::SetBinning(TH3 *h, const TArrayD &binsx, const TArrayD &binsy, const TArrayD &binsz)
{
    MBinning bx;
    MBinning by;
    MBinning bz;
    bx.SetEdges(binsx);
    by.SetEdges(binsy);
    bz.SetEdges(binsz);
    SetBinning(h, &bx, &by, &bz);
}

// --------------------------------------------------------------------------
//
// Applies the binning of a TAxis (eg from a root histogram) to the axis
// of a 1D-histogram
//
void MH::SetBinning(TH1 *h, const TAxis *binsx)
{
    const Int_t nx = binsx->GetNbins();

    TArrayD bx(nx+1);
    for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
    bx[nx] = binsx->GetXmax();

    SetBinning(h, bx);
}

// --------------------------------------------------------------------------
//
// Applies the binnings of the TAxis' (eg from a root histogram) to the
// two axis' of a 2D-histogram
//
void MH::SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy)
{
    const Int_t nx = binsx->GetNbins();
    const Int_t ny = binsy->GetNbins();

    TArrayD bx(nx+1);
    TArrayD by(ny+1);
    for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
    for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
    bx[nx] = binsx->GetXmax();
    by[ny] = binsy->GetXmax();

    SetBinning(h, bx, by);
}

// --------------------------------------------------------------------------
//
// Applies the binnings of the TAxis' (eg from a root histogram) to the
// three axis' of a 3D-histogram
//
void MH::SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
{
    const Int_t nx = binsx->GetNbins();
    const Int_t ny = binsy->GetNbins();
    const Int_t nz = binsz->GetNbins();

    TArrayD bx(nx+1);
    TArrayD by(ny+1);
    TArrayD bz(nz+1);
    for (int i=0; i<nx; i++) bx[i] = binsx->GetBinLowEdge(i+1);
    for (int i=0; i<ny; i++) by[i] = binsy->GetBinLowEdge(i+1);
    for (int i=0; i<nz; i++) bz[i] = binsz->GetBinLowEdge(i+1);
    bx[nx] = binsx->GetXmax();
    by[ny] = binsy->GetXmax();
    bz[nz] = binsz->GetXmax();

    SetBinning(h, bx, by, bz);
}

// --------------------------------------------------------------------------
//
// Applies the binnings of one root-histogram x to another one h
// Both histograms must be of the same type: TH1, TH2 or TH3
//
void MH::SetBinning(TH1 *h, TH1 *x)
{
    if (h->InheritsFrom(TH3::Class()) && x->InheritsFrom(TH3::Class()))
    {
        SetBinning((TH3*)h, x->GetXaxis(), x->GetYaxis(), x->GetZaxis());
        return;
    }
    if (h->InheritsFrom(TH3::Class()) || x->InheritsFrom(TH3::Class()))
        return;
    if (h->InheritsFrom(TH2::Class()) && x->InheritsFrom(TH2::Class()))
    {
        SetBinning((TH2*)h, x->GetXaxis(), x->GetYaxis());
        return;
    }
    if (h->InheritsFrom(TH2::Class()) || x->InheritsFrom(TH2::Class()))
        return;
    if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
    {
        SetBinning(h, x->GetXaxis());
        return;
    }
}

// --------------------------------------------------------------------------
//
// Multiplies all entries in a TArrayD by a float f
//
void MH::ScaleArray(TArrayD &bins, Double_t f)
{
    for (int i=0; i<bins.GetSize(); i++)
        bins[i] *= f;
}

// --------------------------------------------------------------------------
//
// Scales the binning of a TAxis by a float f
//
TArrayD MH::ScaleAxis(TAxis &axe, Double_t f)
{
    TArrayD arr(axe.GetNbins()+1);

    for (int i=1; i<=axe.GetNbins()+1; i++)
        arr[i-1] = axe.GetBinLowEdge(i);

    ScaleArray(arr, f);

    return arr;
}

// --------------------------------------------------------------------------
//
// Scales the binning of one, two or three axis of a histogram by a float f
//
void MH::ScaleAxis(TH1 *h, Double_t fx, Double_t fy, Double_t fz)
{
    if (h->InheritsFrom(TH3::Class()))
    {
        SetBinning((TH3*)h,
                   ScaleAxis(*h->GetXaxis(), fx),
                   ScaleAxis(*h->GetYaxis(), fy),
                   ScaleAxis(*h->GetZaxis(), fz));
        return;
    }

    if (h->InheritsFrom(TH2::Class()))
    {
        SetBinning((TH2*)h,
                   ScaleAxis(*h->GetXaxis(), fx),
                   ScaleAxis(*h->GetYaxis(), fy));
        return;
    }

    if (h->InheritsFrom(TH1::Class()))
        SetBinning(h, ScaleAxis(*h->GetXaxis(), fx));
}

// --------------------------------------------------------------------------
//
// Tries to find a MBinning container with the name "Binning"+name
// in the given parameter list. If it was found it is applied to the
// given histogram. This is only valid for 1D-histograms
//
Bool_t MH::ApplyBinning(const MParList &plist, TString name, TH1 *h)
{
    if (h->InheritsFrom(TH2::Class()) || h->InheritsFrom(TH3::Class()))
    {
        gLog << warn << "MH::ApplyBinning: '" << h->GetName() << "' is not a basic TH1 object... no binning applied." << endl;
        return kFALSE;
    }

    const MBinning *bins = (MBinning*)plist.FindObject("Binning"+name);
    if (!bins)
    {
        gLog << warn << "Object 'Binning" << name << "' [MBinning] not found... no binning applied." << endl;
        return kFALSE;
    }

    SetBinning(h, bins);
    return kTRUE;
}

void MH::FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger)
{
//*-*-*-*-*-*-*-*-*Find reasonable bin values*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*              ==========================

    Double_t dx = 0.1*(xmax-xmin);
    Double_t umin = xmin - dx;
    Double_t umax = xmax + dx;

    if (umin < 0 && xmin >= 0)
        umin = 0;

    if (umax > 0 && xmax <= 0)
        umax = 0;

    Double_t binlow  =0;
    Double_t binhigh =0;
    Double_t binwidth=0;

#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,01)
    TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, nbins, binwidth, "");
#else
    gLog << all << "*********************************************************" << endl;
    gLog << all << "ERROR, because Abelardo removed TGaxis::Optimize from" << endl;
    gLog << all << " MH::FindGoodLimits to be able to compile Mars with newer" << endl;
    gLog << all << " root versions, whatever you are trying to do will fail!!!" << endl;
    gLog << all << "*********************************************************" << endl;
#endif

    if (binwidth <= 0 || binwidth > 1.e+39)
    {
        xmin = -1;
        xmax = 1;
    }
    else
    {
        xmin = binlow;
        xmax = binhigh;
    }

    if (isInteger)
    {
        Int_t ixmin = (Int_t)xmin;
        Int_t ixmax = (Int_t)xmax;
        Double_t dxmin = (Double_t)ixmin;
        Double_t dxmax = (Double_t)ixmax;

        xmin = xmin<0 && xmin!=dxmin ? dxmin - 1 : dxmin;
        xmax = xmax>0 && xmax!=dxmax ? dxmax + 1 : dxmax;

        if (xmin>=xmax)
            xmax = xmin+1;

        Int_t bw = 1 + (Int_t)((xmax-xmin)/nbins);

        nbins = (Int_t)((xmax-xmin)/bw);

        if (xmin+nbins*bw < xmax)
        {
            nbins++;
            xmax = xmin +nbins*bw;
        }
    }

    newbins = nbins;
}

// --------------------------------------------------------------------------
//
//  Returns the lowest entry in a histogram which is greater than gt (eg >0)
//
Double_t MH::GetMinimumGT(const TH1 &h, Double_t gt)
{
    Double_t min = FLT_MAX;

    const TAxis &axe = *((TH1&)h).GetXaxis();

    for (int i=1; i<=axe.GetNbins(); i++)
    {
        Double_t x = h.GetBinContent(i);
        if (gt<x && x<min)
            min = x;
    }
    return min;
}

// --------------------------------------------------------------------------
//
//  Returns the bin center in a logarithmic scale. If the given bin
//  number is <1 it is set to 1. If it is =GetNbins() it is set to
//  GetNbins()
//
Double_t MH::GetBinCenterLog(const TAxis &axe, Int_t nbin)
{
    if (nbin>axe.GetNbins())
        nbin = axe.GetNbins();

    if (nbin<1)
        nbin = 1;

    const Double_t lo = axe.GetBinLowEdge(nbin);
    const Double_t hi = axe.GetBinUpEdge(nbin);

    const Double_t val = log10(lo) + log10(hi);

    return pow(10, val/2);
}

// --------------------------------------------------------------------------
//
// Draws a copy of the two given histograms. Uses title as the pad title.
// Also layout the two statistic boxes and a legend.
//
void MH::DrawCopy(const TH1 &hist1, const TH1 &hist2, const TString title)
{
    //
    // Draw first histogram
    //
    TH1 *h1 = (TH1*)((TH1&)hist1).DrawCopy();
    gPad->Update();

    TPaveText *t = (TPaveText*)gPad->FindObject("title");
    if (t)
    {
        t->SetName((TString)"MHTitle");     // rename object
        t->Clear();                         // clear old lines
        t->AddText((TString)" "+title+" "); // add the new title
        t->SetBit(kCanDelete);              // make sure object is deleted

        //
        // FIXME: This is a stupid workaround to hide the redrawn
        // (see THistPainter::PaintTitle) title
        //
        gPad->Modified();  // indicates a change
        gPad->Update();    // recreates the original title
        t->Pop();          // bring our title on top
    }

    //
    // Rename first statistics box
    //
    TPaveStats &s1 = *(TPaveStats*)gPad->FindObject("stats");
    s1.SetX1NDC(s1.GetX1NDC()-0.01);
    s1.SetName("MHStat");

    //
    // Draw second histogram
    //
    TH1 *h2 = (TH1*)((TH1&)hist2).DrawCopy("sames");
    gPad->Update();

    //
    // Set new position of second statistics box
    //
    TPaveStats &s2 = *(TPaveStats*)gPad->FindObject("stats");
    s2.SetX1NDC(s1.GetX1NDC()-(s2.GetX2NDC()-s2.GetX1NDC())-0.01);
    s2.SetX2NDC(s1.GetX1NDC()-0.01);

    //
    // Draw Legend
    //
    const Int_t    n = s1.GetListOfLines()->GetSize();
    const Double_t h = s1.GetY2NDC()-s1.GetY1NDC();
    TLegend &l = *new TLegend(s1.GetX1NDC(), s1.GetY1NDC()-0.015-h*2/n,
                              s1.GetX2NDC(), s1.GetY1NDC()-0.01
                            );
    l.AddEntry(h1, h1->GetTitle());
    l.AddEntry(h2, h2->GetTitle());
    l.SetTextSize(s2.GetTextSize());
    l.SetTextFont(s2.GetTextFont());
    l.SetBorderSize(s2.GetBorderSize());
    l.SetBit(kCanDelete);

    l.Draw();
}

// --------------------------------------------------------------------------
//
// Draws the two given histograms. Uses title as the pad title.
// Also layout the two statistic boxes and a legend.
//
void MH::Draw(TH1 &hist1, TH1 &hist2, const TString title)
{
    //
    // Draw first histogram
    //
    hist1.Draw();
    gPad->Update();

    TPaveText *t = (TPaveText*)gPad->FindObject("title");
    if (t)
    {
        t->SetName((TString)"MHTitle");     // rename object
        t->Clear();                         // clear old lines
        t->AddText((TString)" "+title+" "); // add the new title
        t->SetBit(kCanDelete);              // make sure object is deleted

        //
        // FIXME: This is a stupid workaround to hide the redrawn
        // (see THistPainter::PaintTitle) title
        //
        gPad->Modified();  // indicates a change
        gPad->Update();    // recreates the original title
        t->Pop();          // bring our title on top
    }

    //
    // Rename first statistics box
    //
    TPaveStats &s1 = *(TPaveStats*)gPad->FindObject("stats");
    s1.SetX1NDC(s1.GetX1NDC()-0.01);
    s1.SetName((TString)"Stat"+hist1.GetTitle());

    //
    // Draw second histogram
    //
    hist2.Draw("sames");

    gPad->Update();

    //
    // Set new position of second statistics box
    //
    TPaveStats &s2 = *(TPaveStats*)gPad->FindObject("stats");
    s2.SetX1NDC(s1.GetX1NDC()-(s2.GetX2NDC()-s2.GetX1NDC())-0.01);
    s2.SetX2NDC(s1.GetX1NDC()-0.01);

    //
    // Draw Legend
    //
    TLegend &l = *new TLegend(s1.GetX1NDC(),
                              s1.GetY1NDC()-0.015-(s1.GetY2NDC()-s1.GetY1NDC())/2,
                              s1.GetX2NDC(),
                              s1.GetY1NDC()-0.01
                             );
    l.AddEntry(&hist1, hist1.GetTitle());
    l.AddEntry(&hist2, hist2.GetTitle());
    l.SetTextSize(s2.GetTextSize());
    l.SetTextFont(s2.GetTextFont());
    l.SetBorderSize(s2.GetBorderSize());
    l.SetBit(kCanDelete);
    l.Draw();
}