#!/usr/bin/env python2
import sys, getopt, os
import array
import datetime
import csv
import math


maxY = 0
#maskcut = "!(mask & (4+8+16+1))" #4=inBleedZone, 8=nearBigEvent, 16=crossTalk, 1=hasNeighbor
#maskcut = "!(mask & (4+8+16+128))" #4=inBleedZone, 8=nearBigEvent, 16=crossTalk, 128=serial register hit
#maskcut = "!(mask & (1+4+8+16+64+128+256+512+1024+4096+8192))" #4=inBleedZone, 8=nearBigEvent, 16=crossTalk, 128=serial register hit
maskbits = 4+8+16 #4=inBleedZone, 8=nearBigEvent, 16=crossTalk
outfilename=""

options, remainder = getopt.gnu_getopt(sys.argv[1:], 'y:o:m:h')

for opt, arg in options:
    if opt == '-y':
        maxY = int(arg)
    elif opt == '-o':
        outfilename = arg
    elif opt == '-m':
        maskbits |= int(arg)
    elif opt == '-h':
        print("\nUsage: "+sys.argv[0]+" <output basename> <root files>")
        print("Arguments: ")
        print("\t-y: max-Y cut to apply to data")
        print("\t-o: basename for output file")
        print("\t-m: additional mask bits")
        print("\n")
        sys.exit(0)

#import ROOT now, since otherwise PyROOT eats the command-line options
#from ROOT import gROOT, gStyle, TFile, TTree, TChain, TCanvas, gDirectory, TH1, TGraph, gPad, TF1, THStack, TLegend, TGraphErrors, TLatex, TEfficiency, TMath
from ROOT import gROOT
gROOT.SetBatch(True)
from ROOT import gStyle
from ROOT import TFile
from ROOT import TTree
from ROOT import TChain
from ROOT import TCanvas
from ROOT import gDirectory
from ROOT import TH1
from ROOT import TGraph
from ROOT import gPad
from ROOT import TF1
from ROOT import THStack
from ROOT import TLegend
from ROOT import TGraphErrors
from ROOT import TLatex
from ROOT import TEfficiency
from ROOT import TMath
from ROOT import TVectorD
from ROOT import nullptr

import skipper_utils

gStyle.SetOptStat(110011)
gStyle.SetOptFit(1)
#gStyle.SetPalette(57)

if (len(remainder)<1):
    print(sys.argv[0]+' <root files>')
    sys.exit()

infiles = remainder
#if len(infiles)>1:
#    infiles.sort(key=skipper_utils.decodeRunnum)

#if (len(sys.argv)<3):
#    print("not enough input args")
#    sys.exit(1)
if outfilename=="":
    basenames = [os.path.splitext(os.path.basename(a))[0] for a in infiles] #drop the directory path and file extension from each filename
    outfilename = "ana_"+os.path.commonprefix(basenames) #take the common leading substring
    print("no output filename supplied, using default: "+outfilename)

csvfile = open(outfilename+'.csv','w')
csvFields = ['runid', 'ltaname', 'readoutStart', 'readoutEnd', 'exposure']
csvDict = {}

data = TChain("calPixTree")

#default dimensions - for prototype CCD
ccdPrescan = 7
ccdCol = 724
ccdRow = 1248
nRow = 700
nCol = 450
binRows = 1
binCols = 1
skipRows = 0
skipCols = 0

numRuns = 0
readoutDays = 1.0
for i in range(0,len(infiles)):
    print(infiles[i])
    data.Add(infiles[i])
    numRuns += 1
    thefile = TFile(infiles[i])
    header = thefile.Get("headerTree_0")
    if i==0: #get CCD dimensions from the first file
        ccdPrescan = int(skipper_utils.getHeaderValue(header, "CCDNPRES"))
        ccdCol = int(skipper_utils.getHeaderValue(header, "CCDNCOL"))
        ccdRow = int(skipper_utils.getHeaderValue(header, "CCDNROW"))
        nRow = int(skipper_utils.getHeaderValue(header, "NROW"))
        nCol = int(skipper_utils.getHeaderValue(header, "NCOL"))
        try:
            binRows = int(skipper_utils.getHeaderValue(header, "NBINROW"))
        except:
            binRows = 1
        try:
            binCols = int(skipper_utils.getHeaderValue(header, "NBINCOL"))
        except:
            binCols = 1
        try:
            skipRows = int(skipper_utils.getHeaderValue(header, "SKIPROW"))
        except:
            skipRows = 0
        try:
            skipCols = int(skipper_utils.getHeaderValue(header, "SKIPCOL"))
        except:
            skipCols = 0
        start = datetime.datetime.strptime(skipper_utils.getHeaderString(header, "DATESTART"),"%Y-%m-%dT%H:%M:%S")
        end = datetime.datetime.strptime(skipper_utils.getHeaderString(header, "DATEEND"),"%Y-%m-%dT%H:%M:%S")
        runID = int(skipper_utils.getHeaderValue(header, "RUNID"))
        try:
            ltaName = skipper_utils.getHeaderValue(header, "LTANAME")
        except:
            ltaName = ""
        csvDict['runid'] = runID
        csvDict['ltaname'] = ltaName
        csvDict['readoutStart'] = start
        csvDict['readoutEnd'] = end
        try:
            exposure = int(skipper_utils.getHeaderValue(header, "EXPOSURE"))
            csvDict['exposure'] = exposure
        except:
            exposure = 0
        readoutTime = end-start
        readoutDays = readoutTime.days + (float(readoutTime.seconds) - exposure)/(24*60*60)
    #ogl = skipper_utils.getHeaderValue(header,"OGAL")
    #swl = skipper_utils.getHeaderValue(header,"SWAL")
    #print("OGL={0}, SWL={1}".format(ogl,swl))

daysPerRow = readoutDays/nRow

if maxY!=0:
    nRow = min(maxY,nRow)

#x=0 and y=0 look weird (negative or extreme values)
#for prototype CCD:
#x=[370,450] is overscan
#x=[1,7] is prescan
#x=[8,369] is active area - matches dimension in paper (362 pixels)
#y=[625,700] is vertical overscan
#y=[1,624] is active area - matches dimension in paper (624 pixels)
#y=624 (top row of active area) seems to have excess charge, so we exclude it
#x=8 (first column of active area) seems to have excess charge, so we exclude it

prescanEdge = (ccdPrescan + 1 - skipCols)/binCols #X of first active pixel
activeEdgeX = (ccdPrescan + ccdCol/2 - skipCols)/binCols #X of last active pixel

#activeEdgeY = (ccdRow/2)/binRows #Y of last active pixel
if nRow*binRows>ccdRow:
    activeEdgeY = ccdRow/binRows
else:
    activeEdgeY = (ccdRow/2)/binRows #Y of last active pixel

maskcut = "!(mask & {0})".format(maskbits)
activecutX = "x>{0} && x<={1}".format(prescanEdge+1, activeEdgeX)
activecutY = "y>0 && y<={0}".format(activeEdgeY)
edgecut = "x!={0} && y!={1}".format(prescanEdge, activeEdgeY)

#pick bin sizes and ranges for efficiency plots
binSizeX = 20
nbinsX = (nCol/binSizeX) + 1
minX = prescanEdge - binSizeX - 0.5
maxX = minX + binSizeX + nCol
binningX = "{0},{1},{2}".format(nbinsX, minX, maxX)

activebinsY = 16 #16 is the common divisor of proto, science, and C sizes
binSizeY = ccdRow/2/activebinsY
datarows = min(nRow, activeEdgeY)
activerows = min(nRow, activeEdgeY)
activebins = activerows/binSizeY
while activebins<20:
    if binSizeY%2==0:
        activebinsY *= 2
        binSizeY /= 2
        activebins = activerows/binSizeY
    else:
        break
minY = 0.5
nbinsY = ((nRow-1)/binSizeY)+1
topbinedgeY = nbinsY*binSizeY + 0.5
binningY = "{0},{1},{2}".format(nbinsY, minY, topbinedgeY)
#topbinedgeY = (databins+1)*binSizeY + 0.5

fitfunc = skipper_utils.poissonFitfunc()

pixmax = 50.0
pixmin = -1.0
pixnbin = 510
pixbinning = "{0},{1},{2}".format(pixnbin,pixmin,pixmax)
binwidth = (pixmax-pixmin)/pixnbin

c = TCanvas("c","c",1200,900);
c.Print(outfilename+".pdf[")

outfile = TFile(outfilename+".root","RECREATE")

latex = TLatex()
latex.SetNDC(True)

c.SetLogy(1)

fitvals_default = {
        "zero" : 0.0,
        "gain" : 1.0,
        "noise" : 0.1
        }

ohdus = []
ohduelists = []

for ohdu in range(0,6):#LTA data has 1 through 4, Monsoon data has 2 through 5
    elistname = "e{0}".format(ohdu)
    if maxY!=0:
        n = data.Draw(">>"+elistname,"ohdu=={0} && {1} && y<={2}".format(ohdu, maskcut, maxY))
    else:
        n = data.Draw(">>"+elistname,"ohdu=={0} && {1}".format(ohdu, maskcut))
    if n>0:
        ohdus.append(ohdu)
        ohduelists.append(gDirectory.Get(elistname))
        csvFields.append('overscan_y{0}'.format(ohdu))
        csvFields.append('overscan_yErr{0}'.format(ohdu))
        csvFields.append('overscan_x{0}'.format(ohdu))
        csvFields.append('overscan_xErr{0}'.format(ohdu))
        csvFields.append('rate{0}'.format(ohdu))
        csvFields.append('rateErr{0}'.format(ohdu))
        csvFields.append('minCharge{0}'.format(ohdu))
        csvFields.append('minChargeErr{0}'.format(ohdu))
        csvFields.append('midCharge{0}'.format(ohdu))
        csvFields.append('midChargeErr{0}'.format(ohdu))
        csvFields.append('maxCharge{0}'.format(ohdu))
        csvFields.append('maxChargeErr{0}'.format(ohdu))

csvwriter = csv.DictWriter(csvfile,fieldnames=csvFields)
#csvwriter.writeheader()
csvwriter.writer.writerow(csvwriter.fieldnames) #hack for python 2.6

arrPoisson = array.array('d')
arrPoissonErr = array.array('d')
arrActiveTime = array.array('d')
arrZero = array.array('d')

htypes = ["prescan","overscan_x","overscan_y","active"]
regioncuts = ["x>0 && x<{0} && y>0".format(prescanEdge),
        "x>{0} && y>0".format(activeEdgeX),
        "x>{0} && x<={1} && y>{2}".format(prescanEdge, activeEdgeX, activeEdgeY),
        "x>{0} && x<={1} && y>0 && y<={2} && y!={2}".format(prescanEdge, activeEdgeX, activeEdgeY)]
#exposures = [0.0, ((ccdCol/2.0)/nCol)*READOUT/nRow, READOUT*((ccdRow/2.0)/nRow), 0.5*READOUT*((ccdRow/2.0)/nRow)+EXPOSURE]

c.Clear()
c.Divide(len(ohdus),4)
histList = [] #not used, but needed so Python doesn't delete the histograms
noise = [0.1] * len(ohdus)
for iHdu in range(0,len(ohdus)):
    data.SetEventList(ohduelists[iHdu])
    for hnum in range(0,len(htypes)):
        c.cd(iHdu + hnum*len(ohdus) +1)
        gPad.SetLogy(1)
        hname = "h{0}{1}".format(ohdus[iHdu],htypes[hnum])
        data.Draw("ePix>>{0}({1}".format(hname, pixbinning),regioncuts[hnum],"colz")
        h = gDirectory.Get(hname)
        h.GetXaxis().SetRangeUser(-1.0,5.0)
        histList.append(h)
        fitvals_thisfit = fitvals_default
        if (htypes[hnum]=="overscan_x" or (htypes[hnum]=="active") and noise[iHdu]==0.1):
            s = h.Fit("gaus","QSN","",-0.5,0.5);
            if (int(s)==0): #good fit
                noise[iHdu] = s.Parameter(2)
                fitvals_thisfit["noise"] = noise[iHdu]
                print("hdu {0} {1}: noise={2}".format(ohdus[iHdu], htypes[hnum], noise[iHdu]))
        s = skipper_utils.fitPeaksPoisson(fitfunc,h,fitvals_thisfit,fixgain=True,fixnoise=(noise[iHdu]!=0.1))
        if (int(s)==0): #good fit
            #arrPoisson.append(s.Parameter(4))
            #arrPoissonErr.append(s.Error(4))
            #arrActiveTime.append(exposures[hnum])
            #arrZero.append(0.0)
            print("hdu {0} {1}: mu={2} pm {3}".format(ohdus[iHdu], htypes[hnum], s.Parameter(4), s.Error(4)))
            latex.DrawLatex(0.3,0.8,"mu={0:.3} \pm {1:.3}".format(s.Parameter(4),s.Error(4)))
            if (htypes[hnum]=="overscan_y"):
                csvDict['overscan_y{0}'.format(ohdus[iHdu])] = s.Parameter(4)/(binRows*binCols)
                csvDict['overscan_yErr{0}'.format(ohdus[iHdu])] = s.Error(4)/(binRows*binCols)
            elif (htypes[hnum]=="overscan_x"):
                csvDict['overscan_x{0}'.format(ohdus[iHdu])] = s.Parameter(4)/(binCols)
                csvDict['overscan_xErr{0}'.format(ohdus[iHdu])] = s.Error(4)/(binCols)
    c.cd(iHdu + (len(htypes)-1)*len(ohdus) +1)
    latex.DrawLatex(0.4,0.025,"OHDU{0}".format(ohdus[iHdu]))


c.cd()
c.Print(outfilename+".pdf");
c.Clear()

#if not badFit:
    #print(arrPoisson)
    #print(arrPoissonErr)
    #print(arrActiveTime)
    #c.SetLogy(0)
    #gStyle.SetStatX(0.4)
    #thegraph = TGraphErrors(len(arrPoisson),arrActiveTime,arrPoisson,arrZero,arrPoissonErr)
    #thegraph.Draw("AP")
    #thegraph.Fit("pol1")
    #thegraph.SetMinimum(0.0)
    #thegraph.SetMarkerSize(10)
    #c.Print(outfilename+".pdf");
#c.SetLogy(1)
#gStyle.SetStatX(0.9)

pixcuts = [("any","ePix>0.6"),
        ("big","ePix>10"),
        ("neg","ePix<-0.7")]

gStyle.SetOptStat(0)
c.Clear()
c.Divide(len(ohdus),4)
for iHdu in range(0,len(ohdus)):
    data.SetEventList(ohduelists[iHdu])
        #data.SetAlias("ele","(ePix-{0})/{1}".format(fitvals[iHdu]["zero"],fitvals[iHdu]["gain"]))

        #gPad.SetLogy(1)
        #data.Draw("ele>>hdata(500,-1,4)","x>8 && y>0 && "+ohducut,"colz")
        #c.Print(outfilename+".pdf");
    c.cd(iHdu +1)
    hname = "h2d{0}_colz".format(ohdus[iHdu])
    binning = "({0},-0.5,{1}-0.5,{2},-0.5,{3}-0.5)".format(nCol/25, nCol, nRow/25, nRow)
    data.Draw("y:x>>"+hname+binning,"x>0 && y>0 && {0} && ePix>0.6".format(edgecut),"colz")

    for hnum in range(0,len(pixcuts)):
        c.cd(iHdu + (hnum+1)*len(ohdus) +1)
        hname = "h2d{0}_{1}".format(ohdus[iHdu],pixcuts[hnum][0])
        binning = "({0},-0.5,{0}-0.5,{1},-0.5,{1}-0.5)".format(nCol, nRow)
        print(data.Draw("y:x>>"+hname+binning,"x>0 && y>0 && {0} && {1}".format(edgecut, pixcuts[hnum][1]),""))
    c.cd(iHdu + 3*len(ohdus) +1)
    latex.DrawLatex(0.4,0.025,"OHDU{0}".format(ohdus[iHdu]))
c.cd()
c.Print(outfilename+".pdf");


c.SetLogy(1)
gStyle.SetOptStat(0)

maxActiveY = min(activeEdgeY,nRow)
pol1 = TF1("f1", "pol1", 0, maxActiveY)


hists = []
c.Clear()
c.Divide(len(ohdus),4)
for iHdu in range(0,len(ohdus)):
    data.SetEventList(ohduelists[iHdu])

    data.Draw(">>elist","x>0 && y>0 && {0} && {1}".format(edgecut, maskcut))
    data.SetEventList(gDirectory.Get("elist"))

    c.cd(iHdu + 0*len(ohdus) +1)
    data.Draw("y:x>>hdenom{0}xy({1},{2})".format(ohdus[iHdu], binningX, binningY),"","colz")
    hDenomXY = gDirectory.Get("hdenom{0}xy".format(ohdus[iHdu]))
    hists.append(hDenomXY)
    hDenomXY.Draw("colz")
    hDenomXY.GetYaxis().SetRangeUser(0,nRow)

    c.cd(iHdu + 1*len(ohdus) +1)
    data.Draw("y:x>>h{0}xy({1},{2})".format(ohdus[iHdu], binningX, binningY),"ePix>0.6","colz")
    hXY = gDirectory.Get("h{0}xy".format(ohdus[iHdu]))
    hists.append(hXY)
    hXY.Draw("colz")
    hXY.GetYaxis().SetRangeUser(0,nRow)

    effXY = TEfficiency(hXY,hDenomXY)
    effXY.SetName("eff{0}xy".format(ohdus[iHdu]))
    hists.append(effXY)
    effXY.Draw("colz")
    gPad.Update()
    effHistXY = effXY.GetPaintedHistogram()
    effHistXY.GetYaxis().SetRangeUser(0,nRow)

    c.cd(iHdu + 2*len(ohdus) +1)
    gPad.SetLogy(1)

    hsX = THStack("hsx{0}".format(ohdus[iHdu]),"")
    data.Draw("x>>hdenom{0}x({1})".format(ohdus[iHdu], binningX),activecutY,"colz")
    hDenomX = gDirectory.Get("hdenom{0}x".format(ohdus[iHdu]))
    hsX.Add(hDenomX)
    data.Draw("x>>h{0}x({1})".format(ohdus[iHdu], binningX),"ePix>0.6 && {0}".format(activecutY),"colz")
    hX = gDirectory.Get("h{0}x".format(ohdus[iHdu]))
    hX.SetLineColor(2)
    hsX.Add(hX)

    hsX.Draw("nostack")
    hists.append(hsX)

    c.cd(iHdu + 3*len(ohdus) +1)
    effX = TEfficiency(hX,hDenomX)
    effX.SetName("eff{0}x".format(ohdus[iHdu]))
    hists.append(effX)
    effX.Draw()
    gPad.Update()
    effGraphX = effX.GetPaintedGraph()
    #maxEff = TMath.MaxElement(effGraphX.GetN(), effGraphX.GetY())
    #effGraphX.GetYaxis().SetRangeUser(0,1.5*maxEff)
    gStyle.SetStatX(0.4)
    effX.Write()

    c.cd(iHdu + 3*len(ohdus) +1)
    latex.DrawLatex(0.4,0.025,"OHDU{0}".format(ohdus[iHdu]))

c.cd()
c.Print(outfilename+".pdf");

hists = []
c.Clear()
c.Divide(len(ohdus),4)
for iHdu in range(0,len(ohdus)):
    data.SetEventList(ohduelists[iHdu])

    data.Draw(">>elist","x>0 && y>0 && {0} && {1}".format(edgecut, maskcut))
    data.SetEventList(gDirectory.Get("elist"))

    c.cd(iHdu + 0*len(ohdus) +1)
    gPad.SetLogz(1)
    data.Draw("ePix:y>>hpix{0}y({1},{2})".format(ohdus[iHdu], binningY, pixbinning),activecutX,"colz")
    hpixY = gDirectory.Get("hpix{0}y".format(ohdus[iHdu]))
    hists.append(hpixY)

    c.cd(iHdu + 1*len(ohdus) +1)
    fitfunc.SetParameter(0,1000)
    fitfunc.FixParameter(1,0)
    fitfunc.FixParameter(2,1)
    fitfunc.FixParameter(3,noise[iHdu])
    #fitfunc.SetParameter(4,abs(hpixY.GetMean(2)))
    fitfunc.SetParameter(4,1.0)
    fitfunc.SetParLimits(4,0.0,2*pixmax)
    #hpixY.FitSlicesY(fitfunc,0,-1,0,"QNL")

    xArr = array.array('d')
    yArr = array.array('d')
    yErrArr = array.array('d')
    for ixbin in range(1,hpixY.GetNbinsX()):
        proj = hpixY.ProjectionY("temp",ixbin,ixbin)
        if proj.Integral()<10: continue
        #if proj.GetMean()>0.9*pixmax: continue
        fitfunc.SetParameter(0,proj.GetEntries()*proj.GetXaxis().GetBinWidth(1))
        #fitfunc.SetParameter(4,0.01)
        fitfunc.SetParameter(4,abs(proj.GetMean()))
        #s = proj.Fit(fitfunc,"QNS")
        s = proj.Fit(fitfunc,"QNSL","",-0.5,pixmax)
        if int(s)!=0: continue
        #print(s.MinFcnValue())
        xArr.append(hpixY.GetXaxis().GetBinCenter(ixbin))
        yArr.append(fitfunc.GetParameter(4))
        yErrArr.append(fitfunc.GetParError(4))
        #print(ixbin,proj.Integral(),yArr[-1],yErrArr[-1])
        #maxY = max(maxY, fitfunc.GetParameter(4)+fitfunc.GetParError(4))
        maxY = max(maxY, fitfunc.GetParameter(4))
    
    #change the axis range after fitting, since FitSlicesY only fits data inside the axis limits
    hpixY.GetXaxis().SetRangeUser(0,nRow)
    hpixY.GetYaxis().SetRangeUser(-1.0,5.0)

    if len(xArr)>2:
        graph=TGraphErrors(len(xArr),xArr,yArr,nullptr,yErrArr)
        graph.Write("hgraph_{0}".format(iHdu))
        hists.append(graph)
        graph.Draw("A*")
        graph.GetYaxis().SetRangeUser(-0.1*maxY,1.5*maxY)


        #start with the first 3 points, and add more points as long as the fit chi2 remains acceptable
        lastgoodpoint = -1
        for lastpoint in range(2, len(xArr)):
            s = graph.Fit(pol1,"QNS","",xArr[0]-0.5,xArr[lastpoint]+0.5)
            if s.Prob()>1e-2:
                lastgoodpoint = lastpoint
            else:
                break

        if lastgoodpoint>0:
            s = graph.Fit(pol1,"QNS","",xArr[0]-0.5,xArr[lastgoodpoint]+0.5)
        #print(s.Chi2())
#
#    #gPad.SetLogy(1)
#    hmuY = gDirectory.Get("hpix{0}y_4".format(ohdus[iHdu]))
#    hists.append(hmuY)
#    if hmuY.GetEntries()>=5:
#        hmuY.Draw()
#        hmuY.GetXaxis().SetRangeUser(0,topbinedgeY+1.0)
#        pol1.SetParameters(0,0)
#        s = hmuY.Fit(pol1,"QNSR")
        #if int(s)==0:
            pol1.DrawCopy("same")
            minCharge = pol1.GetParameter(0)/(binRows*binCols)
            minChargeErr = pol1.GetParError(0)/(binRows*binCols)
            coeffVec = TVectorD(2)
            coeffVec[0] = 1.0
            coeffVec[1] = maxActiveY
            maxCharge = pol1.Eval(maxActiveY)/(binRows*binCols)
            maxChargeErr = math.sqrt(s.GetCovarianceMatrix().Similarity(coeffVec))/(binRows*binCols)
            coeffVec[1] = maxActiveY/2.0
            midCharge = pol1.Eval(maxActiveY/2.0)/(binRows*binCols)
            midChargeErr = math.sqrt(s.GetCovarianceMatrix().Similarity(coeffVec))/(binRows*binCols)
            rate = pol1.GetParameter(1)/daysPerRow/(binRows*binCols)
            rateErr = pol1.GetParError(1)/daysPerRow/(binRows*binCols)
            print("hdu {0}: mu(y=0)={1} pm {2} rate={3} pm {4}".format(ohdus[iHdu], minCharge, minChargeErr, rate, rateErr))
            latex.DrawLatex(0.15,0.7,"#splitline{{#mu(0)={0:.3}#pm{1:.3}}}{{rate={2:.3}#pm{3:.3}/day}}".format(minCharge, minChargeErr, rate, rateErr))
            csvDict['rate{0}'.format(ohdus[iHdu])] = rate
            csvDict['rateErr{0}'.format(ohdus[iHdu])] = rateErr
            csvDict['minCharge{0}'.format(ohdus[iHdu])] = minCharge
            csvDict['minChargeErr{0}'.format(ohdus[iHdu])] = minChargeErr
            csvDict['midCharge{0}'.format(ohdus[iHdu])] = midCharge
            csvDict['midChargeErr{0}'.format(ohdus[iHdu])] = midChargeErr
            csvDict['maxCharge{0}'.format(ohdus[iHdu])] = maxCharge
            csvDict['maxChargeErr{0}'.format(ohdus[iHdu])] = maxChargeErr
            #gStyle.SetStatX(0.4)

    c.cd(iHdu + 2*len(ohdus) +1)
    gPad.SetLogy(1)

    hsY = THStack("hsy{0}".format(ohdus[iHdu]),"")
    data.Draw("y>>hdenom{0}y({1})".format(ohdus[iHdu], binningY),activecutX,"colz")
    hDenomY = gDirectory.Get("hdenom{0}y".format(ohdus[iHdu]))
    hsY.Add(hDenomY)
    data.Draw("y>>h{0}y({1})".format(ohdus[iHdu], binningY),"ePix>0.6 && "+activecutX,"colz")
    hY = gDirectory.Get("h{0}y".format(ohdus[iHdu]))
    hY.SetLineColor(2)
    hsY.Add(hY)
    hists.append(hsY)
    hsY.Draw("nostack")
    hsY.GetXaxis().SetRangeUser(0,topbinedgeY)

    c.cd(iHdu + 3*len(ohdus) +1)
    effY = TEfficiency(hY,hDenomY)
    effY.SetName("eff{0}y".format(ohdus[iHdu]))
    hists.append(effY)
    #effY.Draw()
    #gPad.Update()
    effGraphY = effY.GetPaintedGraph()
    #maxEff = TMath.MaxElement(effGraphY.GetN(), effGraphY.GetY())
    #effGraphY.GetYaxis().SetRangeUser(0,1.5*maxEff)
    #effGraphY.GetXaxis().SetRangeUser(0,topbinedgeY)
    pol1.SetParameters(0,0)
    #s = effY.Fit(pol1,"QSR")
    effY.Write()
    effY.Draw()
    #if int(s)==0:
    #    minCharge = pol1.GetParameter(0)
    #    minChargeErr = pol1.GetParError(0)
    #    rate = pol1.GetParameter(1)/daysPerRow/(binRows*binCols)
    #    rateErr = pol1.GetParError(1)/daysPerRow/(binRows*binCols)
    #    #print "hdu {0}: minMu={1} pm {2} rate={3} pm {4}".format(ohdus[iHdu], minCharge, minChargeErr, rate, rateErr)
    #    #latex.DrawLatex(0.15,0.7,"#splitline{{min #mu={0:.3}#pm{1:.3}}}{{rate={2:.3}#pm{3:.3}/day}}".format(minCharge, minChargeErr, rate, rateErr))
    #    #gStyle.SetStatX(0.4)

    c.cd(iHdu + 3*len(ohdus) +1)
    latex.DrawLatex(0.4,0.025,"OHDU{0}".format(ohdus[iHdu]))

c.cd()
c.Print(outfilename+".pdf");

for key in csvFields:
    if key in csvDict:
        print(key, csvDict[key])
csvwriter.writerow(csvDict)

c.Print(outfilename+".pdf]");
outfile.Write()
outfile.Close()
csvfile.close()