source: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc@ 4017

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

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHEffOnTimeTheta                                                        //
//                                                                          //
//  calculates the effective on time for each bin in Theta                  //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHEffOnTimeTheta.h"

#include <TStyle.h>

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

#include "MTime.h"

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MHEffOnTimeTheta);

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title of the histograms. 
//
MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title)
    : fHEffOn()
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHEffOnTimeTheta";
    fTitle = title ? title : "1-D histogram of Eff On Time";

    // effective on time versus theta
    fHEffOn.SetName("EffOn");
    fHEffOn.SetTitle("Effective On Time vs. Theta");

    fHEffOn.SetDirectory(NULL);

    fHEffOn.SetXTitle("\\Theta [\\circ]");
    fHEffOn.SetYTitle("t_{eff} [s]");

    // chi2/NDF versus theta
    fHChi2.SetName("Chi2/NDF");
    fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Theta");

    fHChi2.SetDirectory(NULL);

    fHChi2.SetXTitle("\\Theta [\\circ]");
    fHChi2.SetYTitle("chi^{2}/NDF");

    // lambda versus theta
    fHLambda.SetName("lambda");
    fHLambda.SetTitle("lambda of OnTimeFit vs. Theta");

    fHLambda.SetDirectory(NULL);

    fHLambda.SetXTitle("\\Theta [\\circ]");
    fHLambda.SetYTitle("\\lambda [Hz]");

    // N0del versus theta
    fHN0del.SetName("N0del");
    fHN0del.SetTitle("N0del of OnTimeFit vs. Theta");

    fHN0del.SetDirectory(NULL);

    fHN0del.SetXTitle("\\Theta [\\circ]");
    fHN0del.SetYTitle("N0del");
}

// -----------------------------------------------------------------------
//
// Calculate the effective on time by fitting the distribution of
// time differences
//
void MHEffOnTimeTheta::Calc(TH2D *hist)
{
    // nbins = number of Theta bins
    const Int_t nbins = hist->GetNbinsY();

    for (int i=1; i<=nbins; i++)
    {
        //        TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
        TH1D &h = *hist->ProjectionX("Calc-theta", i, i, "E");

        // Nmdel = Nm * binwidth,  with Nm = number of observed events
        const Double_t Nmdel = h.Integral("width");
        const Double_t Nm    = h.Integral();
        //        Double_t mean  = h->GetMean();

        //...................................................
        // determine range (yq[0], yq[1]) of time differences 
        // where fit should be performed;
        // require a fraction >=xq[0] of all entries to lie below yq[0]
        //     and a fraction <=xq[1] of all entries to lie below yq[1];  
        // within the range (yq[0], yq[1]) there must be no empty bin;
        // choose pedestrian approach as long as GetQuantiles is not available

        Double_t xq[2] = { 0.15, 0.95 };
        Double_t yq[2];

        // GetQuantiles doesn't seem to be available in root 3.01/06
        // h->GetQuantiles(2,yq,xq);

        const Double_t sumtot = h.Integral();
        const Int_t    jbins  = h.GetNbinsX();

        if (sumtot > 0.0)
        {
            Double_t sum1 = 0.0;
            yq[0] = h.GetBinLowEdge(jbins+1);
            for (int j=1; j<=jbins; j++)
            {
                if (sum1 >= xq[0]*sumtot)
                {
                    yq[0] = h.GetBinLowEdge(j);
                    break;
                }
                sum1 += h.GetBinContent(j);
            }

            Double_t sum2 = 0.0;
            yq[1] = h.GetBinLowEdge(jbins+1);
            for (int j=1; j<=jbins; j++)
            {
                Double_t content = h.GetBinContent(j);
                sum2 += content;
                if (sum2 >= xq[1]*sumtot || content == 0.0)
                {
                    yq[1] = h.GetBinLowEdge(j);
                    break;
                }
            }

            //...................................................

            // parameter 0 = lambda
            // parameter 1 = N0*del        with N0 = ideal number of events
            //                             and del = bin width of time difference
            TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", yq[0], yq[1]);

            func.SetParameter(0, 100); // [Hz]
            func.SetParameter(1, Nmdel);

            func.SetParLimits(0, 0, 1000);    // [Hz]
            func.SetParLimits(1, 0, 10*Nmdel);

            func.SetParName(0, "lambda");
            func.SetParName(1, "Nmdel");

            // options : 0  do not plot the function
            //           I  use integral of function in bin rather than value at bin center
            //           R  use the range specified in the function range
            //           Q  quiet mode
            h.Fit("Poisson", "0IRQ");

            // gPad->SetLogy();
            // gStyle->SetOptFit(1011);
            // h->GetXaxis()->SetTitle("time difference [s]");
            // h->GetYaxis()->SetTitle("Counts");
            // h->DrawCopy();

            // func.SetRange(yq[0], yq[1]); // Range of Drawing
            // func.DrawCopy("same");

            const Double_t lambda = func.GetParameter(0);
            const Double_t N0del  = func.GetParameter(1);
            const Double_t chi2   = func.GetChisquare();
            const Int_t    NDF    = func.GetNDF();

            // was fit successful ?
            if (NDF>0 && chi2<2.5*NDF)
            {
                // the effective on time is Nm/lambda
                fHEffOn.SetBinContent(i, Nm/lambda);

                // plot chi2/NDF of fit
                fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0);

                // lambda of fit
                fHLambda.SetBinContent(i, lambda);

                // N0del of fit
                fHN0del.SetBinContent(i, N0del);

                delete &h;
                continue;
            }
        }

        fHEffOn.SetBinContent (i, 1.e-20);
        fHChi2.SetBinContent  (i, 1.e-20);
        fHLambda.SetBinContent(i, 1.e-20);
        fHN0del.SetBinContent (i, 1.e-20);

        delete &h;
    }
}

// -------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histograms
//
Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist)
{
   const MBinning* bins = (MBinning*)plist->FindObject("BinningTheta");
   if (!bins)
   {
       *fLog << err << dbginf << "BinningTheta [MBinning] not found... aborting." << endl;
       return kFALSE;
   }

   SetBinning(&fHEffOn,  bins);
   SetBinning(&fHChi2,   bins);
   SetBinning(&fHLambda, bins);
   SetBinning(&fHN0del,  bins);

   fHEffOn.Sumw2();
   fHChi2.Sumw2();
   fHLambda.Sumw2();
   fHN0del.Sumw2();

   return kTRUE;
}

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHEffOnTimeTheta::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    pad->Divide(2,2);

    pad->cd(1);
    gPad->SetBorderMode(0);
    fHEffOn.Draw(opt);

    pad->cd(2);
    gPad->SetBorderMode(0);
    fHChi2.Draw(opt);

    pad->cd(3);
    gPad->SetBorderMode(0);
    fHLambda.Draw(opt);

    pad->cd(4);
    gPad->SetBorderMode(0);
    fHN0del.Draw(opt);

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