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

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

#include <TCanvas.h>

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

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

#include "MGeomCam.h"

ClassImp(MHCamEvent);

using namespace std;

// --------------------------------------------------------------------------
//
// Initialize the name and title of the task.
// Resets the sum histogram
//
MHCamEvent::MHCamEvent(const char *name, const char *title)
    : fSum(NULL), fEvt(NULL), fType(0)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHCamEvent";
    fTitle = title ? title : "Average of MCamEvents";
}

// --------------------------------------------------------------------------
//
// 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 << warn << GetDescriptor() << ": No MCamEvent available..." << endl;
    }

    MGeomCam *cam = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!cam)
        *fLog << warn << GetDescriptor() << ": No MGeomCam found." << endl;

    if (fSum)
        delete (fSum);

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

    fSum = new MHCamera(*cam, name+";avg", fTitle);
    fSum->SetYTitle("a.u.");

    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 << "No MCamEvent found..." << endl;
        return kFALSE;
    }
    fSum->AddCamContent(*evt, fType);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Scale the sum container with the number of entries
//
Bool_t MHCamEvent::Finalize()
{
    if (fSum->GetEntries()>0)
        fSum->Scale(1./fSum->GetEntries());
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Print pixel indices which value is higher than mean+sigma*rms
//
void MHCamEvent::PrintOutliers(Float_t sigma) const
{
    const Double_t mean = fSum->GetMean();
    const Double_t rms  = fSum->GetRMS();

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

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

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

// --------------------------------------------------------------------------
//
// Return fSum.
//
TH1 *MHCamEvent::GetHistByName(const TString name)
{
    return fSum;
}

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

    pad->Divide(1,2);

    pad->cd(1);
    gPad->SetBorderMode(0);
    gPad->Divide(1,1);
    gPad->cd(1);
    gPad->SetBorderMode(0);
    fSum->Draw();

    pad->cd(2);
    gPad->SetBorderMode(0);
    fSum->Draw("EPhist");
}