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

# # Introduction
# 
# QA plots for the generic vetex finding and DCA 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"
    qa_file_name_new = "G4sPHENIX_test-tracking_Event1000_Sum16_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 = "G4sPHENIX_test-tracking_Event1000_Sum16_qa.root"
#     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"Openning 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"Openning 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')


# # DCA performance
# 
# Per-track DCA performance

# ## Transverse DCA, $DCA_{r\phi}$
# 
# First without track selection cuts, which included many track that do not have good MVTX matches

# 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    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n    \n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("DCArPhi_pT"), "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("DCArPhi_pT"), "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_DCArPhi") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCArPhi") + 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, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\n\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    if (pt_range.first < 2.0)\n    {\n      h_proj_new->GetXaxis()->SetRangeUser(-.05,.05);\n      h_proj_new->Rebin(5);\n    }\n    else\n    {\n      h_proj_new->GetXaxis()->SetRangeUser(-.01,.01);\n    }\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        "DCA (r #phi) [cm]"));\n    h_proj_new->GetXaxis()->SetNdivisions(5,5);\n    f1 = new TF1("f1","gaus",-.01,.01);\n    h_proj_new->Fit(f1,"qm");\n    sigma = f1->GetParameter(2);\n    sigma_unc = f1->GetParError(2);\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_%d_%d",\n                  h_new->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n      if (pt_range.first < 2.0)\n      {\n\t//h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n\th_proj_ref->Rebin(5);\n      }\n    }\n    \n    DrawReference(h_proj_new, h_proj_ref);\n\n    sprintf(resstr,"#sigma = %.5f #pm %.5f cm", 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  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  TPaveText *pt = new TPaveText(.05,.1,.95,.8);\n  pt->AddText("No cuts");\n  pt->Draw();\n\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n\n\n\n    c1->Draw();\n}\n')


# Next, with cuts on minimal tracker hits as defined in `QAG4SimulationTracking` and denoted on plot

# In[14]:


get_ipython().run_cell_magic('cpp', '', '{\n    const 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  TH2 *h_new2 = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("DCArPhi_pT_cuts"), "TH2");\n  assert(h_new2);\n\n  //  h_new->Rebin(1, 2);\n  //h_new2->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref2 = NULL;\n  if (qa_file_ref)\n  {\n    h_ref2 = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("DCArPhi_pT_cuts"), "TH2");\n    assert(h_ref2);\n\n    //    h_ref->Rebin(1, 2);\n    //h_ref2->Sumw2();\n    h_ref2->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c2 = new TCanvas(TString("QA_Draw_Tracking_DCArPhi2") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCArPhi2") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c2->Divide(4, 2);\n  int idx2 = 1;\n  TPad *p2;\n\n  vector<pair<double, double>> gpt_ranges2{\n      {0, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\n  TF1 *f2 = nullptr;\n  TF1 *fit2 = nullptr;\n  Double_t sigma2 = 0;\n  Double_t sigma_unc2 = 0;\n  char resstr2[500];\n  TLatex *res2 = nullptr;\n  for (auto pt_range : gpt_ranges2)\n  {\n    //cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";\n\n    p2 = (TPad *) c2->cd(idx2++);\n    c2->Update();\n    p2->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new2->GetXaxis()->FindBin(pt_range.first + epsilon);\n    const int bin_end = h_new2->GetXaxis()->FindBin(pt_range.second - epsilon);\n\n    TH1 *h_proj_new2 = h_new2->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new2->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n    if (pt_range.first < 2.0)\n    {\n      h_proj_new2->GetXaxis()->SetRangeUser(-.05,.05);\n      h_proj_new2->Rebin(5);\n    }\n    else\n    {\n      h_proj_new2->GetXaxis()->SetRangeUser(-.01,.01);\n    }\n    h_proj_new2->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));\n    h_proj_new2->GetXaxis()->SetTitle(TString::Format(\n        "DCA (r #phi) [cm]"));\n    h_proj_new2->GetXaxis()->SetNdivisions(5,5);\n    f2 = new TF1("f2","gaus",-.01,.01);\n    h_proj_new2->Fit(f2 , "mq");\n    sigma2 = f2->GetParameter(2);\n    sigma_unc2 = f2->GetParError(2);\n\n    TH1 *h_proj_ref2 = nullptr;\n    if (h_ref2)\n    {\n      h_proj_ref2 =\n          h_ref2->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new2->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n      if (pt_range.first < 2.0)\n      {\n\t//h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n\th_proj_ref2->Rebin(5);\n      }\n    }\n    DrawReference(h_proj_new2, h_proj_ref2);\n\n    sprintf(resstr2,"#sigma = %.5f #pm %.5f cm", sigma2, sigma_unc2);\n    res2 = new TLatex(0.325,0.825,resstr2);\n    res2->SetNDC();\n    res2->SetTextSize(0.05);\n    res2->SetTextAlign(13);\n    res2->Draw();\n  }\n    \n    \n  p2 = (TPad *) c2->cd(idx2++);\n  c2->Update();\n  TPaveText *pt2 = new TPaveText(.05,.1,.95,.8);\n  pt2->AddText("Cuts: MVTX hits>=2, INTT hits>=1,");\n  pt2->AddText("TPC hits>=20");\n  pt2->Draw();\n\n  //SaveCanvas(c2, TString(qa_file_name_new) + TString("_") + TString(c2->GetName()), true);\n\n    c2->Draw();\n}\n')


# ## Longitudinal DCA, $DCA_z$
# 
# Start with the version without track quality cuts

# In[15]:


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    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\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("Truth Track"));\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("Truth Track"));\n  }\n\n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("DCAZ_pT"), "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("DCAZ_pT"), "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_DCAZ") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCAZ") + 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, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\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    if (pt_range.first < 2.0)\n    {\n      h_proj_new->GetXaxis()->SetRangeUser(-.05, .05);\n      h_proj_new->Rebin(5);\n    }\n    else\n    {\n      h_proj_new->GetXaxis()->SetRangeUser(-.01, .01);\n    }\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        "DCA (Z) [cm]"));\n    h_proj_new->GetXaxis()->SetNdivisions(5, 5);\n\n    f1 = new TF1("f1", "gaus", -.01, .01);\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    {\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      if (pt_range.first < 2.0)\n      {\n        //h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n        h_proj_ref->Rebin(5);\n      }\n    }\n    DrawReference(h_proj_new, h_proj_ref);\n    sprintf(resstr, "#sigma = %.5f #pm %.5f cm", 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  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  TPaveText *pt = new TPaveText(.05, .1, .95, .8);\n  pt->AddText("No cuts");\n  pt->Draw();\n\n//  SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    \n    c1->Draw();\n}\n')


# with cuts on minimal hits on each of the trackers

# In[16]:


get_ipython().run_cell_magic('cpp', '', '\n{\n    \n    const 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  // cuts plots\n  TH2 *h_new2 = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("DCAZ_pT_cuts"), "TH2");\n  assert(h_new2);\n\n  //  h_new->Rebin(1, 2);\n  //h_new2->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref2 = NULL;\n  if (qa_file_ref)\n  {\n    h_ref2 = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("DCAZ_pT_cuts"), "TH2");\n    assert(h_ref2);\n\n    //    h_ref->Rebin(1, 2);\n    //h_ref2->Sumw2();\n    h_ref2->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c2 = new TCanvas(TString("QA_Draw_Tracking_DCAZ2") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCAZ2") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c2->Divide(4, 2);\n  int idx2 = 1;\n  TPad *p2;\n\n  vector<pair<double, double>> gpt_ranges2{\n      {0, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\n  TF1 *f2 = nullptr;\n  TF1 *fit2 = nullptr;\n  Double_t sigma2 = 0;\n  Double_t sigma_unc2 = 0;\n  char resstr2[500];\n  TLatex *res2 = nullptr;\n  for (auto pt_range : gpt_ranges2)\n  {\n    //cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";\n\n    p2 = (TPad *) c2->cd(idx2++);\n    c2->Update();\n    p2->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new2->GetXaxis()->FindBin(pt_range.first + epsilon);\n    const int bin_end = h_new2->GetXaxis()->FindBin(pt_range.second - epsilon);\n\n    TH1 *h_proj_new2 = h_new2->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new2->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n    if (pt_range.first < 2.0)\n    {\n      h_proj_new2->GetXaxis()->SetRangeUser(-.05, .05);\n      h_proj_new2->Rebin(5);\n    }\n    else\n    {\n      h_proj_new2->GetXaxis()->SetRangeUser(-.01, .01);\n    }\n    h_proj_new2->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                          ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));\n    h_proj_new2->GetXaxis()->SetTitle(TString::Format(\n        "DCA (Z) [cm]"));\n    h_proj_new2->GetXaxis()->SetNdivisions(5, 5);\n\n    f2 = new TF1("f2", "gaus", -.01, .01);\n    h_proj_new2->Fit(f2,"mq");\n    sigma2 = f2->GetParameter(2);\n    sigma_unc2 = f2->GetParError(2);\n\n    TH1 *h_proj_ref2 = nullptr;\n    if (h_ref2)\n    {\n      h_proj_ref2 =\n          h_ref2->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new2->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n      if (pt_range.first < 2.0)\n      {\n        //h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n        h_proj_ref2->Rebin(5);\n      }\n    }\n    DrawReference(h_proj_new2, h_proj_ref2);\n    sprintf(resstr2, "#sigma = %.5f #pm %.5f cm", sigma2, sigma_unc2);\n    res2 = new TLatex(0.325, 0.825, resstr2);\n    res2->SetNDC();\n    res2->SetTextSize(0.05);\n    res2->SetTextAlign(13);\n    res2->Draw();\n  }\n  p2 = (TPad *) c2->cd(idx2++);\n  c2->Update();\n  TPaveText *pt2 = new TPaveText(.05, .1, .95, .8);\n  pt2->AddText("Cuts: MVTX hits>=2, INTT hits>=1,");\n  pt2->AddText("TPC hits>=20");\n  pt2->Draw();\n\n//  SaveCanvas(c2, TString(qa_file_name_new) + TString("_") + TString(c2->GetName()), true);\n    c2->Draw();\n}\n')


# In[17]:


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


# In[18]:


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    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\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  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_DCA_Resolution") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCA_Resolution") + 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->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -0.01, 50, 0.01,\n                              ";Truth p_{T} [GeV/c];<DCA (r #phi)> #pm #sigma(DCA (r #phi)) [cm]");\n    gPad->SetLeftMargin(.2);\n    frame->GetYaxis()->SetTitleOffset(2);\n    TLine *l = new TLine(0.1, 0, 50, 0);\n    l->SetLineColor(kGray);\n    l->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "DCArPhi_pT_cuts", "TH2");\n    assert(h_QAG4SimulationTracking_DCArPhi);\n\n    h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n    // h_QAG4SimulationTracking_DCArPhi->Draw("colz");\n    TGraphErrors *ge_QAG4SimulationTracking_DCArPhi = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n    ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "DCArPhi_pT_cuts", "TH2");\n      assert(h_QAG4SimulationTracking_DCArPhi);\n\n      h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n      ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n    }\n\n    ge_QAG4SimulationTracking_DCArPhi->SetTitle("DCA (r #phi, #geq 2MVTX, #geq 1INTT, #geq 20TPC) [cm]");\n    DrawReference(ge_QAG4SimulationTracking_DCArPhi, h_ratio_ref, true);\n      \n  }\n\n  {\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -0.01, 50, 0.01,\n                              "DCA (Z) [cm];Truth p_{T} [GeV/c];<DCA (Z)> #pm #sigma(DCA (Z)) [cm]");\n    // gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    //TLine *l = new TLine(0.1, 0, 50, 0);\n    //l->SetLineColor(kGray);\n    //l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "DCAZ_pT_cuts", "TH2");\n    assert(h_QAG4SimulationTracking_DCAZ);\n\n    h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n    TGraphErrors *ge_QAG4SimulationTracking_DCAZ = FitProfile(h_QAG4SimulationTracking_DCAZ);\n    ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    ge_QAG4SimulationTracking_DCAZ->SetTitle("DCA (Z) [cm]");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "DCAZ_pT_cuts", "TH2");\n      assert(h_QAG4SimulationTracking_DCAZ);\n\n      h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCAZ);\n      ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_DCAZ, h_ratio_ref, true);\n  }\n\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# In[19]:


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    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\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  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_DCA_Resolution") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_DCA_Resolution") + 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->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -0.01, 50, 0.01,\n                              ";Truth p_{T} [GeV/c];<DCA (r #phi)> #pm #sigma(DCA (r #phi)) [cm]");\n    gPad->SetLeftMargin(.2);\n    frame->GetYaxis()->SetTitleOffset(2);\n    TLine *l = new TLine(0.1, 0, 50, 0);\n    l->SetLineColor(kGray);\n    l->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "DCArPhi_pT_cuts", "TH2");\n    assert(h_QAG4SimulationTracking_DCArPhi);\n\n    h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n    // h_QAG4SimulationTracking_DCArPhi->Draw("colz");\n    TGraphErrors *ge_QAG4SimulationTracking_DCArPhi = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n    ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "DCArPhi_pT_cuts", "TH2");\n      assert(h_QAG4SimulationTracking_DCArPhi);\n\n      h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n      ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n    }\n\n    ge_QAG4SimulationTracking_DCArPhi->SetTitle("DCA (r #phi, #geq 2MVTX, #geq 1INTT, #geq 20TPC) [cm]");\n    DrawReference(ge_QAG4SimulationTracking_DCArPhi, h_ratio_ref, true);\n      \n      SaveGraphError2CSV(ge_QAG4SimulationTracking_DCArPhi, "QAG4SimulationTracking_DCArPhi");\n  }\n\n  {\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -0.01, 50, 0.01,\n                              "DCA (Z) [cm];Truth p_{T} [GeV/c];<DCA (Z)> #pm #sigma(DCA (Z)) [cm]");\n    // gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    //TLine *l = new TLine(0.1, 0, 50, 0);\n    //l->SetLineColor(kGray);\n    //l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "DCAZ_pT_cuts", "TH2");\n    assert(h_QAG4SimulationTracking_DCAZ);\n\n    h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n    TGraphErrors *ge_QAG4SimulationTracking_DCAZ = FitProfile(h_QAG4SimulationTracking_DCAZ);\n    ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    ge_QAG4SimulationTracking_DCAZ->SetTitle("DCA (Z) [cm]");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "DCAZ_pT_cuts", "TH2");\n      assert(h_QAG4SimulationTracking_DCAZ);\n\n      h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCAZ);\n      ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_DCAZ, h_ratio_ref, true);\n      SaveGraphError2CSV(ge_QAG4SimulationTracking_DCAZ, "QAG4SimulationTracking_DCAZ");\n  }\n\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# # Sigmalized DCA

# In[20]:


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


# In[21]:


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  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("Truth Track"));\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("Truth Track"));\n  }\n    \n    \n    \n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("SigmalizedDCArPhi_pT"), "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("SigmalizedDCArPhi_pT"), "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_SigmalizedDCArPhi") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_SigmalizedDCArPhi") + 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, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\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    h_proj_new->GetXaxis()->SetRangeUser(-5.,5.);\n    h_proj_new->Rebin(5);\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        "Sigmalized DCA (r #phi)"));\n    h_proj_new->GetXaxis()->SetNdivisions(5,5);\n\n    f1 = new TF1("f1","gaus",-4.,4.);\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    {\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      //h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n      h_proj_ref->Rebin(5);\n    }\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  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  TPaveText *pt = new TPaveText(.05,.1,.95,.8);\n  pt->AddText("Cuts: MVTX hits>=2, INTT hits>=1,");\n  pt->AddText("TPC hits>=20");\n  pt->Draw();\n\n  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n  c1->Draw();\n}\n')


# In[22]:


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  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("Truth Track"));\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("Truth Track"));\n  }\n    \n    \n    \n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("SigmalizedDCArPhi_pT"), "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("SigmalizedDCArPhi_pT"), "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_SigmalizedDCArPhi") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_SigmalizedDCArPhi") + 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, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\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    h_proj_new->GetXaxis()->SetRangeUser(-5.,5.);\n    h_proj_new->Rebin(5);\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        "Sigmalized DCA (r #phi)"));\n    h_proj_new->GetXaxis()->SetNdivisions(5,5);\n\n    f1 = new TF1("f1","gaus",-4.,4.);\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    {\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      //h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n      h_proj_ref->Rebin(5);\n    }\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  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  TPaveText *pt = new TPaveText(.05,.1,.95,.8);\n  pt->AddText("Cuts: MVTX hits>=2, INTT hits>=1,");\n  pt->AddText("TPC hits>=20");\n  pt->Draw();\n\n  // SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n  c1->Draw();\n}\n')


# ## Longitudinal DCA, $DCA_z/\sigma[DCA_z]$

# In[23]:


get_ipython().run_cell_magic('cpp', '', '\n{\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  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("Truth Track"));\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("Truth Track"));\n  }\n    \n  TH2 *h_new2 = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("SigmalizedDCAZ_pT"), "TH2");\n  assert(h_new2);\n\n  //  h_new->Rebin(1, 2);\n  //h_new2->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref2 = NULL;\n  if (qa_file_ref)\n  {\n    h_ref2 = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("SigmalizedDCAZ_pT"), "TH2");\n    assert(h_ref2);\n\n    //    h_ref->Rebin(1, 2);\n    //h_ref2->Sumw2();\n    h_ref2->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c2 = new TCanvas(TString("QA_Draw_Tracking_SigmalizedDCAZ") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_SigmalizedDCAZ") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c2->Divide(4, 2);\n  int idx2 = 1;\n  TPad *p2;\n\n  vector<pair<double, double>> gpt_ranges2{\n      {0, 0.5},\n      {0.5, 1},\n      {1, 1.5},\n      {1.5, 2},\n      {2, 4},\n      {4, 16},\n      {16, 40}};\n  TF1 *f2 = nullptr;\n  TF1 *fit2 = nullptr;\n  Double_t sigma2 = 0;\n  Double_t sigma_unc2 = 0;\n  char resstr2[500];\n  TLatex *res2 = nullptr;\n  for (auto pt_range : gpt_ranges2)\n  {\n   // cout << __PRETTY_FUNCTION__ << " process " << pt_range.first << " - " << pt_range.second << " GeV/c";\n\n    p2 = (TPad *) c2->cd(idx2++);\n    c2->Update();\n    p2->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new2->GetXaxis()->FindBin(pt_range.first + epsilon);\n    const int bin_end = h_new2->GetXaxis()->FindBin(pt_range.second - epsilon);\n\n    TH1 *h_proj_new2 = h_new2->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new2->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n    h_proj_new2->GetXaxis()->SetRangeUser(-5.,5.);\n    h_proj_new2->Rebin(5);\n    h_proj_new2->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f GeV/c", pt_range.first, pt_range.second));\n    h_proj_new2->GetXaxis()->SetTitle(TString::Format(\n        "Sigmalized DCA (Z)"));\n    h_proj_new2->GetXaxis()->SetNdivisions(5,5);\n    \n    f2 = new TF1("f2","gaus",-4.,4.);\n    h_proj_new2->Fit(f2, "mq");\n    sigma2 = f2->GetParameter(2);\n    sigma_unc2 = f2->GetParError(2);\n\n    TH1 *h_proj_ref2 = nullptr;\n    if (h_ref2)\n    {\n      h_proj_ref2 =\n          h_ref2->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new2->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n      //h_proj_ref->GetXaxis()->SetRangeUser(-.05,.05);\n      h_proj_ref2->Rebin(5);\n    }\n    DrawReference(h_proj_new2, h_proj_ref2);\n    sprintf(resstr2,"#sigma = %.5f #pm %.5f", sigma2, sigma_unc2);\n    res2 = new TLatex(0.325,0.825,resstr2);\n    res2->SetNDC();\n    res2->SetTextSize(0.05);\n    res2->SetTextAlign(13);\n    res2->Draw();\n  }\n  p2 = (TPad *) c2->cd(idx2++);\n  c2->Update();\n  TPaveText *pt2 = new TPaveText(.05,.1,.95,.8);\n  pt2->AddText("Cuts: MVTX hits>=2, INTT hits>=1,");\n  pt2->AddText("TPC hits>=20");\n  pt2->Draw();\n\n  //SaveCanvas(c2, TString(qa_file_name_new) + TString("_") + TString(c2->GetName()), true);\n    c2->Draw();\n}\n')


# In[ ]:





# ## Sigmalized DCA summaries
# 
# With tracker hit cuts (>=2MVTX, >=1INTT, >=20TPC)
# 

# In[24]:


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


# In[25]:


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    \n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\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  TCanvas *c1 = new TCanvas(TString("QA_Draw_Tracking_SigmalizedDCA_Resolution") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Tracking_SigmalizedDCA_Resolution") + 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->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -2, 50, 2,\n                              ";Truth p_{T} [GeV/c];<Sigmalized DCA (r #phi)> #pm #sigma(Sigmalized DCA (r #phi))");\n    gPad->SetLeftMargin(.2);\n    frame->GetYaxis()->SetTitleOffset(2);\n    TLine *l = new TLine(0.1, 0, 50, 0);\n    l->SetLineColor(kGray);\n    l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n    HorizontalLine(gPad, -1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "SigmalizedDCArPhi_pT", "TH2");\n    assert(h_QAG4SimulationTracking_DCArPhi);\n\n    h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n    // h_QAG4SimulationTracking_DCArPhi->Draw("colz");\n    TGraphErrors *ge_QAG4SimulationTracking_DCArPhi = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n    ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCArPhi = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "SigmalizedDCArPhi_pT", "TH2");\n      assert(h_QAG4SimulationTracking_DCArPhi);\n\n      h_QAG4SimulationTracking_DCArPhi->Rebin2D(20, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCArPhi);\n      ge_QAG4SimulationTracking_DCArPhi->Draw("pe");\n    }\n\n    ge_QAG4SimulationTracking_DCArPhi->SetTitle("DCA_{r#phi}/#sigma[DCA_{r#phi}]");\n    DrawReference(ge_QAG4SimulationTracking_DCArPhi, h_ratio_ref, true);\n  }\n\n  {\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogx();\n    TH1 *frame = p->DrawFrame(0.1, -2, 50, 2,\n                              "DCA_z/#sigma[DCA_z];Truth p_{T} [GeV/c];<Sigmalized DCA (Z)> #pm #sigma(Sigmalized DCA (Z))");\n    // gPad->SetLeftMargin(.2);\n    gPad->SetTopMargin(-1);\n    frame->GetYaxis()->SetTitleOffset(1.7);\n    TLine *l = new TLine(0.1, 0, 50, 0);\n    l->SetLineColor(kGray);\n    l->Draw();\n    HorizontalLine(gPad, 1)->Draw();\n    HorizontalLine(gPad, -1)->Draw();\n\n    TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_new->GetObjectChecked(\n        prefix + "SigmalizedDCAZ_pT", "TH2");\n    assert(h_QAG4SimulationTracking_DCAZ);\n\n    h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n    TGraphErrors *ge_QAG4SimulationTracking_DCAZ = FitProfile(h_QAG4SimulationTracking_DCAZ);\n    ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    ge_QAG4SimulationTracking_DCAZ->SetTitle("DCA_z/#sigma[DCA_z]");\n\n    TGraphErrors *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2 *h_QAG4SimulationTracking_DCAZ = (TH2 *) qa_file_ref->GetObjectChecked(\n          prefix + "SigmalizedDCAZ_pT", "TH2");\n      assert(h_QAG4SimulationTracking_DCAZ);\n\n      h_QAG4SimulationTracking_DCAZ->Rebin2D(40, 1);\n\n      h_ratio_ref = FitProfile(h_QAG4SimulationTracking_DCAZ);\n      ge_QAG4SimulationTracking_DCAZ->Draw("pe");\n    }\n\n    DrawReference(ge_QAG4SimulationTracking_DCAZ, h_ratio_ref, true);\n  }\n\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n}\n')


# # Vertexing

# ## Number of tracks associated
# 
# Tracks assocaition with vertex. `gntrack` are truth tracks associated with the vertex and `ntrack` are reconstructed tracks

# In[26]:


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


# In[27]:


get_ipython().run_cell_magic('cpp', '', '{\n    const char *hist_name_prefix = "QAG4SimulationVertex_SvtxVertexMap";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n    \n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Vertex_nVertex") +\n                                TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Vertex_nVertex") +\n                                TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(2, 2);\n  int idx = 1;\n  TPad *p;\n\n  {\n    static const int nrebin = 1;\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 + "gntracks", "TH1");\n    assert(h_pass);\n\n    h_pass->Rebin(nrebin);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_pass->GetYaxis()->SetTitle("Counts");\n    // h_pass->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ref = NULL;\n    if (qa_file_ref) {\n      h_ref =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "gntracks", "TH1");\n      assert(h_ref);\n\n      h_ref->Rebin(nrebin);\n    }\n\n    h_pass->SetTitle(TString(hist_name_prefix) + ": gntracks");\n\n    DrawReference(h_pass, h_ref, false);\n  }\n\n  {\n    static const int nrebin = 1;\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 + "gntracksmaps", "TH1");\n    assert(h_pass);\n\n    h_pass->Rebin(nrebin);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_pass->GetYaxis()->SetTitle("Counts");\n    // h_pass->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ref = NULL;\n    if (qa_file_ref) {\n      h_ref =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "gntracksmaps", "TH1");\n      assert(h_pass);\n\n      h_ref->Rebin(nrebin);\n    }\n\n    h_pass->SetTitle(TString(hist_name_prefix) + ": gntracksmaps");\n\n    DrawReference(h_pass, h_ref, false);\n  }\n\n  {\n    static const int nrebin = 1;\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 + "ntracks", "TH1");\n    assert(h_pass);\n\n    h_pass->Rebin(nrebin);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_pass->GetYaxis()->SetTitle("Counts");\n    // h_pass->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ref = NULL;\n    if (qa_file_ref) {\n      h_ref =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "ntracks", "TH1");\n      assert(h_pass);\n\n      h_ref->Rebin(nrebin);\n    }\n\n    h_pass->SetTitle(TString(hist_name_prefix) + ": ntracks");\n\n    DrawReference(h_pass, h_ref, false);\n  }\n\n  {\n    static const int nrebin = 1;\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 + "ntracks_cuts", "TH1");\n    assert(h_pass);\n\n    h_pass->Rebin(nrebin);\n\n    //    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_pass->GetYaxis()->SetTitle("Counts");\n    // h_pass->GetYaxis()->SetRangeUser(-0, 1.);\n\n    TH1 *h_ref = NULL;\n    if (qa_file_ref) {\n      h_ref =\n          (TH1 *)qa_file_ref->GetObjectChecked(prefix + "ntracks_cuts", "TH1");\n      assert(h_pass);\n\n      h_ref->Rebin(nrebin);\n    }\n\n    h_pass->SetTitle(TString(hist_name_prefix) + ": ntracks (#geq 2 MVTX)");\n\n    DrawReference(h_pass, h_ref, false);\n  }\n\n // SaveCanvas(c1,\n //            TString(qa_file_name_new) + TString("_") + TString(c1->GetName()),\n //            true);\n    \n    c1->Draw();\n}\n')


# ## Vertex resolution

# In[28]:


get_ipython().run_cell_magic('cpp', '', '{\n    const char *hist_name_prefix = "QAG4SimulationVertex_SvtxVertexMap";\n    TString prefix = TString("h_") + hist_name_prefix + TString("_");\n    \n    \n  // X-direction\n\n  TH2 *h_new = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("vxRes_gvz"), "TH2");\n  assert(h_new);\n\n  // h_new->Rebin2D(1, 5);\n  //h_new->Sumw2();\n  // h_new->GetXaxis()->SetRangeUser(-15,15);\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("vxRes_gvz"), "TH2");\n    assert(h_ref);\n\n    // h_ref->Rebin2D(1, 5);\n    //h_ref->Sumw2();\n    // h_ref->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c1 = new TCanvas(TString("QA_Draw_Vertex_Resolution_x") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Vertex_Resolution_x") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c1->Divide(2,1);\n  int idx = 1;\n  TPad *p;\n\n  vector<pair<double, double>> gvz_ranges{\n      {-10.0, 10.0}};\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 gvz_range : gvz_ranges)\n  {\n   // cout << __PRETTY_FUNCTION__ << " process " << gvz_range.first << " - " << gvz_range.second << " cm";\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(gvz_range.first + epsilon);\n    const int bin_end = h_new->GetXaxis()->FindBin(gvz_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\t// bin_start, bin_end);\n\n    h_proj_new->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f cm - gvz", gvz_range.first, gvz_range.second));\n    h_proj_new->GetXaxis()->SetTitle(TString::Format(\n        "Vertex Resolution (x) [cm]"));\n    h_proj_new->GetXaxis()->SetNdivisions(5,5);\n    h_proj_new->GetXaxis()->SetRangeUser(-0.002,0.002);\n\n    f1 = new TF1("f1","gaus",-.002,.002);\n    h_proj_new->Fit(f1, "qm");\n    sigma = f1->GetParameter(2);\n    sigma_unc = f1->GetParError(2);\n\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_%d_%d",\n                  h_new->GetName(), bin_start, bin_end));\n              // bin_start, bin_end);\n      h_proj_ref->GetXaxis()->SetRangeUser(-10,10);\n    }\n    \n    DrawReference(h_proj_new, h_proj_ref);\n    sprintf(resstr,"#sigma = %.5f #pm %.5f cm", 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  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  gPad->SetLeftMargin(.2);\n  //h_new->GetYaxis()->SetTitleOffset(2);\n  h_new->Draw("colz");\n  \n\n//  SaveCanvas(c1, TString(qa_file_name_new) + TString("_") + TString(c1->GetName()), true);\n    c1->Draw();\n\n  // Y-direction\n\n  TH2 *h_new2 = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("vyRes_gvz"), "TH2");\n  assert(h_new2);\n\n  //  h_new->Rebin(1, 2);\n  //h_new2->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref2 = NULL;\n  if (qa_file_ref)\n  {\n    h_ref2 = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("vyRes_gvz"), "TH2");\n    assert(h_ref2);\n\n    // h_ref->Rebin(1, 2);\n    //h_ref2->Sumw2();\n    // h_ref->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c2 = new TCanvas(TString("QA_Draw_Vertex_Resolution_y") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Vertex_Resolution_y") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c2->Divide(2,1);\n  int idx2 = 1;\n  TPad *p2;\n\n  vector<pair<double, double>> gvz_ranges2{\n      {-10.0, 10.0}};\n  TF1 *f2 = nullptr;\n  TF1 *fit2 = nullptr;\n  Double_t sigma2 = 0;\n  Double_t sigma_unc2 = 0;\n  char resstr2[500];\n  TLatex *res2 = nullptr;\n  for (auto gvz_range : gvz_ranges2)\n  {\n    //cout << __PRETTY_FUNCTION__ << " process " << gvz_range.first << " - " << gvz_range.second << " cm";\n\n    p2 = (TPad *) c2->cd(idx2++);\n    c2->Update();\n    // p->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new2->GetXaxis()->FindBin(gvz_range.first + epsilon);\n    const int bin_end = h_new2->GetXaxis()->FindBin(gvz_range.second - epsilon);\n\n    TH1 *h_proj_new2 = h_new2->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new2->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n\n    h_proj_new2->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f cm - gvz", gvz_range.first, gvz_range.second));\n    h_proj_new2->GetXaxis()->SetTitle(TString::Format(\n        "Vertex Resolution (y) [cm]"));\n    h_proj_new2->GetXaxis()->SetNdivisions(5,5);\n    h_proj_new2->GetXaxis()->SetRangeUser(-0.002,0.002);\n    \n    f2 = new TF1("f2","gaus",-.002,.002);\n    h_proj_new2->Fit(f2,  "qm");\n    sigma2 = f2->GetParameter(2);\n    sigma_unc2 = f2->GetParError(2);\n\n    TH1 *h_proj_ref2 = nullptr;\n    if (h_ref2)\n    {\n      h_proj_ref2 =\n          h_ref2->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new2->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n    }\n    \n    DrawReference(h_proj_new2, h_proj_ref2);\n    sprintf(resstr2,"#sigma = %.5f #pm %.5f cm", sigma2, sigma_unc2);\n    res2 = new TLatex(0.325,0.825,resstr2);\n    res2->SetNDC();\n    res2->SetTextSize(0.05);\n    res2->SetTextAlign(13);\n    res2->Draw();\n  }\n  p2 = (TPad *) c2->cd(idx2++);\n  c2->Update();\n  gPad->SetLeftMargin(.2);\n  //h_new2->GetYaxis()->SetTitleOffset(2);\n  h_new2->Draw("colz");\n\n  //SaveCanvas(c2, TString(qa_file_name_new) + TString("_") + TString(c2->GetName()), true);\n    c2->Draw();\n    \n  // Z-direction\n\n  TH2 *h_new3 = (TH2 *) qa_file_new->GetObjectChecked(\n      prefix + TString("vzRes_gvz"), "TH2");\n  assert(h_new3);\n\n  //  h_new->Rebin(1, 2);\n  //h_new3->Sumw2();\n  //  h_new->Scale(1. / Nevent_new);\n\n  TH2 *h_ref3 = NULL;\n  if (qa_file_ref)\n  {\n    h_ref3 = (TH2 *) qa_file_ref->GetObjectChecked(\n        prefix + TString("vzRes_gvz"), "TH2");\n    assert(h_ref3);\n\n    // h_ref->Rebin(1, 2);\n    //h_ref3->Sumw2();\n    // h_ref->Scale(Nevent_new / Nevent_ref);\n  }\n\n  TCanvas *c3 = new TCanvas(TString("QA_Draw_Vertex_Resolution_z") + TString("_") + hist_name_prefix,\n                            TString("QA_Draw_Vertex_Resolution_z") + TString("_") + hist_name_prefix,\n                            950, 600);\n  c3->Divide(2,1);\n  int idx3 = 1;\n  TPad *p3;\n\n  vector<pair<double, double>> gvz_ranges3{\n      {-10.0, 10.0}};\n  TF1 *f3 = nullptr;\n  TF1 *fit3 = nullptr;\n  Double_t sigma3 = 0;\n  Double_t sigma_unc3 = 0;\n  char resstr3[500];\n  TLatex *res3 = nullptr;\n  for (auto gvz_range : gvz_ranges3)\n  {\n   // cout << __PRETTY_FUNCTION__ << " process " << gvz_range.first << " - " << gvz_range.second << " cm";\n\n    p3 = (TPad *) c3->cd(idx3++);\n    c3->Update();\n    // p->SetLogy();\n\n    const double epsilon = 1e-6;\n    const int bin_start = h_new3->GetXaxis()->FindBin(gvz_range.first + epsilon);\n    const int bin_end = h_new3->GetXaxis()->FindBin(gvz_range.second - epsilon);\n\n    TH1 *h_proj_new3 = h_new3->ProjectionY(\n        TString::Format(\n            "%s_New_ProjX_%d_%d",\n            h_new3->GetName(), bin_start, bin_end),\n        bin_start, bin_end);\n\n    h_proj_new3->SetTitle(TString(hist_name_prefix) + TString::Format(\n                                                         ": %.1f - %.1f cm -gvz", gvz_range.first, gvz_range.second));\n    h_proj_new3->GetXaxis()->SetTitle(TString::Format(\n        "Vertex Resolution (z) [cm]"));\n    h_proj_new3->GetXaxis()->SetNdivisions(5,5);\n    h_proj_new3->GetXaxis()->SetRangeUser(-0.002,0.002);\n    \n    f3 = new TF1("f3","gaus",-.002,.002);\n    h_proj_new3->Fit(f3,  "qm");\n    sigma3 = f3->GetParameter(2);\n    sigma_unc3 = f3->GetParError(2);\n\n    TH1 *h_proj_ref3 = nullptr;\n    if (h_ref3)\n    {\n      h_proj_ref3 =\n          h_ref3->ProjectionY(\n              TString::Format(\n                  "%s_Ref_ProjX_%d_%d",\n                  h_new3->GetName(), bin_start, bin_end),\n              bin_start, bin_end);\n    }\n    \n    DrawReference(h_proj_new3, h_proj_ref3);\n    sprintf(resstr3,"#sigma = %.5f #pm %.5f cm", sigma3, sigma_unc3);\n    res3 = new TLatex(0.325,0.825,resstr3);\n    res3->SetNDC();\n    res3->SetTextSize(0.05);\n    res3->SetTextAlign(13);\n    res3->Draw();\n  }\n  p3 = (TPad *) c3->cd(idx3++);\n  c3->Update();\n  gPad->SetLeftMargin(.2);\n  //h_new3->GetYaxis()->SetTitleOffset(2);\n  h_new3->Draw("colz");\n\n//  SaveCanvas(c3, TString(qa_file_name_new) + TString("_") + TString(c3->GetName()), true);\n    c3->Draw();\n}\n')


# ## Vertex resolution w/ multiplicity dependence
# 
# To be completed

# In[ ]:





# # Summary statistics

# In[29]:


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


# In[30]:


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