chuckie82 / pysingfel Goto Github PK

For the Epix10k detector, the calibration files (e.g. pedestal, pixel_status, etc.) seem to have seven values per pixel:

detector.pedestals.shape: (7,16,352,384)
shot.shape: (16,352,384) # random data frame taken from exp=xcsx39618

Do these values correspond to different camera configurations? However, there appear to be only five configurations (two auto-ranging, three fixed), so I’m not sure which set of calibration constants is the appropriate one to use in different circumstances, e.g. for adding in pedestal-subtracted dark noise.

Due to active development (API changes etc), example scripts quickly become outdated and don't work anymore.

it would be great if we can run example scripts in pytest

@marcgri

I think this could be sped up with memoization. Converting these if-else calls to a dictionary lookup: e.g.
lookup{'CBILE':'CH', 'CBTHR':'CH', 'CBVAL':'CH', 'CBALA':'CH3', add other amino acid combinations, e.g. 'CBLEU':'CH2'}

@irischang020

I noted the recent following change commits 4da07da and b29311b, coming from the geometry file added in commit 95724f3:

     det = ps.PnccdDetector(
-        geom=ex_dir_+'/lcls/amo86615/'
-             'PNCCD::CalibV1/Camp.0:pnCCD.1/geometry/0-end.data',
+        geom=ex_dir_+'/input/geometry/0-end.data',
         beam=beam)

In Python 2 / psana 1, the PNCCD::CalibV1 part might be needed to get the pedestals.
The doc says, in an error message:

Sorry, at present, the package is not very smart. Please specify the path of the detector as deep as the geometry profile.
And example would be like: /reg/d/psdm/amo/experiment_name/calib/group/source/geometry/0-end.data where the '/calib/group/source/geometry/0-end.data' part is essential. The address before that part is not essential and can be replaced with your absolute address or relative address.

In Python 3 / psana 2, the amo86615 part is needed to get the pedestals (through the call to MongoDB).

The version of pysingfel on PyPI has not been updated in 3 years. Should we remove the PyPI project or try to update it?

update only works the first time.
A clean-up is required if we want to make it work properly, and to maybe integrate it with SASE later on.

Add tox.ini for python version

There are hard coded and relative paths that should be removed

Right now, we create an empty particle and then we populate it by reading a PDB file.
We should probably be able to provide the PDB file name to the constructor instead, as already done with HDF5 files.

Try oversampling the detector up to electron sampling, then bin.
That would allow us to stop assuming photons are centered in the middle of the pixel and start doing subpixel resolution stuff.

Hi @chuckie82, this radiationDamage.py script created a long time ago was implemented in SimEx to generate diffraction patterns from the PMI output of XMDYN. When I checked the codes, I found an inconsistency between what we expected and what we got from the code. The formulas were derived from the supplementary material of Chuck's 2016srep paper.

The code excerpt from pysingfel/radiationDamage.py:

        detector_intensity += photon_field*beam.get_photons_per_pulse_per_area()
   detector_intensity *= (det.solid_angle_per_pixel *
                           det.polarization_correction)

to manifest the expected function in the picture.

Do you have any comments on that? Is there a new script to conduct similar computing in this new version of pySingFEL?

By the way, The total_phot was calculated on line 160 but was not called in the code.

@irischang020
separate out smear method out of matricize method

Keep all atoms regardless of occupancy. Currently, throw away <50% occupancy atoms

The photon -> ADU transformation is unclear and requires a precomputed database.

I had a discussion with Haoyuan. He considers that photons produce electrons in a 5-by-5 pixel region around their hit point. The number of electrons is considered to be 130 and they are sampled according to a (truncated) Gaussian distribution in that 5-by-5 window.

Since his sampling method of taking 130 electron for each photon was extremely slow, he exported 1 million of pre-sampled distributions in a numpy dump file (the precomputed database). He then samples one of these 5-by-5 spreading kernels for each photon.
In practice, I think we could use

np.random.multinomial(130, the_5_by_5_distribution_weights,
                      (number_of_photons_in_the_image,))

to generate our kernels.

The sigma value is unclear. The sampling also assumes that the photon is centered on one of the pixels.

To interpret the data, we can use psana's ADU -> photon functionality presented in:
https://lcls-psana.github.io/Detector/index.html#module-AreaDetector

PnccdDetector needs the entire path to the geometry file because it uses some info in the file and some info in the path.
It might be worth uncoupling the two.

One feature that would be nice to have is:
2 displays side-by-side.
We show the particle on the left. We show the diffraction pattern on the right.
And allow the user to rotate the particle and display corresponding diffraction pattern near real-time.

@AntoineDujardin this seems to fail for py2 and py3

(ps-3.1.12) [yoon82@psanagpu101 pysingfel]$ pytest
============================= test session starts ==============================
platform linux -- Python 3.7.6, pytest-5.4.3, py-1.9.0, pluggy-0.13.1
PyQt5 5.12.3 -- Qt runtime 5.12.5 -- Qt compiled 5.12.5
rootdir: /reg/neh/home4/yoon82/Software/pysingfel
plugins: qt-3.3.0, asyncio-0.12.0
collected 96 items / 3 skipped / 93 selected

pysingfel/beam/tests/test_beam_base.py ................. [ 17%]
pysingfel/detector/tests/test_user_defined_detector.py ...... [ 23%]
pysingfel/geometry/tests/test_convert.py ............................... [ 56%]
[ 56%]
pysingfel/geometry/tests/test_generate.py ..FFF [ 61%]
pysingfel/geometry/tests/test_geometry.py ......... [ 70%]
pysingfel/geometry/tests/test_mapping.py ... [ 73%]
pysingfel/geometry/tests/test_slice_.py ........ [ 82%]
pysingfel/solvent_form_factor/tests/form_factor_test.py . [ 83%]
pysingfel/solvent_form_factor/tests/form_factor_variants_test.py .... [ 87%]
pysingfel/solvent_form_factor/tests/libsaxs_compare_test.py ..... [ 92%]
pysingfel/solvent_form_factor/tests/radial_distribution_test.py .... [ 96%]
pysingfel/tests/test_particle.py ... [100%]

=================================== FAILURES ===================================
_____________ TestPointsOn2Sphere.test_points_on_2sphere_moment_1 ______________

self = <test_generate.TestPointsOn2Sphere object at 0x7fe88435de90>

def test_points_on_2sphere_moment_1(self):
    """Test points_on_2sphere for moments of order 1."""
    thresh = 2.3e-3  # Lower limit for original implementation (1000)
    points = self.points.copy()
    self._test_points_on_2sphere_moment_1(points, thresh)
    with pytest.raises(AssertionError):

      self._test_points_on_2sphere_moment_1(points, thresh*.95)

E Failed: DID NOT RAISE <class 'AssertionError'>

pysingfel/geometry/tests/test_generate.py:48: Failed
---------------------------- Captured stdout setup -----------------------------
Deprecation warning: The function points_on_2sphere actually generates points on a 3-sphere, in 4D. Please call points_on_3sphere instead.
_____________ TestPointsOn2Sphere.test_points_on_2sphere_moment_2 ______________

self = <test_generate.TestPointsOn2Sphere object at 0x7fe884332c50>

def test_points_on_2sphere_moment_2(self):
    """Test points_on_2sphere for moments of order 2."""
    thresh_xx = 0.34  # Lower limit for original implementation (1000)
    thresh_xy = 1.16  # Lower limit for original implementation (1000)
    points = self.points.copy()

  self._test_points_on_2sphere_moment_2(points, thresh_xx, thresh_xy)

pysingfel/geometry/tests/test_generate.py:70:

self = <test_generate.TestPointsOn2Sphere object at 0x7fe884332c50>
points = array([[ 0.99115132, 0.07138709, 0.04162661, 0.10387578],
[ 0.99115132, 0.07138709, -0.10184068, -0.0463834... [-0.94078156, -0.12973254, 0.22756036, -0.21521108],
[-0.94078156, -0.3132062 , 0.06721999, 0.1109658 ]])
thresh_xx = 0.34, thresh_xy = 1.16

def _test_points_on_2sphere_moment_2(self, points, thresh_xx, thresh_xy):
    """Utilisty to test points_on_2sphere for moments of order 2."""
    # No need to normalize for double-cover because for even order
    moments_xx = []
    moments_xy = []
    for (i, j) in itertools.product(*(range(4),)*2):
        M2 = np.dot(points[:,i], points[:,j])
        if i == j:
            moments_xx.append(M2)
        else:
            moments_xy.append(M2)
    moments_xx = np.array(moments_xx)

  assert np.all(np.abs(moments_xx-moments_xx.mean()) < thresh_xx)

E AssertionError: assert False
E + where False = <function all at 0x7fe8997c0050>(array([ 4.44651051, 3.06173139, 3.65832287, 11.16656477]) < 0.34)
E + where <function all at 0x7fe8997c0050> = np.all
E + and array([ 4.44651051, 3.06173139, 3.65832287, 11.16656477]) = <ufunc 'absolute'>((array([245.55348949, 246.93826861, 246.34167713, 261.16656477]) - 250.0000000000001))
E + where <ufunc 'absolute'> = np.abs
E + and 250.0000000000001 = <built-in method mean of numpy.ndarray object at 0x7fe8815743f0>()
E + where <built-in method mean of numpy.ndarray object at 0x7fe8815743f0> = array([245.55348949, 246.93826861, 246.34167713, 261.16656477]).mean

pysingfel/geometry/tests/test_generate.py:62: AssertionError
_______________ TestPointsOn2Sphere.test_points_on_2sphere_angle _______________

self = <test_generate.TestPointsOn2Sphere object at 0x7fe8842ac6d0>

def test_points_on_2sphere_angle(self):
    """Test points_on_2sphere for angle between elements."""
    thresh = 0.105  # Lower limit for original implementation (1000)
    points = self.points.copy()
    self._test_points_on_2sphere_angle(points, thresh)
    with pytest.raises(AssertionError):

      self._test_points_on_2sphere_angle(points, thresh*.95)

E Failed: DID NOT RAISE <class 'AssertionError'>

pysingfel/geometry/tests/test_generate.py:90: Failed
=========================== short test summary info ============================
FAILED pysingfel/geometry/tests/test_generate.py::TestPointsOn2Sphere::test_points_on_2sphere_moment_1
FAILED pysingfel/geometry/tests/test_generate.py::TestPointsOn2Sphere::test_points_on_2sphere_moment_2
FAILED pysingfel/geometry/tests/test_generate.py::TestPointsOn2Sphere::test_points_on_2sphere_angle
=================== 3 failed, 93 passed, 3 skipped in 19.95s ===================

I will look into it soon, but put this here in case someone else runs into it. It seems that generate_new_sample_state() in the FXSExperiment class sometimes results in a Recursion Error. Example shown below for pysingfel installed while on ps-4.0.9 and using LCLS-II py3 kernel in Jupyter

The code:

input_dir='../input'
pdbfile=input_dir+'/pdb/2cex.pdb'
beamfile=input_dir+'/beam/amo86615.beam'
geom=input_dir+'/lcls/amo86615/PNCCD::CalibV1/Camp.0:pnCCD.1/geometry/0-end.data'

beam = ps.Beam(beamfile)
increase_factor = 1e2
beam.set_photons_per_pulse(increase_factor * beam.get_photons_per_pulse())

det = ps.PnccdDetector(geom=geom, beam=beam)
increase_factor = -0.5
det.distance = increase_factor * det.distance

particle = ps.Particle()
particle.read_pdb(pdbfile, ff='WK')

experiment = ps.FXSExperiment(det, beam, [particle], 100)

image = experiment.generate_image()

The error:

/cds/home/f/fpoitevi/Toolkit/pysingfel/pysingfel/particlePlacement.py:60: RuntimeWarning: invalid value encountered in true_divide
  direction = direction/np.linalg.norm(direction)
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-6-5bf0731ec4f4> in <module>
----> 1 image = experiment.generate_image()
      2 fig = plt.figure(figsize=(4,3), dpi=180)
      3 viz = ps.Visualizer(experiment, diffraction_rings="auto", log_scale=True)
      4 viz.imshow(image)

~/Toolkit/pysingfel/pysingfel/experiment/base.py in generate_image(self, return_orientation)
     36             return self.det.assemble_image_stack(img_stack), orientation
     37         else:
---> 38             img_stack = self.generate_image_stack()
     39             return self.det.assemble_image_stack(img_stack)
     40 

~/Toolkit/pysingfel/pysingfel/experiment/base.py in generate_image_stack(self, return_photons, return_intensities, return_orientation, always_tuple)
     62         beam_spectrum = self.beam.generate_new_state()
     63 
---> 64         sample_state = self.generate_new_sample_state()
     65 
     66         intensities_stack = 0.

~/Toolkit/pysingfel/pysingfel/experiment/fxs.py in generate_new_sample_state(self)
     37         particle_dict = {self.particles[i]: n_particles for i, n_particles in enumerate(particle_distribution)}
     38         print("particle_dict =", particle_dict)
---> 39         part_states, part_positions = distribute_particles(particle_dict, self.beam.get_focus()[0]/2, jet_radius=1e-4, gamma=0.5)
     40         part_states = np.array(part_states)
     41         for i in range(self.n_particle_kinds):

~/Toolkit/pysingfel/pysingfel/particlePlacement.py in distribute_particles(particles, beam_focus_radius, jet_radius, gamma)
     65     checkList = [overlap[i][j] for i in range(N) for j in range(N) if j > i]
     66     if any(item == True for item in checkList):
---> 67         distribute_particles(particles, beam_focus_radius, jet_radius, gamma)
     68     return state, coords
     69 

... last 1 frames repeated, from the frame below ...

~/Toolkit/pysingfel/pysingfel/particlePlacement.py in distribute_particles(particles, beam_focus_radius, jet_radius, gamma)
     65     checkList = [overlap[i][j] for i in range(N) for j in range(N) if j > i]
     66     if any(item == True for item in checkList):
---> 67         distribute_particles(particles, beam_focus_radius, jet_radius, gamma)
     68     return state, coords
     69 

RecursionError: maximum recursion depth exceeded while calling a Python object

Has the convention of what’s considered a positive versus negative detector distance changed? A PnCCD detector initialized from the current example .geom file (amo86615/PNCCD::CalibV1/Camp.0:pnCCD.1/geometry/0-end.data) has a negative distance:

detector = ps.PnccdDetector(geom = geom_file, beam = beam)
print(f"Detector distance is: {detector.distance:.2f} m")
Detector distance is: -0.58 m

Looking through example scripts, it looks like this detector distance was used previously. However, when computing the intensity field using detector.get_intensity_field, it seems that a positive detector distance yields the correct diffraction pattern (centrosymmetric, with intensity decreasing at high q) while the negative distance does not. @fredericpoitevin mentioned that orientation conventions in psana might have changed recently (thanks for the tip!). I find similar results for the PnCCDDetector class and a non-psana style SimpleSquareDetector, but I'm not sure if psana could affect conventions elsewhere in the code.

On the other hand, it seems like this is reversed for the calculated Compton scattering (calculate_compton in diffraction.py). In this case, a negative detector distance results in decreasing intensity as a function of q:

In addition, it looks like computing the intensity field using either calculate_molecular_form_factor_square in diffraction.py or detector.get_intensity_field (which calls the gpu code) yields different results, as shown below. In this case the patterns were calculated with a positive detector distance.

For completeness, here are the results using a negative detector distance:

Create a function that changes the PnCCD detector distance.
Update pixel_position accordingly and propagate the change (initialize_pixels_with_beam).

For the same beam, detector, and particle in the same orientation, I anticipated experiment.generate_image and detector.get_photons to yield the same result within Poisson noise/quantization error. However, they appear to be transposes, as shown in the attached. I think experiment.generate_images rotates the particle while detector.get_photons rotates the diffraction volume, but I don’t understand why this would swap the axes of the resulting image.

There is a 41MB file at examples/scripts/imStack.hdf5, which should be removed & purged.

In the form factor tables, the variant O2p is key both to O and OH.
Since it's currently a suite of conditional statements, only the first one is used, but that might hide a bug.

Do we have a source for these coefficients?

@AntoineDujardin
It would be nice to have a UserDefinedDetector example for people using pysingfel without psana.

I couldn't find an example notebook:

@AntoineDujardin
Install pysingfel with
pip install -e .
and define pysingfel.dir to eliminate need to define ROOT_DIR and sys.path in the code

This instruction would need to go on top of the README

It seems that unit tests in test_diffraction.py and for the detector class are failing due to the missing geometry file: ‘PNCCD::CalibV1/Camp.0:pnCCD.1/geometry/0-end.data'. Tests that rely on a user-defined detector pass.