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

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

#include <TRandom.h>

#include <TLatex.h>

#include "MMath.h"

#include "MLogManip.h"

ClassImp(MAlphaFitter);

using namespace std;

void MAlphaFitter::Clear(Option_t *o)
{
    fSignificance=0;
    fEventsExcess=0;
    fEventsSignal=0;
    fEventsBackground=0;

    fChiSqSignal=0;
    fChiSqBg=0;
    fIntegralMax=0;

    fCoefficients.Reset();
}

// --------------------------------------------------------------------------
//
// 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)
{
    Clear();
    if (h.GetEntries()==0)
        return kFALSE;

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

Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Bool_t paint)
{
    /*
     Clear();
     if (hon.GetEntries()==0)
     return kFALSE;
     */

    TH1D h(hon);
    h.Add(&hof, -1);

    MAlphaFitter fit(*this);
    fit.SetPolynomOrder(1);

    if (!fit.Fit(h, paint))
        return kFALSE;

    fChiSqSignal = fit.GetChiSqSignal();
    fChiSqBg     = fit.GetChiSqBg();
    fCoefficients = fit.GetCoefficients();

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

    fIntegralMax      = hon.GetBinLowEdge(bin+1);
    fEventsBackground = hof.Integral(0, bin);
    fEventsSignal     = hon.Integral(0, bin);
    fEventsExcess     = fEventsSignal-fEventsBackground;
    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);

    if (fEventsExcess<0)
        fEventsExcess=0;
/*
    TF1 func("", "gaus(0)+pol0(3)", 0., 90.);

    const Double_t A  = fEventsSignal/bin;
    const Double_t dA = TMath::Abs(A);
    func.SetParLimits(0, -dA*4, dA*4);
    func.SetParLimits(2, 0, 90);
    func.SetParLimits(3, -dA, dA);

    func.SetParameter(0, A);
    func.FixParameter(1, 0);
    func.SetParameter(2, fSigMax*0.75);
    func.SetParameter(3, 0);

    // 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
    TH1D h(hon);
    h.Add(&hof, -1);
    h.Fit(&func, "NQI", "", 0, 90);

    fChiSqSignal = func.GetChisquare()/func.GetNDF();

    const Int_t bin1 = h.GetXaxis()->FindFixBin(func.GetParameter(2)*2);

    fChiSqBg = 0;
    for (int i=bin1; i<=h.GetNbinsX(); i++)
    {
        const Float_t val = h.GetBinContent(i);
        fChiSqBg = val*val;
    }
    if (fChiSqBg>0)
        fChiSqBg /= h.GetNbinsX()+1-bin1;

    fCoefficients.Set(func.GetNpar(), func.GetParameters());

    // ------------------------------------
    if (paint)
    {
        func.SetLineColor(kBlue);
        func.SetLineWidth(2);
        func.Paint("same");
    }
    // ------------------------------------
  */
    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}/ndf=%.1f  \\chi_{s}^{2}/ndf=%.1f  c_{0}=%.1f)",
                           fSignificance, GetGausSigma(),
                           (int)fEventsExcess, fIntegralMax,
                           fChiSqBg, fChiSqSignal, fCoefficients[3]));

    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.fScaleMin     = fScaleMin;
    f.fScaleMax     = fScaleMax;
    f.fPolynomOrder = fPolynomOrder;
    f.fScaleMode    = fScaleMode;
    f.fScaleUser    = fScaleUser;
    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;
    *fLog << " ...scale mode: ";
    switch (fScaleMode)
    {
    case kNone:        *fLog << "none.";         break;
    case kEntries:     *fLog << "entries.";      break;
    case kIntegral:    *fLog << "integral.";     break;
    case kOffRegion:   *fLog << "off region.";   break;
    case kLeastSquare: *fLog << "least square."; break;
    case kUserScale:   *fLog << "user def (" << fScaleUser << ")"; break;
    }
    *fLog << endl;

    if (TString(o).Contains("result"))
    {
        *fLog << "Result:" << endl;
        *fLog << " - Significance           " << fSignificance << endl;
        *fLog << " - Excess Events          " << fEventsExcess << endl;
        *fLog << " - Signal Events          " << fEventsSignal << endl;
        *fLog << " - Background Events      " << fEventsBackground << endl;
        *fLog << " - Chi^2/ndf (Signal)     " << fChiSqSignal << endl;
        *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl;
        *fLog << " - Signal integrated      " << fIntegralMax << "°" << endl;
    }
}

Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
{
    const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}

Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
{
    const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}

Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
{
    const TString name(Form("TempAlpha%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, -1, 9999, -1, 9999, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}

Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    Scale(*h0, *h1);

    const Bool_t rc = Fit(*h1, *h0, paint);

    delete h0;
    delete h1;

    return rc;
}

Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));

    TH1D *h1 = hon.ProjectionZ(name1, bin, bin, -1, 9999, "E");
    TH1D *h0 = hof.ProjectionZ(name0, bin, bin, -1, 9999, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    Scale(*h0, *h1);

    const Bool_t rc = Fit(*h1, *h0, paint);

    delete h0;
    delete h1;

    return rc;
}

Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));

    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E");
    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    Scale(*h0, *h1);

    const Bool_t rc = Fit(*h1, *h0, paint);

    delete h0;
    delete h1;

    return rc;
}

void MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
{
    Float_t scaleon = 1;
    Float_t scaleof = 1;
    switch (fScaleMode)
    {
    case kNone:
        return;

    case kEntries:
        scaleon = on.GetEntries();
        scaleof = of.GetEntries();
        break;

    case kIntegral:
        scaleon = on.Integral();
        scaleof = of.Integral();
        break;

    case kOffRegion:
        {
            const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
            const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
            scaleon = on.Integral(min, max);
            scaleof = of.Integral(min, max);
        }
        break;

    case kUserScale:
        scaleon = fScaleUser;
        break;

    default:
        return;
    }

    if (scaleof!=0)
        of.Scale(scaleon/scaleof);
}