Merge pull request #130 from slaclab/beamtime_september_1_2024

merging in changes from beamtime september 1 2024

calibrationSuite/fitFunctions.py

@@ -16,6 +16,7 @@
 
 logger = logging.getLogger(__name__)
 
+sqrt2pi = np.sqrt(2*np.pi)
 
 def linear(x, a, b):
     return a * x + b
@@ -34,20 +35,20 @@ def gaussian(x, a, mu, sigma):
 
 
 def gaussianArea(a, sigma):
-    return a * sigma * 6.28
+    return a * sigma * sqrt2pi
 
 
 def estimateGaussianParametersFromUnbinnedArray(flatData):
     sigma = flatData.std()
     entries = len(flatData)
     ## will crash if sigma is 0
-    return entries / (sigma * 6.28), flatData.mean(), sigma
+    return entries / (sigma * sqrt2pi), flatData.mean(), sigma
 
 
 def estimateGaussianParametersFromXY(x, y):
     mean, sigma = getHistogramMeanStd(x, y)
     ## will crash if sigma is 0
-    return sum(y) / (sigma * 6.28), mean, sigma
+    return sum(y) / (sigma * sqrt2pi), mean, sigma
 
 
 def getHistogramMeanStd(binCenters, counts):

config_files/epixMSuiteConfig.py

@@ -12,27 +12,28 @@
 # experimentHash = {'exp':'mfxx1005021', 'location':'MfxEndstation', 'fluxSource':'MFX-USR-DIO', 'fluxChannels':[11], 'fluxSign':-1}
 # experimentHash = {'exp':'rixc00121', 'location':'RixEndstation',
 singlePixelArray = []
-for i in range(0, 3):
-    singlePixelArray.append([i, 10, 10])
-    singlePixelArray.append([i, 180, 10])
-    singlePixelArray.append([i, 10, 200])
-    singlePixelArray.append([i, 180, 200])
-    singlePixelArray.append([i, 10, 380])
-    singlePixelArray.append([i, 180, 380])
+##for i in [0, 2, 3]:##range(0, 3):
+for i in [2, 3]:##range(0, 3):
+    singlePixelArray.append([i, 55, 10])
+    singlePixelArray.append([i, 120, 10])
+    singlePixelArray.append([i, 55, 200])
+    singlePixelArray.append([i, 120, 200])
+    singlePixelArray.append([i, 55, 300])
+    singlePixelArray.append([i, 120, 300])
 
 experimentHash = {
     "detectorType": "epixm",
     "detectorVersion":1,## new firmware
     "exp": "rixx1005922",
     "location": "RixEndstation",
-    "analyzedModules": [0, 2, 3],
+    "analyzedModules": [2, 3],
     "seedCut": 40,  ## pure guess
     "neighborCut": 10,  ##pure guess
-    # "fluxSource": "MfxDg1BmMon",
-    "fluxSource": "MfxDg2BmMon",
-    "fluxChannels": [15],
-    "fluxSign": 1,  ## for dg2
-    # "fluxSign": -1,
+    "fluxSource": "MfxDg1BmMon",
+    #"fluxSource": "MfxDg2BmMon",
+    "fluxChannels": [15],## or 11 if we see saturation
+    #"fluxSign": 1,  ## for dg2
+    "fluxSign": -1,
     "singlePixels": singlePixelArray,
     # 'ROIs':['module0', 'module2', 'module4', 'module6', 'module10','module12', 'module14']
     # 'ROIs':['roiFromSwitched_e557_rmfxx1005021']

standalone_scripts/analyze_npy.py

@@ -35,22 +35,29 @@ def __init__(self):
             "g1r2": 5,
             "g0max": 6,
             "g1min": 7,
+            "min": 8,
+            "max": 9,
+            "g0Ped": 10,
+            "g1Ped": 11,
+            "g0Gain": 12,
+            "offset": 13
         }
         self.dataRanges = {  ##"g0slope":[0,1000],
-            "g0slope": [0, 5],
+            ##"g0slope": [-5, 5],
+            "g0slope": [-1, 0], ##just xtalk...
             "g0intercept": [0, 10000],
             "g0r2": [0.9, 1.0],
-            "g1slope": [0, 0.1],
+            "g1slope": [-.1, 0.1],
             "g1intercept": [0, 10000],
             "g1r2": [0.9, 1.0],
-            "g0max": [10000, 16384],
-            "g1min": [0, 16384],
+            "g0max": [10000, 32000],
+            "g1min": [0, 32000],
         }
 
 
 class AnalyzeOneScan(object):
-    def __init__(self, scanObj, statsArray, data, label, ratio=False):
-        self.data = data
+    def __init__(self, scanObj, statsArray, data, module, label, ratio=False):
+        self.data = data[module]
         self.label = label
         self.ratio = ratio
         self.statsArray = statsArray
@@ -61,7 +68,10 @@ def __init__(self, scanObj, statsArray, data, label, ratio=False):
     def analyze(self):
         if self.ratio:
             for array in self.statsArray:
-                self.plotRatio(*tuple(array))
+                try:
+                    self.plotRatio(*tuple(array))
+                except:
+                    print("problem with", array, "in plotRatio, probably divide by zero")
         else:
             for array in self.statsArray:
                 if len(array) == 2:
@@ -209,12 +219,13 @@ def analyze(self):
 
 if __name__ == "__main__":
     f = sys.argv[1]
+    module = eval(sys.argv[2])
     statsArray = None
     plainStatsArray = None
     ratioStatsArray = None
 
     try:
-        statsArray = [sys.argv[2].split(",")]
+        statsArray = [sys.argv[3].split(",")]
     except:
         pass
 
@@ -251,9 +262,9 @@ def analyze(self):
                 plainStatsArray = statsArray
 
         if plainStatsArray is not None:
-            a = AnalyzeOneScan(scanObj, plainStatsArray, data, label)
+            a = AnalyzeOneScan(scanObj, plainStatsArray, data, module, label)
             a.analyze()
 
         if ratioStatsArray is not None:
-            a = AnalyzeOneScan(scanObj, ratioStatsArray, data, label, ratio=True)
+            a = AnalyzeOneScan(scanObj, ratioStatsArray, data, module, label, ratio=True)
             a.analyze()

suite_scripts/LinearityPlotsParallelSlice.py

@@ -193,7 +193,7 @@ def analyze_h5(self, dataFile, label):
 
     def analyze_h5_slice(self, dataFile, label):
         module = 2
-        nModules = 3
+        nModules = 4
         data = h5py.File(dataFile)
         fluxes = data["fluxes"][()]
         pixels = data["slice"][()]
@@ -204,7 +204,7 @@ def analyze_h5_slice(self, dataFile, label):
                 (nModules, rows, cols, 13)
             )  ##g0 slope, intercept, r2; g1 x3; max, min, g0Ped, g1Ped, g0Gain, g1Gain, offset
 
-        for module in [1, 2]:
+        for module in lpp.analyzedModules:
             for i in range(rows):
                 # so we can end early on testing runs:
                 if self.special is not None and "testing" in self.special and i >= self.maxNevents:
@@ -243,7 +243,7 @@ def analyze_h5_slice(self, dataFile, label):
                             )
                             if True:
                                 if i % 10 == 0 and j % 10 == 0:
-                                    print(i, j, fitPar, r2, 0)
+                                    print(module, i, j, fitPar, r2, 0)
                             ##np.save("temp_r%dc%d_x.py" %(i,j), fluxes[g0])
                             ##np.save("temp_r%dc%d_y.py" %(i,j), y)
                             ##np.save("temp_r%dc%d_func.py" %(i,j), fitFunc)
@@ -277,7 +277,9 @@ def analyze_h5_slice(self, dataFile, label):
                                 logger.info("empty profile for %d, %d" % (i, j))
                         if x is not None:
                             fitPar, covar, fitFunc, r2 = fitFunctions.fitLinearUnSaturatedData(x, y)
-                            print(i, j, fitPar, r2, 1)
+                            if True:
+                                if i % 10 == 0 and j % 10 == 0:
+                                    print(module, i, j, fitPar, r2, 1)
                             self.fitInfo[module, i, j, 3:5] = fitPar
                             self.fitInfo[module, i, j, 5] = r2
                             self.fitInfo[module, i, j, 7] = y_g1_min

suite_scripts/simplePhotonCounter.py

@@ -63,6 +63,8 @@
 
         if spc.fakePedestal is not None:
             rawFrames = spc.getRawData(evt)
+            if rawFrames is None:
+                continue
             frames = rawFrames.astype("float") - spc.fakePedestal
             frames /= gain  ## psana may not have the right default
         else:

suite_scripts/simpleTestScript.py

@@ -0,0 +1,30 @@
+from psana import *
+import sys
+
+##N.b. - this is psana2
+
+exp = 'rixc00122'
+run = eval(sys.argv[1])
+ds = DataSource(exp=exp,run=run,
+                intg_det='archon', max_events=666666)
+
+smd = ds.smalldata(filename="foo.h5")
+print(help(smd.event))
+
+myrun = next(ds.runs())
+det = myrun.Detector('archon')
+
+nGood = 0
+nBad = 0
+for nevt, evt in enumerate(myrun.events()):
+    ##print(dir(evt))
+    ##print(help(evt))
+    if det.raw.raw(evt) is None:
+        if nBad<3: print("event %d is None" %(nevt))
+        nBad += 1
+    else:
+        if nGood<3: print("good event")
+        nGood += 1
+
+
+print("good vs evil:", nGood, nBad)

tests/test_FitFunctions.py

@@ -63,7 +63,7 @@ def test_gaussian(a, mu, sigma, expected):
     np.testing.assert_allclose(result, expected)
 
 
-@pytest.mark.parametrize("a, sigma, expected", [(10, 0.5, 31.4)])
+@pytest.mark.parametrize("a, sigma, expected", [(10, 0.5, 12.533141373155)])
 def test_gaussianArea(a, sigma, expected):
     result = gaussianArea(a, sigma)
     assert np.isclose(result, expected)
@@ -73,7 +73,7 @@ def test_estimateGaussianParametersFromUnbinnedArray():
     flatData = np.array([1, 2, 3, 4, 5])
     result_amp, result_mean, result_sigma = estimateGaussianParametersFromUnbinnedArray(flatData)
     expected_amp, expected_mean, expected_sigma = (
-        0.5629831060402448,
+        1.4104739588693909,
         3.0,
         1.4142135623730951,
     )
@@ -87,7 +87,7 @@ def test_estimateGaussianParametersFromXY():
     y = np.array([2, 4, 6, 4, 2])
     result_amp, result_mean, result_sigma = estimateGaussianParametersFromXY(x, y)
     expected_amp, expected_mean, expected_sigma = (
-        2.482238418490429,
+        6.21889469048404,
         3.0,
         1.1547005383792515,
     )