source: trunk/MagicSoft/Mars/mimage/MHVsSize.cc@ 6979

/* ======================================================================== *\
!
! *
! * 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@astro.uni-wuerzburg.de>
!   Author(s): Wolfgang Wittek  2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHVsSize
//
// This class contains histograms for the source independent image parameters
//
/////////////////////////////////////////////////////////////////////////////
#include "MHVsSize.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 "MHillasSrc.h"
#include "MHillasExt.h"
#include "MNewImagePar.h"
#include "MGeomCam.h"
#include "MBinning.h"

ClassImp(MHVsSize);

using namespace std;

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

    fLength.SetNameTitle("Length", "Length vs. Size");
    fWidth.SetNameTitle( "Width",  "Width vs. Size");
    fDist.SetNameTitle(  "Dist",   "Dist vs. Size");
    fConc1.SetNameTitle( "Conc1",  "Conc1 vs. Size");
    fArea.SetNameTitle(  "Area",   "Area vs. Size");
    fM3Long.SetNameTitle("M3Long", "M3Long vs. Size");

    fLength.SetDirectory(NULL);
    fWidth.SetDirectory(NULL);
    fDist.SetDirectory(NULL);
    fConc1.SetDirectory(NULL);
    fArea.SetDirectory(NULL);
    fM3Long.SetDirectory(NULL);

    fLength.SetXTitle("Size [phe]");
    fWidth.SetXTitle("Size [phe]");
    fDist.SetXTitle("Size [phe]");
    fConc1.SetXTitle("Size [phe]");
    fArea.SetXTitle("Size [phe]");
    fM3Long.SetXTitle("Size [phe]");

    fLength.SetYTitle("Length [mm]");
    fWidth.SetYTitle("Width [mm]");
    fDist.SetYTitle("Distance [mm]");
    fConc1.SetYTitle("Conc1 [ratio]");
    fArea.SetYTitle("Conc1 [mm^{2}]");
    fM3Long.SetYTitle("M3Long [mm]");

    MBinning binse, binsl, binsd, binsc, binsa, binsm;
    binse.SetEdgesLog( 50,   10, 1e5);
    binsl.SetEdges(   100,    0, 296.7/2);
    binsd.SetEdges(   100,    0, 600);
    binsc.SetEdgesLog(100, 1e-5, 5e-3);
    binsa.SetEdges(   100,    0, 445*45);
    binsm.SetEdges(   100, -445, 445);

    MH::SetBinning(&fLength, &binse, &binsl);
    MH::SetBinning(&fWidth,  &binse, &binsl);
    MH::SetBinning(&fDist,   &binse, &binsd);
    MH::SetBinning(&fConc1,  &binse, &binsc);
    MH::SetBinning(&fArea,   &binse, &binsa);
    MH::SetBinning(&fM3Long, &binse, &binsm);

    fLength.UseCurrentStyle();
    fWidth.UseCurrentStyle();
    fDist.UseCurrentStyle();
    fConc1.UseCurrentStyle();
    fArea.UseCurrentStyle();
    fM3Long.UseCurrentStyle();
}

// --------------------------------------------------------------------------
//
// 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 MHVsSize::SetupFill(const MParList *plist)
{
    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);
    }

    fHillas = (MHillas*)plist->FindObject("MHillas");
    if (!fHillas)
    {
        *fLog << err << "MHillas not found... abort." << endl;
        return kFALSE;
    }

    fHillasExt = (MHillasExt*)plist->FindObject("MHillasExt");
    if (!fHillasExt)
    {
        *fLog << err << "MHillasExt not found... abort." << endl;
        return kFALSE;
    }

    fNewImagePar = (MNewImagePar*)plist->FindObject("MNewImagePar");
    if (!fNewImagePar)
    {
        *fLog << err << "MNewImagePar not found... abort." << endl;
        return kFALSE;
    }

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

    const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
    if (!bins)
    {
        float r = fGeomCam ? fGeomCam->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 MHVsSize::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 MHVsSize::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, 1, scale);
    MH::ScaleAxis(&fWidth,  1, scale);
    MH::ScaleAxis(&fDist,   1, scale);
    MH::ScaleAxis(&fM3Long, 1, scale);
    MH::ScaleAxis(&fArea,   1, scale*scale);

    if (mmscale)
    {
        fLength.SetYTitle("Length [mm]");
        fWidth.SetYTitle("Width [mm]");
        fDist.SetYTitle("Distance [mm]");
        fArea.SetYTitle("Area [mm^{2}]");
        fM3Long.SetYTitle("M3Long [mm]");
    }
    else
    {
        fLength.SetYTitle("Length [\\circ]");
        fWidth.SetYTitle("Width [\\circ]");
        fDist.SetYTitle("Distance [\\circ]");
        fArea.SetYTitle("Area [\\circ^{2}]");
        fM3Long.SetYTitle("M3Long [\\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 MHVsSize::Fill(const MParContainer *par, const Stat_t w)
{
    const MHillasSrc *src = dynamic_cast<const MHillasSrc*>(par);
    if (!src)
    {
        *fLog << err << "MHVsSize::Fill: Wrong argument... abort." << endl;
        return kFALSE;
    }

    const Double_t scale = fUseMmScale ? 1 : fMm2Deg;

    fLength.Fill(fHillas->GetSize(), scale*fHillas->GetLength(),     w);
    fWidth.Fill( fHillas->GetSize(), scale*fHillas->GetWidth(),      w);
    fDist.Fill(  fHillas->GetSize(), scale*src->GetDist(),           w);
    fConc1.Fill( fHillas->GetSize(), scale*fNewImagePar->GetConc1(), w);
    fArea.Fill(  fHillas->GetSize(), scale*scale*fHillas->GetArea(), w);
    fM3Long.Fill(fHillas->GetSize(), scale*fHillasExt->GetM3Long()*TMath::Sign(1.0f, src->GetCosDeltaAlpha()), w);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// 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 MHVsSize::Draw(Option_t *o)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("");

    TString opt(o);
    opt.ToLower();

    // FIXME: If same-option given make two independant y-axis!
    const Bool_t same = opt.Contains("same");

    if (!same)
        pad->Divide(3,2);
    else
    {
        fLength.SetName("LengthSame");
        fWidth.SetName("WidthSame");
        fDist.SetName("DistSame");
        fConc1.SetName("Conc1Same");
        fArea.SetName("AreaSame");
        fM3Long.SetName("M3LongSame");

        fLength.SetDirectory(0);
        fWidth.SetDirectory(0);
        fDist.SetDirectory(0);
        fConc1.SetDirectory(0);
        fArea.SetDirectory(0);
        fM3Long.SetDirectory(0);

        fDist.SetMarkerColor(kBlue);
        fConc1.SetMarkerColor(kBlue);
        fWidth.SetMarkerColor(kBlue);
        fLength.SetMarkerColor(kBlue);
        fArea.SetMarkerColor(kBlue);
        fM3Long.SetMarkerColor(kBlue);
    }

    pad->cd(1);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    RemoveFromPad("LengthSame");
    fLength.Draw(same?"same":"");

    pad->cd(2);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    RemoveFromPad("WidthSame");
    fWidth.Draw(same?"same":"");

    pad->cd(3);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    RemoveFromPad("DistSame");
    fDist.Draw(same?"same":"");

    pad->cd(4);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    RemoveFromPad("AreaSame");
    fArea.Draw(same?"same":"");

    pad->cd(5);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    RemoveFromPad("M3LongSame");
    fM3Long.Draw(same?"same":"");

    pad->cd(6);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    RemoveFromPad("Conc1Same");
    fConc1.Draw(same?"same":"");
}