Skip to content

Commit

Permalink
STY: Apply black and ruff again after rebase
Browse files Browse the repository at this point in the history
  • Loading branch information
nstelter-slac committed Jan 30, 2024
1 parent 526c211 commit 464791e
Show file tree
Hide file tree
Showing 7 changed files with 187 additions and 86 deletions.
41 changes: 30 additions & 11 deletions scripts/AnalyzeH5.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,14 @@ def __init__(self, outputDir, className, run, camera, label):
self.camera = camera
self.label = label


# Setup logging.
# Log file gets appended to each new run, can manually delete for fresh log.
# Could change so makes new unique log each run or overwrites existing log.
logging.basicConfig(
filename='analyze_h5.log',
level=logging.INFO, # For full logging set to INFO which includes ERROR logging too
format='%(asctime)s - %(levelname)s - %(message)s' # levelname is log severity (ERROR, INFO, etc)
filename="analyze_h5.log",
level=logging.INFO, # For full logging set to INFO which includes ERROR logging too
format="%(asctime)s - %(levelname)s - %(message)s", # levelname is log severity (ERROR, INFO, etc)
)


Expand Down Expand Up @@ -61,9 +62,15 @@ def __init__(self):
self.sliceEdges = args.slice_edges.split(",")
self.sliceEdges = [int(curr) for curr in self.sliceEdges]
self.nBins = 100
self.shiftEnergy = False if args.shift_energy_bits == None else True
self.analysisNum = 1 if args.analysis_mode == None else int(args.analysis_mode)
self.fileNameInfo = FileNamingInfo(args.path, self.__class__.__name__, args.run, 0, args.label,)
self.shiftEnergy = False if args.shift_energy_bits is None else True
self.analysisNum = 1 if args.analysis_mode is None else int(args.analysis_mode)
self.fileNameInfo = FileNamingInfo(
args.path,
self.__class__.__name__,
args.run,
0,
args.label,
)
print("Output dir: " + self.fileNameInfo.outputDir)
logging.info("Output dir: " + self.fileNameInfo.outputDir)

Expand Down Expand Up @@ -203,24 +210,36 @@ def plot_gain_distribution(self, fitInfo, fileInfo, fitIndex):
plt.savefig(fileName)

def analyzeSimpleClusters(self, clusters):
energy = clusters[:, :, 0] #.flatten()
energy = clusters[:, :, 0] # .flatten()
if self.shiftEnergy:
energy *= 2 # temporary, due to bit shift
rows = self.sliceEdges[0]
cols = self.sliceEdges[1]
fitInfo = np.zeros((rows, cols, 4)) # mean, std, mu, sigma
fitInfo = np.zeros((rows, cols, 4)) # mean, std, mu, sigma
fitIndex = 0

self.plot_overall_energy_distribution(energy, self.fileNameInfo)

print("Analysis Mode: " + str(self.analysisNum))
logging.info("Analysis Mode: " + str(self.analysisNum))
if self.analysisNum == 1:
fitIndex = 2
fitInfo = pixelAnalysis.analysis_one(clusters, energy, rows, cols, fitInfo, self.lowEnergyCut, self.highEnergyCut, self.fileNameInfo)
fitInfo = pixelAnalysis.analysis_one(
clusters, energy, rows, cols, fitInfo, self.lowEnergyCut, self.highEnergyCut, self.fileNameInfo
)
else:
fitIndex = 3
fitInfo = pixelAnalysis.analysis_two(clusters, self.nBins, self.sliceCoordinates, rows, cols, fitInfo, self.lowEnergyCut, self.highEnergyCut, self.fileNameInfo)
fitInfo = pixelAnalysis.analysis_two(
clusters,
self.nBins,
self.sliceCoordinates,
rows,
cols,
fitInfo,
self.lowEnergyCut,
self.highEnergyCut,
self.fileNameInfo,
)

self.save_fit_information(fitInfo, rows, cols, self.fileNameInfo)
self.plot_gain_distribution(fitInfo, self.fileNameInfo, fitIndex)
Expand Down
26 changes: 15 additions & 11 deletions scripts/ancillaryMethods.py
Original file line number Diff line number Diff line change
@@ -1,43 +1,47 @@
import numpy as np
from scipy.stats import binned_statistic


def makeProfile(x, y, bins, range=None, spread=False):
## NaN for empty bins are suppressed
## using mean root(N) for non-empty bins to calculate 0 var weights
##
## spread=True to return standard deviation instead of standard error

meansObj = binned_statistic(x, [y, y**2], bins=bins, range=range, statistic='mean')
meansObj = binned_statistic(x, [y, y**2], bins=bins, range=range, statistic="mean")
means, means2 = meansObj.statistic
countsObj = binned_statistic(x, [y, y**2], bins=bins, range=(0,1), statistic='count')
countsObj = binned_statistic(x, [y, y**2], bins=bins, range=(0, 1), statistic="count")
bin_N = countsObj.statistic
yErr = np.sqrt(means2 - means**2)
if not spread:
root_N = np.sqrt(bin_N)
root_N[root_N==0] = root_N[root_N>0].mean()
yErr = yErr/root_N
root_N[root_N == 0] = root_N[root_N > 0].mean()
yErr = yErr / root_N
##yErr = yErr.clip(0, 6666666.)
bin_edges = means_result.bin_edges
bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.
usefulBins = bin_N>0
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0
usefulBins = bin_N > 0
return bin_centers[usefulBins], means[usefulBins], yErr[usefulBins]


def plotProfile(x, y, yErr):
plt.errorbar(x=x, y=y, yerr=yErr, linestyle='none', marker='.')
plt.errorbar(x=x, y=y, yerr=yErr, linestyle="none", marker=".")


def selectedClusters(clusters, row, col, lowEnerygCut, highEnergyCut, nPixelCut=4, isSquare=1):
pass


def goodClusters(clusters, row, col, nPixelCut=4, isSquare=None):
##print(clusters)
pixelRowCol = np.bitwise_and((clusters[:,:,1] == row),
(clusters[:,:,2] == col))
pixelRowCol = np.bitwise_and((clusters[:, :, 1] == row), (clusters[:, :, 2] == col))
if isSquare is None:
small = clusters[:,:,3]<nPixelCut
small = clusters[:, :, 3] < nPixelCut
else:
small = np.bitwise_and((clusters[:,:,3]<nPixelCut), (clusters[:,:,4]==isSquare))
small = np.bitwise_and((clusters[:, :, 3] < nPixelCut), (clusters[:, :, 4] == isSquare))
return clusters[np.bitwise_and(small, pixelRowCol)]


def getClusterEnergies(clusters):
##print(clusters)
return clusters[:, 0]
59 changes: 38 additions & 21 deletions scripts/fitFunctions.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,46 +6,57 @@
def linear(x, a, b):
return a * x + b


def saturatedLinear(x, a, b, c, d):
return [a * X + b if X < c else d for X in x]


def saturatedLinearB(x, a, b, d):
return [a * X + b if a * X + b < d else d for X in x]


def gaussian(x, a, mu, sigma):
return a * np.exp(-((x - mu) ** 2) / (2 * sigma**2))


def gaussianArea(a, sigma):
return a * sigma * 6.28


def estimateGaussianParameters(flatData):
return flatData.mean(), flatData.std()


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 * 6.28), 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 * 6.28), mean, sigma


def getHistogramMeanStd(binCenters, counts):
mean = np.average(binCenters, weights=counts)
var = np.average((binCenters - mean) ** 2, weights=counts)
return mean, np.sqrt(var)


def getBinCentersFromNumpyHistogram(bins):
return (bins[1:] + bins[:-1]) / 2.0


def getRestrictedHistogram(bins, counts, x0, x1):
cut = np.where(np.logical_and(bins >= x0, bins <= x1))
x = bins[cut]
y = counts[cut]
return x, y


# unused? and won't work in current state i think...
"""
def getGaussianFitFromHistogram(binCenters, counts, x0=None, x1=None):
Expand All @@ -71,19 +82,22 @@ def fitNorm(data):
mean, std = norm.fit(data)
return mean, std


def twoGaussSilvermanModeTest(x0, x1):
a = np.random.normal(0, 1, 1000)
b = np.random.normal(x0, 1, 500)
c = np.random.normal(x1, 1, 500)
d = np.append(b, c)
if False:
import matplotlib.pyplot as plt

plt.hist(d, 100)
##plt.hist(a, 100)
plt.show()
print(a.mean(), a.std(), d.mean(), d.std())
print(x0, x1, bw_silverman(a), bw_silverman(d))
return d.std()/a.std(), bw_silverman(d)/bw_silverman(a)
return d.std() / a.std(), bw_silverman(d) / bw_silverman(a)


def testSilvermanModeTest():
a = np.linspace(0, 3, 10)
Expand All @@ -95,66 +109,69 @@ def testSilvermanModeTest():
S.append(r1)

import matplotlib.pyplot as plt

plt.plot(noise, S)
plt.title("Two gaussian test of Silverman's test separating peaks")
plt.xlabel("two peak rms/one peak rms")
plt.ylabel("two peak Silverman/one peak Silverman")
plt.show()

a = np.random.normal(0, 100, 1000)
b = np.random.normal(0, 100, 10000)
c = np.where(b!=50)
c = np.where(b != 50)
print("rms for a gap in gaussian, notched gaussian:", a.std(), b[c][:1000].std())
print("Silverman for a gap in gaussian, notched gaussian:", bw_silverman(a), bw_silverman(b[c][:1000]))
print("rms for a gap in gaussian/notched gaussian:", a.std()/b[c][:1000].std())
print("Silverman for a gap in gaussian/notched gaussian:", bw_silverman(a)/bw_silverman(b[c][:1000]))
print("rms for a gap in gaussian/notched gaussian:", a.std() / b[c][:1000].std())
print("Silverman for a gap in gaussian/notched gaussian:", bw_silverman(a) / bw_silverman(b[c][:1000]))



def missingBinTest(binCenters, counts):
mu, sigma = getHistogramMeanStd(binCenters, counts)
binCenters, counts = getRestrictedHistogram(binCenters, counts, mu-sigma, mu+sigma)
binCenters, counts = getRestrictedHistogram(binCenters, counts, mu - sigma, mu + sigma)
##print(len(b), len(c))
n = len(counts)
if n>=10:
step = int(n/10)
if n >= 10:
step = int(n / 10)
else:
step = n

cuts = range(1, 6)
missingBins = np.zeros(len(cuts))

for i in range(step):
i0 = step*i
i1 = min(step*(i+1), n)
i0 = step * i
i1 = min(step * (i + 1), n)
med = np.median(counts[i0:i1])
mean = counts[i0:i1].mean()
rootN = np.sqrt(mean)
for k in cuts:
for j in range(i0, i1):
if counts[j]<med-k*rootN:
if counts[j] < med - k * rootN:
missingBins[k] += 1
##print(n, step, i, i0, i1, k, j, binCenters[j], counts[j], med-k*rootN)
return missingBins


def testMissingBinTest():
nSamples = 5000
a = np.random.normal(0, 100, nSamples)
b = np.random.normal(0, 100, nSamples*2)
b = np.random.normal(0, 100, nSamples * 2)
missingValue = 30
c = np.where(np.logical_or(b>(missingValue+1), b<missingValue))
c = np.where(np.logical_or(b > (missingValue + 1), b < missingValue))
b = b[c][:nSamples]
##print(50 in a, 50 in b)
ha, bins = np.histogram(a, 600, [-300, 300])
hb,_ = np.histogram(b, 600, [-300, 300])
hb, _ = np.histogram(b, 600, [-300, 300])
binCenters = getBinCentersFromNumpyHistogram(bins)
##print(binCenters, ha, hb)
if True:
import matplotlib.pyplot as plt

plt.stairs(ha, bins)
plt.stairs(hb, bins, color='b')
plt.stairs(hb, bins, color="b")
plt.show()

mbt_a = missingBinTest(binCenters, ha)
mbt_b = missingBinTest(binCenters, hb)
print("per n sigma check for gap in gaussian, notched gaussian:", mbt_a, mbt_b)
print((range(1,6)*mbt_a).sum(), (range(1,6)*mbt_b).sum())
print((range(1, 6) * mbt_a).sum(), (range(1, 6) * mbt_b).sum())
Loading

0 comments on commit 464791e

Please sign in to comment.