source: trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc@ 1823

/* ======================================================================== *\
!
! *
! * 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): A. Moralejo 3/2003  moralejo@pd.infn.it
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHMcCT1CollectionArea                                                   //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHMcCT1CollectionArea.h" 

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

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

ClassImp(MHMcCT1CollectionArea);

// --------------------------------------------------------------------------
//
//  Creates the three necessary histograms:
//   - selected showers (input)
//   - all showers (input)
//   - collection area (result)
//
MHMcCT1CollectionArea::MHMcCT1CollectionArea(const char *name, const char *title)
{ 
    //   initialize the histogram for the distribution r vs E
    //
    //   we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders
    //   of magnitude) and for each order we take 10 subdivision --> 45 xbins
    //   we set the theta range from 12.5 to 48 deg, with 6 bins.
    //
    fName  = name  ? name  : "MHMcCT1CollectionArea";
    fTitle = title ? title : "Collection Area vs. Energy";

    MBinning binsx;
    MBinning binsy;

    binsx.SetEdgesLog(45, 10, 30000);
    const Double_t yedge[7] = {12.5, 17.5, 23.5, 29.5, 35.5, 42., 48.};
    const TArrayD yed(7,yedge);
    binsy.SetEdges(yed);

    fHistAll = new TH2D;
    fHistSel = new TH2D;
    fHistCol = new TH2D;

    MH::SetBinning(fHistAll, &binsx, &binsy);
    MH::SetBinning(fHistSel, &binsx, &binsy);
    MH::SetBinning(fHistCol, &binsx, &binsy);


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

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

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

    fHistAll->SetXTitle("E [GeV]");
    fHistAll->SetYTitle("theta [deg]");
    fHistAll->SetZTitle("N");

    fHistSel->SetXTitle("E [GeV]");
    fHistSel->SetYTitle("theta [deg]");
    fHistSel->SetYTitle("N");

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

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

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains the selected showers
//
void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t theta)
{
    fHistSel->Fill(energy, theta);
}

// --------------------------------------------------------------------------
//
// Draw the histogram with all showers
//
void MHMcCT1CollectionArea::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 MHMcCT1CollectionArea::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 *MHMcCT1CollectionArea::DrawClone(Option_t* option) const
{
    TCanvas *c = MH::MakeDefCanvas(fHistCol);

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

    fHistCol->DrawCopy(option);

    gPad->SetLogx();

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

    return c;
}

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

    fHistCol->Draw(option);

    gPad->SetLogx();

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

//
//  Calculate the Efficiency (collection area) for the CT1 MC sample
//  and set the 'ReadyToSave' flag
//
void MHMcCT1CollectionArea::CalcEfficiency()
{
  //
  // Here we estimate the total number of showers in each energy bin
  // from the known the energy range and spectral index of the generated 
  // showers. This procedure is intended for the CT1 MC files. The total 
  // number of generated events, collection area, spectral index etc will be 
  // set here by hand, so make sure that the MC sample you are using is the 
  // right one (check all these quantities in your files and compare with
  // is written below. In some theta bins, there are two different 
  // productions, with different energy limits but with the same spectral 
  // slope. We account for this when calculating the original number of 
  // events in each energy bin.
  //

  for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
    {
      // This theta is not exactly the one of the MC events, just about 
      // the same:
      Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);

      Float_t emin1, emax1, emin2, emax2;
      Float_t index, expo, k1, k2;
      Float_t numevts1, numevts2;
      Float_t r1, r2;        // Impact parameter range (on ground).

      emin1 = 0;  emax1 = 0; emin2 = 0;  emax2 = 0;
      expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.;
      numevts1 = 0; numevts2 = 0;

      if (theta > 14 && theta < 16)   // 15 deg
        {
          r1 = 0.;
          r2 = 250.;     //meters
          emin1 = 300.;
          emax1 = 400.;  // Energies in GeV.
          emin2 = 400.;
          emax2 = 30000.;
          numevts1 = 4000.;
          numevts2 = 25740.;
        }
      else if (theta > 20 && theta < 21)  // 20.5 deg 
        {
          r1 = 0.;
          r2 = 263.;     //meters
          emin1 = 300.;
          emax1 = 400.;  // Energies in GeV.
          emin2 = 400.;
          emax2 = 30000.;
          numevts1 = 6611.;
          numevts2 = 24448.;
        }
      else if (theta > 26 && theta < 27)  // 26.5 degrees
        {
          r1 = 0.;
          r2 = 290.;     //meters
          emin1 = 300.;
          emax1 = 400.;  // Energies in GeV.
          emax2 = emax1;          emin2 = 400.;
          emax2 = 30000.;
          numevts1 = 4000.;
          numevts2 = 26316.;
        }
      else if (theta > 32 && theta < 33)  // 32.5 degrees
        {
          r1 = 0.;
          r2 = 350.;     //meters
          emin1 = 300.;
          emax1 = 30000.;  // Energies in GeV.
          emax2 = emax1;
          numevts1 = 33646.;
        }
      else if (theta > 38 && theta < 39)  // 38.75 degrees
        {
          r1 = 0.;
          r2 = 380.;     //meters
          emin1 = 300.;
          emax1 = 30000.;  // Energies in GeV.
          emax2 = emax1;
          numevts1 = 38415.;
        }
      else if (theta > 44 && theta < 46)  // 45 degrees
        {
          r1 = 0.;
          r2 = 565.;     //meters
          emin1 = 300.;
          emax1 = 50000.;  // Energies in GeV.
          emax2 = emax1;
          numevts1 = 30197.;
        }

      index = 1.5;     // Differential spectral Index.
      expo = 1.-index;
      k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo));
      k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo));

      for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
        {
          const Float_t e1 = fHistAll->GetXaxis()->GetBinLowEdge(i);
          const Float_t e2 = fHistAll->GetXaxis()->GetBinLowEdge(i+1);

          if (e1 < emin1 || e2 > emax2)
            continue;

          Float_t events;

          if (e2 <= emax1)
            events = k1 * (pow(e2, expo) - pow(e1, expo));
          else if (e1 >= emin2)
            events = k2 * (pow(e2, expo) - pow(e1, expo));
          else
            events = 
              k1 * (pow(emax1, expo) - pow(e1, expo))+
              k2 * (pow(e2, expo) - pow(emin2, expo));

          fHistAll->SetBinContent(i, thetabin, events);
          fHistAll->SetBinError(i, thetabin, sqrt(events));
        }

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

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

      for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
        {
          const Float_t Na = fHistAll->GetBinContent(ix,thetabin);

          if (Na <= 0)
            continue;

          const Float_t Ns = fHistSel->GetBinContent(ix,thetabin);

          // 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(theta*TMath::Pi()/180.);

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

    SetReadyToSave();
}