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

/* ======================================================================== *\
!
! *
! * 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>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  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 "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"

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

#include "MSigmabar.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 : Sigmabar, \\Theta");
    fSigmaTheta.SetXTitle("\\Theta [\\circ]");
    fSigmaTheta.SetYTitle("Sigmabar");

    fSigmaPixTheta.SetDirectory(NULL);
    fSigmaPixTheta.SetName("3D-ThetaPixSigma");
    fSigmaPixTheta.SetTitle("3D : \\Theta, pixel no., Sigma");
    fSigmaPixTheta.SetXTitle("\\Theta [\\circ]");
    fSigmaPixTheta.SetYTitle("pixel number");
    fSigmaPixTheta.SetZTitle("Sigma");

    fDiffPixTheta.SetDirectory(NULL);
    fDiffPixTheta.SetName("3D-ThetaPixDiff");
    fDiffPixTheta.SetTitle("3D : \\Theta, pixel, Sigma^2-Sigmabar^2");
    fDiffPixTheta.SetXTitle("\\Theta [\\circ]");
    fDiffPixTheta.SetYTitle("pixel number");
    fDiffPixTheta.SetZTitle("Sigma^2-sigmabar^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 << dbginf << "MHSigmaTheta::SetupFill : MMcEvt not found... aborting." << endl;
       return kFALSE;
     }
  
   fPed = (MPedestalCam*)plist->FindObject("MPedestalCam");
   if (!fPed)
     {
       *fLog << dbginf << "MHSigmaTheta::SetupFill : MPedestalCam not found... aborting." << endl;
       return kFALSE;
     }
  
   fCam = (MGeomCam*)plist->FindObject("MGeomCam");
   if (!fCam)
     {
       *fLog << dbginf << "MHSigmaTheta::SetupFill : MGeomCam not found (no geometry information available)... aborting." << endl;
       return kFALSE;
     }
  
   fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
   if (!fEvt)
     {
       *fLog << dbginf << "MHSigmaTheta::SetupFill : MCerPhotEvt not found... aborting." << endl;
       return kFALSE;
     }

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

   // Get Theta Binning  
   MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta");
   if (!binstheta)
     {
       *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningTheta not found... aborting." << endl;
       return kFALSE;      
     }

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

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


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

   // Get binning for pixel number
   UInt_t npix = fPed->GetSize();
   MBinning binspix("BinningPixel");
   MBinning* binspixel = &binspix;
   binspixel->SetEdges(npix, -0.5, ((float)npix)-0.5 ); 

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

   return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Fill the histograms
//
Bool_t MHSigmaTheta::Fill(const MParContainer *par)
{
  //*fLog << "entry Fill" << endl;

    Double_t Theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
    Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
    fSigmaTheta.Fill(Theta, mySig);

    //*fLog << "Theta, mySig = " << Theta << ",  " << mySig << endl;

    const UInt_t npix = fEvt->GetNumPixels();
    for (UInt_t i=0; i<npix; i++)
    {
      MCerPhotPix cerpix = fEvt->operator[](i);
      if (!cerpix.IsPixelUsed())
        continue;

      if (cerpix.GetNumPhotons() == 0)
        continue;

      Int_t j = cerpix.GetPixId();
      const MPedestalPix pix = fPed->operator[](j);

      Double_t Sigma = pix.GetMeanRms();
      Double_t Area  = fCam->GetPixRatio(j);

      fSigmaPixTheta.Fill(Theta, (Double_t)j, Sigma);

      Double_t Diff = Sigma*Sigma/Area - mySig*mySig;
      fDiffPixTheta.Fill (Theta, (Double_t)j, Diff);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Plot the results
//
Bool_t MHSigmaTheta::Finalize()
{
    DrawClone();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Draw a copy of the histogram
//
TObject *MHSigmaTheta::DrawClone(Option_t *opt) 
{
    TCanvas &c = *MakeDefCanvas("SigmaThetaPlot", "Sigmabar vs. Theta",
                                 900, 900);
    c.Divide(3, 3);

    gROOT->SetSelectedPad(NULL);

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

    c.cd(1);
    h = ((TH2*)&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->DrawCopy(opt);
    h->SetBit(kCanDelete);;
    gPad->SetLogy();

    c.cd(4);
    h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E");
    h->SetDirectory(NULL);
    h->SetTitle("Distribution of Sigmabar");
    h->SetXTitle("Sigmabar");
    h->SetYTitle("No.of events");

    h->DrawCopy(opt);
    h->SetBit(kCanDelete);;

    c.cd(7);
    ((TH2*)&fSigmaTheta)->DrawCopy(opt);

    //--------------------------------------------------------------------
    // draw the 3D histogram : Theta, pixel, Sigma^2-Sigmabar^2

    TH2D *l;

    c.cd(2);
    l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx");
    l->SetDirectory(NULL);
    l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
    l->SetXTitle("\\Theta [\\circ]");
    l->SetYTitle("Sigma^2-Sigmabar^2");

    l->DrawCopy("box");
    l->SetBit(kCanDelete);;

    c.cd(5);
    l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy");
    l->SetDirectory(NULL);
    l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
    l->SetXTitle("pixel");
    l->SetYTitle("Sigma^2-Sigmabar^2");

    l->DrawCopy("box");
    l->SetBit(kCanDelete);;

    c.cd(8);
    ((TH2*)&fDiffPixTheta)->DrawCopy(opt);


    //--------------------------------------------------------------------
    // draw the 3D histogram : Theta, pixel, Sigma

    TH2D *k;

    c.cd(3);
    k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx");
    k->SetDirectory(NULL);
    k->SetTitle("Sigma vs. \\Theta (all pixels)");
    k->SetXTitle("\\Theta [\\circ]");
    k->SetYTitle("Sigma");

    k->DrawCopy("box");
    k->SetBit(kCanDelete);;

    c.cd(6);
    k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy");
    k->SetDirectory(NULL);
    k->SetTitle("Sigma vs. pixel number (all \\Theta)");
    k->SetXTitle("pixel");
    k->SetYTitle("Sigma");

    k->DrawCopy("box");
    k->SetBit(kCanDelete);;

    c.cd(9);
    ((TH2*)&fSigmaPixTheta)->DrawCopy(opt);

    //--------------------------------------------------------------------
    c.Modified();
    c.Update();

    return &c;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHSigmaTheta::Draw(Option_t *opt)
{
    if (!gPad)
        MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta",
                       600, 600);

    TH1D *h;

    gPad->Divide(2,2);

    gPad->cd(1);
    h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E");
    h->SetTitle("Distribution of \\Theta");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("No.of events");

    h->Draw(opt);
    h->SetBit(kCanDelete);;
    gPad->SetLogy();

    gPad->cd(2);
    h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E");
    h->SetTitle("Distribution of Sigmabar");
    h->SetXTitle("Sigmabar");
    h->SetYTitle("No.of events");

    h->Draw(opt);
    h->SetBit(kCanDelete);;

    gPad->cd(3);
    fSigmaTheta.DrawCopy(opt);



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