From f840578f1706d4d75f1d3d31fa57c1974e2fee13 Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Thu, 24 Oct 2024 16:07:57 -0700 Subject: [PATCH 1/2] add non-psana class and support; smarter common mode correction --- calibrationSuite/argumentParser.py | 3 + calibrationSuite/basicSuiteScript.py | 35 ++- calibrationSuite/nonPsanaBase.py | 223 +++++++++++++++++++ calibrationSuite/psana2Base.py | 19 +- suite_scripts/AnalyzeH5.py | 11 + suite_scripts/EventScanParallelSlice.py | 8 +- suite_scripts/LinearityPlotsParallelSlice.py | 90 ++++++-- suite_scripts/SimpleClustersParallelSlice.py | 11 +- suite_scripts/simpleTestScript.py | 3 +- 9 files changed, 359 insertions(+), 44 deletions(-) create mode 100644 calibrationSuite/nonPsanaBase.py diff --git a/calibrationSuite/argumentParser.py b/calibrationSuite/argumentParser.py index d748545..879a127 100644 --- a/calibrationSuite/argumentParser.py +++ b/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") diff --git a/calibrationSuite/basicSuiteScript.py b/calibrationSuite/basicSuiteScript.py index 8c0e012..d3a0deb 100755 --- a/calibrationSuite/basicSuiteScript.py +++ b/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 @@ -254,18 +262,23 @@ def colCommonModeCorrection(self, frame, arbitraryCut=1000): for c in range(self.detectorInfo.nCols): rowOffset = 0 for b in range(0, self.detectorInfo.nBanksRow): - try: + if True:##try: + testPixels = np.s_[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c] + relevantPixels = frame[testPixels] cut: + raise Exception("overcorrection: colCM, cut:", colCM, cut) + frame[testPixels] -= colCM + else:##except Exception: print("colCM problem") logger.error("colCM problem") print(frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank], c) diff --git a/calibrationSuite/nonPsanaBase.py b/calibrationSuite/nonPsanaBase.py new file mode 100644 index 0000000..613b993 --- /dev/null +++ b/calibrationSuite/nonPsanaBase.py @@ -0,0 +1,223 @@ +import sys, os +import numpy as np +import h5py +##from calibrationSuite.argumentParser import ArgumentParser +##import importlib.util +from calibrationSuite.psanaCommon import PsanaCommon +## actually not psana, should rename + + +## need this for DataSource.smalldata in analysis script +from calibrationSuite.psana2Base import *##DataSource +## 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) diff --git a/calibrationSuite/psana2Base.py b/calibrationSuite/psana2Base.py index 09112f0..e642e3b 100755 --- a/calibrationSuite/psana2Base.py +++ b/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,13 @@ 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 +211,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): diff --git a/suite_scripts/AnalyzeH5.py b/suite_scripts/AnalyzeH5.py index a04cea5..06d4439 100644 --- a/suite_scripts/AnalyzeH5.py +++ b/suite_scripts/AnalyzeH5.py @@ -114,6 +114,8 @@ 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)) @@ -307,6 +309,15 @@ def histogramAndFitGaussian(self, ax, energies, nBins): 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() ah5.getFiles() diff --git a/suite_scripts/EventScanParallelSlice.py b/suite_scripts/EventScanParallelSlice.py index b9204f1..3064245 100644 --- a/suite_scripts/EventScanParallelSlice.py +++ b/suite_scripts/EventScanParallelSlice.py @@ -273,12 +273,12 @@ def analyze_h5(self, dataFile, label): ##print(frames[tuple(esp.singlePixels[2])]) if esp.special is not None and "rowCommonMode" in esp.special: - frames = np.array([esp.rowCommonModeCorrection3d(frames)]) + frames = esp.rowCommonModeCorrection3d(frames, 3.0) if esp.special is not None and "colCommonMode" in esp.special: - frames = np.array([esp.colCommonModeCorrection3d(frames)]) + frames = esp.colCommonModeCorrection3d(frames, 3.0) if esp.special is not None and "regionCommonMode" in esp.special: ##oldFrames = frames - frames = np.array([esp.regionCommonModeCorrection(frames, esp.regionSlice, 666)]) + frames = esp.regionCommonModeCorrection(frames, esp.regionSlice, 3.0) ##print(frames-oldFrames) else: if esp.special is not None and "regionCommonMode" in esp.special: @@ -286,7 +286,7 @@ def analyze_h5(self, dataFile, label): eventNumbers.append(nevt) for i, roi in enumerate(esp.ROIs): - m = frames[roi == 1].mean() + m = frames[roi > 0].mean() roiMeans[i].append(m) for i, roi in enumerate(esp.singlePixels): diff --git a/suite_scripts/LinearityPlotsParallelSlice.py b/suite_scripts/LinearityPlotsParallelSlice.py index 4b223f2..2df4fd1 100755 --- a/suite_scripts/LinearityPlotsParallelSlice.py +++ b/suite_scripts/LinearityPlotsParallelSlice.py @@ -10,6 +10,7 @@ import logging import os import sys +import copy import h5py import matplotlib.pyplot as plt @@ -374,15 +375,23 @@ def analyze_h5_slice(self, dataFile, label): print("not doing Kaz events") logger.info("not doing Kaz events") + commonModeCorrectionForZongde = False + if lpp.special is not None and "colCommonMode" in lpp.special: + commonModeCorrectionForZongde = True + lpp.setupPsana() try: size = comm.Get_size() # noqa: F821 except Exception: size = 1 - smd = lpp.ds.smalldata( - filename="%s/%s_%s_c%d_r%d_n%d.h5" % (lpp.outputDir, lpp.className, lpp.label, lpp.camera, lpp.run, size) - ) + + filename="%s/%s_%s_c%d_r%d_n%d.h5" % (lpp.outputDir, lpp.className, lpp.label, lpp.camera, lpp.run, size) + if lpp.psanaType == 1: + smd = lpp.ds.small_data(filename, gather_interval=100) + else: + smd = lpp.get_smalldata(filename) + nGoodEvents = 0 fluxes = [] ## common for all ROIs @@ -400,7 +409,8 @@ def analyze_h5_slice(self, dataFile, label): if doFast: ec = lpp.getEventCodes(evt) beamEvent = lpp.isBeamEvent(evt) - if ec[137]: + ##print(ec, beamEvent) + if not ec == [] and ec[137]: lpp.flux = lpp._getFlux(evt) ## fix this lpp.fluxTS = lpp.getTimestamp(evt) ##print("found flux", lpp.flux) @@ -409,10 +419,18 @@ def analyze_h5_slice(self, dataFile, label): ##frames = lpp.det.raw(evt)&0x3fff ##elif ec[281]: elif beamEvent: - lpp.framesTS = lpp.getTimestamp(evt) + try: + lpp.framesTS = lpp.getTimestamp(evt) + except: ## for nonPsana + lpp.framesTS = nevt + lpp.fluxTS = nevt rawFrames = lpp.getRawData(evt, gainBitsMasked=False) - frames = lpp.getCalibData(evt) + try: + frames = lpp.getCalibData(evt) + except: + frames = rawFrames ## for nonPsana stuff else: + print("something problematic with event code going on") continue else: lpp.flux = lpp._getFlux(evt) ## fix this @@ -451,12 +469,30 @@ def analyze_h5_slice(self, dataFile, label): ##print('nSwitched: %d' %(nSwitched)) continue + if commonModeCorrectionForZongde: + ## find switched pixels, strip gain bits,## make a copy + ## subtract pedestal from non-switched + ## (later make this apply to all, probably) + ## do common mode correction + ## cast back to int in case it matters + switchedPixels = lpp.getSwitchedPixels(rawFrames) + rawFrames = rawFrames & lpp.gainBitsMask + ##no_cm_rawFrames = copy.copy(rawFrames) + rawFrames = rawFrames.astype("float") + rawFrames[~switchedPixels] -= lpp.fakePedestal[~switchedPixels] + rawFrames = lpp.colCommonModeCorrection3d(rawFrames, cut=1000, switchedPixels=switchedPixels) + rawFrames = rawFrames.round().astype("int32") + roiMeans = [] for i, roi in enumerate(lpp.ROIs): ##m = np.multiply(roi, frames).mean() - m = frames[roi == 1].mean() - roiMeans.append(m) - + try: + m = frames[roi == 1].mean() + roiMeans.append(m) + except: + ## probably have a bad config pointing to inapposite ROIs + pass + if doKazFlux: rf = rawFrames[tuple(lpp.singlePixels[0])] if not (flux < 20000 and rf >= lpp.g0cut and rf > 2950): @@ -468,23 +504,34 @@ def analyze_h5_slice(self, dataFile, label): singlePixelData = [] for j, p in enumerate(lpp.singlePixels): - singlePixelData.append([int(rawFrames[tuple(p)] >= lpp.g0cut), rawFrames[tuple(p)] & lpp.gainBitsMask]) - + if commonModeCorrectionForZongde: + singlePixelData.append([switchedPixels[tuple(p)], rawFrames[tuple(p)]]) + else: + singlePixelData.append([int(rawFrames[tuple(p)] >= lpp.g0cut), rawFrames[tuple(p)] & lpp.gainBitsMask]) + eventDict = { "fluxes": flux, - "rois": np.array(roiMeans), + "rois": np.array(roiMeans), "pixels": np.array(singlePixelData), "slice": rawFrames[lpp.regionSlice], } - smd.event( - evt, - eventDict, - fluxes=flux, - rois=np.array(roiMeans), - pixels=np.array(singlePixelData), - slice=rawFrames[lpp.regionSlice], - ) + if lpp.psanaType == 0: + eventDict['fluxes'] = nevt + smd.event( + nevt, + eventDict) + elif lpp.psanaType == 2: + smd.event( + evt, + eventDict, + ##fluxes=flux, + ##rois=np.array(roiMeans), + ##pixels=np.array(singlePixelData), + ##slice=rawFrames[lpp.regionSlice], + ) + else: + smd.event(eventDict) nGoodEvents += 1 if nGoodEvents % 100 == 0: @@ -569,6 +616,7 @@ def analyze_h5_slice(self, dataFile, label): lpp.plotDataROIs(roiMeans, fluxes, "roiTest") """ - smd.done() + if lpp.psanaType == 2: + smd.done() lpp.dumpEventCodeStatistics() diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index 321875a..56a4d56 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -139,10 +139,11 @@ def analyze_h5(self, dataFile, label): print("analyzed modules:", sic.analyzedModules) ## move this to psana setup size = 666 filename = "%s/%s_%s_c%d_r%d_n%d.h5" % (sic.outputDir, sic.className, sic.label, sic.camera, sic.run, size) - if sic.psanaType == 1: - smd = sic.ds.small_data(filename=filename, gather_interval=100) - else: + if sic.psanaType == 2: smd = sic.get_smalldata(filename=filename) + else: + smd = sic.ds.small_data(filename=filename, gather_interval=100) + ## 50x50 pixels, 3x3 clusters, 10% occ., 2 sensors maxClusters = 10000 ##int(50 * 50 / 3 / 3 * 0.1 * 2) @@ -359,11 +360,11 @@ def analyze_h5(self, dataFile, label): ## np.save("%s/eventNumbers_c%d_r%d_%s.npy" %(sic.outputDir, sic.camera, sic.run, sic.exp), np.array(eventNumbers)) ## sic.plotData(roiMeans, pixelValues, eventNumbers, "foo") - if sic.psanaType == 1 or smd.summary: + if sic.psanaType < 2 or smd.summary: ## guess at desired psana1 behavior - no smd.summary there ## maybe check smd.rank == 0? sumhSum = smd.sum(hSum) - if sic.psanaType == 1: + if sic.psanaType < 2: smd.save({"energyHistogram": sumhSum}) smd.save(sic.configHash) else: diff --git a/suite_scripts/simpleTestScript.py b/suite_scripts/simpleTestScript.py index 5967c8e..6977042 100644 --- a/suite_scripts/simpleTestScript.py +++ b/suite_scripts/simpleTestScript.py @@ -9,7 +9,7 @@ ds = DataSource(exp=exp, run=run, intg_det="archon", max_events=666666) # noqa: F405 smd = ds.smalldata(filename="foo.h5") -print(help(smd.event)) +##print(help(smd.event)) myrun = next(ds.runs()) det = myrun.Detector("archon") @@ -29,4 +29,5 @@ nGood += 1 + print("good vs evil:", nGood, nBad) From b4dd4dfad391b37dc21093f8126092432580aabd Mon Sep 17 00:00:00 2001 From: philiph-slac Date: Thu, 24 Oct 2024 17:14:11 -0700 Subject: [PATCH 2/2] lint changes --- calibrationSuite/basicSuiteScript.py | 16 ++-- calibrationSuite/nonPsanaBase.py | 93 ++++++++++---------- calibrationSuite/psana2Base.py | 11 ++- suite_scripts/AnalyzeH5.py | 21 +++-- suite_scripts/LinearityPlotsParallelSlice.py | 21 ++--- suite_scripts/SimpleClustersParallelSlice.py | 1 - suite_scripts/simpleTestScript.py | 1 - suite_scripts/slurmSubModulesAndRegions.py | 2 +- 8 files changed, 82 insertions(+), 84 deletions(-) diff --git a/calibrationSuite/basicSuiteScript.py b/calibrationSuite/basicSuiteScript.py index d3a0deb..db77c0a 100755 --- a/calibrationSuite/basicSuiteScript.py +++ b/calibrationSuite/basicSuiteScript.py @@ -23,7 +23,7 @@ if os.getenv("foo") == "0": print("non-psana") - from calibrationSuite.nonPsanaBase import PsanaBase + from calibrationSuite.nonPsanaBase import PsanaBase elif os.getenv("foo") == "1": print("psana1") from calibrationSuite.psana1Base import PsanaBase @@ -100,7 +100,7 @@ def getRawData(self, evt, gainBitsMasked=True, negativeGain=False): return None nZero = 0 - if not 'skipZeroCheck' in dir(self): + if not "skipZeroCheck" in dir(self): nZero = frames.size - np.count_nonzero(frames) try: dz = self.nZero - nZero @@ -253,7 +253,7 @@ def rowCommonModeCorrection(self, frame, arbitraryCut=1000): colOffset += self.detectorInfo.nColsPerBank return frame - def colCommonModeCorrection(self, frame, cut=1000, switchedPixels = None): + 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 @@ -262,15 +262,13 @@ def colCommonModeCorrection(self, frame, cut=1000, switchedPixels = None): for c in range(self.detectorInfo.nCols): rowOffset = 0 for b in range(0, self.detectorInfo.nBanksRow): - if True:##try: + try: testPixels = np.s_[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank, c] - relevantPixels = frame[testPixels] cut: raise Exception("overcorrection: colCM, cut:", colCM, cut) frame[testPixels] -= colCM - else:##except Exception: + except Exception: print("colCM problem") logger.error("colCM problem") print(frame[rowOffset : rowOffset + self.detectorInfo.nRowsPerBank], c) diff --git a/calibrationSuite/nonPsanaBase.py b/calibrationSuite/nonPsanaBase.py index 613b993..cb83ec2 100644 --- a/calibrationSuite/nonPsanaBase.py +++ b/calibrationSuite/nonPsanaBase.py @@ -1,18 +1,22 @@ -import sys, os -import numpy as np +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 -## actually not psana, should rename - -## need this for DataSource.smalldata in analysis script -from calibrationSuite.psana2Base import *##DataSource ## do we actually need this? -class PsanaBase(PsanaCommon): ## not actually psana, should rename - def __init__(self, analysisType='scan'): + +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) @@ -41,49 +45,46 @@ def setupNonPsana(self): ## temp ## nonXtcFiles = self.args.nonXtcFiles if nonXtcFiles is not None: - self.file = None ## nonPsanaBase only used as substitute xtc reader + 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 + 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 + 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(',')]) + 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.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 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.data = self.h5File["frames"][()] self.setDataArray(self.data) def getStepInfoFromFile(self): @@ -99,27 +100,27 @@ def setDataArray(self, data): self.dataIter = iter(self.data) print("have set up dataIter") self.myrun = MyRun(self.data) - + def getStepGen(self): - self.stepIndex = -1 ## needed? + self.stepIndex = -1 ## needed? return self.stepGen() - def setupPsana(self): ## should rename to setupBase or whatever + 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 + 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.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 - + self.skipZeroCheck = True ## some Alex data gets entirely lost else + if False: self.getDataFromFile() self.getStepInfoFromFile() - + def stepGen(self): while True: try: @@ -144,7 +145,7 @@ def getEventCodes(self, foo): return [283] def plainGetRawData(self, evt): - return evt.astype('uint16') + return evt.astype("uint16") def getFlux(self, foo): return 1 @@ -176,15 +177,17 @@ def __init__(self, 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) - + return self.array + + class Step(object): def __init__(self, array, nStep): self.array = array self.nStep = nStep def events(self): - return(self.array) + return self.array + if __name__ == "__main__": npb = NonPsanaBase() @@ -195,18 +198,18 @@ def events(self): 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)) + 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)) + print("step %d" % (nStep)) for nEvt, evt in enumerate(step.events()): print(nEvt, evt) diff --git a/calibrationSuite/psana2Base.py b/calibrationSuite/psana2Base.py index e642e3b..9bc5adc 100755 --- a/calibrationSuite/psana2Base.py +++ b/calibrationSuite/psana2Base.py @@ -87,7 +87,7 @@ def setupPsana(self): 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: @@ -116,13 +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"]] + 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, - detectors=detectors + exp=self.exp, + run=run, + intg_det=self.experimentHash["detectorType"], + max_events=self.maxNevents, + detectors=detectors, ) ##, dir=tmpDir) return ds diff --git a/suite_scripts/AnalyzeH5.py b/suite_scripts/AnalyzeH5.py index 06d4439..2a7b445 100644 --- a/suite_scripts/AnalyzeH5.py +++ b/suite_scripts/AnalyzeH5.py @@ -117,8 +117,8 @@ def analyze(self): 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 @@ -308,16 +308,15 @@ 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") - - + 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() ah5.getFiles() diff --git a/suite_scripts/LinearityPlotsParallelSlice.py b/suite_scripts/LinearityPlotsParallelSlice.py index 2df4fd1..3b57618 100755 --- a/suite_scripts/LinearityPlotsParallelSlice.py +++ b/suite_scripts/LinearityPlotsParallelSlice.py @@ -7,10 +7,10 @@ ## may be copied, modified, propagated, or distributed except according to ## the terms contained in the LICENSE.txt file. ############################################################################## +import copy import logging import os import sys -import copy import h5py import matplotlib.pyplot as plt @@ -386,13 +386,12 @@ def analyze_h5_slice(self, dataFile, label): except Exception: size = 1 - filename="%s/%s_%s_c%d_r%d_n%d.h5" % (lpp.outputDir, lpp.className, lpp.label, lpp.camera, lpp.run, size) + filename = "%s/%s_%s_c%d_r%d_n%d.h5" % (lpp.outputDir, lpp.className, lpp.label, lpp.camera, lpp.run, size) if lpp.psanaType == 1: smd = lpp.ds.small_data(filename, gather_interval=100) else: smd = lpp.get_smalldata(filename) - nGoodEvents = 0 fluxes = [] ## common for all ROIs roiMeans = [[] for i in lpp.ROIs] @@ -421,14 +420,14 @@ def analyze_h5_slice(self, dataFile, label): elif beamEvent: try: lpp.framesTS = lpp.getTimestamp(evt) - except: ## for nonPsana + except: ## for nonPsana lpp.framesTS = nevt lpp.fluxTS = nevt rawFrames = lpp.getRawData(evt, gainBitsMasked=False) try: frames = lpp.getCalibData(evt) except: - frames = rawFrames ## for nonPsana stuff + frames = rawFrames ## for nonPsana stuff else: print("something problematic with event code going on") continue @@ -492,7 +491,7 @@ def analyze_h5_slice(self, dataFile, label): except: ## probably have a bad config pointing to inapposite ROIs pass - + if doKazFlux: rf = rawFrames[tuple(lpp.singlePixels[0])] if not (flux < 20000 and rf >= lpp.g0cut and rf > 2950): @@ -508,19 +507,17 @@ def analyze_h5_slice(self, dataFile, label): singlePixelData.append([switchedPixels[tuple(p)], rawFrames[tuple(p)]]) else: singlePixelData.append([int(rawFrames[tuple(p)] >= lpp.g0cut), rawFrames[tuple(p)] & lpp.gainBitsMask]) - + eventDict = { "fluxes": flux, - "rois": np.array(roiMeans), + "rois": np.array(roiMeans), "pixels": np.array(singlePixelData), "slice": rawFrames[lpp.regionSlice], } if lpp.psanaType == 0: - eventDict['fluxes'] = nevt - smd.event( - nevt, - eventDict) + eventDict["fluxes"] = nevt + smd.event(nevt, eventDict) elif lpp.psanaType == 2: smd.event( evt, diff --git a/suite_scripts/SimpleClustersParallelSlice.py b/suite_scripts/SimpleClustersParallelSlice.py index 56a4d56..ce3941a 100644 --- a/suite_scripts/SimpleClustersParallelSlice.py +++ b/suite_scripts/SimpleClustersParallelSlice.py @@ -144,7 +144,6 @@ def analyze_h5(self, dataFile, label): else: smd = sic.ds.small_data(filename=filename, gather_interval=100) - ## 50x50 pixels, 3x3 clusters, 10% occ., 2 sensors maxClusters = 10000 ##int(50 * 50 / 3 / 3 * 0.1 * 2) if sic.seedCut is not None: diff --git a/suite_scripts/simpleTestScript.py b/suite_scripts/simpleTestScript.py index 6977042..4bac51f 100644 --- a/suite_scripts/simpleTestScript.py +++ b/suite_scripts/simpleTestScript.py @@ -29,5 +29,4 @@ nGood += 1 - print("good vs evil:", nGood, nBad) diff --git a/suite_scripts/slurmSubModulesAndRegions.py b/suite_scripts/slurmSubModulesAndRegions.py index 85c6707..83a836b 100644 --- a/suite_scripts/slurmSubModulesAndRegions.py +++ b/suite_scripts/slurmSubModulesAndRegions.py @@ -1,5 +1,5 @@ -from subprocess import call import time +from subprocess import call from calibrationSuite.detectorInfo import DetectorInfo