source: trunk/MagicSoft/Mars/mhbase/MBinning.cc@ 5687

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

//////////////////////////////////////////////////////////////////////////////
//
//  MBinning
//
// This is a MParCOntainer storing a binning for a histogram. Doing this
// you are able to distribute a single binning to several histograms
// in your parameter list.
//
// In some classes the title of the container is used to set the
// axis-title of the corresponding axis in your histogram.
//
// For all the features supported see the function descriptions in
//  MBinning and MH
//
//////////////////////////////////////////////////////////////////////////////
#include "MBinning.h"

#include <ctype.h>      // tolower
#include <fstream>

#include <TH1.h>        // InheritsFrom

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

#include "MH.h"

ClassImp(MBinning);

using namespace std;

const TString MBinning::gsDefName  = "MBinning";
const TString MBinning::gsDefTitle = "Container describing the binning of an axis";

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title only. Typically you won't
// need to change this.
//
MBinning::MBinning(const char *name, const char *title)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    SetEdges(10, 0, 1);
    fType = kIsDefault;
}

// --------------------------------------------------------------------------
//
// Instantiate MBinning with nbins number of bins between lo (lower edge)
// and hi (upper edge), name name and title title.
//
MBinning::MBinning(Int_t nbins, Axis_t lo, Axis_t hi, const char *name, const char *opt, const char *title)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    SetEdges(nbins, lo, hi, opt);
}

// --------------------------------------------------------------------------
//
// Setup the edges stored in MBinning from the TAxis axe
//
void MBinning::SetEdges(const TAxis &axe)
{
    const TArrayD &arr = *axe.GetXbins();
    if (arr.GetSize()>0)
    {
        SetEdges(arr);
        return;
    }

    SetEdges(axe.GetNbins(), axe.GetXmin(), axe.GetXmax());
}

// --------------------------------------------------------------------------
//
// Add a new upper edge to the edges stored in MBinning. The new upper
// edge must be greater than the current greatest. Using this you can
// enhance a histogram bin-by-bin, eg:
//   TH1F h("", "", 2, 0, 1);
//   MBinning b;
//   b.SetEdges(h);
//   b.AddEdge(2);
//   b.Apply(h);
//   b.AddEdge(3);
//   b.Apply(h);
//   [...]
//
void MBinning::AddEdge(Axis_t up)
{
    const Int_t n = fEdges.GetSize();

    if (up<=fEdges[n-1])
    {
        *fLog << warn << dbginf << "WARNING - New upper edge not greater than old upper edge... ignored." << endl;
        return;
    }

    fEdges.Set(n+1);
    fEdges[n] = up;

    fType = kIsUserArray;
}

void MBinning::RemoveFirstEdge()
{
    const Int_t n = fEdges.GetSize();
    for (int i=0; i<n-1; i++)
        fEdges[i] = fEdges[i+1];
    fEdges.Set(n-1);
}

void MBinning::RemoveLastEdge()
{
    fEdges.Set(fEdges.GetSize()-1);
}

// --------------------------------------------------------------------------
//
// Set the edges in MBinning from a histogram-axis. Which axis is
// specified by axis ('x', 'y', 'z')
//
void MBinning::SetEdges(const TH1 &h, const Char_t axis)
{
    switch (tolower(axis))
    {
    case 'x':
        SetEdges(*h.GetXaxis());
        return;
    case 'y':
        SetEdges(*h.GetYaxis());
        return;
    case 'z':
        SetEdges(*h.GetZaxis());
        return;
    default:
        *fLog << warn << "MBinning::SetEdges: Axis '" << axis << "' unknown... using x." << endl;
        SetEdges(*h.GetXaxis());
    }
}

// --------------------------------------------------------------------------
//
// Specify the number of bins <nbins> (not the number of edges), the
// lowest <lo> and highest <up> Edge (of your histogram)
//
void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up)
{
    const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
    fEdges.Set(nbins+1);
    for (int i=0; i<=nbins; i++)
        fEdges[i] = binsize*i + lo;

    fType = kIsLinear;
}

// --------------------------------------------------------------------------
//
// Calls SetEdgesLog if opt contains "log"
// Calls SetEdgesCos if opt contains "cos"
// Calls SetEdges in all other cases
//
void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up, const char *opt)
{
    const TString o(opt);
    if (o.Contains("log", TString::kIgnoreCase))
    {
        SetEdgesLog(nbins, lo, up);
        return;
    }
    if (o.Contains("cos", TString::kIgnoreCase))
    {
        SetEdgesCos(nbins, lo, up);
        return;
    }
    SetEdges(nbins, lo, up);
}

// --------------------------------------------------------------------------
//
// Specify the number of bins <nbins> (not the number of edges), the
// lowest <lo> and highest <up> Edge (of your histogram)
//
void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up)
{
    // if (lo==0) ...

    const Double_t binsize = log10(up/lo)/nbins;
    fEdges.Set(nbins+1);
    for (int i=0; i<=nbins; i++)
        fEdges[i] = pow(10, binsize*i) * lo;

    fType = kIsLogarithmic;
}

// --------------------------------------------------------------------------
//
// Specify the number of bins <nbins> (not the number of edges), the
// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
//
void MBinning::SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up)
{
    // if (lo==0) ...
    const Axis_t ld = lo/kRad2Deg;
    const Axis_t ud = up/kRad2Deg;

    const Double_t binsize = (cos(ld)-cos(ud))/nbins;
    fEdges.Set(nbins+1);
    for (int i=0; i<=nbins; i++)
        fEdges[i] = acos(cos(ld)-binsize*i)*kRad2Deg;

    fType = kIsCosinic;
}

// --------------------------------------------------------------------------
//
// Apply this binning to the given histogram.
// (By definition this works only for 1D-histograms. For 2D- and 3D-
//  histograms use MH::SetBinning directly)
//
void MBinning::Apply(TH1 &h) const
{
    if (h.InheritsFrom("TH2") || h.InheritsFrom("TH3"))
    {
        *fLog << warn << "MBinning::Apply: '" << h.GetName() << "' is not a basic TH1 object... no binning applied." << endl;
        return;
    }

    MH::SetBinning(&h, this);
}

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MBinning::StreamPrimitive(ofstream &out) const
{
    out << "   MBinning " << GetUniqueName();
    if (fName!=gsDefName || fTitle!=gsDefTitle)
    {
        out << "(\"" << fName << "\"";
        if (fTitle!=gsDefTitle)
            out << ", \"" << fTitle << "\"";
        out <<")";
    }
    out << ";" << endl;

    if (IsDefault())
        return;

    if (IsLinear() || IsLogarithmic() || IsCosinic())
    {
        out << "   " << GetUniqueName() << ".SetEdges";
        if (IsLogarithmic())
            out << "Log";
        if (IsCosinic())
            out << "Cos";
        out << "(" << GetNumBins() << ", " << GetEdgeLo() << ", " << GetEdgeHi() << ");" << endl;
        return;
    }

    out << "   {" << endl;
    out << "      TArrayD dummy;" << endl;
    for (int i=0; i<GetNumEdges(); i++)
        out << "      dummy[" << i << "]=" << GetEdges()[i] << ";" << endl;
    out << "      " << GetUniqueName() << ".SetEdges(dummy);" << endl;
    out << "   }" << endl;
}