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

# # Introduction
# 
# QA plots for the single particle in calorimeters and reconstructed as calorimetric jet objects

# 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]:


# initialization

display(Markdown('pyROOT env check:'))

import ROOT

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


# In[4]:


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


# In[5]:


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}"))

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[6]:


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[7]:


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[8]:


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


# ## Inputs and file checks

# In[9]:


qa_file_name_new = os.getenv("qa_file_name_new")
if qa_file_name_new is None:
    qa_file_name_new = "G4sPHENIX_pi+_pT30_Sum15_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_pi+_pT30_Sum15_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[10]:


# 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[11]:


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')


# In[12]:


get_ipython().run_cell_magic('cpp', '-d', '\n// Jet utility functions\n\n\nvoid QA_Draw_Jet_Spectrum(\n    const char *jet = "h_QAG4SimJet_AntiKt_Truth_r07")\n{\n\n  // obtain normalization\n  double Nevent_new = 1;\n  double Nevent_ref = 1;\n\n  if (qa_file_new)\n  {\n    TH1D *h_norm = (TH1D *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Normalization"), "TH1D");\n    assert(h_norm);\n\n    Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));\n  }\n  if (qa_file_ref)\n  {\n    TH1D *h_norm = (TH1D *) qa_file_ref->GetObjectChecked(\n        TString(jet) + TString("_Normalization"), "TH1D");\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_Jet_Spectrum_") + TString(jet),\n                            TString("QA_Draw_Jet_Spectrum_") + TString(jet), 950, 600);\n  c1->Divide(4, 2);\n  int idx = 1;\n  TPad *p;\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_eta"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_eta"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_phi"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_phi"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_Et"), "TH1F");\n    assert(h_new);\n\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_Et"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_Mass"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_Mass"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogx();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_CompSize"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_CompSize"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_CEMC_Ratio"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_CEMC_Ratio"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  p->SetLogy();\n\n  {\n    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n        TString(jet) + TString("_Leading_CEMC_HCalIN_Ratio"), "TH1F");\n    assert(h_new);\n\n    h_new->Rebin(2);\n    // h_new->Sumw2();\n    h_new->Scale(1. / Nevent_new);\n\n    TH1F *h_ref = NULL;\n    if (qa_file_ref)\n    {\n      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + TString("_Leading_CEMC_HCalIN_Ratio"), "TH1F");\n      assert(h_ref);\n\n      h_ref->Rebin(2);\n      h_ref->Scale(1. / Nevent_ref);\n    }\n\n    h_new->GetYaxis()->SetTitleOffset(1.5);\n    h_new->GetYaxis()->SetTitle("Count / event / bin");\n    //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n    DrawReference(h_new, h_ref);\n  }\n\n  if (TString(jet).Contains("Truth"))\n  {\n    // truth jets\n\n    p = (TPad *) c1->cd(idx++);\n    c1->Update();\n    p->SetLogy();\n\n    {\n      TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(\n          TString(jet) + TString("_Leading_Leakage_Ratio"), "TH1F");\n      assert(h_new);\n\n      h_new->Rebin(2);\n      // h_new->Sumw2();\n      h_new->Scale(1. / Nevent_new);\n\n      TH1F *h_ref = NULL;\n      if (qa_file_ref)\n      {\n        h_ref = (TH1F *) qa_file_ref->GetObjectChecked(\n            TString(jet) + TString("_Leading_Leakage_Ratio"), "TH1F");\n        assert(h_ref);\n\n        h_ref->Rebin(2);\n        h_ref->Scale(1. / Nevent_ref);\n      }\n\n      h_new->GetYaxis()->SetTitleOffset(1.5);\n      h_new->GetYaxis()->SetTitle("Count / event / bin");\n      //      h_new->GetXaxis()->SetRangeUser(-0, .1);\n\n      DrawReference(h_new, h_ref);\n    }\n  }\n\n//  PutInputFileName(c1, .04, qa_file_name_new, qa_file_name_ref);\n    c1->Draw();\n    \n    //  SaveCanvas(c1, TString(qa_file_name_new) + TString(c1->GetName()), false);\n    \n}\n')


# # $R=0.7$ Jets spectrums
# 
# Check the particular set of jets with $R=0.7$ parameters

# ## Truth jet spectrum
# 
# Jet the details of the truth jet. Since only one particle is used in this simulation, the number of component is 1 and jet mass is 0.

# In[13]:


get_ipython().run_cell_magic('cpp', '', '\nQA_Draw_Jet_Spectrum("h_QAG4SimJet_AntiKt_Truth_r07");\n')


# ## Reco tower jet spectrum
# 
# Jet the details of the reco tower jet in response to this single particle

# In[14]:


get_ipython().run_cell_magic('cpp', '', '\nQA_Draw_Jet_Spectrum("h_QAG4SimJet_AntiKt_Tower_r07");\n')


# # Jet-Truth matching
# 
# Bunch of plots matching reco jet to truth jet

# In[15]:


get_ipython().run_cell_magic('cpp', '', '\n//! jet family, tower, cluster or track jets \nconst char *jet_family = "AntiKt_Tower";\n\n//! drawing energy range\nconst double min_Et = 10;\nconst double max_Et = 80;\n')


# In[16]:


get_ipython().run_cell_magic('cpp', '-d', '\n// utility of drawing a matching pair\n\nvector<TGraphErrors *>\nQA_Draw_Jet_TruthMatching(const char *jet =\n                              "h_QAG4SimJet_AntiKt_Truth_r07_AntiKt_Tower_r07")\n{\n\n  vector<TGraphErrors *> resolution_collections;\n\n  TCanvas *c1 = new TCanvas(\n      TString("QA_Draw_Jet_TruthMatching_") + TString(jet),\n      TString("QA_Draw_Jet_TruthMatching_") + TString(jet), 950, 600);\n  c1->Divide(3, 2);\n  int idx = 1;\n  TPad *p;\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *proj_new = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_dPhi", "TH2F");\n\n    assert(proj_new);\n\n    proj_new->Rebin2D(1, 5);\n\n    TGraphErrors *ge = FitProfile(proj_new);\n\n    proj_new->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    proj_new->GetYaxis()->SetTitleOffset(1.5);\n    proj_new->Draw("COLZ");\n\n    TGraphErrors *ge_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *proj_ref = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_dPhi", "TH2F");\n      assert(proj_ref);\n      proj_ref->Rebin2D(1, 5);\n      TGraphErrors *ge_ref = FitProfile(proj_ref);\n    }\n    DrawReference(ge, ge_ref);\n\n    resolution_collections.push_back(ge);\n  }\n  TLine *l = new TLine(min_Et, 0, max_Et, 00);\n  l->Draw();\n    \n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *proj_new = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_dEta", "TH2F");\n\n    assert(proj_new);\n\n    proj_new->Rebin2D(1, 5);\n    TGraphErrors *ge = FitProfile(proj_new);\n\n    proj_new->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    proj_new->GetYaxis()->SetTitleOffset(1.5);\n    proj_new->Draw("COLZ");\n    proj_new->SetTitle(jet);\n\n    TGraphErrors *ge_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *proj_ref = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_dEta", "TH2F");\n      assert(proj_ref);\n      proj_ref->Rebin2D(1, 5);\n      TGraphErrors *ge_ref = FitProfile(proj_ref);\n    }\n    DrawReference(ge, ge_ref);\n\n    resolution_collections.push_back(ge);\n  }\n  l = new TLine(min_Et, 0, max_Et, 00);\n  l->Draw();\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *proj_new = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_dE", "TH2F");\n\n    assert(proj_new);\n\n    //    proj_new->Rebin2D(1,5);\n\n    TGraphErrors *ge = FitProfile(proj_new);\n\n    proj_new->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    proj_new->GetYaxis()->SetTitleOffset(1.5);\n    proj_new->Draw("COLZ");\n    proj_new->SetTitle(jet);\n\n    TGraphErrors *ge_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *proj_ref = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_dE", "TH2F");\n      assert(proj_ref);\n      proj_ref->Rebin2D(1, 5);\n      TGraphErrors *ge_ref = FitProfile(proj_ref);\n    }\n    DrawReference(ge, ge_ref);\n\n    resolution_collections.push_back(ge);\n  }\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *proj_new = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_dEt", "TH2F");\n\n    assert(proj_new);\n\n    //    proj_new->Rebin2D(1,5);\n    TGraphErrors *ge = FitProfile(proj_new);\n\n    proj_new->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    proj_new->GetYaxis()->SetTitleOffset(1.5);\n    proj_new->Draw("COLZ");\n    proj_new->SetTitle(jet);\n\n    TGraphErrors *ge_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *proj_ref = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_dEt", "TH2F");\n      assert(proj_ref);\n      proj_ref->Rebin2D(1, 5);\n      TGraphErrors *ge_ref = FitProfile(proj_ref);\n    }\n    DrawReference(ge, ge_ref);\n\n    resolution_collections.push_back(ge);\n  }\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *h2 = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_Count_Truth_Et", "TH2F");\n    assert(h2);\n\n    TH1 *h_norm = h2->ProjectionX(\n        TString(jet) + "_Matching_Count_Truth_Et" + "_All", 1, 1);\n    //          TH1 * h_pass = h2->ProjectionX(\n    //              TString(jet) + "_Matching_Count_Truth_Et" + "_Matched", 2, 2);// inclusive match\n    TH1 *h_pass = h2->ProjectionX(\n        TString(jet) + "_Matching_Count_Truth_Et" + "_Matched", 3, 3);  // unique match\n    assert(h_norm);\n    assert(h_pass);\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.2);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *h2 = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_Count_Truth_Et", "TH2F");\n      assert(h2);\n      TH1 *h_norm = h2->ProjectionX(\n          TString(jet) + "_Matching_Count_Truth_Et" + "_All", 1, 1);\n      //          TH1 * h_pass = h2->ProjectionX(\n      //              TString(jet) + "_Matching_Count_Truth_Et" + "_Matched", 2, 2);// inclusive match\n      TH1 *h_pass = h2->ProjectionX(\n          TString(jet) + "_Matching_Count_Truth_Et" + "_Matched", 3, 3);  // unique match\n      assert(h_norm);\n      assert(h_pass);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    DrawReference(h_ratio, h_ratio_ref, true);\n\n    resolution_collections.push_back(new TGraphErrors(h_ratio));\n  }\n\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  {\n    TH2F *h2 = (TH2F *) qa_file_new->GetObjectChecked(\n        TString(jet) + "_Matching_Count_Reco_Et", "TH2F");\n    assert(h2);\n\n    TH1 *h_norm = h2->ProjectionX(\n        TString(jet) + "_Matching_Count_Reco_Et" + "_All", 1, 1);\n    //      TH1 * h_pass = h2->ProjectionX(\n    //          TString(jet) + "_Matching_Count_Reco_Et" + "_Matched", 2, 2); // inclusive match\n    TH1 *h_pass = h2->ProjectionX(\n        TString(jet) + "_Matching_Count_Reco_Et" + "_Matched", 3, 3);  // unique match\n    assert(h_norm);\n    assert(h_pass);\n    TH1 *h_ratio = GetBinominalRatio(h_pass, h_norm);\n\n    h_ratio->GetXaxis()->SetRangeUser(min_Et, max_Et);\n    h_ratio->GetYaxis()->SetTitle("Reconstruction Purity");\n    h_ratio->GetYaxis()->SetRangeUser(-0, 1.2);\n\n    TH1 *h_ratio_ref = NULL;\n    if (qa_file_ref)\n    {\n      TH2F *h2 = (TH2F *) qa_file_ref->GetObjectChecked(\n          TString(jet) + "_Matching_Count_Reco_Et", "TH2F");\n      assert(h2);\n\n      TH1 *h_norm = h2->ProjectionX(\n          TString(jet) + "_Matching_Count_Reco_Et" + "_All", 1, 1);\n      //      TH1 * h_pass = h2->ProjectionX(\n      //          TString(jet) + "_Matching_Count_Reco_Et" + "_Matched", 2, 2); // inclusive match\n      TH1 *h_pass = h2->ProjectionX(\n          TString(jet) + "_Matching_Count_Reco_Et" + "_Matched", 3, 3);  // unique match\n      assert(h_norm);\n      assert(h_pass);\n      h_ratio_ref = GetBinominalRatio(h_pass, h_norm);\n    }\n\n    DrawReference(h_ratio, h_ratio_ref, true);\n\n    resolution_collections.push_back(new TGraphErrors(h_ratio));\n  }\n\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n\n//  PutInputFileName(c1, .04, qa_file_name_new, qa_file_name_ref);\n    c1->Draw();\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString(c1->GetName()), false);\n    \n  return resolution_collections;\n}\n')


# In[17]:


get_ipython().run_cell_magic('cpp', '', '{\n//      qa_file_new->GetListOfKeys().Print()\n//     qa_file_new->ls();\n//    QA_Draw_Jet_TruthMatching("h_QAG4SimJet_AntiKt_Truth_r07_AntiKt_Tower_r07");\n\n}\n')


# ## Kinematic dependent Jet matching plots for various radius

# In[18]:


get_ipython().run_cell_magic('cpp', '', '\n// buffer for results\nvector<float> vec_radius;\nvector<TGraphErrors *> vec_phi_res;\nvector<TGraphErrors *> vec_eta_res;\nvector<TGraphErrors *> vec_e_res;\nvector<TGraphErrors *> vec_et_res;\nvector<TGraphErrors *> vec_reco_eff;\nvector<TGraphErrors *> vec_purity;\n\n{\n    \n  // list and process all jets\n  TList *hist_key_list = qa_file_new->GetListOfKeys();\n    TString s_re_fullname = Form(\n        "h_QAG4SimJet_.*_r[0-9]*_%s_r[0-9]*_Matching_Count_Truth_Et",\n        jet_family);  // regular expression for search\n    TRegexp re_fullname(s_re_fullname, false);\n    \n  cout <<"Seraching for jet pairs in "<<hist_key_list->GetSize()<<" QA histograms with pattern "<<s_re_fullname<<endl;\n  for (int i = 0; i < hist_key_list->GetSize(); ++i)\n  {\n    TString key_name = hist_key_list->At(i)->GetName();\n\n    if (key_name.Index(re_fullname) == kNPOS)\n      continue;\n\n    //      cout << " key_name = " << key_name << endl;\n    TString jet_pair_name = key_name(0,\n                                     key_name.Length() - TString("_Matching_Count_Truth_Et").Length());  // remove suffix\n\n    // cout << "Processing jet_pair_name = " << jet_pair_name << endl;\n\n    //get jet radius\n    TRegexp re_jetradius("_r[0-9]*", false);\n    Ssiz_t index_radius = key_name.Index(re_jetradius);             // first radius\n    index_radius = key_name.Index(re_jetradius, index_radius + 1);  // second radius\n    assert(index_radius != kNPOS);\n    float radius = 0;\n    sscanf(key_name(index_radius, 100).Data(), "_r%f", &radius);\n    //      cout << " index_radius = " << index_radius << endl;\n    assert(radius != 0);\n    radius /= 10;  // jet radius convention in DST names\n\n    cout << "QA_Draw_Jet_Summary - process jet pair " << jet_pair_name\n         << " with radius = " << radius << endl;\n\n    vector<TGraphErrors *> resolution_efficiency_summary(\n        QA_Draw_Jet_TruthMatching(jet_pair_name));\n\n    //save results\n    vec_radius.push_back(radius);\n    vec_phi_res.push_back(resolution_efficiency_summary[0]);\n    vec_eta_res.push_back(resolution_efficiency_summary[1]);\n    vec_e_res.push_back(resolution_efficiency_summary[2]);\n    vec_et_res.push_back(resolution_efficiency_summary[3]);\n    vec_reco_eff.push_back(resolution_efficiency_summary[4]);\n    vec_purity.push_back(resolution_efficiency_summary[5]);\n\n  }\n    \n}\n')


# ## Summary plot as function of jet parameter

# In[19]:


get_ipython().run_cell_magic('cpp', '', '\n{\n    \n  // plot\n  TCanvas *c1 = new TCanvas(\n      TString("QA_Draw_Jet_Summary_") + TString(jet_family),\n      TString("QA_Draw_Jet_Summary_") + TString(jet_family), 950, 600);\n  c1->Divide(3, 2);\n  int idx = 1;\n  TPad *p;\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  TH1 *h_frame =\n      p->DrawFrame(min_Et, -.1, max_Et, .1,\n                   TString(jet_family) + " #phi Reconstruction;E_{T, Truth} (GeV);#phi_{Reco} - #phi_{Truth} (rad)");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  TLine *l = new TLine(min_Et, 0, max_Et, 0);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  TLegend *legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_phi_res[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  h_frame =\n      p->DrawFrame(min_Et, -.1, max_Et, .1,\n                   TString(jet_family) + " #eta Reconstruction;E_{T, Truth} (GeV);#eta_{Reco} - #eta_{Truth}");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  l = new TLine(min_Et, 0, max_Et, 0);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_eta_res[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  h_frame = p->DrawFrame(min_Et, 0, max_Et, 2,\n                         TString(jet_family) + " Jet Energy Reconstruction;E_{Truth} (GeV);E_{Reco} / E_{Truth}");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_e_res[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  h_frame =\n      p->DrawFrame(min_Et, 0, max_Et, 2,\n                   TString(jet_family) + " Jet E_{T} Reconstruction;E_{T, Truth} (GeV);E_{T, Reco} / E_{T, Truth}");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_et_res[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  h_frame = p->DrawFrame(min_Et, 0, max_Et, 1.2,\n                         TString(jet_family) + " Reco Efficiency;E_{T, Truth} (GeV);Reco efficiency");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_reco_eff[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  // ------------------------------------\n  p = (TPad *) c1->cd(idx++);\n  c1->Update();\n  //  p->SetLogz();\n\n  h_frame = p->DrawFrame(min_Et, 0, max_Et, 1.2,\n                         TString(jet_family) + " Reconstruction Purity;E_{T, Reco} (GeV);Reconstruction Purity");\n  //  h_frame->GetYaxis()->SetTitleOffset(1.01);\n  l = new TLine(min_Et, 1, max_Et, 1);\n  l->Draw();\n  p->SetGridx(0);\n  p->SetGridy(0);\n  legend = new TLegend(0.7, 0.2, .95, 0.5);\n  legend->SetFillColor(kWhite);\n  legend->SetFillStyle(1001);\n  legend->SetLineWidth(2);\n  legend->SetLineColor(kBlack);\n  legend->SetLineStyle(kSolid);\n  for (int i = 0; i < vec_radius.size(); ++i)\n  {\n    const float radius = vec_radius[i];\n\n    TGraphErrors *ge = vec_purity[i];\n    assert(ge);\n    ge = new TGraphErrors(*ge);  // make a copy\n\n    ge->SetLineColor(i + 2);    // automatic color scheme from ROOT\n    ge->SetMarkerColor(i + 2);  // automatic color scheme from ROOT\n    for (int idata = 0; idata < ge->GetN(); ++idata)\n    {\n      (ge->GetX())[idata] += i * 0.5;  // shift x a little bit\n      (ge->GetEX())[idata] = 0;        // no x error bar\n    }\n    ge->Draw("p E l");\n    legend->AddEntry(ge, Form("r = %.1f", radius), "elp");\n  }\n  legend->Draw();\n\n  //  PutInputFileName(c1, .03, qa_file_name_new, qa_file_name_ref);\n    c1->Draw();\n  //SaveCanvas(c1, TString(qa_file_name_new) + TString(c1->GetName()), true);\n\n}\n')


# In[ ]:





# # Summary statistics

# In[20]:


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


# In[21]:


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