In This Package:

Tools.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# Python tools for working with calibration parameters
# 

__all__ = [ "PmtTableWriterAlg", "mergeTables", "parseTable" ]

# Load GaudiPython
from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl
import GaudiKernel.SystemOfUnits as units
import math

# Shortcuts to DayaBay classes
Detector = gbl.DayaBay.Detector
Site = gbl.Site
DetectorId = gbl.DetectorId
Context = gbl.Context
ServiceMode = gbl.ServiceMode
PmtCalibData = gbl.DayaBay.PmtCalibData
AdPmtSensor = gbl.DayaBay.AdPmtSensor

# Generate PMT calibration table
class PmtTableWriterAlg(DybPythonAlg):
    "Generate a calibration table from calibration run statistics"
    def __init__(self,name):
        DybPythonAlg.__init__(self,name)
        # Properties
        self.StatsPath = "/file0/pmtCalibLeadingEdge"
        self.OutputTableName = "pmtCalibTable.txt"
        self.IsDryAd = False
        self.SetDefaultCalib = True
        self.SourcePosition = [0.,0.,0.]
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.info("initializing")

        # Initialize services
        self.cableSvc = self.svc('ICableSvc','StaticCableSvc')
        if self.cableSvc == None:
            self.error("Failed to get StaticCableSvc")
            return FAILURE
        self.statsSvc = self.svc('IStatisticsSvc','StatisticsSvc')
        if self.statsSvc == None:
            self.error("Failed to get StatisticsSvc")
            return FAILURE
        self.pmtGeomSvc = self.svc('IPmtGeomInfoSvc','PmtGeomInfoSvc')
        if self.pmtGeomSvc == None:
            self.error("Failed to get PmtGeomInfoSvc")
            return FAILURE
        return status

    def execute(self):
        self.info("executing")
        return SUCCESS
        
    def finalize(self):
        self.info("finalizing")
        # Generate the PMT calibration table, and write to file
        # Step 1: initialize the table
        pmtCalibList = self.initializeCalibTable()
        self.info("I initialized the table") #tmp
        # Step 2: incorporate data using the current calibration statistics
        status = self.updateCalibration(pmtCalibList)
        if not status.isSuccess(): return status
        # Step 3: Print out the table
        self.writeCalibTable(pmtCalibList)
        status = DybPythonAlg.finalize(self)
        return status

    def initializeCalibTable(self):
        # Initialize table of PMT calibrations
        context = Context()
        context.SetSite(Site.kAll)
        context.SetDetId(DetectorId.kAll)
        svcMode = ServiceMode(context, 0)  # Dummy ServiceMode
        adPmts = self.cableSvc.adPmtSensors(svcMode)
        poolPmts = self.cableSvc.poolPmtSensors(svcMode)
        adPmtCalibList = []
        for adPmt in adPmts:
            pmtCalib = PmtCalibData()
            pmtCalib.m_pmtId = adPmt
            if self.SetDefaultCalib: self.initializeCalib(pmtCalib)
            adPmtCalibList.append( pmtCalib )
        poolPmtCalibList = []
        for poolPmt in poolPmts:
            pmtCalib = PmtCalibData()
            pmtCalib.m_pmtId = poolPmt
            if self.SetDefaultCalib: self.initializeCalib(pmtCalib)
            poolPmtCalibList.append( pmtCalib )
        pmtCalibList = adPmtCalibList + poolPmtCalibList
        return pmtCalibList

    def initializeCalib(self, pmtCalib):
        # Initialize using an approximate default calibration
        pmtCalib.m_status = PmtCalibData.kGood
        pmtCalib.m_speHigh = 20.0
        pmtCalib.m_sigmaSpeHigh = 0.0
        pmtCalib.m_speLow = 2.0
        pmtCalib.m_timeOffset = 0.0
        pmtCalib.m_timeSpread = 0.0
        pmtCalib.m_efficiency = 1.0
        pmtCalib.m_prePulseProb = 0.0
        pmtCalib.m_afterPulseProb = 0.0
        pmtCalib.m_darkRate = 0.0

    def updateCalibration(self, pmtCalibList):
        # Process the calibration statistics and update the
        # calibration parameters
        detDirs = self.statsSvc.getSubFolders(self.StatsPath)
        for detDir in detDirs:
            detector = Detector(detDir)
            if( detector.site() == Site.kUnknown
                or detector.detectorId() == DetectorId.kUnknown ):
                self.error("Unknown detector: "+detDir)
                return FAILURE
            else:
                self.info("Processing: "+detDir)
            # Find the PMT IDs that are in the current detector
            pmtCalibs = []
            for pmtCalib in pmtCalibList:
                pmtId = pmtCalib.m_pmtId
                if (pmtId.site() == detector.site()
                    and pmtId.detectorId() == detector.detectorId()):
                    pmtCalibs.append(pmtCalib)
            if detector.isAD():
                status = self.updateAdCalibration(detector, pmtCalibs)
                if not status.isSuccess(): return status
            elif detector.isWaterShield():
                status = self.updateWaterShieldCalibration(detector, pmtCalibs)
                if not status.isSuccess(): return status
        return SUCCESS

    def updateAdCalibration(self, detector, pmtCalibList):
        print "I am in updateAdCalibration" #tmp
        # Loop over the pmts and process calibration stats
        meanTimeOffset = 0.0
        nMeanTimeOffset = 0
        relativeEfficiency = {}
        for ring in range(1,9): relativeEfficiency[ring] = []
        for pmtCalib in pmtCalibList:

            pmtId = pmtCalib.m_pmtId
            if pmtId.ring() < 1: continue  # Skip 2-inch calib pmts
            # Get the calibration histograms for this pmt
            [occupancyH, tdcOffsetH, adcMeanH, adcSigmaH] = self.getAdCalibHists(pmtId)
            if occupancyH==None or tdcOffsetH==None or adcMeanH==None:
                self.error("Failed to get calibration histograms for "
                           +detector.detName())
                return FAILURE
            occupancy = occupancyH.GetBinContent(pmtId.column())
            occupancyUncert = occupancyH.GetBinError(pmtId.column())
            tdcOffset = tdcOffsetH.GetBinContent(pmtId.column())
            tdcOffsetUncert = tdcOffsetH.GetBinError(pmtId.column())
            adcMean       = adcMeanH.GetBinContent(pmtId.column())            
            adcMeanUncert = adcMeanH.GetBinError(pmtId.column())
            adcReso       = adcSigmaH.GetBinContent(pmtId.column())
            adcResoUncert = adcSigmaH.GetBinError(pmtId.column())
            # Update PMT status
            pmtCalib.m_status = PmtCalibData.kUnknown
            if( occupancy < 0.01 or occupancyUncert > 0.1 ):
                pmtCalib.m_status += PmtCalibData.kBadEfficiency
            if( tdcOffsetUncert > 20.0 ):
                pmtCalib.m_status += PmtCalibData.kBadTiming
            if( adcMean < 10.0 or adcMean > 50.0 or adcMeanUncert > 4.0 ):
                pmtCalib.m_status += PmtCalibData.kBadGain
            # If no problems found, PMT is good
            if( pmtCalib.m_status == PmtCalibData.kUnknown ):
                pmtCalib.m_status = PmtCalibData.kGood
            # Convert time offset to nanoseconds
            tdcToNs = 1.5625
            lightTravelTime = self.timeToAdPmt(self.SourcePosition, pmtId)
            timeOffset = -tdcOffset * tdcToNs - lightTravelTime
            timeSpread = tdcOffsetUncert*tdcToNs
            pmtCalib.m_timeOffset = timeOffset
            pmtCalib.m_timeSpread = timeSpread
            # wangzhe: comment out temporarily
            # pmtCalib.m_efficiency = -math.log(1 - occupancy)
            # wz
            pmtCalib.m_speHigh = adcMean
            pmtCalib.m_sigmaSpeHigh = adcReso
            # FIXME: need dedicated Low-gain calibration:
            # According to DocDB 5038, for dry run it is set to 19.5
            pmtCalib.m_speLow = adcMean / 19.5
            # Update mean time
            if not (pmtCalib.m_status & PmtCalibData.kBadTiming):
                meanTimeOffset += timeOffset
                nMeanTimeOffset += 1
            # Update mean efficiency by ring
            if not (pmtCalib.m_status & PmtCalibData.kBadEfficiency):
                relativeEfficiency[pmtId.ring()].append(pmtCalib.m_efficiency)
        # Adjust timing offset so mean over detector is zero
        if nMeanTimeOffset != 0:
            meanTimeOffset /= nMeanTimeOffset
        for pmtCalib in pmtCalibList:
            pmtId = pmtCalib.m_pmtId
            if pmtId.ring() < 1: continue  # Skip 2-inch calib pmts
            pmtCalib.m_timeOffset -= meanTimeOffset
        # Adjust the relative efficiency according to the mean for each ring
        meanEfficiency = {}
        for ring in relativeEfficiency.keys():
            relEff = relativeEfficiency[ring]
            meanEfficiency[ring] = 0
            if len(relEff) > 0:
                meanEfficiency[ring] = sum(relEff)/len(relEff)
        for pmtCalib in pmtCalibList:
            pmtId = pmtCalib.m_pmtId
            if pmtId.ring() < 1: continue  # Skip 2-inch calib pmts
            if meanEfficiency[pmtId.ring()] != 0: #jpochoa
                pmtCalib.m_efficiency /= meanEfficiency[pmtId.ring()]
            else:
                pmtCalib.m_efficiency = -1

        
                
        return SUCCESS

    def getAdCalibHists(self, pmtId):
        # Return the calibration histograms for this PMT
        path = self.StatsPath + "/" + pmtId.detName()
        path += "/ring_%d" % pmtId.ring()
        occupancyH = self.statsSvc.get(path + "/occupancy")
        tdcOffsetH = self.statsSvc.get(path + "/tdcOffset")
        if self.StatsPath=="/file0/pmtCalibLeadingEdge":
            adcMeanH =   self.statsSvc.get(path + "/adcMean")
        else:
            adcMeanH = self.statsSvc.get(path + "/adcMedian")
        adcSigmaH =  self.statsSvc.get(path + "/adcSigma")
        return [occupancyH, tdcOffsetH, adcMeanH, adcSigmaH]
    
    def getPmtDistance(self, position, pmtId):
        "Return the distance from position to the pmt with ID pmtID"
        # No Site.kSAB in current geometry tree. Change it to kDayaBay to work around
        site = pmtId.site()
        if site == Site.kSAB:
            site = Site.kDayaBay
        tmpPmtId = AdPmtSensor( pmtId.ring(), pmtId.column(), site, pmtId.detectorId() )
        pmtGeomInfo = self.pmtGeomSvc.get( tmpPmtId.fullPackedData() )
        Rpmt = pmtGeomInfo.localPosition()
        dR = [position[0] - Rpmt.x(),
              position[1] - Rpmt.y(),
              position[2] - Rpmt.z()]
        return math.sqrt(dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2])

    def timeToAdPmt(self, position, pmtId):
        "Calculate the light travel time from source to pmt"
        dr = self.getPmtDistance(position, pmtId)
        lightSpeedInVac = 299.79 * units.mm/units.ns
        liquidIndexOfRefrac = 1.45
        lightSpeed = lightSpeedInVac
        if not self.IsDryAd:
            lightSpeed /= liquidIndexOfRefrac
        return dr / lightSpeed


    def updateWaterShieldCalibration(self, detector, pmtCalibList):
        # Update the calib data for this water shield based on calibration data
        self.error("Calibration of water shield not implemented")
        self.error("I am ignoring this for now") #jpochoa
        return SUCCESS #jpochoa
            
    def calibTableLine(self, pmtCalib):
        # Write a line for the calibration table
        pmtLabel = None
        pmtId = pmtCalib.m_pmtId
        if pmtId.isAD():
            pmtLabel = "%s-ring%02d-column%02d" % (pmtId.detName(),
                                                   pmtId.ring(),
                                                   pmtId.column())
        elif pmtId.isWaterShield():
            inward = "in"
            if not pmtId.inwardFacing(): inward = "out"
            pmtLabel = "%s-wall%02d-spot%02d-%s" % (pmtId.detName(),
                                                    pmtId.wallNumber(),
                                                    pmtId.wallSpot(),
                                                    inward)
        line=("%d %30s %2d %7.3f %5.3f %5.3f %7.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n"
              % (pmtCalib.m_pmtId.fullPackedData(),
                 pmtLabel,
                 pmtCalib.m_status,
                 pmtCalib.m_speHigh,
                 pmtCalib.m_sigmaSpeHigh,
                 pmtCalib.m_speLow,
                 pmtCalib.m_timeOffset / units.ns,
                 pmtCalib.m_timeSpread / units.ns,
                 pmtCalib.m_efficiency,
                 pmtCalib.m_prePulseProb,
                 pmtCalib.m_afterPulseProb,
                 pmtCalib.m_darkRate / units.hertz
                 )
              )
        return line

    def writeCalibTable(self, pmtCalibList):
        # Write a calibration table using the list of pmt calibrations
        tableLines = []
        tableLines.append("# Table of Calibration constants for each channel\n")
        tableLines.append("# [pmtID] [description] [status] [speHigh] [sigmaSpe] [speLow] [timeOffset] [timeSpread] [efficiency] [prePulse] [afterPulse] [darkRate]\n")

        # sort them by id, otherwise it is out of control
        idCalib = {}
        for pmtCalib in pmtCalibList:
            idCalib[ pmtCalib.m_pmtId.fullPackedData() ] = pmtCalib

        for id in sorted( idCalib.keys() ):
            pmtCalib = idCalib[id]
            line = self.calibTableLine(pmtCalib)
            tableLines.append(line)
        outFile = open(self.OutputFileName,'w')
        outFile.writelines(tableLines)
        outFile.close()        


# Merge calibration tables
def mergeTables(primaryFile, changeMap, outputFile):
    "Simple Non-NuWa python function for merging calibration text tables"
    inputTable = parseTable(primaryFile)
    # Update table
    for changeFile in changeMap.keys():
        changeTable = parseTable(changeFile)
        changeFields = changeMap[changeFile]
        for fieldName in changeFields:
            fieldName = fieldName.strip()
            if not inputTable["columns"].has_key(fieldName):
                print "Field \"",fieldName,"\" doesn't exist in input table."
                return
            inputTable["columns"][fieldName] = changeTable["columns"][fieldName]
    # Write output table comments and header
    outputTableLines = []
    for comment in inputTable["comments"]:
        outputTableLines.append(comment)
    outputTableLines.append(inputTable["header"])
    # determine field lengths
    fieldLength = {}
    for fieldName in inputTable["fields"]:
        fieldLength[fieldName] = max( [len(data) for data in
                                       inputTable["columns"][fieldName] ] )
    # write data lines
    for rowIdx in range(inputTable["numberOfRows"]):
        line = ""
        for fieldName in inputTable["fields"]:
            format = "%"+str(fieldLength[fieldName])+"s " 
            line += format % inputTable["columns"][fieldName][rowIdx]
        outputTableLines.append(line+"\n")
    open(outputFile,"w").writelines(outputTableLines)

def parseTable(filename):
    "Simple text table parser"
    tableLines = open(filename).readlines()
    fieldNames = []
    comments = []
    headerLine = None
    columns = {}
    nRows = 0
    for line in tableLines:
        data = line.strip()
        # Catch comments / header
        if len(data)>0 and data[0]=='#':
            comment = data[1:].strip()
            if len(comment)>0 and comment[0] =='[':
                # Found header
                if len(fieldNames) > 0:
                    print "Headers already defined by line \"",line,"\""
                    return None
                tokens = comment.split(']')
                if len(tokens[-1].strip())==0: tokens = tokens[:-1]
                for token in tokens:
                    fieldStart = token.rfind('[')
                    if fieldStart<0:
                        print "Invalid token \"",token,"\" in line \"",line,"\""
                        return None
                    fieldName = token[fieldStart+1:].strip()
                    fieldNames.append(fieldName)
                    columns[fieldName] = []
                headerLine = line
            else:
                comments.append(line)
            continue
        # Data line
        fields = data.split()
        if len(fields)==0: continue  # Empty line
        if len(fields) != len(columns):
            print "Invalid number of fields in line \"",line,"\""
            return None
        for fieldIdx in range(len(fields)):
            columns[fieldNames[fieldIdx]].append( fields[fieldIdx] )
        nRows += 1
    return { "fields":fieldNames,
             "comments":comments,
             "header":headerLine,
             "columns":columns,
             "numberOfRows":nRows}

---