// -----------------------------------------------------------------------------
// ----- pl_over_MCPVReco.C                                                -----
// -----                                                                   -----
// ----- created by C. Simon on 2018-10-10                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------


void pl_over_MCPVReco()
{
  gStyle->SetPalette(1);
  gStyle->SetOptStat(kTRUE);
  gStyle->SetOptFit(kTRUE);

  gROOT->cd();
  gROOT->SetDirLevel(1);

  TH1* h1(NULL);
  TH2* h2(NULL);
  TEfficiency* heff(NULL);
  TF1* f1(NULL);
  TString hname;

  TCanvas* can = new TCanvas("can", "MC PV reconstruction overview", 0, 0, 1500, 800);
  can->Divide(6, 3);

  can->cd(1);
  hname = "hPVResXAll";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(7);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaX) [cm]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(13);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }

  can->cd(2);
  hname = "hPVResYAll";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(8);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaY) [cm]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(14);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }

  can->cd(3);
  hname = "hPVResTAll";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(9);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaT) [ns]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(15);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }

  can->cd(4);
  hname = "hPVResXRef";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(10);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaX) [cm]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(16);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }

  can->cd(5);
  hname = "hPVResYRef";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(11);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaY) [cm]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(17);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }

  can->cd(6);
  hname = "hPVResTRef";
  h2 = dynamic_cast<TH2*>(gROOT->FindObjectAny(hname));
  if(h2) {
    h2->Draw("colz");

    can->cd(12);
    h1 = h2->ProjectionY(Form("%s_py", hname.Data()), 1, h2->GetNbinsX());
    gPad->SetLogy();
    h1->Draw();
    h1->Fit("gaus", "Q");
    gPad->Update();
    TPaveStats* stats = dynamic_cast<TPaveStats*>(h1->FindObject("stats"));
    stats->SetOptFit(10001);
    gPad->Modified();

    TH1* h1hs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sHistSigma", hname.Data())));
    if(!h1hs)
    {
      h1hs = new TH1F(Form("%sHistSigma", hname.Data()),
                      Form("#color[602]{hist}/#color[2]{gaus} residual StdDev; %s;  #sigma(#DeltaT) [ns]", h2->GetXaxis()->GetTitle()),
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    TH1* h1gs = dynamic_cast<TH1*>(gROOT->FindObjectAny(Form("%sGausSigma", hname.Data())));
    if(!h1gs)
    {
      h1gs = new TH1F(Form("%sGausSigma", hname.Data()),
                      ";;",
                      h2->GetNbinsX(), h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax());
    }

    for(Int_t iBin = 1; iBin <= h2->GetNbinsX(); iBin++)
    {
      h1 = h2->ProjectionY(Form("%s_py%d", hname.Data(), iBin), iBin, iBin);

      h1hs->SetBinContent(iBin, h1->GetStdDev(1));
      h1hs->SetBinError(iBin, h1->GetStdDevError(1));

      if(h1->GetEntries())
      {
        if(0 == static_cast<Int_t>(h1->Fit("gaus", "Q0")))
        {
          h1gs->SetBinContent(iBin, (h1->GetFunction("gaus"))->GetParameter("Sigma"));
          h1gs->SetBinError(iBin, (h1->GetFunction("gaus"))->GetParError(2));
        }
      }

      delete h1;
    }

    can->cd(18);
    h1hs->SetStats(0);
    h1hs->Draw("PE1");
    h1gs->SetLineColor(kRed);
    h1gs->Draw("PE1same");
  }
  else {
    cout << hname << " not found" << endl;
  }


  can->SaveAs("pl_over_MCPVReco.pdf");
}