source: trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc@ 2160

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

#include "MLog.h"
#include "MLogManip.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 2 Gev to 20000 GeV (in log 4 orders
    //   of magnitude) and for each order we take 25 subdivision --> 100 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";

    MBinning binsx;
    MBinning binsy;
    binsx.SetEdgesLog(100, 2., 20000);
    binsy.SetEdges   (50,  0,   500);

    fHistAll = new TH2D;
    fHistSel = new TH2D;
    fHistCol = new TH1D;
 
    MH::SetBinning(fHistAll, &binsx, &binsy);
    MH::SetBinning(fHistSel, &binsx, &binsy);

    fHistCol->SetName(fName);
    fHistAll->SetName("AllEvents");
    fHistSel->SetName("SelectedEvents");

    fHistCol->SetTitle(fTitle);
    fHistAll->SetTitle("All showers - Radius vs Energy distribution");
    fHistSel->SetTitle("Selected showers - Radius vs Energy distribution");

    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(Double_t energy, Double_t radius)
{
    fHistAll->Fill(energy, radius);
}

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains the selected showers
//
void MHMcCollectionArea::FillSel(Double_t energy, Double_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();
}

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

    fHistCol->Draw(option);

    pad->SetLogx();

    pad->Modified();
    pad->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();

    TH1D &pally = *hall.ProjectionY();

    Double_t max = 0.;
    for (Int_t i = hall.GetNbinsY(); i > 0; i--)
      if (pally.GetBinContent(i) > 0)
        {
          max = pally.GetBinLowEdge(i+1);
          break;
        }

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

    delete &psel;
    delete &pall;
    delete &pally;

    SetReadyToSave();
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t emin, Float_t emax, Float_t index)
{
    // 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 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:  TO BE FIXED! Impact parameter shoud be
    // read from run header, but it is not yet in!!
    //
    const Float_t r1 = 0;
    const Float_t r2 = 300;

    *fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of 300 meters for the MC events. Check that this is the true one!" <<endl<<endl;
    const Float_t area = TMath::Pi() * (r2*r2 - r1*r1);

    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 efferr = sqrt((1.-eff)*Ns)/Na;

        fHistCol->SetBinContent(ix, eff*area);
        fHistCol->SetBinError(ix, efferr*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 efferr = h.GetCellError(ix, iy);

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

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

    SetReadyToSave();
}