#!/usr/bin/env python
# coding: utf-8

# # Introduction
# 
# QA plots for the generic tracking performance

# In[1]:


# imports to write dynamic markdown contents
import os
from IPython.display import display, Markdown, Latex
from IPython.display import HTML


# In[2]:


# turn off/on code for the result HTML page

display(Markdown('*For the result HTML page:* '))
    
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 $('div.jp-CodeMirrorEditor').hide(); // for newer versions of nbviewer.org
 } else {
 $('div.input').show();
 $('div.jp-CodeMirrorEditor').show(); // for newer versions of nbviewer.org
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')


# In[3]:


import os.path

# readme file of the macros, available if run under JenkinsCI
# https://github.com/sPHENIX-Collaboration/utilities/blob/master/jenkins/built-test/test-tracking-qa.sh 

macro_markdown = 'Fun4All-macros-README.md'


if os.path.isfile(macro_markdown) :
    with open(macro_markdown, 'r') as file:
        display(Markdown(file.read()))


# ## `pyROOT` env check

# In[4]:


import ROOT

OFFLINE_MAIN = os.getenv("OFFLINE_MAIN")
if OFFLINE_MAIN is not None:
    display(Markdown(f"via sPHENIX software distribution at `{OFFLINE_MAIN}`"))


# ## Plotting source code

# In[5]:


import subprocess

try:
    git_url = \
        subprocess.run(['git','remote','get-url','origin'], stdout=subprocess.PIPE)\
        .stdout.decode('utf-8').strip()\
        .replace('git@github.com:','https://github.com/')

    display(Markdown(f"View the source code repository at {git_url}"))
except: # catch *all* exceptions
    #     well do nothing
    pass


# ## JenkinsCI information (if available)

# In[6]:


display(Markdown('Some further details about the QA run, if executed under the Jenkins CI:'))

checkrun_repo_commit = os.getenv("checkrun_repo_commit")
if checkrun_repo_commit is not None:
    display(Markdown(f"* The commit being checked is {checkrun_repo_commit}"))
        
ghprbPullLink = os.getenv("ghprbPullLink")
if ghprbPullLink is not None:
    display(Markdown(f"* Link to the pull request: {ghprbPullLink}"))

BUILD_URL =  os.getenv("BUILD_URL")
if BUILD_URL is not None:
    display(Markdown(f"* Link to the build: {BUILD_URL}"))

git_url_macros =  os.getenv("git_url_macros")
sha_macros =  os.getenv("sha_macros")
if git_url_macros is not None:
    display(Markdown(f"* Git repo for macros: {git_url_macros} , which merges `{sha_macros}` and the QA tracking branch"))

RUN_ARTIFACTS_DISPLAY_URL = os.getenv("RUN_ARTIFACTS_DISPLAY_URL")
if RUN_ARTIFACTS_DISPLAY_URL is not None:
    display(Markdown(f"* Download the QA ROOT files: {RUN_ARTIFACTS_DISPLAY_URL}"))

JENKINS_URL  = os.getenv("JENKINS_URL")
if JENKINS_URL is not None:
    display(Markdown(f"Automatically generated by [sPHENIX Jenkins continuous integration]({JENKINS_URL}) [![sPHENIX](https://raw.githubusercontent.com/sPHENIX-Collaboration/utilities/master/jenkins/material/sphenix-logo-white-bg-72p.png)](https://www.sphenix.bnl.gov/web/) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; [![jenkins.io](https://raw.githubusercontent.com/sPHENIX-Collaboration/utilities/master/jenkins/material/jenkins_logo_title-72p.png)](https://jenkins.io/)"))


# 
# # Initialization

# In[7]:


get_ipython().run_cell_magic('cpp', '-d', '\n#include "QA_Draw_Utility.C"\n\n#include <sPhenixStyle.C>\n\n#include <TFile.h>\n#include <TLine.h>\n#include <TString.h>\n#include <TTree.h>\n#include <cassert>\n#include <cmath>\n')


# In[8]:


get_ipython().run_cell_magic('cpp', '', '\nSetsPhenixStyle();\nTVirtualFitter::SetDefaultFitter("Minuit2");\n\n// test sPHENIX lib load\n// gSystem->Load("libg4eval.so");\n\n// test libs\n// gSystem->ListLibraries();\n')


# In[9]:


get_ipython().run_line_magic('jsroot', 'on')


# ## Inputs and file checks

# In[10]:


qa_file_name_new = os.getenv("qa_file_name_new")
if qa_file_name_new is None:
    qa_file_name_new = "G4sPHENIX_test-tracking-low-occupancy-qa_Event100_Sum10_qa.root"
    display(Markdown(f"`qa_file_name_new` env not set. use the default `qa_file_name_new={qa_file_name_new}`"))

qa_file_name_ref = os.getenv("qa_file_name_ref")
if qa_file_name_ref is None:
    qa_file_name_ref = "reference/G4sPHENIX_test-tracking-low-occupancy-qa_Event100_Sum10_qa.root"
    display(Markdown(f"`qa_file_name_ref` env not set. use the default `qa_file_name_ref={qa_file_name_ref}`"))
elif qa_file_name_ref == 'None':
    qa_file_name_ref = None
    display(Markdown(f"`qa_file_name_ref` = None and we are set to not to use the reference histograms"))
    


# In[11]:


# qa_file_new = ROOT.TFile.Open(qa_file_name_new);

# assert qa_file_new.IsOpen()
# qa_file_new.ls()
display(Markdown(f"Opening QA file at `{qa_file_name_new}`"))
ROOT.gInterpreter.ProcessLine(f"TFile *qa_file_new = new TFile(\"{qa_file_name_new}\");")
ROOT.gInterpreter.ProcessLine(f"const char * qa_file_name_new = \"{qa_file_name_new}\";")

if qa_file_name_ref is not None:
#     qa_file_ref = ROOT.TFile.Open(qa_file_name_ref);

#     assert qa_file_ref.IsOpen()
    display(Markdown(f"Opening QA reference file at `{qa_file_name_ref}`"))
    ROOT.gInterpreter.ProcessLine(f"TFile *qa_file_ref = new TFile(\"{qa_file_name_ref}\");")
    ROOT.gInterpreter.ProcessLine(f"const char *  qa_file_name_ref = \"{qa_file_name_ref}\";")
else:
    ROOT.gInterpreter.ProcessLine(f"TFile *qa_file_ref = nullptr;")
    ROOT.gInterpreter.ProcessLine(f"const char *  qa_file_name_ref = nullptr;")


# In[12]:


get_ipython().run_cell_magic('cpp', '', '\nif (qa_file_new == nullptr) \n{\n    cout <<"Error, can not open QA root file"<<qa_file_name_new<<endl;\n    exit(1);\n}\n\n// list inputs histograms if needed\n// qa_file_new ->ls();\n\n//TFile *qa_file_ref = NULL;\n//if (qa_file_name_ref)\n//{\n//    qa_file_ref = new TFile(qa_file_name_ref);\n//  \n//    if (qa_file_ref == nullptr) \n//    {\n//        cout <<"Error, can not open QA root file"<<qa_file_name_ref<<endl;\n//        exit(1);\n//    }\n//}\n')


# # Truth matched efficiency and $p_T$ resolution

# ## Match truth track to reco tracks
# 
# Reco track matched with truth track. Please note the efficiency only require truth-reco matching in unique pairs, no cut on the reco quality is required, as defined in [`QAG4SimulationTracking::process_event(PHCompositeNode * )`](https://sphenix-collaboration.github.io/doxygen/de/dbf/classQAG4SimulationTracking.html#ac2aeab5238fcf828f70c61b4976b57fe)

# In[13]:


get_ipython().run_cell_magic('cpp', '', '{\n    const char *hist_name_prefix = "QAG4SimulationTracking";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n    TCanvas *rawResCan = new TCanvas(TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,\n                                     TString("QA_Draw_pTResolution") + TString("_") + hist_name_prefix,\n                                     900,500);\n    rawResCan->Divide(2,2);\n    \n    TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new =\n        (TH2F *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",\n                                             "TH2");\n    assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new);\n    \n    h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetDirectory(nullptr);\n    \n    TGraphErrors *newScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 1);\n    TGraphErrors *newRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new, false, 2);\n    newScale->SetMarkerStyle(20);\n    newScale->SetMarkerColor(kBlack);\n    newScale->SetMarkerSize(1.0);\n    newRes->SetMarkerStyle(20);\n    newRes->SetMarkerColor(kBlack);\n    newRes->SetMarkerSize(1.0);\n    \n    rawResCan->cd(3);\n    newScale->GetXaxis()->SetTitle("Truth p_{T} [GeV]");\n    newScale->GetYaxis()->SetTitle("#LTp_{T}^{reco}/p_{T}^{true}#GT");\n    newScale->GetYaxis()->SetRangeUser(0.9,1.1);\n    newScale->Draw("ap");\n    \n    rawResCan->cd(4);\n    newRes->GetXaxis()->SetTitle("Truth p_{T} [GeV]");\n    newRes->GetYaxis()->SetTitle("#sigma(p_{T}^{reco}/p_{T}^{true})");\n    newRes->GetYaxis()->SetRangeUser(0,0.08);\n    newRes->Draw("ap");\n    \n    rawResCan->cd(2);\n    gPad->SetRightMargin(0.2);\n    gPad->SetLogx();\n    gPad->SetLogz();\n    h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->SetTitle("New pT Spectrum");\n    h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_new->Draw("colz");\n    \n    TLatex newl; \n    newl.SetTextSize(0.05); \n    newl.SetNDC();\n    newl.SetTextColor(kBlack);\n    newl.DrawLatex(0.45,0.96,"New");\n    \n    if (qa_file_ref) {\n      TH2F *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref =\n          (TH2F *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",\n                                               "TH2");\n      assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref);\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetDirectory(nullptr);\n        \n      TGraphErrors *refScale = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 1);\n      TGraphErrors *refRes = FitResolution(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref, false, 2);\n      refScale->SetMarkerStyle(24);\n      refScale->SetMarkerColor(kRed);\n      refScale->SetMarkerSize(1.0);\n      refRes->SetMarkerStyle(24);\n      refRes->SetMarkerColor(kRed);\n      refRes->SetMarkerSize(1.0);\n        \n      TLegend *ptresleg = new TLegend(0.4,0.77,0.5,0.9);\n      ptresleg->AddEntry(refScale,"Reference","P");\n      ptresleg->AddEntry(newScale,"New","P"); \n        \n      rawResCan->cd(3);\n      refScale->Draw("psame");\n      ptresleg->Draw("same");\n       \n      rawResCan->cd(4);\n      refRes->Draw("psame");\n      ptresleg->Draw("same"); \n        \n      rawResCan->cd(1);\n  \n      gPad->SetRightMargin(0.2);\n      gPad->SetLogx();\n      gPad->SetLogz();\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->SetTitle("Reference pT Spectrum");\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen_ref->Draw("colz");\n      TLatex refl; \n      refl.SetTextSize(0.05); \n      refl.SetNDC();\n      refl.SetTextColor(kBlack);\n      refl.DrawLatex(0.45,0.96,"Reference");\n  \n    }\n        \n    rawResCan->Draw();\n    \n   \n}\n\n\n')


# In[14]:


get_ipython().run_line_magic('jsroot', 'off')


# In[15]:


get_ipython().run_cell_magic('cpp', '', '\n{\n  //base histogram from the reco module name \n  const char *hist_name_prefix = "QAG4SimulationTracking";\n  TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatchingOverview") +\n                                TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_TruthMatchingOverview") +\n                                TString("_") + hist_name_prefix,\n                            900, 400);\n  c1->Divide(3, 1);\n  int idx = 1;\n  TPad *p;\n\n  {\n    static const int nrebin = 5;\n\n    p = (TPad *)c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    p->SetGridy();\n\n    TH1 *h_pass =\n        (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_pTGen", "TH1");\n    TH1 *h_norm =\n        (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_pTGen", "TH1");\n    assert(h_norm);\n    assert(h_pass);\n      \n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n\n    h_norm->Rebin(nrebin);\n    h_pass->Rebin(nrebin);\n\n    TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_ratio->GetYaxis()->SetTitle("Reco efficiency");\n    h_ratio->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref) {\n      TH1 *h_pass =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_pTGen", "TH1");\n      TH1 *h_norm =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_pTGen", "TH1");\n      assert(h_norm);\n      assert(h_pass);\n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n      h_norm->Rebin(nrebin);\n      h_pass->Rebin(nrebin);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");\n\n    DrawReference(h_ratio, h_ratio_ref, false);\n  }\n\n  {\n    static const int nrebin = 4;\n\n    p = (TPad *)c1->cd(idx++);\n    c1->Update();\n    // p->SetLogx();\n    p->SetGridy();\n\n    TH1 *h_pass =\n        (TH1 *)qa_file_new->GetObjectChecked(prefix + "nReco_etaGen", "TH1");\n    TH1 *h_norm =\n        (TH1 *)qa_file_new->GetObjectChecked(prefix + "nGen_etaGen", "TH1");\n    assert(h_norm);\n    assert(h_pass);\n\n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n    h_norm->Rebin(nrebin);\n    h_pass->Rebin(nrebin);\n\n    TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n\n    h_ratio->GetXaxis()->SetRangeUser(-1.1, 1.1);\n    h_ratio->GetYaxis()->SetTitle("Reco efficiency");\n    h_ratio->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref) {\n      TH1 *h_pass =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nReco_etaGen", "TH1");\n      TH1 *h_norm =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "nGen_etaGen", "TH1");\n      assert(h_norm);\n      assert(h_pass);\n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n      h_norm->Rebin(nrebin);\n      h_pass->Rebin(nrebin);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    h_ratio->SetTitle(TString(hist_name_prefix) + ": Tracking Efficiency");\n\n    DrawReference(h_ratio, h_ratio_ref, false);\n  }\n\n  {\n    p = (TPad *)c1->cd(idx++);\n    c1->Update();\n    //    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,\n                              "Mean and sigma, p_{T,reco}/p_{T,truth};Truth p_{T} [GeV/c];<p_{T,reco}/p_{T,truth}> #pm #sigma(p_{T,reco}/p_{T,truth})");\n    //gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    //TLine *l = new TLine(0, 1, 50, 1);\n    //l->SetLineColor(kGray);\n    //l->Draw();\n    HorizontalLine( gPad, 1 )->Draw();\n\n    TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =\n        (TH2 *)qa_file_new->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",\n                                             "TH2");\n    assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);\n\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);\n    h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);\n    TGraphErrors *ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen =\n        FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);\n    ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");\n    ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetTitle(\n        "Mean and sigma, p_{T,reco}/p_{T,truth}");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref) {\n      TH2 *h_QAG4SimulationTracking_pTRecoGenRatio_pTGen =\n          (TH2 *)qa_file_ref->GetObjectChecked(prefix + "pTRecoGenRatio_pTGen",\n                                               "TH2");\n      assert(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);\n\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->SetDirectory(nullptr);\n      h_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Rebin2D(16, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_pTRecoGenRatio_pTGen);\n      ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen, h_ratio_ref,\n                  true);\n            \n      SaveGraphError2CSV(ge_QAG4SimulationTracking_pTRecoGenRatio_pTGen, "QAG4SimulationTracking_pTRecoGenRatio_pTGen");\n  }\n\n  //SaveCanvas(c1,\n  //           TString(qa_file_name_new) + TString("_") + TString(c1->GetName()),\n  //           true);\n    \n  c1->Draw();\n}\n')


# # $p_T$ resolution and lineshape

# In[16]:


get_ipython().run_line_magic('jsroot', 'on')


# In[17]:


get_ipython().run_cell_magic('cpp', '', '\n{\n    const char *hist_name_prefix = "QAG4SimulationTracking";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n    \n    \n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("pTRecoGenRatio_pTGen"), "TH2");\n  assert(h_new);\n\n  //  h_new->Rebin(1, 2);\n  //h_new->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref = NULL;\n  if (qa_file_ref)\n  {\n    h_ref = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("pTRecoGenRatio_pTGen"), "TH2");\n    assert(h_ref);\n\n    //    h_ref->Rebin(1, 2);\n    //h_ref->Sumw2();\n    h_ref->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_pTRatio") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(4, 2);\n  int idx = 1;\n  TPad *p;\n\n  vector<pair<double, double>> gpt_ranges{\n      {0, 1},\n      {1, 5},\n      {5, 10},\n      {10, 15},\n      {15, 50}};\n  TF1 *f1 = nullptr;\n  TF1 *fit = nullptr;\n  Double_t sigma = 0;\n  Double_t sigma_unc = 0;\n  char resstr[500];\n  TLatex *res = nullptr;\n  for (auto pt_range : gpt_ranges)\n  {\n    //cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";\n\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new->GetXaxis()->FindBin(pt_range.first + epsilon);\n    const int bin_end = h_new->GetXaxis()->FindBin(pt_range.second - epsilon);\n\n    TH1 *h_proj_new = h_new->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n\n    h_proj_new->GetXaxis()->SetRangeUser(.7, 1.3);\n    h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));\n    h_proj_new->GetXaxis()->SetTitle(TString::Format(\n        "Reco p_{T}/Truth p_{T}"));\n\n    f1 = new TF1("f1", "gaus", -.85, 1.15);\n    h_proj_new->Fit(f1, "mq");\n    sigma = f1->GetParameter(2);\n    sigma_unc = f1->GetParError(2);\n\n    TH1 *h_proj_ref = nullptr;\n    if (h_ref)\n      h_proj_ref =\n          h_ref->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n\n    DrawReference(h_proj_new, h_proj_ref);\n    sprintf(resstr, "#sigma = %.5f #pm %.5f", sigma, sigma_unc);\n    res = new TLatex(0.325, 0.825, resstr);\n    res->SetNDC();\n    res->SetTextSize(0.05);\n    res->SetTextAlign(13);\n    res->Draw();\n  }\n\n // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# In[18]:


get_ipython().run_line_magic('jsroot', 'off')


# In[19]:


get_ipython().run_cell_magic('cpp', '', '\n{\nconst char *hist_name_prefix = "QAG4SimulationTracking";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\n    //cout << "Open new QA file " << qa_file_new->GetName() << endl;\n\n    TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n  if (qa_file_ref)\n  {\n   // cout << "Open ref QA file " << qa_file_ref->GetName() << endl;\n    TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n    \n    \n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_RecoTruthMatching") +\n                                TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_RecoTruthMatching") +\n                                TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(2, 1);\n  int idx = 1;\n  TPad *p;\n\n  {\n    static const int nrebin = 5;\n\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    p->SetGridy();\n\n    TH1 *h_pass =\n        (TH1 *) qa_file_new->GetObjectChecked(prefix + "nGen_pTReco", "TH1");\n    TH1 *h_norm =\n        (TH1 *) qa_file_new->GetObjectChecked(prefix + "nReco_pTReco", "TH1");\n    assert(h_norm);\n    assert(h_pass);\n      \n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n\n    h_norm->Rebin(nrebin);\n    h_pass->Rebin(nrebin);\n\n    TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_ratio->GetYaxis()->SetTitle("Tracking Purity");\n    h_ratio->GetYaxis()->SetRangeUser(-0, 1.1);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH1 *h_pass =\n          (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nGen_pTReco", "TH1");\n      TH1 *h_norm =\n          (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nReco_pTReco", "TH1");\n      assert(h_norm);\n      assert(h_pass);\n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n      h_norm->Rebin(nrebin);\n      h_pass->Rebin(nrebin);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    h_ratio->SetTitle("Tracking Purity (matched truth-reco pairs)");\n\n    DrawReference(h_ratio, h_ratio_ref, false);\n  }\n\n  {\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    //    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,\n                              "Mean and sigma p_{Tmatched}/p_{Treco};Reco p_{T} [GeV/c];<p_{T,matched}/p_{T,reco}> #pm #sigma(p_{T,matched}/p_{T,reco})");\n    // gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    // TLine *l = new TLine(0, 1, 50, 1);\n    // l->SetLineColor(kGray);\n    // l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n        (TH2 *) qa_file_new->GetObjectChecked(\n            prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");\n    assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);\n    h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);\n\n    TGraphErrors *ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n        FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n    ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");\n    ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetTitle(\n        "Mean and sigma p_{Tmatched}/p_{Treco}");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n          (TH2 *) qa_file_ref->GetObjectChecked(\n              prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");\n      assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);\n\n      h_ratio_ref =\n          FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n      ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco,\n                  h_ratio_ref, true);\n\n  }\n\n    c1->Draw();\n}\n')


# In[20]:


get_ipython().run_cell_magic('cpp', '', '\n{\nconst char *hist_name_prefix = "QAG4SimulationTracking";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\n    //cout << "Open new QA file " << qa_file_new->GetName() << endl;\n\n    TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n  if (qa_file_ref)\n  {\n   // cout << "Open ref QA file " << qa_file_ref->GetName() << endl;\n    TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n    \n    \n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_RecoTruthMatching") +\n                                TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_RecoTruthMatching") +\n                                TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(2, 1);\n  int idx = 1;\n  TPad *p;\n\n  {\n    static const int nrebin = 5;\n\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    p->SetGridy();\n\n    TH1 *h_pass =\n        (TH1 *) qa_file_new->GetObjectChecked(prefix + "nGen_pTReco", "TH1");\n    TH1 *h_norm =\n        (TH1 *) qa_file_new->GetObjectChecked(prefix + "nReco_pTReco", "TH1");\n    assert(h_norm);\n    assert(h_pass);\n      \n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n\n    h_norm->Rebin(nrebin);\n    h_pass->Rebin(nrebin);\n\n    TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_ratio->GetYaxis()->SetTitle("Tracking Purity");\n    h_ratio->GetYaxis()->SetRangeUser(-0, 1.1);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH1 *h_pass =\n          (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nGen_pTReco", "TH1");\n      TH1 *h_norm =\n          (TH1 *) qa_file_ref->GetObjectChecked(prefix + "nReco_pTReco", "TH1");\n      assert(h_norm);\n      assert(h_pass);\n      h_norm->SetDirectory(nullptr);\n      h_pass->SetDirectory(nullptr);\n      h_norm->Rebin(nrebin);\n      h_pass->Rebin(nrebin);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    h_ratio->SetTitle("Tracking Purity (matched truth-reco pairs)");\n\n    DrawReference(h_ratio, h_ratio_ref, false);\n  }\n\n  {\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    //    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0, .9, 50, 1.1,\n                              "Mean and sigma p_{Tmatched}/p_{Treco};Reco p_{T} [GeV/c];<p_{T,matched}/p_{T,reco}> #pm #sigma(p_{T,matched}/p_{T,reco})");\n    // gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    // TLine *l = new TLine(0, 1, 50, 1);\n    // l->SetLineColor(kGray);\n    // l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n        (TH2 *) qa_file_new->GetObjectChecked(\n            prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");\n    assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);\n    h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);\n\n    TGraphErrors *ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n        FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n    ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");\n    ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetTitle(\n        "Mean and sigma p_{Tmatched}/p_{Treco}");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco =\n          (TH2 *) qa_file_ref->GetObjectChecked(\n              prefix + "pTRecoTruthMatchedRatio_pTReco", "TH2");\n      assert(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->SetDirectory(nullptr);\n      h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Rebin2D(16, 1);\n\n      h_ratio_ref =\n          FitProfile(h_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco);\n      ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_pTRecoTruthMatchedRatio_pTReco,\n                  h_ratio_ref, true);\n  }\n\n    c1->Draw();\n}\n')


# # Tracker hit checks
# 
# Hits per tracker and layer

# In[21]:


get_ipython().run_line_magic('jsroot', 'on')


# In[22]:


get_ipython().run_cell_magic('cpp', '', '{\n    const char *hist_name_prefix = "QAG4SimulationTracking";\n      TString prefix = TString("h_") + hist_name_prefix + TString("_");\n\n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n    \n    \n\n  if (qa_file_new)\n  {\n    TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n  if (qa_file_ref)\n  {\n    TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n\n  //MVTX, INTT, TPC\n  vector<TString> detectors{"MVTX", "INTT", "TPC"};\n  vector<int> eff_ncluster_cuts{2, 2, 40};\n  vector<double> ncluster_spectrum_pt_cuts{2, 2, 2};\n  vector<TH2 *> h_pass_detectors(3, nullptr);\n  static const int nrebin = 5;\n\n  h_pass_detectors[0] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nMVTX_nReco_pTGen", "TH1") ;\n  h_pass_detectors[1] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nINTT_nReco_pTGen", "TH1") ;\n  h_pass_detectors[2] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nTPC_nReco_pTGen", "TH1") ;\n\n  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n      prefix + "nGen_pTGen", "TH1") ;\n  assert(h_norm);\n    h_norm->SetDirectory(nullptr);\n  h_norm->Rebin(nrebin);\n\n  vector<TH2 *> h_pass_detectors_ref(3, nullptr);\n  TH1 *h_norm_ref = nullptr;\n  if (qa_file_ref)\n  {\n    h_pass_detectors_ref[0] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nMVTX_nReco_pTGen", "TH1") ;\n    h_pass_detectors_ref[1] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nINTT_nReco_pTGen", "TH1") ;\n    h_pass_detectors_ref[2] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nTPC_nReco_pTGen", "TH1") ;\n\n    h_norm_ref = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + "nGen_pTGen", "TH1") ;\n    h_norm_ref->SetDirectory(nullptr);\n    h_norm_ref->Rebin(nrebin);\n\n  }\n\n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(3, 2);\n  TPad *p;\n\n  for (int i = 0; i < 3; ++i)\n  {\n    TString detector = detectors[i];\n    TH2 *h_pass_detector = h_pass_detectors[i];\n    TH2 *h_pass_detector_ref = h_pass_detectors_ref[i];\n    assert(h_pass_detector);\n\n    {\n      p = (TPad *) c1->cd(i + 1);\n      c1->Update();\n      p->SetLogy();\n\n      const int bin_start = h_pass_detector->GetXaxis()->FindBin(ncluster_spectrum_pt_cuts[i]);\n\n      TH1 *h_pass_detector_ncluster = h_pass_detector->ProjectionY(\n          TString(h_pass_detector->GetName()) + "_nCluster_new",\n          bin_start);\n      TH1 *h_pass_detector_ncluster_ref = nullptr;\n      if (h_pass_detector_ref)\n      {\n        h_pass_detector_ncluster_ref = h_pass_detector_ref->ProjectionY(\n            TString(h_pass_detector_ref->GetName()) + "_nCluster_ref",\n            bin_start);\n      }\n\n      h_pass_detector_ncluster->SetTitle(TString(hist_name_prefix) + ": " + detector + Form(" n_{Cluster} | p_{T} #geq %.1fGeV/c", ncluster_spectrum_pt_cuts[i]));\n      h_pass_detector_ncluster->SetYTitle("# of reconstructed track");\n      DrawReference(h_pass_detector_ncluster, h_pass_detector_ncluster_ref, false);\n    }\n\n    {\n      p = (TPad *) c1->cd(i + 3 + 1);\n      c1->Update();\n      p->SetLogx();\n      p->SetGridy();\n\n      const int bin_start = h_pass_detector->GetYaxis()->FindBin(eff_ncluster_cuts[i]);\n      TH1 *h_pass = h_pass_detector->ProjectionX(\n          TString(h_pass_detector->GetName()) + "_nReco_new",\n          bin_start);\n\n      assert(h_pass);\n        h_pass->SetDirectory(nullptr);\n      h_pass->Rebin(nrebin);\n\n      TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n      h_ratio->GetYaxis()->SetTitle("Reco efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));\n      h_ratio->GetYaxis()->SetRangeUser(-0, 1.);\n      //\n      TH1 *h_ratio_ref = NULL;\n      if (h_pass_detector_ref)\n      {\n        TH1 *h_pass = h_pass_detector_ref->ProjectionX(\n            TString(h_pass_detector->GetName()) + "_nReco_ref",\n            bin_start);\n\n        assert(h_pass);\n        h_pass->SetDirectory(nullptr);\n      h_pass->Rebin(nrebin);\n\n        h_ratio_ref = GetBinominalRatio(h_pass, h_norm_ref);\n      }\n      //\n      h_ratio->SetTitle("Tracking efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));\n      DrawReference(h_ratio, h_ratio_ref, false);\n    }\n  }\n\n  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# In[23]:


get_ipython().run_cell_magic('cpp', '', '{\n    const char *hist_name_prefix = "QAG4SimulationTracking";\n      TString prefix = TString("h_") + hist_name_prefix + TString("_");\n\n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n    \n    \n\n  if (qa_file_new)\n  {\n    TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n  if (qa_file_ref)\n  {\n    TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("Normalization"), "TH1");\n    assert(h_norm);\n\n    Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n\n  //MVTX, INTT, TPC\n  vector<TString> detectors{"MVTX", "INTT", "TPC"};\n  vector<int> eff_ncluster_cuts{2, 2, 40};\n  vector<double> ncluster_spectrum_pt_cuts{2, 2, 2};\n  vector<TH2 *> h_pass_detectors(3, nullptr);\n  static const int nrebin = 5;\n\n  h_pass_detectors[0] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nMVTX_nReco_pTGen", "TH1") ;\n  h_pass_detectors[1] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nINTT_nReco_pTGen", "TH1") ;\n  h_pass_detectors[2] = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + "nTPC_nReco_pTGen", "TH1") ;\n\n  TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(\n      prefix + "nGen_pTGen", "TH1") ;\n  assert(h_norm);\n    h_norm->SetDirectory(nullptr);\n  h_norm->Rebin(nrebin);\n\n  vector<TH2 *> h_pass_detectors_ref(3, nullptr);\n  TH1 *h_norm_ref = nullptr;\n  if (qa_file_ref)\n  {\n    h_pass_detectors_ref[0] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nMVTX_nReco_pTGen", "TH1") ;\n    h_pass_detectors_ref[1] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nINTT_nReco_pTGen", "TH1") ;\n    h_pass_detectors_ref[2] = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + "nTPC_nReco_pTGen", "TH1") ;\n\n    h_norm_ref = (TH1 *) qa_file_ref->GetObjectChecked(\n        prefix + "nGen_pTGen", "TH1") ;\n    h_norm_ref->SetDirectory(nullptr);\n    h_norm_ref->Rebin(nrebin);\n\n  }\n\n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_TruthMatching_NumOfClusters") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(3, 2);\n  TPad *p;\n\n  for (int i = 0; i < 3; ++i)\n  {\n    TString detector = detectors[i];\n    TH2 *h_pass_detector = h_pass_detectors[i];\n    TH2 *h_pass_detector_ref = h_pass_detectors_ref[i];\n    assert(h_pass_detector);\n\n    {\n      p = (TPad *) c1->cd(i + 1);\n      c1->Update();\n      p->SetLogy();\n\n      const int bin_start = h_pass_detector->GetXaxis()->FindBin(ncluster_spectrum_pt_cuts[i]);\n\n      TH1 *h_pass_detector_ncluster = h_pass_detector->ProjectionY(\n          TString(h_pass_detector->GetName()) + "_nCluster_new",\n          bin_start);\n      TH1 *h_pass_detector_ncluster_ref = nullptr;\n      if (h_pass_detector_ref)\n      {\n        h_pass_detector_ncluster_ref = h_pass_detector_ref->ProjectionY(\n            TString(h_pass_detector_ref->GetName()) + "_nCluster_ref",\n            bin_start);\n      }\n\n      h_pass_detector_ncluster->SetTitle(TString(hist_name_prefix) + ": " + detector + Form(" n_{Cluster} | p_{T} #geq %.1fGeV/c", ncluster_spectrum_pt_cuts[i]));\n      h_pass_detector_ncluster->SetYTitle("# of reconstructed track");\n      DrawReference(h_pass_detector_ncluster, h_pass_detector_ncluster_ref, false);\n    }\n\n    {\n      p = (TPad *) c1->cd(i + 3 + 1);\n      c1->Update();\n      p->SetLogx();\n      p->SetGridy();\n\n      const int bin_start = h_pass_detector->GetYaxis()->FindBin(eff_ncluster_cuts[i]);\n      TH1 *h_pass = h_pass_detector->ProjectionX(\n          TString(h_pass_detector->GetName()) + "_nReco_new",\n          bin_start);\n\n      assert(h_pass);\n        h_pass->SetDirectory(nullptr);\n      h_pass->Rebin(nrebin);\n\n      TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n      h_ratio->GetYaxis()->SetTitle("Reco efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));\n      h_ratio->GetYaxis()->SetRangeUser(-0, 1.);\n      //\n      TH1 *h_ratio_ref = NULL;\n      if (h_pass_detector_ref)\n      {\n        TH1 *h_pass = h_pass_detector_ref->ProjectionX(\n            TString(h_pass_detector->GetName()) + "_nReco_ref",\n            bin_start);\n\n        assert(h_pass);\n        h_pass->SetDirectory(nullptr);\n      h_pass->Rebin(nrebin);\n\n        h_ratio_ref = GetBinominalRatio(h_pass, h_norm_ref);\n      }\n      //\n      h_ratio->SetTitle("Tracking efficiency | " + detector + Form(" n_{Cluster} #geq %d", eff_ncluster_cuts[i]));\n      DrawReference(h_ratio, h_ratio_ref, false);\n    }\n  }\n\n  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# ## Hits per layer

# In[24]:


get_ipython().run_cell_magic('cpp', '-d', '\n\nnamespace\n{\n  // Normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  void GetNormalization(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)\n  {\n    if (qa_file_new)\n    {\n      TH1 *h_norm = (TH1 *) qa_file_new->GetObjectChecked(prefix + TString("Normalization"), "TH1");\n      assert(h_norm);\n      Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));\n    }\n\n    if (qa_file_ref)\n    {\n      TH1 *h_norm = (TH1 *) qa_file_ref->GetObjectChecked(prefix + TString("Normalization"), "TH1");\n      assert(h_norm);\n      Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin(tag));\n    }\n  }\n\n  void Draw(TFile *qa_file_new, TFile *qa_file_ref, const TString &prefix, const TString &tag)\n  {\n    auto h_new = static_cast<TH1 *>(qa_file_new->GetObjectChecked(prefix + tag, "TH1"));\n    assert(h_new);\n    //h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1 *h_ref = nullptr;\n    if (qa_file_ref)\n    {\n      h_ref = static_cast<TH1 *>(qa_file_ref->GetObjectChecked(prefix + tag, "TH1"));\n      assert(h_ref);\n      //h_ref->Sumw2();\n      h_ref->Scale(1.0 / Nevent_ref);\n    }\n\n    DrawReference(h_new, h_ref);\n    HorizontalLine(gPad, 1)->Draw();\n  }\n\n}  // namespace\n')


# In[25]:


get_ipython().run_cell_magic('cpp', '', '{\n     const char *hist_name_prefix = "QAG4SimulationTracking";\n      TString prefix = TString("h_") + hist_name_prefix + TString("_");\n\n  auto c1 = new TCanvas(TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,\n                        TString("QA_Draw_Tracking_nClus_Layer") + TString("_") + hist_name_prefix,\n                        950, 600);\n\n  c1->Divide(2, 1);\n  c1->cd(1);\n  GetNormalization(qa_file_new, qa_file_ref, prefix, "Truth Track");\n  Draw(qa_file_new, qa_file_ref, prefix, "nClus_layerGen");\n\n  c1->cd(2);\n  GetNormalization(qa_file_new, qa_file_ref, prefix, "Reco Track");\n  Draw(qa_file_new, qa_file_ref, prefix, "nClus_layer");\n\n  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# 

# # Upsilon reconstruction
# 
# One $\Upsilon(1S) \rightarrow e^+ e^-$ is embedded. Here is its reco results

# In[26]:


get_ipython().run_cell_magic('cpp', '-d', '\n// square\ninline double square( const double& x ) { return x*x; }\n\n// christal ball function for momentum resolution fit\nDouble_t CBcalc(Double_t *xx, Double_t *par)\n{\n  // Crystal Ball fit to one state\n  double f;\n  const double x = xx[0];\n\n  // The four parameters (alpha, n, x_mean, sigma) plus normalization (N) are:\n   \n  const double alpha = par[0];\n  const double n = par[1];\n  const double x_mean = par[2];\n  const double sigma = par[3];\n  const double N = par[4];\n\n  // dimensionless variable\n  const double t = (x-x_mean)/sigma;\n\n  // The Crystal Ball function is:\n   \n  if( t > -alpha)\n  {\n    return N * exp( -square(t)/2 );\n  }\n  else\n  {\n    const double A = pow( (n/TMath::Abs(alpha)),n) * exp(-square(alpha)/2.0);\n    const double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);\n    return N * A * pow(B - t, -n);\n  }\n\n}\n\n// christal ball function for Upsilon fits\nDouble_t CBcalc2(Double_t *xx, Double_t *par)\n{\n  // Crystal Ball fit to one state\n  const double x = xx[0];\n\n  /* The parameters are: \n   * N the normalization\n   * x_mean, sigma\n   * alpha_left, n_left\n   * alpha_right, n_right\n   */\n  const double N = par[0];\n  const double x_mean = par[1];\n  const double sigma = par[2];\n  const double t = (x-x_mean)/sigma;    \n\n  // left tail\n  const double alpha_left = std::abs(par[3]);\n  const double n_left = par[4];\n\n  // right tail\n  const double alpha_right = std::abs(par[5]);\n  const double n_right = par[6];\n    \n  // The Crystal Ball function is:   \n  if( t < -alpha_left )\n  {      \n    const double A = pow( (n_left/TMath::Abs(alpha_left)),n_left) * exp(-square(alpha_left)/2.0);\n    const double B = n_left/std::abs(alpha_left) - std::abs(alpha_left);\n    return N * A * pow(B - t, -n_left);\n  } else if( t > alpha_right ) {\n     const double A = pow( (n_right/TMath::Abs(alpha_right)),n_right) * exp(-square(alpha_right)/2.0);\n     const double B = n_right/std::abs(alpha_right) - std::abs(alpha_right);\n     return N * A * pow(B + t, -n_right);  \n  } else {\n      return N * exp( -square(t)/2);\n  }\n}\n')


# In[27]:


get_ipython().run_cell_magic('cpp', '', '\n{\n  const char *hist_name_prefix = "QAG4SimulationUpsilon";\n  TString prefix = TString("h_") + hist_name_prefix + TString("_");\n\n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if ( qa_file_new->GetObjectChecked(\n        prefix + TString("pTRecoGenRatio_pTGen"), "TH2")\n     == nullptr )\n    {\n        cout <<"QAG4SimulationUpsilon is not enabled. Skip...."<<endl;\n    }\n    else\n    {\n\n      TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,\n                                TString("QA_Draw_Tracking_UpsilonOverview") + TString("_") + hist_name_prefix,\n                                950, 600);\n      c1->Divide(2, 1);\n      int idx = 1;\n      TPad *p;\n\n      {\n        p = (TPad *) c1->cd(idx++);\n        c1->Update();\n        p->SetLogy();\n\n        TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n            prefix + TString("pTRecoGenRatio_pTGen"), "TH2");\n        assert(h_new);\n\n        //  h_new->Rebin(1, 2);\n        //h_new->Sumw2();\n        //  h_new->Scale(1. / Nevent_new);\n\n        TH2 *h_ref = NULL;\n        if (qa_file_ref)\n        {\n          h_ref = (TH2 *) qa_file_ref->GetObjectChecked(\n              prefix + TString("pTRecoGenRatio_pTGen"), "TH2");\n          assert(h_ref);\n\n          //    h_ref->Rebin(1, 2);\n          //h_ref->Sumw2();\n          h_ref->Scale(Nevent_new / Nevent_ref);\n        }\n\n        TH1 *h_proj_new = h_new->ProjectionY(\n            TString::Format(\n                "%s_New_ProjX",\n                h_new->GetName()));\n\n        h_proj_new->GetXaxis()->SetRangeUser(0, 1.3);\n        h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                             ": Electron lineshape"));\n        h_proj_new->GetXaxis()->SetTitle(TString::Format(\n            "Reco p_{T}/Truth p_{T}"));\n\n        TF1 *f_eLineshape = new TF1("f_eLineshape", CBcalc, 7, 11, 5);\n        f_eLineshape->SetParameter(0, 1.0);\n        f_eLineshape->SetParameter(1, 1.0);\n        f_eLineshape->SetParameter(2, 0.95);\n        f_eLineshape->SetParameter(3, 0.08);\n        f_eLineshape->SetParameter(4, 20.0);\n\n        f_eLineshape->SetParNames("alpha","n","m","sigma","N");\n        f_eLineshape->SetLineColor(kRed);\n        f_eLineshape->SetLineWidth(3);\n        f_eLineshape->SetLineStyle(kSolid);\n        f_eLineshape->SetNpx(1000);\n\n        h_proj_new->Fit(f_eLineshape);\n\n        TH1 *h_proj_ref = nullptr;\n        if (h_ref)\n        {\n          h_proj_ref =\n              h_ref->ProjectionY(\n                  TString::Format(\n                      "%s_Ref_ProjX",\n                      h_new->GetName()));\n        }\n        TF1 *f_eLineshape_ref = new TF1("f_eLineshape_ref", CBcalc, 7, 11, 5);\n        f_eLineshape_ref->SetParameter(0, 1.0);\n        f_eLineshape_ref->SetParameter(1, 1.0);\n        f_eLineshape_ref->SetParameter(2, 0.95);\n        f_eLineshape_ref->SetParameter(3, 0.08);\n        f_eLineshape_ref->SetParameter(4, 20.0);\n\n        f_eLineshape_ref->SetParNames("alpha","n","m","sigma","N");\n        f_eLineshape_ref->SetLineColor(kRed);\n        f_eLineshape_ref->SetLineWidth(3);\n        f_eLineshape_ref->SetLineStyle(kSolid);\n\n        h_proj_ref->Fit(f_eLineshape_ref);\n\n\n        DrawReference(h_proj_new, h_proj_ref);\n        f_eLineshape->Draw("same");\n\n        char resstr_1[500];\n        sprintf(resstr_1,"#sigma_{dp/p} = %.2f #pm %.2f %%", f_eLineshape->GetParameter(3)*100, f_eLineshape->GetParError(3)*100);\n        TLatex *res_1 = new TLatex(0.2,0.75,resstr_1);\n        res_1->SetNDC();\n        res_1->SetTextSize(0.05);\n        res_1->SetTextAlign(13);\n        res_1->Draw();\n\n        char resstr_2[500];\n        sprintf(resstr_2,"#sigma_{dp/p,ref} = %.2f #pm %.2f %%", f_eLineshape_ref->GetParameter(3)*100, f_eLineshape_ref->GetParError(3)*100);\n        TLatex *res_2 = new TLatex(0.2,0.7,resstr_2);\n        res_2->SetNDC();\n        res_2->SetTextSize(0.05);\n        res_2->SetTextAlign(13);\n        res_2->Draw();\n      }\n\n      {\n        p = (TPad *) c1->cd(idx++);\n        c1->Update();\n    //    p->SetLogy();\n\n        TH1 *h_new = (TH1 *) qa_file_new->GetObjectChecked(\n            prefix + TString("nReco_Pair_InvMassReco"), "TH1");\n        assert(h_new);\n\n        // h_new->Rebin(2);\n        // h_new->Sumw2();\n        // h_new->Scale(1. / Nevent_new);\n\n        TF1 *f1S = new TF1("f1S", CBcalc2, 7, 11, 7);\n        f1S->SetParameter(0, 50.0);\n        f1S->SetParameter(1, 9.46);\n        f1S->SetParameter(2, 0.08);\n        f1S->SetParameter(3, 1.0);\n        f1S->SetParameter(4, 3.0);\n        f1S->SetParameter(5, 1.0);\n        f1S->SetParLimits(3, 0.120, 10);\n        f1S->SetParLimits(4, 1.05, 10);\n        f1S->SetParameter(6, 3.0);\n        f1S->SetParLimits(5, 0.120, 10);\n        f1S->SetParLimits(6, 1.05, 10);\n \n        f1S->SetParNames("N", "m", "#sigma", "#alpha_{left}","n_{left}","#alpha_{right}","#sigma_{right}");\n        f1S->SetLineColor(kRed);\n        f1S->SetLineWidth(3);\n        f1S->SetLineStyle(kSolid);\n        f1S->SetNpx(1000);\n\n        h_new->Fit(f1S);\n\n        TH1 *h_ref = NULL;\n        if (qa_file_ref)\n        {\n          h_ref = (TH1 *) qa_file_ref->GetObjectChecked(\n              prefix + TString("nReco_Pair_InvMassReco"), "TH1");\n          assert(h_ref);\n\n          // h_ref->Rebin(2);\n          // h_ref->Sumw2();\n          // h_ref->Scale(Nevent_new / Nevent_ref);\n        }\n\n        h_new->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                        ": #Upsilon #rightarrow e^{+}e^{-} lineshape"));\n        h_new->GetXaxis()->SetRangeUser(7, 10);\n\n        TF1 *f1S_ref = new TF1("f1S_ref", CBcalc2, 7, 11, 7);\n        f1S_ref->SetParameter(0, 50.0);\n        f1S_ref->SetParameter(1, 9.46);\n        f1S_ref->SetParameter(2, 0.08);\n        f1S_ref->SetParameter(3, 1.0);\n        f1S_ref->SetParameter(4, 3.0);\n        f1S_ref->SetParameter(5, 1.0);\n        f1S_ref->SetParLimits(3, 0.120, 10);\n        f1S_ref->SetParLimits(4, 1.05, 10);\n        f1S_ref->SetParameter(6, 3.0);\n        f1S_ref->SetParLimits(5, 0.120, 10);\n        f1S_ref->SetParLimits(6, 1.05, 10);\n \n        f1S_ref->SetParNames("N", "m", "#sigma", "#alpha_{left}","n_{left}","#alpha_{right}","#sigma_{right}");\n        f1S_ref->SetLineColor(kRed);\n        f1S_ref->SetLineWidth(3);\n        f1S_ref->SetLineStyle(kSolid);\n\n        h_ref->Fit(f1S_ref);\n\n        DrawReference(h_new, h_ref, false);\n        f1S->Draw("same");\n\n        // cout << "f1S pars " <<  f1S->GetParameter(3) << "   " << f1S->GetParError(3) << endl;\n\n        char resstr_3[500];\n        sprintf(resstr_3,"#sigma_{1S} = %.1f #pm %.1f MeV", f1S->GetParameter(2)*1000, f1S->GetParError(2)*1000);\n        TLatex *res_3 = new TLatex(0.2,0.75,resstr_3);\n        res_3->SetNDC();\n        res_3->SetTextSize(0.05);\n        res_3->SetTextAlign(13);\n        res_3->Draw();\n\n        char resstr_4[500];\n        sprintf(resstr_4,"#sigma_{1S,ref} = %.1f #pm %.1f MeV", f1S_ref->GetParameter(2)*1000, f1S_ref->GetParError(2)*1000);\n        TLatex *res_4 = new TLatex(0.2,0.7,resstr_4);\n        res_4->SetNDC();\n        res_4->SetTextSize(0.05);\n        res_4->SetTextAlign(13);\n        res_4->Draw();\n          \n        ofstream fcsv;\n        fcsv.open ("Upsilon_mean.csv");\n        fcsv<<"Upsilon mean (GeV)"<<endl;\n        fcsv<<f1S->GetParameter(1)<<endl;\n        fcsv.close();\n        fcsv.open ("Upsilon_width.csv");\n        fcsv<<"Upsilon width (GeV)"<<endl;\n        fcsv<<f1S->GetParameter(2)<<endl;\n        fcsv.close();\n        fcsv.open ("Upsilon_count.csv");\n        fcsv<<"Upsilon reconstructed"<<endl;\n        fcsv<<f1S->GetParameter(0)<<endl;\n        fcsv.close();\n      }\n\n      // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n\n      c1 -> Draw();\n    }// if checks\n}\n')


# # Summary statistics

# In[28]:


get_ipython().run_cell_magic('cpp', '', '\nKSTestSummary::getInstance()->make_summary_txt("QA-tracking.txt");\n')


# In[29]:


get_ipython().run_cell_magic('cpp', '', '\nKSTestSummary::getInstance()->make_summary_TCanvas() -> Draw();\n')


# In[ ]:





# In[ ]: