source: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc@ 2951

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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHSigmaTheta (extension of Robert's MHSigmabarTheta)                    //
//                                                                          //
//  calculates - the 2D-histogram   sigmabar vs. Theta, and                 //
//             - the 3D-histogram   sigma, pixel no., Theta                 //
//             - the 3D-histogram   (sigma^2-sigmabar^2), pixel no., Theta  //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHSigmaTheta.h"

#include <TCanvas.h>

#include "MTime.h"
#include "MMcEvt.hxx"

#include "MBinning.h"
#include "MParList.h"
#include "MSigmabar.h"

#include "MGeomCam.h"

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

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

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

ClassImp(MHSigmaTheta);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title of the histogram.
//
MHSigmaTheta::MHSigmaTheta(const char *name, const char *title)
{
    fName  = name  ? name  : "MHSigmaTheta";
    fTitle = title ? title : "2D histogram sigmabar vs. Theta";

    fSigmaTheta.SetDirectory(NULL);
    fSigmaTheta.SetName("2D-ThetaSigmabar");
    fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta");
    fSigmaTheta.SetXTitle("\\Theta [\\circ]");
    fSigmaTheta.SetYTitle("Sigmabar");

    fSigmaPixTheta.SetDirectory(NULL);
    fSigmaPixTheta.SetName("3D-ThetaPixSigma");
    fSigmaPixTheta.SetTitle("3D: \\Theta, Pixel Id, \\sigma");
    fSigmaPixTheta.SetXTitle("\\Theta [\\circ]");
    fSigmaPixTheta.SetYTitle("Pixel Id");
    fSigmaPixTheta.SetZTitle("\\sigma");

    fDiffPixTheta.SetDirectory(NULL);
    fDiffPixTheta.SetName("3D-ThetaPixDiff");
    fDiffPixTheta.SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}");
    fDiffPixTheta.SetXTitle("\\Theta [\\circ]");
    fDiffPixTheta.SetYTitle("Pixel Id");
    fDiffPixTheta.SetZTitle("{\\sigma}^{2}-\\bar{\\sigma}^{2}");

    // Set default binning
    // FIXME: Maybe ist's necessary to adapt the value to the
    // Magic default values
    MBinning binsb;
    MBinning binst;
    MBinning binsd;
    binsd.SetEdges(100, -5, 20);
    binsb.SetEdges(100, 0, 5);
    binst.SetEdgesCos(10, 0, 90);

    MBinning binspix("BinningPixel");
    binspix.SetEdges(578, -0.5, 577.5);

    SetBinning(&fSigmaTheta,    &binst, &binsb);
    SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb);
    SetBinning(&fDiffPixTheta,  &binst, &binspix, &binsd);
}

// --------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histogram
//
Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
{
    fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
    if (!fMcEvt)
        *fLog << warn << "MMcEvt not found... aborting." << endl;

    fPed = (MPedPhotCam*)plist->FindObject("MPedPhotCam");
    if (!fPed)
    {
        *fLog << err << "MPedPhotCam not found... aborting." << endl;
        return kFALSE;
    }

    fCam = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

    fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
    if (!fEvt)
    {
        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar");
    if (!fSigmabar)
    {
        *fLog << err << "MSigmabar not found... aborting." << endl;
        return kFALSE;
    }

    // Get Theta Binning
    MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
    if (!binstheta)
    {
        *fLog << warn << "Object 'BinningTheta' [MBinning] not found... no binning applied." << endl;
        return kTRUE;
    }

    // Get Sigmabar binning
    MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar", "MBinning");
    if (!binssigma)
        *fLog << warn << "Object 'BinningSigmabar' [MBinning] not found... no binning applied." << endl;

    // Get binning for (sigma^2-sigmabar^2)
    MBinning* binsdiff  = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
    if (!binsdiff)
        *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... no binning applied." << endl;

    //FIXME: Missing: Apply new binning to only one axis!

    // Get binning for pixel number
    const UInt_t npix1 = fPed->GetSize()+1;

    MBinning binspix("BinningPixel");
    binspix.SetEdges(npix1, -0.5, npix1-0.5);

    // Set binnings in histograms
    if (binssigma)
    {
        SetBinning(&fSigmaTheta, binstheta, binssigma);
        SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
    }

    if (binsdiff)
        SetBinning(&fDiffPixTheta,  binstheta, &binspix, binsdiff);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Fill the histograms
//
Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
{
    Double_t theta = fMcEvt ? fMcEvt->GetTelescopeTheta()*kRad2Deg : 0;
    Double_t mysig = fSigmabar->Calc(*fCam, *fPed, *fEvt);

    fSigmaTheta.Fill(theta, mysig);

    const UInt_t npix = fEvt->GetNumPixels();

    for (UInt_t i=0; i<npix; i++)
    {
        MCerPhotPix cerpix = (*fEvt)[i];
        if (!cerpix.IsPixelUsed())
            continue;

        const Int_t id = cerpix.GetPixId();
        const MPedPhotPix &pix = (*fPed)[id];

        // ratio is the area of pixel 0 
        //          divided by the area of the current pixel
        const Double_t ratio = fCam->GetPixRatio(id);
        const Double_t sigma = pix.GetRms();

        fSigmaPixTheta.Fill(theta, (Double_t)id, sigma);

        const Double_t diff = sigma*sigma*ratio - mysig*mysig;
        fDiffPixTheta.Fill(theta, (Double_t)id, diff);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Update the projections and (if possible) set log scales before painting
//
void MHSigmaTheta::Paint(Option_t *opt)
{
    TVirtualPad *padsave = gPad;

    TH1D* h;

    padsave->cd(1);
    if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
    {
        ProjectionX(*h, fSigmaTheta);
        if (h->GetEntries()!=0)
            gPad->SetLogy();
    }

    padsave->cd(4);
    if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
        ProjectionY(*h, fSigmaTheta);

    gPad = padsave;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHSigmaTheta::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("");

    pad->Divide(3, 2);

    // draw the 2D histogram Sigmabar versus Theta
    TH1 *h;

    pad->cd(1);
    gPad->SetBorderMode(0);
    h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E");
    h->SetDirectory(NULL);
    h->SetTitle("Distribution of \\Theta");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("No.of events");
    h->Draw(opt);
    h->SetBit(kCanDelete);

    pad->cd(2);
    gPad->SetBorderMode(0);
    h = fDiffPixTheta.Project3D("zx");
    h->SetDirectory(NULL);
    h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(3);
    gPad->SetBorderMode(0);
    h = fSigmaPixTheta.Project3D("zx");
    h->SetDirectory(NULL);
    h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("\\sigma_{ped}");
    h->Draw("box");
    h->SetBit(kCanDelete);

    //pad->cd(7);
    //gPad->SetBorderMode(0);
    //h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E");
    //h->SetDirectory(NULL);
    //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
    //h->SetXTitle("\\bar{\\sigma}_{ped}");
    //h->SetYTitle("No.of events");
    //h->Draw(opt);
    //h->SetBit(kCanDelete);

    pad->cd(5);
    gPad->SetBorderMode(0);
    h = fDiffPixTheta.Project3D("zy");
    h->SetDirectory(NULL);
    h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all  \\Theta)");
    h->SetXTitle("Pixel Id");
    h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(6);
    gPad->SetBorderMode(0);
    h = fSigmaPixTheta.Project3D("zy");
    h->SetDirectory(NULL);
    h->SetTitle("\\sigma_{ped} vs. pixel Id (all  \\Theta)");
    h->SetXTitle("Pixel Id");
    h->SetYTitle("\\sigma_{ped}");
    h->Draw("box");
    h->SetBit(kCanDelete);

    pad->cd(4);
    fSigmaTheta.Draw(opt);

    //pad->cd(8);
    //fDiffPixTheta.Draw(opt);

    //pad->cd(9);
    //fSigmaPixTheta.Draw(opt);
}