source: trunk/MagicSoft/Mars/mhist/MHStarMap.cc@ 1889

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

/////////////////////////////////////////////////////////////////////////////
//
// MHStarMap
//
// This class contains a 2-dimensional histogram. It should show some
// kind of star map. The algorith which calculates the star map
// from the Hillas parameters (Fill) can be enhanced.
//
/////////////////////////////////////////////////////////////////////////////

#include "MHStarMap.h"

#include <TStyle.h>   // gStyle
#include <TColor.h>   // SetRGB
#include <TCanvas.h>  // TCanvas

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

#include "MParList.h"

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

ClassImp(MHStarMap);

// --------------------------------------------------------------------------
//
// Create the star map histogram
//
MHStarMap::MHStarMap(const char *name, const char *title)
    : fMm2Deg(-1), fUseMmScale(kTRUE)
{
    //
    //  default constructor
    //  creates an a list of histograms for all pixels and both gain channels
    //

    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHStarMap" ;
    fTitle = title ? title : "Container for a Star Map" ;

    *fLog << warn << "WARNING - Using MHStarMap doesn't take care of the Source Position!" << endl;

    //
    //   loop over all Pixels and create two histograms
    //   one for the Low and one for the High gain
    //   connect all the histogram with the container fHist
    //
    fStarMap = new TH2F("StarMap", " 2D Hillas Star Map ",
                        150, -300, 300,
                        150, -300, 300);

    fStarMap->SetDirectory(NULL);

    fStarMap->SetXTitle("x [mm]");
    fStarMap->SetYTitle("y [mm]");
    fStarMap->SetZTitle("Counts");
}

// --------------------------------------------------------------------------
//
// delete the histogram instance
//
MHStarMap::~MHStarMap()
{
    delete fStarMap;
}

// --------------------------------------------------------------------------
//
// Setup the Binning for the histograms automatically if the correct
// instances of MBinning (with the names 'BinningCamera')
// 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 MHStarMap::SetupFill(const MParList *plist)
{
    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    if (geom)
    {
        fMm2Deg = geom->GetConvMm2Deg();
        fUseMmScale = kFALSE;

        fStarMap->SetXTitle("x [\\circ]");
        fStarMap->SetYTitle("y [\\circ]");
    }

    const MBinning *bins = (MBinning*)plist->FindObject("BinningStarMap");
    if (!bins)
    {
        float r = geom ? geom->GetMaxRadius() : 600;
        r *= 5./6;
        if (!fUseMmScale)
            r *= fMm2Deg;

        MBinning b;
        b.SetEdges(100, -r, r);
        SetBinning(fStarMap, &b, &b);
    }
    else
        SetBinning(fStarMap, bins, bins);

    if (!geom)
    {
        *fLog << warn << dbginf << "No Camera Geometry available. Using mm-scale for histograms." << endl;
        return kTRUE;
    }

    return kTRUE;
}

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

    const float delta = h.GetDelta();

    const float m    = tan(delta);
    const float cosd = 1.0/sqrt(1.0+m*m);
    const float sind = sqrt(1.0-cosd*cosd);


    float t = h.GetMeanY() - m*h.GetMeanX();

    if (!fUseMmScale)
        t *= fMm2Deg;

    // get step size ds along the main axis of the ellipse
    TAxis &axe = *fStarMap->GetXaxis();
    const int   N    = axe.GetNbins();
    const float xmin = axe.GetBinLowEdge(1);
    const float xmax = axe.GetBinLowEdge(N+1);
    const float ds   = (xmax-xmin) / 200.0;

    if (m>-1 && m<1)
    {
        float dx = ds * cosd;
        float  x = xmin + dx/2.0;
        int   N1 = (int) ((xmax-xmin)/dx+1.0);

        for (int i=0; i<N1; i++)
        {
            const float y = m*x+t;
            fStarMap->Fill(x, y);
            x += dx;
        }
    }
    else
    {
        TAxis &axe = *fStarMap->GetYaxis();
        const int   M    = axe.GetNbins();
        const float ymin = axe.GetBinLowEdge(1);
        const float ymax = axe.GetBinLowEdge(M+1);

        float dy = ds * sind;
        float  y = ymin + dy/2.0;
        int   M1 = (int) ((ymax-ymin)/dy+1.0);

        for (int i=0; i<M1; i++)
        {
            const float x = (y-t)/m;
            fStarMap->Fill(x, y);
            y += dy;
        }
    }

    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 MHStarMap::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 MHStarMap::SetMmScale(Bool_t mmscale)
{
    if (fUseMmScale == mmscale)
        return;

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

    if (fUseMmScale)
    {
        fStarMap->SetXTitle("x [mm]");
        fStarMap->SetYTitle("y [mm]");

        fStarMap->Scale(1./fMm2Deg);
    }
    else
    {
        fStarMap->SetXTitle("x [\\circ]");
        fStarMap->SetYTitle("y [\\circ]");

        fStarMap->Scale(1./fMm2Deg);
    }

    fUseMmScale = mmscale;
}

// --------------------------------------------------------------------------
//
// Set the palette you wanna use:
//  - you could set the root "Pretty Palette Violet->Red" by
//    gStyle->SetPalette(1, 0), but in some cases this may look
//    confusing
//  - The maximum colors root allowes us to set by ourself
//    is 50 (idx: 51-100). This colors are set to a grayscaled
//    palette
//  - the number of contours must be two less than the number
//    of palette entries
//
void MHStarMap::PrepareDrawing() const
{
    const Int_t numg = 32; // number of gray scaled colors
    const Int_t numw = 32; // number of white

    Int_t palette[numg+numw];

    //
    // The first half of the colors are white.
    // This is some kind of optical background supression
    //
    gROOT->GetColor(51)->SetRGB(1, 1, 1);

    Int_t i;
    for (i=0; i<numw; i++)
        palette[i] = 51;

    //
    // now the (gray) scaled part is coming
    //
    for (;i<numw+numg; i++)
    {
        const Float_t gray = 1.0-(float)(i-numw)/(numg-1.0);

        gROOT->GetColor(52+i)->SetRGB(gray, gray, gray);
        palette[i] = 52+i;
    }

    //
    // Set the palette and the number of contour levels
    //
    gStyle->SetPalette(numg+numw, palette);
    fStarMap->SetContour(numg+numw-2);
}


// --------------------------------------------------------------------------
//
// Draw clones of the  histograms, so that the object can be deleted
// and the histogram is still visible in the canvas.
// The cloned object is deleted together with the canvas if the canvas is
// destroyed. If you want to handle destroying the canvas you can get a
// pointer to it from this function
//
TObject *MHStarMap::DrawClone(Option_t *opt) const
{
    TCanvas *c=MakeDefCanvas(fStarMap, 500, 500);

    //
    // This is necessary to get the expected bahviour of DrawClone
    //
    gROOT->SetSelectedPad(NULL);

    PrepareDrawing();

    fStarMap->DrawCopy("colz");

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

    return c;
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draw the histogram into it.
// Be careful: The histogram belongs to this object and won't get deleted
// together with the canvas.
//
void MHStarMap::Draw(Option_t *)
{
    if (!gPad)
        MakeDefCanvas(fStarMap, 500, 500);

    PrepareDrawing();

    fStarMap->Draw("colz");

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