source: trunk/MagicSoft/Mars/mhflux/MHPhi.cc@ 7574

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

/////////////////////////////////////////////////////////////////////////////
//
// MHPhi
//
// Plot delta phi the angle between the reconstructed shower origin and
// the source position.
//
// More detail can be found at:
// http://www.astro.uni-wuerzburg.de/results/ringmethod/
//
////////////////////////////////////////////////////////////////////////////
#include "MHPhi.h"

#include <TCanvas.h>
#include <TVector2.h>
#include <TLatex.h>
#include <TLine.h>
#include <TMarker.h>

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

#include "MParList.h"

#include "MHillas.h"
#include "MSrcPosCam.h"
#include "MParameters.h"
#include "MGeomCam.h"
#include "MBinning.h"
#include "MMath.h"

ClassImp(MHPhi);

using namespace std;

// --------------------------------------------------------------------------
//
// Setup histograms 
//
MHPhi::MHPhi(const char *name, const char *title)
: fHillas(0), fSrcPos(0), fDisp(0)//, fOnOffMode(kTRUE), fIsOffLoop(kFALSE)
{
    fName  = name  ? name  : "MHPhi";
    fTitle = title ? title : "Graphs for rate data";

    // Init Graphs
    fHPhi.SetNameTitle("Phi", "\\Delta\\Phi-Distribution");

    fHPhi.SetXTitle("\\Delta\\Phi [\\circ]");
    fHPhi.SetYTitle("Counts");

    fHPhi.SetMinimum(0);
    fHPhi.SetDirectory(0);
    fHPhi.Sumw2();
    fHPhi.SetBit(TH1::kNoStats);
    fHPhi.SetMarkerStyle(kFullDotMedium);

    fHPhi.GetYaxis()->SetTitleOffset(1.2);

    /*
    fNameParameter = "Disp";

    fHist.SetNameTitle("Phi", "\\Delta\\Phi-Distribution");
    fHist.SetZTitle("\\Delta\\Phi [\\circ]");
    fHist.SetDirectory(NULL);

    // Main histogram
    fHistTime.SetName("Phi");
    fHistTime.SetXTitle("\\Delta\\Phi [\\circ]");
    fHistTime.SetDirectory(NULL);

    MBinning binsa, binse, binst;
    //binsa.SetEdges(75, 0, 1.5);
    //binsa.SetEdges(arr);
    binse.SetEdgesLog(15, 10, 100000);
    binst.SetEdgesASin(67, -0.005, 0.665);
    //binsa.Apply(fHistTime);

    MH::SetBinning(&fHist, &binst, &binse, &binsa);
     */
}

/*
Double_t MHPhi::GetVal() const
{
    const Dopuble_t disp = static_cast<const MParameterD*>(fParameter)->GetVal();

    const TVector2 pos = fHillas->GetMean()*fConvMm2Deg + fHillas->GetNormAxis()*disp;
    const TVector2 src = fSrcPos->GetXY()*fConvMm2Deg;

    // Calculate radial distance.
    const Double_t d = pos.Mod() - src.Mod();

    if (d<-fThetaCut*0.913 || d>fThetaCut)
        return kTRUE;

    const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();
    const Double_t absd  = TMath::Abs(delta)

    return fHistOff ? absd : 180-absd;
}
*/

// --------------------------------------------------------------------------
//
// Setup the Binning for the histograms automatically if the correct
// instances of MBinning
//
Bool_t MHPhi::SetupFill(const MParList *plist)
{
    fHillas = (MHillas*)plist->FindObject("MHillas");
    if (!fHillas)
    {
        *fLog << err << "MHillas not found... abort." << endl;
        return kFALSE;
    }
    fSrcPos = (MSrcPosCam*)plist->FindObject("MSrcPosCam");
    if (!fSrcPos)
    {
        *fLog << err << "MSrcPosCam not found... abort." << endl;
        return kFALSE;
    }
    fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD");
    if (!fDisp)
    {
        *fLog << err << "Disp [MParameterD] not found... abort." << endl;
        return kFALSE;
    }

    MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!geom)
    {
        *fLog << err << "MGeomCam not found... abort." << endl;
        return kFALSE;
    }

    fConvMm2Deg = geom->GetConvMm2Deg();

    fNumBinsSignal = 2;
    fThetaCut      = 0.21/1.2;
    fDistSrc       = 0.4;
    //fIsOffLoop = !fIsOffLoop;

    const Double_t w  = TMath::ATan(fThetaCut/fDistSrc);
    const Float_t  sz = TMath::RadToDeg()*w/fNumBinsSignal;
    const Int_t    n  = TMath::Nint(TMath::Ceil(180/sz));

    MBinning(n, 0, n*sz).Apply(fHPhi);
    /*

    // Get Histogram binnings
    MBinning binst, binse;
    binst.SetEdges(fHist, 'x');
    binse.SetEdges(fHist, 'y');

    MBinning binsa(n, 0, n*sz);

    // Apply binning
    binsa.Apply(fHistTime);
    MH::SetBinning(&fHist, &binst, &binse, &binsa);

    // Remark: Binnings might be overwritten in MHAlpha::SetupFill
    return MHAlpha::SetupFill(pl);
     */
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histograms with data from a MMuonCalibPar and
// MMuonSearchPar container.
//
Bool_t MHPhi::Fill(const MParContainer *par, const Stat_t weight)
{
    const TVector2 pos = fHillas->GetMean()*fConvMm2Deg + fHillas->GetNormAxis()*fDisp->GetVal();
    const TVector2 src = fSrcPos->GetXY()*fConvMm2Deg;

    // Calculate radial distance.
    const Double_t d = pos.Mod() - src.Mod();

    if (d<-fThetaCut*0.913 || d>fThetaCut)
        return kTRUE;

    const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();

    fHPhi.Fill(TMath::Abs(delta), weight);

    // const Double_t absd = TMath::Abs(delta)
    // fHPhi.Fill(fHistOff ? absd : 180-absd, weight);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// This displays the TGraph like you expect it to be (eg. time on the axis)
// It should also make sure, that the displayed time always is UTC and
// not local time...
//
void MHPhi::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("combine");

    fHPhi.Draw();

    AppendPad(opt);
}

void MHPhi::Paint(Option_t *o)
{
    //TString opt(o);
    //opt.ToLower();

    // if (opt=="combine" && fHistOff)
    // {
    //    fHPhi.Add(fHist, fHistOff);
    //    return;
    // }

    const Bool_t wobble = TString(o).Contains("anticut", TString::kIgnoreCase);

    const Double_t cut  = fHPhi.GetBinLowEdge(fNumBinsSignal+1);

    const Int_t maxbin  = wobble ? fHPhi.GetXaxis()->FindFixBin(180-cut)-1 : fHPhi.GetNbinsX();
    const Double_t cut2 = wobble ? fHPhi.GetBinLowEdge(maxbin+1) : 180;

    const Double_t sig  = fHPhi.Integral(1, fNumBinsSignal);
    const Double_t bg   = fHPhi.Integral(1+fNumBinsSignal, maxbin);

    const Double_t f    = cut/(cut2-cut);

    const Double_t S0   = MMath::SignificanceLiMaSigned(sig, bg*f);
    const Double_t S    = MMath::SignificanceLiMaSigned(sig, bg, f);

    const TString  fmt  = Form("\\sigma_{L/M}=%.1f (\\sigma_{0}=%.1f)  \\Delta\\Phi_{on}<%.1f\\circ  \\Delta\\Phi_{off}<%.1f\\circ  E=%.0f  B=%.0f  f=%.2f",
                               S, S0, cut, cut2, sig-bg*f, bg*f, f);

    const Double_t b = bg             *f/fNumBinsSignal;
    const Double_t e = TMath::Sqrt(bg)*f/fNumBinsSignal;

    TLatex text(0.27, 0.94, fmt);
    text.SetBit(TLatex::kTextNDC);
    text.SetTextSize(0.035);
    text.Paint();

    TLine line;
    line.SetLineColor(14);
    line.PaintLine(cut, gPad->GetUymin(), cut, gPad->GetUymax());
    if (maxbin<fHPhi.GetNbinsX())
        line.PaintLine(cut2, gPad->GetUymin(), cut2, gPad->GetUymax());
    line.SetLineColor(kBlue);
    line.PaintLine(0, b, cut, b);
    line.PaintLine(cut/2, b-e, cut/2, b+e);
    line.SetLineStyle(kDashed);
    line.PaintLine(cut, b, fHPhi.GetXaxis()->GetXmax(), b);

    TMarker m;
    m.SetMarkerColor(kBlue);
    m.SetMarkerStyle(kFullDotMedium);
    m.PaintMarker(cut/2, b);
}