source: trunk/MagicSoft/Mars/mhist/MHMcEnergy.cc@ 863

/* ======================================================================== *\
!
! *
! * 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): Javier Lopez 05/2001 (jlopez@ifae.es)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MHMcEnergy
//
// This class holds the information (histogram and fit function)
// about the energy threshold for a particular trigger condition.
//
////////////////////////////////////////////////////////////////////////////
#include "MHMcEnergy.h" 

#include <iostream.h>

#include <TH1.h> 
#include <TF1.h> 
#include <TCanvas.h>
#include <TPaveLabel.h> 

ClassImp(MHMcEnergy);

// -------------------------------------------------------------------------
//
//  Default Constructor.
//
MHMcEnergy::MHMcEnergy(const UInt_t idx, const char *name, const char *title)
{ 
    char aux[15]="MHMcEnergy";

    if (idx>0)
        sprintf(aux+10, ";%i", idx);
    *fName  = name  ? name  : aux;
    *fTitle = title ? title : "Container for a MC energy histogram" ;

    //  - we initialize the histogram and the gaus function
    //  - we have to give diferent names for the diferent histograms because
    //    root don't allow us to have diferent histograms with the same name

    strcpy(aux, "fLogEnergy");
    if (idx>0)
        sprintf(aux+10, ";%i", idx);
    fLogEner = new TF1(aux, "gaus", 1., 3.);

    strcpy(aux, "hLogEnergy");
    if (idx>0)
        sprintf(aux+10, ";%i", idx);
    hLogEner = new TH1F(aux, "", 40, 0.5, 4.5);
    hLogEner->SetXTitle("log(E) [GeV]");
    hLogEner->SetYTitle("dN/dE");
}

//-------------------------------------------------------------------------
//
//  Defualt Destructor
//
MHMcEnergy::~MHMcEnergy()
{
    delete hLogEner;
    delete fLogEner;
}

//--------------------------------------------------------------------------
//
//  Fill the histogram with the log10 of the energy for triggered events.
//
void MHMcEnergy::Fill(Float_t log10E, Float_t w)
{
    hLogEner->Fill(log10E, w);
}

// -------------------------------------------------------------------------
//
// Fitting function
//
void MHMcEnergy::Fit(Axis_t xxmin, Axis_t xxmax)
{
    //
    // 0: don't draw the function (it is drawn together with the histogram)
    // Q: quiet mode
    //
    hLogEner->Fit(fLogEner->GetName(), "Q0", "", xxmin, xxmax);
}

// ------------------------------------------------------------------------
// 
// Drawing function. It creates its own canvas.
//
void MHMcEnergy::Draw(Option_t *option)
{
    char text[50];
    sprintf(text, "Energy Threshold = %4.1f +- %4.1f GeV",
            GetThreshold(), GetThresholdErr());

    const Float_t min = hLogEner->GetMinimum();
    const Float_t max = hLogEner->GetMaximum();
    const Float_t sum = min+max;

    TCanvas *c=new TCanvas("Energy Distribution","Energy distribution for triggered events");

    hLogEner->Draw(option);

    TPaveLabel* label = new TPaveLabel(2.2, 0.75*sum, 4.4, 0.90*sum, text);
    label->SetFillColor(10);
    label->SetTextSize(0.3);
    label->Draw();

    c->Modified();
    c->Update();
}

// --------------------------------------------------------------------------
//
// Set the number of bins in the histogran.
//
void MHMcEnergy::SetNumBins(Int_t nbins)
{
    hLogEner->SetBins(nbins, 0.5, 4.5);
}
// --------------------------------------------------------------------------
//
// Write the threshold and its error in the standard output
//
void MHMcEnergy::Print(Option_t*)
{
    cout << "Threshold: " << GetThreshold() << " +- " << GetThresholdErr() << endl;
}

// -------------------------------------------------------------------------
//
//  Return the threshold
//
Float_t MHMcEnergy::GetThreshold() const
{
    const Float_t p1 = fLogEner->GetParameter(1);

    return pow(10, p1);
}

// -------------------------------------------------------------------------
//
// Return the error of the threshold.
//
Float_t MHMcEnergy::GetThresholdErr() const
{
    const Float_t lg10  = log(10);
    const Float_t p1    = fLogEner->GetParameter(1);
    const Float_t p1err = fLogEner->GetParError(1);

    // The error has into accuont the error in the fit
    return pow(10, p1) * p1err * lg10;
}

// -------------------------------------------------------------------------
//
// Return the peak of the fitted gaussan function.
//
Float_t MHMcEnergy::GetGaussPeak() const
{
    return fLogEner->GetParameter(1);
}

// -------------------------------------------------------------------------
//
// Return the sigma of the fitted gaussan function.
//
Float_t MHMcEnergy::GetGaussSigma() const
{
    return fLogEner->GetParameter(2);
}