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

/* ======================================================================== *\
!
! *
! * 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.    //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MH.h"

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TGaxis.h>
#include <TCanvas.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 << GetDescriptor() << ": Fill not overloaded! Can't be used!" << endl;
    return kFALSE;
}

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

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
}

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
}

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
}

void MH::SetBinning(TH1 *h, const TArrayD &binsx)
{
    MBinning bx;
    bx.SetEdges(binsx);
    SetBinning(h, &bx);
}

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);
}

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);
}

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);
}

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);
}

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);
}

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(TH2::Class()) && x->InheritsFrom(TH2::Class()))
    {
        SetBinning((TH2*)h, x->GetXaxis(), x->GetYaxis());
        return;
    }
    if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
    {
        SetBinning(h, x->GetXaxis());
        return;
    }
}

void MH::ScaleArray(TArrayD &bins, Double_t f)
{
    for (int i=0; i<bins.GetSize(); i++)
        bins[i] *= 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;
}

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));
}

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;

    Int_t n=0;
    Double_t binlow  =0;
    Double_t binhigh =0;
    Double_t binwidth=0;
    TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, n, binwidth, "");

    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;
}

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;
}