diff --git a/scripts/ancillaryMethods.py b/scripts/ancillaryMethods.py new file mode 100644 index 0000000..470c01a --- /dev/null +++ b/scripts/ancillaryMethods.py @@ -0,0 +1,43 @@ +import numpy as np +from scipy.stats import binned_statistic + +def makeProfile(x, y, bins, range=None, spread=False): + ## NaN for empty bins are suppressed + ## using mean root(N) for non-empty bins to calculate 0 var weights + ## + ## spread=True to return standard deviation instead of standard error + + meansObj = binned_statistic(x, [y, y**2], bins=bins, range=range, statistic='mean') + means, means2 = meansObj.statistic + countsObj = binned_statistic(x, [y, y**2], bins=bins, range=(0,1), statistic='count') + bin_N = countsObj.statistic + yErr = np.sqrt(means2 - means**2) + if not spread: + root_N = np.sqrt(bin_N) + root_N[root_N==0] = root_N[root_N>0].mean() + yErr = yErr/root_N + ##yErr = yErr.clip(0, 6666666.) + bin_edges = means_result.bin_edges + bin_centers = (bin_edges[:-1] + bin_edges[1:])/2. + usefulBins = bin_N>0 + return bin_centers[usefulBins], means[usefulBins], yErr[usefulBins] + +def plotProfile(x, y, yErr): + plt.errorbar(x=x, y=y, yerr=yErr, linestyle='none', marker='.') + +def selectedClusters(clusters, row, col, lowEnerygCut, highEnergyCut, nPixelCut=4, isSquare=1): + pass + +def goodClusters(clusters, row, col, nPixelCut=4, isSquare=None): + ##print(clusters) + pixelRowCol = np.bitwise_and((clusters[:,:,1] == row), + (clusters[:,:,2] == col)) + if isSquare is None: + small = clusters[:,:,3] lowEnergyCut + ) + pixelEcut = np.bitwise_and( + pixelEcut0, energy < highEnergyCut + ) + nPixelClusters = (pixelEcut > 0).sum() + + mean = std = mu = sigma = 0 + + # Select energy vals that passed cut conditions + # (selects elements from energy where corresponding element in pixelEcut is True) + pixelE = energy[pixelEcut > 0] + + if nPixelClusters > 5: # only do analysis if enough pixels + print("pixel %d,%d has %d photons" % (i, j, nPixelClusters)) + logger.info("pixel %d,%d has %d photons" % (i, j, nPixelClusters)) + ax = plt.subplot() + y, bin_edges, _ = ax.hist(pixelE, 100) + bins = (bin_edges[:-1] + bin_edges[1:]) / 2 + # print(y, bins) + mean, std = fitFunctions.estimateGaussianParameters(pixelE) + try: + # Set maxfev arg > 800?? (fails to find optimal params for some clusters) + popt, pcov = curve_fit(fitFunctions.gaussian, bins, y, [3, mean, std]) + + mu = popt[1] + sigma = popt[2] + fitInfo[i, j] = (mean, std, popt[1], popt[2]) + fittedFunc = fitFunctions.gaussian(bins, *popt) + #ax.plot(bins, fittedFunc, color="b") + except Exception as e: + print(f"An exception occurred: {e}") + logger.error(f"An exception occurred: {e}") + pass + + ax.set_xlabel("energy (keV)") + ax.set_title("pixel %d,%d in slice, small cluster cuts" % (i, j)) + plt.figtext(0.7, 0.8, "%d entries" % (nPixelClusters)) + plt.figtext(0.7, 0.75, "mu %0.2f" % (mu)) + plt.figtext(0.7, 0.7, "sigma %0.2f" % (sigma)) + fileName = "%s/%s_r%d_c%d_r%d_c%d_%s_E.png" % (fileInfo.outputDir, fileInfo.className, fileInfo.run, fileInfo.camera, i, j, fileInfo.label) + logger.info("Writing plot: " + fileName) + plt.savefig(fileName) + plt.close() + return fitInfo + +#fitIndex=3 +def analysis_two(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) + goodClusters = ancillaryMethods.goodClusters(clusters, i, j, nPixelCut=4, isSquare=1) + if len(goodClusters) <5: + print("too few clusters in slice pixel %d, %d: %d" %(i,j, len(goodClusters))) + continue + + energies = ancillaryMethods.getClusterEnergies(goodClusters) + photonEcut = np.bitwise_and(energies>lowEnergyCut, energies0).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 = 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)) + 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 = 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 + +def sliceToDetector(sliceRow, sliceCol, sliceCoordinates): + return sliceRow + sliceCoordinates[0][0], sliceCol + sliceCoordinates[1][0]