/** CbmAnaDielectronTaskDraw.cxx
 * @author Elena Lebedeva <e.lebedeva@gsi.de>
 * @since 2011
 * @version 1.0
 **/

#include "CbmAnaDielectronTaskDraw.h"


#include "CbmDrawHist.h"
#include "CbmHistManager.h"
#include "std/utils/CbmLitUtils.h"

#include <string>
#include <iostream>
#include <iomanip>

#include <boost/assign/list_of.hpp>
#include "CbmAnaPTree.h"

#include "TText.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TMath.h"
#include "TKey.h"
#include "TClass.h"
#include "TF1.h"
#include "TEllipse.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TSystem.h"

ClassImp(CbmAnaDielectronTaskDraw);

using namespace std;
using namespace lit;
using boost::assign::list_of;

void CbmAnaDielectronTaskDraw::DrawHistFromFile(
      const string& fileName,
      const string& outputDir,
      Bool_t useMvd,
      Bool_t drawSig)
{
   SetDefaultDrawStyle();
   fOutputDir = outputDir;
   fUseMvd = useMvd;
   fDrawSignificance = drawSig;

   fPt = new CbmAnaPTree();

   fHM = new CbmHistManager();
   TFile* file = new TFile(fileName.c_str());
   fHM->ReadFromFile(file);

   Int_t fNofEvents = (Int_t)H1("fh_event_number")->GetEntries();
   cout << "File name = " << fileName << endl;
   cout << "Number of events = " << fNofEvents<< endl;

   fHM->ScaleByPattern(".*", 1./fNofEvents);
   SOverBgAll();
   RebinMinvHist();
   DrawPtYDistributionAll();
   DrawPtYEfficiencyAll();
   DrawMomentumDistributionAll();
   DrawMomentumEfficiencyAll();
   DrawMotherPdg();
   DrawSourcesBgPairsAll();
   DrawGammaVertex();
   DrawCutDistributions();
   DrawMinvForEachAnalysisStep();
   DrawMinvSandBgAll();
   DrawBGSourceTracks();
   DrawMinvPtAll();
   DrawBgSourcesVsMomentum();
   SaveCanvasToImage();

   string qaFile = fOutputDir + "lmvm_results.json";
   fPt->Write(qaFile);
}

void CbmAnaDielectronTaskDraw::RebinMinvHist()
{
   int nRebin = 10;
   for (int i = 0; i < CbmAnaLmvmNames::fNofAnaSteps; i++){
      fHM->Rebin("fh_signal_minv_" + CbmAnaLmvmNames::fAnaSteps[i], nRebin);
      fHM->Rebin("fh_bg_minv_" + CbmAnaLmvmNames::fAnaSteps[i], nRebin);
      fHM->Rebin("fh_pi0_minv_" + CbmAnaLmvmNames::fAnaSteps[i], nRebin);
      fHM->Rebin("fh_eta_minv_" + CbmAnaLmvmNames::fAnaSteps[i], nRebin);
   }
}

TH1D* CbmAnaDielectronTaskDraw::H1(
      const string& name)
{
   return (TH1D*) fHM->H1(name);
}

TH2D* CbmAnaDielectronTaskDraw::H2(
      const string& name)
{
   return (TH2D*) fHM->H1(name);
}

TCanvas* CbmAnaDielectronTaskDraw::CreateCanvas(
      const string& name,
      const string& title,
      int width,
      int height)
{
   TCanvas *c = new TCanvas(name.c_str(), title.c_str(), width, height);
   fCanvas.push_back(c);
   return c;
}

void CbmAnaDielectronTaskDraw::DrawEfficiencyOnHist(
      TH1* h1,
      TH1* h2,
      Double_t xPos,
      Double_t yPos)
{
   string effTxt = "";
   if (h2->GetEntries() != 0){
      effTxt = lit::NumberToString<Double_t>(((Double_t)h1->GetEntries()/h2->GetEntries()*100.), 1);
   }
   TText *t = new TText(xPos, yPos, effTxt.c_str());
   t->SetTextSize(0.1);
   t->Draw();
}

void CbmAnaDielectronTaskDraw::DrawTextOnHist(
      const string& text,
      Double_t x1,
      Double_t y1,
      Double_t x2,
      Double_t y2)
{
   TLegend* leg = new TLegend(x1, y1, x2, y2);
   leg->AddEntry(new TH2D(), text.c_str(), "");
   leg->SetFillColor(kWhite);
   leg->SetFillStyle(0);
   leg->SetBorderSize(0);
   leg->Draw();
}

TH2D* CbmAnaDielectronTaskDraw::DivideH2D(
      TH2D* h1,
      TH2D* h2)
{
    Int_t nBinsX = h1->GetNbinsX();
    Double_t minX = h1->GetXaxis()->GetXmin();
    Double_t maxX = h1->GetXaxis()->GetXmax();
    Int_t nBinsY = h1->GetNbinsY();
    Double_t minY = h1->GetYaxis()->GetXmin();
    Double_t maxY = h1->GetYaxis()->GetXmax();

    string hname = string(h1->GetName()) + "_divide";
    TH2D* h3 = new TH2D(hname.c_str(), h1->GetTitle(), nBinsX, minX, maxX, nBinsY, minY, maxY);
    //h1->Sumw2();
    //h2->Sumw2();
    h3->Divide(h1,h2,100.,1.,"B");
    return h3;
}

TH1D* CbmAnaDielectronTaskDraw::DivideH1D(
      TH1D* h1,
      TH1D* h2)
{
    Int_t nBins = h1->GetNbinsX();
    Double_t min = h1->GetXaxis()->GetXmin();
    Double_t max = h1->GetXaxis()->GetXmax();
    string hname = string(h1->GetName()) + "_divide";
    TH1D* h3 = new TH1D(hname.c_str(), h1->GetTitle(), nBins, min, max);
    h3->Sumw2();
    h3->Divide(h1, h2, 100.,1., "B");
    return h3;
}

TH1D* CbmAnaDielectronTaskDraw::CreateSignificanceH1D(
      TH1D* s,
      TH1D* bg,
      const string& name,
      const string& option)
{
    Int_t nBins = s->GetNbinsX();
    Double_t xmin = s->GetXaxis()->GetXmin();
    Double_t xmax = s->GetXaxis()->GetXmax();
    TH1D* hsig = new TH1D(name.c_str(), name.c_str(), nBins, xmin, xmax);
    hsig->GetXaxis()->SetTitle(s->GetXaxis()->GetTitle());

    Double_t sumSignal = 0., sumBg = 0., significance = 0.;
    // "right" - one wants to reject right part of the histogram.
    // value > cut -> reject
    if(option == "right"){
        for (Int_t i = 1; i <= nBins; i++) {
            sumSignal = s->Integral(1, i, "width");
            sumBg = bg->Integral(1, i, "width");
            Double_t prov = TMath::Sqrt(sumSignal+sumBg);
            if (prov != 0. ) significance = sumSignal / prov;
            else significance = 0.;
            hsig->SetBinContent(i, significance);
        }
    // "left" - one wants to reject left part of the histogram.
    // value < cut -> reject
    } else if (option == "left"){
        for (Int_t i = nBins; i >= 1; i--) {
            sumSignal = s->Integral(i, nBins,"width");
            sumBg = bg->Integral(i, nBins,"width");
            Double_t prov = TMath::Sqrt(sumSignal+sumBg);
            if (prov != 0. ) significance = sumSignal / prov;
            else significance = 0.;
            hsig->SetBinContent(i, significance);
        }
    }
    return hsig;
}

TH2D* CbmAnaDielectronTaskDraw::CreateSignificanceH2D(
      TH2D* signal,
      TH2D* bg,
      const string& name,
      const string& title)
{

   Double_t xmin = 1.0;
   Double_t xmax = 5.0;
   Double_t ymin = 1.0;
   Double_t ymax = 5.0;
   Double_t delta = 0.1;
   Int_t nStepsX = (Int_t)( (xmax - xmin) / delta );
   Int_t nStepsY = (Int_t)( (ymax - ymin) / delta );
   Int_t nBinsX = signal->GetNbinsX();
   Int_t nBinsY = signal->GetNbinsY();

   TH2D* hsig = new TH2D(name.c_str(), title.c_str(), nStepsX, xmin, xmax, nStepsY, ymin, ymax);

   Double_t sumSignal = 0;
   Double_t sumBg = 0;

   Int_t binX = 1;
   for (Double_t xcut = xmin; xcut <= xmax; xcut+=delta, binX++) {
      Int_t binY = 1;
      cout << "x " << xcut << endl;
      for (Double_t ycut = ymin; ycut <= ymax; ycut+=delta, binY++){
         sumSignal = 0;
         sumBg = 0;
         for (Int_t ix = 1; ix <= nBinsX; ix++){
            for (Int_t iy = 1; iy <=nBinsY; iy++){
               Double_t xcenter = signal->GetXaxis()->GetBinCenter(ix);
               Double_t ycenter = signal->GetYaxis()->GetBinCenter(iy);
               Double_t val = -1 * (ycut/xcut)*xcenter + ycut;

               if (!(xcenter < xcut && ycenter < val)) {
                  sumSignal += signal->GetBinContent(ix,iy);
                  sumBg += bg->GetBinContent(ix,iy);
               }
            }
         }
         Double_t prov = TMath::Sqrt(sumSignal+sumBg);
         Double_t significance = 0.;
         if (prov != 0) significance = sumSignal / prov;
         hsig->SetBinContent(binX, binY, significance);
      }
   }
   return hsig;
}

void CbmAnaDielectronTaskDraw::SOverBg(
      AnalysisSteps step)
{
   TH1D* s = H1("fh_signal_minv_" + CbmAnaLmvmNames::fAnaSteps[step]);
   TH1D* bg = H1("fh_bg_minv_" + CbmAnaLmvmNames::fAnaSteps[step]);
   TH2D* pty = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[step]);
   TH2D* ptymc = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[0]);

   TH1D* sClone = (TH1D*)s->Clone();
   sClone->Fit("gaus", "Q");

   Double_t mean = sClone->GetFunction("gaus")->GetParameter("Mean");
   Double_t sigma = sClone->GetFunction("gaus")->GetParameter("Sigma");

   Int_t minInd = s->FindBin(mean - 2.*sigma);
   Int_t maxInd = s->FindBin(mean + 2.*sigma);

   Double_t sumSignal = 0.;
   Double_t sumBg = 0.;
   for (Int_t i = minInd + 1; i <= maxInd - 1; i++){
      sumSignal += s->GetBinContent(i);
      sumBg += bg->GetBinContent(i);
   }
   fPt->Put("sbg_"+CbmAnaLmvmNames::fAnaSteps[step], sumSignal/sumBg);
   fPt->Put("eff_"+CbmAnaLmvmNames::fAnaSteps[step], (Double_t)pty->GetEntries()/ptymc->GetEntries()*100.);
   fPt->Put("signal_minv_mean_"+CbmAnaLmvmNames::fAnaSteps[step], 1000.*mean);
   fPt->Put("signal_minv_rms_"+CbmAnaLmvmNames::fAnaSteps[step], 1000.*sigma);
}

void CbmAnaDielectronTaskDraw::SOverBgAll()
{
   TCanvas *c1 = CreateCanvas("lmvm_signal_fitting", "lmvm_signal_fitting", 600, 600);
   SOverBg(kReco);
   SOverBg(kChi2Prim);
   SOverBg(kElId);
   SOverBg(kGammaCut);
   if (fUseMvd) SOverBg(kMvd1Cut);
   if (fUseMvd) SOverBg(kMvd2Cut);
   SOverBg(kStCut);
   SOverBg(kTtCut);
   SOverBg(kPtCut);
}

void CbmAnaDielectronTaskDraw::DrawPtYDistribution(
      int step)
{
   TH2D* h = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[step]);
   TH2D* hmc = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[0]);

   DrawH2(h, kLinear, kLinear, kLinear, "COLZ");
   DrawEfficiencyOnHist(h, hmc, 0.2, 1.8);
   DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.78, 0.70, 0.9);
}

void CbmAnaDielectronTaskDraw::DrawPtYDistributionAll()
{
   Int_t hi = 1;
   TCanvas *c = CreateCanvas("lmvm_pty", "lmvm_pty", 750, 1000);
   c->Divide(3, 4);
   for (int step = 0; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c->cd(hi++);
      DrawPtYDistribution(step);
   }
}

void CbmAnaDielectronTaskDraw::DrawPtYEfficiency(
      int step)
{
   TH2D* h = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[step]);
   // efficiency is normalized to the previous step (step - 1)
   TH2D* hmc = H2("fh_signal_pty_" + CbmAnaLmvmNames::fAnaSteps[step - 1]);

   TH2D* eff = DivideH2D(h, hmc);
   eff->GetZaxis()->SetTitle("Efficiency [%]");
   DrawH2(eff);
   DrawEfficiencyOnHist(h, hmc, 0.2, 1.8);
   DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.78, 0.70, 0.9);
}

void CbmAnaDielectronTaskDraw::DrawPtYEfficiencyAll()
{
   Int_t hi = 1;
   TCanvas *c = CreateCanvas("lmvm_pty_efficiency", "lmvm_pty_efficiency", 750, 1000);
   c->Divide(3,4);
   for (int step = kAcc; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c->cd(hi++);
      DrawPtYEfficiency(step);
   }
}

void CbmAnaDielectronTaskDraw::DrawMomentumDistributionAll()
{
    TCanvas *c = CreateCanvas("lmvm_signal_momentum_distribution", "lmvm_signal_momentum_distribution", 600, 600);
    Draw1DHistoForEachAnalysisStep("fh_signal_mom", true);
}

void CbmAnaDielectronTaskDraw::DrawMomentumEfficiencyAll()
{
//EFFICIENCY vs. MOMENTUM
/*    TCanvas *c5 = CreateCanvas("signal_momentum_efficiency","signal_momentum_efficiency", 600, 600);
    Draw1DHistoForEachAnalysisStep(
          NULL,
          NULL, //DivideHisto1D(H1("fh_acc_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_reco_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_chi_prim_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_el_id_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_gammacut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_dstscut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_dsts2cut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_stcut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_ttcut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_ptcut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_anglecut_signal_mom"), H1("fh_mc_signal_mom")),
          DivideHisto1D(H1("fh_apmcut_signal_mom"), H1("fh_mc_signal_mom")),
          false
    );*/
}

void CbmAnaDielectronTaskDraw::DrawMotherPdg()
{
   TCanvas *c = CreateCanvas("lmvm_mother_pdg", "lmvm_mother_pdg", 500, 500);
   DrawH1(list_of(H1("fh_mc_mother_pdg"))(H1("fh_acc_mother_pdg")),
         list_of("MC")("acc"), kLinear, kLog, true, 0.7, 0.7, 0.99, 0.99);
}

void CbmAnaDielectronTaskDraw::Draw1DSourceTypes(
      const string& hName)
{
   vector<TH1*> h;
   vector<string> hLegend;
   for (int i = 0; i < CbmAnaLmvmNames::fNofSourceTypes; i++){
      string fullName = hName+"_"+CbmAnaLmvmNames::fSourceTypes[i];
      h.push_back( H1(fullName) );
      h[i]->SetLineWidth(2);
      h[i]->SetLineColor(CbmAnaLmvmNames::fSourceTypesColor[i]);
      h[i]->Scale(1. / h[i]->Integral());
      hLegend.push_back( CbmAnaLmvmNames::fSourceTypesLatex[i] );
   }
   DrawH1(h, hLegend, kLinear, kLog, true, 0.90, 0.7, 0.99, 0.99);
}

void CbmAnaDielectronTaskDraw::Draw1DCut(
      const string& hName,
      const string& sigOption)
{
   Int_t w = 600;
   Int_t h = 600;
   if (fDrawSignificance) w = 1200;
   TCanvas *c = CreateCanvas( ("lmvm_" + hName).c_str(), ("lmvm_" + hName).c_str(), w, h);
   if (fDrawSignificance){
      c->Divide(2,1);
      c->cd(1);
   }
   Draw1DSourceTypes(hName);
   if (fDrawSignificance){
      c->cd(2);
      string sName = hName+"_"+CbmAnaLmvmNames::fSourceTypes[kSignal];
      string bgName = hName+"_"+CbmAnaLmvmNames::fSourceTypes[kBg];
      TH1D* sign = CreateSignificanceH1D(H1(sName), H1(bgName), hName+"_significance", sigOption);
      DrawH1(sign);
   }
}

void CbmAnaDielectronTaskDraw::DrawCutDistributions()
{
   Draw1DCut("fh_richann", "left");
   Draw1DCut("fh_trdann", "left");
   Draw2DCut("fh_tofm2");

   Draw1DCut("fh_chi2prim", "right");
   Draw1DCut("fh_pt", "left");
   Draw1DCut("fh_mom", "left");
   Draw1DCut("fh_chi2sts", "right");

   Draw2DCut("fh_stcut");
   Draw2DCut("fh_ttcut");
   Draw2DCut("fh_apmcut");
   Draw2DCut("fh_apcut");
   if (fUseMvd) {
      Draw2DCut("fh_mvd1cut");
      Draw2DCut("fh_mvd2cut");
   }
}


void CbmAnaDielectronTaskDraw::DrawSourcesBgPairs(
      int step,
      bool inPercent)
{
   TH2D* h = (TH2D*)H2("fh_source_pair_" + CbmAnaLmvmNames::fAnaSteps[step])->Clone();
   gStyle->SetPaintTextFormat("4.1f");
   string labels[3] = {"#gamma", "#pi^{0}", "oth"};
   for (Int_t i = 1; i <= 3; i++){
      h->GetYaxis()->SetBinLabel(i, labels[i-1].c_str());
      h->GetXaxis()->SetBinLabel(i, labels[i-1].c_str());
   }
   h->SetMarkerColor(0);
   h->SetMarkerSize(3);
   if (inPercent) {
      h->Scale(100. / h->Integral());
      h->GetZaxis()->SetTitle("[%]");
   } else {
      h->Scale(100.);
      h->GetZaxis()->SetTitle("Number of pairs/event x10^{-2}");
   }
   DrawH2(h, kLinear, kLinear, kLinear, "text COLZ");
   h->GetXaxis()->SetLabelSize(0.1);
   h->GetYaxis()->SetLabelSize(0.1);
   DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.90, 0.70, 0.99);
}

void CbmAnaDielectronTaskDraw::DrawSourcesBgPairsAll()
{
   Int_t hi = 1;
   TCanvas *c1 = CreateCanvas("lmvm_bg_sources_pair_abs", "lmvm_bg_sources_pair_abs", 900, 900);
   c1->Divide(3,3);
   for (int step = kReco; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c1->cd(hi++);
      DrawSourcesBgPairs(step, false);
   }

   hi = 1;
   TCanvas *c2 = CreateCanvas("lmvm_bg_sources_pair_percent", "lmvm_bg_sources_pair_percent", 900, 900);
   c2->Divide(3,3);
   for (int step = kReco; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c2->cd(hi++);
      DrawSourcesBgPairs(step, true);
   }
}

void CbmAnaDielectronTaskDraw::Draw2DCut(
      const string& hist)
{
   string sName = hist + "_"+ CbmAnaLmvmNames::fSourceTypes[kSignal];
   string bgName = hist + "_"+ CbmAnaLmvmNames::fSourceTypes[kBg];
   TCanvas *c = CreateCanvas(("lmvm_" + hist).c_str(), ("lmvm_" + hist).c_str(), 1200, 600);
   c->Divide(4,2);
   c->cd(1);
   DrawH2(H2( sName ));
   DrawTextOnHist(CbmAnaLmvmNames::fSourceTypesLatex[kSignal], 0.5, 0.89, 0.6, 0.99);
   c->cd(2);
   DrawH2(H2( bgName ));
   DrawTextOnHist(CbmAnaLmvmNames::fSourceTypesLatex[kBg], 0.5, 0.89, 0.6, 0.99);
   c->cd(3);
   DrawH2(H2( hist + "_"+ CbmAnaLmvmNames::fSourceTypes[kPi0] ));
   DrawTextOnHist(CbmAnaLmvmNames::fSourceTypesLatex[kPi0], 0.5, 0.89, 0.6, 0.99);
   c->cd(4);
   DrawH2(H2( hist + "_"+ CbmAnaLmvmNames::fSourceTypes[kGamma] ));
   DrawTextOnHist(CbmAnaLmvmNames::fSourceTypesLatex[kGamma], 0.5, 0.89, 0.6, 0.99);
   // Create and draw X projection
   c->cd(5);
   TH1D* pxS = H2(sName)->ProjectionX();
   TH1D* pxBg = H2(bgName)->ProjectionX();
   DrawH1(list_of(pxS)(pxBg), list_of("S")("BG"), kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99);
   // Create and draw Y projection
   c->cd(6);
   TH1D* signX = CreateSignificanceH1D(pxS, pxBg, hist + "_px_significance", "left");
   DrawH1(signX);
   c->cd(7);
   TH1D* pyS = H2(sName)->ProjectionY();
   TH1D* pyBg = H2(bgName)->ProjectionY();
   DrawH1(list_of(pyS)(pyBg), list_of("S")("BG"), kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99);
   c->cd(8);
   TH1D* signY = CreateSignificanceH1D(pyS, pyBg, hist + "_py_significance", "left");
   DrawH1(signY);
   //c->cd(9);
  //TH2D* fh_significance = CalculateSignificance2D(fh_stcut_signal, fh_stcut_bg, "stcut_2dsignificance", "significance");
  //fh_significance->Draw("COLZ");
}

void CbmAnaDielectronTaskDraw::DrawGammaVertex()
{
   TCanvas *c = CreateCanvas("lmvm_gamma_vertex","lmvm_gamma_vertex", 1200, 400);
   c->Divide(3,1);
   c->cd(1);
   DrawH2(H2("fh_mc_vertex_gamma_xz"));
   gPad->SetLogz(true);
   c->cd(2);
   DrawH2(H2("fh_mc_vertex_gamma_yz"));
   gPad->SetLogz(true);
   c->cd(3);
   DrawH2(H2("fh_mc_vertex_gamma_xy"));
   gPad->SetLogz(true);
}

void CbmAnaDielectronTaskDraw::Draw1DHistoForEachAnalysisStep(
      const string& hist,
      Bool_t logy)
{
   string drOpt = ""; // First histogram will be drawn without options.
   Double_t min = 9999999999.;
   TH1D* firstHisto = NULL;
   TLegend* leg = new TLegend(0.80, 0.32, 0.99, 0.99);
   for (int step = 0; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      string fullName = hist + "_" + CbmAnaLmvmNames::fAnaSteps[step];
      TH1D* h = H1(fullName);
      if (h == NULL || h->GetEntries() <= 0) continue;
      leg->AddEntry(h, CbmAnaLmvmNames::fAnaStepsLatex[step].c_str(), "l");
      int color = CbmAnaLmvmNames::fAnaStepsColor[step];
      DrawH1( h, kLinear, kLinear, drOpt, color);
      if (firstHisto == NULL) firstHisto = h;
      drOpt = "same"; // If the histogram is drawn, next histograms must be drawn with "same" option.
      if (min > h->GetMinimum()) { min = h->GetMinimum(); }
   }
   if (min <= 0.0) min = 1e-7;
   if (firstHisto != NULL) firstHisto->SetMinimum(min);

   leg->SetFillColor(kWhite);
   leg->Draw();

   gPad->SetGridx(true);
   gPad->SetGridy(true);
   gPad->SetLogy(logy);
}

void CbmAnaDielectronTaskDraw::DrawMinvForEachAnalysisStep()
{
   TCanvas *c1 = CreateCanvas("lmvm_minv_for_each_analysis_step_s_bg", "lmvm_minv_for_each_analysis_step_s_bg", 1200, 600);
   c1->Divide(2,1);
   c1->cd(1);
   Draw1DHistoForEachAnalysisStep("fh_signal_minv", true);
   c1->cd(2);
   Draw1DHistoForEachAnalysisStep("fh_bg_minv", true);

   TCanvas *c2 = CreateCanvas("lmvm_minv_for_each_analysis_step_pi0_eta", "lmvm_minv_for_each_analysis_step_pi0_eta", 1200, 600);
   c2->Divide(2,1);
   c2->cd(1);
   Draw1DHistoForEachAnalysisStep("fh_pi0_minv", true);
   c2->cd(2);
   Draw1DHistoForEachAnalysisStep("fh_eta_minv", true);
}

void CbmAnaDielectronTaskDraw::DrawMinvSandBg(
      int step)
{
   TH1* s1 = H1("fh_signal_minv_" + CbmAnaLmvmNames::fAnaSteps[step]);
   TH1* bg1 = H1("fh_bg_minv_" + CbmAnaLmvmNames::fAnaSteps[step]);

   TH1D* s = (TH1D*)s1->Clone();
   TH1D* bg = (TH1D*)bg1->Clone();
   TH1D* sbg = (TH1D*)bg->Clone();
   sbg->Add(s);
   sbg->SetMinimum(1e-8);

   DrawH1(list_of(sbg)(bg)(s), list_of("")("")(""), kLinear, kLog, false, 0,0,0,0);
   s->SetFillColor(kRed);
   s->SetLineColor(kBlack);
   s->SetLineWidth(1);
   s->SetLineStyle(CbmDrawingOptions::MarkerStyle(1));
   bg->SetFillColor(kYellow - 10);
   bg->SetLineColor(kBlack);
   bg->SetLineWidth(1);
   bg->SetLineStyle(CbmDrawingOptions::MarkerStyle(1));
   sbg->SetFillColor(kBlue);
   sbg->SetLineColor(kBlack);
   sbg->SetLineWidth(1);
   sbg->SetLineStyle(CbmDrawingOptions::MarkerStyle(1));
   s->SetMarkerStyle(1);
   bg->SetMarkerStyle(1);
   sbg->SetMarkerStyle(1);

   DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.65, 0.78, 0.85, 0.9);
}

void CbmAnaDielectronTaskDraw::DrawMinvSandBgAll()
{
   Int_t hi = 1;
   TCanvas *c = CreateCanvas("lmvm_minv_both_s_and_bg", "lmvm_minv_both_s_and_bg", 900, 900);
   c->Divide(3,3);
   for (int step = kReco; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c->cd(hi++);
      DrawMinvSandBg(step);
   }
}

void CbmAnaDielectronTaskDraw::DrawBGSourceTracks()
{
   gStyle->SetPaintTextFormat("4.1f");
   TCanvas *c1 = CreateCanvas("lmvm_nof_tracks_per_event", "lmvm_nof_tracks_per_event", 1200, 400);
   c1->Divide(3,1);
   c1->cd(1);
   DrawH1( H1("fh_nof_bg_tracks"), kLinear, kLog, "hist text0");

   c1->cd(2);
   DrawH1( H1("fh_nof_el_tracks"), kLinear, kLog);

   c1->cd(3);
   TH1D* purity = new TH1D("purity","purity;Analysis steps;Purity",
         CbmAnaLmvmNames::fNofAnaSteps, 0., CbmAnaLmvmNames::fNofAnaSteps);
   purity->Divide(H1("fh_nof_bg_tracks"), H1("fh_nof_el_tracks"));
   DrawH1( purity, kLinear, kLog, "hist text0");

   TCanvas *c2 = CreateCanvas("lmvm_nof_tracks_2d_abs", "lmvm_nof_tracks_2d_abs", 900, 600);
   TH2D* source_tracks_abs = (TH2D*)H2("fh_source_tracks")->Clone();
   source_tracks_abs->SetStats(false);
   source_tracks_abs->Scale(100);
   source_tracks_abs->GetZaxis()->SetTitle("Tracks per event x10^{-2}");
   source_tracks_abs->GetXaxis()->SetRange(kReco, CbmAnaLmvmNames::fNofAnaSteps);
   DrawH2(source_tracks_abs, kLinear, kLinear, kLog, "text COLZ");

   TCanvas *c3 = CreateCanvas("lmvm_nof_tracks_2d_percent", "lmvm_nof_tracks_2d_percent", 900, 600);
   TH2D* source_tracks_perc = (TH2D*)H2("fh_source_tracks")->Clone();
   source_tracks_perc->SetStats(false);
   Int_t nBinsX = source_tracks_perc->GetNbinsX();
   Int_t nBinsY = source_tracks_perc->GetNbinsY();
   for (Int_t x = 1; x <= nBinsX; x++){
      Double_t scale =100./(H1("fh_nof_bg_tracks")->GetBinContent(x));
      for (Int_t y = 1; y <= nBinsY; y++){
         Double_t val = scale * source_tracks_perc->GetBinContent(x,y);
         source_tracks_perc->SetBinContent(x,y,val);
      }
   }
   source_tracks_perc->GetZaxis()->SetTitle("[%]");
   source_tracks_perc->GetXaxis()->SetLabelSize(0.06);
   source_tracks_perc->GetYaxis()->SetLabelSize(0.06);
   source_tracks_perc->SetMarkerColor(kBlack);
   source_tracks_perc->SetMarkerSize(1.8);
   source_tracks_perc->GetXaxis()->SetRange(kReco, CbmAnaLmvmNames::fNofAnaSteps);
   DrawH2(source_tracks_perc, kLinear, kLinear, kLinear, "text COLZ");

   Int_t ny = 6;
   string yLabels[6] = {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "oth"};
   for (Int_t y = 1; y <= ny; y++){
      source_tracks_perc->GetYaxis()->SetBinLabel(y, yLabels[y-1].c_str());
      source_tracks_abs->GetYaxis()->SetBinLabel(y, yLabels[y-1].c_str());
   }
   SetAnalysisStepLabels(H1("fh_nof_bg_tracks"));
   SetAnalysisStepLabels(H1("fh_nof_el_tracks"));
   SetAnalysisStepLabels(purity);
   SetAnalysisStepLabels(source_tracks_perc);
   SetAnalysisStepLabels(source_tracks_abs);
}

void CbmAnaDielectronTaskDraw::SetAnalysisStepLabels(
      TH1* h)
{
   h->GetXaxis()->SetLabelSize(0.06);
   for (Int_t x = 1; x <= CbmAnaLmvmNames::fNofAnaSteps; x++){
      h->GetXaxis()->SetBinLabel(x, CbmAnaLmvmNames::fAnaStepsLatex[x-1].c_str());
   }
}

void CbmAnaDielectronTaskDraw::DrawMinvPtAll()
{
   Int_t hi = 1;
   TCanvas *c = CreateCanvas("lmvm_fh_signal_minv_pt", "lmvm_fh_signal_minv_pt", 750, 1000);
   c->Divide(3, 4);
   for (int step = 0; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c->cd(hi++);

      TH2D* h = H2("fh_signal_minv_pt_" + CbmAnaLmvmNames::fAnaSteps[step]);
      DrawH2(h, kLinear, kLinear, kLinear, "COLZ");
      DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.78, 0.70, 0.9);
   }
}

void CbmAnaDielectronTaskDraw::DrawBgSourcesVsMomentum()
{
   int hi = 1;
   TCanvas *c1 = CreateCanvas("lmvm_fh_source_mom","lmvm_fh_source_mom", 900, 900);
   c1->Divide(3,3);
   for (Int_t step = kReco; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c1->cd (hi++);
      Draw1DSourceTypes("fh_source_mom_" + CbmAnaLmvmNames::fAnaSteps[step]);
      DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.90, 0.70, 0.99);
   }

   hi = 1;
   TCanvas *c2 = CreateCanvas("lmvm_fh_opening_angle","lmvm_fh_opening_angle", 900, 900);
   c2->Divide(3,3);
   for (Int_t step = kReco; step < CbmAnaLmvmNames::fNofAnaSteps; step++){
      if ( !fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) continue;
      c2->cd (hi++);
      Draw1DSourceTypes("fh_opening_angle_" + CbmAnaLmvmNames::fAnaSteps[step]);
      DrawTextOnHist(CbmAnaLmvmNames::fAnaStepsLatex[step], 0.50, 0.90, 0.70, 0.99);
   }
}

void CbmAnaDielectronTaskDraw::SaveCanvasToImage()
{
   for (int i = 0; i < fCanvas.size(); i++){
      lit::SaveCanvasAsImage(fCanvas[i], fOutputDir);
   }
}