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

/* ======================================================================== *\
!
! *
! * 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 1D-histogram   ID of blind pixels
//             - the 1D-histogram   no.of pixels in MCerPhotEvt, and        //
//             - 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";

    fNpix.SetDirectory(NULL);
    fNpix.SetName("1D-Npix");
    fNpix.SetTitle("1D : no.of used pixels in MCerPhotEvt");
    fNpix.SetXTitle("no. of used pixels in MCerPhotEvt");
    fNpix.SetYTitle("Counts");

    fBlindId.SetDirectory(NULL);
    fBlindId.SetName("1D-BlindId");
    fBlindId.SetTitle("1D : Id of blind pixel");
    fBlindId.SetXTitle("Id of blind pixel");
    fBlindId.SetYTitle("Counts");

    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 << "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(&fNpix,    &binspix);
   SetBinning(&fBlindId, &binspix);
   SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
   SetBinning(&fDiffPixTheta,  binstheta, &binspix, 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);

    //*fLog << "MHSigmaTheta::Fill; theta, mySig = " << theta << ",  "
    //      << mySig << endl;

    fSigmaTheta.Fill(theta, mySig);

    const UInt_t npix = fEvt->GetNumPixels();

    Int_t npixused=0;
    for (UInt_t i=0; i<npix; i++)
    {
      MCerPhotPix cerpix = fEvt->operator[](i);
      if (!cerpix.IsPixelUsed())
      {
        fBlindId.Fill( cerpix.GetPixId() );        
        continue;
      }
      npixused++;

      /*      
      if (cerpix.GetNumPhotons() == 0)
        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);

      //*fLog << "theta, id, sigma, diff, area  = " << theta << ",  " << id 
      //      << ",  "
      //      << sigma << ",  " << diff << ",  " << area << endl;
    }
    fNpix.Fill(npixused);

    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);

    c.cd(8);
    ((TH1*)&fBlindId)->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.cd(9);
    ((TH1*)&fNpix)->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();
}
// --------------------------------------------------------------------------