source: trunk/MagicSoft/Mars/mhist/MHCamEvent.cc@ 7867

/* ======================================================================== *\
!
! *
! * 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, 12/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHCamEvent
//
// A histogram to store the sum of camera events. This can be photons,
// currents or enything else derived from MCamEvent
//
// Setup
// =====
//
// To count how often a certain pixel is above or below a threshold do:
//    MHCamEvent::SetThreshold(5.5);  // Default=LowerBound
//    MHCamEvent::SetThreshold(5.5, MHCamEvent::kIsUpperBound);
//
//
// Axis titles
// ===========
//
// 1) If no other title is given 'a.u.' is used.
// 2) If the title of MHCamEvent is different from the default,
//    it is used as histogram title. You can use this to set the
//    axis title, too. For more information see TH1::SetTitle, eg.
//       SetTitle("MyHist;;y[cm];Counts");
//    Make sure to set an empty x-axis title.
//
// For example:
//   MHCamEvent myhist("Titele;;y [cm]");
//
/////////////////////////////////////////////////////////////////////////////
#include "MHCamEvent.h"

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

#include "MParList.h"
#include "MCamEvent.h"
#include "MHCamera.h"

#include "MGeomCam.h"

ClassImp(MHCamEvent);

using namespace std;

const TString MHCamEvent::gsDefName  = "MHCamEvent";
const TString MHCamEvent::gsDefTitle = "Average of MCamEvents";

// --------------------------------------------------------------------------
//
// Initialize the name and title of the task
//
void MHCamEvent::Init(const char *name, const char *title)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    fNameGeom = "MGeomCam";
}

// --------------------------------------------------------------------------
//
// Initialize the name and title of the task. Set fType to 0
//
MHCamEvent::MHCamEvent(const char *name, const char *title)
    : fSum(NULL), fEvt(NULL), fType(0), fThreshold(0), fUseThreshold(0)
{
    Init(name, title);
}

// --------------------------------------------------------------------------
//
// Initialize the name and title of the task. Set fType to type
//
MHCamEvent::MHCamEvent(Int_t type, const char *name, const char *title)
    : fSum(NULL), fEvt(NULL), fType(type), fThreshold(0), fUseThreshold(0)
{
    Init(name, title);
}

// --------------------------------------------------------------------------
//
// Delete the corresponding camera display if available
//
MHCamEvent::~MHCamEvent()
{
    if (fSum)
        delete fSum;
}

// --------------------------------------------------------------------------
//
// Get the event (MCamEvent) the histogram might be filled with. If
// it is not given, it is assumed, that it is filled with the argument
// of the Fill function.
// Looks for the camera geometry MGeomCam and resets the sum histogram.
//
Bool_t MHCamEvent::SetupFill(const MParList *plist)
{
    fEvt = (MCamEvent*)plist->FindObject(fNameEvt, "MCamEvent");
    if (!fEvt)
    {
        if (!fNameEvt.IsNull())
        {
            *fLog << err << GetDescriptor() << ": No " << fNameEvt <<" [MCamEvent] available..." << endl;
            return kFALSE;
        }
        *fLog << inf << GetDescriptor() << ": Assuming to get MCamEvent from Fill()..." << endl;
    }

    MGeomCam *cam = (MGeomCam*)plist->FindObject(fNameGeom, "MGeomCam");
    if (!cam)
    {
        *fLog << err << fNameGeom << " [MGeomCam] not found... abort." << endl;
        return kFALSE;
    }

    // Delete a possible old histogram from a previous loop
    if (fSum)
        delete (fSum);

    // combine name
    const TString name = fNameEvt.IsNull() ? fName : fNameEvt;

    // create and setup MHCamera
    fSum = new MHCamera(*cam, name+";avg");
    if (fTitle!=gsDefTitle)
        fSum->SetTitle(fTitle);
    if (!fTitle.Contains(";"))
        fSum->SetYTitle("a.u.");
    fSum->SetBit(MHCamera::kProfile);

    fSum->SetXTitle("Pixel Idx");

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histograms with data from a MCamEvent-Container.
//
Bool_t MHCamEvent::Fill(const MParContainer *par, const Stat_t w)
{
    const MCamEvent *evt = par ? dynamic_cast<const MCamEvent*>(par) : fEvt;
    if (!evt)
    {
        *fLog << err << dbginf << "Got no MCamEvent as argument of Fill()..." << endl;
        return kFALSE;
    }
    if (fUseThreshold)
        fSum->CntCamContent(*evt, fThreshold, fType, fUseThreshold>0);
    else
        fSum->AddCamContent(*evt, fType);
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Take the mean of the sum histogram and print all pixel indices
// which are above sum+s*rms
//
void MHCamEvent::PrintOutliers(Float_t s) const
{
    const Double_t mean = fSum->GetMean();
    const Double_t rms  = fSum->GetRMS();

    *fLog << all << underline << GetDescriptor() << ": Mean=" << mean << ", Rms=" << rms << endl;

    for (UInt_t i=0; i<fSum->GetNumPixels(); i++)
    {
        if (!fSum->IsUsed(i))
            continue;

        if ((*fSum)[i+1]>mean+s*rms)
            *fLog << "Contents of Pixel-Index #" << i << ": " << (*fSum)[i+1] << " > " << s << "*rms" << endl;
    }
}

// --------------------------------------------------------------------------
//
// Return fSum for "sum" and fRms for "rms"
//
TH1 *MHCamEvent::GetHistByName(const TString name) const
{
    return fSum;
}

void MHCamEvent::Paint(Option_t *)
{
    TVirtualPad *pad = gPad;

    pad->cd(2);
    if (gPad->FindObject(Form("%s;proj", fName.Data())))
    {
        TH1 *h=fSum->Projection(Form("%s;proj", fName.Data()));
        if (h->GetMaximum()>0)
            gPad->SetLogy();
    }

    pad->cd(5);
    if (gPad->FindObject(Form("%s;rad", fName.Data())))
        fSum->RadialProfile(Form("%s;rad", fName.Data()));

    pad->cd(6);
    if (gPad->FindObject(Form("%s;az", fName.Data())))
        fSum->AzimuthProfile(Form("%s;az", fName.Data()));

    pad->cd(4);
    gPad->cd(1);
    MHCamera *cam = (MHCamera*)gPad->FindObject(Form("%s;err", fName.Data()));
    if (cam)
        cam->SetCamContent(*fSum, 1);
}

void MHCamEvent::Draw(Option_t *)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    const Int_t col = pad->GetFillColor();
    pad->SetBorderMode(0);

    AppendPad();

    TString name = Form("%s_1", pad->GetName());
    TPad *p = new TPad(name,name,0.005, 0.5, 0.66, 0.995,col,0,0);
    p->SetNumber(1);
    p->Draw();
    p->cd();
    fSum->Draw("EPhist");

    pad->cd();
    name = Form("%s_2", pad->GetName());
    p = new TPad(name,name,0.66, 0.5, 0.995, 0.995,col,0,0);
    p->SetNumber(2);
    p->Draw();
    p->cd();
    TH1 *h = fSum->Projection(Form("%s;proj", fName.Data()), 50);
    h->SetTitle("Projection");
    h->SetBit(kCanDelete);
    h->Draw();

    pad->cd();
    name = Form("%s_3", pad->GetName());
    p = new TPad(name,name,0.005, 0.005, 3./8, 0.5,col,0,0);
    p->SetNumber(3);
    p->Draw();
    p->cd();
    fSum->Draw();

    pad->cd();
    name = Form("%s_4", pad->GetName());
    p = new TPad(name,name,3./8, 0.005, 6./8-0.005, 0.5,col,0,0);
    p->SetNumber(4);
    p->Draw();
    p->cd();

    MHCamera *cam = new MHCamera(*fSum->GetGeometry());
    cam->SetName(Form("%s;err", fName.Data()));
    cam->SetTitle("Sqrt(Variance)");
    cam->SetYTitle(fSum->GetYaxis()->GetTitle());
    cam->SetCamContent(*fSum, 1);
    cam->SetBit(kCanDelete);
    cam->Draw();

    pad->cd();
    name = Form("%s_5", pad->GetName());
    p = new TPad(name,name,6./8,0.25,0.995,0.5,col,0,0);
    p->SetNumber(5);
    p->Draw();
    p->cd();
    h = (TH1*)fSum->RadialProfile(Form("%s;rad", fName.Data()), 20);
    h->SetTitle("Radial Profile");
    h->SetBit(kCanDelete|TH1::kNoStats);
    h->Draw();

    pad->cd();
    name = Form("%s_6", pad->GetName());
    p = new TPad(name,name,6./8,0.005,0.995,0.25,col,0,0);
    p->SetNumber(6);
    p->Draw();
    p->cd();
    h = (TH1*)fSum->AzimuthProfile(Form("%s;az", fName.Data()), 30);
    h->SetTitle("Azimuth Profile");
    h->SetBit(kCanDelete|TH1::kNoStats);
    h->Draw();
}