#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Author: Elena Graverini
# @Date:   2015-10-27 18:26:49
# @Last Modified by:   Elena Graverini
# @Last Modified time: 2017-05-03 18:58:44

import os
from math import sqrt
import ROOT as r

location = {'TT': "/disk/data1/hep/elena/data/ST/Aging_Tuples/TT/TTaU/",
            'IT': "/disk/data1/hep/elena/data/ST/Aging_Tuples/IT/T3X2/"}
out_location = os.getcwd() + '/trackSel'
os.system('mkdir -p %s' % out_location)

# Load fill numbers
macros = os.path.expandvars('$CCEHOME/macros/CCEScan/')
with open(macros + 'Fills.dat', 'rb') as f:
    fills = f.read().splitlines()

fills = [3478, 4518, 4643, 4856, 5162, 5448]
dets = ['TT', 'IT']
val = {'TT': 5, 'IT': 7}
layer = {'TT': 'TTaU', 'IT': 'T3X2'}
cut = {'TT': 12, 'IT': 20}
stat = 10000000

r.gROOT.SetBatch(True)
r.gROOT.ProcessLine('.X %s/../include/lhcbstyle.C' % os.getcwd())
# r.gStyle.SetPadRightMargin(0.16)
# r.gStyle.SetPadLeftMargin(r.gStyle.GetPadLeftMargin() * 0.7)
r.gStyle.SetPadLeftMargin(0.20)
r.gStyle.SetPadRightMargin(0.25)
r.gStyle.SetPadBottomMargin(0.15)
r.gStyle.SetPadTopMargin(0.05)


def sqrt_of_2d_histogram(h):
    new = h.Clone()
    nbins = h.GetNbinsX() * h.GetNbinsY()
    for i in range(nbins):
        new.SetBinContent(i, sqrt(h.GetBinContent(i)))
    return new


def s_over_b(s, b):
    h = s.Clone()
    res = h.Divide(b)
    res.name = 's_over_b'
    return res


def s_over_sqrtb(s, b):
    h = s.Clone()
    sqrtb = sqrt_of_2d_histogram(b)
    res = h.Divide(sqrtb)
    res.name = 's_over_sqrtb'
    return res


def punzi_significance(s, b):
    h = s.Clone()
    sb = s.Add(b)
    sqrt_sb = sqrt_of_2d_histogram(sb)
    res = h.Divide(sqrt_sb)
    res.name = 'punzi_significance'
    return res


def put_fill(c, fill):
    label = r.TPaveText( 0.74 - r.gStyle.GetPadRightMargin(), 0.87 - r.gStyle.GetPadTopMargin(),
                         0.95 - r.gStyle.GetPadRightMargin(), 0.95 - r.gStyle.GetPadTopMargin(), "BRNDC")
    label.SetFillColor(0)
    label.SetTextAlign(12)
    label.SetBorderSize(0)
    label.SetTextFont(132)
    label.SetTextSize(0.06)
    label.SetTextAlign(22)
    label.AddText('Fill %s' % fill)
    label.Draw('same')
    c.Modified()
    c.Update()
    return c


def printSelection(f, t_sig, det, fill, stat):
    r.gStyle.SetPadLeftMargin(0.20)
    r.gStyle.SetPadRightMargin(0.25)
    r.gStyle.SetPadBottomMargin(0.15)
    r.gStyle.SetPadTopMargin(0.05)
    print t_sig.Draw("GhostP : TrChi2/TrNDoF>>hsig(100,0., 5,100, 0.,1.)", "val%s>%s" % (val[det], cut[det]), "goff", stat, 0)
    h_sig = r.gDirectory.Get('hsig')
    print t_sig.Draw("GhostP : TrChi2/TrNDoF>>hnoise(100,0., 5,100, 0.,1.)", "val%s<%s" % (val[det], cut[det]), "goff", stat, 0)
    h_noise = r.gDirectory.Get('hnoise')
    # c = r.TCanvas()
    csaver = []
    h_s_over_b = s_over_b(h_sig, h_noise)
    h_s_over_b.GetZaxis().SetRangeUser(0, 10)
    h_s_over_b.GetZaxis().SetTitle('S/B ratio')
    h_s_over_sqrtb = s_over_sqrtb(h_sig, h_noise)
    h_s_over_sqrtb.GetZaxis().SetRangeUser(0, 50)
    h_s_over_sqrtb.GetZaxis().SetTitle('S/#sqrt{B}')
    h_punzi_significance = punzi_significance(h_sig, h_noise)
    h_punzi_significance.GetZaxis().SetRangeUser(0, 10)
    h_punzi_significance.GetZaxis().SetTitle('S/#sqrt{S+B}')
    for h in [h_s_over_b, h_s_over_sqrtb, h_punzi_significance]:
        csaver.append(r.TCanvas())
        csaver[-1].cd()
        h.GetXaxis().SetTitle(r'#chi^{2}_{track}/ndf')
        h.GetYaxis().SetTitle(r'Ghost Prob.')
        h.Draw('colz')
        csaver[-1].Update()
        csaver[-1], line = putLine(c, det)
        csaver[-1] = put_fill(fill)
        csaver[-1].SaveAs('%s/%s_%s_fill%s.pdf' % (out_location, det, h.name, fill))
        h.Delete()
    # h_sig.Divide(h_noise)
    # h_sig.GetZaxis().SetRangeUser(0, 10)
    # h_sig.SetTitle('TTaU S/B ratio for fill %s' % fill)
    # h_sig.GetXaxis().SetTitle(r'#chi^{2}_{track}/ndf')
    # h_sig.GetYaxis().SetTitle(r'Ghost Prob.')
    # h_sig.GetZaxis().SetTitle(r'S/B ratio')
    # h_sig.Draw("colz")
    # c.Update()
    # c, line = putLine(c, det)
    # c = put_fill(c, fill)
    # c.SaveAs('%s/%s_fill%s.pdf' % (out_location, det, fill))
    h_sig.Delete()
    h_noise.Delete()
    return f, t_sig


def printADC(f, t_sig, det, fill, stat, step):
    if step > 10:
        print t_sig.Draw("val%s>>hsig%s(210,-35, 175)" % (val[det], step), "odinStep==%s" % step, "goff")  # , stat, 0)
    else:
        print t_sig.Draw("val%s>>hsig%s(210,-35, 175)" % (val[det], step), "odinStep==%s" % step, "goff", stat, 0)
    h_sig = r.gDirectory.Get('hsig%s' % step)
    r.gROOT.ProcessLine('.X %s/../include/lhcbstyle.C' % os.getcwd())
    r.gStyle.SetPadLeftMargin(0.15)
    c = r.TCanvas()
    h_sig.GetXaxis().SetTitle('Signal height [ADC value]')
    h_sig.GetYaxis().SetTitle('Fraction of hits / ADC value')
    h_sig.DrawNormalized()
    c.Update()
    label = r.TPaveText( 0.74 - r.gStyle.GetPadRightMargin(), 0.80 - r.gStyle.GetPadTopMargin(),
                         0.95 - r.gStyle.GetPadRightMargin(), 0.95 - r.gStyle.GetPadTopMargin(), "BRNDC")
    label.SetFillColor(0)
    label.SetTextAlign(12)
    label.SetBorderSize(0)
    label.SetTextFont(132)
    label.SetTextSize(0.06)
    label.SetTextAlign(22)
    label.AddText('Fill %s\n' % fill)
    from checkPulseData import get_vmap
    v_bias = get_vmap(det)[int(step) / 6]
    label.AddText('V_{bias} = %s V' % v_bias)
    label.Draw('same')
    line = r.TLine()
    line.SetLineWidth(3)
    line.SetLineColor(r.kRed)
    line.SetLineStyle(7)
    line.DrawLine(float(cut[det]), 0.0, float(cut[det]), r.gPad.GetUymax())
    c.Modified()
    c.Update()
    c.SaveAs('%s/%s_odinStep%s_fill%s.pdf' % (out_location, det, step, fill))
    h_sig.Delete()
    return f, t_sig


def putLine(c, det):
    line = r.TLine()
    line.SetLineWidth(4)
    line.SetLineColor(r.kWhite)
    c.cd()
    if 'TT' in det:
        line.DrawLine(0.0, 0.01, 1.6, 0.16)
        line.DrawLine(3.0, 0.07, 1.6, 0.16)
        line.DrawLine(3.0, 0.07, 3.0, 0.0)
    else:
        line.DrawLine(0.0, 0.0, 3.5, 0.35)
        line.DrawLine(3.5, 0.35, 5.0, 0.20)
    return c, line


def printByStep(fill, stat):
    f = r.TFile(location_tt + '%s.root' % fill, 'read')
    t_sig = f.Get('STADCTrackMonitor/HitInfo/TTaU')

    # Look at Landau distribution
    for step in range(67):
        c = r.TCanvas()
        t_sig.Draw('val5>>hnew(70,-30,90)', 'odinStep==%s' % step)
        c.SaveAs()

    draw_sig = "GhostP : TrChi2/TrNDoF>>hsig(100, 0., 10., 100, 0., 1.)"
    draw_bg = "GhostP : TrChi2/TrNDoF>>hnoise(100, 0., 10., 100, 0., 1.)"
    f.Close()


if __name__ == '__main__':
    for det in dets:
        for fill in fills:
            f = r.TFile(location[det] + '%s.root' % fill, 'read')
            t_sig = f.Get('STADCTrackMonitor/HitInfo/%s' % layer[det])

            f, t_sig = printSelection(f, t_sig, det, fill, stat)
            for step in [0, 12, 48, 54, 60]:
                f, t_sig = printADC(f, t_sig, det, fill, stat, step)

            f.Close()