From 8c599f18c4fc7d32cd04c9b99daf11d000d7ff45 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Fri, 19 Apr 2024 18:29:48 -0700 Subject: [PATCH 1/6] Make cluster analysis much faster, better --- calibrationSuite/ancillaryMethods.py | 19 ++++ calibrationSuite/detectorInfo.py | 2 +- suite_scripts/AnalyzeH5.py | 107 +++++++++++-------- suite_scripts/CalcNoiseAndMean.py | 17 ++- suite_scripts/SimpleClustersParallelSlice.py | 2 +- 5 files changed, 96 insertions(+), 51 deletions(-) diff --git a/calibrationSuite/ancillaryMethods.py b/calibrationSuite/ancillaryMethods.py index 6c3fe35..f0c2f8e 100644 --- a/calibrationSuite/ancillaryMethods.py +++ b/calibrationSuite/ancillaryMethods.py @@ -52,8 +52,27 @@ def plotProfile(x, y, yErr): def selectedClusters(clusters, row, col, lowEnerygCut, highEnergyCut, nPixelCut=4, isSquare=1): pass +def getEnergeticClusters(clusters): + ## expects [events, maxClusters, nClusterElements] + ## returns [nEnergeticClusters, nClusterElements] + return clusters[clusters[:,:,0]>0] + +def getSmallSquareClusters(clusters, nPixelCut=4): + smallClusters = getSmallClusters(clusters, nPixelCut=4) + return getSquareClusters(smallClusters) + +def getSmallClusters(clusters, nPixelCut=4): + return clusters[clusters[:,4] < nPixelCut] + +def getSquareClusters(clusters): + return clusters[clusters[:,5]==1] + +def getMatchedClusters(clusters, m, row, col): + matched = np.bitwise_and.reduce([(clusters[:,1]==m), (clusters[:,2]==row), clusters[:,3]==col]) + return clusters[matched] def goodClusters(clusters, module, row, col, nPixelCut=4, isSquare=None): + ## this is too slow mCut = clusters[:,:,1] == module pixelRowCol = np.bitwise_and((clusters[:, :, 2] == row), (clusters[:, :, 3] == col)) if isSquare is None: diff --git a/calibrationSuite/detectorInfo.py b/calibrationSuite/detectorInfo.py index 0403d35..f3f3f21 100644 --- a/calibrationSuite/detectorInfo.py +++ b/calibrationSuite/detectorInfo.py @@ -79,7 +79,7 @@ def setup_epixM(self, version=0): self.negativeGain = True self.aduPerKeV = 666 self.seedCut = 4 ## hopefully can be lowered - self.neighborCut = 0.5 ## ditto + self.neighborCut = 0.25 ## ditto def setup_rixsCCD(self, mode='1d', version=0): self.nTestPixelsPerBank = 36 diff --git a/suite_scripts/AnalyzeH5.py b/suite_scripts/AnalyzeH5.py index 8d5b979..24a3dea 100644 --- a/suite_scripts/AnalyzeH5.py +++ b/suite_scripts/AnalyzeH5.py @@ -78,7 +78,13 @@ def analyze(self): def clusterAnalysis(self): clusters = None energyHist = None - clusters = np.concatenate([h5["clusterData"][()] for h5 in self.h5Files]) + clusters = np.concatenate( + [ancillaryMethods.getEnergeticClusters(h5["clusterData"][()]) + for h5 in self.h5Files]) + ## getEnergeticClusters + ## makes [events, maxClusters, nClusterElements] -> [m, n] + ## makes memory violation a bit less likely too + try: energyHist = np.concatenate(energyHist, h5["energyHistogram"][()]) except: @@ -129,8 +135,18 @@ def analyzeSimpleClusters(self, clusters): ## should make a dict for this info and use below instead of indices ax = plt.subplot() - energy = clusters[:, :, 0] ##.flatten() - maximumModule = int(clusters[:,:,1].max()) + energy = clusters[:, 0] ##.flatten() + maximumModule = int(clusters[:,1].max()) + analyzedModules = np.unique(clusters[:,1]).astype('int') + print("analyzing modules", analyzedModules) + + ##rows = self.sliceEdges[0] + ##cols = self.sliceEdges[1] + ## doesn't exist in h5 yet so calculate dumbly instead + rows = int(clusters[:,2].max())+1 + cols = int(clusters[:,3].max())+1 + print("appear to have a slice with %d rows, %d cols" %(rows, cols)) + print("mean energy above 0:" + str(energy[energy > 0].mean())) logger.info("mean energy above 0:" + str(energy[energy > 0].mean())) @@ -148,49 +164,48 @@ def analyzeSimpleClusters(self, clusters): logger.info("Wrote file: " + figFileName) plt.close() - rows = self.sliceEdges[0] - cols = self.sliceEdges[1] - m = 1## temp hack, Kaz's favorite asic, off by 1 - fitInfo = np.zeros((maximumModule, 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, m, i, j, nPixelCut=3, isSquare=1) - if len(currGoodClusters) < 5: - print("too few clusters in slice pixel %d, %d, %d: %d" % (m, i, j, len(currGoodClusters))) - logger.info("too few clusters in slice pixel %d, %d, %d: %d" % (m, 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,%d has about %d photons" % (m, i, j, nPixelClusters)) - logger.info("pixel %d,%d,%d has about %d photons" % (m, 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) - fitInfo[m, i, j] = mean, std, area, mu, sigma - - ax.set_xlabel("energy (keV)") - ax.set_title("pixel %d,%d,%d, small cluster cuts" % (m, 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_m%d_r%d_c%d_%s_E.png" % ( - self.outputDir, - self.__class__.__name__, - self.run, - self.camera, - m, - detRow, - detCol, - self.label, - ) - plt.savefig(figFileName) - logger.info("Wrote file: " + figFileName) - plt.close() + fitInfo = np.zeros((maximumModule+1, rows, cols, 5)) ## mean, std, area, mu, sigma + smallSquareClusters = ancillaryMethods.getSmallSquareClusters(clusters, nPixelCut=3) + for m in analyzedModules: + for i in range(rows): + for j in range(cols): + detRow, detCol = self.sliceToDetector(i, j) + ax = plt.subplot() + currGoodClusters = ancillaryMethods.getMatchedClusters(smallSquareClusters, m, i, j) + if len(currGoodClusters) < 5: + print("too few clusters in slice pixel %d, %d, %d: %d" % (m, i, j, len(currGoodClusters))) + logger.info("too few clusters in slice pixel %d, %d, %d: %d" % (m, 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,%d has about %d photons" % (m, i, j, nPixelClusters)) + logger.info("pixel %d,%d,%d has about %d photons" % (m, 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) + fitInfo[m, i, j] = mean, std, area, mu, sigma + + ax.set_xlabel("energy (keV)") + ax.set_title("pixel %d,%d,%d, small cluster cuts" % (m, 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_m%d_r%d_c%d_%s_E.png" % ( + self.outputDir, + self.__class__.__name__, + self.run, + self.camera, + m, + detRow, + detCol, + self.label, + ) + plt.savefig(figFileName) + logger.info("Wrote file: " + figFileName) + plt.close() npyFileName = "%s/%s_r%d_c%d_%s_fitInfo.npy" % ( self.outputDir, diff --git a/suite_scripts/CalcNoiseAndMean.py b/suite_scripts/CalcNoiseAndMean.py index 2a23f69..9d1fa20 100755 --- a/suite_scripts/CalcNoiseAndMean.py +++ b/suite_scripts/CalcNoiseAndMean.py @@ -59,11 +59,19 @@ def __init__(self): if frames is None: continue if "noCommonMode" in cn.special: - frames = cn.noCommonModeCorrection(frames[0]) + frames = cn.noCommonModeCorrection(frames) elif "rowCommonMode" in cn.special: - frames = cn.rowCommonModeCorrection(frames[0], 2.0) + if self.fakePedestal is None: + print("row common mode needs reasonable pedestal") + raise Exception + frames = frames - cn.fakePedestal + frames = cn.rowCommonModeCorrection3d(frames, 2.0) elif "colCommonMode" in cn.special: - frames = cn.colCommonModeCorrection(frames[0], 2.0) + if self.fakePedestal is None: + print("col common mode needs reasonable pedestal") + raise Exception + frames = frames - cn.fakePedestal + frames = cn.colCommonModeCorrection3d(frames, 2.0) elif "regionCommonMode" in cn.special: frames = cn.regionCommonModeCorrection(frames[0], cn.regionSlice, commonModeCut) else: @@ -107,6 +115,9 @@ def __init__(self): else: pass + if cn.fakePedestal is not None: + cn.label += "_fakePdestal" + meanRmsFileName = "%s/CalcNoiseAndMean_%s_rms_r%d_step%s.npy" % (cn.outputDir, cn.label, cn.run, nstep) np.save(meanRmsFileName, noise) meanFileName = "%s/CalcNoiseAndMean_mean_%s_r%d_step%s.npy" % (cn.outputDir, cn.label, cn.run, nstep) diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index c9de2bb..d87c83e 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -128,7 +128,7 @@ def analyze_h5(self, dataFile, label): smd = sic.ds.smalldata(filename="%s/%s_c%d_r%d_n%d.h5" % (sic.outputDir, sic.className, sic.camera, sic.run, size)) ## 50x50 pixels, 3x3 clusters, 10% occ., 2 sensors - maxClusters = 1000##int(50 * 50 / 3 / 3 * 0.1 * 2) + maxClusters = 10000##int(50 * 50 / 3 / 3 * 0.1 * 2) try: seedCut = sic.detectorInfo.seedCut except: From 64ed343b86ed6f0c91301e496a41bec94dca7339 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Fri, 26 Apr 2024 09:15:43 -0700 Subject: [PATCH 2/6] updates for faster single photon analysis, 3d common mode --- calibrationSuite/ancillaryMethods.py | 16 ++++-- calibrationSuite/basicSuiteScript.py | 22 ++++---- calibrationSuite/detectorInfo.py | 6 +-- suite_scripts/AnalyzeH5.py | 55 +++++++++++--------- suite_scripts/LinearityPlotsParallelSlice.py | 10 +++- suite_scripts/SimpleClustersParallelSlice.py | 8 +-- suite_scripts/archonSuiteConfig.py | 37 +++++++++++++ suite_scripts/runAnalyzeH5.py | 3 +- 8 files changed, 110 insertions(+), 47 deletions(-) create mode 100755 suite_scripts/archonSuiteConfig.py diff --git a/calibrationSuite/ancillaryMethods.py b/calibrationSuite/ancillaryMethods.py index f0c2f8e..a9b51ba 100644 --- a/calibrationSuite/ancillaryMethods.py +++ b/calibrationSuite/ancillaryMethods.py @@ -67,9 +67,19 @@ def getSmallClusters(clusters, nPixelCut=4): def getSquareClusters(clusters): return clusters[clusters[:,5]==1] -def getMatchedClusters(clusters, m, row, col): - matched = np.bitwise_and.reduce([(clusters[:,1]==m), (clusters[:,2]==row), clusters[:,3]==col]) - return clusters[matched] +def getMatchedClusters(clusters, dimension, n): + if dimension == 'column': + return clusters[(clusters[:,3]==n)] + if dimension == 'row': + return clusters[(clusters[:,2]==n)] + if dimension == 'module': + return clusters[(clusters[:,1]==n)] + return None + + +##def getMatchedClusters(clusters, m, row, col): +## matched = np.bitwise_and.reduce([(clusters[:,1]==m), (clusters[:,2]==row), clusters[:,3]==col]) +## return clusters[matched] def goodClusters(clusters, module, row, col, nPixelCut=4, isSquare=None): ## this is too slow diff --git a/calibrationSuite/basicSuiteScript.py b/calibrationSuite/basicSuiteScript.py index 8de264a..753930b 100755 --- a/calibrationSuite/basicSuiteScript.py +++ b/calibrationSuite/basicSuiteScript.py @@ -262,10 +262,12 @@ def regionCommonModeCorrection(self, frame, region, arbitraryCut=1000): def rowCommonModeCorrection3d(self, frames, arbitraryCut=1000): for module in self.analyzedModules: frames[module] = self.rowCommonModeCorrection(frames[module], arbitraryCut) - + return frames + def colCommonModeCorrection3d(self, frames, arbitraryCut=1000): for module in self.analyzedModules: frames[module] = self.colCommonModeCorrection(frames[module], arbitraryCut) + return frames def rowCommonModeCorrection(self, frame, arbitraryCut=1000): ## this takes a 2d object @@ -273,21 +275,21 @@ def rowCommonModeCorrection(self, frame, arbitraryCut=1000): for r in range(self.detectorInfo.nRows): colOffset = 0 - ##for b in range(0, self.detNbanks): - for b in range(0, 2): + for b in range(0, self.detectorInfo.nBanksCol): + ##for b in range(0, 2): try: rowCM = np.median( - frame[r, colOffset : colOffset + self.detColsPerBank][ - frame[r, colOffset : colOffset + self.detColsPerBank] < arbitraryCut + frame[r, colOffset : colOffset + self.detectorInfo.nColsPerBank][ + frame[r, colOffset : colOffset + self.detectorInfo.nColsPerBank] < arbitraryCut ] ) - frame[r, colOffset : colOffset + self.detColsPerBank] -= rowCM + frame[r, colOffset : colOffset + self.detectorInfo.nColsPerBank] -= rowCM except: rowCM = -666 print("rowCM problem") logger.error("rowCM problem") - print(frame[r, colOffset : colOffset + self.detColsPerBank]) - colOffset += self.detColsPerBank + print(frame[r, colOffset : colOffset + self.detectorInfo.nColsPerBank]) + colOffset += self.detectorInfo.nColsPerBank return frame def colCommonModeCorrection(self, frame, arbitraryCut=1000): @@ -295,9 +297,9 @@ def colCommonModeCorrection(self, frame, arbitraryCut=1000): ## cut keeps photons in common mode - e.g. set to <<1 photon ##rand = np.random.random() - for c in range(self.detCols): + for c in range(self.detectorInfo.nCols): rowOffset = 0 - for b in range(0, self.detNbanksCol): + for b in range(0, self.detectorInfo.nBanksRow): try: colCM = np.median( frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c][ diff --git a/calibrationSuite/detectorInfo.py b/calibrationSuite/detectorInfo.py index f3f3f21..0fc5bad 100644 --- a/calibrationSuite/detectorInfo.py +++ b/calibrationSuite/detectorInfo.py @@ -52,7 +52,7 @@ def setup_epixhr(self, version=0): self.nRows = 288 self.nCols = 284 self.nColsPerBank = 96 - self.detNbanks = int(self.nCols / self.nColsPerBank) + self.nBanksRow = int(self.nCols / self.nColsPerBank) self.nBanksCol = 2 self.nRowsPerBank = int(self.nRows / self.nBanksCol) @@ -69,10 +69,10 @@ def setup_epixM(self, version=0): ## per module (aka asic) self.nRows = 192 self.nBanksRow = 4 - self.nRowsPerBank = self.nRows/self.nBanksRow + self.nRowsPerBank = int(self.nRows/self.nBanksRow) self.nCols = 384 self.nBanksCol = 6 - self.nColsPerBank = self.nCols/self.nBanksCol + self.nColsPerBank = int(self.nCols/self.nBanksCol) self.preferredCommonMode = 'rowCommonMode'## guess self.clusterShape = [3,3] self.gainMode = None ## may want to know about default, softHigh, softLow diff --git a/suite_scripts/AnalyzeH5.py b/suite_scripts/AnalyzeH5.py index 24a3dea..974b860 100644 --- a/suite_scripts/AnalyzeH5.py +++ b/suite_scripts/AnalyzeH5.py @@ -12,7 +12,6 @@ import calibrationSuite.fitFunctions as fitFunctions import calibrationSuite.ancillaryMethods as ancillaryMethods -##import seaborn as sns import matplotlib.pyplot as plt from matplotlib.ticker import AutoMinorLocator @@ -32,7 +31,6 @@ class AnalyzeH5(object): def __init__(self): - print("in init") args = ArgumentParser().parse_args() @@ -146,7 +144,9 @@ def analyzeSimpleClusters(self, clusters): rows = int(clusters[:,2].max())+1 cols = int(clusters[:,3].max())+1 print("appear to have a slice with %d rows, %d cols" %(rows, cols)) - + self.sliceCoordinates = [[0, rows], [0, cols]] ## temp - get from h5 + self.sliceEdges = [rows, cols] + print("mean energy above 0:" + str(energy[energy > 0].mean())) logger.info("mean energy above 0:" + str(energy[energy > 0].mean())) @@ -164,14 +164,17 @@ def analyzeSimpleClusters(self, clusters): logger.info("Wrote file: " + figFileName) plt.close() + verbose = False fitInfo = np.zeros((maximumModule+1, rows, cols, 5)) ## mean, std, area, mu, sigma smallSquareClusters = ancillaryMethods.getSmallSquareClusters(clusters, nPixelCut=3) for m in analyzedModules: + modClusters = ancillaryMethods.getMatchedClusters(smallSquareClusters, 'module', m) for i in range(rows): + rowModClusters = ancillaryMethods.getMatchedClusters(modClusters, 'row', i) + for j in range(cols): detRow, detCol = self.sliceToDetector(i, j) - ax = plt.subplot() - currGoodClusters = ancillaryMethods.getMatchedClusters(smallSquareClusters, m, i, j) + currGoodClusters = ancillaryMethods.getMatchedClusters(rowModClusters, 'column', j) if len(currGoodClusters) < 5: print("too few clusters in slice pixel %d, %d, %d: %d" % (m, i, j, len(currGoodClusters))) logger.info("too few clusters in slice pixel %d, %d, %d: %d" % (m, i, j, len(currGoodClusters))) @@ -184,29 +187,32 @@ def analyzeSimpleClusters(self, clusters): photonRegion = energies[photonEcut] mean = photonRegion.mean() std = photonRegion.std() + ax = plt.subplot() a, mu, sigma = self.histogramAndFitGaussian(ax, energies, self.nBins) area = fitFunctions.gaussianArea(a, sigma) fitInfo[m, i, j] = mean, std, area, mu, sigma - - ax.set_xlabel("energy (keV)") - ax.set_title("pixel %d,%d,%d, small cluster cuts" % (m, 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_m%d_r%d_c%d_%s_E.png" % ( - self.outputDir, - self.__class__.__name__, - self.run, - self.camera, - m, - detRow, - detCol, - self.label, - ) - plt.savefig(figFileName) - logger.info("Wrote file: " + figFileName) + if i%13==1 and j%17==1: + ## don't save a zillion plots when analyzing full array + ## should add singlePixelArray here + ax.set_xlabel("energy (keV)") + ax.set_title("pixel %d,%d,%d, small cluster cuts" % (m, 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_m%d_r%d_c%d_%s_E.png" % ( + self.outputDir, + self.__class__.__name__, + self.run, + self.camera, + m, + detRow, + detCol, + self.label, + ) + plt.savefig(figFileName) + logger.info("Wrote file: " + figFileName) plt.close() - + npyFileName = "%s/%s_r%d_c%d_%s_fitInfo.npy" % ( self.outputDir, self.__class__.__name__, @@ -235,7 +241,6 @@ def analyzeSimpleClusters(self, clusters): 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]) diff --git a/suite_scripts/LinearityPlotsParallelSlice.py b/suite_scripts/LinearityPlotsParallelSlice.py index e3324f5..87d77b2 100755 --- a/suite_scripts/LinearityPlotsParallelSlice.py +++ b/suite_scripts/LinearityPlotsParallelSlice.py @@ -80,11 +80,17 @@ def plotAutorangingData_profile(self, g0s, g1s, g0Fluxes, g1Fluxes, label, parti if partialMode != 1 and len(g0s[i]) > 0: x = g0Fluxes[i] y = g0s[i] - sns.regplot(x=x, y=y, x_bins=40, marker=".", ax=ax, order=order) ##add fit_reg=None for no plot + ##yMaxPlot = y.max() + sns.regplot(x=x, y=y, x_bins=40, marker=".", ax=ax, order=order, truncate=True) ##add fit_reg=None for no plot + ## truncate added to keep epixM plot lines from messing with y limits if partialMode != 0 and len(g1s[i]) > 0: x = g1Fluxes[i] y = g1s[i] - sns.regplot(x=x, y=y, x_bins=40, marker=".", ax=ax, order=order) ##add fit_reg=None for no plot + ##try:## using truncate instead + ##yMaxPlot = max(y.max(), yMaxPlot) + ##except: + ##yMaxPlot = y.max() + sns.regplot(x=x, y=y, x_bins=40, marker=".", ax=ax, order=order, truncate=True) ##add fit_reg=None for no plot plt.xlabel("wave8 flux (ADU)") if "raw" in label: plt.ylabel("Red medium, blue low (ADU)") diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index d87c83e..dfc31c8 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -125,7 +125,7 @@ def analyze_h5(self, dataFile, label): sys.exit(0) sic.setupPsana() - smd = sic.ds.smalldata(filename="%s/%s_c%d_r%d_n%d.h5" % (sic.outputDir, sic.className, sic.camera, sic.run, size)) + smd = sic.ds.smalldata(filename="%s/%s_%s_c%d_r%d_n%d.h5" % (sic.outputDir, sic.className, sic.label, sic.camera, sic.run, size)) ## 50x50 pixels, 3x3 clusters, 10% occ., 2 sensors maxClusters = 10000##int(50 * 50 / 3 / 3 * 0.1 * 2) @@ -138,6 +138,8 @@ def analyze_h5(self, dataFile, label): except: neighborCut = 0.5 + print("using seed, neighbor cuts", seedCut, neighborCut) + sic.clusterElements = ["energy", "module", "row", "col", "nPixels", "isSquare"] nClusterElements = len(sic.clusterElements) @@ -211,7 +213,7 @@ def analyze_h5(self, dataFile, label): if sic.special is not None: if "regionCommonMode" in sic.special: - frame = sic.regionCommonModeCorrection(frame, sic.regionSlice, 2.0) + frames = sic.regionCommonModeCorrection(frames, sic.regionSlice, 2.0) if "rowCommonMode" in sic.special: frames = sic.rowCommonModeCorrection3d(frames, 2.0) if "colCommonMode" in sic.special: @@ -219,7 +221,7 @@ def analyze_h5(self, dataFile, label): if frames is None: print("common mode killed frames???") - continue + raise Exception flux = sic.flux if sic.useFlux and flux is None: diff --git a/suite_scripts/archonSuiteConfig.py b/suite_scripts/archonSuiteConfig.py new file mode 100755 index 0000000..5b33989 --- /dev/null +++ b/suite_scripts/archonSuiteConfig.py @@ -0,0 +1,37 @@ +############################################################################## +## This file is part of 'SLAC Beamtime Calibration Suite'. +## It is subject to the license terms in the LICENSE.txt file found in the +## top-level directory of this distribution and at: +## https://confluence.slac.stanford.edu/display/ppareg/LICENSE.html. +## No part of 'SLAC Beamtime Calibration Suite', including this file, +## may be copied, modified, propagated, or distributed except according to +## the terms contained in the LICENSE.txt file. +############################################################################## +import numpy as np + +##experimentHash = {'exp':'mfxx1005021', 'location':'MfxEndstation', 'fluxSource':'MFX-USR-DIO', 'fluxChannels':[11], 'fluxSign':-1} +##experimentHash = {'exp':'rixc00121', 'location':'RixEndstation', +experimentHash = { + "detectorType": "archon", + "exp": "rixc00121", + "location": "RixEndstation", + "seedCut":30,##from old scripts + ##"fluxSource": "MfxDg1BmMon", + ## 'fluxChannels':[9, 10, 11, 12], + ##"fluxChannels": [15], + ## 'fluxChannels':[8, 14], + ##"fluxSign": -1, + "singlePixels": [ + [0, 10], + [0, 100], + [0, 1000], + [0, 2000], + ], + # 'ROIs':['module0', 'module2', 'module4', 'module6', 'module10','module12', 'module14'] + # 'ROIs':['roiFromSwitched_e557_rmfxx1005021'] + ## 'ROIs':['allHRasicPixels', 'goodboxROI']#'roiAbove7k_raw_r123'] + ##"ROIs": ["../data/XavierV4_2", "../data/OffXavierV4_2"], + ##"regionSlice": np.s_[270:288], +} +##more complex approach allowing run ranges +##fluxHash = {1:['MFX-USR-DIO', 11]} diff --git a/suite_scripts/runAnalyzeH5.py b/suite_scripts/runAnalyzeH5.py index 9338cf1..305a61c 100644 --- a/suite_scripts/runAnalyzeH5.py +++ b/suite_scripts/runAnalyzeH5.py @@ -23,8 +23,9 @@ for r in runs: globString = "%s/%s*r%s*.h5" % (basePath, analysisType, r) f = glob.glob(globString) + f0 = [x for x in f if 'part' not in x] try: - command += f[-1] + "," + command += f0[-1] + "," except: print("could not find %s" % (globString)) From 00ce6af11fdb16034804cba7af02a2b54e85cfb8 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Wed, 8 May 2024 13:53:15 -0700 Subject: [PATCH 3/6] Linearity plot update for Bhavna --- suite_scripts/LinearityPlotsParallelSlice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/suite_scripts/LinearityPlotsParallelSlice.py b/suite_scripts/LinearityPlotsParallelSlice.py index 87d77b2..35e5130 100755 --- a/suite_scripts/LinearityPlotsParallelSlice.py +++ b/suite_scripts/LinearityPlotsParallelSlice.py @@ -107,7 +107,7 @@ def plotAutorangingData_profile(self, g0s, g1s, g0Fluxes, g1Fluxes, label, parti ##lot of hackiness here plt.title(label + "_profile") - figFileName = "%s/%s_p%d_r%d_c%d_%s_profile.png" % (self.outputDir, self.className, i, self.run, self.camera, label) + figFileName = "%s/%s_p%d_mod%d_row%d_col%d_r%d_c%d_%s_profile.png" % (self.outputDir, self.className, i, p[0], p[1], p[2], self.run, self.camera, label) plt.savefig(figFileName) plt.close() From 632ade7aa874da03ba41e8af5873f21e04e1f630 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Fri, 10 May 2024 12:08:09 -0700 Subject: [PATCH 4/6] Last changes to support analyses for beamtime_april_15_2024 data --- suite_scripts/EventScanParallelSlice.py | 21 +++++++++++++++++++- suite_scripts/SimpleClustersParallelSlice.py | 15 +++++++++++++- suite_scripts/TimeScanParallelSlice.py | 21 +++++++++++++------- 3 files changed, 48 insertions(+), 9 deletions(-) diff --git a/suite_scripts/EventScanParallelSlice.py b/suite_scripts/EventScanParallelSlice.py index e13718e..6b168c0 100644 --- a/suite_scripts/EventScanParallelSlice.py +++ b/suite_scripts/EventScanParallelSlice.py @@ -197,7 +197,7 @@ def analyze_h5(self, dataFile, label): except: skip_283_check = False ## for running at MFX - h5FileName = "%s/%s_c%d_r%d_n%d.h5" % (esp.outputDir, esp.className, esp.camera, esp.run, size) + h5FileName = "%s/%s_c%d_r%d_%s_n%d.h5" % (esp.outputDir, esp.className, esp.camera, esp.run, esp.label, size) smd = esp.ds.smalldata(filename=h5FileName) esp.nGoodEvents = 0 @@ -215,11 +215,30 @@ def analyze_h5(self, dataFile, label): ##print(ec) continue frames = esp.getRawData(evt, gainBitsMasked=True) + zeroLowGain = False + zeroHighGain = False + if esp.special: + if "zeroLowGain" in esp.special: + zeroLowGain = True + elif "zeroHighGain" in esp.special: + zeroHighGain = True + if zeroLowGain or zeroHighGain: + frames = esp.getRawData(evt, gainBitsMasked=False) + if zeroLowGain: + gainFilter = frames>=esp.g0cut + else: + gainFilter = not (frames>=esp.g0cut) + frames[gainFilter] = 0 + frames = frames & esp.gainBitsMask + if frames is None: ##print("no frame") continue if esp.fakePedestal is not None: frames = frames.astype("float") - esp.fakePedestal + if zeroLowGain or zeroHighGain: + frames[gainFilter] = 0 + if esp.special is not None and "crazyModeForDionisio" in esp.special: print("crazy mode for Dionisio") plt.imshow(np.vstack(frames[1:3].clip(-500, 500))) diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index dfc31c8..8602025 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -191,15 +191,28 @@ def analyze_h5(self, dataFile, label): if not sic.fakeBeamCode and not sic.isBeamEvent(evt): continue - rawFrames = sic.getRawData(evt) + zeroLowGain = False + if sic.special and 'zeroLowGain' in sic.special: + rawFrames = sic.getRawData(evt, gainBitsMasked=False) + else: + rawFrames = sic.getRawData(evt) + + if rawFrames is None: continue + if zeroLowGain: g1 = rawFrames>=sic.g0cut frames = None if sic.fakePedestal is not None: frames = rawFrames.astype("float") - sic.fakePedestal + if zeroLowGain: frames[g1] = 0 elif pedestal is not None: frames = rawFrames.astype("float") - pedestal + if zeroLowGain: frames[g1] = 0 + ##else: + ##print("something is probably wrong, need a pedestal to cluster") + ##sys.exit(0) + if frames is not None and gain is not None: if sic.special is not None and "addFakePhotons" in sic.special: frames, nAdded = sic.addFakePhotons(frames, 0.01, 666*10, 10) diff --git a/suite_scripts/TimeScanParallelSlice.py b/suite_scripts/TimeScanParallelSlice.py index 034d5df..a5539ca 100755 --- a/suite_scripts/TimeScanParallelSlice.py +++ b/suite_scripts/TimeScanParallelSlice.py @@ -120,16 +120,16 @@ def analyze_h5(self, dataFile, norm, label): import h5py a = h5py.File(dataFile)[norm] - delays = np.array([k for k in a.keys()]) - print(delays) + delays = np.array([int(eval(k)) for k in a.keys()]) + print('delays:', delays) ##delays = delays.astype('float').astype("int") ##delays = delays.astype('float').astype("int") delays.sort() d = np.array([a[str(k)] for k in delays]) delays = np.array([d for d in delays]) - delays = [eval(d)/1000. for d in delays] + delays = [d/1000. for d in delays] ##delays /= 1000. - print(delays) + print('scaled delays', delays) runString = "_r%d" % (self.run) if norm != "slice": @@ -146,13 +146,15 @@ def analyze_h5(self, dataFile, norm, label): tsp = TimeScanParallel() print("have built a", tsp.className, "class") logger.info("have built a" + tsp.className + "class") - if tsp.file is not None: + print('tsp.file', tsp.file) + fileMadeByScript = tsp.file.split('/')[-1].startswith(tsp.className) + if tsp.file is not None and fileMadeByScript:##and tsp.psanaType != 0: ## added type for rogue tsp.analyze_h5(tsp.file, 'means', tsp.label) ## tsp.analyze_h5(tsp.file, 'ratios', tsp.label) tsp.analyze_h5(tsp.file, "slice", tsp.label) print("done with standalone analysis of %s, exiting" % (tsp.file)) logger.info("done with standalone analysis of %s, exiting" % (tsp.file)) - sys.exit() + sys.exit(0) tsp.setupPsana() tsp.use_281_for_old_data = False @@ -161,7 +163,12 @@ def analyze_h5(self, dataFile, norm, label): tsp.use_281_for_old_data = True print("using all event code 281 frames for old data") logger.info("using all event code 281 frames for old data") - + + if 'size' in dir(): ## check for rogue + print("size is", size) + else: + size = 666 + h5FileName = "%s/%s_%s_c%d_r%d_n%d.h5" % (tsp.outputDir, tsp.className, tsp.label, tsp.camera, tsp.run, size) smd = tsp.ds.smalldata(filename=h5FileName) From 7db0917e6927b83c767bd466aef320bfb4779ec8 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Fri, 10 May 2024 13:47:42 -0700 Subject: [PATCH 5/6] fixed a bug, improved a loop --- suite_scripts/EventScanParallelSlice.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/suite_scripts/EventScanParallelSlice.py b/suite_scripts/EventScanParallelSlice.py index 6b168c0..1f57da2 100644 --- a/suite_scripts/EventScanParallelSlice.py +++ b/suite_scripts/EventScanParallelSlice.py @@ -112,6 +112,7 @@ def plotData(self, rois, pixels, eventNumbers, dPulseId, label): ax.hist(pixels[i], 100, range=[pixels[i].min().astype("int"), pixels[i].max().astype("int")]) plt.xlabel("Pixel ADU") plt.title("Event scan projection of pixel %d" % (i)) + plt.yscale('log') pltFileName = "%s/%s_r%d_c%d_%s_pixel%d_hist.png" % ( self.outputDir, self.__class__.__name__, @@ -197,6 +198,14 @@ def analyze_h5(self, dataFile, label): except: skip_283_check = False ## for running at MFX + zeroLowGain = False + zeroHighGain = False + if esp.special: + if "zeroLowGain" in esp.special: + zeroLowGain = True + elif "zeroHighGain" in esp.special: + zeroHighGain = True + h5FileName = "%s/%s_c%d_r%d_%s_n%d.h5" % (esp.outputDir, esp.className, esp.camera, esp.run, esp.label, size) smd = esp.ds.smalldata(filename=h5FileName) @@ -215,22 +224,15 @@ def analyze_h5(self, dataFile, label): ##print(ec) continue frames = esp.getRawData(evt, gainBitsMasked=True) - zeroLowGain = False - zeroHighGain = False - if esp.special: - if "zeroLowGain" in esp.special: - zeroLowGain = True - elif "zeroHighGain" in esp.special: - zeroHighGain = True if zeroLowGain or zeroHighGain: frames = esp.getRawData(evt, gainBitsMasked=False) if zeroLowGain: gainFilter = frames>=esp.g0cut else: - gainFilter = not (frames>=esp.g0cut) + gainFilter = ~ (frames>=esp.g0cut) frames[gainFilter] = 0 frames = frames & esp.gainBitsMask - + if frames is None: ##print("no frame") continue From 185f0f87c28dc3d9237e87e5c8b888ea541b2e8b Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Fri, 10 May 2024 13:51:16 -0700 Subject: [PATCH 6/6] fixed a bug, improved a loop --- suite_scripts/SimpleClustersParallelSlice.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index 8602025..d03df98 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -184,15 +184,19 @@ def analyze_h5(self, dataFile, label): print("sic.det:", sic.det) pass + + zeroLowGain = False + if sic.special and 'zeroLowGain' in sic.special: + zeroLowGain = True + hSum = None for nevt, evt in enumerate(evtGen): if evt is None: continue if not sic.fakeBeamCode and not sic.isBeamEvent(evt): continue - - zeroLowGain = False - if sic.special and 'zeroLowGain' in sic.special: + + if zeroLowGain: rawFrames = sic.getRawData(evt, gainBitsMasked=False) else: rawFrames = sic.getRawData(evt)