diff --git a/calibrationSuite/pixelAnalysis.py b/calibrationSuite/pixelAnalysis.py new file mode 100644 index 0000000..2cef6fd --- /dev/null +++ b/calibrationSuite/pixelAnalysis.py @@ -0,0 +1,90 @@ +import numpy as np +import logging +import matplotlib.pyplot as plt +from calibrationSuite.fitFunctions import * +from calibrationSuite.ancillaryMethods import * +from scipy.optimize import curve_fit + + +logger = logging.getLogger(__name__) + +# add more analysis here... + +# fitIndex=3 +def analysis_one(clusters, nBins, sliceCoordinates, rows, cols, fitInfo, lowEnergyCut, highEnergyCut, fileInfo): + fitInfo = np.zeros((rows, cols, 5)) ## mean, std, area, mu, sigma + for i in range(rows): + for j in range(cols): + ax = plt.subplot() + + detRow, detCol = sliceToDetector(i, j, sliceCoordinates) + currGoodClusters = goodClusters(clusters, i, j, nPixelCut=4, isSquare=1) + if len(currGoodClusters) < 5: + print("too few clusters in slice pixel %d, %d: %d" % (i, j, len(currGoodClusters))) + continue + + energies = getClusterEnergies(currGoodClusters) + photonEcut = np.bitwise_and(energies > lowEnergyCut, energies < highEnergyCut) + nPixelClusters = (photonEcut > 0).sum() + print("pixel %d,%d has about %d photons" % (i, j, nPixelClusters)) + logger.info("pixel %d,%d has %d photons" % (i, j, nPixelClusters)) + + photonRegion = energies[photonEcut] + mean = photonRegion.mean() + std = photonRegion.std() + a, mu, sigma = histogramAndFitGaussian(ax, energies, nBins) + area = gaussianArea(a, sigma) + + ax.set_xlabel("energy (keV)") + ax.set_title("pixel %d,%d, small cluster cuts" % (detRow, detCol)) + plt.figtext(0.7, 0.8, "%d entries (peak)" % (area)) + plt.figtext(0.7, 0.75, "mu %0.2f" % (mu)) + plt.figtext(0.7, 0.7, "sigma %0.2f" % (sigma)) + fileNamePlot = "%s/%s_r%d_c%d_r%d_c%d_%s_E.png" % ( + fileInfo.outputDir, + fileInfo.className, + fileInfo.run, + fileInfo.camera, + detRow, + detCol, + fileInfo.label, + ) + logger.info("Writing plot: " + fileNamePlot) + plt.savefig(fileNamePlot) + plt.close() + + fileNameNpy = "%s/%s_r%d_c%d_r%d_c%d_%s_fitInfo.npy" % ( + fileInfo.outputDir, + fileInfo.className, + fileInfo.run, + fileInfo.camera, + detRow, + detCol, + fileInfo.label, + ) + logger.info("Writing npy: " + fileNameNpy) + np.save(fileNameNpy, fitInfo) + + fitInfo[i, j] = mean, std, area, mu, sigma + return fitInfo + + +# Helpers +def histogramAndFitGaussian(ax, energies, nBins): + y, bin_edges, _ = ax.hist(energies, nBins) + bins = (bin_edges[:-1] + bin_edges[1:]) / 2 + ##print(y, bins) + a, mean, std = estimateGaussianParametersFromUnbinnedArray(energies) + try: + popt, pcov = curve_fit(gaussian, bins, y, [a, mean, std]) + popt[1] + popt[2] + fittedFunc = gaussian(bins, *popt) + ax.plot(bins, fittedFunc, color="b") + return popt + except Exception: + return 0, 0, 0 + + +def sliceToDetector(sliceRow, sliceCol, sliceCoordinates): + return sliceRow + sliceCoordinates[0][0], sliceCol + sliceCoordinates[1][0] \ No newline at end of file diff --git a/suite_scripts/AnalyzeH5.py b/suite_scripts/AnalyzeH5.py index 5d756fb..d264c8c 100644 --- a/suite_scripts/AnalyzeH5.py +++ b/suite_scripts/AnalyzeH5.py @@ -1,75 +1,102 @@ import h5py +import argparse +import logging import numpy as np -import calibrationSuite.fitFunctions as fitFunctions -import calibrationSuite.ancillaryMethods as ancillaryMethods - -##import seaborn as sns import matplotlib.pyplot as plt from matplotlib.ticker import AutoMinorLocator -##import sys -import argparse +import calibrationSuite.fitFunctions as fitFunctions +import calibrationSuite.ancillaryMethods as ancillaryMethods +from calibrationSuite.pixelAnalysis import * from calibrationSuite.argumentParser import ArgumentParser -import logging - # for logging from current file logger = logging.getLogger(__name__) - import calibrationSuite.loggingSetup as ls import os - # log to file named .log currFileName = os.path.basename(__file__) ls.setupScriptLogging(currFileName[:-3] + ".log", logging.ERROR) # change to logging.INFO for full logging output +class FileNamingInfo: + def __init__(self, outputDir, className, run, camera, label): + self.outputDir = outputDir + self.className = className + self.run = run + self.camera = camera + self.label = label + + +# Setup logging. +# Log file gets appended to each new run, can manually delete for fresh log. +# Could change so makes new unique log each run or overwrites existing log. +logging.basicConfig( + filename="analyze_h5.log", + level=logging.INFO, # For full logging set to INFO which includes ERROR logging too + format="%(asctime)s - %(levelname)s - %(message)s", # levelname is log severity (ERROR, INFO, etc) +) + + +# Analysis of Hierarchical Data Format files (.h5 files) class AnalyzeH5(object): def __init__(self): - print("in init") - args = ArgumentParser().parse_args() - - self.run = args.run - self.files = args.files.replace(" ", "") - print(self.files) + self.outputDir = args.path logging.info("Output dir: " + self.outputDir) + self.label = args.label self.camera = 0 + self.run = args.run - def getFiles(self): - fileNames = self.files.split(",") - self.h5Files = [] - for f in fileNames: - print(f) - self.h5Files.append(h5py.File(f)) + if args.files is None: + print("No input files specified, quitting...") + logging.info("No input files specified, quitting...") + exit(1) + self.files = args.files.replace(" ", "") - def identifyAnalysis(self): - try: - self.analysisType = self.h5Files[0]["analysisType"] - self.sliceCoordinates = self.h5Files[0]["sliceCoordinates"][()] - except: - ## do something useful here, maybe - self.analysisType = None - ## but for now - self.analysisType = "cluster" - self.sliceCoordinates = [[270, 288], [59, 107]] + self.lowEnergyCut = 4 # fix - should be 0.5 photons or something + self.highEnergyCut = 30 # fix - should be 1.5 photons or something + self.sliceEdges = None + if args.slice_edges is not None: + self.sliceEdges = args.slice_edges.split(",") + self.sliceEdges = [int(curr) for curr in self.sliceEdges] + else: self.sliceEdges = [288 - 270, 107 - 59] + self.nBins = 100 + self.shiftEnergy = False if args.shift_energy_bits is None else True + self.analysisNum = 1 if args.analysis_mode is None else int(args.analysis_mode) - def sliceToDetector(self, sliceRow, sliceCol): - return sliceRow + self.sliceCoordinates[0][0], sliceCol + self.sliceCoordinates[1][0] + self.fileNameInfo = FileNamingInfo( + args.path, + self.__class__.__name__, + args.run, + 0, + args.label, + ) - def analyze(self): - if self.analysisType == "cluster": - self.clusterAnalysis() - else: + print("Output dir: " + self.fileNameInfo.outputDir) + logging.info("Output dir: " + self.fileNameInfo.outputDir) + + + def getFiles(self): + fileNames = self.files.split(",") + self.h5Files = [] + for currName in fileNames: + currFile = h5py.File(currName) + self.h5Files.append(currFile) + print("Using input file: " + currFile.filename) + logging.info("Using input file: " + currFile.filename) + + def clusterAnalysis(self): + + if self.analysisType != "cluster": # for now, until we have more types print("unknown analysis type %s" % (self.analysisType)) logging.info("unknown analysis type %s" % (self.analysisType)) + return - def clusterAnalysis(self): - clusters = None - energyHist = None + # energyHist = None clusters = np.concatenate([h5["clusterData"][()] for h5 in self.h5Files]) try: energyHist = np.concatenate(energyHist, h5["energyHistogram"][()]) @@ -79,10 +106,20 @@ def clusterAnalysis(self): self.nBins = 100 self.lowEnergyCut = 1 ## fix - should be 0.5 photons or something self.highEnergyCut = 10 ## fix - should be 1.5 photons or something - ##tmp - npyFileName = "%s/r%d_clusters.npy" % (self.outputDir, self.run) + + # concat never works here since h5 undefined + try: + # meant to do similar thing as clusters above? + pass # np.concatenate(energyHist, h5["energyHistogram"][()]) + # self.plotEnergyHist(energyHist, self.fileNameInfo) + except Exception as e: + print(f"An exception occurred: {e}") + logging.error(f"An exception occurred: {e}") + pass + + npyFileName = "%s/r%d_clusters.npy" % (self.fileNameInfo.outputDir, self.fileNameInfo.run) np.save(npyFileName, clusters) - logger.info("Wrote file: " + npyFileName) + logging.info("Wrote file: " + npyFileName) self.analyzeSimpleClusters(clusters) @@ -91,6 +128,12 @@ def clusterAnalysis(self): _, bins = np.histogram(energyHist, 250, [-5, 45]) plt.hist(bins[:-1], bins, weights=energyHist) ##, log=True) + + ''' comment-out for now ... + def plotEnergyHist(self, energyHist, fileInfo): + _, bins = np.histogram(energyHist, 250, [-5, 45]) + + plt.hist(bins[:-1], bins, weights=energyHist) # , log=True) plt.grid(which="major", linewidth=0.5) plt.title = "All pixel energies in run after common mode correction" plt.xlabel = "energy (keV)" @@ -115,84 +158,79 @@ def clusterAnalysis(self): np.save(npyFileName, energyHist) logger.info("Wrote file: " + npyFileName) plt.close() + ''' def analyzeSimpleClusters(self, clusters): + energy = clusters[:, :, 0] # .flatten() + if self.shiftEnergy: + energy *= 2 # temporary, due to bit shift + rows = self.sliceEdges[0] + cols = self.sliceEdges[1] + fitInfo = np.zeros((rows, cols, 4)) # mean, std, mu, sigma + fitIndex = 0 + + self.plot_overall_energy_distribution(energy, self.fileNameInfo) + + print("Analysis Mode: " + str(self.analysisNum)) + logging.info("Analysis Mode: " + str(self.analysisNum)) + if self.analysisNum == 1: + fitIndex = 3 + fitInfo = analysis_one( + clusters, + self.nBins, + self.sliceCoordinates, + rows, + cols, + fitInfo, + self.lowEnergyCut, + self.highEnergyCut, + self.fileNameInfo, + ) + #else if self.analysisNum == 2: + # ... + + self.save_fit_information(fitInfo, rows, cols, self.fileNameInfo) + self.plot_gain_distribution(fitInfo, self.fileNameInfo, fitIndex) + + def plot_overall_energy_distribution(self, energy, fileInfo): ax = plt.subplot() - energy = clusters[:, :, 0] ##.flatten() - ##energy *= 2 ## temporary, due to bit shift - ##print(energy[energy>0][666:777]) - print("mean energy above 0:" + str(energy[energy > 0].mean())) - logger.info("mean energy above 0:" + str(energy[energy > 0].mean())) - foo = ax.hist(energy[energy > 0], 100) + # print(energy[energy>0][666:777]) + print("mean energy above 0:", energy[energy > 0].mean()) + logging.info("mean energy above 0:" + str(energy[energy > 0].mean())) + + ax.hist(energy[energy > 0], 100) # 100 bins plt.xlabel = "energy (keV)" plt.title = "All pixels" - figFileName = "%s/%s_r%d_c%d_%s_E.png" % ( - self.outputDir, - self.__class__.__name__, - self.run, - self.camera, - self.label, + + fileName = "%s/%s_r%d_c%d_%s_E.png" % ( + fileInfo.outputDir, + fileInfo.className, + fileInfo.run, + fileInfo.camera, + fileInfo.label, ) - plt.savefig(figFileName) - logger.info("Wrote file: " + figFileName) + logging.info("Writing plot: " + fileName) + plt.savefig(fileName) plt.close() - rows = self.sliceEdges[0] - cols = self.sliceEdges[1] - fitInfo = np.zeros((rows, cols, 5)) ## mean, std, area, mu, sigma - for i in range(rows): - for j in range(cols): - detRow, detCol = self.sliceToDetector(i, j) - ax = plt.subplot() - currGoodClusters = ancillaryMethods.goodClusters(clusters, i, j, nPixelCut=4, isSquare=1) - if len(currGoodClusters) < 5: - print("too few clusters in slice pixel %d, %d: %d" % (i, j, len(currGoodClusters))) - logger.info("too few clusters in slice pixel %d, %d: %d" % (i, j, len(currGoodClusters))) - continue - energies = ancillaryMethods.getClusterEnergies(currGoodClusters) - photonEcut = np.bitwise_and(energies > self.lowEnergyCut, energies < self.highEnergyCut) - nPixelClusters = (photonEcut > 0).sum() - print("pixel %d,%d has about %d photons" % (i, j, nPixelClusters)) - logger.info("pixel %d,%d has about %d photons" % (i, j, nPixelClusters)) - photonRegion = energies[photonEcut] - mean = photonRegion.mean() - std = photonRegion.std() - a, mu, sigma = self.histogramAndFitGaussian(ax, energies, self.nBins) - area = fitFunctions.gaussianArea(a, sigma) - ax.set_xlabel("energy (keV)") - ax.set_title("pixel %d,%d, small cluster cuts" % (detRow, detCol)) - plt.figtext(0.7, 0.8, "%d entries (peak)" % (area)) - plt.figtext(0.7, 0.75, "mu %0.2f" % (mu)) - plt.figtext(0.7, 0.7, "sigma %0.2f" % (sigma)) - figFileName = "%s/%s_r%d_c%d_r%d_c%d_%s_E.png" % ( - self.outputDir, - self.__class__.__name__, - self.run, - self.camera, - detRow, - detCol, - self.label, - ) - plt.savefig(figFileName) - logger.info("Wrote file: " + figFileName) - plt.close() - - npyFileName = "%s/%s_r%d_c%d_r%d_c%d_%s_fitInfo.npy" % ( - self.outputDir, - self.__class__.__name__, - self.run, - self.camera, - detRow, - detCol, - self.label, - ) - np.save(npyFileName, fitInfo) - logger.info("Wrote file: " + npyFileName) - fitInfo[i, j] = mean, std, area, mu, sigma - - gains = fitInfo[:, :, 3] - goodGains = gains[np.bitwise_and(gains > 0, gains < 15)] + def save_fit_information(self, fitInfo, rows, cols, fileInfo): + fileName = "%s/%s_r%d_c%d_r%d_c%d_%s_fitInfo.npy" % ( + fileInfo.outputDir, + fileInfo.className, + fileInfo.run, + fileInfo.camera, + rows - 1, + cols - 1, + fileInfo.label, + ) + logging.info("Writing npy: " + fileName) + np.save(fileName, fitInfo) + + def plot_gain_distribution(self, fitInfo, fileInfo, fitIndex): + gains = fitInfo[:, :, fitIndex] + goodGains = gains[np.bitwise_and(gains > 0, gains < 30)] + ax = plt.subplot() ax.hist(goodGains, 100) ax.set_xlabel("energy (keV)") @@ -206,24 +244,24 @@ def analyzeSimpleClusters(self, clusters): ) plt.savefig(figFileName) - def histogramAndFitGaussian(self, ax, energies, nBins): - y, bin_edges, _ = ax.hist(energies, nBins) - bins = (bin_edges[:-1] + bin_edges[1:]) / 2 - ##print(y, bins) - a, mean, std = fitFunctions.estimateGaussianParametersFromUnbinnedArray(energies) - try: - popt, pcov = fitFunctions.curve_fit(fitFunctions.gaussian, bins, y, [a, mean, std]) - mu = popt[1] - sigma = popt[2] - fittedFunc = fitFunctions.gaussian(bins, *popt) - ax.plot(bins, fittedFunc, color="b") - return popt - except: - return 0, 0, 0 - if __name__ == "__main__": + print("Starting new run!") + logging.info("Starting new run!") ah5 = AnalyzeH5() ah5.getFiles() - ah5.identifyAnalysis() - ah5.analyze() + + # identify type of analysis + + # script fails earlier if not >= 1 h5 file + if "analysisType" in self.h5Files[0]: + self.analysisType = self.h5Files[0]["analysisType"] + # '[()]' gets us the data and not a reference + self.sliceCoordinates = self.h5Files[0]["sliceCoordinates"][()] + else: + # do something useful here, maybe + # but for now + self.analysisType = "cluster" + self.sliceCoordinates = [[270, 288], [59, 107]] + + ah5.clusterAnalysis() \ No newline at end of file