source: trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc@ 2894

/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek 5/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHEffOnTime                                                             //
//                                                                          //
//  calculates the effective on time for each bin in the variable Var       //
//             (Var may be time or Theta)                                   //
//                                                                          //
//  Important : The sample of events used for this should be as close       //
//              to the sample of recorded events as possible.               //
//              This means that NO additional cuts (be it quality cuts or   //
//              gamma/hadron cuts) should be applied (see MAGIC-TDAS 02-02).//
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MHEffOnTime.h"

#include <TStyle.h>

#include <TMinuit.h>
#include <TFitter.h>

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

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

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

#include "MHTimeDiffTheta.h"

ClassImp(MHEffOnTime);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. It sets the name of the variable ("Time" or "Theta")
//                      and the units of the variable ("[s]" or "[\\circ])")
//
MHEffOnTime::MHEffOnTime(const char *varname, const char *unit)
    : fHEffOn(), fHProb(), fHLambda(), fHRdead()
{
    fVarname  = varname;
    fUnit     = unit;

    TString strg = fVarname + " " + fUnit;

    //
    //   set the name and title of this object
    //
    fName  = TString("MHEffOnTime-")+fVarname;
    fTitle = "1-D histogram of Eff On Time";

    // effective on time versus Var
    fHEffOn.SetName("EffOn");
    fHEffOn.SetTitle(TString("T_{on, eff} vs. ")+fVarname);

    fHEffOn.SetDirectory(NULL);

    fHEffOn.SetXTitle(strg);
    fHEffOn.SetYTitle("T_{on, eff} [s]");

    // chi2-probability versus Var
    fHProb.SetName("Chi2-prob");
    TString strg3("\\chi^{2}-prob of OnTimeFit vs. ");
    strg3 += fVarname;
    fHProb.SetTitle(strg3);

    fHProb.SetDirectory(NULL);

    fHProb.SetXTitle(strg);
    fHProb.SetYTitle("\\chi^{2}-probability");

    // lambda versus Var
    fHLambda.SetName("lambda");
    fHLambda.SetTitle(TString("\\lambda from OnTimeFit vs. ")+fVarname);

    fHLambda.SetDirectory(NULL);

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

    // Rdead versus Var
    // Rdead is the fraction from real time lost by the dead time
    fHRdead.SetName("Rdead");

    fHRdead.SetTitle(TString("Rdead of OnTimeFit vs. ")+fVarname);

    fHRdead.SetDirectory(NULL);

    fHRdead.SetXTitle(strg);
    fHRdead.SetYTitle("Rdead");

}

// --------------------------------------------------------------------
//
//  determine range (yq[0], yq[1]) of time differences
//  where fit should be performed;
//  require a fraction >=min of all entries to lie below yq[0]
//      and a fraction <=max of all entries to lie below yq[1];
//  within the range (yq[0], yq[1]) there must be no empty bin;
//
void MHEffOnTime::GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const
{
#if ROOT_VERSION_CODE < ROOT_VERSION(3,02,07)
    // WOrkaround for missing const qualifier of TH1::Integral
    TH1 &h1 = (TH1&)h;

    // choose pedestrian approach as long as GetQuantiles is
    // not available
    const Int_t    jbins = h1.GetNbinsX();
    const Double_t Nm    = h1.Integral();

    const Double_t xq[2] = { 0.15*Nm, 0.98*Nm };

    yq[0] = yq[1] = h1.GetBinLowEdge(jbins+1);

    for (int j=1; j<=jbins; j++)
        if (h1.Integral(2, j) >= xq[0])
        {
            yq[0] = h1.GetBinLowEdge(j);
            break;
        }

    for (int j=1; j<=jbins; j++)
        if (h1.Integral(1, j) >= xq[1] || h1.GetBinContent(j)==0)
        {
            yq[1] = h1.GetBinLowEdge(j);
            break;
        }
#else
    // GetQuantiles doesn't seem to be available in root 3.01/06
    Double_t xq[2] = { min, max };
    ((TH1&)h).GetQuantiles(2, yq, xq);
#endif
}

void MHEffOnTime::DrawBin(TH1 &h, Int_t i) const
{
    TString strg1 = fVarname+"-bin #";
    strg1 += i;

    new TCanvas(strg1, strg1);

    gPad->SetLogy();
    gStyle->SetOptFit(1011);

    TString name="Bin_";
    name += i;

    h.SetName(name);
    h.SetXTitle("\\Delta t [s]");
    h.SetYTitle("Counts");
    h.DrawCopy();
}

Bool_t MHEffOnTime::CalcResults(const TF1 &func, Double_t Nm, Int_t i)
{
    const Double_t lambda = func.GetParameter(0);
    const Double_t N0     = func.GetParameter(1);
    const Double_t prob   = func.GetProb();
    const Int_t    NDF    = func.GetNDF();

    Double_t xmin, xmax;
    ((TF1&)func).GetRange(xmin, xmax);

    *fLog << inf;
    *fLog << "Fitted bin #" << i << " from " << xmin << " to " << xmax;
    *fLog << ",  got: lambda=" << lambda << "Hz N0=" << N0 << endl;

    if (prob<=0.01)
    {
        *fLog << warn << "WARN - Fit bin#" << i << " gives:";
        *fLog << " Chi^2-Probab(" << prob << ")<0.01";
        *fLog << " NoFitPts=" << func.GetNumberFitPoints();
        *fLog << " Chi^2=" << func.GetChisquare();
    }

    // was fit successful ?
    if (NDF<=0 || /*prob<=0.001 ||*/ lambda<=0 || N0<=0)
    {
        *fLog << warn << dbginf << "Fit failed bin #" << i << ": ";
        if (NDF<=0)
            *fLog << " NDF(Number of Degrees of Freedom)=0";
        if (lambda<=0)
            *fLog << " Parameter#0(lambda)=0";
        if (N0<=0)
            *fLog << " Parameter#1(N0)=0";
        *fLog << endl;

        return kFALSE;
    }

    //
    // -------------- start error calculation ----------------
    //
    Double_t emat[2][2];
    gMinuit->mnemat(&emat[0][0], 2);

    //
    // Rdead : fraction of real time lost by the dead time
    // kappa = number of observed events (Nm) divided by
    //         the number of genuine events (N0)
    // Teff  : effective on-time
    //
    const Double_t Teff  = Nm/lambda;
    const Double_t kappa = Nm/N0;
    const Double_t Rdead = 1.0 - kappa;

    const Double_t dldl   = emat[0][0];
    const Double_t dN0dN0 = emat[1][1];

    const Double_t dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
    const Double_t dl = sqrt(dldl);
    const Double_t dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
    //
    // -------------- end error calculation ----------------
    //

    // the effective on time is Nm/lambda
    fHEffOn.SetBinContent(i,  Teff);
    fHEffOn.SetBinError  (i, dTeff);

    // plot chi2-probability of fit
    fHProb.SetBinContent(i, prob);

    // lambda from fit
    fHLambda.SetBinContent(i, lambda);
    fHLambda.SetBinError  (i,     dl);

    // Rdead from fit
    fHRdead.SetBinContent(i, Rdead);
    fHRdead.SetBinError  (i,dRdead);

    return kTRUE;
}

void MHEffOnTime::ResetBin(Int_t i)
{
    fHEffOn.SetBinContent (i, 1.e-20);
    fHProb.SetBinContent  (i, 1.e-20);
    fHLambda.SetBinContent(i, 1.e-20);
    fHRdead.SetBinContent (i, 1.e-20);
}

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

    for (int i=1; i<=nbins; i++)
    {
        TH1D &h = *((TH2D*)hist)->ProjectionX(TString("Calc-")+fVarname,
                                              i, i, "E");
        if (draw)
            DrawBin(h, i);

        ResetBin(i);

        // Nmdel = Nm * binwidth,  with Nm = number of observed events
        const Double_t Nm    = h.Integral();
        const Double_t Nmdel = h.Integral("width");
        if (Nm <= 0)
        {
            *fLog << warn << dbginf << "Nm<0 for bin #" << i << endl;
            delete &h;
            continue;
        }

        Double_t yq[2];
        GetQuantiles(h, 0.15, 0.95, yq);

        //
        // Setup Poisson function for the fit:
        // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
        //
        TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
        func.SetParNames("lambda", "N0", "del");

        func.SetParameter(0, 1);
        func.SetParameter(1, Nm);
        func.FixParameter(2, Nmdel/Nm);

        // For me (TB) it seems that setting parameter limits isn't necessary

        // For a description of the options see TH1::Fit
        h.Fit("Poisson", "0IRQ");

        // Calc and fill results of fit into the histograms
        const Bool_t rc = CalcResults(func, Nm, i);

        // draw the distribution of time differences if requested
        if (rc && draw)
            func.DrawCopy("same");

        delete &h;
    }
}

// -------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histograms
//
Bool_t MHEffOnTime::SetupFill(const MParList *plist)
{
    TString bn = "Binning"+fVarname;

    const MBinning* binsVar = (MBinning*)plist->FindObject(bn);

    if (!binsVar)
    {
        *fLog << err << dbginf << bn << " [MBinning] not found... aborting." << endl;
        return kFALSE;
    }

    SetBinning(&fHEffOn,  binsVar);
    SetBinning(&fHProb,   binsVar);
    SetBinning(&fHLambda, binsVar);
    SetBinning(&fHRdead,  binsVar);

    fHEffOn.Sumw2();
    fHProb.Sumw2();
    fHLambda.Sumw2();
    fHRdead.Sumw2();

    return kTRUE;
}

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHEffOnTime::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);
    fHProb.Draw(opt);

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

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

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

void MHEffOnTime::Calc(const MHTimeDiffTheta &hist, const Bool_t Draw)
{
    Calc(hist.GetHist(), Draw);
}