source: trunk/MagicSoft/Mars/mhist/MHHadroness.cc@ 1474

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

/////////////////////////////////////////////////////////////////////////////
//
// MHHadroness
//
// This is histogram is a way to evaluate the quality of a gamma/hadron
// seperation method. It is filled from a MHadroness 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 vs. the hadroness
//    - red:   acceptance of non gammas 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 "MHHadroness.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 "MHadroness.h"

#include "MMcEvt.hxx"

ClassImp(MHHadroness);

// --------------------------------------------------------------------------
//
// Setup histograms, nbins is the number of bins used for the evaluation.
// The default is 100 bins.
//
MHHadroness::MHHadroness(Int_t nbins, const char *name, const char *title)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHHadroness";
    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("Hadroness");
    fPhness->SetXTitle("Hadroness");
    fGhness->SetYTitle("Counts");
    fPhness->SetYTitle("Counts");
    
    fIntGhness = new TH1D("AccGammas",  "Acceptance", nbins, 0, 1);
    fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
    fIntGhness->SetXTitle("Hadroness");
    fIntPhness->SetXTitle("Hadroness");
    fIntGhness->SetYTitle("Acceptance");
    fIntPhness->SetYTitle("Acceptance");

    fQfac = new TH1D("Qfac", "Naive Quality factor", nbins, 0, 1);
    fQfac->SetXTitle("Hadroness");
    fQfac->SetYTitle("Quality");

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

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

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

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

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

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

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

     fHist.Sumw2();
     */

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the Hadroness from a MHadroness container into the corresponding
// histogram dependant on the particle id.
//
Bool_t MHHadroness::Fill(const MParContainer *par)
{
    if (fMcEvt->GetPartId()==kGAMMA)
        fGhness->Fill(fHadroness->GetHadroness());
    else
        fPhness->Fill(fHadroness->GetHadroness());

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// 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 MHHadroness::Finalize()
{
    Int_t n = fGhness->GetNbinsX();

    fGraph->Set(n);

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

    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)));
        fQfac->SetBinContent(i, ig/sqrt(ip));
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Search the corresponding points for the given hadron acceptance (acchad)
// and interpolate the to pointsd (linear)
//
Double_t MHHadroness::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;
}

// --------------------------------------------------------------------------
//
// 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 MHHadroness::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Hadroness histograms:" << endl;
    *fLog << "---------------------" << endl;
    *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
    *fLog << "Acc Gammas at  1% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.01)*100) << "%" << endl;
    *fLog << "Acc Gammas at  2% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.02)*100) << "%" << endl;
    *fLog << "Acc Gammas at  5% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.05)*100) << "%" << endl;
    *fLog << "Acc Gammas at 10% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.1)*100) << "%" << endl;
    *fLog << "Acc Gammas at 20% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.2)*100) << "%" << endl;
    *fLog << "Acc Gammas at 30% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.3)*100) << "%" << endl;
    *fLog << "Acc Gammas at 40% Hadron-acc: " << Form("%3.0f", GetGammaAcceptance(0.4)*100) << "%" << endl;
    *fLog << "Acc Gammas at 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;
    const Int_t q = GetQfac()->GetMaximumBin();
    *fLog << "Maximum Q-Factor: " << GetQfac()->GetMaximum() << " @ H=";
    *fLog << GetQfac()->GetBinCenter(q) << endl;;
    *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(q)*100) << "%, ";
    *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(q)*100) << "%" << endl;
    *fLog << endl;
}

// --------------------------------------------------------------------------
//
// Draw clone of all histograms. (For the Meaning see class description)
//
TObject *MHHadroness::DrawClone(Option_t *opt) const
{
    TCanvas &c = *MakeDefCanvas("Hadroness", fTitle);
    c.Divide(2, 2);

    gROOT->SetSelectedPad(NULL);

    c.cd(1);
    gStyle->SetOptStat(10);
    Getghness()->DrawCopy();
    Getphness()->SetLineColor(kRed);
    Getphness()->DrawCopy("same");
    c.cd(4);
    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();
    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();

    return &c;
}

// --------------------------------------------------------------------------
//
// Draw all histograms. (For the Meaning see class description)
//
void MHHadroness::Draw(Option_t *)
{
    if (!gPad)
        MakeDefCanvas("Hadroness", fTitle);

    gPad->Divide(2, 2);

    gPad->cd(1);
    gStyle->SetOptStat(10);
    Getghness()->Draw();
    Getphness()->SetLineColor(kRed);
    Getphness()->Draw("same");
    gPad->cd(4);
    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();
    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();
}