source: trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc@ 1978

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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  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, Int_t nbins, Axis_t minEnergy, Axis_t maxEnergy)
{ 
  //
  //   nbins, minEnergy, maxEnergy defaults:
  //   we set the energy range from 100 Gev to 30000 GeV (in log, 3.5 orders
  //   of magnitude) and for each order we take 10 subdivisions --> 35 xbins
  //   we set the theta range from 12.5 to 48 deg, with 6 bins (the latter
  //   choice has been done to make the bin centers as close as possible to 
  //   the actual zenith angles in the CT1 MC sample).
  //

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

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

  SetBins(nbins, minEnergy, maxEnergy);

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

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

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

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

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

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

// --------------------------------------------------------------------------
//
// Set binning of histograms. Binning of energy axis can be changed, that
//   of theta axis is always the same (to match the CT1 MC sample theta dist):
//
void MHMcCT1CollectionArea::SetBins(Int_t nbins, Axis_t minEnergy, Axis_t maxEnergy)
{
  MBinning binsx;
  binsx.SetEdges(nbins, minEnergy, maxEnergy);

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

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


// --------------------------------------------------------------------------
//
// 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(log10(energy), theta);
}

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

  fHistAll->Draw(option);

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

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

  return c;
}

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

  fHistCol->Draw(option);

  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 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
          const Float_t e2 = pow(10.,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)
            {
              //
              // If energy is large, this case means that no or very few events
              // were generated at this energy bin. In this case we assign it 
              // the effective area of the bin below it in energy. If energy is
              // below 1E4, it means that no events triggered -> eff area = 0
              //

              if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.)
                {
                  fHistCol->SetBinContent(ix, thetabin, fHistCol->GetBinContent(ix-1, thetabin));
                  fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin));
                }
              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();
}