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

/* ======================================================================== *\
!
! *
! * 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 <TCanvas.h>
#include <TPaveText.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 "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;
}

// --------------------------------------------------------------------------
//
// Calculate Significance as
// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
// 
Double_t MHAlpha::Significance(Double_t s, Double_t b)
{
    const Double_t k = b==0 ? 0 : s/b;
    const Double_t f = s+k*k*b;

    return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
}

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

        alpha = hil->GetAlpha();
    }

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Update the projections
//
void MHAlpha::Paint(Option_t *opt)
{
    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...
    const Double_t alpha0 = h->GetBinContent(1);
    if (alpha0==0)
        return;

    // 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);
    h->Fit(&func, "N0Q", "", bgmin, bgmax);
    func.SetRange(0, 90);
    func.SetLineColor(kRed);
    func.Paint("same");

    //h4a.Fill(func.GetChisquare());
    //h5a.Fill(func.GetProb());

    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);
    h->Fit(&func, "N0Q", "", 0, sigmax);
    func.SetLineColor(kGreen);
    func.Paint("same");
}

// --------------------------------------------------------------------------
//
// 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::Finalize()
{
    //*fLog << dbg << "Fitting..." << endl;

    Float_t sigint=fSigInt;
    Float_t sigmax=fSigMax;
    Float_t bgmin=fBgMin;
    Float_t bgmax=fBgMax;
    Byte_t polynom=fPolynom;

    //                      xmin, xmax, npar
    //TF1 func("MyFunc", fcn, 0, 90, 6);
    // Implementing the function yourself is only about 5% faster
    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    TArrayD maxpar(func.GetNpar());

    //*fLog << "Setup pars..." << endl;

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

    //    TStopwatch clk;
    //    clk.Start();
   /*
    *fLog << inf;
    *fLog << "Signal fit:     alpha < " << sigmax << endl;
    *fLog << "Integration:    alpha < " << sigint << endl;
    *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
    *fLog << "Polynom order:  " << (int)polynom << endl;
    *fLog << "Fitting False Source Plot..." << flush;
    */
    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)
    {
        *fLog << warn << "Histogram empty." << endl;
        return kTRUE;
    }

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

    //*fLog << "Release pars..." << endl;

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

    //*fLog << dbg << "Fit 1" << endl;
    //h->Fit(&func, "N0Q", "", bgmin, bgmax);
    h->Fit(&func, "N0Q", "", bgmin, bgmax);
    //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;

    //h4a.Fill(func.GetChisquare());
    //h5a.Fill(func.GetProb());

    //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    //func.SetParLimits(2, 5, 90);

    //*fLog << dbg << "Change params..." << endl;

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

    //*fLog << dbg << "Fit 2..." << endl;
    //h->Fit(&func, "N0Q", "", 0, sigmax);
    h->Fit(&func, "N0Q", "", 0, sigmax);
    //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;

    //*fLog << dbg << "Copy par..." << endl;
    TArrayD p(func.GetNpar(), func.GetParameters());

    // Fill results into some histograms
    /*
    h0a.Fill(p[0]);
    h0b.Fill(p[3]);
    h1.Fill(p[1]);
    h2.Fill(p[2]);
    if (polynom>1)
        h3.Fill(p[5]);
    h4b.Fill(func.GetChisquare());
    //h5b.Fill(func.GetProb());
    */
    // Implementing the integral as analytical function
    // gives the same result in the order of 10e-5
    // and it is not faster at all...

    //*fLog << dbg << "Int1 par..." << endl;
    const Double_t s = func.Integral(0, sigint)/alphaw;

    func.SetParameter(0, 0);
    //func.SetParameter(2, 1);

    //*fLog << dbg << "Int2 par..." << endl;
    const Double_t b   = func.Integral(0, sigint)/alphaw;
    const Double_t sig = Significance(s, b);

    //*fLog << dbg << "Set Val "<<sig<<" ..." << fResult << endl;
    if (fResult)
        fResult->SetVal(sig);

    *fLog << inf << "Found Significance: " << sig << "  (width=";
    *fLog << func.GetParameter(2) << "°, excess=" << s-b << ")" << endl;

    //*fLog << dbg << "Done." << endl;
    return kTRUE;
    /*
    if (sig!=0)
        h6.Fill(sig);

    if (sig<maxs)
    {
        maxs = sig;
        maxx = ix;
        maxy = iy;
        maxpar = p;
        }
       */
/*

    *fLog << inf << "done." << endl;

    h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);

    //hists->SetMinimum(GetMinimumGT(*hists));
    histb->SetMinimum(GetMinimumGT(*histb));

//    clk.Stop();
//    clk.Print();
    TCanvas *c=new TCanvas;

    c->Divide(3,2, 0, 0);
    c->cd(1);
    gPad->SetBorderMode(0);
    hists->Draw("colz");
    hists->SetBit(kCanDelete);
    c->cd(2);
    gPad->SetBorderMode(0);
    hist->Draw("colz");
    hist->SetBit(kCanDelete);
    c->cd(3);
    gPad->SetBorderMode(0);
    histb->Draw("colz");
    histb->SetBit(kCanDelete);
    c->cd(4);
    gPad->Divide(1,3, 0, 0);
    TVirtualPad *p=gPad;
    p->SetBorderMode(0);
    p->cd(1);
    gPad->SetBorderMode(0);
    h0b.DrawCopy();
    h0a.DrawCopy("same");
    p->cd(2);
    gPad->SetBorderMode(0);
    h3.DrawCopy();
    p->cd(3);
    gPad->SetBorderMode(0);
    h2.DrawCopy();
    c->cd(6);
    gPad->Divide(1,2, 0, 0);
    TVirtualPad *q=gPad;
    q->SetBorderMode(0);
    q->cd(1);
    gPad->SetBorderMode(0);
    gPad->SetBorderMode(0);
    h4b.DrawCopy();
    h4a.DrawCopy("same");
    h6.DrawCopy("same");
    q->cd(2);
    gPad->SetBorderMode(0);
    //h5b.DrawCopy();
    //h5a.DrawCopy("same");
    c->cd(5);
    gPad->SetBorderMode(0);
    if (maxx>0 && maxy>0)
    {
        const char *title = Form(" \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f) ",
                                 fHist.GetXaxis()->GetBinCenter(maxx),
                                 fHist.GetYaxis()->GetBinCenter(maxy), maxs);

        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
        result->SetDirectory(NULL);
        result->SetNameTitle("Result \\alpha", title);
        result->SetBit(kCanDelete);
        result->SetXTitle("\\alpha [\\circ]");
        result->SetYTitle("Counts");
        result->Draw();

        TF1 f1("", func.GetExpFormula(), 0, 90);
        TF1 f2("", func.GetExpFormula(), 0, 90);
        f1.SetParameters(maxpar.GetArray());
        f2.SetParameters(maxpar.GetArray());
        f2.FixParameter(0, 0);
        f2.FixParameter(1, 0);
        f2.FixParameter(2, 1);
        f1.SetLineColor(kGreen);
        f2.SetLineColor(kRed);

        f2.DrawCopy("same");
        f1.DrawCopy("same");

        TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC");
        leg->SetBorderSize(1);
        leg->SetTextSize(0.04);
        leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
        //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
        leg->AddLine(0, 0.65, 0, 0.65);
        leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
        leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
        leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
        leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
        leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
        if (polynom>1)
            leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
        leg->SetBit(kCanDelete);
        leg->Draw();

        q->cd(2);

        TGraph *g = new TGraph;
        g->SetPoint(0, 0, 0);

        Int_t max=0;
        Float_t maxsig=0;
        for (int i=1; i<89; i++)
        {
            const Double_t s = f1.Integral(0, (float)i);
            const Double_t b = f2.Integral(0, (float)i);

            const Double_t sig = Significance(s, b);

            g->SetPoint(g->GetN(), i, sig);

            if (sig>maxsig)
            {
                max = i;
                maxsig = sig;
            }
        }

        g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
        g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
        g->GetHistogram()->SetYTitle("Significance");
        g->SetBit(kCanDelete);
        g->Draw("AP");

        leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
        leg->SetBorderSize(1);
        leg->SetTextSize(0.1);
        leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
        leg->SetBit(kCanDelete);
        leg->Draw();
    }*/
}

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