source: trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc@ 3913

/* ======================================================================== *\
!
! *
! * 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, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Markus Gaug, 02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//
//  MGCamDisplays
//
//  Graphical interfaces to display the camera with fits and projections
//
/////////////////////////////////////////////////////////////////////////////
#include "MGCamDisplays.h"

#include <TStyle.h>
#include <TCanvas.h>

#include "MH.h"
#include "MHCamera.h"
#include "MGeomCam.h"
#include "TVirtualPad.h"
#include "TProfile.h"
#include "TF1.h"

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

#include "MStatusDisplay.h"

ClassImp(MGCamDisplays);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MGCamDisplays::MGCamDisplays() 
{
}

// --------------------------------------------------------------------------
//
// Draw the MHCamera into the MStatusDisplay: 
// 
// 1) Draw it as histogram (MHCamera::DrawCopy("hist")
// 2) Draw it as a camera, with MHCamera::SetPrettyPalette() set. 
// 3) If "rad" is not zero, draw its values vs. the radius from the camera center. 
//    (DrawRadialProfile())
// 4) Depending on the variable "fit", draw the values projection on the y-axis
//    (DrawProjection()):
//    0: don't draw
//    1: Draw fit to Single Gauss (for distributions flat-fielded over the whole camera)
//    2: Draw and fit to Double Gauss (for distributions different for inner and outer pixels)
//    3: Draw and fit to Triple Gauss (for distributions with inner, outer pixels and outliers)
//    4: Draw and fit to Polynomial grade 0: (for the probability distributions)
//    >4: Draw and don;t fit.
//
void MGCamDisplays::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, 
                            const Int_t fit, const Int_t rad, TObject *notify)
{

    c.cd(x);
    gPad->SetBorderMode(0);
    gPad->SetTicks();
    MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
    obj1->SetDirectory(NULL);

    if (notify)
      obj1->AddNotify(notify);

    c.cd(x+y);
    gPad->SetBorderMode(0);
    obj1->SetPrettyPalette();
    obj1->Draw();

    if (rad)
      {
        c.cd(x+2*y);
        gPad->SetBorderMode(0);
        gPad->SetTicks();
        DrawRadialProfile(obj1);
      }
    

    if (!fit)
        return;

    c.cd(rad ? x+3*y : x+2*y);
    gPad->SetBorderMode(0);
    gPad->SetTicks();
    DrawProjection(obj1, fit);
}

// --------------------------------------------------------------------------
//
// Draw a projection of MHCamera onto the y-axis values. Depending on the 
// variable fit, the following fits are performed:
//
// 1: Single Gauss (for distributions flat-fielded over the whole camera)
// 2: Double Gauss (for distributions different for inner and outer pixels)
// 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
// 4: flat         (for the probability distributions)
// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are 
//        drawn separately, for inner and outer pixels.
// 5: Fit Inner and Outer pixels separately by a single Gaussian 
//                 (only for MAGIC cameras)
// 6: Fit Inner and Outer pixels separately by a single Gaussian and display 
//                 additionally the two camera halfs separately (for MAGIC camera)
//
//
void MGCamDisplays::DrawProjection(MHCamera *obj, Int_t fit) const
{
  
  TArrayI inner(1);
  inner[0] = 0;
  
  TArrayI outer(1);
  outer[0] = 1;
          
  if (fit==5 || fit==6)
    {
      
      if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
        {
          TArrayI s0(6);
          s0[0] = 6;
          s0[1] = 1;
          s0[2] = 2;
          s0[3] = 3;
          s0[4] = 4;
          s0[5] = 5;

          TArrayI s1(3);
          s1[0] = 6;
          s1[1] = 1;
          s1[2] = 2;
          
          TArrayI s2(3);
          s2[0] = 3;
          s2[1] = 4;
          s2[2] = 5;

          gPad->Clear();
          TVirtualPad *pad = gPad;
          pad->Divide(2,1);
          
          TH1D *inout[2];
          inout[0] = obj->ProjectionS(s0, inner, "Inner");
          inout[1] = obj->ProjectionS(s0, outer, "Outer");
          
          inout[0]->SetDirectory(NULL);
          inout[1]->SetDirectory(NULL);
          
          for (int i=0; i<2; i++)
            {
              pad->cd(i+1);
              inout[i]->SetLineColor(kRed+i);
              inout[i]->SetBit(kCanDelete);
              inout[i]->Draw();
              inout[i]->Fit("gaus","Q");

              if (fit == 6)
                {
                  TH1D *half[2];
                  half[0] = obj->ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
                  half[1] = obj->ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
                  
                  for (int j=0; j<2; j++)
                    {
                      half[j]->SetLineColor(kRed+i+j);
                      half[j]->SetDirectory(0);
                      half[j]->SetBit(kCanDelete);
                      half[j]->Draw("same");
                    }
                }
              
            }
          
          gLog << all << obj->GetName() 
               << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ",
                       inout[0]->GetFunction("gaus")->GetParameter(1),"+-",
                       inout[0]->GetFunction("gaus")->GetParError(1));
          gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
                       inout[1]->GetFunction("gaus")->GetParameter(1),"+-",
                       inout[1]->GetFunction("gaus")->GetParError(1)) << endl;
          gLog << all << obj->GetName()
               << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ",
                       inout[0]->GetFunction("gaus")->GetParameter(2),"+-",
                       inout[0]->GetFunction("gaus")->GetParError(2));
          gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
                       inout[1]->GetFunction("gaus")->GetParameter(2),"+-",
                       inout[1]->GetFunction("gaus")->GetParError(2)) << endl;

        }
      return;
    }
  
  TH1D *obj2 = (TH1D*)obj->Projection(obj->GetName());
  obj2->SetDirectory(0);
  obj2->Draw();
  obj2->SetBit(kCanDelete);
  
  
  if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    {
      TArrayI s0(3);
      s0[0] = 6;
      s0[1] = 1;
      s0[2] = 2;
      
      TArrayI s1(3);
      s1[0] = 3;
      s1[1] = 4;
      s1[2] = 5;
      
      
      TH1D *halfInOut[4];
      
      // Just to get the right (maximum) binning
      halfInOut[0] = obj->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
      halfInOut[1] = obj->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
      halfInOut[2] = obj->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
      halfInOut[3] = obj->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
      
      for (int i=0; i<4; i++)
        {
          halfInOut[i]->SetLineColor(kRed+i);
          halfInOut[i]->SetDirectory(0);
          halfInOut[i]->SetBit(kCanDelete);
          halfInOut[i]->Draw("same");
        }
    }
  
  const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
  const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
  const Double_t integ = obj2->Integral("width")/2.5;
  const Double_t mean  = obj2->GetMean();
  const Double_t rms   = obj2->GetRMS();
  const Double_t width = max-min;
  
  const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
    "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
  
  const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
    "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
    "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
  TF1 *f=0;
  switch (fit)
    {
    case 1:
      f = new TF1("sgaus", "gaus(0)", min, max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("Area", "#mu", "#sigma");
      f->SetParameters(integ/rms, mean, rms);
      f->SetParLimits(0, 0,   integ);
      f->SetParLimits(1, min, max);
      f->SetParLimits(2, 0,   width/1.5);
      
      obj2->Fit(f, "QLR");
      break;
      
    case 2:
      f = new TF1("dgaus",dgausformula.Data(),min,max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
      f->SetParameters(integ,(min+mean)/2.,width/4.,
                       integ/width/2.,(max+mean)/2.,width/4.);
      // The left-sided Gauss
      f->SetParLimits(0,integ-1.5      , integ+1.5);
      f->SetParLimits(1,min+(width/10.), mean);
      f->SetParLimits(2,0              , width/2.);
      // The right-sided Gauss
      f->SetParLimits(3,0   , integ);
      f->SetParLimits(4,mean, max-(width/10.));
      f->SetParLimits(5,0   , width/2.);
      obj2->Fit(f,"QLRM");
      break;
      
    case 3:
      f = new TF1("tgaus",tgausformula.Data(),min,max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
                     "A_{2}","#mu_{2}","#sigma_{2}",
                     "A_{3}","#mu_{3}","#sigma_{3}");
      f->SetParameters(integ,(min+mean)/2,width/4.,
                       integ/width/3.,(max+mean)/2.,width/4.,
                       integ/width/3.,mean,width/2.);
      // The left-sided Gauss
      f->SetParLimits(0,integ-1.5,integ+1.5);
      f->SetParLimits(1,min+(width/10.),mean);
      f->SetParLimits(2,width/15.,width/2.);
      // The right-sided Gauss
      f->SetParLimits(3,0.,integ);
      f->SetParLimits(4,mean,max-(width/10.));
      f->SetParLimits(5,width/15.,width/2.);
      // The Gauss describing the outliers
      f->SetParLimits(6,0.,integ);
      f->SetParLimits(7,min,max);
      f->SetParLimits(8,width/4.,width/1.5);
      obj2->Fit(f,"QLRM");
      break;
      
    case 4:
      obj2->Fit("pol0", "Q");
      obj2->GetFunction("pol0")->SetLineColor(kYellow);
      break;
      
    case 9:
      break;
      
    default:
      obj2->Fit("gaus", "Q");
      obj2->GetFunction("gaus")->SetLineColor(kYellow);
      break;
    }
}

// --------------------------------------------------------------------------
//
// Draw a projection of MHCamera vs. the radius from the central pixel. 
//
// The inner and outer pixels are drawn separately, both fitted by a polynomial
// of grade 1.
//
void MGCamDisplays::DrawRadialProfile(MHCamera *obj) const
{
  
  TProfile *obj2 = (TProfile*)obj->RadialProfile(obj->GetName());
  obj2->SetDirectory(0);
  obj2->Draw();
  obj2->SetBit(kCanDelete);
  
  if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    {

      TArrayI s0(6);
      s0[0] = 1;
      s0[1] = 2;
      s0[2] = 3;
      s0[3] = 4;
      s0[4] = 5;
      s0[5] = 6;

      TArrayI inner(1);
      inner[0] = 0;
      
      TArrayI outer(1);
      outer[0] = 1;
      
      // Just to get the right (maximum) binning
      TProfile *half[2];
      half[0] = obj->RadialProfileS(s0, inner,Form("%s%s",obj->GetName(),"Inner"));
      half[1] = obj->RadialProfileS(s0, outer,Form("%s%s",obj->GetName(),"Outer"));
      
      for (Int_t i=0; i<2; i++)
        {
          Double_t min = obj->GetGeomCam().GetMinRadius(i);
          Double_t max = obj->GetGeomCam().GetMaxRadius(i);

          half[i]->SetLineColor(kRed+i);
          half[i]->SetDirectory(0);
          half[i]->SetBit(kCanDelete);
          half[i]->Draw("same");
          half[i]->Fit("pol1","Q","",min,max);
          half[i]->GetFunction("pol1")->SetLineColor(kRed+i);
          half[i]->GetFunction("pol1")->SetLineWidth(1);
        }
    }
}