source: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc@ 5672

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

//////////////////////////////////////////////////////////////////////////////
//
// MAlphaFitter
//
// 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 "MAlphaFitter.h"

#include <TF1.h>
#include <TH1.h>

#include <TLatex.h>

#include "MMath.h"

#include "MLogManip.h"

ClassImp(MAlphaFitter);

using namespace std;

// --------------------------------------------------------------------------
//
// 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(hist,    0, 1.25*sigint)
//   b = int(pol2(3), 0, 1.25*sigint)
//
// The Significance is calculated using the Significance() member
// function.
//
Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
{
    Double_t sigmax=fSigMax;
    Double_t bgmin =fBgMin;
    Double_t bgmax =fBgMax;

    //*fLog << inf << "Fit: " << sigmax << " " << fSigInt << "  " << bgmin << " " << bgmax << endl;

    //TF1 fFunc("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90);

    //fFunc->SetParameters(fCoefficients.GetArray());

    fFunc->FixParameter(1, 0);
    fFunc->FixParameter(4, 0);
    fFunc->SetParLimits(2, 0, 90);
    fFunc->SetParLimits(3, -1, 1);

    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
    fFunc->FixParameter(0, 0);
    fFunc->FixParameter(2, 1);
    fFunc->ReleaseParameter(3);

    for (int i=5; i<fFunc->GetNpar(); i++)
        fFunc->ReleaseParameter(i);

    // options : N  do not store the function, do not draw
    //           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(fFunc, "NQI", "", bgmin, bgmax);

    fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();

    // ------------------------------------
    if (paint)
    {
        fFunc->SetRange(0, 90);
        fFunc->SetLineColor(kRed);
        fFunc->SetLineWidth(2);
        fFunc->Paint("same");
    }
    // ------------------------------------

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

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

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

    // options : N  do not store the function, do not draw
    //           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(fFunc, "NQI", "", 0, sigmax);

    fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
    fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());

    //const Bool_t ok = NDF>0 && chi2<2.5*NDF;

    // ------------------------------------
    if (paint)
    {
        fFunc->SetLineColor(kGreen);
        fFunc->SetLineWidth(2);
        fFunc->Paint("same");
    }
    // ------------------------------------

     //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw;
     fFunc->SetParameter(0, 0);
     fFunc->SetParameter(2, 1);
     //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw;
     //fSignificance = MMath::SignificanceLiMaSigned(s, b);


    const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);

    fIntegralMax      = h.GetBinLowEdge(bin+1);
    fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
    fEventsSignal     = h.Integral(0, bin);
    fEventsExcess     = fEventsSignal-fEventsBackground;
    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);

    if (fEventsExcess<0)
        fEventsExcess=0;

    return kTRUE;
}

void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
{
    TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f  \\omega=%.1f\\circ  E=%d  (\\alpha<%.1f\\circ)  (\\chi_{b}^{2}=%.1f  \\chi_{s}^{2}=%.1f)",
                           fSignificance, GetGausSigma(),
                           (int)fEventsExcess, fIntegralMax,
                           fChiSqBg, fChiSqSignal));

    text.SetBit(TLatex::kTextNDC);
    text.SetTextSize(size);
    text.Paint();
}

void MAlphaFitter::Copy(TObject &o) const
{
    MAlphaFitter &f = static_cast<MAlphaFitter&>(o);

    f.fSigInt       = fSigInt;
    f.fSigMax       = fSigMax;
    f.fBgMin        = fBgMin;
    f.fBgMax        = fBgMax;
    f.fPolynomOrder = fPolynomOrder;
    f.fCoefficients.Set(fCoefficients.GetSize());
    f.fCoefficients.Reset();

    TF1 *fcn = f.fFunc;
    f.fFunc = new TF1(*fFunc);
    gROOT->GetListOfFunctions()->Remove(f.fFunc);
    f.fFunc->SetName("Dummy");
    delete fcn;
}

void MAlphaFitter::Print(Option_t *o) const
{
    *fLog << inf << GetDescriptor() << ": Fitting..." << endl;
    *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
    *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
    *fLog << " ...polynom order " << fPolynomOrder << endl;
}