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

/* ======================================================================== *\
!
! *
! * 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 "MPedestalCam.h"
#include "MPedestalPix.h"

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

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


ClassImp(MHSigmaTheta);


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

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

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

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

// --------------------------------------------------------------------------
//
// 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 << err << "MMcEvt not found... aborting." << endl;
       return kFALSE;
     }
  
   fPed = (MPedestalCam*)plist->FindObject("MPedestalCam");
   if (!fPed)
     {
       *fLog << err << "MPedestalCam 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");
   if (!binstheta)
     {
       *fLog << err << "BinningTheta [MBinning] not found... aborting." << endl;
       return kFALSE;      
     }

   // Get Sigmabar binning  
   MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar");
   if (!binssigma)
     {
       *fLog << err << "BinningSigmabar [MBinning] not found... aborting." << endl;
       return kFALSE;      
     }

   // Get binning for (sigma^2-sigmabar^2)  
   MBinning* binsdiff  = (MBinning*)plist->FindObject("BinningDiffsigma2");
   if (!binsdiff)
     {
       *fLog << err << "BinningDiffsigma2 [MBinning] not found... aborting." << endl;
       return kFALSE;      
     }


   // Set binnings in histograms
   SetBinning(&fSigmaTheta, binstheta, binssigma);

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

   MBinning binspix("BinningPixel");
   binspix.SetEdges(npix+1, -0.5, 0.5+npix );

   SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
   SetBinning(&fDiffPixTheta,  binstheta, &binspix, binsdiff);

   return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Fill the histograms
//
Bool_t MHSigmaTheta::Fill(const MParContainer *par)
{
    Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
    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 MPedestalPix &pix = (*fPed)[id];

      const Double_t sigma = pix.GetMeanRms();
      const Double_t area  = fCam->GetPixRatio(id);


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

      const Double_t diff = sigma*sigma/area - mySig*mySig;
      fDiffPixTheta.Fill(theta, (Double_t)id, diff);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Draw a copy of the histogram
//
TObject *MHSigmaTheta::DrawClone(Option_t *opt) 
{
    return MH::DrawClone(opt, 900, 900);
}

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

    AppendPad("");

    pad->Divide(3, 3);

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

    pad->cd(1);
    gPad->SetLogy();
    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("box");
    h->SetBit(kCanDelete);

    pad->cd(2);
    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->SetBit(kCanDelete);

    pad->cd(3);
    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(4);
    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);
    h = fDiffPixTheta.Project3D("zy");
    h->SetDirectory(NULL);
    h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)");
    h->SetXTitle("Id");
    h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}");
    h->Draw("box");
    h->SetBit(kCanDelete);

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

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

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

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