source: trunk/MagicSoft/Mars/mhist/MHHillas.cc@ 1524

/* ======================================================================== *\
!
! *
! * 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  2001 <mailto:tbretz@uni-sw.gwdg.de>
!              Wolfgang Wittek  2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHHillas
//
// This class contains histograms for the source independent image parameters
//
/////////////////////////////////////////////////////////////////////////////

#include "MHHillas.h"

#include <math.h>

#include <TH1.h>
#include <TH2.h>
#include <TPad.h>
#include <TStyle.h>
#include <TCanvas.h>

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

#include "MParList.h"

#include "MHillas.h"
#include "MGeomCam.h"
#include "MBinning.h"


ClassImp(MHHillas);

// --------------------------------------------------------------------------
//
// Setup four histograms for Width, Length
//
MHHillas::MHHillas(const char *name, const char *title)
    : fMm2Deg(1), fUseMmScale(kTRUE)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHHillas";
    fTitle = title ? title : "Source independent image parameters";

    fLength  = new TH1F("Length",  "Length of Ellipse",               100,   0, 296.7);
    fWidth   = new TH1F("Width",   "Width of Ellipse",                100,   0, 296.7);
    fDistC   = new TH1F("DistC",   "Distance from center of camera",  100,   0, 445);
    fDelta   = new TH1F("Delta",   "Angle (Main axis - x-axis)",      101, -90,  90);
    fUsedPix = new TH1F("UsedPix", "Number of used pixels",           150,   0, 150);
    fCorePix = new TH1F("CorePix", "Number of core pixels",           150,   0, 150);

    fLength->SetLineColor(kBlue);
    fCorePix->SetLineColor(kRed);
    fUsedPix->SetLineColor(kGreen);

    fLength->SetDirectory(NULL);
    fWidth->SetDirectory(NULL);
    fDistC->SetDirectory(NULL);
    fDelta->SetDirectory(NULL);
    fUsedPix->SetDirectory(NULL);
    fCorePix->SetDirectory(NULL);

    fLength->SetXTitle("Length [mm]");
    fWidth->SetXTitle("Width [mm]");
    fDistC->SetXTitle("Distance [mm]");
    fDelta->SetXTitle("Delta [\\circ]");
    fUsedPix->SetXTitle("Number of Pixels");
    fCorePix->SetXTitle("Number of Pixels");

    fLength->SetYTitle("Counts");
    fWidth->SetYTitle("Counts");
    fDistC->SetYTitle("Counts");
    fDelta->SetYTitle("Counts");
    fUsedPix->SetYTitle("Counts");
    fCorePix->SetYTitle("Counts");

    MBinning bins;
    bins.SetEdgesLog(50, 1, 1e7);

    fSize  = new TH1F;
    fSize->SetName("Size");
    fSize->SetTitle("Number of Photons");
    fSize->SetDirectory(NULL);
    fSize->SetXTitle("Size");
    fSize->SetYTitle("Counts");
    fSize->GetXaxis()->SetTitleOffset(1.2);
    fSize->GetXaxis()->SetLabelOffset(-0.015);

    bins.Apply(*fSize);

    fCenter = new TH2F("Center", "Center of Ellipse", 51, -445, 445, 51, -445, 445);
    fCenter->SetDirectory(NULL);
    fCenter->SetXTitle("x [mm]");
    fCenter->SetYTitle("y [mm]");
    fCenter->SetZTitle("Counts");
}

// --------------------------------------------------------------------------
//
// Delete the histograms
//
MHHillas::~MHHillas()
{
    delete fLength;
    delete fWidth;

    delete fDistC;
    delete fDelta;

    delete fSize;
    delete fCenter;

    delete fUsedPix;
    delete fCorePix;
}

// --------------------------------------------------------------------------
//
// Setup the Binning for the histograms automatically if the correct
// instances of MBinning (with the names 'BinningWidth' and 'BinningLength')
// are found in the parameter list
// Use this function if you want to set the conversion factor which
// is used to convert the mm-scale in the camera plain into the deg-scale
// used for histogram presentations. The conversion factor is part of
// the camera geometry. Please create a corresponding MGeomCam container.
//
Bool_t MHHillas::SetupFill(const MParList *plist)
{
    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!geom)
        *fLog << warn << GetDescriptor() << ": No Camera Geometry available. Using mm-scale for histograms." << endl;
    else
    {
        fMm2Deg = geom->GetConvMm2Deg();
        SetMmScale(kFALSE);
    }

    ApplyBinning(*plist, "Width",  fWidth);
    ApplyBinning(*plist, "Length", fLength);
    ApplyBinning(*plist, "Dist",   fDistC);
    ApplyBinning(*plist, "Delta",  fDelta);
    ApplyBinning(*plist, "Size",   fSize);
    ApplyBinning(*plist, "Pixels", fUsedPix);
    ApplyBinning(*plist, "Pixels", fCorePix);

    const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
    if (!bins)
    {
        float r = geom ? geom->GetMaxRadius() : 600;
        r *= 0.9;
        if (!fUseMmScale)
            r *= fMm2Deg;

        MBinning b;
        b.SetEdges(61, -r, r);
        SetBinning(fCenter, &b, &b);
    }
    else
        SetBinning(fCenter, bins, bins);


    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Use this function to setup your own conversion factor between degrees
// and millimeters. The conversion factor should be the one calculated in
// MGeomCam. Use this function with Caution: You could create wrong values
// by setting up your own scale factor.
//
void MHHillas::SetMm2Deg(Float_t mmdeg)
{
    if (mmdeg<0)
    {
        *fLog << warn << dbginf << "Warning - Conversion factor < 0 - nonsense. Ignored." << endl;
        return;
    }

    if (fMm2Deg>=0)
        *fLog << warn << dbginf << "Warning - Conversion factor already set. Overwriting" << endl;

    fMm2Deg = mmdeg;
}

// --------------------------------------------------------------------------
//
// With this function you can convert the histogram ('on the fly') between
// degrees and millimeters.
//
void MHHillas::SetMmScale(Bool_t mmscale)
{
    if (fUseMmScale == mmscale)
        return;

    if (fMm2Deg<0)
    {
        *fLog << warn << dbginf << "Warning - Sorry, no conversion factor for conversion available." << endl;
        return;
    }

    const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg;
    MH::ScaleAxis(fLength, scale);
    MH::ScaleAxis(fWidth,  scale);
    MH::ScaleAxis(fDistC,  scale);
    MH::ScaleAxis(fCenter, scale, scale);

    if (mmscale)
    {
        fLength->SetXTitle("Length [mm]");
        fWidth->SetXTitle("Width [mm]");
        fDistC->SetXTitle("Distance [mm]");
        fCenter->SetXTitle("x [mm]");
        fCenter->SetYTitle("y [mm]");
    }
    else
    {
        fLength->SetXTitle("Length [\\circ]");
        fWidth->SetXTitle("Width [\\circ]");
        fDistC->SetXTitle("Distance [\\circ]");
        fCenter->SetXTitle("x [\\circ]");
        fCenter->SetYTitle("y [\\circ]");
    }

    fUseMmScale = mmscale;
}

// --------------------------------------------------------------------------
//
// Fill the histograms with data from a MHillas-Container.
// Be careful: Only call this with an object of type MHillas
//
Bool_t MHHillas::Fill(const MParContainer *par)
{
    const MHillas &h = *(MHillas*)par;

    const Double_t d = sqrt(h.GetMeanX()*h.GetMeanX() + h.GetMeanY()*h.GetMeanY());
    const Double_t scale = fUseMmScale ? 1 : fMm2Deg;

    fLength ->Fill(scale*h.GetLength());
    fWidth  ->Fill(scale*h.GetWidth());
    fDistC  ->Fill(scale*d);
    fCenter ->Fill(scale*h.GetMeanX(), scale*h.GetMeanY());
    fDelta  ->Fill(kRad2Deg*h.GetDelta());
    fSize   ->Fill(h.GetSize());
    fUsedPix->Fill(h.GetNumUsedPixels());
    fCorePix->Fill(h.GetNumCorePixels());

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Setup a inversed deep blue sea palette for the fCenter histogram.
//
void MHHillas::SetColors() const
{
    // FIXME: This must be redone each time the canvas is repainted....
    gStyle->SetPalette(51, NULL);
    Int_t c[50];
    for (int i=0; i<50; i++)
        c[49-i] = gStyle->GetColorPalette(i);
    gStyle->SetPalette(50, c);
}

// --------------------------------------------------------------------------
//
// Draw clones of four histograms. So that the object can be deleted
// and the histograms are still visible in the canvas.
// The cloned object are deleted together with the canvas if the canvas is
// destroyed. If you want to handle dostroying the canvas you can get a
// pointer to it from this function
//
TObject *MHHillas::DrawClone(Option_t *opt) const
{
    TCanvas *c = MakeDefCanvas(this, 720, 810);
    c->Divide(2,3);

    gROOT->SetSelectedPad(NULL);

    c->cd(1);
    DrawCopy(*fWidth, *fLength, "Width / Length");

    c->cd(2);
    gPad->SetLogx();
    fSize->DrawCopy();

    c->cd(3);
    DrawCopy(*fCorePix, *fUsedPix, "Number of core/used Pixels");

    c->cd(4);
    fDelta->DrawCopy();

    c->cd(5);
    fDistC->DrawCopy();

    c->cd(6);
    SetColors();
    fCenter->DrawCopy("colz");

    c->Modified();
    c->Update();

    return c;
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the four histograms into it.
// Be careful: The histograms belongs to this object and won't get deleted
// together with the canvas.
//
void MHHillas::Draw(Option_t *)
{
    if (!gPad)
        MakeDefCanvas(this, 720, 810);

    gPad->Divide(2,3);

    gPad->cd(1);
    MH::Draw(*fWidth, *fLength, "Width / Length");

    gPad->cd(2);
    gPad->SetLogx();
    fSize->Draw();

    gPad->cd(3);
    MH::Draw(*fCorePix, *fUsedPix, "Number of core/used Pixels");

    gPad->cd(4);
    fDelta->Draw();

    gPad->cd(5);
    fDistC->Draw();

    gPad->cd(6);
    SetColors();
    fCenter->Draw("colz");

    gPad->Modified();
    gPad->Update();
}