#include "TCanvas.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TH1D.h"
#include "TH2D.h"
#include "THStack.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TPaveText.h"
#include "TPie.h"
#include "TROOT.h"
#include "TString.h"
#include "TStyle.h"
#include "TSystem.h"

#include <algorithm>
#include <iostream>
#include <sstream>

using namespace std;

#include "scaleHists.C"

#include "/u/ekrebs/workspace/cbmroot/JUL13/dielectron/CbmMvdAnaNames.h"

const Int_t nofMinStationsReco = 3;

//TString sSaveDir("/u/ekrebs/docs/talks/Cbm_Dilepton_2014_02_21");
TString sInDir("./histos_jun14");
//TString sInDir("./");
TString sAna("reco");
TString sMvd("mvd3");
TString sField("1.0");
TString sFieldCutoff("");
TString sLvm("rho0");
TString sDecay("epem");
TString sEnergy("25gev");
TString sSignalType(Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));

const Int_t nofCuts = 7;
const Int_t cutTypes[nofCuts] = {
    CbmMvdAnaNames::kCutNone,
    CbmMvdAnaNames::kCutGamma,
    CbmMvdAnaNames::kCutChi2Primary,
    CbmMvdAnaNames::kCutMvdPrimary,
    CbmMvdAnaNames::kCutTT, CbmMvdAnaNames::kCutPt,
    CbmMvdAnaNames::kCutKfp
};
const Color_t cutTypesColors[nofCuts] = {
    kBlack,
    kGray+1,
    kRed,
    kMagenta,
    kGreen+2,
    kBlue,
    kViolet
};

const char* cutTypesNames[CbmMvdAnaNames::kNofCutTypes] = {
    "PID",
    "m_{#gamma} cut",
    "Primary #chi^{2} track cut",
    "1st MVD cut",
    "Track topology cut", "p_{t} cut",
    "KFParticle cut"
};

// THStack must only contain histograms of TH1 and must have been
// drawn already.
void setStackAxis(THStack *hStack) {
    TH1 *h = (TH1*)hStack->GetStack()->At(0);
    hStack->Draw("nostack hist");
    hStack->GetXaxis()->SetTitle    (h->GetXaxis()->GetTitle());
    hStack->GetXaxis()->SetTitleFont(h->GetXaxis()->GetTitleFont());
    hStack->GetXaxis()->SetTitleSize(h->GetXaxis()->GetTitleSize());
    hStack->GetXaxis()->SetLabelFont(h->GetXaxis()->GetLabelFont());
    hStack->GetXaxis()->SetLabelSize(h->GetXaxis()->GetLabelSize());

    //hStackTrueZCut[iSource]->GetYaxis()->SetTitleOffset(1.3);
    hStack->GetYaxis()->SetTitle    (h->GetYaxis()->GetTitle());
    hStack->GetYaxis()->SetTitleFont(h->GetYaxis()->GetTitleFont());
    hStack->GetYaxis()->SetTitleSize(h->GetYaxis()->GetTitleSize());
    hStack->GetYaxis()->SetLabelFont(h->GetYaxis()->GetLabelFont());
    hStack->GetYaxis()->SetLabelSize(h->GetYaxis()->GetLabelSize());
}


Int_t getNevents(TFile *file) {
    Int_t nEvents = 0;
    file->cd("/");
    TH1F* h_events = (TH1F*)gDirectory->Get("fh_reco_event_number;1");
    if(h_events) {
	nEvents = h_events->GetEntries();
    } else {
	cout<<"Event counter histogram not found."<<endl;
	nEvents = 0;
    }
    return nEvents;
}

void drawMcGammaConv(TFile *file) {
    if(!file->cd("GammaConversions")) {
        cout<<"Directory Gamma Conversions not found"<<endl;
        return;
    }

    TCanvas *cLowZ  = new TCanvas("c_gamma_conv_zr", "Gamma Conversions",800,800);
    TCanvas *cHighZ = new TCanvas("c_gamma_conv_zr_high", "Gamma Conversions",800,800);

    TPaveText *ptSts  = new TPaveText(0.,  330., 150.,  430.);
    ptSts->AddText("MVD, STS");
    TPaveText *ptRich = new TPaveText(200., 500., 300.,  600.);
    ptRich->AddText("RICH");
    TPaveText *ptTrd  = new TPaveText(550., 850., 650.,  950.);
    ptTrd->AddText("TRD");
    TPaveText *ptTof  = new TPaveText(950., 950,  1050., 1050.);
    ptTof->AddText("TOF");

    cLowZ->cd();

    TH1F *fh_low = (TH1F*)gDirectory->Get("fh_gamma_conv_zr;1");
    cLowZ->SetLogz();
    fh_low->SetTitle("Origin of #gamma-Conversions");
    fh_low->GetXaxis()->SetRangeUser(0., 120.);
    fh_low->GetXaxis()->SetNoExponent(kTRUE);
    fh_low->GetYaxis()->SetRangeUser(-150., 150.);
    fh_low->GetYaxis()->SetTitleOffset(1.35);
    fh_low->GetYaxis()->SetNoExponent(kTRUE);
    fh_low->Draw("colz");

    TPaveText *ptMvdLow  = new TPaveText(8.,  50., 23.,  70.);
    ptMvdLow->AddText("MVD");
    ptMvdLow->SetFillColor(kYellow-6);
    ptMvdLow->Draw();

    TPaveText *ptStsLow  = new TPaveText(50.,  50., 65.,  70.);
    ptStsLow->AddText("STS");
    ptStsLow->SetFillColor(kYellow-6);
    ptStsLow->Draw();


    cHighZ->cd();
    TH1F *fh_high = (TH1F*)gDirectory->Get("fh_gamma_conv_zr_high;1");
    cHighZ->SetLogz();
    fh_high->SetTitle("Origin of #gamma-Conversions");
    fh_high->GetXaxis()->SetRangeUser(0., 1200.);
    fh_high->GetXaxis()->SetNoExponent(kTRUE);
    fh_high->GetYaxis()->SetRangeUser(-1100., 1100.);
    fh_high->GetYaxis()->SetTitleOffset(1.35);
    fh_high->GetYaxis()->SetNoExponent(kTRUE);
    fh_high->Draw("colz");

    ptSts->Draw();
    ptRich->Draw();
    ptTrd->Draw();
    ptTof->Draw();


    file->cd("/");
}

void drawMcDetHistShort(TFile *file) {
    if(!file->cd("Detectors")) {
        cout<<"Directory Detector not found"<<endl;
        return;
    }

    TCanvas *cDetHistShortSignal = new TCanvas("cDetHistShortSignal", "Detector History",800,800);
    TCanvas *cDetHistShortBg     = new TCanvas("cDetHistShortSignalBg", "Detector History",800,800);
    TCanvas *cDetHistShortPi0    = new TCanvas("cDetHistShortSignalPi0", "Detector History",800,800);
    TCanvas *cDetHistShortGamma  = new TCanvas("cDetHistShortSignalGamma", "Detector History",800,800);

    cDetHistShortSignal->cd();
    TH1F *fh_ee_mc_det_short_signal = (TH1F*)gDirectory->Get("fh_ee_mc_det_short_signal;1");
    Double_t n = fh_ee_mc_det_short_signal->GetEntries();
    fh_ee_mc_det_short_signal->SetTitle("Signal");
    fh_ee_mc_det_short_signal->Scale(100. / n);
    fh_ee_mc_det_short_signal->GetXaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_signal->GetYaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_signal->GetYaxis()->LabelsOption("v");
    fh_ee_mc_det_short_signal->Draw("colztext");

    cDetHistShortBg->cd();
    TH1F *fh_ee_mc_det_short_bg = (TH1F*)gDirectory->Get("fh_ee_mc_det_short_bg;1");
    n = fh_ee_mc_det_short_bg->GetEntries();
    fh_ee_mc_det_short_bg->SetTitle("#eta");
    fh_ee_mc_det_short_bg->Scale(100. / n);
    fh_ee_mc_det_short_bg->GetXaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_bg->GetYaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_bg->Draw("colztext");

    cDetHistShortPi0->cd();
    TH1F *fh_ee_mc_det_short_pi0 = (TH1F*)gDirectory->Get("fh_ee_mc_det_short_pi0;1");
    n = fh_ee_mc_det_short_pi0->GetEntries();
    fh_ee_mc_det_short_pi0->SetTitle("#pi^{0}");
    fh_ee_mc_det_short_pi0->Scale(100. / n);
    fh_ee_mc_det_short_pi0->GetXaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_pi0->GetYaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_pi0->Draw("colztext");

    cDetHistShortGamma->cd();
    TH1F *fh_ee_mc_det_short_gamma = (TH1F*)gDirectory->Get("fh_ee_mc_det_short_gamma;1");
    n = fh_ee_mc_det_short_gamma->GetEntries();
    fh_ee_mc_det_short_gamma->SetTitle("#gamma-Conversion");
    fh_ee_mc_det_short_gamma->Scale(100. / n);
    fh_ee_mc_det_short_gamma->GetXaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_gamma->GetYaxis()->SetBinLabel(4, "Segment");
    fh_ee_mc_det_short_gamma->Draw("colztext");

    TString sHistPrefix("fh_ee_mc_det_short");
    if(sFieldCutoff.Length() > 0) {
	sHistPrefix.Append("_");
	sHistPrefix.Append(sFieldCutoff.Data());
    }

    //cDetHistShortSignal->SaveAs(Form("%s/%s_signal_%sfield_%s.png", sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cDetHistShortBg    ->SaveAs(Form("%s/%s_bg_%sfield_%s.png",     sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cDetHistShortGamma ->SaveAs(Form("%s/%s_gamma_%sfield_%s.png",  sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cDetHistShortPi0   ->SaveAs(Form("%s/%s_pi0_%sfield_%s.png",    sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));

    file->cd("/");
}

void drawMcTrackFragsMvdSts(TFile *file) {
    if(!file->cd("Detectors")) {
        cout<<"Directory Detector not found"<<endl;
        return;
    }

    const Color_t gainFillColor = kBlue;
    const Style_t gainFillStyle = 3004;
    const Float_t barWidth = 0.9;
    const Float_t barOffset = 0.1;

    TCanvas *cTrackFragMvdStsSignal = new TCanvas("cTrackFragMvdStsSignal", "Number of traversed MVD and STS stations for track fragments",800,800);
    TCanvas *cTrackFragMvdStsBg = new TCanvas("cTrackFragMvdStsBg", "Number of traversed MVD and STS stations for track fragments",800,800);
    TCanvas *cTrackFragMvdStsPi0 = new TCanvas("cTrackFragMvdStsPi0", "Number of traversed MVD and STS stations for track fragments",800,800);
    TCanvas *cTrackFragMvdStsGamma = new TCanvas("cTrackFragMvdStsGamma", "Number of traversed MVD and STS stations for track fragments",800,800);

    cTrackFragMvdStsSignal->cd();
    TH1F *fh_e_mc_det_nmvdsts_notreco_signal = (TH1F*)gDirectory->Get("fh_e_mc_det_nmvdsts_notreco_signal;1");
    Double_t n = fh_e_mc_det_nmvdsts_notreco_signal->GetEntries();
    fh_e_mc_det_nmvdsts_notreco_signal->SetTitle("Signal");
    fh_e_mc_det_nmvdsts_notreco_signal->GetXaxis()->SetTitle("Number of MVD+STS Stations");
    fh_e_mc_det_nmvdsts_notreco_signal->GetYaxis()->SetTitleOffset(1.35);
    //fh_e_mc_det_nmvdsts_notreco_signal->Scale(1. / (n));

    Int_t nBins = fh_e_mc_det_nmvdsts_notreco_signal->GetNbinsX()-1;
    TH1F *fh_nmvdsts_ts_signal = new TH1F("fh_nmvdsts_ts_signal", "", nBins-2, 2.5, nBins+.5);
    for(Int_t i = nofMinStationsReco; i < nBins; i++) {
	fh_nmvdsts_ts_signal->SetBinContent(i-nofMinStationsReco, fh_e_mc_det_nmvdsts_notreco_signal->GetBinContent(i));
    }
    fh_nmvdsts_ts_signal->SetFillColor(gainFillColor);
    fh_nmvdsts_ts_signal->SetFillStyle(gainFillStyle);
    fh_nmvdsts_ts_signal->SetBarWidth(barWidth);
    fh_nmvdsts_ts_signal->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_signal->SetLineColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_signal->SetFillColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_signal->SetFillStyle(0);
    fh_e_mc_det_nmvdsts_notreco_signal->SetBarWidth(barWidth);
    fh_e_mc_det_nmvdsts_notreco_signal->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_signal->Draw("bar");
    fh_nmvdsts_ts_signal->Draw("Bsames");


    cTrackFragMvdStsBg->cd();
    TH1F *fh_e_mc_det_nmvdsts_notreco_bg = (TH1F*)gDirectory->Get("fh_e_mc_det_nmvdsts_notreco_bg;1");
    n = fh_e_mc_det_nmvdsts_notreco_bg->GetEntries();
    fh_e_mc_det_nmvdsts_notreco_bg->SetTitle("#eta");
    fh_e_mc_det_nmvdsts_notreco_bg->GetXaxis()->SetTitle("Number of MVD+STS Stations");
    fh_e_mc_det_nmvdsts_notreco_bg->GetYaxis()->SetTitleOffset(1.35);
    //fh_e_mc_det_nmvdsts_notreco_bg->Scale(1. / (n));

    nBins = fh_e_mc_det_nmvdsts_notreco_bg->GetNbinsX()-1;
    TH1F *fh_nmvdsts_ts_bg = new TH1F("fh_nmvdsts_ts_bg", "", nBins-2, 2.5, nBins+0.5);
    for(Int_t i = nofMinStationsReco; i < nBins; i++) {
	fh_nmvdsts_ts_bg->SetBinContent(i-nofMinStationsReco, fh_e_mc_det_nmvdsts_notreco_bg->GetBinContent(i));
    }
    fh_nmvdsts_ts_bg->SetFillColor(gainFillColor);
    fh_nmvdsts_ts_bg->SetFillStyle(gainFillStyle);
    fh_nmvdsts_ts_bg->SetBarWidth(barWidth);
    fh_nmvdsts_ts_bg->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_bg->SetLineColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_bg->SetFillColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_bg->SetFillStyle(0);
    fh_e_mc_det_nmvdsts_notreco_bg->SetBarWidth(barWidth);
    fh_e_mc_det_nmvdsts_notreco_bg->SetBarOffset(barOffset);


    fh_e_mc_det_nmvdsts_notreco_bg->Draw("bar");
    fh_nmvdsts_ts_bg->Draw("Bsames");



    cTrackFragMvdStsPi0->cd();
    TH1F *fh_e_mc_det_nmvdsts_notreco_pi0 = (TH1F*)gDirectory->Get("fh_e_mc_det_nmvdsts_notreco_pi0;1");
    n = fh_e_mc_det_nmvdsts_notreco_pi0->GetEntries();
    fh_e_mc_det_nmvdsts_notreco_pi0->SetTitle("#pi^{0}-Dalitz");
    fh_e_mc_det_nmvdsts_notreco_pi0->GetXaxis()->SetTitle("Number of MVD+STS Stations");
    fh_e_mc_det_nmvdsts_notreco_pi0->GetYaxis()->SetTitleOffset(1.35);
    //fh_e_mc_det_nmvdsts_notreco_pi0->Scale(1. / (n));

    TH1F *fh_nmvdsts_ts_pi0 = new TH1F("fh_nmvdsts_ts_pi0", "", nBins-2, 2.5, nBins+.5);
    for(Int_t i = nofMinStationsReco; i < nBins; i++) {
	fh_nmvdsts_ts_pi0->SetBinContent(i-nofMinStationsReco, fh_e_mc_det_nmvdsts_notreco_pi0->GetBinContent(i));
    }
    fh_nmvdsts_ts_pi0->SetFillColor(gainFillColor);
    fh_nmvdsts_ts_pi0->SetFillStyle(gainFillStyle);
    fh_nmvdsts_ts_pi0->SetBarWidth(barWidth);
    fh_nmvdsts_ts_pi0->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_pi0->SetLineColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_pi0->SetFillColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_pi0->SetFillStyle(0);
    fh_e_mc_det_nmvdsts_notreco_pi0->SetBarWidth(barWidth);
    fh_e_mc_det_nmvdsts_notreco_pi0->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_pi0->Draw("bar");
    fh_nmvdsts_ts_pi0->Draw("Bsames");


    cTrackFragMvdStsGamma->cd();
    TH1F *fh_e_mc_det_nmvdsts_notreco_gamma = (TH1F*)gDirectory->Get("fh_e_mc_det_nmvdsts_notreco_gamma;1");
    n = fh_e_mc_det_nmvdsts_notreco_gamma->GetEntries();
    fh_e_mc_det_nmvdsts_notreco_gamma->SetTitle("#gamma-Conversions");
    fh_e_mc_det_nmvdsts_notreco_gamma->GetXaxis()->SetTitle("Number of MVD+STS Stations");
    fh_e_mc_det_nmvdsts_notreco_gamma->GetYaxis()->SetTitleOffset(1.35);
    //fh_e_mc_det_nmvdsts_notreco_gamma->Scale(1. / (n));

    TH1F *fh_nmvdsts_ts_gamma = new TH1F("fh_nmvdsts_ts_gamma", "", nBins-2, 2.5, nBins+.5);
    for(Int_t i = nofMinStationsReco; i < nBins; i++) {
	fh_nmvdsts_ts_gamma->SetBinContent(i-nofMinStationsReco, fh_e_mc_det_nmvdsts_notreco_gamma->GetBinContent(i));
    }
    fh_nmvdsts_ts_gamma->SetFillColor(gainFillColor);
    fh_nmvdsts_ts_gamma->SetFillStyle(gainFillStyle);
    fh_nmvdsts_ts_gamma->SetBarWidth(barWidth);
    fh_nmvdsts_ts_gamma->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_gamma->SetLineColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_gamma->SetFillColor(kWhite);
    fh_e_mc_det_nmvdsts_notreco_gamma->SetFillStyle(0);
    fh_e_mc_det_nmvdsts_notreco_gamma->SetBarWidth(barWidth);
    fh_e_mc_det_nmvdsts_notreco_gamma->SetBarOffset(barOffset);

    fh_e_mc_det_nmvdsts_notreco_gamma->Draw("bar");
    fh_nmvdsts_ts_gamma->Draw("Bsames");

    TString sHistPrefix("fh_e_mc_det_nmvdsts_notreco");
    if(sFieldCutoff.Length() > 0) {
	sHistPrefix.Append("_");
	sHistPrefix.Append(sFieldCutoff.Data());
    }

    //cTrackFragMvdStsSignal->SaveAs(Form("%s/%s_signal_%sfield_%s.png", sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cTrackFragMvdStsBg    ->SaveAs(Form("%s/%s_bg_%sfield_%s.png",     sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cTrackFragMvdStsGamma ->SaveAs(Form("%s/%s_gamma_%sfield_%s.png",  sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));
    //cTrackFragMvdStsPi0   ->SaveAs(Form("%s/%s_pi0_%sfield_%s.png",    sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));


    file->cd("/");
}

void drawMcMomAccepted(TFile *file) {
    if(!file->cd("MomentumDistribution")) {
        cout<<"Directory MomentumDistribution not found"<<endl;
        return;
    }

    TCanvas *c = new TCanvas("momAccepted", "Momentrum Distribution", 800, 800);
    c->cd();

    TH1F *fh_e_mc_mom_accepted_signal = (TH1F*)gDirectory->Get("fh_e_mc_mom_Accepted_signal;1");
    fh_e_mc_mom_accepted_signal->SetLineColor(kBlue);
    fh_e_mc_mom_accepted_signal->GetXaxis()->SetRangeUser(0., 2.0);
    fh_e_mc_mom_accepted_signal->GetXaxis()->SetTitleOffset(1.35);

    TH1F *fh_e_mc_mom_accepted_pi0 = (TH1F*)gDirectory->Get("fh_e_mc_mom_Accepted_pi0;1");
    fh_e_mc_mom_accepted_pi0->SetLineColor(kRed);
    fh_e_mc_mom_accepted_pi0->SetTitle("Momentum distribution of e^{+}/e^{-}.");

    TH1F *fh_e_mc_mom_accepted_gamma = (TH1F*)gDirectory->Get("fh_e_mc_mom_Accepted_gamma;1");
    fh_e_mc_mom_accepted_gamma->SetLineColor(kMagenta);

    fh_e_mc_mom_accepted_pi0->Draw();
    fh_e_mc_mom_accepted_gamma->Draw("same");
    fh_e_mc_mom_accepted_signal->Draw("same");

    //gPad->Update(); // Needed to get stats window.
    //TPaveStats *stats = (TPaveStats*)fh_e_mc_mom_accepted_pi0->FindObject("stats");
    //stats->SetOptStat(0000);

    TLegend* legend = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
    TLegendEntry* legentry = legend->AddEntry(fh_e_mc_mom_accepted_signal,"#rho^{0} #rightarrow e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_e_mc_mom_accepted_pi0,"#pi^{0} #rightarrow #gamma e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_e_mc_mom_accepted_gamma,"#gamma #rightarrow e^{+} e^{-}","lpf");
    legend->SetFillColor(kWhite);
    legend->Draw();


    file->cd("/");
}

void drawMcThetaAccepted(TFile *file) {
    if(!file->cd("OpeningAngle")) {
        cout<<"Directory OpeningAngle not found"<<endl;
        return;
    }

    TCanvas *c = new TCanvas("thetaAccepted", "Opening Angle", 800, 800);
    c->cd()->SetLogy();

    TH1F *fh_ee_mc_angle_accepted_signal = (TH1F*)gDirectory->Get("fh_ee_mc_angle_Accepted_signal;1");
    fh_ee_mc_angle_accepted_signal->SetLineColor(kBlue);

    TH1F *fh_ee_mc_angle_accepted_pi0 = (TH1F*)gDirectory->Get("fh_ee_mc_angle_Accepted_pi0;1");
    fh_ee_mc_angle_accepted_pi0->SetLineColor(kRed);
    fh_ee_mc_angle_accepted_pi0->SetTitle("Opening angles of e^{+}/e^{-}.");

    TH1F *fh_ee_mc_angle_accepted_gamma = (TH1F*)gDirectory->Get("fh_ee_mc_angle_Accepted_gamma;1");
    fh_ee_mc_angle_accepted_gamma->SetLineColor(kMagenta);

    fh_ee_mc_angle_accepted_pi0->Draw();
    fh_ee_mc_angle_accepted_gamma->Draw("same");
    fh_ee_mc_angle_accepted_signal->Draw("same");

    //gPad->Update(); // Needed to get stats window.
    //TPaveStats *stats = (TPaveStats*)fh_ee_mc_angke_accepted_pi0->FindObject("stats");
    //stats->SetOptStat(0000);

    TLegend* legend = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
    TLegendEntry* legentry = legend->AddEntry(fh_ee_mc_angle_accepted_signal,"#rho^{0} #rightarrow e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_ee_mc_angle_accepted_pi0,"#pi^{0} #rightarrow #gamma e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_ee_mc_angle_accepted_gamma,"#gamma #rightarrow e^{+} e^{-}","lpf");
    legend->SetFillColor(kWhite);
    legend->Draw();


    file->cd("/");
}

void drawMcHitTopology(TFile *file) {
    if(!file->cd("PairDistanceInMVD")) {
        cout<<"Directory PairDistanceInMVD not found"<<endl;
        return;
    }

    TH1F *fh_mc_dist_pair_signal = (TH1F*)gDirectory->Get("fh_ee_mc_mvd1_dist_signal;1");
    TH1F *fh_mc_dist_pair_gamma  = (TH1F*)gDirectory->Get("fh_ee_mc_mvd1_dist_gamma;1");

    if(!file->cd("Distance to NN in MVD")) {
        cout<<"Directory Distance to NN in MVD not found"<<endl;
        return;
    }

    TCanvas *cDistNN = new TCanvas("McHitTopology", "MC Hit Topology Cut", 800, 800);
    cDistNN->cd();

    TH1F *fh_mc_dist_nn_signal = (TH1F*)gDirectory->Get("fh_e_mc_mvd1_dist_nn_signal;1");
    TH1F *fh_mc_dist_nn_gamma  = (TH1F*)gDirectory->Get("fh_e_mc_mvd1_dist_nn_gamma;1");

    fh_mc_dist_pair_signal->SetLineColor(kOrange);
    fh_mc_dist_pair_gamma->SetLineColor(kBlue);

    fh_mc_dist_nn_signal->SetLineColor(kRed);
    fh_mc_dist_nn_gamma->SetLineColor(kBlack);

    fh_mc_dist_pair_gamma->SetTitle("");
    fh_mc_dist_pair_gamma->Draw();
    fh_mc_dist_pair_signal->Draw("sames");
    fh_mc_dist_nn_signal->Draw("sames");
    fh_mc_dist_nn_gamma->Draw("sames");

    TLegend* legend = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
    TLegendEntry* legentry = NULL;
    legentry = legend->AddEntry(fh_mc_dist_pair_gamma,"#gamma #rightarrow e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_mc_dist_nn_gamma,"e_{#gamma} to NN", "lpf");
    legentry = legend->AddEntry(fh_mc_dist_pair_signal,"#rho^{0} #rightarrow e^{+} e^{-}","lpf");
    legentry = legend->AddEntry(fh_mc_dist_nn_signal,"e_{#rho^{0}} to NN","lpf");
    legend->SetFillColor(kWhite);
    legend->Draw();

    TString sHistPrefix("fh_ee_dist");
    //cDistNN->SaveAs(Form("%s/%s_signal_%sfield_%s.png", sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));

    file->cd("/");
}


void drawTrackTopology(TFile *file) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("TrackTopologyCut")) {
        cout<<"Directory HitTopologyCut not found"<<endl;
        return;
    }


    TString histNamePrefix("fh_cut_track_topology");

    TString histName = histNamePrefix + "_signal";
    TH2F *fh_tt_signal = (TH2F*)gDirectory->Get(histName.Data());
    if(!fh_tt_signal) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = histNamePrefix + "_gamma";
    TH2F *fh_tt_gamma = (TH2F*)gDirectory->Get(histName.Data());
    if(!fh_tt_gamma) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = histNamePrefix + "_pi0";
    TH2F *fh_tt_pi0 = (TH2F*)gDirectory->Get(histName.Data());
    if(!fh_tt_pi0) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = histNamePrefix + "_bg";
    TH2F *fh_tt_other = (TH2F*)gDirectory->Get(histName.Data());
    if(!fh_tt_other) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    TH2F *fh_tt_bg(fh_tt_gamma);
    fh_tt_bg->Add(fh_tt_pi0, fh_tt_gamma);
    fh_tt_bg->Add(fh_tt_bg);

    scaleHist2DSignal(fh_tt_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));
    scaleHist2D(fh_tt_pi0, nEvents);
    scaleHist2D(fh_tt_bg, nEvents);

    TCanvas *cTopoSignal = new TCanvas("cTopoSignal", "Track Topology Signal", 800, 800);
    cTopoSignal->cd();
    //fh_tt_signal->GetYaxis()->SetTitleOffset(1.15);
    fh_tt_signal->GetXaxis()->SetTitle("#sqrt{p_{e} p_{reco}} [#frac{GeV}{c}]");
    fh_tt_signal->Draw("colz");

    TCanvas *cTopoPi0 = new TCanvas("cTopoPi0", "Track Topology Pi0", 800, 800);
    cTopoPi0->cd();
    //fh_tt_pi0->GetYaxis()->SetTitleOffset(1.15);
    fh_tt_pi0->GetXaxis()->SetTitle("#sqrt{p_{e} p_{reco}} [#frac{GeV}{c}]");
    fh_tt_pi0->Draw("colz");

    TCanvas *cTopoBg = new TCanvas("cTopoBg", "Track Topology Background", 800, 800);
    cTopoBg->cd();
    //fh_tt_bg->GetYaxis()->SetTitleOffset(1.15);
    fh_tt_bg->GetXaxis()->SetTitle("#sqrt{p_{e} p_{reco}} [#frac{GeV}{c}]");
    fh_tt_bg->Draw("colz");

    file->cd("/");
}



void drawHitTopology(TFile *file) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("HitTopologyCut")) {
        cout<<"Directory HitTopologyCut not found"<<endl;
        return;
    }

    TH1F *fh_dist_pair_signal = (TH1F*)gDirectory->Get("fh_e_mvd_dist_pair_signal;1");
    if(!fh_dist_pair_signal) {
	cout<<"Histogram fh_e_mvd_dist_pair_signal;1 not found."<<endl;
        return;
    }
    TH1F *fh_dist_pair_gamma  = (TH1F*)gDirectory->Get("fh_e_mvd_dist_pair_gamma;1");
    if(!fh_dist_pair_gamma) {
	cout<<"Histogram fh_e_mvd_dist_pair_gamma;1 not found."<<endl;
        return;
    }

    TH1F *fh_dist_nn_signal = (TH1F*)gDirectory->Get("fh_e_mvd_dist_nn_signal;1");
    if(!fh_dist_nn_signal) {
	cout<<"Histogram fh_e_mvd_dist_nn_signal;1 not found."<<endl;
        return;
    }
    TH1F *fh_dist_nn_gamma  = (TH1F*)gDirectory->Get("fh_e_mvd_dist_nn_gamma;1");
    if(!fh_dist_nn_gamma) {
	cout<<"Histogram fh_e_mvd_dist_nn_gamma;1 not found."<<endl;
        return;
    }

    fh_dist_pair_signal->SetLineColor(kOrange);
    fh_dist_pair_gamma->SetLineColor(kBlue);

    fh_dist_nn_signal->SetLineColor(kRed);
    fh_dist_nn_gamma->SetLineColor(kBlack);

    TCanvas *cDistNN = new TCanvas("HitTopology", "Hit Topology Cut", 800, 800);
    cDistNN->SetLogy();

    scaleHistSignal(fh_dist_pair_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));
    scaleHist(fh_dist_pair_gamma, nEvents);

    scaleHistSignal(fh_dist_nn_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));
    scaleHist(fh_dist_nn_gamma, nEvents);

    fh_dist_pair_gamma->SetTitle("");
    fh_dist_nn_gamma->GetYaxis()->SetRangeUser(1.e-5, 5.);
    fh_dist_nn_gamma->GetYaxis()->SetTitle("dN/d_{MVD} [1/cm]");
    fh_dist_nn_gamma->GetYaxis()->SetTitleOffset(1.2);
    fh_dist_nn_gamma->Draw("hist");
    fh_dist_nn_signal->Draw("hist sames");
    //fh_dist_pair_gamma->Draw("hist sames");
    //fh_dist_pair_signal->Draw("hist sames");

    TLegend* legend = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
    TLegendEntry* legentry = NULL;
    legentry = legend->AddEntry(fh_dist_nn_gamma,"e_{#gamma} to NN", "lpf");
    legentry = legend->AddEntry(fh_dist_nn_signal,"e_{#rho^{0}} to NN","lpf");
    //legentry = legend->AddEntry(fh_dist_pair_gamma,"e_{#gamma} Pair","lpf");
    //legentry = legend->AddEntry(fh_dist_pair_signal,"e_{#rho^{0}} Pair","lpf");
    legend->SetFillColor(kWhite);
    legend->Draw();

    TString sHistPrefix("fh_ee_dist");
    //cDistNN->SaveAs(Form("%s/%s_signal_%sfield_%s.png", sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));


    TH2F *fh_dist_nn_mom_signal = (TH2F*)gDirectory->Get("fh_e_mvd_dist_nn_mom_signal;1");
    if(!fh_dist_nn_mom_signal) {
	cout<<"Histogram fh_e_mvd_dist_nn_mom_signal;1 not found."<<endl;
        return;
    }

    TH2F *fh_dist_nn_mom_gamma = (TH2F*)gDirectory->Get("fh_e_mvd_dist_nn_mom_gamma;1");
    if(!fh_dist_nn_mom_gamma) {
	cout<<"Histogram fh_e_mvd_dist_nn_mom_gamma;1 not found."<<endl;
        return;
    }

    TH2F *fh_dist_nn_mom_pi0 = (TH2F*)gDirectory->Get("fh_e_mvd_dist_nn_mom_pi0;1");
    if(!fh_dist_nn_mom_pi0) {
	cout<<"Histogram fh_e_mvd_dist_nn_mom_pi0;1 not found."<<endl;
        return;
    }

    TH2F *fh_dist_nn_mom_bg(fh_dist_nn_mom_gamma);
    fh_dist_nn_mom_bg->Add(fh_dist_nn_mom_pi0, fh_dist_nn_mom_gamma);

    scaleHist2DSignal(fh_dist_nn_mom_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));
    scaleHist2D(fh_dist_nn_mom_bg, nEvents);

    TCanvas *cHitTopoSignal = new TCanvas("cHitTopoSignal", "Hit Topology Signal", 800, 800);
    cHitTopoSignal->cd();
    fh_dist_nn_mom_signal->GetYaxis()->SetTitleOffset(1.15);
    fh_dist_nn_mom_signal->Draw("colz");

    TCanvas *cHitTopoBg = new TCanvas("cHitTopoBg", "Hit Topology Background", 800, 800);
    cHitTopoBg->cd();
    fh_dist_nn_mom_bg->GetYaxis()->SetTitleOffset(1.15);
    fh_dist_nn_mom_bg->Draw("colz");

    file->cd("/");
}

void drawHitMismatches(TFile *file) {
    Int_t nEvents = getNevents(file);

    TString dirName = "HitsInStsTrack";
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
	return;
    }

    TString histName("");

    TH1F *fh_mismatch[CbmMvdAnaNames::kNofSourceTypes];
    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        histName = Form("fh_e_mcid_mvd_hits_mvdstsmatch_%s",
                        CbmMvdAnaNames::kSourceTypesNames[iSource]);
        fh_mismatch[iSource] = (TH1F*)gDirectory->Get(histName.Data());
        if(!fh_mismatch[iSource]) {
            cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
            continue;
        }
    }

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        TH1F *h = fh_mismatch[iSource];
        if(!h) continue;
        TCanvas *c = new TCanvas(Form("MVD mismatches %s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                 Form("c_mvd_mismatch_%s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                 800, 800);
        c->cd();
        for(Int_t iBin = 1; iBin <= h->GetXaxis()->GetNbins(); iBin++) {
            if(iBin % 2 == 1) {
                Int_t mvd = 1 +  (iBin - 1) / 2;
                h->GetXaxis()->SetBinLabel(iBin, Form("Same MVD %i", mvd));
            }
        }
        h->GetYaxis()->SetTitle("Yield per event");
        h->GetXaxis()->SetTitle("");
        h->Scale(1./(Double_t)nEvents);
        h->Draw("");
    }
}

void drawTrackTypes(TFile *file, Int_t cut, Bool_t bDetailedSegs) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("TrackTypes")) {
        cout<<"Directory TrackTypes not found"<<endl;
        return;
    }

    TH2F* fh_ee_tracks[CbmMvdAnaNames::kNofSourceTypes];
    TH1F* fh_tracks[CbmMvdAnaNames::kNofSourceTypes];

    TString sDetSeg("");
    if(bDetailedSegs) {
        sDetSeg = "det_";
    }

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        fh_ee_tracks[iSource] = (TH2F*)gDirectory->Get(Form("fh_ee_track_%stypes_%s_%s;1",
                                                            sDetSeg.Data(),
                                                            CbmMvdAnaNames::kSourceTypesNames[iSource],
                                                            CbmMvdAnaNames::kCutTypesNames[cut]));

        fh_tracks[iSource] = new TH1F(Form("fh_%s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                      CbmMvdAnaNames::kSourceTypesNames[iSource],
                                      //CbmMvdAnaNames::kNofTrackDetTypesReco,
                                      //0.5, CbmMvdAnaNames::kNofTrackDetTypesReco+0.5);
                                      fh_ee_tracks[iSource]->GetNbinsX(),
                                      0.5, fh_ee_tracks[iSource]->GetNbinsX()+0.5);
        for(Int_t iBin = 1; iBin <= fh_tracks[iSource]->GetNbinsX(); iBin++) {
            //fh_tracks[iSource]->GetXaxis()->SetBinLabel(iBin, CbmMvdAnaNames::kTrackDetTypesNamesReco[iBin-1]);
            fh_tracks[iSource]->GetXaxis()->SetBinLabel(iBin, fh_ee_tracks[iSource]->GetXaxis()->GetBinLabel(iBin));
            fh_tracks[iSource]->SetBinContent(iBin, 0.);
        }
    }

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        TH2F *fh_ee = fh_ee_tracks[iSource];
	cout<<"Source: "<<CbmMvdAnaNames::kSourceTypesNames[iSource]<<endl;
        Float_t nIdentified = 0;
        if(!fh_ee) continue;
        for(Int_t iBin = 1; iBin <= fh_ee->GetNbinsX(); iBin++) {
            Float_t sum = fh_ee->GetBinContent(iBin) +
		//fh_ee->GetBinContent(iBin, CbmMvdAnaNames::kStsRich-1);
		fh_ee->GetBinContent(iBin, fh_ee->GetNbinsY());
            if(iBin != fh_ee->GetNbinsX()) {
                //sum += fh_ee->GetBinContent(CbmMvdAnaNames::kStsRich-1, iBin);
                sum += fh_ee->GetBinContent(fh_ee->GetNbinsX(), iBin);
	    }
            nIdentified += sum;
	    sum = sum / (Float_t)nEvents;
	    fh_tracks[iSource]->SetBinContent(iBin, sum);
	    if(bDetailedSegs) {
		cout<<"Track type "<<CbmMvdAnaNames::kTrackDetTypesNames[iBin]
		    <<": "<<sum<<endl;
	    } else {
		cout<<"Track type "<<CbmMvdAnaNames::kTrackDetTypesNames[iBin]
		    <<": "<<sum<<endl;
	    }
        }
	nIdentified = nIdentified / (Float_t)nEvents;
	cout<<"Total number of identified electrons: "<<nIdentified<<endl;
    }

    TCanvas *cTracks = new TCanvas("Tracks", "Tracks", 800, 800);
    cTracks->Divide(2,2);

    for(Int_t i = 0; i < CbmMvdAnaNames::kNofSourceTypes; i++) {
        cTracks->cd(i+1);
        fh_tracks[i]->Draw();
    }

    TString sHistPrefix("");
    if(bDetailedSegs) {
	sHistPrefix = "fh_ee_track_det_types";
    } else {
	sHistPrefix = "fh_ee_track_types";
    }

    //cTracks->SaveAs(Form("%s/%s_%sfield_%s.png", sSaveDir.Data(), sHistPrefix.Data(), sField.Data(), sMvd.Data()));

    TCanvas *cPieTracks = new TCanvas("PieTracks", "PieTracks", 600, 600);
    cPieTracks->cd();

    TH1F *h = fh_tracks[CbmMvdAnaNames::kSourcePi0];
    const Int_t n = CbmMvdAnaNames::kNofTrackTypes;
    Float_t val[n];
    const char *lab[n];

    val[0] = h->GetBinContent(3);
    val[1] = h->GetBinContent(4);
    val[2] = h->GetBinContent(1);
    val[3] = h->GetBinContent(2);

    lab[0] = h->GetXaxis()->GetBinLabel(3);
    lab[1] = h->GetXaxis()->GetBinLabel(4);
    lab[2] = h->GetXaxis()->GetBinLabel(1);
    lab[3] = h->GetXaxis()->GetBinLabel(2);

//    for(Int_t i = 0; i < n; i++) {
//        if(i == n-1) {
//            val[i] = h->GetBinContent(1);
//            lab[i] = h->GetXaxis()->GetBinLabel(1);
//        } else {
//            val[i] = h->GetBinContent(i+2);
//            lab[i] = h->GetXaxis()->GetBinLabel(i+2);
//        }
//    }

    TString sPieTitle("");
    if(sMvd == "nomvd") {
        sPieTitle = "#pi^{0}-Dalitz, No MVD stations";
    }
    if(sMvd == "mvd2") {
        sPieTitle = "#pi^{0}-Dalitz, 2 MVD stations";
    }
    if(sMvd == "mvd3") {
        sPieTitle = "#pi^{0}-Dalitz, 3 MVD stations";
    }
    TPie *pie_tracks = new TPie(h);
    //TPie *pie_tracks = new TPie("pie", sPieTitle.Data(), n, &val[0], &col[0], &lab[0]);
    pie_tracks->SetAngularOffset(90.-((h->GetBinContent(1)+h->GetBinContent(2))/h->Integral())*360.);
    pie_tracks->SetTitle(sPieTitle.Data());
    //pie_tracks->SetAngularOffset(90.);
    pie_tracks->SetRadius(0.25);
    pie_tracks->SetEntryFillColor(0, kGray+1);
    pie_tracks->SetEntryFillColor(1, kRed-4);
    pie_tracks->SetEntryFillColor(2, kBlue);
    pie_tracks->SetEntryFillColor(3, kGreen+1);
//    pie_tracks->SetLabelsOffset(1.1);
    TLegend *pie_legend = pie_tracks->MakeLegend();
    pie_tracks->Draw("r");
    pie_legend->Draw("same");
}

void drawTrueBgCuts(TFile *file) {
    Int_t nEvents = getNevents(file);
    string dirName = "InvariantMass";
    if(!file->cd(dirName.c_str())) {
        cout<<"Directory "<<dirName.c_str()<<" not found."<<endl;
        return;
    }
    string histName ="";

    Color_t color_true_bg[CbmMvdAnaNames::kNofCutTypes];
    color_true_bg[CbmMvdAnaNames::kCutNone]        = kBlack;
    color_true_bg[CbmMvdAnaNames::kCutGamma]       = kGray+1;
    color_true_bg[CbmMvdAnaNames::kCutChi2Primary] = kRed;
    color_true_bg[CbmMvdAnaNames::kCutMvdPrimary]  = kMagenta;
    color_true_bg[CbmMvdAnaNames::kCutD]           = kGreen;
    color_true_bg[CbmMvdAnaNames::kCutHT]          = kCyan+2;
    color_true_bg[CbmMvdAnaNames::kCutTT]          = kGreen+2;
    color_true_bg[CbmMvdAnaNames::kCutPt]          = kBlue;
    color_true_bg[CbmMvdAnaNames::kCutKfp]         = kViolet;

    Bool_t bDrawCut[CbmMvdAnaNames::kNofCutTypes];
    bDrawCut[CbmMvdAnaNames::kCutNone]        = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutGamma]       = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutChi2Primary] = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutMvdPrimary]  = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutD]           = kFALSE;
    bDrawCut[CbmMvdAnaNames::kCutHT]          = kFALSE;
    //if(sMvd == "nomvd") {
    //    bDrawCut[CbmMvdAnaNames::kCutHT]      = kFALSE;
    //}
    bDrawCut[CbmMvdAnaNames::kCutTT]      = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutPt]      = kTRUE;
    bDrawCut[CbmMvdAnaNames::kCutKfp]     = kTRUE;

    const char* CutTypesNames[CbmMvdAnaNames::kNofCutTypes];
    CutTypesNames[CbmMvdAnaNames::kCutNone]        = "No cut";
    CutTypesNames[CbmMvdAnaNames::kCutGamma]       = "#gamma cut";
    CutTypesNames[CbmMvdAnaNames::kCutChi2Primary] = "Primary #chi^{2} track cut";
    CutTypesNames[CbmMvdAnaNames::kCutMvdPrimary]  = "MVD station cut";
    CutTypesNames[CbmMvdAnaNames::kCutD]           = "distance cut MVD";
    CutTypesNames[CbmMvdAnaNames::kCutHT]          = "Hit topology cut";
    CutTypesNames[CbmMvdAnaNames::kCutTT]          = "Track topology cut";
    CutTypesNames[CbmMvdAnaNames::kCutPt]          = "p_{t} cut";
    CutTypesNames[CbmMvdAnaNames::kCutKfp]         = "KFParticle cut";


    //TH1F* fh_signal = (TH1F*)gDirectory->Get("fh_ee_minv_Global Track_signal;1");
    //if(!fh_signal) {
    //    cout<<"Histogram fh_ee_minv_Global Track_signal;1 not found"<<endl;
    //    return;
    //}
    //fh_signal->SetLineColor(color_signal);
    //scaleHistSignal(fh_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));

    TH1F* fh_signal[CbmMvdAnaNames::kNofCutTypes];

    for(Int_t iCut = 0; iCut < CbmMvdAnaNames::kNofCutTypes; iCut++) {
        fh_signal[iCut] = (TH1F*)gDirectory->Get(Form("fh_ee_minv_signal_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[iCut]));
        if(!fh_signal[iCut]) {
	    cout<<"Histogram signal for cut "<<iCut<<" not found."<<endl;
            continue;
	}
	fh_signal[iCut]->SetLineColor(color_true_bg[iCut]);
        scaleHistSignal(fh_signal[iCut], nEvents, sSignalType.Data());
    }

    TH1F* fh_true_bg[CbmMvdAnaNames::kNofCutTypes];
    TH1F* fh_true_bg_signal[CbmMvdAnaNames::kNofCutTypes];

    for(Int_t iCut = 0; iCut < CbmMvdAnaNames::kNofCutTypes; iCut++) {
        fh_true_bg[iCut] = (TH1F*)gDirectory->Get(Form("fh_ee_minv_true_bg_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[iCut]));
        if(!fh_true_bg[iCut]) {
	    cout<<"Histogram true bg for cut "<<iCut<<" not found."<<endl;
            continue;
	}
	fh_true_bg[iCut]->SetLineColor(color_true_bg[iCut]);
        scaleHist(fh_true_bg[iCut], nEvents);

        fh_true_bg_signal[iCut] = (TH1F*)gDirectory->Get(Form("fh_ee_minv_true_bg_signal_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[iCut]));
        if(!fh_true_bg_signal[iCut]) {
	    cout<<"Histogram true bg_signal for cut "<<iCut<<" not found."<<endl;
            continue;
	}
	fh_true_bg_signal[iCut]->SetLineColor(color_true_bg[iCut]);
        scaleHistSignal(fh_true_bg_signal[iCut], nEvents, sSignalType.Data());
        //scaleHist(fh_true_bg_signal[iCut], nEvents);

        if(fh_true_bg[iCut] && fh_true_bg_signal[iCut]) {
            cout<<"adding bg hists"<<endl;
            fh_true_bg[iCut]->Add(fh_true_bg_signal[iCut]);
        }

    }

    TH1F *fh_pi0, *fh_gamma, *fh_eta;
    histName = "fh_ee_minv_Global Track_pi0;1";
    fh_pi0 = (TH1F*)gDirectory->Get(histName.c_str());
    if(!fh_pi0) {
        cout<<"Histogram "<<dirName.c_str()<<"/"<<histName.c_str()<<" not found."<<endl;
    } else {
        scaleHist(fh_pi0, nEvents);
    }

    histName = "fh_ee_minv_Global Track_gamma;1";
    fh_gamma = (TH1F*)gDirectory->Get(histName.c_str());
    if(!fh_gamma) {
        cout<<"Histogram "<<dirName.c_str()<<"/"<<histName.c_str()<<" not found."<<endl;
    } else {
        scaleHist(fh_gamma, nEvents);
    }
    histName = "fh_ee_minv_Global Track_eta;1";
    fh_eta = (TH1F*)gDirectory->Get(histName.c_str());
    if(!fh_eta) {
        cout<<"Histogram "<<dirName.c_str()<<"/"<<histName.c_str()<<" not found."<<endl;
    } else {
        scaleHist(fh_eta, nEvents);
    }

    TCanvas *cTrueBg = new TCanvas("cTrueBgCuts", "True Combinatorial Background", 800, 800);
    cTrueBg->cd();
    cTrueBg->SetLogy();

    fh_true_bg[0]->GetYaxis()->SetRangeUser(5.e-5, 60.);
    fh_true_bg[0]->GetYaxis()->SetTitle("#frac{dN}{dM_{ee}} [(GeV/c^{2})^{-1}]");
    fh_true_bg[0]->Draw("hist");
    if(fh_signal[0]) {
        fh_signal[0]->Draw("hist sames");
    }
    for(Int_t iCut = 1; iCut < CbmMvdAnaNames::kNofCutTypes; iCut++) {
        if(bDrawCut[iCut]) {
            if(fh_signal[iCut]) {
                fh_signal[iCut]->Draw("hist sames");
            }
            if(fh_true_bg[iCut]) {
                fh_true_bg[iCut]->Draw("hist sames");
            }
        }
    }
    //fh_signal->DrawCopy("hist sames");

    TLegend* legend = new TLegend(0.62,0.65,0.96,0.95, "","brNDC");
    TLegendEntry* legentry = NULL;
    //legentry = legend->AddEntry(fh_signal, "#rho^{0} #rightarrow e^{+} e^{-}","lpf");
    for(Int_t iCut = 0; iCut < CbmMvdAnaNames::kNofCutTypes; iCut++) {
        if(bDrawCut[iCut] && fh_true_bg[iCut]) {
            legentry = legend->AddEntry(fh_true_bg[iCut],
                                        CutTypesNames[iCut], "lpf");
        }
    }
    legend->SetFillColor(kWhite);
    legend->Draw();

    TCanvas *cTrueBgNoCut = new TCanvas("cTrueBgNoCut", "cTrueBgNoCut", 800, 800);
    cTrueBgNoCut->cd();
    cTrueBgNoCut->SetLogy();
    THStack *hstackTrueBgNoCut = new THStack("hstackTrueBgNoCut", "hstackTrueBgNoCut");
    if(fh_true_bg[0]) {
        hstackTrueBgNoCut->Add(fh_true_bg[0]);
    }
    if(fh_signal[0]) {
        fh_signal[0]->SetLineColor(kMagenta);
        hstackTrueBgNoCut->Add(fh_signal[0]);
    }
    if(fh_pi0) {
        fh_pi0->SetLineColor(kRed);
        hstackTrueBgNoCut->Add(fh_pi0);
    }
    if(fh_gamma) {
        fh_gamma->SetLineColor(kBlue);
        hstackTrueBgNoCut->Add(fh_gamma);
    }
    if(fh_eta) {
        fh_eta->SetLineColor(kGreen+2);
        hstackTrueBgNoCut->Add(fh_eta);
    }
    hstackTrueBgNoCut->Draw("hist nostack");
    cTrueBgNoCut->Update();
    hstackTrueBgNoCut->GetYaxis()->SetRangeUser(1.e-5, 60.);

    hstackTrueBgNoCut->GetXaxis()->SetTitle(fh_true_bg[0]->GetXaxis()->GetTitle());
    hstackTrueBgNoCut->GetXaxis()->SetTitleOffset(1.2);
    hstackTrueBgNoCut->GetXaxis()->SetTitleFont(fh_true_bg[0]->GetXaxis()->GetTitleFont());
    hstackTrueBgNoCut->GetXaxis()->SetTitleSize(fh_true_bg[0]->GetXaxis()->GetTitleSize());
    hstackTrueBgNoCut->GetXaxis()->SetLabelFont(fh_true_bg[0]->GetXaxis()->GetLabelFont());
    hstackTrueBgNoCut->GetXaxis()->SetLabelSize(fh_true_bg[0]->GetXaxis()->GetLabelSize());

    hstackTrueBgNoCut->GetYaxis()->SetTitle("#frac{dN}{dM_{ee}} [(GeV/c^{2})^{-1}]");
    hstackTrueBgNoCut->GetYaxis()->SetTitleOffset(1.1);
    hstackTrueBgNoCut->GetYaxis()->SetTitleFont(fh_true_bg[0]->GetYaxis()->GetTitleFont());
    hstackTrueBgNoCut->GetYaxis()->SetTitleSize(fh_true_bg[0]->GetYaxis()->GetTitleSize());
    hstackTrueBgNoCut->GetYaxis()->SetLabelFont(fh_true_bg[0]->GetYaxis()->GetLabelFont());
    hstackTrueBgNoCut->GetYaxis()->SetLabelSize(fh_true_bg[0]->GetYaxis()->GetLabelSize());

    TLegend* legendSource = new TLegend(0.62,0.65,0.96,0.95, "","brNDC");
    TLegendEntry *legentrySource = NULL;
    legentrySource = legendSource->AddEntry(fh_true_bg[0], "True Background", "lpf");
    legentrySource = legendSource->AddEntry(fh_signal[0],  "Signal", "lpf");
    legentrySource = legendSource->AddEntry(fh_pi0,        "#pi^{0}-Dalitz", "lpf");
    legentrySource = legendSource->AddEntry(fh_gamma,      "#gamma-Conversion", "lpf");
    legentrySource = legendSource->AddEntry(fh_eta,        "#eta-Dalitz", "lpf");
    legendSource->Draw();
}

void drawTrueBgTracks(TFile *file, Bool_t mc=kFALSE) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("InvariantMass")) {
        if(!file->cd("Invariantmass")) { // stupid typo in analysis task
            return;
        }
    }

    TString smc("");
    if(mc) smc = "_mc";

    Color_t color_signal = kMagenta;
    Color_t color_true_bg[CbmMvdAnaNames::kNofBgTypes];
    color_true_bg[CbmMvdAnaNames::kBgAll]    = kBlack;
    color_true_bg[CbmMvdAnaNames::kBgNoTs]   = kOrange;;
    color_true_bg[CbmMvdAnaNames::kBgNoTf]   = kGreen;
    color_true_bg[CbmMvdAnaNames::kBgNoTsTf] = kCyan;
    color_true_bg[CbmMvdAnaNames::kBgMinv]   = kBlue;
    color_true_bg[CbmMvdAnaNames::kBgPt]     = kViolet;

    TH1F* fh_signal = (TH1F*)gDirectory->Get(Form("fh_ee%s_minv_Global Track_signal;1", smc.Data()));
    if(!fh_signal) {
	cout<<"Histogram fh_ee_minv_Global Track_signal;1 not found"<<endl;
        return;
    }
    fh_signal->SetLineColor(color_signal);
    scaleHistSignal(fh_signal, nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));

    TH1F* fh_true_bg[CbmMvdAnaNames::kNofBgTypes];
    TH1F* fh_true_bg_signal[CbmMvdAnaNames::kNofBgTypes];

    for(Int_t iCut = 0; iCut < CbmMvdAnaNames::kNofBgTypes; iCut++) {
        fh_true_bg[iCut] = (TH1F*)gDirectory->Get(Form("fh_ee%s_minv_true_bg_%s;1",
                                                       smc.Data(),
                                                       CbmMvdAnaNames::kBgTypesNames[iCut]));
        if(!fh_true_bg[iCut]) {
	    cout<<"Histogram true bg for background type "<<iCut<<" not found."<<endl;
            return;
	}
	fh_true_bg[iCut]->SetLineColor(color_true_bg[iCut]);
        scaleHist(fh_true_bg[iCut], nEvents);

        fh_true_bg_signal[iCut] = (TH1F*)gDirectory->Get(Form("fh_ee%s_minv_true_bg_signal_%s;1",
                                                       smc.Data(),
                                                       CbmMvdAnaNames::kBgTypesNames[iCut]));
        if(!fh_true_bg_signal[iCut]) {
	    cout<<"Histogram true bg_signal for background type "<<iCut<<" not found."<<endl;
            return;
	}
	fh_true_bg_signal[iCut]->SetLineColor(color_true_bg[iCut]);
        scaleHistSignal(fh_true_bg_signal[iCut], nEvents, sSignalType);

        if(fh_true_bg[iCut] && fh_true_bg_signal[iCut]) {
            fh_true_bg[iCut]->Add(fh_true_bg_signal[iCut]);
        }
    }

    TString sCanvasTitle("True Combinatorial Background");
    if(mc) sCanvasTitle.Append(" MC");
    else sCanvasTitle.Append(" Reconstructed");
    TCanvas *cTrueBg = new TCanvas(Form("cTrueBgTracks%s", smc.Data()),
                                   sCanvasTitle.Data(),
                                        800, 800);
    cTrueBg->cd();
    cTrueBg->SetLogy();

    fh_true_bg[0]->GetYaxis()->SetRangeUser(1.e-5, 40.);
    fh_true_bg[0]->Draw("hist");
    for(Int_t iCut = 1; iCut < CbmMvdAnaNames::kNofBgTypes; iCut++) {
	fh_true_bg[iCut]->Draw("hist sames");
    }
    fh_signal->DrawCopy("hist sames");

    TLegend* legend = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
    TLegendEntry* legentry = NULL;
    legentry = legend->AddEntry(fh_signal, "#rho^{0} #rightarrow e^{+} e^{-}","lpf");
    for(Int_t iCut = 0; iCut < CbmMvdAnaNames::kNofBgTypes; iCut++) {
	legentry = legend->AddEntry(fh_true_bg[iCut],
				    CbmMvdAnaNames::kBgTypesNames[iCut], "lpf");
    }
    legend->SetFillColor(kWhite);
    legend->Draw();

}

Bool_t getCutEfficiencies(TFile *file, TH1D *&h_cut_effs_signal,
                          TH1D *&h_cut_effs_bg, TH1D *&h_cut_sb) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("InvariantMass")) {
        cout<<"Directory InvariantMass not found"<<endl;
        return kFALSE;
    }



    TH1D* fh_signal[nofCuts];

    for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
        fh_signal[iCut] = (TH1D*)gDirectory->Get(Form("fh_ee_minv_signal_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[cutTypes[iCut]]));
        if(!fh_signal[iCut]) {
	    cout<<"Histogram signal for cut "<<iCut<<" not found."<<endl;
            continue;
	}
        scaleHistSignal(fh_signal[iCut], nEvents, Form("%s.%s.%s", sLvm.Data(), sDecay.Data(), sEnergy.Data()));
    }

    TH1D* fh_true_bg[nofCuts];
    TH1D* fh_true_bg_signal[nofCuts];

    for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
        fh_true_bg[iCut] = (TH1D*)gDirectory->Get(Form("fh_ee_minv_true_bg_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[cutTypes[iCut]]));
        if(!fh_true_bg[iCut]) {
	    cout<<"Histogram true bg for cut "<<iCut<<" not found."<<endl;
            continue;
	}
        scaleHist(fh_true_bg[iCut], nEvents);

        fh_true_bg_signal[iCut] = (TH1D*)gDirectory->Get(Form("fh_ee_minv_true_bg_signal_cut_%s;1",
                                                       CbmMvdAnaNames::kCutTypesNames[cutTypes[iCut]]));
        if(!fh_true_bg_signal[iCut]) {
	    cout<<"Histogram true bg_signal for cut "<<iCut<<" not found."<<endl;
            continue;
	}
        scaleHistSignal(fh_true_bg_signal[iCut], nEvents, sSignalType);

        if(fh_true_bg[iCut] && fh_true_bg_signal[iCut]) {
            fh_true_bg[iCut]->Add(fh_true_bg_signal[iCut]);
        }

    }

    h_cut_effs_signal = new TH1D("h_cut_effs_signal", "Cut Efficiencies",
                                 nofCuts, 0, nofCuts);

    h_cut_effs_bg     = new TH1D("h_cut_effs_bg", "Cut Efficiencies",
                                 nofCuts, 0, nofCuts);

    h_cut_sb          = new TH1D("h_cut_sb", "Signal to Background Ratio",
                                 nofCuts, 0, nofCuts);

    Float_t minMom = 0.55;
    Float_t maxMom = 1.2;

    Float_t denomSignal = 1.;
    if(fh_signal[0]) {
        Int_t minBin = fh_signal[0]->GetXaxis()->FindBin(minMom);
        Int_t maxBin = fh_signal[0]->GetXaxis()->FindBin(maxMom);
        denomSignal = fh_signal[0]->Integral(minBin, maxBin);
    }

    Float_t denomBg     = 1.;
    if(fh_true_bg[0]) {
        Int_t minBin = fh_true_bg[0]->GetXaxis()->FindBin(minMom);
        Int_t maxBin = fh_true_bg[0]->GetXaxis()->FindBin(maxMom);
        denomBg = fh_true_bg[0]->Integral(minBin, maxBin);
    }

    for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
        Double_t nSignal = 0.F;
        Double_t nBg = 1.F;
        if(fh_signal[iCut]) {
            Int_t minBin = fh_signal[iCut]->GetXaxis()->FindBin(minMom);
            Int_t maxBin = fh_signal[iCut]->GetXaxis()->FindBin(maxMom);
            nSignal      = fh_signal[iCut]->Integral(minBin, maxBin);
            cout<<"signal: minBin "<<minBin<<", maxBin "<<maxBin<<endl;
            cout<<"signal: n "<<nSignal<<", denom "<<denomSignal<<endl;
            Double_t effSignal = nSignal  / denomSignal;
            h_cut_effs_signal->SetBinContent(iCut+1, effSignal);
        } else {
            h_cut_effs_signal->SetBinContent(iCut+1, h_cut_effs_signal->GetBinContent(iCut));
        }
        if(fh_true_bg[iCut]) {
            Int_t minBin = fh_true_bg[iCut]->GetXaxis()->FindBin(minMom);
            Int_t maxBin = fh_true_bg[iCut]->GetXaxis()->FindBin(maxMom);
            nBg          = fh_true_bg[iCut]->Integral(minBin, maxBin);
            Double_t effBg     = nBg / denomBg;
            cout<<"bg: minBin "<<minBin<<", maxBin "<<maxBin<<endl;
            cout<<"bg: n "<<nBg<<", denom "<<denomBg<<endl;
            h_cut_effs_bg    ->SetBinContent(iCut+1, effBg);
        } else {
            h_cut_effs_bg->SetBinContent(iCut+1, h_cut_effs_bg->GetBinContent(iCut));
        }
	if(fh_signal[iCut] && fh_true_bg[iCut]) {
	    if(nBg > 0.) {
		Float_t sb = nSignal / nBg;
		h_cut_sb->SetBinContent(iCut+1, sb);
	    }
        } else {
            h_cut_sb->SetBinContent(iCut+1, h_cut_sb->GetBinContent(iCut));
        }
    }

    for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
        h_cut_effs_signal->GetXaxis()->SetBinLabel(iCut+1, cutTypesNames[iCut]);
    }
    h_cut_effs_signal->SetLineColor(kRed);
    h_cut_effs_bg    ->SetLineColor(kBlack);
    h_cut_effs_signal->GetYaxis()->SetRangeUser(1.e-3, 1.3);
    h_cut_effs_signal->GetYaxis()->SetTitle("Cut Efficiency");

    for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
        h_cut_sb->GetXaxis()->SetBinLabel(iCut+1, cutTypesNames[iCut]);
        h_cut_sb->GetYaxis()->SetTitle("S/B");
    }

    file->cd("/");

    return kTRUE;
}

void drawCutEfficiencies(TFile *file) {

    TH1D *h_cut_effs_signal = NULL;
    TH1D *h_cut_effs_bg     = NULL;
    TH1D *h_cut_sb          = NULL;

    if(!getCutEfficiencies(file, h_cut_effs_signal, h_cut_effs_bg,
                          h_cut_sb)) {
        return;
    }

    TCanvas *cCutEff = new TCanvas("cCutEff", "Cut Efficiencies", 800, 800);
    cCutEff->cd();
    cCutEff->SetLogy();

    if(h_cut_effs_signal)
        h_cut_effs_signal->Draw("");
    else
        cout<<"Histogram with cut efficiencies for signal does not exist."<<endl;
    if(h_cut_effs_bg)
        h_cut_effs_bg    ->Draw("same");

    TLegend* legend = new TLegend(0.65,0.83,0.95,0.99, "","brNDC");
    TLegendEntry* legentry = NULL;
    legentry = legend->AddEntry(h_cut_effs_signal, "Signal","lpf");
    legentry = legend->AddEntry(h_cut_effs_bg, "Background","lpf");
    legend->Draw("same");

    TCanvas *cCutSb = new TCanvas("cCutSb", "Signal to Background", 800, 800);
    cCutSb->cd();

    if(h_cut_sb) {
        h_cut_sb->SetMarkerStyle(kFullSquare);
        h_cut_sb->Draw("lp");
    }

}

void drawMinStations(TFile *file) {
    Int_t nEvents = getNevents(file);

    TString dirName = "HitsInStsTrack";
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
	return;
    }

    TString histName("");

    histName = "fh_e_mvdsts_hits_min_station_all_signal";
    TH1F *fh_min_station_signal = (TH1F*)gDirectory->Get(histName.Data());
    if(!fh_min_station_signal) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = "fh_e_mvdsts_hits_min_station_sec_gamma";
    TH1F *fh_min_station_sec_gamma = (TH1F*)gDirectory->Get(histName.Data());
    if(!fh_min_station_sec_gamma) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = "fh_e_mvdsts_hits_min_station_all_gamma";
    TH1F *fh_min_station_all_gamma = (TH1F*)gDirectory->Get(histName.Data());
    if(!fh_min_station_all_gamma) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    histName = "fh_e_mvdsts_hits_min_station_prim_gamma";
    TH1F *fh_min_station_prim_gamma = (TH1F*)gDirectory->Get(histName.Data());
    if(!fh_min_station_prim_gamma) {
	cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
        return;
    }

    TString histFileReco = file->GetName();
    Double_t br   = getDecayBR(histFileReco);
    Double_t mult = getDecayMultiplicity(histFileReco);


    TGaxis::SetMaxDigits(3);

    TCanvas *cMinSignal = new TCanvas("cMinSignal", "First Stations Signal", 800, 800);
    cMinSignal->cd();

    fh_min_station_signal->Scale(br * mult / (Double_t)nEvents);
    fh_min_station_signal->GetYaxis()->SetTitle("Yield per Event");
    fh_min_station_signal->GetYaxis()->SetTitleOffset(1.1);
    fh_min_station_signal->Draw("");

    TCanvas *cMinGamma = new TCanvas("cMinGamma", "First Stations Gammas",
                                    800, 800);
    cMinGamma->cd();

    fh_min_station_all_gamma->Scale(1. / (Double_t)nEvents);
    fh_min_station_all_gamma->GetYaxis()->SetTitle("Yield per Event");
    fh_min_station_all_gamma->GetYaxis()->SetTitleOffset(1.1);
    fh_min_station_all_gamma->SetLineColor(kBlack);
    fh_min_station_all_gamma->SetLineWidth(2);

    fh_min_station_prim_gamma->Scale(1. / (Double_t)nEvents);
    fh_min_station_prim_gamma->GetYaxis()->SetTitle("Yield per Event");
    fh_min_station_prim_gamma->GetYaxis()->SetTitleOffset(1.1);
    fh_min_station_prim_gamma->SetLineColor(kRed);
    fh_min_station_prim_gamma->SetLineWidth(2);

    fh_min_station_sec_gamma->Scale(1. / (Double_t)nEvents);
    fh_min_station_sec_gamma->GetYaxis()->SetTitle("Yield per Event");
    fh_min_station_sec_gamma->GetYaxis()->SetTitleOffset(1.1);
    fh_min_station_sec_gamma->SetLineColor(kBlue);
    fh_min_station_sec_gamma->SetLineWidth(2);

    THStack *hStackGamma = new THStack("hStackGamma", "hStackGamma");
    hStackGamma->Add(fh_min_station_all_gamma);
    hStackGamma->Add(fh_min_station_prim_gamma);
    hStackGamma->Add(fh_min_station_sec_gamma);

    hStackGamma->Draw("nostack hist");
    setStackAxis(hStackGamma);
    hStackGamma->Draw("nostack hist");

    TLegend *legendGamma = new TLegend(0.50,0.72,0.95,0.95, "","brNDC");
    legendGamma->AddEntry(fh_min_station_all_gamma,
                          Form("#gamma-conversions: %.2f",
                               fh_min_station_all_gamma->GetEntries() / nEvents),"lpf");
    legendGamma->AddEntry(fh_min_station_prim_gamma,
                          Form("Primary #gamma-conversions: %.2f",
                               fh_min_station_prim_gamma->GetEntries() / nEvents),"lpf");
    legendGamma->AddEntry(fh_min_station_sec_gamma,
                          Form("Secondary #gamma-conversions: %.2f",
                               fh_min_station_sec_gamma->GetEntries() / nEvents),"lpf");
    legendGamma->Draw("same");
}

void drawPrimaryVertexCut(TFile *file) {
    Int_t nEvents = getNevents(file);
    if(!file->cd("PrimaryTrackCut")) {
        cout<<"Directory PrimaryTrackCut not found"<<endl;
	return;
    }

    TH1F *fh_true_z_cut[CbmMvdAnaNames::kNofSourceTypes][nofCuts];
    TString histName = "";

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
            histName = Form("fh_e_vertex_true_z_%s_cut_%s",
                            CbmMvdAnaNames::kSourceTypesNames[iSource],
                            CbmMvdAnaNames::kCutTypesNames[cutTypes[iCut]]);
            fh_true_z_cut[iSource][iCut] = (TH1F*)gDirectory->Get(histName.Data());
            if(!fh_true_z_cut[iSource][iCut]) {
                cout<<"Histogram "<<histName.Data()<<" not found."<<endl;
            }
        }
    }

    THStack *hStackTrueZCut[CbmMvdAnaNames::kNofSourceTypes];
    TLegend *legendTrueZCut = new TLegend(0.58,0.82,0.98,0.97, "","brNDC");

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        hStackTrueZCut[iSource] = new THStack(Form("hStackTrueZCut%s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                     Form("hStackTrueZCut%s", CbmMvdAnaNames::kSourceTypesNames[iSource]));
        for(Int_t iCut = 0; iCut < nofCuts; iCut++) {
            TH1F *h = fh_true_z_cut[iSource][iCut];
            if(!h) continue;
            if(iSource == CbmMvdAnaNames::kSourceSignal) {
                scaleHistSignal(h, nEvents, sSignalType.Data());
            } else {
                scaleHist(h, nEvents);
            }
            h->SetLineColor(cutTypesColors[iCut]);
            if(cutTypes[iCut] <= CbmMvdAnaNames::kCutMvdPrimary &&
               cutTypes[iCut] != CbmMvdAnaNames::kCutGamma) {
                hStackTrueZCut[iSource]->Add(h);
                if(iSource == 0) {
                    legendTrueZCut->AddEntry(h, cutTypesNames[iCut],"lpf");
                }
            }
        }
    }

    TCanvas *cTrueZCut[CbmMvdAnaNames::kNofSourceTypes];
    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
        cTrueZCut[iSource] = new TCanvas(Form("cTrueZCut%s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                        Form("cTrueZCut%s", CbmMvdAnaNames::kSourceTypesNames[iSource]),
                                        800, 800);
        cTrueZCut[iSource]->cd();
        cTrueZCut[iSource]->SetLogy();
        setStackAxis(hStackTrueZCut[iSource]);
        hStackTrueZCut[iSource]->GetYaxis()->SetTitle("#frac{dN}{dz} [cm]^{-1}");
        //hStackTrueZCut[iSource]->GetYaxis()->SetTitleOffset(1.3);

        legendTrueZCut->Draw("same");
    }
/*
    TH1F* fh_true_z_all_gamma = (TH1F*)gDirectory->Get("fh_e_vertex_true_z_all_gamma;1");
    if(!fh_true_z_all_gamma) {
	cout<<"Histogram fh_e_vertex_true_z_all_gamma;1 not found"<<endl;
        return;
    }
    fh_true_z_all_gamma->Scale(1.F/(Float_t)nEvents);

    TH1F* fh_true_z_secondary_gamma = (TH1F*)gDirectory->Get("fh_e_vertex_true_z_cut_chi2prim_gamma;1");
    if(!fh_true_z_secondary_gamma) {
	cout<<"Histogram fh_e_vertex_true_z_cut_chi2prim_gamma;1 not found"<<endl;
        return;
    }
    fh_true_z_secondary_gamma->Scale(1.F/(Float_t)nEvents);

    TH1F *fh_true_z_chi2prim_gamma = new TH1F("fh_e_vertex_true_z_prim_gamma", "fh_e_vertex_true_z_prim_gamma",
					      fh_true_z_secondary_gamma->GetNbinsX(),
					      fh_true_z_secondary_gamma->GetXaxis()->GetXmin(),
					      fh_true_z_secondary_gamma->GetXaxis()->GetXmax());
    // Correct stupid inverse logic for filling histograms in analysis task.
    for(Int_t i = 1; i <= fh_true_z_secondary_gamma->GetNbinsX(); i++) {
	Double_t primaries = fh_true_z_all_gamma->GetBinContent(i) -
	    fh_true_z_secondary_gamma->GetBinContent(i);
	fh_true_z_chi2prim_gamma->SetBinContent(i, primaries);
    }

    TCanvas *cPrimCutZ = new TCanvas("cPrimCutZ", "Primary Track Cut", 800, 800);
    cPrimCutZ->cd();
    cPrimCutZ->SetLogy();

    fh_true_z_all_gamma->SetLineColor(kRed);
    fh_true_z_all_gamma->GetYaxis()->SetTitle("Yield per event");
    fh_true_z_all_gamma->GetYaxis()->SetTitleOffset(1.2);
    fh_true_z_all_gamma->SetTitle("MC Vertex of reconstructed e^{+}/e^{-} #gamma-conversions");
    fh_true_z_all_gamma->Draw("");
    //fh_true_z_secondary_gamma->Draw("sames");
    fh_true_z_chi2prim_gamma->SetLineColor(kBlue);
    fh_true_z_chi2prim_gamma->Draw("sames");

    TLegend* legendZ = new TLegend(0.58,0.77,0.98,0.92, "","brNDC");
    legendZ->AddEntry(fh_true_z_all_gamma, "before primary track cut","lpf");
    legendZ->AddEntry(fh_true_z_chi2prim_gamma, "after primary track cut","lpf");
    legendZ->Draw("same");


    TH1F* fh_vertex_fit_chi2_all_gamma = (TH1F*)gDirectory->Get("fh_e_vertex_fit_chi2_all_gamma;1");
    if(!fh_vertex_fit_chi2_all_gamma) {
	cout<<"Histogram fh_e_vertex_fit_chi2_all_gamma;1 not found"<<endl;
        return;
    }
    fh_vertex_fit_chi2_all_gamma->Scale(1.F/(Float_t)nEvents);

    TH1F* fh_vertex_fit_chi2_trueprim_gamma = (TH1F*)gDirectory->Get("fh_e_vertex_fit_chi2_trueprim_gamma;1");
    if(!fh_vertex_fit_chi2_trueprim_gamma) {
	cout<<"Histogram fh_vertex_fit_chi2_trueprim_gamma;1 not found"<<endl;
        return;
    }
    fh_vertex_fit_chi2_trueprim_gamma->Scale(1.F/(Float_t)nEvents);

    TCanvas *cPrimCutChi2 = new TCanvas("cPrimCutChi2", "Primary Track Cut", 800, 800);
    cPrimCutChi2->cd();

    fh_vertex_fit_chi2_all_gamma->SetLineColor(kRed);
    fh_vertex_fit_chi2_all_gamma->GetXaxis()->SetRangeUser(0., 15.);
    fh_vertex_fit_chi2_all_gamma->GetXaxis()->SetTitle("#chi^{2}");
    fh_vertex_fit_chi2_all_gamma->GetYaxis()->SetTitle("Yield per event");
    fh_vertex_fit_chi2_all_gamma->GetYaxis()->SetTitleOffset(1.4);
    fh_vertex_fit_chi2_all_gamma->GetYaxis()->SetRangeUser(0., 0.55);
    fh_vertex_fit_chi2_all_gamma->SetTitle("Vertex fit of #gamma-Conversions");
    fh_vertex_fit_chi2_all_gamma->Draw("");
    //fh_true_z_secondary_gamma->Draw("sames");
    fh_vertex_fit_chi2_trueprim_gamma->SetLineColor(kBlue);
    fh_vertex_fit_chi2_trueprim_gamma->Draw("sames");

    TLegend* legendChi2 = new TLegend(0.62,0.78,0.97,0.98, "","brNDC");
    legendChi2->AddEntry(fh_vertex_fit_chi2_all_gamma, "all tracks","lpf");
    legendChi2->AddEntry(fh_vertex_fit_chi2_trueprim_gamma, "true primary tracks","lpf");
    legendChi2->Draw("same");
    */
}

void drawFindGammaComb(TFile *file) {

    //Int_t nEvents = getNevents(file);
    Int_t nEvents = 800000;

    TString dirName("Gamma");
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        //return;
    }

    TString histFileName("");

    histFileName = "h_gamma_comb_chi2_true;1";
    TH1F* h_gamma_comb_chi2_true = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_comb_chi2_true) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_gamma_comb_chi2_false;1";
    TH1F* h_gamma_comb_chi2_false = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_comb_chi2_false) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_gamma_comb_nhits_true;1";
    TH1F* h_gamma_comb_nhits_true = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_comb_nhits_true) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_gamma_comb_nhits_false;1";
    TH1F* h_gamma_comb_nhits_false = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_comb_nhits_false) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_gamma_comb_corr;1";
    TH1F* h_gamma_comb_corr = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_comb_corr) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    TCanvas *cCombChi2 = new TCanvas("cCombChi2", "Shi2", 800, 800);
    cCombChi2->cd();
    h_gamma_comb_chi2_true->SetLineColor(kBlue);
    h_gamma_comb_chi2_true->Scale(1.F / (Float_t)nEvents);
    h_gamma_comb_chi2_true->GetXaxis()->SetRangeUser(-1., 15.);
    h_gamma_comb_chi2_true->GetXaxis()->SetTitle("#chi^{2}");
    h_gamma_comb_chi2_true->GetYaxis()->SetTitle("Yield per event");
    h_gamma_comb_chi2_true->GetYaxis()->SetTitleOffset(1.4);
    h_gamma_comb_chi2_true->GetYaxis()->SetTitleColor(h_gamma_comb_chi2_true->GetLineColor());
    h_gamma_comb_chi2_true->GetYaxis()->SetLabelColor(h_gamma_comb_chi2_true->GetLineColor());
    h_gamma_comb_chi2_true->Draw("");

    cCombChi2->Update();

    h_gamma_comb_chi2_false->Scale(1.F / (Float_t)nEvents);
    Float_t rightmax = 1.1*h_gamma_comb_chi2_false->GetMaximum();
    Float_t scale    = gPad->GetUymax()/rightmax;
    h_gamma_comb_chi2_false->SetLineColor(kRed);
    h_gamma_comb_chi2_false->Scale(scale);
    h_gamma_comb_chi2_false->Draw("same");
    //draw an axis on the right side
    TGaxis*axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
        		     gPad->GetUxmax(),gPad->GetUymax(),
        		     0,rightmax,510,"+L");
    axis->SetLineColor(kRed);
    axis->SetLabelColor(kRed);
    axis->SetTitle("Yield per event");
    axis->SetTitleColor(h_gamma_comb_chi2_false->GetLineColor());
    axis->SetTitleFont(h_gamma_comb_chi2_true->GetYaxis()->GetTitleFont());
    axis->SetTitleSize(h_gamma_comb_chi2_true->GetYaxis()->GetTitleSize());
    axis->SetLabelFont(h_gamma_comb_chi2_true->GetYaxis()->GetLabelFont());
    axis->SetLabelSize(h_gamma_comb_chi2_true->GetYaxis()->GetLabelSize());
    axis->Draw();

    TLegend* legendChi2 = new TLegend(0.52,0.78,0.87,0.93, "","brNDC");
    legendChi2->AddEntry(h_gamma_comb_chi2_true, "True Hypotheses","lpf");
    legendChi2->AddEntry(h_gamma_comb_chi2_false, "False Hypotheses","lpf");
    legendChi2->Draw("same");


    TCanvas *cCombHits = new TCanvas("cCombHits", "Hits", 800, 800);
    cCombHits->cd();
    h_gamma_comb_nhits_true->Scale(1.F / (Float_t)h_gamma_comb_nhits_true->GetEntries());
    h_gamma_comb_nhits_true->SetLineColor(kBlue);
    h_gamma_comb_nhits_true->GetXaxis()->SetTitle("Number of Hits");
    h_gamma_comb_nhits_true->GetYaxis()->SetTitle("Fraction");
    h_gamma_comb_nhits_true->GetYaxis()->SetTitleOffset(1.4);
    h_gamma_comb_nhits_true->GetYaxis()->SetRangeUser(0., 0.5);
    h_gamma_comb_nhits_true->Draw("");

    h_gamma_comb_nhits_false->Scale(1.F / (Float_t)h_gamma_comb_nhits_false->GetEntries());

    h_gamma_comb_nhits_false->SetLineColor(kRed);
    h_gamma_comb_nhits_false->Draw("same");

    TLegend* legendHits = new TLegend(0.62,0.78,0.97,0.93, "","brNDC");
    legendHits->AddEntry(h_gamma_comb_nhits_true, "True Hypotheses","lpf");
    legendHits->AddEntry(h_gamma_comb_nhits_false, "False Hypotheses","lpf");
    legendHits->Draw("same");


    TCanvas *cCombCorr = new TCanvas("cCombCorr", "Hypothesis Selection", 800, 800);
    cCombCorr->cd();
    h_gamma_comb_corr->Scale(1.F / (Float_t)nEvents);
    h_gamma_comb_corr->GetYaxis()->SetTitle("Yield per event");
    h_gamma_comb_corr->GetYaxis()->SetTitleOffset(1.1);
    h_gamma_comb_corr->Draw("");
}

void drawFindGammaInput(TFile *file) {

    //Int_t nEvents = getNevents(file);
    Int_t nEvents = 800000;

    TString dirName("Gamma");
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        //return;
    }

    TString histFileName("");

    histFileName = "h_mom_gamma_cand;1";
    TH1F* h_mom_in = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_mom_in) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_angle_gamma_cand;1";
    TH1F* h_angle_in = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_mom_in) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_minv_gamma_cand;1";
    TH1F* h_minv_in = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_mom_in) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dx;1";
    TH1F* h_seed_dx = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dx) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dy;1";
    TH1F* h_seed_dy = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dy) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dz;1";
    TH1F* h_seed_dz = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dz) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dtx;1";
    TH1F* h_seed_dtx = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dtx) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dty;1";
    TH1F* h_seed_dty = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dty) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_seed_error_dp;1";
    TH1F* h_seed_dp = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_seed_dp) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_used_hits;1";
    TH1F* h_used_hits = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_used_hits) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_gamma_vertices;1";
    TH1F* h_gamma_vertices = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_gamma_vertices) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    TCanvas *cMomIn = new TCanvas("cMomIn", "Input momentum", 800, 800);
    cMomIn->cd();
    h_mom_in->Scale(1.F / (Float_t)nEvents);
    h_mom_in->GetYaxis()->SetTitle("Yield per event");
    h_mom_in->GetYaxis()->SetTitleOffset(1.4);
    h_mom_in->Draw();

    TCanvas *cAngleIn = new TCanvas("cAngleIn", "Input opening angle", 800, 800);
    cAngleIn->cd();
    h_angle_in->Scale(1.F / (Float_t)nEvents);
    h_angle_in->GetYaxis()->SetTitle("Yield per event");
    h_angle_in->GetYaxis()->SetTitleOffset(1.4);
    h_angle_in->Draw();

    TCanvas *cMinvIn = new TCanvas("cMinvIn", "Input invariant mass", 800, 800);
    cMinvIn->cd();
    h_minv_in->Scale(1.F / (Float_t)nEvents);
    h_minv_in->GetYaxis()->SetTitle("Yield per event");
    h_minv_in->GetYaxis()->SetTitleOffset(1.4);
    h_minv_in->Draw();

    TCanvas *cSeedDx = new TCanvas("cSeedDx", "Seed Error dx", 800, 800);
    cSeedDx->cd();
    h_seed_dx->Scale(1.F / (Float_t)nEvents);
    h_seed_dx->GetYaxis()->SetTitle("Yield per event");
    h_seed_dx->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dx->Draw();

    TCanvas *cSeedDy = new TCanvas("cSeedDy", "Seed Error dy", 800, 800);
    cSeedDy->cd();
    h_seed_dy->Scale(1.F / (Float_t)nEvents);
    h_seed_dy->GetYaxis()->SetTitle("Yield per event");
    h_seed_dy->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dy->Draw();

    TCanvas *cSeedDz = new TCanvas("cSeedDz", "Seed Error dz", 800, 800);
    cSeedDz->cd();
    h_seed_dz->Scale(1.F / (Float_t)nEvents);
    h_seed_dz->GetYaxis()->SetTitle("Yield per event");
    h_seed_dz->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dz->Draw();

    TCanvas *cSeedDtx = new TCanvas("cSeedDtx", "Seed Error dtx", 800, 800);
    cSeedDtx->cd();
    h_seed_dtx->Scale(1.F / (Float_t)nEvents);
    h_seed_dtx->GetYaxis()->SetTitle("Yield per event");
    h_seed_dtx->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dtx->Draw();

    TCanvas *cSeedDty = new TCanvas("cSeedDty", "Seed Error dty", 800, 800);
    cSeedDty->cd();
    h_seed_dty->Scale(1.F / (Float_t)nEvents);
    h_seed_dty->GetYaxis()->SetTitle("Yield per event");
    h_seed_dty->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dty->Draw();

    TCanvas *cSeedDp = new TCanvas("cSeedDp", "Seed Error dp", 800, 800);
    cSeedDp->cd();
    h_seed_dp->Scale(1.F / (Float_t)nEvents);
    h_seed_dp->GetYaxis()->SetTitle("Yield per event");
    h_seed_dp->GetYaxis()->SetTitleOffset(1.3);
    h_seed_dp->Draw();

    TCanvas *cUsedHits = new TCanvas("cUsedHits", "Used Hits", 800, 800);
    cUsedHits->cd();
    Int_t usedMvd = h_used_hits->GetBinContent(1);
    Int_t unusedMvd = h_used_hits->GetBinContent(2);
    Int_t sumMvd = usedMvd + unusedMvd;
    Int_t usedSts = h_used_hits->GetBinContent(3);
    Int_t unusedSts = h_used_hits->GetBinContent(4);
    Int_t sumSts = usedSts + unusedSts;
    h_used_hits->SetBinContent(1, (Double_t)usedMvd / (Double_t)sumMvd);
    h_used_hits->SetBinContent(2, (Double_t)unusedMvd / (Double_t)sumMvd);
    h_used_hits->SetBinContent(3, (Double_t)usedSts / (Double_t)sumSts);
    h_used_hits->SetBinContent(4, (Double_t)unusedSts / (Double_t)sumSts);
    h_used_hits->GetYaxis()->SetTitle("Fraction");
    h_used_hits->GetYaxis()->SetRangeUser(0., 1.05);
    h_used_hits->Draw();

    TCanvas *cGammaVertices = new TCanvas("cGammaVertices", "Gamma Vertices",
					  800, 800);
    cGammaVertices->cd();
    h_gamma_vertices->Draw();

}

void drawFindGammaOutput(TFile *file) {

    //Int_t nEvents = getNevents(file);
    Int_t nEvents = 800000;

    TString dirName("Gamma");
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        //return;
    }

    TString histFileName("");

    histFileName = "h_e_mom_eff_corr_hits;1";
    TH1F* h_eff_hits = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_eff_hits) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_e_mom_nof_tracks;1";
    TH1F* h_nof_tracks = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_nof_tracks) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_e_track_types;1";
    TH1F* h_e_track_types = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_e_track_types) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    histFileName = "h_e_mom_eff_find_tracks;1";
    TH1F* h_eff_tracks = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_eff_tracks) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
        return;
    }

    const Int_t nMvd = 5;
    Color_t color_hit_types[nMvd] = { kViolet, kBlue, kGreen+2, kOrange, kRed };
    TH1F *h_hit_types_mvd[nMvd] = { NULL };
    Int_t iMvd = 0;
    Bool_t bFoundHist = kTRUE;
    Double_t ymax = 0.;
    while(bFoundHist && iMvd <= nMvd) {
	// Names start with 1
	histFileName = Form("h_e_hit_types_mvd_%i;1", iMvd+1);
	h_hit_types_mvd[iMvd] = (TH1F*)gDirectory->Get(histFileName.Data());
	if(h_hit_types_mvd[iMvd]) {
            h_hit_types_mvd[iMvd]->Scale(100./h_hit_types_mvd[iMvd]->GetEntries());
	    ymax = std::max(ymax, h_hit_types_mvd[iMvd]->GetMaximum()*1.1);
	    iMvd++;
	} else {
	    bFoundHist = kFALSE;
	}
    }

    if(iMvd > 0) {
	TCanvas *cHitTypeMvd = new TCanvas("cHitTypes", "Hit Topology", 800, 800);
	cHitTypeMvd->cd();
	TLegend* legendHitTypes = new TLegend(0.62,0.72,0.92,0.88, "","brNDC");
	for(Int_t i = 0; i <= iMvd; i++) {
	    if(h_hit_types_mvd[i]) {
		h_hit_types_mvd[i]->SetLineColor(color_hit_types[i]);
		h_hit_types_mvd[i]->SetLineWidth(2);
		if(i == 0) {
		    h_hit_types_mvd[i]->GetYaxis()->SetRangeUser(0., 105.);
		    h_hit_types_mvd[i]->GetYaxis()->SetTitle("Fraction of Hits [%]");
		    h_hit_types_mvd[i]->GetYaxis()->SetTitleOffset(1.2);
		    h_hit_types_mvd[i]->Draw("");
		} else {
		    h_hit_types_mvd[i]->Draw("sames");
		}
		//
		switch(i) {
		case 0:
		    legendHitTypes->AddEntry(h_hit_types_mvd[i], "1st MVD", "lpf");
		    break;
		case 1:
		    legendHitTypes->AddEntry(h_hit_types_mvd[i], "2nd MVD", "lpf");
		    break;
		case 2:
		    legendHitTypes->AddEntry(h_hit_types_mvd[i], "3rd MVD", "lpf");
		    break;
		default:
		    legendHitTypes->AddEntry(h_hit_types_mvd[i], Form("%ith MVD", i+1), "lpf");
		    break;
		}
	    }
	}
	legendHitTypes->Draw("same");
    }

    const Int_t nSts = 8;
    TH1F *h_hit_types_sts[nSts] = { NULL };
    TCanvas *cHitTypesSts = new TCanvas("cHitTypesSts", "Hit Types in Sts",
					1200, 600);
    cHitTypesSts->Divide(nSts/2 ,2);

    for(Int_t iSts = 0; iSts < nSts; iSts++) {
	histFileName = Form("h_e_hit_types_sts_%i;1", iSts+1);
	h_hit_types_sts[iSts] = (TH1F*)gDirectory->Get(histFileName.Data());
        TH1F *h = h_hit_types_sts[iSts];
	if(h) {
	    h->Scale(1.F / (Float_t)h->GetEntries());
	    h->SetTitle(Form("STS %i", iSts+1));
            h->GetYaxis()->SetRangeUser(0., 1.05);
	    h->GetYaxis()->SetTitle("Fraction of Hits");
	    h->GetYaxis()->SetTitleOffset(1.2);
	    cHitTypesSts->cd(iSts+1);
            h->Draw("");
	} else {
	    cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	}

    }

    TCanvas *cEffHits = new TCanvas("cEffHits", "Efficiency MVD Hit Assignment", 800, 800);
    cEffHits->cd();
    h_eff_hits->GetYaxis()->SetTitle("Efficiency");
    h_eff_hits->GetYaxis()->SetTitleOffset(1.1);
    h_eff_hits->GetYaxis()->SetRangeUser(0., 1.);
    h_eff_hits->Draw();

    TCanvas *cNofTracks = new TCanvas("cNofTracks", "Input Tracks", 800, 800);
    cNofTracks->cd();
    h_nof_tracks->Scale(1.F / (Float_t)nEvents);
    h_nof_tracks->GetYaxis()->SetTitle("Yield per event");
    h_nof_tracks->Draw();

    TCanvas *cTrackTypes = new TCanvas("cTrackTypes", "Output Tracks", 800, 800);
    cTrackTypes->cd();
    h_e_track_types->Scale(1.F / (Float_t)h_e_track_types->GetEntries());
    h_e_track_types->GetYaxis()->SetTitle("Fraction of total tracks");
    h_e_track_types->GetYaxis()->SetTitleOffset(1.3);
    h_e_track_types->Draw();

    TCanvas *cEffTracks = new TCanvas("cEffTracks", "Efficiency Track Finding", 800, 800);
    cEffTracks->cd();
    h_eff_tracks->GetYaxis()->SetTitle("Efficiency");
    h_eff_tracks->GetYaxis()->SetTitleOffset(1.1);
    h_eff_tracks->GetYaxis()->SetRangeUser(0., 1.);
    h_eff_tracks->Draw();
}

void drawKfpInput(TFile *file) {
    Int_t nEvents = getNevents(file);
    TString dirName("KFParticle");
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        return;
    }

    TString histFileName("");

    TH1F *h_kfp_mc_m[CbmMvdAnaNames::kNofSourceTypes];

    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
	histFileName = Form("fh_kfp_true_pair_mc_m_%s;1", +
			    CbmMvdAnaNames::kSourceTypesNames[iSource]);
	h_kfp_mc_m[iSource] = (TH1F*)gDirectory->Get(histFileName.Data());
	if(!h_kfp_mc_m[iSource]) {
	    cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	    return;
	}
    }

    histFileName = "fh_kfp_true_pair_mc_z";
    TH1F *h_kfp_mc_z = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_mc_z) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    Float_t nPairs[CbmMvdAnaNames::kNofSourceTypes];
    Float_t sumPairs = 0.F;
    for(Int_t iSource = 0; iSource < CbmMvdAnaNames::kNofSourceTypes; iSource++) {
	nPairs[iSource] = (Float_t)h_kfp_mc_m[iSource]->GetEntries() / (Float_t)nEvents;
        sumPairs += nPairs[iSource];
    }
    nPairs[CbmMvdAnaNames::kSourceSignal] *= getDecayBR(sSignalType) * getDecayMultiplicity(sSignalType);

    TH1F *h_pairs = new TH1F("h_pairs", "", 4, -0.5, 3.5);
    h_pairs->GetXaxis()->SetBinLabel(1, "Signal");
    h_pairs->GetXaxis()->SetBinLabel(2, "#pi^{0}");
    h_pairs->GetXaxis()->SetBinLabel(3, "#gamma");
    h_pairs->GetXaxis()->SetBinLabel(4, "#eta");

    h_pairs->GetYaxis()->SetTitle("Yield per event");

    h_pairs->SetBinContent(1, nPairs[CbmMvdAnaNames::kSourceSignal]);
    h_pairs->SetBinContent(2, nPairs[CbmMvdAnaNames::kSourcePi0]);
    h_pairs->SetBinContent(3, nPairs[CbmMvdAnaNames::kSourceGamma]);
    h_pairs->SetBinContent(4, nPairs[CbmMvdAnaNames::kSourceEta]);


    TCanvas *cPairs = new TCanvas("cPairs", "Dilepton Sources", 800, 800);
    cPairs->cd();
    cPairs->SetLogy();
    h_pairs->Draw();

    TCanvas *cPairsZ = new TCanvas("cPairsZ", "Dilepton Sources", 800, 800);
    cPairsZ->cd();
    cPairsZ->SetLogy();
    h_kfp_mc_z->Scale(1.F / (Float_t)nEvents);
    h_kfp_mc_z->GetYaxis()->SetTitle("Yield per event");
    //scaleHist(h_kfp_mc_z, nEvents);
    //h_kfp_mc_z->GetYaxis()->SetTitle("dN/dz [1/cm]");
    h_kfp_mc_z->Draw("hist");

    histFileName = "fh_kfp_pair_type";
    TH1F *h_kfp_pair_type = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_pair_type) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }
    h_kfp_pair_type->Scale(1.F/(Float_t)nEvents);
    h_kfp_pair_type->GetYaxis()->SetTitle("Yield per event");
    h_kfp_pair_type->GetYaxis()->SetTitleOffset(1.5);
    TCanvas *cPairType = new TCanvas("cPairType", "cPairType", 800, 800);
    cPairType->cd();
    h_kfp_pair_type->Draw("hist");
}

void drawKfpObservables(TFile *file) {
    Int_t nEvents = getNevents(file);
    TString dirName("KFParticle");

    const Int_t nObs = 6;
    string fObs[nObs] = {"angle", "chi2", "m", "merr", "z", "pt" };

    string histPrefix = "fh_kfp_";
    const Int_t nTypes = 3;
    string fTypes[nTypes] = {"fake_pairs", "prim_pairs", "sec_pairs" };
    Color_t fTypeCols[nTypes] = { kRed, kBlue, kViolet };
    string fTypeLegend[nTypes] = { "Fake pairs", "Primary pairs", "Secondary pairs" };
    Bool_t fTypeDraw[nTypes] = { kTRUE, kTRUE, kTRUE };
    string fObsYAxisTitle[nObs] = { "Yield per event", "Yield per event",
    "Yield per event", "Yield per event",
    "Yield per event", "Yield per event" };


/*
    string histPrefix = "fh_kfp_true_pair_";
    const Int_t nTypes = 4;
    string fTypes[nTypes] = {"signal", "gamma", "pi0", "eta" };
    Color_t fTypeCols[nTypes] = { kMagenta, kBlue, kRed, kGreen+2 };
    string fTypeLegend[nTypes] = { "Signal", "#gamma", "#pi^{0}-Dalitz", "#eta" };
    Bool_t fTypeDraw[nTypes] = { kTRUE, kTRUE, kTRUE, kTRUE };
    string fObsYAxisTitle[nObs] = { "dN/d#theta [1/#circ]", "dN/d#chi^{2}",
    "dN/dm_{inv} [GeV/c^{2}]^{-1}", "dN/d#Delta m_{inv} [GeV/c^{2}]^{-1}",
    "dN/dz [1/cm]", "dN/dp_{t} [GeV/c]^{-1}" };
*/

    const Int_t nHypos = 5;
    string fHypoDir[nHypos] = { "IsPrimary", "NoVtxGuess", "Pion", "VtxGuess",
                                "IdealVtxGuess"};

    TH1F *fh_kfp[nObs][nHypos][nTypes];
    // pt same for every hypothesis.
    TH1F *fh_kfp_pt[nTypes];

    TString histFileName("");

    for(Int_t iObs = 0; iObs < nObs; iObs++) {
	for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
            TString hypoDirName(Form("/%s/Hypo%s", dirName.Data(), fHypoDir[iHypo].c_str()));
	    if(!file->cd(hypoDirName.Data())) {
		cout<<"Directory "<<hypoDirName.Data()<<" not found."<<endl;
		return;
	    }
	    for(Int_t iType = 0; iType < nTypes; iType++) {
		histFileName = Form("%s%s_hypo%i_%s;1",
                                    histPrefix.c_str(),
				    fObs[iObs].c_str(), iHypo,
				    fTypes[iType].c_str());
		fh_kfp[iObs][iHypo][iType] = (TH1F*)gDirectory->Get(histFileName.Data());
		if(!fh_kfp[iObs][iHypo][iType]) {
		    cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
		    continue;
		}
	    }
            //file->cd("/");
	}
    }

    for(Int_t iType = 0; iType < nTypes; iType++) {
        if(!file->cd(dirName.Data())) {
            cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
            return;
        }
        histFileName = Form("fh_kfp_single_pt_%s", fTypes[iType].c_str());
        fh_kfp_pt[iType] = (TH1F*)gDirectory->Get(histFileName.Data());
        if(!fh_kfp_pt[iType]) {
            cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
            continue;
        }
    }

    TCanvas *c[nObs][nHypos];
    THStack *hstack[nObs][nHypos];

    for(Int_t iObs = 0; iObs < nObs; iObs++) {
	for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
	    c[iObs][iHypo] = new TCanvas(Form("c%sHypo%s", fObs[iObs].c_str(), fHypoDir[iHypo].c_str()),
					 Form("c%sHypo%s", fObs[iObs].c_str(), fHypoDir[iHypo].c_str()),
                                         800, 800);
	    hstack[iObs][iHypo] = new THStack(Form("stack%sHypo%s", fObs[iObs].c_str(), fHypoDir[iHypo].c_str()),
					      Form("stack%sHypo%s", fObs[iObs].c_str(), fHypoDir[iHypo].c_str()));
	}
    }

    TLegend* legendType = new TLegend(0.60,0.72,0.95,0.92, "","brNDC");
    TLegendEntry* legentry = NULL;
    for(Int_t iType = 0; iType < nTypes; iType++) {
        if(fh_kfp[0][0][iType]) {
	    if(fTypeDraw[iType]) {
		legentry = legendType->AddEntry(fh_kfp[0][0][iType],
						fTypeLegend[iType].c_str(),"lpf");
            }
        }
    }

    for(Int_t iObs = 0; iObs < nObs; iObs++) {
	for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
	    c[iObs][iHypo]->cd();
	    c[iObs][iHypo]->SetLogy();
	    TH1F *h = NULL;
	    for(Int_t iType = 0; iType <nTypes; iType++) {
		h = fh_kfp[iObs][iHypo][iType];
		if(!h) {
		    continue;
		}
		if(fTypeDraw[iType]) {
		    if(fTypes[iType].find("signal") != std::string::npos) {
                        scaleHistSignal(h, nEvents, sSignalType.Data());
		    } else {
                        scaleHist(h, nEvents);
		    }
                    //h->GetYaxis()->SetTitle("Yield per event");
		    h->SetLineColor(fTypeCols[iType]);
		    hstack[iObs][iHypo]->Add(h);
		    //h->Draw("hist same");
		}
	    }
            hstack[iObs][iHypo]->Draw("nostack hist");
            if(h) {
                hstack[iObs][iHypo]->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
                hstack[iObs][iHypo]->GetXaxis()->SetTitleFont(h->GetXaxis()->GetTitleFont());
                hstack[iObs][iHypo]->GetXaxis()->SetTitleSize(h->GetXaxis()->GetTitleSize());
                hstack[iObs][iHypo]->GetXaxis()->SetLabelFont(h->GetXaxis()->GetLabelFont());
                hstack[iObs][iHypo]->GetXaxis()->SetLabelSize(h->GetXaxis()->GetLabelSize());

                hstack[iObs][iHypo]->GetYaxis()->SetTitle(fObsYAxisTitle[iObs].c_str());
                hstack[iObs][iHypo]->GetYaxis()->SetTitleOffset(1.3);
                hstack[iObs][iHypo]->GetYaxis()->SetTitleFont(h->GetYaxis()->GetTitleFont());
                hstack[iObs][iHypo]->GetYaxis()->SetTitleSize(h->GetYaxis()->GetTitleSize());
                hstack[iObs][iHypo]->GetYaxis()->SetLabelFont(h->GetYaxis()->GetLabelFont());
                hstack[iObs][iHypo]->GetYaxis()->SetLabelSize(h->GetYaxis()->GetLabelSize());
            }
            legendType->Draw("same");
	}
    }

    TCanvas *cPt = new TCanvas("cPtSingle", "cPtSingle", 800, 800);
    THStack *hstackPt = new THStack("stackPt", "stackPt");
    for(Int_t iType = 0; iType < nTypes; iType++) {
        TH1F *h = fh_kfp_pt[iType];
        if(!h) {
            continue;
        }
        if(fTypeDraw[iType]) {
            if(fTypes[iType].find("signal") != std::string::npos) {
                scaleHistSignal(h, nEvents, sSignalType.Data());
            } else {
                scaleHist(h, nEvents);
                h->GetYaxis()->SetTitle("Yield per event");
            }
            h->SetLineColor(fTypeCols[iType]);
            hstackPt->Add(h);
        }
    }
    cPt->cd();
    cPt->SetLogy();
    hstackPt->Draw("nostack hist");
    TH1F *h = fh_kfp_pt[0];
    if(h) {
        hstackPt->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
        hstackPt->GetXaxis()->SetTitleFont(h->GetXaxis()->GetTitleFont());
        hstackPt->GetXaxis()->SetTitleSize(h->GetXaxis()->GetTitleSize());
        hstackPt->GetXaxis()->SetLabelFont(h->GetXaxis()->GetLabelFont());
        hstackPt->GetXaxis()->SetLabelSize(h->GetXaxis()->GetLabelSize());

        hstackPt->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());
        hstackPt->GetYaxis()->SetTitleFont(h->GetYaxis()->GetTitleFont());
        hstackPt->GetYaxis()->SetTitleSize(h->GetYaxis()->GetTitleSize());
        hstackPt->GetYaxis()->SetLabelFont(h->GetYaxis()->GetLabelFont());
        hstackPt->GetYaxis()->SetLabelSize(h->GetYaxis()->GetLabelSize());
    }
    legendType->Draw("same");

    TH1F *fh_kfp_r[nTypes];
    for(Int_t iType = 0; iType < nTypes; iType++) {
        if(!file->cd(dirName.Data())) {
            cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
            return;
        }
        histFileName = Form("fh_kfp_dist_r_first_hits_%s", fTypes[iType].c_str());
        fh_kfp_r[iType] = (TH1F*)gDirectory->Get(histFileName.Data());
        if(!fh_kfp_r[iType]) {
            cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
            continue;
        }
    }
    TCanvas *cDistR = new TCanvas("cDistR", "cDistR", 800, 800);
    THStack *hstackR = new THStack("stackR", "stackR");
    for(Int_t iType = 0; iType < nTypes; iType++) {
        h = fh_kfp_r[iType];
        if(!h) {
            continue;
        }
        if(fTypeDraw[iType]) {
            if(fTypes[iType].find("signal") != std::string::npos) {
                scaleHistSignal(h, nEvents, sSignalType.Data());
            } else {
                scaleHist(h, nEvents);
                h->GetYaxis()->SetTitle("Yield per event");
            }
            h->SetLineColor(fTypeCols[iType]);
            hstackR->Add(h);
        }
    }
    cDistR->cd();
    cDistR->SetLogy();
    hstackR->Draw("nostack hist");
    h = fh_kfp_r[0];
    if(h) {
        hstackR->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
        hstackR->GetXaxis()->SetTitleFont(h->GetXaxis()->GetTitleFont());
        hstackR->GetXaxis()->SetTitleSize(h->GetXaxis()->GetTitleSize());
        hstackR->GetXaxis()->SetLabelFont(h->GetXaxis()->GetLabelFont());
        hstackR->GetXaxis()->SetLabelSize(h->GetXaxis()->GetLabelSize());

        hstackR->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());
        hstackR->GetYaxis()->SetTitleFont(h->GetYaxis()->GetTitleFont());
        hstackR->GetYaxis()->SetTitleSize(h->GetYaxis()->GetTitleSize());
        hstackR->GetYaxis()->SetLabelFont(h->GetYaxis()->GetLabelFont());
        hstackR->GetYaxis()->SetLabelSize(h->GetYaxis()->GetLabelSize());
    }
    legendType->Draw("same");


}

void drawKfpFakePairCut(TFile *file) {
    Int_t nEvents = getNevents(file);
    TString dirName("KFParticle");


    const Int_t nHypos = 5;
    string fHypoDir[nHypos] = { "IsPrimary", "NoVtxGuess", "Pion", "VtxGuess",
                                "IdealVtxGuess"};
    TCanvas *cFakePairCut[nHypos];
    TH1F *fh_kfp_fake_pair_cut[nHypos][2];
    TCanvas *cFakePairCutSb[nHypos];
    TH1F *fh_kfp_sb_fake_pair_cut[nHypos];
    for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
        cFakePairCut[iHypo] = new TCanvas(Form("cFakePairCut%s", fHypoDir[iHypo].c_str()),
                                          Form("cFakePairCut%s", fHypoDir[iHypo].c_str()),
                                          800, 800);

        TString hypoDirName(Form("/%s/Hypo%s", dirName.Data(), fHypoDir[iHypo].c_str()));
        if(!file->cd(hypoDirName.Data())) {
            cout<<"Directory "<<hypoDirName.Data()<<" not found."<<endl;
            return;
        }

        TString histFileName = Form("%s_hypo%i_%s;1",
                                    "fh_kfp_nof",
                                    iHypo,
                                    "fake_pairs");
        fh_kfp_fake_pair_cut[iHypo][0] = (TH1F*)gDirectory->Get(histFileName.Data());
        if(!fh_kfp_fake_pair_cut[iHypo][0]) {
            cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
            continue;
        }
        histFileName = Form("%s_hypo%i_%s;1",
                            "fh_kfp_nof",
                            iHypo,
                            "true_pairs");
        fh_kfp_fake_pair_cut[iHypo][1] = (TH1F*)gDirectory->Get(histFileName.Data());
        if(!fh_kfp_fake_pair_cut[iHypo][1]) {
            cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
            continue;
        }

        cFakePairCutSb[iHypo] = new TCanvas(Form("cFakePairCutSb%s", fHypoDir[iHypo].c_str()),
                                            Form("cFakePairCutSb%s", fHypoDir[iHypo].c_str()),
                                            800, 800);
        fh_kfp_sb_fake_pair_cut[iHypo] = new TH1F(Form("fh_kfp_sb_fake_pair_cut_hypo%i", iHypo),
                                                  Form("fh_kfp_sb_fake_pair_cut_hypo%i", iHypo),
                                                  fh_kfp_fake_pair_cut[iHypo][0]->GetNbinsX(),
                                                  fh_kfp_fake_pair_cut[iHypo][0]->GetXaxis()->GetXmin(),
                                                  fh_kfp_fake_pair_cut[iHypo][0]->GetXaxis()->GetXmax());
        fh_kfp_sb_fake_pair_cut[iHypo]->Divide(fh_kfp_fake_pair_cut[iHypo][1],
                                               fh_kfp_fake_pair_cut[iHypo][0]);

        fh_kfp_fake_pair_cut[iHypo][0]
            ->Scale(1.F / fh_kfp_fake_pair_cut[iHypo][0]->GetBinContent(1));
        //fh_kfp_fake_pair_cut[iHypo][0]->Scale(1.F / (Float_t)nEvents);

        //fh_kfp_fake_pair_cut[iHypo][1]->Scale(1.F / (Float_t)nEvents);
        //Double_t scale_true = fh_kfp_fake_pair_cut[iHypo][0]->GetBinContent(1) /
        //                      fh_kfp_fake_pair_cut[iHypo][1]->GetBinContent(1);
        fh_kfp_fake_pair_cut[iHypo][1]
            ->Scale(1.F / fh_kfp_fake_pair_cut[iHypo][1]->GetBinContent(1));

    }

    TLegend* legendPairCut = new TLegend(0.60,0.83,0.95,0.98, "","brNDC");
    TLegendEntry *legentry = NULL;
    legentry = legendPairCut->AddEntry(fh_kfp_fake_pair_cut[0][0],
                                       "Fake Pairs","lpf");
    legentry = legendPairCut->AddEntry(fh_kfp_fake_pair_cut[0][1],
                                       "True Pairs","lpf");

    for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
        if(fh_kfp_fake_pair_cut[iHypo][0] && fh_kfp_fake_pair_cut[iHypo][1]) {
            cFakePairCut[iHypo]->cd();
            cFakePairCut[iHypo]->SetLogy();
            fh_kfp_fake_pair_cut[iHypo][0]->SetLineColor(kRed);
            fh_kfp_fake_pair_cut[iHypo][0]->GetYaxis()->SetTitle("Fraction");
            fh_kfp_fake_pair_cut[iHypo][0]->GetYaxis()->SetRangeUser(5.e-6, 1.5);
            //fh_kfp_fake_pair_cut[iHypo][0]->GetYaxis()->SetTitle("Yield per event");
            //fh_kfp_fake_pair_cut[iHypo][0]->GetYaxis()->SetTitleOffset(1.1);
            //fh_kfp_fake_pair_cut[iHypo][0]->GetYaxis()->SetRangeUser(1.e-2, 5.e3);
            fh_kfp_fake_pair_cut[iHypo][0]->Draw("hist");
            fh_kfp_fake_pair_cut[iHypo][1]->SetLineColor(kBlue);
            fh_kfp_fake_pair_cut[iHypo][1]->Draw("hist same");
            legendPairCut->Draw("same");
        }
        if(fh_kfp_sb_fake_pair_cut[iHypo]) {
            cFakePairCutSb[iHypo]->cd();
            cFakePairCutSb[iHypo]->SetLogy();

            //cFakePairCut[iHypo]->Update();
            //Float_t rightmax = 1.1*fh_kfp_fake_pair_cut[iHypo][0]->GetBinContent(1);
            // scale    = gPad->GetUymax()/rightmax;

            fh_kfp_sb_fake_pair_cut[iHypo]->SetMarkerStyle(kFullCircle);
            for(Int_t iBin = 1; iBin <= fh_kfp_sb_fake_pair_cut[iHypo]->GetNbinsX(); iBin++) {
                fh_kfp_sb_fake_pair_cut[iHypo]
                    ->GetXaxis()->SetBinLabel(iBin, fh_kfp_fake_pair_cut[iHypo][0]->GetXaxis()->GetBinLabel(iBin));
            }
                fh_kfp_sb_fake_pair_cut[iHypo]
                    ->GetYaxis()->SetTitle("S/B");
            //fh_kfp_sb_fake_pair_cut[iHypo]->SetLineColor(kViolet);
            //fh_kfp_sb_fake_pair_cut[iHypo]->Scale(scale);
            fh_kfp_sb_fake_pair_cut[iHypo]->Draw("pl");

            //TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
            //                          gPad->GetUxmax(),gPad->GetUymax(),
            //                          0,rightmax,510,"+L");
            //axis->Draw();
        }
    }

    TH1F    *fh_kfp_cut[nHypos] = { NULL };
    TCanvas *cKfpCut[nHypos]   = { NULL };
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        return;
    }
    for(Int_t iHypo = 0; iHypo < nHypos; iHypo++) {
        cKfpCut[iHypo] = new TCanvas(Form("cKfpCut%s", fHypoDir[iHypo].c_str()),
                                     Form("cKfpCut%s", fHypoDir[iHypo].c_str()),
                                     800, 800);
        TString histFileName = Form("fh_kfp_pair_cut_%i;1",
                                    iHypo);
        fh_kfp_cut[iHypo] = (TH1F*)gDirectory->Get(histFileName.Data());
        if(!fh_kfp_cut[iHypo]) {
            cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
            continue;
        }
        scaleHist(fh_kfp_cut[iHypo], nEvents);
        cKfpCut[iHypo]->cd();
        fh_kfp_cut[iHypo]->GetYaxis()->SetTitle("Yield per event");
        fh_kfp_cut[iHypo]->Draw("hist");
    }
}

void drawKfpHypos(TFile *file) {
    Int_t nEvents = getNevents(file);
    TString dirName("KFParticle");
    if(!file->cd(dirName.Data())) {
        cout<<"Directory "<<dirName.Data()<<" not found."<<endl;
        return;
    }

    TString histFileName("");

    histFileName = "fh_kfp_correct_hypo;1";
    TH1F *h_kfp_correct_hypo = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_correct_hypo) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    histFileName = "fh_kfp_found_correct_hypo;1";
    TH1F *h_kfp_found_correct_hypo = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_correct_hypo) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    histFileName = "fh_kfp_shift_primsec_m_fake_pairs;1";
    TH1F *h_kfp_shift_m_fake = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_shift_m_fake) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    histFileName = "fh_kfp_shift_primsec_m_prim_pairs;1";
    TH1F *h_kfp_shift_m_prim = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_shift_m_prim) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    histFileName = "fh_kfp_shift_primsec_m_sec_pairs;1";
    TH1F *h_kfp_shift_m_sec = (TH1F*)gDirectory->Get(histFileName.Data());
    if(!h_kfp_shift_m_sec) {
	cout<<"Histogram "<<histFileName.Data()<<" not found"<<endl;
	return;
    }

    TCanvas *cCorrHypo = new TCanvas("cCorrHypo", "Correct Hypothesis", 800, 800);
    cCorrHypo->cd();
    cCorrHypo->SetLogy();
    h_kfp_correct_hypo->Scale(1.F / (Float_t)nEvents);
    h_kfp_correct_hypo->GetXaxis()->SetBinLabel(1, "Primary Electron");
    h_kfp_correct_hypo->GetXaxis()->SetBinLabel(2, "Secondary Electron");
    h_kfp_correct_hypo->GetXaxis()->SetBinLabel(3, "#pi");
    h_kfp_correct_hypo->GetXaxis()->SetBinLabel(4, "Proton");
    h_kfp_correct_hypo->GetXaxis()->SetBinLabel(5, "Other Particle");
    h_kfp_correct_hypo->GetYaxis()->SetTitle("Yield per event");
    h_kfp_correct_hypo->Draw();

    TCanvas *cFoundCorrHypo = new TCanvas("cFoundCorrHypo",
					  "Found Correct Hypothesis", 800, 800);
    cFoundCorrHypo->cd();
    h_kfp_found_correct_hypo->Scale(1.F / (Float_t)(h_kfp_found_correct_hypo->GetEntries())) ;
    h_kfp_found_correct_hypo->GetYaxis()->SetTitle("Fraction");
    h_kfp_found_correct_hypo->GetYaxis()->SetRangeUser(0., 1.);
    h_kfp_found_correct_hypo->Draw();

    TCanvas *cShiftM = new TCanvas("cShiftM", "Mass shift", 800, 800);
    cShiftM->cd();
    cShiftM->SetLogy();
    h_kfp_shift_m_fake->Scale(1.F / (Float_t)nEvents);
    h_kfp_shift_m_fake->SetLineColor(kRed);
    //h_kfp_shift_m_fake->Draw();
    h_kfp_shift_m_prim->Scale(1.F / (Float_t)nEvents);
    h_kfp_shift_m_prim->SetLineColor(kBlue);
    //h_kfp_shift_m_prim->Draw("sames");
    h_kfp_shift_m_sec->Scale(1.F / (Float_t)nEvents);
    h_kfp_shift_m_sec->SetLineColor(kViolet);
    //h_kfp_shift_m_sec->Draw("sames");

    THStack *hStackShiftM = new THStack("hs", "");
    hStackShiftM->Add(h_kfp_shift_m_fake);
    hStackShiftM->Add(h_kfp_shift_m_prim);
    hStackShiftM->Add(h_kfp_shift_m_sec);
    hStackShiftM->Draw("nostack");

    hStackShiftM->GetXaxis()->SetTitle(h_kfp_shift_m_fake->GetXaxis()->GetTitle());
    hStackShiftM->GetXaxis()->SetTitleOffset(1.2);
    hStackShiftM->GetXaxis()->SetTitleFont(h_kfp_shift_m_fake->GetXaxis()->GetTitleFont());
    hStackShiftM->GetXaxis()->SetTitleSize(h_kfp_shift_m_fake->GetXaxis()->GetTitleSize());
    hStackShiftM->GetXaxis()->SetLabelFont(h_kfp_shift_m_fake->GetXaxis()->GetLabelFont());
    hStackShiftM->GetXaxis()->SetLabelSize(h_kfp_shift_m_fake->GetXaxis()->GetLabelSize());

    hStackShiftM->GetYaxis()->SetTitle("Yield per event");
    hStackShiftM->GetYaxis()->SetTitleOffset(1.1);
    hStackShiftM->GetYaxis()->SetTitleFont(h_kfp_shift_m_fake->GetYaxis()->GetTitleFont());
    hStackShiftM->GetYaxis()->SetTitleSize(h_kfp_shift_m_fake->GetYaxis()->GetTitleSize());
    hStackShiftM->GetYaxis()->SetLabelFont(h_kfp_shift_m_fake->GetYaxis()->GetLabelFont());
    hStackShiftM->GetYaxis()->SetLabelSize(h_kfp_shift_m_fake->GetYaxis()->GetLabelSize());

    TLegend* legend = new TLegend(0.62,0.78,0.97,0.93, "","brNDC");
    legend->AddEntry(h_kfp_shift_m_fake, "Fake Pairs","lpf");
    legend->AddEntry(h_kfp_shift_m_prim, "Primary Pairs","lpf");
    legend->AddEntry(h_kfp_shift_m_sec, "Secondary Pairs","lpf");
    legend->Draw("same");

}

void drawHists() {

    TString fileNameMc(Form("%s/analysis.mvd.%s.%s_%sfield_%s.%s.%s.root",
                          sInDir.Data(), "mc", sEnergy.Data(), sField.Data(), sMvd.Data(), sLvm.Data(), sDecay.Data()));
    TString fileNameReco(Form("%s/analysis.mvd.%s.%s_%sfield_%s.%s.%s.root",
                          sInDir.Data(), "reco", sEnergy.Data(), sField.Data(), sMvd.Data(),  sLvm.Data(), sDecay.Data()));
    TString fileNameGamma(Form("%s/analysis.mvd.%s.%s_%sfield_%s.%s.%s.root",
                          sInDir.Data(), "gamma", sEnergy.Data(), sField.Data(), sMvd.Data(),  sLvm.Data(), sDecay.Data()));

    //gROOT->SetStyle("Default");
    gROOT->SetStyle("Plain");

    // n = 1 the name of histogram is printed
    // e = 1 the number of entries
    // m = 1 the mean value m = 2 the mean and mean error values
    // r = 1 the root mean square (RMS) r = 2 the RMS and RMS error
    // u = 1 the number of underflows
    // o = 1 the number of overflows
    // i = 1 the integral of bins s = 1 the skewness s = 2 the skewness and the skewness error
    // k = 1 the kurtosis k = 2 the kurtosis and the kurtosis error
    // 1 -> lower case letter, 2 -> upper case letter
    //gStyle->SetOptStat("mr");
    gStyle->SetOptStat("");
    //gStyle->SetPaintTextFormat("4.2f %%");

    //TGaxis::SetMaxDigits(2);

    TFile *histFileMc = NULL;
    if(!gSystem->IsFileInIncludePath(fileNameMc.Data())) {
        cout<<"Input file "<<fileNameMc.Data()<<" not found."<<endl;
    } else {
        histFileMc   = TFile::Open(fileNameMc.Data());
    }

    if(histFileMc) {
        gFile=histFileMc;
        //drawGammaConv(histFileMc);
        //drawDetHistShort(histFileMc);
        //drawTrackFragsMvdSts(histFileMc);
        //drawMomAccepted(histFileMc);
        //drawThetaAccepted(histFileMc);
        //drawMcHitTopology(histFileMc);
        //drawTrueBgTracks(histFileMc, kTRUE);
    } else {
        cout<<"File "<<fileNameMc<<" not found."<<endl;
    }

    Int_t cut = CbmMvdAnaNames::kCutChi2Primary;

    TFile *histFileReco = NULL;
    if(!gSystem->IsFileInIncludePath(fileNameReco.Data())) {
        cout<<"Input file "<<fileNameReco.Data()<<" not found."<<endl;
    } else {
        histFileReco = TFile::Open(fileNameReco.Data());
    }

    if(histFileReco) {
        gFile = histFileReco;
        //drawHitMismatches(histFileReco);
	//drawHitTopology(histFileReco);
        //drawTrackTopology(histFileReco);
        //drawTrackTypes(histFileReco, cut, kFALSE);
        //drawTrueBgTracks(histFileReco, kFALSE);
	//drawTrueBgCuts(histFileReco);
        //drawCutEfficiencies(histFileReco);
        drawMinStations(histFileReco);
	//drawPrimaryVertexCut(histFileReco);

	//drawKfpInput(histFileReco);
        //drawKfpObservables(histFileReco);
        //drawKfpFakePairCut(histFileReco);
        //drawKfpHypos(histFileReco);
    } else {
        cout<<"File "<<fileNameReco<<" not found."<<endl;
    }


    TFile *histFileGamma = NULL;
    if(!gSystem->IsFileInIncludePath(fileNameGamma.Data())) {
        cout<<"Input file "<<fileNameGamma.Data()<<" not found."<<endl;
    } else {
        histFileGamma = TFile::Open(fileNameGamma.Data());
    }

    if(histFileGamma) {
	gFile = histFileGamma;
        //drawFindGammaInput(histFileGamma);
        //drawFindGammaOutput(histFileGamma);
        //drawFindGammaComb(histFileGamma);
    }
}