source: trunk/MagicSoft/Mars/mhist/MHAlpha.cc@ 4896

/* ======================================================================== *\
!
! *
! * 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, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHAlpha
//
// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
// check the result when it is filles in the MStatusDisplay
// For more information see MHFalseSource::FitSignificance
//
// For convinience (fit) the output significance is stored in a
// container in the parlisrt
//
// PRELIMINARY!
//
//////////////////////////////////////////////////////////////////////////////
#include "MHAlpha.h"

#include <TF1.h>
#include <TH2.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TStopwatch.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MHillasSrc.h"
#include "MTime.h"
#include "MObservatory.h"
#include "MPointingPos.h"
#include "MAstroSky2Local.h"
#include "MStatusDisplay.h"
#include "MParameters.h"
#include "MHMatrix.h"

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

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

ClassImp(MHAlpha);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHAlpha::MHAlpha(const char *name, const char *title)
    : fAlphaCut(12.5), fBgMean(55), fResult(0), fMatrix(0)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHAlpha";
    fTitle = title ? title : "Alpha plot";

    fHist.SetDirectory(NULL);

    fHist.SetName("Alpha");
    fHist.SetTitle("Alpha");
    fHist.SetXTitle("|\\alpha| [\\circ]");
    fHist.SetYTitle("Counts");
    fHist.UseCurrentStyle();

    MBinning binsa;
    binsa.SetEdges(18, 0, 90);
    binsa.Apply(fHist);

    fSigInt=15;
    fSigMax=75;
    fBgMin=45;
    fBgMax=90;
    fPolynom=1;
}

Bool_t MHAlpha::SetupFill(const MParList *pl)
{
    fHist.Reset();

    fResult = (MParameterD*)pl->FindObject("Significance", "MParameterD");
    if (!fResult)
        *fLog << warn << "Significance [MParameterD] not found... ignored." << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
{
    Double_t alpha;

    if (fMatrix)
        alpha = (*fMatrix)[fMap];
    else
    {
        const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
        if (!par)
        {
            *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
            return kFALSE;
        }

        alpha = hil->GetAlpha();
    }

    fHist.Fill(TMath::Abs(alpha), w);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Update the projections
//
void MHAlpha::Paint(Option_t *opt)
{
    DoFit();
}

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

    fHist.Draw();

    AppendPad("");
}

// --------------------------------------------------------------------------
//
// This is a preliminary implementation of a alpha-fit procedure for
// all possible source positions. It will be moved into its own
// more powerfull class soon.
//
// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
// or
//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
//
// Parameter [1] is fixed to 0 while the alpha peak should be
// symmetric around alpha=0.
//
// Parameter [4] is fixed to 0 because the first derivative at
// alpha=0 should be 0, too.
//
// In a first step the background is fitted between bgmin and bgmax,
// while the parameters [0]=0 and [2]=1 are fixed.
//
// In a second step the signal region (alpha<sigmax) is fittet using
// the whole function with parameters [1], [3], [4] and [5] fixed.
//
// The number of excess and background events are calculated as
//   s = int(0, sigint, gaus(0)+pol2(3))
//   b = int(0, sigint,         pol2(3))
//
// The Significance is calculated using the Significance() member
// function.
//
Bool_t MHAlpha::DoFit(Double_t &sig, Bool_t paint)
{
    Float_t sigint=fSigInt;
    Float_t sigmax=fSigMax;
    Float_t bgmin=fBgMin;
    Float_t bgmax=fBgMax;
    Byte_t polynom=fPolynom;

    // Implementing the function yourself is only about 5% faster
    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
    TArrayD maxpar(func.GetNpar());

    func.FixParameter(1, 0);
    func.FixParameter(4, 0);
    func.SetParLimits(2, 0, 90);
    func.SetParLimits(3, -1, 1);

    TH1 *h=&fHist;
    const Double_t alpha0 = h->GetBinContent(1);
    const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);

    // Check for the regios which is not filled...
    if (alpha0==0)
        return kFALSE; //*fLog << warn << "Histogram empty." << endl;

    // First fit a polynom in the off region
    func.FixParameter(0, 0);
    func.FixParameter(2, 1);
    func.ReleaseParameter(3);

    for (int i=5; i<func.GetNpar(); i++)
        func.ReleaseParameter(i);

    h->Fit(&func, "N0Q", "", bgmin, bgmax);

    const Double_t chisq1 = func.GetChisquare()/func.GetNDF();

    // ------------------------------------
    if (paint)
    {
        func.SetRange(0, 90);
        func.SetLineColor(kRed);
        func.Paint("same");
    }
    // ------------------------------------

    func.ReleaseParameter(0);
    //func.ReleaseParameter(1);
    func.ReleaseParameter(2);
    func.FixParameter(3, func.GetParameter(3));
    for (int i=5; i<func.GetNpar(); i++)
        func.FixParameter(i, func.GetParameter(i));

    // Do not allow signals smaller than the background
    const Double_t A  = alpha0-func.GetParameter(3);
    const Double_t dA = TMath::Abs(A);
    func.SetParLimits(0, -dA*4, dA*4);
    func.SetParLimits(2, 0, 90);

    // Now fit a gaus in the on region on top of the polynom
    func.SetParameter(0, A);
    func.SetParameter(2, sigmax*0.75);

    h->Fit(&func, "N0Q", "", 0, sigmax);

    const Double_t chisq2    = func.GetChisquare()/func.GetNDF();
    const Double_t sigmagaus = func.GetParameter(2);

    // ------------------------------------
    if (paint)
    {
        func.SetLineColor(kGreen);
        func.Paint("same");
    }
    // ------------------------------------

    const Double_t s   = func.Integral(0, sigint)/alphaw;
    func.SetParameter(0, 0);
    func.SetParameter(2, 1);
    const Double_t b   = func.Integral(0, sigint)/alphaw;

    sig = MMath::SignificanceLiMaSigned(s, b);

    // ------------------------------------
    // This is always draw one update too late... why?
    //fHist.SetTitle(Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ)  \\omega=%.1f\\circ  E=%d  (\\chi_{b}^{2}/N=%.1f  \\chi_{s}^{2}/N=%.1f)",
    //                    sig, fSigInt, sigmagaus, (int)(s-b), chisq1, chisq2));

    if (paint)
    {
        TLatex text(0.13, 0.94, Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ)  \\omega=%.1f\\circ  E=%d  (\\chi_{b}^{2}/N=%.1f  \\chi_{s}^{2}/N=%.1f)",
                                     sig, fSigInt, sigmagaus, (int)(s-b), chisq1, chisq2));
        text.SetBit(TLatex::kTextNDC);
        text.SetTextSize(0.04);
        text.Paint();
    }
    // ------------------------------------

    return kTRUE;
}

Bool_t MHAlpha::Finalize()
{
    Double_t sig = 0;

    if (!DoFit(sig))
    {
        *fLog << warn << "Histogram empty." << endl;
        return kFALSE;
    }

    if (fResult)
        fResult->SetVal(sig);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of
// MMcEvt. This function adds all necessary columns to the
// given matrix. Afterward you should fill the matrix with the corresponding
// data (eg from a file by using MHMatrix::Fill). If you now loop
// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
// will take the values from the matrix instead of the containers.
//
void MHAlpha::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
        return;

    fMatrix = mat;

    fMap = fMatrix->AddColumn("MHillasSrc.fAlpha");
}

void MHAlpha::StopMapping()
{
    fMatrix = NULL; 
}