source: trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc@ 1610

/* ======================================================================== *\
!
! *
! * 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/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHMcCollectionArea                                                      //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHMcCollectionArea.h" 

#include <TH2.h>
#include <TCanvas.h>

#include "MH.h"
#include "MHMcEfficiency.h"
#include "MHMcEnergyImpact.h"

ClassImp(MHMcCollectionArea);

// --------------------------------------------------------------------------
//
//  Creates the three necessary histograms:
//   - selected showers (input)
//   - all showers (input)
//   - collection area (result)
//
MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title)
{ 
    //   initialize the histogram for the distribution r vs E
    //
    //   we set the energy range from 1 Gev to 10000 GeV (in log 5 orders
    //   of magnitude) and for each order we take 10 subdivision --> 50 xbins
    //
    //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins

  
    fName  = name  ? name  : "MHMcCollectionArea";
    fTitle = title ? title : "Collection Area vs. Energy";

    const Int_t  nbins = 50;
    const Float_t maxx = 5;
    const Float_t maxy = 500;

    Float_t *binsx = new Float_t[nbins+1];
    for (int i=0; i<nbins+1; i++)
        binsx[i] = pow(10, maxx*i/nbins);

    Float_t *binsy = new Float_t[nbins+1];
    for (int i=0; i<nbins+1; i++)
        binsy[i] = maxy*i/nbins;

    fHistAll = new TH2D("AllEvents", "All showers - Radius vs Energy distribution",
                        nbins, binsx, nbins, binsy);
    fHistSel = new TH2D("SelectedEvents", "Selected showers - Radius vs Energy distribution",
                        nbins, binsx, nbins, binsy);

    delete binsx;
    delete binsy;

    fHistCol = new TH1D;
    fHistCol->SetName(fName);
    fHistCol->SetTitle(fTitle);

    fHistAll->SetDirectory(NULL);
    fHistSel->SetDirectory(NULL);
    fHistCol->SetDirectory(NULL);

    fHistAll->SetXTitle("E [GeV]");
    fHistAll->SetYTitle("r [m]");
    fHistAll->SetZTitle("N");

    fHistSel->SetXTitle("E [GeV]");
    fHistSel->SetYTitle("r [m]");
    fHistSel->SetYTitle("N");

    fHistCol->SetXTitle("E [GeV]");
    fHistCol->SetYTitle("A [m^{2}]");
}

// --------------------------------------------------------------------------
//
// Delete the three histograms
//
MHMcCollectionArea::~MHMcCollectionArea()
{
    delete fHistAll;
    delete fHistSel;
    delete fHistCol;
}

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains all showers
//
void MHMcCollectionArea::FillAll(Float_t energy, Float_t radius)
{
    fHistAll->Fill(energy, radius);
}

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains the selected showers
//
void MHMcCollectionArea::FillSel(Float_t energy, Float_t radius)
{
    fHistSel->Fill(energy, radius);
}

// --------------------------------------------------------------------------
//
// Draw the histogram with all showers
//
void MHMcCollectionArea::DrawAll(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistAll);

    fHistAll->Draw(option);

    gPad->SetLogx();

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

// --------------------------------------------------------------------------
//
// Draw the histogram with the selected showers only.
//
void MHMcCollectionArea::DrawSel(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistSel);

    fHistSel->Draw(option);

    gPad->SetLogx();

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

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the histogram into it.
// Be careful: The histogram belongs to this object and won't get deleted
// together with the canvas.
//
TObject *MHMcCollectionArea::DrawClone(Option_t* option) const
{
    TCanvas *c = MH::MakeDefCanvas(fHistCol);

    //
    // This is necessary to get the expected bahviour of DrawClone
    //
    gROOT->SetSelectedPad(NULL);

    fHistCol->DrawCopy(option);

    gPad->SetLogx();

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

    return c;
}

void MHMcCollectionArea::Draw(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistCol);

    fHistCol->Draw(option);

    gPad->SetLogx();

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

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::CalcEfficiency()
{
    Calc(*fHistSel, *fHistAll);
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall)
{
    Calc((TH2D&)*mcsel.GetHist(), (TH2D&)*mcall.GetHist());
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::Calc(TH2D &hsel, TH2D &hall)
{
    MH::SetBinning(fHistCol, hsel.GetXaxis());

    fHistCol->Sumw2();

    TH1D &psel = *hsel.ProjectionX();
    TH1D &pall = *hall.ProjectionX();

    const Double_t max = psel.GetYaxis()->GetXmax();

    fHistCol->Divide(&psel, &pall, TMath::Pi()*max*max, 1);
    fHistCol->SetEntries(hsel.GetEntries());

    delete &psel;
    delete &pall;

    SetReadyToSave();
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t angle)
{
    // Here we estimate the total number of showers in each energy bin
    // known the energy range and spectral index of the generated showers
    // (set now by hand since the info is not available in the header!)

    TH1D &histsel = *fHistSel->ProjectionX();

    TH1D histall;

    TAxis &xaxis = *histsel.GetXaxis();

    MH::SetBinning(&histall, &xaxis);
    MH::SetBinning(fHistCol, &xaxis);

    const Float_t emin  = 10.;
    const Float_t emax  = 30000.;  // Energies in GeV.
    const Float_t index = 2.6;     // Differential spectral Index.

    const Float_t expo = 1.-index;

    const Float_t k = (Float_t)numevts / (pow(emax,expo) - pow(emin,expo));

    const Int_t nbinx = xaxis.GetNbins();

    for (Int_t i=1; i<=nbinx; i++)
    {
        const Float_t e1 = histall.GetBinLowEdge(i);
        const Float_t e2 = histall.GetBinLowEdge(i+1);

        if (e1 < emin || e2 > emax)
            continue;

        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));

        histall.SetBinContent(i, events);
        histall.SetBinError(i, sqrt(events));

    }

    // -----------------------------------------------------------

    // Impact parameter range:
    const Float_t r1 = 0;
    const Float_t r2 = 400;

    const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);

    angle *= TMath::Pi()/180;

    for (Int_t ix=1; ix<=nbinx; ix++)
    {
        const Float_t Na = histall.GetBinContent(ix);

        if (Na <= 0)
            continue;

        const Float_t Ns = histsel.GetBinContent(ix);

        // Since Na is an estimate of the total number of showers generated
        // in the energy bin, it may happen that Ns (triggered showers) is
        // larger than Na. In that case, the bin is skipped:

        if (Na < Ns)
            continue;

        const Double_t eff = Ns/Na;

        const Double_t err = sqrt((1.-eff)*Ns)/Na;

        const Float_t area = dr * cos(angle);

        fHistCol->SetBinContent(ix, eff*area);
        fHistCol->SetBinError(ix, err*area);
    }

    delete &histsel;

    SetReadyToSave();
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::Calc(const MHMcEfficiency &heff)
{
    //
    //  now calculate the Collection area for different
    //  energies
    //
    TH2D &h = (TH2D&)*heff.GetHist();

    MH::SetBinning(fHistCol, h.GetXaxis());

    const Int_t nbinx = h.GetXaxis()->GetNbins();
    const Int_t nbiny = h.GetYaxis()->GetNbins();

    for (Int_t ix=1; ix<=nbinx; ix++)
    {
        Double_t errA = 0;
        Double_t colA = 0;

        for (Int_t iy=1; iy<=nbiny; iy++)
        {
            TAxis *yaxis = h.GetYaxis();

            const Double_t r1  = yaxis->GetBinLowEdge(iy);
            const Double_t r2  = yaxis->GetBinLowEdge(iy+1);

            const Double_t A   = TMath::Pi() * (r2*r2 - r1*r1);

            const Double_t eff = h.GetCellContent(ix, iy);
            const Double_t err = h.GetCellError(ix, iy);

            colA += eff*A;
            errA += A*A * err*err;
        }

        fHistCol->SetBinContent(ix, colA);
        fHistCol->SetBinError(ix, sqrt(errA));
    }

    SetReadyToSave();
}