Merge analyzeNpyScan into deveopment

calibrationSuite/argumentParser.py

@@ -55,6 +55,9 @@ def __init__(self):
         self.parser.add_argument(
             "-f", "--files", type=str, default=None, help="run analysis on file or comma-separated files"
         )
+        self.parser.add_argument(
+            "-F", "--nonXtcFiles", type=str, default=None, help="run initial analysis on file or comma-separated files"
+        )
         self.parser.add_argument("-cf", "--configFile", type=str, help="config file path, can be relative or full path")
         self.parser.add_argument("-L", "--label", type=str, help="analysis label")
         self.parser.add_argument("-t", "--threshold", help="threshold (ADU or keV or wave8) depending on --detObj")

calibrationSuite/basicSuiteScript.py

@@ -21,7 +21,10 @@
 
 logger = logging.getLogger(__name__)
 
-if os.getenv("foo") == "1":
+if os.getenv("foo") == "0":
+    print("non-psana")
+    from calibrationSuite.nonPsanaBase import PsanaBase
+elif os.getenv("foo") == "1":
     print("psana1")
     from calibrationSuite.psana1Base import PsanaBase
 else:
@@ -96,7 +99,9 @@ def getRawData(self, evt, gainBitsMasked=True, negativeGain=False):
         if frames is None:
             return None
 
-        nZero = frames.size - np.count_nonzero(frames)
+        nZero = 0
+        if not "skipZeroCheck" in dir(self):
+            nZero = frames.size - np.count_nonzero(frames)
         try:
             dz = self.nZero - nZero
             if dz != 0:
@@ -171,6 +176,9 @@ def sliceToDetector(self, sliceRow, sliceCol):  ## cp from AnalyzeH5: import?
     def getNswitchedPixels(self, data, region=None):
         return ((data >= self.g0cut) * 1).sum()
 
+    def getSwitchedPixels(self, data, region=None):
+        return data >= self.g0cut
+
     def dumpEventCodeStatistics(self):
         print(
             "have counted %d run triggers, %d DAQ triggers, %d beam events"
@@ -217,9 +225,9 @@ def rowCommonModeCorrection3d(self, frames, arbitraryCut=1000):
             frames[module] = self.rowCommonModeCorrection(frames[module], arbitraryCut)
         return frames
 
-    def colCommonModeCorrection3d(self, frames, arbitraryCut=1000):
+    def colCommonModeCorrection3d(self, frames, cut=1000, switchedPixels=None):
         for module in self.analyzedModules:
-            frames[module] = self.colCommonModeCorrection(frames[module], arbitraryCut)
+            frames[module] = self.colCommonModeCorrection(frames[module], cut, switchedPixels[module])
         return frames
 
     def rowCommonModeCorrection(self, frame, arbitraryCut=1000):
@@ -245,7 +253,7 @@ def rowCommonModeCorrection(self, frame, arbitraryCut=1000):
                 colOffset += self.detectorInfo.nColsPerBank
         return frame
 
-    def colCommonModeCorrection(self, frame, arbitraryCut=1000):
+    def colCommonModeCorrection(self, frame, cut=1000, switchedPixels=None):
         ## this takes a 2d frame
         ## cut keeps photons in common mode - e.g. set to <<1 photon
 
@@ -255,16 +263,19 @@ def colCommonModeCorrection(self, frame, arbitraryCut=1000):
             rowOffset = 0
             for b in range(0, self.detectorInfo.nBanksRow):
                 try:
-                    colCM = np.median(
-                        frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c][
-                            frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c] < arbitraryCut
-                        ]
-                    )
+                    testPixels = np.s_[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c]
+                    relevantPixels = frame[testPixels] < cut
+                    if switchedPixels is not None:
+                        ##print(testPixels, relevantPixels)
+                        relevantPixels = np.bitwise_and(relevantPixels, ~switchedPixels[testPixels])
+                    colCM = np.median(frame[testPixels][relevantPixels])
                     if not np.isnan(colCM):  ## if no pixels < cut we get nan
                         if False:
                             if c < 100:
                                 self.commonModeVals.append(colCM)
-                        frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c] -= colCM
+                        if colCM > cut:
+                            raise Exception("overcorrection: colCM, cut:", colCM, cut)
+                        frame[testPixels] -= colCM
                 except Exception:
                     print("colCM problem")
                     logger.error("colCM problem")

calibrationSuite/nonPsanaBase.py

@@ -0,0 +1,226 @@
+import os
+import sys
+
+import h5py
+import numpy as np
+
+## actually not psana, should rename
+## need this for DataSource.smalldata in analysis script
+from calibrationSuite.psana2Base import *  ##DataSource
+
+##from calibrationSuite.argumentParser import ArgumentParser
+##import importlib.util
+from calibrationSuite.psanaCommon import PsanaCommon
+
+## do we actually need this?
+
+
+class PsanaBase(PsanaCommon):  ## not actually psana, should rename
+    def __init__(self, analysisType="scan"):
+        super().__init__(analysisType)
+        commandUsed = sys.executable + " " + " ".join(sys.argv)
+        ##logger.info("Ran with cmd: " + commandUsed)
+
+        self.psanaType = 0
+        print("in nonPsanaBase")
+        ##logger.info("in nonPsanaBase")
+        self.myrun = None
+        self.setupNonPsana()
+
+    def setupNonPsana(self):
+        ## temp hack ##
+        ##self.args = ArgumentParser().parse_args()
+        ## now done in PsanaCommon
+        defaultConfigFileName = "suiteConfig.py"
+        secondaryConfigFileName = defaultConfigFileName if self.args.configFile is None else self.args.configFile
+        # secondaryConfigFileName is returned if env var not set
+        configFileName = os.environ.get("SUITE_CONFIG", secondaryConfigFileName)
+        config = self.importConfigFile(configFileName)
+        if config is None:
+            print("\ncould not find or read config file: " + configFileName)
+            print("please set SUITE_CONFIG env-var or use the '-cf' cmd-line arg to specify a valid config file")
+            print("exiting...")
+            sys.exit(1)
+        self.experimentHash = config.experimentHash
+        ## temp ##
+        nonXtcFiles = self.args.nonXtcFiles
+        if nonXtcFiles is not None:
+            self.file = None  ## nonPsanaBase only used as substitute xtc reader
+            print(nonXtcFiles)
+            self.setupNonXtcFiles(nonXtcFiles)
+
+        if "run" not in dir(self) or self.run is None:
+            self.run = 666666
+
+        self.ignoreEventCodeCheck = True
+
+    def setupNonXtcFiles(self, nonXtcFiles):
+        if nonXtcFiles.endswith("h5"):
+            self.h5File = h5py.File(nonXtcFiles)  ## need concat
+            self.getDataFromFile()
+            self.getStepInfoFromFile()
+            self.npyFile = False
+        elif nonXtcFiles.endswith("npy"):
+            if "," in nonXtcFiles:
+                self.data = np.concatenate([np.load(f) for f in nonXtcFiles.split(",")])
+            else:
+                self.data = np.load(nonXtcFiles)
+                fixAlexOrdering = True
+                if fixAlexOrdering:
+                    self.data = np.array([np.array([self.data[:, :, i]]) for i in range(self.data.shape[2])])
+
+            self.setDataArray(self.data)
+            self.h5File = False
+        else:
+            raise Exception("non-psana file formats h5 and npy only")
+
+    ##    def importConfigFile(self, file_path):
+    ##        if not os.path.exists(file_path):
+    ##            print(f"The file '{file_path}' does not exist")
+    ##            return None
+    ##        spec = importlib.util.spec_from_file_location("config", file_path)
+    ##        config_module = importlib.util.module_from_spec(spec)
+    ##        spec.loader.exec_module(config_module)
+    ##        return config_module
+
+    def getDataFromFile(self):
+        self.data = self.h5File["frames"][()]
+        self.setDataArray(self.data)
+
+    def getStepInfoFromFile(self):
+        try:
+            self.step_value = self.h5File["step_value"][()]
+        except:
+            print("could not find step_value in h5 file foo")
+            self.step_value = [None]
+
+    def setDataArray(self, data):
+        self.data = data
+        print("data shape:", self.data.shape)
+        self.dataIter = iter(self.data)
+        print("have set up dataIter")
+        self.myrun = MyRun(self.data)
+
+    def getStepGen(self):
+        self.stepIndex = -1  ## needed?
+        return self.stepGen()
+
+    def setupPsana(self):  ## should rename to setupBase or whatever
+        print("in fake setupPsana; get info from npy or h5 and fake stuff")
+        if False:
+            self.ds = DataSource(files=[nonXtcFiles])  ## to make smalldata happy
+            print("data source is", self.ds)
+        else:
+            self.ds = psana.DataSource(exp="rixx1005922", run=6)  ##66)
+            ##self.ds = psana.MPIDataSource('exp=detdaq23:run=1')##66)
+            ## get event to make smd call happy
+            ##evt = next(self.ds.events())
+        self.skipZeroCheck = True  ## some Alex data gets entirely lost else
+
+        if False:
+            self.getDataFromFile()
+            self.getStepInfoFromFile()
+
+    def stepGen(self):
+        while True:
+            try:
+                self.stepIndex += 1
+                step = Step(next(self.dataIter), self.stepIndex)
+                yield step
+            except:
+                ##raise StopIteration
+                return
+            ##print("in stepGen, stepIndex:", self.stepIndex)
+
+    def getScanValue(self, step, useStringInfo=False):
+        if useStringInfo:
+            ##print("no idea how to use stringInfo, giving up")
+            ##sys.exit(0)
+            pass
+        ##print(self.step_value[step.nStep])
+        ##print(type(self.step_value[step.nStep]))
+        return self.step_value[step.nStep]
+
+    def getEventCodes(self, foo):
+        return [283]
+
+    def plainGetRawData(self, evt):
+        return evt.astype("uint16")
+
+    def getFlux(self, foo):
+        return 1
+
+    def get_smalldata(self, filename):  ##, gather_interval=100):
+        try:
+            return self.ds.smalldata(filename=filename)  ##, gather_interval=gather_interval)
+        except Exception:
+            print("can't make smalldata - is datasource defined?")
+        return None
+
+
+class MyRun(object):
+    def __init__(self, array):
+        ## it would probably be smarter to have a reader (?) class instead
+        ## of loading everything into memory
+        ## could do
+        ## with h5py.File('TID.data') as data:
+        ##     frames = data['frames']
+        ##     eventN = 0
+        ##     while True:
+        ##         try:
+        ##             event = frames[eventN++, :...]
+        ##             yield event
+        ##         except:
+        ##             return
+        self.array = array
+
+    def events(self):
+        ## it would be smarter to yield the next event instead of filling RAM
+        ##return(np.vstack(self.array))
+        return self.array
+
+
+class Step(object):
+    def __init__(self, array, nStep):
+        self.array = array
+        self.nStep = nStep
+
+    def events(self):
+        return self.array
+
+
+if __name__ == "__main__":
+    npb = NonPsanaBase()
+
+    print("try running over 2, 3, 2 step data with step calls")
+    ##npb.setDataArray(np.load("../data/testStepData.npy"))
+    npb.setDataArray(np.load("testStepData.npy"))
+    stepGen = npb.getStepGen()
+    ##stepGen = npb.dataIter doesn't work because we want step.events below
+    for nStep, step in enumerate(stepGen):
+        print("step %d" % (nStep))
+        for nEvt, evt in enumerate(step.events()):
+            print("nEvt, evt:", nEvt, evt)
+
+    print("end up with stepIndex:", npb.stepIndex)
+
+    if True:
+        print("try running over 1, 3, 2 non-step data with step calls")
+        npb.setDataArray(np.load("../data/testNoStepData.npy"))
+        stepGen = npb.getStepGen()
+        for nStep, step in enumerate(stepGen):
+            print("step %d" % (nStep))
+            for nEvt, evt in enumerate(step.events()):
+                print(nEvt, evt)
+
+    print("try running over 2, 3, 2 non-step data with event call")
+    npb.setDataArray(np.load("../data/testStepData.npy"))
+    evtGen = npb.myrun.events()
+    for nEvt, evt in enumerate(evtGen):
+        print(nEvt, evt)
+
+    print("try running over 1, 3, 2 non-step data with event call")
+    npb.setDataArray(np.load("../data/testNoStepData.npy"))
+    evtGen = npb.myrun.events()
+    for nEvt, evt in enumerate(evtGen):
+        print(nEvt, evt)

calibrationSuite/psana2Base.py

@@ -81,7 +81,13 @@ def setupPsana(self):
             ##raise Exception
         ## could set to None and reset with first frame I guess, or does the det object know?
 
-        self.timing = self.myrun.Detector("timing")
+        try:
+            self.timing = self.myrun.Detector("timing")
+        except:
+            self.timing = None
+            ## timing makes us look at every single event
+            ## we only want this in MFX+RIX mode
+
         self.desiredCodes = {"120Hz": 272, "4kHz": 273, "5kHz": 274}
 
         try:
@@ -110,8 +116,16 @@ def get_ds(self, run=None):
         if run is None:
             run = self.run
         ##tmpDir = '/sdf/data/lcls/ds/rix/rixx1005922/scratch/xtc'
+        detectors = [self.experimentHash["detectorType"]]
+        if not self.ignoreEventCodeCheck:
+            detectors.append("timing")
+            detectors.append("MfxDg1BmMon")
         ds = psana.DataSource(
-            exp=self.exp, run=run, intg_det=self.experimentHash["detectorType"], max_events=self.maxNevents
+            exp=self.exp,
+            run=run,
+            intg_det=self.experimentHash["detectorType"],
+            max_events=self.maxNevents,
+            detectors=detectors,
         )  ##, dir=tmpDir)
         return ds
 
@@ -200,9 +214,13 @@ def getFlux(self, evt):
         return self.flux
 
     def getEventCodes(self, evt):
+        if self.timing is None:
+            return []
         return self.timing.raw.eventcodes(evt)
 
     def getPulseId(self, evt):
+        if self.timing is None:
+            return []
         return self.timing.raw.pulseId(evt)
 
     def isKicked(self, evt):

suite_scripts/AnalyzeH5.py

@@ -114,9 +114,11 @@ def getRowsColsFromSliceCoordinates(self):
     def analyze(self):
         if self.analysis == "cluster":
             self.clusterAnalysis()
+        if self.analysis == "linearity":
+            self.linearityAnalysis()
         else:
-            print("unknown analysis type %s" % (self.analysisType))
-            logging.info("unknown analysis type %s" % (self.analysisType))
+            print("unknown analysis type %s" % (self.analysis))
+            logging.info("unknown analysis type %s" % (self.analysis))
 
     def clusterAnalysis(self):
         clusters = None
@@ -306,6 +308,14 @@ def histogramAndFitGaussian(self, ax, energies, nBins):
         except Exception:
             return 0, 0, 0
 
+    def linearityAnalysis(self):
+        fluxes = np.concatenate([h5["fluxes"][()] for h5 in self.h5Files])
+        rois = np.concatenate([h5["rois"][()] for h5 in self.h5Files])
+        someSinglePixels = np.concatenate([h5["pixels"][()] for h5 in self.h5Files])
+        slice = np.concatenate([h5["slice"][()] for h5 in self.h5Files])
+        self.singlePixels = [[0, 66, 66], [0, 6, 6]]
+        print("load temporary fake pixels")
+
 
 if __name__ == "__main__":
     ah5 = AnalyzeH5()