Pure Python Histogramming (slower, but easier)

This script lives in /reg/g/psdm/tutorials/examplePython/npHist.py and demonstrates how to create and access data in a histogram, for the case where you can save all the histogram data in memory:

from psana import *
import numpy as np
ds = DataSource('exp=xpptut15:run=54:smd')
det = Detector('cspad')

pixel0_vals = []
for nevent,evt in enumerate(ds.events()):
    if nevent>=4: break
    calib_array = det.calib(evt)
    pixel0_vals.append(calib_array[0,0,0])
hist,bin_edges = np.histogram(pixel0_vals,bins=5,range=(-5.0,5.0))

import matplotlib.pyplot as plt
plt.bar(bin_edges[:-1],hist)
plt.show()

If you want to have the ability to "update" an existing histogram you can use this pattern which lives in /reg/g/psdm/tutorials/examplePython/hist.py:

from psana import *
from pypsalg.Histogram import hist1d, hist2d
ds = DataSource('exp=xpptut15:run=54:smd')
det = Detector('cspad',ds.env())

hist = hist1d(5,-5.,5.0)
for nevent,evt in enumerate(ds.events()):
    if nevent>=4: break
    calib_array = det.calib(evt)
    hist.fill(calib_array[0,0,0])

import matplotlib.pyplot as plt
plt.bar(hist.xaxis.values(),hist.data)
plt.show()

Python Histogramming via C++ (faster, but a little more complex)

These above "updating histogram" classes are written in python and can be slow.  We also have faster histogram classes written in C++ (callable from python) but they can be more difficult to use because they must be passed the correct types (python "floats", that can be created using python's "float()" function).  These classes can be used with the following:

from pypsalg import Hist1D, Hist2D

The main functions are "fill" (fills histograms with data), "get" (returns histogram contents as array), and "xaxis"/"yaxis".  Help for all these can be seen in ipython with commands like:

Hist1D.fill?

"fill" has several different options for filling individual numbers, arrays of numbers, and varying weights.

An example script is in  /reg/g/psdm/tutorials/examplePython/hist_c.py:

from psana import *
from pypsalg import Hist1D
ds = DataSource('exp=xpptut15:run=54:smd')
det = Detector('cspad',ds.env())
hist = Hist1D(5,-5.,5.0)
for nevent,evt in enumerate(ds.events()):
    if nevent>=4: break
    calib_array = det.calib(evt)
    hist.fill(float(calib_array[0,0,0]),1.0)
import matplotlib.pyplot as plt
plt.bar(hist.xaxis(),hist.get())
plt.show()
  • No labels