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

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

ClassImp(MHEffOnTime);


// --------------------------------------------------------------------------
//
// 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);
    strg += "  ";
    strg += fUnit;

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

    // effective on time versus Var
    fHEffOn.SetName("EffOn");
    TString strg2("effective On Time vs. ");
    strg2 += fVarname;
    fHEffOn.SetTitle(strg2);

    fHEffOn.SetDirectory(NULL);

    fHEffOn.SetXTitle(strg);
    fHEffOn.SetYTitle("effective ontime [s]");

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

    fHProb.SetDirectory(NULL);

    fHProb.SetXTitle(strg);
    fHProb.SetYTitle("chi2-probability");

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

    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");
    TString strg5("Rdead of  OnTimeFit vs. ");
    strg5 += fVarname;
    fHRdead.SetTitle(strg5);

    fHRdead.SetDirectory(NULL);

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

}

// -----------------------------------------------------------------------
//
// Calculate the effective on time by fitting the distribution of
// time differences
//
void MHEffOnTime::Calc(TH2D *hist, const Bool_t Draw)
{
     // Draw != 0   means the distribution of time differences should be drawn

    // nbins = number of Var bins
    const Int_t nbins = hist->GetNbinsY();

    // start  of 'for' loop ---------------------------
    for (int i=1; i<=nbins; i++)
    {
      TString strg0("Calc-");
      strg0 += fVarname;
      TH1D &h = *hist->ProjectionX(strg0, 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();

        //...................................................
        // 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 Int_t    jbins  = h.GetNbinsX();

        fHEffOn.SetBinContent (i, 1.e-20);
        fHProb.SetBinContent  (i, 1.e-20);
        fHLambda.SetBinContent(i, 1.e-20);
        fHRdead.SetBinContent (i, 1.e-20);

        // start  of 'if' loop ---------------------------
        if (Nm > 0.0)
        {
            // del is bin width in time difference
            const Double_t del = Nmdel/Nm;

            Double_t sum1 = 0.0;
            yq[0]  = h.GetBinLowEdge(jbins+1);
            for (int j=1; j<=jbins; j++)
            {
                if (sum1 >= xq[0]*Nm)
                {
                    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++)
            {
                const Double_t content = h.GetBinContent(j);
                sum2 += content;
                if (sum2 >= xq[1]*Nm || content == 0.0)
                {
                    yq[1] = h.GetBinLowEdge(j);
                    break;
                }
            }

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

            // parameter 0 = lambda
            // parameter 1 = N0           with N0 = ideal number of events
            // parameter 2 = del (fixed)  with del = bin width of time difference

            TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
            func.FixParameter(2, del); 

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

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

            func.SetParName(0, "lambda");
            func.SetParName(1, "N0");
            func.SetParName(2, "del");

            // options : 0  (=zero) 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");

            const Double_t lambda = func.GetParameter(0);
            const Double_t N0     = func.GetParameter(1);
            // const Double_t chi2   = func.GetChisquare();
            const Double_t Prob   = func.GetProb();
            const Int_t    NDF    = func.GetNDF();

            // was fit successful ?
            if (NDF>0  &&  Prob>0.01  &&  lambda>0.0  &&  N0>0.0)
            {
                // start error calculation .....................

                Double_t emat[2][2];
                // Double_t eplus, eminus, eparab, gcc;
                gMinuit->mnemat(&emat[0][0], 2);

                // for (int m=0; m<2; m++)
                // {
                //   *fLog << "emat[" << m << "][n] = "; 
                //   for (int n=0; n<2; n++)
                //   {
                //     *fLog << emat[m][n] << "  ";
                //   }
                //   *fLog << endl;
                // }

                // *fLog << "eplus, eminus, eparab, gcc :" << endl; 
                // for (int k=0; k<2; k++)
                // {
                //   gMinuit->mnerrs(k, eplus, eminus, eparab, gcc);
                //   *fLog << eplus << " " << eminus << " " << eparab << " "
                //         << gcc << endl;
                // }

                // 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

                Double_t Rdead, kappa, Teff, dTeff, dldl, dl, dRdead;
                Double_t dN0dN0;
                //Double_t dldN0,lN0corr;
                Teff  = Nm/lambda;
                kappa = Nm/N0;
                Rdead = 1.0 - kappa;

                dldl   = emat[0][0];
                dN0dN0 = emat[1][1];
                // dldN0  = emat[0][1];

                dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
                dl = sqrt(dldl);
                dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);

                // lN0corr = dldN0 / sqrt(dldl * dN0dN0);
                // *fLog << "MHEffOnTime : correlation between lambda and N0 = "
                //       << lN0corr << endl;

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

            //........................
            // draw the distribution of time differences if requested
            //
            if (Draw == kTRUE)
              {
                char txt[100];
                TString strg1(fVarname);
                strg1 += "-bin %d";
                sprintf(txt, strg1, i);
                new TCanvas(txt, txt);
                // TCanvas &c = *MakeDefCanvas(txt, txt);
                // gROOT->SetSelectedPad(NULL);

                gPad->SetLogy();
                gStyle->SetOptFit(1011);
                // gStyle->SetOptStat(1);

                h.SetName(txt);
                h.SetXTitle("time difference [s]");
                h.SetYTitle("Counts");
                h.DrawCopy();

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

                // c.Modified();
                // c.Update();
              }
            //........................
        }
        // end  of 'if' loop ---------------------------

        delete &h;
    }
    // end  of 'for' loop ---------------------------
    //    gStyle->SetOptStat(1111);

}

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

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

    if (!binsVar)
    {
        *fLog << err << dbginf << strg0 
              << " [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;
}

// -------------------------------------------------------------------------
//
// Dummy Fill
//
Bool_t MHEffOnTime::Fill(const MParContainer *par)
{
  return kTRUE;
}

// -------------------------------------------------------------------------
//
// Draw a copy of the histogram
//
TObject *MHEffOnTime::DrawClone(Option_t *opt) const
{
    TString strg0("EffOnTime");
    strg0 += fVarname;    
    TString strg1("Results from on time fit vs. ");
    strg1 += fVarname;    

    TCanvas &c = *MakeDefCanvas(strg0, strg1);
    c.Divide(2, 2);

    gROOT->SetSelectedPad(NULL);

    c.cd(1);
    ((TH2*)&fHEffOn)->DrawCopy(opt);

    c.cd(2);
    ((TH2*)&fHProb)->DrawCopy(opt);

    c.cd(3);
    ((TH2*)&fHLambda)->DrawCopy(opt);

    c.cd(4);
    ((TH2*)&fHRdead)->DrawCopy(opt);

    c.Modified();
    c.Update();

    return &c;
}

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHEffOnTime::Draw(Option_t *opt) 
{
    TString strg0("EffOnTime");
    strg0 += fVarname;    
    TString strg1("Results from on time fit vs. ");
    strg1 += fVarname;    

    if (!gPad)
        MakeDefCanvas(strg0, strg1);

    gPad->Divide(2,2);

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

    gPad->cd(2);
    fHProb.Draw(opt);

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

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

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