source: trunk/MagicSoft/Mars/mhist/MHHadronness.cc@ 1587

/* ======================================================================== *\
!
! *
! * 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, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Abelardo Moralejo <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHHadronness
//
// This is histogram is a way to evaluate the quality of a gamma/hadron
// seperation method. It is filled from a MHadronness container, which
// stores a hadroness for the current event. The Value must be in the
// range [0,1]. To fill the histograms correctly the information
// whether it is a gamma or hadron (not a gamma) must be available from
// a MMcEvt container.
//
// In the constructor you can change the number of used bns for the
// evaluation.
//
// The meaning of the histograms (Draw, DrawClone) are the following:
//  * Upper Left Corner:
//    - black: histogram of all hadronesses for gammas
//    - red:   histogram of all hadronesses for non gammas
//  * Upper Right Corner:
//    - black: acceptance of gammas (Ag) vs. the hadroness
//    - red:   acceptance of non gammas (Ah) vs. the hadroness
//    - blue:  2D distance of (acceptance_hadrons, acceptances_gammas)
//             to optimum (0, 1)
//             1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
//  * Bottom Left Corner:
//    Naive quality factor: Ag/sqrt(Ah)
//  * Bottom Right Corner:
//    - black: Acceprtance Gammas vs. Acceptance Hadrons
//    - blue cross: minimum of distance to (0, 1)
//
////////////////////////////////////////////////////////////////////////////
#include "MHHadronness.h"

#include <TPad.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMarker.h>

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

#include "MParList.h"
#include "MBinning.h"
#include "MHadronness.h"

#include "MMcEvt.hxx"

ClassImp(MHHadronness);

// --------------------------------------------------------------------------
//
// Setup histograms, nbins is the number of bins used for the evaluation.
// The default is 100 bins.
//
MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHHadronness";
    fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms";

    fGraph = new TGraph;
    fGraph->SetTitle("Acceptance Gammas vs. Hadrons");
    fGraph->SetMaximum(1);

    fGhness = new TH1D("Ghness", "Hadronness", nbins, 0, 1);
    fPhness = new TH1D("Phness", "Hadronness", nbins, 0, 1);
    fGhness->SetXTitle("Hadronness");
    fPhness->SetXTitle("Hadronness");
    fGhness->SetYTitle("Counts");
    fPhness->SetYTitle("Counts");
    
    fIntGhness = new TH1D("AccGammas",  "Acceptance", nbins, 0, 1);
    fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
    fIntGhness->SetXTitle("Hadronness");
    fIntPhness->SetXTitle("Hadronness");
    fIntGhness->SetYTitle("Acceptance");
    fIntPhness->SetYTitle("Acceptance");

    /*
    fQfac = new TH1D("Qfac", "Naive Quality factor", nbins, 0, 1);
    fQfac->SetXTitle("Hadronness");
    fQfac->SetYTitle("Quality");
    */
    fQfac = new TGraph;
    fQfac->SetTitle(" Naive Quality factor ");

    fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
    fMinDist->SetXTitle("Hadronness");
    fMinDist->SetYTitle("Distance");

    //fQfac->SetDirectory(NULL);
    fGhness->SetDirectory(NULL);
    fPhness->SetDirectory(NULL);
    fMinDist->SetDirectory(NULL);
    fIntGhness->SetDirectory(NULL);
    fIntPhness->SetDirectory(NULL);
}

// --------------------------------------------------------------------------
//
// Delete the histograms.
//
MHHadronness::~MHHadronness()
{
    delete fGhness;
    delete fIntGhness;
    delete fPhness;
    delete fIntPhness;
    delete fQfac;
    delete fMinDist;
    delete fGraph;
}

// --------------------------------------------------------------------------
//
// Setup Filling of the histograms. It needs:
//  MMcEvt and MHadronness
//
Bool_t MHHadronness::SetupFill(const MParList *plist)
{
    fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
    if (!fMcEvt)
    {
        *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
        return kFALSE;
    }

    fHadronness = (MHadronness*)plist->FindObject("MHadronness");
    if (!fHadronness)
    {
        *fLog << err << dbginf << "MHadronness not found... aborting." << endl;
        return kFALSE;
    }

    /*
     MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
     if (!bins)
     {
     *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
     return kFALSE;
     }

     SetBinning(&fHist, binsalpha, binsenergy, binstheta);

     fHist.Sumw2();
     */

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the Hadronness from a MHadronness container into the corresponding
// histogram dependant on the particle id.
//
Bool_t MHHadronness::Fill(const MParContainer *par)
{
    if (fMcEvt->GetPartId()==kGAMMA)
        fGhness->Fill(fHadronness->GetHadronness());
    else
        fPhness->Fill(fHadronness->GetHadronness());

    return kTRUE;
}

Float_t MHHadronness::GetQ05() const
{
    Int_t n = fQfac->GetN();

    Double_t val1x=0;
    Double_t val2x=1;

    Double_t val1y=0;
    Double_t val2y=0;

    for (Int_t i=1; i<=n; i++)
    {
        Double_t x, y;

        fQfac->GetPoint(i, x, y);

        if (x<0.5 && x>val1x)
        {
            val1x = x;
            val1y = y;
        }

        if (x>0.5 && x<val2x)
        {
            val2x = x;
            val2y = y;
        }
    }

    return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
}

// --------------------------------------------------------------------------
//
// Finalize the histograms:
//  - integrate the hadroness histograms --> acceptance
//  - fill the Minimum Distance histogram (formular see class description)
//  - fill the Quality histogram (formular see class description)
//
Bool_t MHHadronness::Finalize()
{
    Int_t n = fGhness->GetNbinsX();

    fGraph->Set(n);
    fQfac->Set(n);

    const Stat_t sumg = fGhness->Integral(1, n+1);
    const Stat_t sump = fPhness->Integral(1, n+1);

    Float_t max=0;

    for (Int_t i=1; i<=n; i++)
    {
        const Stat_t ip = fPhness->Integral(1, i)/sump;
        const Stat_t ig = fGhness->Integral(1, i)/sumg;

        fIntPhness->SetBinContent(i, ip);
        fIntGhness->SetBinContent(i, ig);

        fGraph->SetPoint(i, ip, ig);

        if (ip<=0)
            continue;

        fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));

        Double_t val = ig/sqrt(ip);
        fQfac->SetPoint(i, ig, val);

        if (val>max)
            max = val;
    }

    fQfac->SetMaximum(max*1.05);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Search the corresponding points for the given hadron acceptance (acchad)
// and interpolate the tow points (linear)
//
Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
{
    const Int_t n = fGraph->GetN();
    const Double_t *x = fGraph->GetX();
    const Double_t *y = fGraph->GetY();

    Int_t i = 0;
    while (i<n && x[i]<acchad)
        i++;

    if (i==0 || i==n)
        return 0;

    if (i==n-1)
        i--;

    const Double_t x1 = x[i-1];
    const Double_t y1 = y[i-1];

    const Double_t x2 = x[i];
    const Double_t y2 = y[i];

    return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
}

// --------------------------------------------------------------------------
//
// Search the corresponding points for the given gamma acceptance (accgam)
// and interpolate the tow points (linear)
//
Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
{
    const Int_t n = fGraph->GetN();
    const Double_t *x = fGraph->GetX();
    const Double_t *y = fGraph->GetY();

    Int_t i = 0;
    while (i<n && y[i]<accgam)
        i++;

    if (i==0 || i==n)
        return 0;

    if (i==n-1)
        i--;

    const Double_t x1 = y[i-1];
    const Double_t y1 = x[i-1];

    const Double_t x2 = y[i];
    const Double_t y2 = x[i];

    return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
}

// --------------------------------------------------------------------------
//
// Print the corresponding Gammas Acceptance for a hadron acceptance of
// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
// acceptance and the corresponding acceptances and hadroness value is
// printed, together with the maximum Q-factor.
//
void MHHadronness::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Hadronness histograms:" << endl;
    *fLog << "---------------------" << endl;

    if (fGraph->GetN()==0)
    {
        *fLog << " <No Entries>" << endl;
        return;
    }

    *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
    *fLog << "Acc Gammas @  1% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.01)*100) << "%" << endl;
    *fLog << "Acc Gammas @  2% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.02)*100) << "%" << endl;
    *fLog << "Acc Gammas @  5% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.05)*100) << "%" << endl;
    *fLog << "Acc Gammas @ 10% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.1)*100) << "%" << endl;
    *fLog << "Acc Gammas @ 20% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.2)*100) << "%" << endl;
    *fLog << "Acc Gammas @ 30% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.3)*100) << "%" << endl;
    *fLog << "Acc Gammas @ 40% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.4)*100) << "%" << endl;
    *fLog << "Acc Gammas @ 50% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.5)*100) << "%" << endl;

    const Int_t h = fMinDist->GetMaximumBin();
    *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
    *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
    *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;

    *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
    *fLog << "  Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
    *fLog << endl;
}

// --------------------------------------------------------------------------
//
// Draw clone of all histograms. (For the Meaning see class description)
//
TObject *MHHadronness::DrawClone(Option_t *opt) const
{
    if (fGraph->GetN()==0)
        return NULL;

    TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
    c.Divide(2, 2);

    gROOT->SetSelectedPad(NULL);

    c.cd(1);
    gStyle->SetOptStat(10);
    Getghness()->DrawCopy();
    Getphness()->SetLineColor(kRed);
    Getphness()->DrawCopy("same");
    c.cd(2);
    gStyle->SetOptStat(0);
    Getighness()->DrawCopy();
    Getiphness()->SetLineColor(kRed);
    Getiphness()->DrawCopy("same");
    fMinDist->SetLineColor(kBlue);
    fMinDist->DrawCopy("same");
    c.cd(3);
    //GetQfac()->DrawCopy();
    TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
    g2.SetBit(kCanDelete);
    gPad->Modified();
    gPad->Update();
    if (g2.GetHistogram())
    {
        g2.GetXaxis()->SetRangeUser(0, 1);
        g2.GetXaxis()->SetTitle("Acceptance Gammas");
        g2.GetYaxis()->SetTitle("Quality");
        g2.SetMarkerStyle(kFullDotSmall);
        g2.Draw("P");

        gPad->Modified();
        gPad->Update();
    }
    c.cd(4);
    gPad->Modified();
    gPad->Update();
    TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
    g.SetBit(kCanDelete);
    gPad->Modified();
    gPad->Update();
    if (g.GetHistogram())
    {
        g.GetXaxis()->SetRangeUser(0, 1);
        g.GetXaxis()->SetTitle("Acceptance Hadrons");
        g.GetYaxis()->SetTitle("Acceptance Gammas");
        g.SetMarkerStyle(kFullDotSmall);
        g.Draw("P");

        gPad->Modified();
        gPad->Update();
    }
    const Int_t h = fMinDist->GetMaximumBin();
    TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
                             fIntGhness->GetBinContent(h), kStar);
    m->SetMarkerColor(kBlue);
    m->SetBit(kCanDelete);
    m->Draw();
    return &c;
}

// --------------------------------------------------------------------------
//
// Draw all histograms. (For the Meaning see class description)
//
void MHHadronness::Draw(Option_t *)
{
   if (fGraph->GetN()==0)
        return;

   if (!gPad)
        MakeDefCanvas("Hadronness", fTitle);

    gPad->Divide(2, 2);

    gPad->cd(1);
    gStyle->SetOptStat(10);
    Getghness()->Draw();
    Getphness()->SetLineColor(kRed);
    Getphness()->Draw("same");
    gPad->cd(2);
    gStyle->SetOptStat(0);
    Getighness()->Draw();
    Getiphness()->SetLineColor(kRed);
    Getiphness()->Draw("same");
    fMinDist->SetLineColor(kBlue);
    fMinDist->Draw("same");
    gPad->cd(3);
    //GetQfac()->Draw();
    fQfac->Draw("A*");
    gPad->Modified();
    gPad->Update();
    if (fQfac->GetHistogram())
    {
        fQfac->GetXaxis()->SetRangeUser(0, 1);
        fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
        fQfac->GetYaxis()->SetTitle("Quality");
        fQfac->SetMarkerStyle(kFullDotSmall);
        fQfac->Draw("P");

        gPad->Modified();
        gPad->Update();
    }
    gPad->cd(4);
    gPad->Modified();
    gPad->Update();
    fGraph->Draw("AC");
    gPad->Modified();
    gPad->Update();
    if (fGraph->GetHistogram())
    {
        fGraph->GetXaxis()->SetRangeUser(0, 1);
        fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
        fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
        fGraph->SetMarkerStyle(kFullDotSmall);
        fGraph->Draw("P");
        gPad->Modified();
        gPad->Update();
    }
    const Int_t h = fMinDist->GetMaximumBin();
    TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
                             fIntGhness->GetBinContent(h), kStar);
    m->SetMarkerColor(kBlue);
    m->SetBit(kCanDelete);
    m->Draw();
}