source: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc@ 6938

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

//////////////////////////////////////////////////////////////////////////////
//
//  MHCollectionArea
//
//////////////////////////////////////////////////////////////////////////////
#include "MHCollectionArea.h" 

#include <TCanvas.h>

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

#include "MH.h"
#include "MBinning.h"

#include "MMcEvt.hxx"
#include "MMcRunHeader.hxx"
#include "MMcCorsikaRunHeader.h"

#include "MParList.h"
#include "MParameters.h"

ClassImp(MHCollectionArea);

using namespace std;

// --------------------------------------------------------------------------
//
//  Creates the three necessary histograms:
//   - selected showers (input)
//   - all showers (input)
//   - collection area (result)
//
MHCollectionArea::MHCollectionArea(const char *name, const char *title)
  : fMcEvt(0), fEnergy(0), fMcAreaRadius(300.), fIsExtern(kFALSE)
{ 
    //   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  : "MHCollectionArea";
    fTitle = title ? title : "Collection Area vs. Energy/Theta";

    fHistSel.SetName("SelEvts");
    fHistSel.SetTitle("Number of Events after cuts");
    fHistSel.SetXTitle("\\Theta [deg]");
    fHistSel.SetYTitle("E [GeV]");
    fHistSel.SetDirectory(NULL);
    fHistSel.UseCurrentStyle();

    fHistAll.SetName("AllEvts");
    fHistAll.SetTitle("Number of events produced");
    fHistAll.SetXTitle("\\Theta [deg]");
    fHistAll.SetYTitle("E [GeV]");
    fHistAll.SetDirectory(NULL);
    fHistAll.UseCurrentStyle();

    fHEnergy.SetName("CollEnergy");
    fHEnergy.SetTitle("Collection Area vs Energy");
    fHEnergy.SetXTitle("E [GeV]");
    fHEnergy.SetYTitle("A [m^{2}]");
    fHEnergy.SetDirectory(NULL);
    fHEnergy.UseCurrentStyle();

    MBinning binsa, binse, binst;
    binse.SetEdgesLog(15, 10, 1000000);
    binst.SetEdgesCos(50, 0, 60);

    binse.Apply(fHEnergy);

    MH::SetBinning(&fHistSel, &binst, &binse);
    MH::SetBinning(&fHistAll, &binst, &binse);
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
//  flag
//
void MHCollectionArea::CalcEfficiency()
{
    TH1D *hsel = fHistSel.ProjectionY();
    TH1D *hall = fHistAll.ProjectionY();

    const Int_t nbinx = hsel->GetNbinsX();

    // -----------------------------------------------------------
    //
    // 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 = fMcAreaRadius;

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

    for (Int_t ix=1; ix<=nbinx; ix++)
    {
        const Float_t Na = hall->GetBinContent(ix);

        if (Na <= 0)
            continue;

        const Float_t Ns = hsel->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 = TMath::Sqrt((1.-eff)*Ns)/Na;
        
        const Float_t col_area  =  eff * total_area;
        const Float_t col_area_error =  efferr * total_area;

        fHEnergy.SetBinContent(ix, col_area);
        fHEnergy.SetBinError(ix, col_area_error);
    }

    delete hsel;
    delete hall;
}

Bool_t MHCollectionArea::SetupFill(const MParList *pl)
{
    fHistSel.Reset();
    if (!fIsExtern)
        fHistAll.Reset();

    fMcEvt = (MMcEvt*)pl->FindObject("MMcEvt");
    if (!fMcEvt)
    {
        *fLog << err << "MMcEvt not found... abort." << endl;
        return kFALSE;
    }

    fEnergy = (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
    if (!fEnergy)
    {
        *fLog << err << "MEnergyEst [MParameterD] not found... abort." << endl;
        return kFALSE;
    }

    MBinning binst, binse;
    binst.SetEdges(fHistAll, 'x');
    binse.SetEdges(fHistAll, 'y');

    //if (!fIsExtern)
    {
        MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
        if (bins)
            binst.SetEdges(*bins);

        bins = (MBinning*)pl->FindObject("BinningEnergy", "MBinning");
        if (bins)
            binse.SetEdges(*bins);
    }

    binse.Apply(fHEnergy);

    MH::SetBinning(&fHistSel, &binst, &binse);
    MH::SetBinning(&fHistAll, &binst, &binse);

    return kTRUE;
}

Bool_t MHCollectionArea::ReInit(MParList *plist)
{
    if (fIsExtern)
        return kTRUE;

    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    if (!runheader)
    {
        *fLog << err << "MMcRunHeader not found... abort." << endl;
        return kFALSE;
    }

    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;

 /*
    if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
    {
        *fLog << warn << "Warning - Read files have different TelesTheta (";
        *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl;
    }
    fTheta = runheader->GetTelesTheta();
  */
    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
    fCorsikaVersion = runheader->GetCorsikaVersion();

    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;

    MMcCorsikaRunHeader *crh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
    if (!crh)
    {
        *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
        return kFALSE;
    }

    //
    // Calculate approximately the original number of events in each
    // energy bin:
    //
    const Float_t emin = crh->GetELowLim();
    const Float_t emax = crh->GetEUppLim();
    const Float_t expo = 1 + crh->GetSlopeSpec();
    const Float_t k = runheader->GetNumSimulatedShowers() /
        (pow(emax,expo) - pow(emin,expo));

    const Int_t nbiny = fHistAll.GetNbinsY();

    TAxis &axe = *fHistAll.GetYaxis();
    for (Int_t i = 1; i <= nbiny; i++)
    {
        const Float_t e1 = axe.GetBinLowEdge(i);
        const Float_t e2 = axe.GetBinLowEdge(i+1);

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

        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
        //
        // We fill the i-th energy bin, with the total number of events
        // Second argument of Fill would be impact parameter of each
        // event, but we don't really need it for the collection area,
        // so we just put a dummy value (1.)
        //

        const Float_t energy = (e1+e2)/2.;
        fHistAll.Fill(20, energy, events);
        // you have  MMcRunHeader.fShowerThetaMin and MMcRunHeader.fShowerThetaMax
    }

    return kTRUE;
}

void MHCollectionArea::Paint(Option_t *option)
{
    TVirtualPad *pad = gPad;

    TH1 *h=0;

    pad->cd(1);
    //if (gPad->FindObject("ProjAllX"))
    //    fHistAll.ProjectionX("ProjAllX", -1, 9999, "E");
    if (gPad->FindObject("ProjSelX"))
        fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");

    pad->cd(2);
    if (gPad->FindObject("ProjAllY"))
        h=fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
    if (gPad->FindObject("ProjSelY"))
        fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");

    if (h->GetMaximum()>0)
    {
        gPad->SetLogx();
        gPad->SetLogy();
    }
    pad->cd(4);
    CalcEfficiency();
    if (fHEnergy.GetMaximum()>0)
    {
        gPad->SetLogx();
        gPad->SetLogy();
    }
}

void MHCollectionArea::Draw(Option_t *option)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);

    // Do the projection before painting the histograms into
    // the individual pads
    AppendPad();

    pad->SetBorderMode(0);
    pad->Divide(2,2);

    TH1D *h=0;

    if (fHistSel.GetNbinsX()>1)
    {
        pad->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        /*
        h = fHistAll.ProjectionX("ProjAllX", -1, 9999, "E");
        h->SetXTitle("\\Theta [\\circ]");
        h->SetDirectory(NULL);
        h->SetLineColor(kGreen);
        h->SetBit(kCanDelete);
        h->Draw();
        */
        h = fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");
        h->SetXTitle("\\Theta [\\circ]");
        h->SetDirectory(NULL);
        h->SetLineColor(kRed);
        h->SetBit(kCanDelete);
        h->Draw("hist"/*"same"*/);
    }
    else
        delete pad->GetPad(1);

    if (fHistSel.GetNbinsY()>1)
    {
        pad->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();

        h = fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
        h->SetDirectory(NULL);
        h->SetLineColor(kGreen);
        h->SetXTitle("E [GeV]");
        h->SetBit(kCanDelete);
        h->Draw();

        h = fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");
        h->SetDirectory(NULL);
        h->SetLineColor(kRed);
        h->SetBit(kCanDelete);
        h->Draw("same");
    }
    else
        delete pad->GetPad(2);

    if (fHEnergy.GetNbinsX()>1)
    {
        pad->cd(4);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        fHEnergy.Draw();
    }
    else
        delete pad->GetPad(4);
/*
    if (fHistAll.GetNbinsY()>1)
    {
        pad->cd(4);
    }
    else
        delete pad->GetPad(4);
  */
}

Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
{
    const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/;
    const Double_t theta  = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();

    //*fLog << energy << " " << theta << endl;

    fHistSel.Fill(theta, energy);

    return kTRUE;
}