test/plottingTools/systematicErrors.py

0001 import numpy as np
0002 import ROOT
0003
0004 DIR_BOTH = 0
0005 DIR_UP = 1
0006 DIR_DOWN = -1
0007
0008 NUM_SECTORS = 68
0009 NUM_SECTORS_Y = 14
0010 # systematics has:
0011         #   a dict with with coordinate names "X", "Y" as keys
0012         #       - each value of these keys is a list/an array of systematic errors for each sector
0013         #       - so the list has the length of the number of sectors for that coordinate
0014         #       - these errors are quadratically added
0015         #   direction which can be DIR_DOWN, DIR_BOTH, or DIR_UP, depending on whether it adds only down, symmetric or up
0016         #   isRelative which is a flag telling whether the error is relative to the APE value or absolute
0017
0018 class SystematicErrors:
0019     def __init__(self):
0020         self.X = np.zeros(NUM_SECTORS)
0021         self.Y = np.zeros(NUM_SECTORS)
0022         # is not made seperately for X and Y. If this is wanted, make two separate objects
0023         self.isRelative = np.zeros(NUM_SECTORS, dtype=int)
0024         self.direction  = np.empty(NUM_SECTORS, dtype=int)
0025         self.direction.fill(DIR_BOTH) # just so it is clear that
0026
0027     def __getitem__(self, key):
0028         return getattr(self, key)
0029
0030     def getXFromList(self, X, startat=0):
0031         for i,x in enumerate(X):
0032             self.X[i+startat] = x
0033
0034     def getYFromList(self, Y, startat=0):
0035         for i,y in enumerate(Y):
0036             self.Y[i+startat] = y
0037
0038     # each line has the structure: xerr yerr isrel direction
0039     def write(self, fileName):
0040         with open(fileName, "w") as fi:
0041             for x, y, rel, direc in zip(self.X, self.X, self.isRelative, self.direction):
0042                     fi.write("{} {} {} {}".format(x, y, rel, direc))
0043
0044     def read(self, fileName):
0045         with open(fileName, "r") as fi:
0046             sector = 0
0047             for line in fi:
0048                 x, y, rel, direc = line.rstrip().split(" ")
0049                 self.X[sector] = float(x)
0050                 self.Y[sector] = float(y)
0051                 self.isRelative[sector] = int(rel)
0052                 self.direction[sector] = int(direc)
0053                 sector += 1
0054         return self
0055
0056 # difference between ape values in each sector
0057 # returns a SystematicErrors object with values
0058 def apeDifference(minuend, subtrahend):
0059     fileA = ROOT.TFile(minuend, "READ")
0060     fileB = ROOT.TFile(subtrahend, "READ")
0061     apeTreeA_X = fileA.Get("iterTreeX")
0062     apeTreeA_X.SetDirectory(0)
0063     apeTreeB_X = fileB.Get("iterTreeX")
0064     apeTreeB_X.SetDirectory(0)
0065     apeTreeA_Y = fileA.Get("iterTreeY")
0066     apeTreeA_Y.SetDirectory(0)
0067     apeTreeB_Y = fileB.Get("iterTreeY")
0068     apeTreeB_Y.SetDirectory(0)
0069
0070     fileA.Close()
0071     fileB.Close()
0072
0073     # get to last iteration of each tree
0074     apeTreeA_X.GetEntry(apeTreeA_X.GetEntries()-1)
0075     apeTreeB_X.GetEntry(apeTreeB_X.GetEntries()-1)
0076     apeTreeA_Y.GetEntry(apeTreeA_Y.GetEntries()-1)
0077     apeTreeB_Y.GetEntry(apeTreeB_Y.GetEntries()-1)
0078
0079     difference = SystematicErrors()
0080     isRel = 0
0081     direc = 0
0082
0083     for sector in range(1, NUM_SECTORS+1):
0084         name = "Ape_Sector_{}".format(sector)
0085
0086         diffX = abs(getattr(apeTreeA_X, name) - getattr(apeTreeB_X, name))
0087         difference.X[sector-1] = diffX
0088         if sector <= NUM_SECTORS_Y:
0089             diffY = abs(getattr(apeTreeA_Y, name) - getattr(apeTreeB_Y, name))
0090             difference.Y[sector-1] = diffY
0091         difference.isRel[sector-1] = isRel
0092         difference.direction[sector-1] = direc
0093
0094     return difference
0095
0096
0097 # inFile is allData.root, not allData_iterationApe.root
0098 # returns two arrays with values in x and y
0099 def numberOfHits(inFileName):
0100     inFile = ROOT.TFile(inFileName, "READ")
0101     num_x = np.zeros(NUM_SECTORS, dtype=int)
0102     num_y = np.zeros(NUM_SECTORS, dtype=int)
0103     for sector in range(1, NUM_SECTORS+1):
0104         xhist = inFile.Get("ApeEstimator1/Sector_{}/Results/h_ResX".format(sector))
0105         num_x[sector-1] = xhist.GetEntries()
0106         if sector <= NUM_SECTORS_Y:
0107             yhist = inFile.Get("ApeEstimator1/Sector_{}/Results/h_ResY".format(sector))
0108             num_y[sector-1] = yhist.GetEntries()
0109     inFile.Close()
0110     return num_x, num_y
0111
0112 def main():
0113     pass
0114
0115 if __name__ == "__main__":
0116     main()