makeBrightnessMap.py

Go to the documentation of this file.

#!/usr/bin/env python
import ROOT
import numpy as np
import argparse
import bisect
import math

def findClosestVal(val, quantList):
    n = len( quantList )
    i = bisect.bisect( quantList, val )
    if i == 0:
        return 0
    if i == n:
        return n-1
    if math.fabs( val - quantList[i-1] ) < math.fabs( val - quantList[i] ):
        return i-1
    else:
        return i

if __name__=='__main__':

    parser = argparse.ArgumentParser(description='Process fiber brightness files')
    parser.add_argument('-fd', '--far-detector', action='store_true')
    parser.add_argument('-nd', '--near-detector', action='store_true')
    args = parser.parse_args()

    if args.far_detector and args.near_detector:
        print "Can only set --fd or --nd!"
        exit(1)

    if not args.far_detector and not args.near_detector:
        print "Must set either --fd or --nd!"
        exit(1)

    if args.far_detector:
        fIn = open("FDModuleData.txt", "r")
        fOutName = "fdBrightness.root"
    if args.near_detector:
        fIn = open("NDModuleData.txt", "r")
        fOutName = "ndBrightness.root"

    planeDict = {}
    plane_moduleDict = {}
    brightnessLst = []
    grBright = ROOT.TGraph()

    maxModule = 0
    maxCell = 0
    maxPlane = 0

    maxBrightness = -9999
    minBrightness =  9999
    for line in fIn:
        splitLine = line.split()
        plane = int(splitLine[0])
        cell = int(splitLine[1])
        brightness = float(splitLine[2])
        module = cell/32

        if brightness > 20 and brightness < 180:
            brightnessLst.append(brightness/100.)

        if plane > maxPlane:
            maxPlane = plane
        if cell > maxCell:
            maxCell = cell
        if module > maxModule:
            maxModule = module
        if brightness < minBrightness:
            minBrightness = brightness
        if brightness > maxBrightness:
            maxBrightness = brightness

    print "min brightness:", minBrightness, "max brightness:", maxBrightness

    brightness1D = ROOT.TH1D("brightness1D", "Brightness Level;Brightness;Cells", 1200, -1, 5)
    brightnessValue = ROOT.TH1D("brightnessValue", "Brightness Level Bin Values;Brightness Bin Number;Values", 9, 0, 9)
    brightness2D = ROOT.TH2D("brightness2D", ";Plane;Cell", maxPlane+1, 0, maxPlane+1, maxCell+1, 0, maxCell+1)
    brightness2DX = ROOT.TH2D("brightness2DX", ";Plane;Cell", maxPlane/2+1, 0, maxPlane/2+1, maxCell+1, 0, maxCell+1)
    brightness2DY = ROOT.TH2D("brightness2DY", ";Plane;Cell", maxPlane/2+1, 0, maxPlane/2+1, maxCell+1, 0, maxCell+1)

    brightness2D_module = ROOT.TH2D("brightness2D_module", ";Plane;Module", maxPlane+1, 0, maxPlane+1, maxModule+1, 0, maxModule+1)
    brightness2D_bin = ROOT.TH2C("brightness2D_bin", "Brightness Level Bin Map;Plane;Cell", maxPlane+1, 0, maxPlane+1, maxCell+1, 0, maxCell+1)
    brightness_byplane = ROOT.TH1D("brightness_plane", ";Plane", maxPlane+1, 0, maxPlane+1);

    print "maxplane: {}".format(maxPlane)
    print "maxcell: {}".format(maxCell)

    brightnessLst.sort()
    brightnessArr = np.array(brightnessLst)

    #quantList = []
    #for i in [-2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5]:
    #    quantList.append( 100*0.5*(1.0 + ROOT.TMath.Erf(i/ROOT.TMath.Sqrt(2.0))) )

    quantList = [10, 20, 30, 40, 50, 60, 70, 80, 90]
    quantVals = np.percentile(brightnessArr, quantList)

    brightnessValue.SetBinContent(1, 1)
    for i in range(len(quantVals)):
        brightnessValue.SetBinContent(i+1, quantVals[i])

    for i in xrange(len(quantList)):
        print str(quantList[i])+"%:", quantVals[i]

    for iX in xrange(brightness2D.GetXaxis().GetNbins()):
        for iY in xrange(brightness2D.GetYaxis().GetNbins()):
            brightness2D.SetBinContent(iX+1, iY+1, 1.0)

    for iX in xrange(brightness2D_module.GetXaxis().GetNbins()):
        for iY in xrange(brightness2D_module.GetYaxis().GetNbins()):
            brightness2D_module.SetBinContent(iX+1, iY+1, 1.0)

    for iX in xrange(brightness2D_bin.GetXaxis().GetNbins()):
        for iY in xrange(brightness2D_bin.GetYaxis().GetNbins()):
            brightness2D_bin.SetBinContent(iX+1, iY+1, 1)

    position = fIn.seek(0, 0);
    for line in fIn:
        splitLine = line.split()
        plane = int(splitLine[0])
        cell = int(splitLine[1])
        brightness = float(splitLine[2])
        module = cell/32

        # Fill brightness 1D
        brightness1D.Fill(brightness/100.0)

        # Fill brightness 2D, 2DX, 2DY
        if brightness > 20 and brightness < 180:
            iquantBrightness = findClosestVal(brightness/100., quantVals)
            quantBrightness = quantVals[iquantBrightness]
            brightness2D.SetBinContent(plane+1, cell+1, quantBrightness)

            if plane%2 == 0:
                brightness2DX.SetBinContent((plane/2)+1, cell+1, quantBrightness)
            else:
                brightness2DY.SetBinContent((plane/2)+1, cell+1, quantBrightness)
            if plane in planeDict:
                planeDict[plane][0] += 1
                planeDict[plane][1] += brightness
            else:
                planeDict[plane] = [1, brightness]
            if (plane, module) in plane_moduleDict:
                plane_moduleDict[(plane,module)][0] += 1
                plane_moduleDict[(plane,module)][1] += brightness
            else:
                plane_moduleDict[(plane,module)] = [1, brightness]

    # Fill brightness by plane
    for key in planeDict:
        grBright.SetPoint(key, key, planeDict[key][1]/planeDict[key][0])
        brightness_byplane.SetBinContent(key+1, planeDict[key][1]/(100.0*planeDict[key][0]));

    # Fill brightness by module
    for key in plane_moduleDict:
        plane = key[0]
        module = key[1]
        count = float(plane_moduleDict[key][0])
        totBrightness = float(plane_moduleDict[key][1])/100.
        brightness2D_module.SetBinContent(plane+1, module+1, totBrightness/count)

    # Fill brightness by bin, set cell brightness level to module average if it's not in [20, 180]
    position = fIn.seek(0, 0);
    for line in fIn:
        splitLine = line.split()
        plane = int(splitLine[0])
        cell = int(splitLine[1])
        brightness = float(splitLine[2])
        module = cell/32

        if brightness > 20 and brightness < 180:
            iquantBrightness = findClosestVal(brightness/100., quantVals)
            brightness2D_bin.SetBinContent(plane+1, cell+1, iquantBrightness+1)
        else:
            iquantBrightness = findClosestVal(brightness_byplane.GetBinContent(plane+1, module+1), quantVals)
            brightness2D_bin.SetBinContent(plane+1, cell+1, iquantBrightness+1)

    fOut = ROOT.TFile(fOutName, "recreate")
    fOut.WriteTObject(brightness2D, "BrightnessByCell")
    fOut.WriteTObject(brightnessValue, "BrightnessValue")
    fOut.WriteTObject(brightness2D_module, "BrightnessByModule")
    fOut.WriteTObject(brightness_byplane, "BrightnessByPlane")
    fOut.WriteTObject(brightness2D_bin, "BrightnessByBin")