source: trunk/MagicSoft/Mars/mhist/MHCamera.cc@ 2236

/* ======================================================================== *\
!
! *
! * 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, 05/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer, 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHCamera
//
// Camera Display, based on a TH1D. Pleas be carefull using the
// underlaying TH1D.
//
// To change the scale to a logarithmic scale SetLogy() of the Pad.
//
////////////////////////////////////////////////////////////////////////////
#include "MHCamera.h"

#include <fstream>
#include <iostream>

#include <TBox.h>
#include <TArrow.h>
#include <TLatex.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TArrayF.h>
#include <TRandom.h>
#include <TPaveText.h>
#include <TClonesArray.h>

#include "MH.h"
#include "MHexagon.h"

#include "MGeomPix.h"
#include "MGeomCam.h"

#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

#include "MPedestalPix.h"
#include "MPedestalCam.h"

#include "MCurrents.h"
#include "MCamEvent.h"

#include "MImgCleanStd.h"

#define kItemsLegend 48 // see SetPalette(1,0)

ClassImp(MHCamera);

using namespace std;

// ------------------------------------------------------------------------
//
//  Default Constructor. To be used by the root system ONLY.
//
MHCamera::MHCamera() : TH1D(), fGeomCam(NULL), fColors(kItemsLegend)
{
    SetDirectory(NULL);

    fNotify  = NULL;

#if ROOT_VERSION_CODE < ROOT_VERSION(3,01,06)
    SetPalette(1, 0);
#else
    SetPalette(51, 0);
#endif
}

// ------------------------------------------------------------------------
//
//  Constructor. Makes a clone of MGeomCam. Removed the TH1D from the
// current directory. Calls Sumw2(). Set the histogram line color
// (for error bars) to Green and the marker style to kFullDotMedium.
//
MHCamera::MHCamera(const MGeomCam &geom, const char *name, const char *title)
: TH1D(name, title, geom.GetNumPixels(), -0.5, geom.GetNumPixels()-0.5),
fUsed(geom.GetNumPixels()), fColors(kItemsLegend)
{
    fGeomCam = (MGeomCam*)geom.Clone();

    SetDirectory(NULL);
    Sumw2();

    SetLineColor(kGreen);
    SetMarkerStyle(kFullDotMedium);
    SetXTitle("Pixel Index");

    fNotify = new TList;

    //
    //  create the hexagons of the display
    //
    // root 3.02
    //  * base/inc/TObject.h:
    //    register BIT(8) as kNoContextMenu. If an object has this bit set it will
    //    not get an automatic context menu when clicked with the right mouse button.

    //
    // Construct all hexagons. Use new-operator with placement
    //
    for (Int_t i=0; i<fNcells-2; i++)
        ResetUsed(i);

#if ROOT_VERSION_CODE < ROOT_VERSION(3,01,06)
    SetPalette(1, 0);
#else
    SetPalette(51, 0);
#endif
}

// ------------------------------------------------------------------------
//
// Destructor. Deletes the cloned fGeomCam and the notification list.
//
MHCamera::~MHCamera()
{
    if (fGeomCam)
        delete fGeomCam;
    if (fNotify)
        delete fNotify;
}

// ------------------------------------------------------------------------
//
// Taken from TH1D::Fill(). Uses the argument directly as bin index.
// Doesn't increment the number of entries.
//
//   -*-*-*-*-*-*-*-*Increment bin with abscissa X by 1*-*-*-*-*-*-*-*-*-*-*
//                   ==================================
//
//    if x is less than the low-edge of the first bin, the Underflow bin is incremented
//    if x is greater than the upper edge of last bin, the Overflow bin is incremented
//
//    If the storage of the sum of squares of weights has been triggered,
//    via the function Sumw2, then the sum of the squares of weights is incremented
//    by 1 in the bin corresponding to x.
//
//   -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t MHCamera::Fill(Axis_t x)
{

#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
   if (fBuffer) return BufferFill(x,1);
#endif
   const Int_t bin = (Int_t)x+1;
   AddBinContent(bin);
   if (fSumw2.fN)
       fSumw2.fArray[bin]++;

   if (bin<=0 || bin>fNcells-2)
       return -1;

   fTsumw++;
   fTsumw2++;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   return bin;
}

// ------------------------------------------------------------------------
//
// Taken from TH1D::Fill(). Uses the argument directly as bin index.
// Doesn't increment the number of entries.
//
//   -*-*-*-*-*-*Increment bin with abscissa X with a weight w*-*-*-*-*-*-*-*
//               =============================================
//
//    if x is less than the low-edge of the first bin, the Underflow bin is incremented
//    if x is greater than the upper edge of last bin, the Overflow bin is incremented
//
//    If the storage of the sum of squares of weights has been triggered,
//    via the function Sumw2, then the sum of the squares of weights is incremented
//    by w^2 in the bin corresponding to x.
//
//   -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Int_t MHCamera::Fill(Axis_t x, Stat_t w)
{
#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
   if (fBuffer) return BufferFill(x,w);
#endif
   const Int_t bin = (Int_t)x+1;
   AddBinContent(bin, w);
   if (fSumw2.fN)
       fSumw2.fArray[bin] += w*w;

   if (bin<=0 || bin>fNcells-2)
       return -1;

   const Stat_t z = (w > 0 ? w : -w);
   fTsumw   += z;
   fTsumw2  += z*z;
   fTsumwx  += z*x;
   fTsumwx2 += z*x*x;
   return bin;
}

// ------------------------------------------------------------------------
//
// Return the minimum contents of all pixels (if all is set, otherwise
// only of all 'used' pixels), fMinimum if fMinimum set
//
Double_t MHCamera::GetMinimum(Bool_t all) const
{
    if (fMinimum != -1111)
        return fMinimum;

    Double_t minimum=FLT_MAX;

    if (all)
    {
        for (Int_t idx=0; idx<fNcells-2; idx++)
            if (fArray[idx+1] < minimum)
                minimum = fArray[idx+1];
    }
    else
    {
        for (Int_t idx=0; idx<fNcells-2; idx++)
            if (IsUsed(idx) && fArray[idx+1] < minimum)
                minimum = fArray[idx+1];
    }
    return minimum;
}

// ------------------------------------------------------------------------
//
// Return the maximum contents of all pixels (if all is set, otherwise
// only of all 'used' pixels), fMaximum if fMaximum set
//
Double_t MHCamera::GetMaximum(Bool_t all) const
{
    if (fMaximum != -1111)
        return fMaximum;

    Double_t maximum=-FLT_MAX;
    if (all)
    {
        for (Int_t idx=0; idx<fNcells-2; idx++)
            if (fArray[idx+1] > maximum)
                maximum = fArray[idx+1];
    }
    else
    {
        for (Int_t idx=0; idx<fNcells-2; idx++)
            if (IsUsed(idx) && fArray[idx+1] > maximum)
                maximum = fArray[idx+1];
    }
    return maximum;
}

// ------------------------------------------------------------------------
//
// Call this function to draw the camera layout into your canvas.
// Setup a drawing canvas. Add this object and all child objects
// (hexagons, etc) to the current pad. If no pad exists a new one is
// created.
//
// To draw a camera into its own pad do something like:
//
// TCanvas *c = new TCanvas;
// c->Divide(2,1);
// MGeomCamMagic m;
// MHCamera *d=new MHCamera(&m);
// d->FillRandom();
// c->cd(1);
// gPad->SetBorderMode(0);
// gPad->Divide(1,1);
// gPad->cd(1);
// d->Draw();
// d->SetBit(kCanDelete);
//
void MHCamera::Draw(Option_t *option)
{
    // root 3.02:
    // gPad->SetFixedAspectRatio()
    Int_t col = 16;

    if (gPad)
        col = gPad->GetFillColor();

    TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas("CamDisplay", "Mars Camera Display", 656, 600);
    pad->SetBorderMode(0);
    pad->SetFillColor(col);

    AppendPad(option);
}

// ------------------------------------------------------------------------
//
// Resizes the current pad so that the camera is displayed in its
// correct aspect ratio
//
void MHCamera::SetRange()
{
    const Float_t range = fGeomCam->GetMaxRadius();

    //
    // Maintain aspect ratio
    //
    const float ratio = 1.15;

    //
    // Calculate width and height of the current pad in pixels
    //
    Float_t w = gPad->GetWw();
    Float_t h = gPad->GetWh()*ratio;

    //
    // This prevents the pad from resizing itself wrongly
    //
    if (gPad->GetMother() != gPad)
    {
        w *= gPad->GetMother()->GetAbsWNDC();
        h *= gPad->GetMother()->GetAbsHNDC();
    }

    //
    // Set Range (coordinate system) of pad
    //
    gPad->Range(-range, -range, (2*ratio-1)*range, range);

    //
    // Resize Pad to given ratio
    //
    if (h<w)
        gPad->SetPad((1.-h/w)/2, 0, (h/w+1.)/2, 1);
    else
        gPad->SetPad(0, (1.-w/h)/2, 1, (w/h+1.)/2);
}

// ------------------------------------------------------------------------
//
// Updates the pixel colors and paints the pixels
//
void MHCamera::Update(Bool_t islog, Bool_t isbox, Bool_t iscol)
{
    Double_t min = GetMinimum(kFALSE);
    Double_t max = GetMaximum(kFALSE);
    if (min==FLT_MAX)
    {
        min = 0;
        max = 1;
    }

    if (min==max)
        max += 1;

    UpdateLegend(min, max, islog);

    MHexagon hex;
    for (Int_t i=0; i<fNcells-2; i++)
    {
        if (IsUsed(i) && iscol)
            hex.SetFillColor(GetColor(fArray[i+1], min, max, islog));
        else
            hex.SetFillColor(10);

        MGeomPix &pix = (*fGeomCam)[i];
        if (!isbox)
            hex.PaintHexagon(pix.GetX(), pix.GetY(), pix.GetD());
        else
            if (IsUsed(i))
            {
                Float_t size = pix.GetD()*(fArray[i+1]-min)/(max-min);
                if (size>pix.GetD())
                    size=pix.GetD();
                hex.PaintHexagon(pix.GetX(), pix.GetY(), size);
            }
    }
}

// ------------------------------------------------------------------------
//
// Print minimum and maximum
//
void MHCamera::Print(Option_t *) const
{
    cout << "Minimum: " << GetMinimum();
    if (fMinimum==-1111)
        cout << " <autoscaled>";
    cout << endl;
    cout << "Maximum: " << GetMaximum();
    if (fMaximum==-1111)
        cout << " <autoscaled>";
    cout << endl;
}

// ------------------------------------------------------------------------
//
// Paint the y-axis title
//
void MHCamera::PaintAxisTitle()
{
    Float_t fRange = fGeomCam->GetMaxRadius();

    TLatex *ptitle = new TLatex(1.2*fRange, .97*fRange, GetYaxis()->GetTitle());

    ptitle->SetTextSize(0.03);
    ptitle->SetTextAlign(33);

    // box with the histogram title
    ptitle->SetTextColor(gStyle->GetTitleTextColor());
#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
    ptitle->SetTextFont(gStyle->GetTitleFont(""));
#endif
    ptitle->Paint();
}

// ------------------------------------------------------------------------
//
// Paint the histogram title
//
void MHCamera::PaintTitle()
{
//    *-*-*-*-*-*-*-*-*-*Draw the histogram title*-*-*-*-*-*-*-*-*-*-*-*-*
//                       ========================
   //if (Hoption.Same) return;
#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
    if (TestBit(kNoTitle))
        return;
#endif

    const Int_t nt = strlen(GetTitle());

    TPaveText *title = (TPaveText*)gPad->FindObject("title");
    if (nt == 0 || gStyle->GetOptTitle() <= 0)
    {
        if (title)
            delete title;
        return;
    }

    Double_t ht = gStyle->GetTitleH();
    Double_t wt = gStyle->GetTitleW();

    if (ht <= 0)
        ht = 0.05;
    if (wt <= 0)
    {
        TLatex l;
        l.SetTextSize(ht);
        l.SetTitle(GetTitle());
        Double_t wndc = l.GetXsize()/(gPad->GetX2() - gPad->GetX1());
        wt = TMath::Min(0.7, 0.02+wndc);
    }
    if (title)
    {
        TText *t0 = (TText*)title->GetLine(0);
        if (t0)
        {
            if (!strcmp(t0->GetTitle(), GetTitle()))
                return;

            t0->SetTitle(GetTitle());
            if (wt > 0)
                title->SetX2NDC(title->GetX1NDC()+wt);
        }
        return;
    }

    TPaveText *ptitle = new TPaveText(
                                      gStyle->GetTitleX(),
                                      gStyle->GetTitleY()-ht,
                                      gStyle->GetTitleX()+wt,
                                      gStyle->GetTitleY(),"blNDC");

    // box with the histogram title
#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
    ptitle->SetFillColor(gStyle->GetTitleFillColor());
    ptitle->SetTextFont(gStyle->GetTitleFont(""));
    if (gStyle->GetTitleFont("")%10 > 2)
        ptitle->SetTextSize(gStyle->GetTitleFontSize());
#endif
    ptitle->SetFillStyle(gStyle->GetTitleStyle());
    ptitle->SetName("title");
    ptitle->SetBorderSize(gStyle->GetTitleBorderSize());
    ptitle->SetTextColor(gStyle->GetTitleTextColor());
    ptitle->AddText(GetTitle());
    ptitle->SetBit(kCanDelete);
    ptitle->Draw();
    ptitle->Paint();
}

// ------------------------------------------------------------------------
//
// Paints the camera.
//
void MHCamera::Paint(Option_t *o)
{
    TString opt(o);
    opt.ToLower();

    if (opt.Contains("hist"))
    {
        opt.ReplaceAll("hist", "");

        Int_t mode = gStyle->GetOptStat();
        TVirtualPad *save = gPad;
        gPad=NULL;
        gStyle->SetOptStat(1000011);
        gPad=save;
        TH1D::Paint(o);
        gPad=NULL;
        gStyle->SetOptStat(mode);
        gPad=save;
        return;
    }

    // Maintain aspect ratio
    SetRange();

    Bool_t isbox = opt.Contains("box");
    Bool_t iscol = isbox ? !opt.Contains("nocol") : 1;

    // Update Contents of the pixels and paint legend
    Update(gPad->GetLogy(), isbox, iscol);

    // Paint primitives (pixels, color legend, photons, ...)
    PaintTitle();
    PaintAxisTitle();
}

// ------------------------------------------------------------------------
//
// With this function you can change the color palette. For more
// information see TStyle::SetPalette. Only palettes with 50 colors
// are allowed.
// In addition you can use SetPalette(52, 0) to create an inverse
// deep blue sea palette
//
void MHCamera::SetPalette(Int_t ncolors, Int_t *colors)
{
    //
    // If not enough colors are specified skip this.
    //
    if (ncolors>1 && ncolors<50)
    {
        cout << "MHCamera::SetPalette: Only default palettes with 50 colors are allowed... ignored." << endl;
        return;
    }

    //
    // If ncolors==52 create a reversed deep blue sea palette
    //
    if (ncolors==52)
    {
        gStyle->SetPalette(51, NULL);
        TArrayI c(kItemsLegend);
        for (int i=0; i<kItemsLegend; i++)
            c[kItemsLegend-i-1] = gStyle->GetColorPalette(i);
        gStyle->SetPalette(kItemsLegend, c.GetArray());
    }
    else
        gStyle->SetPalette(ncolors, colors);

    fColors.Set(kItemsLegend);
    for (int i=0; i<kItemsLegend; i++)
        fColors[i] = gStyle->GetColorPalette(i);
}


void MHCamera::SetPrettyPalette()
{
    if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
        SetPalette(1, 0);
}

void MHCamera::SetDeepBlueSeaPalette()
{
    if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
        SetPalette(51, 0);
}

void MHCamera::SetInvDeepBlueSeaPalette()
{
    if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
        SetPalette(52, 0);
}

void MHCamera::DrawPixelIndices()
{
    // FIXME: Is this correct?
    for (int i=0; i<kItemsLegend; i++)
        fColors[i] = 16;

    if (!gPad)
        Draw();

    TText txt;
    txt.SetTextFont(122);
    txt.SetTextAlign(22);   // centered/centered

    for (Int_t i=0; i<fNcells-2; i++)
    {
        TString num;
        num += i;

        const MGeomPix &h = (*fGeomCam)[i];
        TText *nt = txt.DrawText(h.GetX(), h.GetY(), num);
        nt->SetTextSize(0.3*h.GetD()/fGeomCam->GetMaxRadius());
    }
}

void MHCamera::DrawSectorIndices()
{
    for (int i=0; i<kItemsLegend; i++)
        fColors[i] = 16;

    if (!gPad)
        Draw();

    TText txt;
    txt.SetTextFont(122);
    txt.SetTextAlign(22);   // centered/centered

    for (Int_t i=0; i<fNcells-2; i++)
    {
        TString num;
        num += (*fGeomCam)[i].GetSector();

        const MGeomPix &h = (*fGeomCam)[i];
        TText *nt = txt.DrawText(h.GetX(), h.GetY(), num);
        nt->SetTextSize(0.3*h.GetD()/fGeomCam->GetMaxRadius());
    }
}

// ------------------------------------------------------------------------
//
// Call this function to add a MCamEvent on top of the present contents.
// Only 'used' pixels are added.
//
void MHCamera::AddCamContent(const MCamEvent &event, Int_t type)
{
    // FIXME: Security check missing!
    for (Int_t idx=0; idx<fNcells-2; idx++)
    {
        Double_t val=0;
        if (event.GetPixelContent(val, idx, *fGeomCam, type) && !IsUsed(idx))
            SetUsed(idx);

        Fill(idx, val); // FIXME: Slow!
    }
    fEntries++;
}

// ------------------------------------------------------------------------
//
// Call this function to add a MHCamera on top of the present contents.
// Only 'used' pixels are added.
// Type:
//  0) bin content
//  1) errors
//  2) rel. errors
//
void MHCamera::AddCamContent(const MHCamera &d, Int_t type)
{
    if (fNcells!=d.fNcells)
        return;

    // FIXME: Security check missing!
    for (Int_t idx=0; idx<fNcells-2; idx++)
        if (d.IsUsed(idx))
            SetUsed(idx);

    switch (type)
    {
    case 1:
        for (Int_t idx=0; idx<fNcells-2; idx++)
            Fill(idx, d.GetBinError(idx+1));
        break;
    case 2:
        for (Int_t idx=0; idx<fNcells-2; idx++)
            if (d.GetBinContent(idx+1)!=0)
                Fill(idx, fabs(d.GetBinError(idx+1)/d.GetBinContent(idx+1)));
        break;
    default:
        for (Int_t idx=0; idx<fNcells-2; idx++)
            Fill(idx, d.GetBinContent(idx+1));
        break;
    }
    fEntries++;
}

// ------------------------------------------------------------------------
//
// Call this function to add a TArrayD on top of the present contents.
// Only 'used' pixels are added.
//
void MHCamera::AddCamContent(const TArrayD &event, const TArrayC *used)
{
    if (event.GetSize()!=fNcells-2)
        return;

    if (used && used->GetSize()!=fNcells-2)
        return;

    for (Int_t idx=0; idx<fNcells-2; idx++)
    {
        Fill(idx, const_cast<TArrayD&>(event)[idx]); // FIXME: Slow!

        if (used && (*used)[idx])
            SetUsed(idx);
    }
    fEntries++;
}

// ------------------------------------------------------------------------
//
// Call this function to add a MCamEvent on top of the present contents.
// 1 is added to each pixel if the contents of MCamEvent>threshold
//
void MHCamera::CntCamContent(const MCamEvent &event, Double_t threshold, Int_t type)
{
    // FIXME: Security check missing!
    for (Int_t idx=0; idx<fNcells-2; idx++)
    {
        Double_t val=0;
        if (event.GetPixelContent(val, idx, *fGeomCam, type) && !IsUsed(idx))
            SetUsed(idx);

        if (val>threshold)
            Fill(idx);
    }
    fEntries++;
}

// ------------------------------------------------------------------------
//
// Call this function to add a TArrayD on top of the present contents.
// 1 is added to each pixel if the contents of MCamEvent>threshold
//
void MHCamera::CntCamContent(const TArrayD &event, Double_t threshold, Bool_t ispos)
{
    if (event.GetSize()!=fNcells-2)
        return;

    for (Int_t idx=0; idx<fNcells-2; idx++)
    {
        if (const_cast<TArrayD&>(event)[idx]>threshold)
            Fill(idx);

        if (!ispos || fArray[idx+1]>0)
            SetUsed(idx);
    }
    fEntries++;
}

// ------------------------------------------------------------------------
//
// Fill the pixels with random contents.
//
void MHCamera::FillRandom()
{
    Reset();

    // FIXME: Security check missing!
    for (Int_t idx=0; idx<fNcells-2; idx++)
    {
        Fill(idx, gRandom->Uniform()*fGeomCam->GetPixRatio(idx));
        SetUsed(idx);
    }
    fEntries=1;
}


// ------------------------------------------------------------------------
//
// Fill the colors in respect to the cleaning levels
//
void MHCamera::FillLevels(const MCerPhotEvt &event, Float_t lvl1, Float_t lvl2)
{
    SetCamContent(event, 2);

    for (Int_t i=0; i<fNcells-2; i++)
    {
        if (!IsUsed(i))
            continue;

        if (fArray[i+1]>lvl1)
            fArray[i+1] = 0;
        else
            if (fArray[i+1]>lvl2)
                fArray[i+1] = 1;
            else
                fArray[i+1] = 2;
    }
}

// ------------------------------------------------------------------------
//
// Fill the colors in respect to the cleaning levels
//
void MHCamera::FillLevels(const MCerPhotEvt &event, const MImgCleanStd &clean)
{
    FillLevels(event, clean.GetCleanLvl1(), clean.GetCleanLvl2());
}

// ------------------------------------------------------------------------
//
// Reset the all pixel colors to a default value
//
void MHCamera::Reset(Option_t *opt)
{
    TH1::Reset(opt);

    for (Int_t i=0; i<fNcells-2; i++)
    {
        fArray[i+1]=0;
        ResetUsed(i);
    }
    fArray[0]         = 0;
    fArray[fNcells-1] = 0;
} 

// ------------------------------------------------------------------------
//
//  Here we calculate the color index for the current value.
//  The color index is defined with the class TStyle and the
//  Color palette inside. We use the command gStyle->SetPalette(1,0)
//  for the display. So we have to convert the value "wert" into
//  a color index that fits the color palette.
//  The range of the color palette is defined by the values fMinPhe
//  and fMaxRange. Between this values we have 50 color index, starting
//  with 0 up to 49.
//
Int_t MHCamera::GetColor(Float_t val, Float_t min, Float_t max, Bool_t islog)
{
    if (TMath::IsNaN(val)) // FIXME: gLog!
    {
        cout << "MHCamera::GetColor: " << GetName() << " <" << GetTitle() << "> - Color for val=NaN (Not a Number) requested... set val=min." << endl;
        val = min;
    }

    //
    //   first treat the over- and under-flows
    //
    const Int_t maxcolidx = kItemsLegend-1;

    if (val >= max)
        return fColors[maxcolidx];

    if (val <= min)
        return fColors[0];

    //
    // calculate the color index
    //
    Float_t ratio;
    if (islog && min>0)
        ratio = log10(val/min) / log10(max/min);
    else
        ratio = (val-min) / (max-min);

    const Int_t colidx = (Int_t)(ratio*maxcolidx + .5);
    return fColors[colidx];
}

// ------------------------------------------------------------------------
//
//  Change the text on the legend according to the range of the Display
//
void MHCamera::UpdateLegend(Float_t minphe, Float_t maxphe, Bool_t islog)
{
    const Float_t range = fGeomCam->GetMaxRadius();

    const Float_t H = 0.9*range;
    const Float_t h = 2./kItemsLegend;
    const Float_t offset = 0.04*range;

    const Float_t w = range/sqrt((float)(fNcells-2));

    TBox newbox;
    TText newtxt;
    newtxt.SetTextSize(0.025);
    newtxt.SetTextAlign(12);
#if ROOT_VERSION_CODE > ROOT_VERSION(3,01,06)
    newtxt.SetBit(/*kNoContextMenu|*/kCannotPick);
    newbox.SetBit(/*kNoContextMenu|*/kCannotPick);
#endif

    for (Int_t i=0; i<kItemsLegend+1; i+=3)
    {
        const Float_t pos = (Float_t)i/kItemsLegend;

        Float_t val;
        if (islog && minphe>0)
            val = pow(10, log10(maxphe/minphe)*pos) * minphe;
        else
            val = minphe + pos * (maxphe-minphe);

        newtxt.PaintText(range+1.5*w, H*(i*h-1)-offset, Form(val<1e6?"%5.1f":"%5.1e", val));
    }

    for (Int_t i=0; i<kItemsLegend; i++)
    {
        newbox.SetFillColor(fColors[i]);
        newbox.PaintBox(range, H*(i*h-1)-offset, range+w, H*((i+1)*h-1)-offset);
    }

    TArrow arr;
    arr.PaintArrow(-range*.9, -range*.9, -range*.6, -range*.9, 0.025);
    arr.PaintArrow(-range*.9, -range*.9, -range*.9, -range*.6, 0.025);

    TString text;
    text += (int)(range*.3);
    text += "mm";

    TText newtxt2;
    newtxt2.SetTextSize(0.04);
    newtxt2.PaintText(-range*.85, -range*.85, text);

    text = "";
    text += (float)((int)(range*.3*fGeomCam->GetConvMm2Deg()*10))/10;
    text += "\\circ";
    text = text.Strip(TString::kLeading);

    TLatex latex;
    latex.PaintLatex(-range*.85, -range*.75, 0, 0.04, text);
}

// ------------------------------------------------------------------------
//
// Save primitive as a C++ statement(s) on output stream out
//
void MHCamera::SavePrimitive(ofstream &out, Option_t *opt)
{
    cout << "MHCamera::SavePrimitive: Must be rewritten!" << endl;
    /*
    if (!gROOT->ClassSaved(TCanvas::Class()))
        fDrawingPad->SavePrimitive(out, opt);

    out << "   " << fDrawingPad->GetName() << "->SetWindowSize(";
    out << fDrawingPad->GetWw() << "," << fDrawingPad->GetWh() << ");" << endl;
    */
}

// ------------------------------------------------------------------------
//
// compute the distance of a point (px,py) to the Camera
// this functions needed for graphical primitives, that
// means without this function you are not able to interact
// with the graphical primitive with the mouse!!!
//
// All calcutations are done in pixel coordinates
//
Int_t MHCamera::DistancetoPrimitive(Int_t px, Int_t py)
{
    if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
        return TH1D::DistancetoPrimitive(px, py);

    for (Int_t i=0; i<fNcells-2; i++)
    {
        MHexagon hex((*fGeomCam)[i]);
        if (hex.DistancetoPrimitive(px, py)==0)
            return 0;
    }
    return 999999;
}

// ------------------------------------------------------------------------
//
// Function introduced  (31-01-03)  WILL BE REMOVED IN THE FUTURE! DON'T
// USE IT!
//
void MHCamera::SetPix(const Int_t idx, const Int_t color, Float_t min, Float_t max)
{
    fArray[idx+1] = color;
    SetUsed(idx);
}

// ------------------------------------------------------------------------
//
// Function introduced  (31-01-03)  WILL BE REMOVED IN THE FUTURE! DON'T
// USE IT!
//
Int_t MHCamera::GetPixelIndex(Int_t px, Int_t py) const
{
    Int_t i;
    for (i=0; i<fNcells-2; i++)
    {
        MHexagon hex((*fGeomCam)[i]);
        if (hex.DistancetoPrimitive(px, py)>0)
            continue;

        return i;
    }
    return -1;
}

// ------------------------------------------------------------------------
//
// Returns string containing info about the object at position (px,py).
// Returned string will be re-used (lock in MT environment).
//
char *MHCamera::GetObjectInfo(Int_t px, Int_t py) const
{
    if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
        return TH1D::GetObjectInfo(px, py);

    static char info[128];

    const Int_t idx=GetPixelIndex(px, py);

    if (idx<0)
        return TObject::GetObjectInfo(px, py);

    sprintf(info, "Software Pixel Index: %d (Hardware Id=%d)", idx, idx+1);
    return info;
}

// ------------------------------------------------------------------------
//
// Execute a mouse event on the camera
//
void MHCamera::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
    if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
    {
        TH1D::ExecuteEvent(event, px, py);
        return;
    }
    //if (event==kMouseMotion && fStatusBar)
    //    fStatusBar->SetText(GetObjectInfo(px, py), 0);
    if (event!=kButton1Down)
        return;

    const Int_t idx = GetPixelIndex(px, py);
    if (idx<0)
        return;

    cout << GetTitle() << " <" << GetName() << ">" << endl;
    cout << "Software Pixel Index: " << idx << endl;
    cout << "Hardware Pixel Id:    " << idx+1 << endl;
    cout << "Contents:             " << fArray[idx+1] << "  <";
    cout << (IsUsed(idx)?"on":"off");
    cout << ">" << endl;

    if (fNotify && fNotify->GetSize()>0)
        new TCanvas;
    fNotify->ForEach(MCamEvent, DrawPixelContent)(idx);
}