Skip to content

Commit

Permalink
Merge pull request #141 from slaclab/analyzeNpyScan
Browse files Browse the repository at this point in the history
Merge analyzeNpyScan into deveopment
  • Loading branch information
nstelter-slac authored Oct 29, 2024
2 parents 45d1694 + b4dd4df commit edf8c52
Show file tree
Hide file tree
Showing 10 changed files with 357 additions and 44 deletions.
3 changes: 3 additions & 0 deletions calibrationSuite/argumentParser.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
33 changes: 22 additions & 11 deletions calibrationSuite/basicSuiteScript.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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):
Expand All @@ -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

Expand All @@ -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")
Expand Down
226 changes: 226 additions & 0 deletions calibrationSuite/nonPsanaBase.py
Original file line number Diff line number Diff line change
@@ -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)
22 changes: 20 additions & 2 deletions calibrationSuite/psana2Base.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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):
Expand Down
14 changes: 12 additions & 2 deletions suite_scripts/AnalyzeH5.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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()
Expand Down
Loading

0 comments on commit edf8c52

Please sign in to comment.