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

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

#include "MHMcEnergy.h" 

#include <iostream.h>

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

ClassImp(MHMcEnergy)

MHMcEnergy::MHMcEnergy(const UInt_t idx, const char *name, const char *title)
{ 
    //
    //   default constructor
    //
    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, "", 100, 0.5, 4.5);
    hLogEner->SetXTitle("log(E) [GeV]");
    hLogEner->SetYTitle("dN/dE");
    //hLogEner->SetBins(60);
}

MHMcEnergy::~MHMcEnergy()
{
    delete hLogEner;
    delete fLogEner;
}

void MHMcEnergy::Fill(Float_t log10E, Float_t w)
{
    hLogEner->Fill(log10E, w);
}

void MHMcEnergy::Fit(Axis_t xxmin, Axis_t xxmax)
{
    //
    // 0: don't draw the function (it is drawn together with the histogram)
    // +: add these function to the list of fits. Don't delete the last fit.
    //
    // FIXME: R means: use the range specified in the function (xxmin, xxmax are ignored!)
    //        Q means: quiet (why?)
    //
    //
    hLogEner->Fit(fLogEner->GetName(), "Q0+", "", xxmin, xxmax);
}

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

void MHMcEnergy::Print(Option_t*)
{
    cout << "Threshold: " << GetThreshold() << " +- " << GetThresholdErr() << endl;
}

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

    return pow(10, p1);
}

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

    return pow(10, p1) * p1err * lg10;
}

Float_t MHMcEnergy::GetGaussPeak() const
{
    return fLogEner->GetParameter(1);
}

Float_t MHMcEnergy::GetGaussSigma() const
{
    return fLogEner->GetParameter(2);
}