Merge pull request #87 from slaclab/beamtime_april_15_2024

philip development work for may beamtime

calibrationSuite/ancillaryMethods.py

@@ -52,8 +52,37 @@ 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, 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
     mCut = clusters[:,:,1] == module
     pixelRowCol = np.bitwise_and((clusters[:, :, 2] == row), (clusters[:, :, 3] == col))
     if isSquare is None:

calibrationSuite/basicSuiteScript.py

@@ -262,42 +262,44 @@ 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
         ## cut keeps photons in common mode - e.g. set to <<1 photon
 
         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):
         ## this takes a 2d frame
         ## 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][

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,17 +69,17 @@ 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
         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

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()
 
@@ -78,7 +76,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 +133,20 @@ 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))
+        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()))
 
@@ -148,50 +164,55 @@ 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()
+        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)
+                    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)))
+                        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()
+                    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
+                    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__,
@@ -220,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])

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)

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,7 +198,15 @@ 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)
+    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)
 
     esp.nGoodEvents = 0
@@ -215,11 +224,23 @@ def analyze_h5(self, dataFile, label):
                 ##print(ec)
                 continue
         frames = esp.getRawData(evt, gainBitsMasked=True)
+        if zeroLowGain or zeroHighGain:
+            frames = esp.getRawData(evt, gainBitsMasked=False)
+            if zeroLowGain:
+                gainFilter = frames>=esp.g0cut
+            else:
+                gainFilter = ~ (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)))

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)")
@@ -101,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()