The Detectors
=========================

The :code:`scatman.Detectors` submodule contains a set of objects (at the moment 3) that aim to simulate the behavior of different kinds of CDI detectors. 

All detectors derive from an *abstract* base class :code:`scatman.Detectors.Detector`. Being an *abstract* class, it cannot be used directly, i.e. it is not possible to instantiate an object of that type:

.. code::
    
    >>> import scatman
    >>> MyDetector = scatman.Detectors.Detector()
    Traceback (most recent call last):
        File "<stdin>", line 1, in <module>
    TypeError: scatman.Detectors.Detector: No constructor defined!
    
    
For the moment, the following detectors are implemented:

:MSFT: It is actually not a real detector. It just provides the diffraction pattern as it is computed via the MSFT.
:Ideal: The Ideal detector is the one that perfectly detects each photon arriving on it. It basically adds just the photon statistics (Poisson) to the MSFT result.
:MCP: The Micro-Channel Plate simulation involves many features. For the moment, just its non-linear response function is included in the simulation. Things like the 8 degrees hole will be added in future.



.. _coordinates_description:

Detector Coordinates
---------------------

The output of the MSFT gives the result for an ideal *curved* detector, i.e. a detector where each *idealized* pixel has the same cross section and the same distance from the interaction zone between the sample and the beam. 

It is often necessary to translate the pixel position in the diffraction pattern, basically its **x** and **y** coordinates, into a physical quantity of interest. Usually, common physical quantities are the *transfer momentum* and the *scattering angle*. So, it is often desiderable that the single pixel size is a constant amount of *transfer momentum* or *scattering angle*. 

However, the raw MSFT result provides a scattering image whose pixels are equally spaced on the plane orthogonal to the beam propagation direction. In practice, this means that each pixel has a constant size :math:`S_p^\text{proj}` that is a fraction of the maximum *transfer momentum* projected on the *xy* plane, i.e.:

.. math::

    S_p^\text{proj} = \frac{|k_\text{in}| \sin(\theta_\text{max})}{\frac{N_p}{2}}

where :math:`|k_\text{in}|` is the radiation wavenumber, :math:`\theta_\text{max}` is the maximum scattering angle and :math:`N_p` is the linear amount of pixels of the detector.

It is however possible to perform a *change of coordinates* and obtain different representations. It's worth noting that this operation is done by **interpolating** the raw MSFT data, and the result must be thus considered just an approximation of the correct result.

An useful coordinates system is the one where each pixel has a size :math:`S_p^\text{mom}` that corresponds to a constant unit of the *transfer momentum*, that is:

.. math::

    S_p^\text{mom} = \frac{2 |k_\text{in}| \sin \left( \frac{\theta_\text{max}}{2} \right)}{\frac{N_p}{2}}
    
Another one is, instead, the case in which the pixel size :math:`S_p^\text{ang}` has a constant unit of scattering angle, which means:

.. math::

    S_p^\text{ang} = \frac{\theta_\text{max}}{\frac{N_p}{2}}

The last case is when the pixel size :math:`S_p^\text{det}` of the pattern directly corresponds to a pixel of the real detector, that is:

.. math::

    S_p^\text{ang} = \frac{\tan(\theta_\text{max})}{\frac{N_p}{2}}
    
The Scatman is capable of providing these different representations of the diffraction pattern, by setting a proper value for the parameter :code:`coordinates` when a detector is defined. Its value can be:

- :code:`'projection'`: keep the *natural* MSFT coordinates system. No interpolation is performed
- :code:`'momentum'`: Interpolate the MSFT result, such that the pixel size is a constant fraction of the *transfer momentum*.
- :code:`'angle'`: Interpolate the MSFT result, such that the pixel size is a constant fraction of the *scattering angle*.
- :code:`'detector'`: Interpolate the MSFT result, such that the pixel size of the output corresponds to the pixel size of a real detector. **In this case, the intensity of the scattering is scaled too**, to match what is acquired by a *real detector*, where the cross section of the pixels varies as function of the distance from the detector center.
    
The following scheme summarizes the three options for the choice of the coordinates:

.. figure:: figures/coordinates.png
    :width: 80%
    :align: center



.. _MCP_response_description:


Modelling the photons statistics
--------------------------------

The amount of scattered photons depends on the radiation absorbed by the sample. In particular, here we assume that the amount of absorbed radiation is equal to the scattered one. 
The distribution if the power density absorbed by the sample :math:`H_a(x,y)` at coordinates  :math:`x,y`  is computed through the Lamber-Beer law (with the radiation propagating along the :math:`z` direction):

.. math:: 

    H_a(x,y) = H_i \left [ 1 - e^{2 k_0 \int dz \beta(x,y,z)} \right ]
    
where :math:`H_i` is the incident power density, assumed to be uniform on the :math:`x,y` plane, and :math:`\beta(x,y,z)` is the absorption coefficient of the sample material.
Then, the total absorbed power :math:`W_a` is:

.. math:: 

    W_a = \iint dx dy H_a(x,y) = H_i  \iint dx dy \left [ 1 - e^{2 k_0 \int dz \beta(x,y,z)} \right ]
    
As stated before, the Scatman assumes a completely elasitc scattering, which implies that the total scattered power :math:`W_s` is equal to the total absorbed one :math:`W_a`.
Moreover, being the power directly proportional to the amount of photons, the equation can be rewritten as:

.. math:: 

    N_s = N_a = \rho_i  \iint dx dy \left [ 1 - e^{2 k_0 \int dz \beta(x,y,z)} \right ]

where  :math:`N_s` and :math:`N_a` is the amount of scattered and absorbed photons, respectively, and :math:`\rho_i` is the incoming photon density (i.e. photons per area unit).

The final assumption made by the current implementation is that the :math:`N_s` photons are all scattered towards the detector.



The following example deals with the simulation of a realistic diffraction pattern, i.e. including photons statistics, through the use of the :code:`scatman.Detectors.Ideal` detector.
To simulate photons statistics, the photon density :math:`\rho` has to be computed, and depends on the experimental parameters. Here, we simulate an experiment with radiation wavelength :math:`\lambda=40\text{nm}=31\text{eV}` and incoming power density :math:`H = 3 \cdot 10^{14} \frac{\text{W}}{\text{cm}^2} =  3 \cdot 10^{18} \frac{\text{W}}{\text{m}^2}`. The duration of the laser shot is :math:`t = 80\text{fs} = 8\cdot 10^{-14}\text{s}`.
The photon density can be calculated with the following formula:

.. math::

    \rho = \frac{H \cdot t}{q \cdot E} = 4.84\cdot 10^{22} \frac{\text{photons}}{\text{m}^2} = 4.84\cdot 10^{4} \frac{\text{photons}}{\text{nm}^2} 
    
where :math:`q=1.6\cdot 10^{-19} \text{C}` is the elementary charge.

The following code is a simulation of a spherical sample at these experimental conditions, with a detector covering up to 30 degrees scattering angle. 

.. code::

    import scatman

    scatman.set_experiment(wavelength=40., angle=30, resolution=512, photon_density=4.84e4);

    mydetector = scatman.Detectors.Ideal()

    myshape = scatman.Shapes.SphericalHarmonics(a=400, coefficients=[[0,0,1]], delta = 0.01, beta=0.02)

    mypattern = mydetector.acquire(myshape)


.. figure:: figures/example_ideal.png
    :width: 60%
    :align: center


Simulation of the MCP response
-------------------------------

Modelling the MCP response function (i.e. finding out which is the MCP output at a given input) is all but trivial. 
given :math:`I_\text{ideal}` the input intensity value of the scattered field, the MCP output response :math:`I_\text{MCP}` can be modeled with a function :math:`R` such that :math:`I_\text{MCP}=R(I_\text{ideal})`.

In the current implementation, the response function is defined as it follows:

.. math::
    
    R(x)=
    \begin{cases}
      G \cdot x^E, & \text{if}\ G \cdot x^E < L \cdot S \\
      S \cdot \left [ (1-L) \cdot (1 - e^{\frac{1}{L-1} \cdot (\frac{G \cdot x^E}{S} - L)}) + L \right ] , & \text{otherwise}
    \end{cases}
    
where:
    
    - :math:`G` refers to what I called :code:`gain`
    - :math:`E` refers to what I called :code:`exponent`, usually :math:`< 1` (the MCP response is usually sub-linear)
    - :math:`S` is the :code:`saturation` value of the MCP output
    - :math:`L` is the (sub) :code:`linearity` level, i.e. the fraction of the saturation value up to which the response follows the power-law model

All these parameters can be set when defining the :meth:`scatman.Detectors.MCP` .    

Despite this function looks quite messy, it is just a power law that at some point becomes an exponential that goes to the saturation level. The complicated expression is due to the fact that, at the point where there is the transition between the sub-linear behavior and the exponential one, the continuity of both the function and its first derivative are imposed.

To get a better idea, the following plot depicts the meaning of the parameters in a graphical way:

.. figure:: figures/MCP_response.png
    :width: 80%
    :align: center
    
    

.. .. autoclass:: scatman.Detectors.Detector
..     :members:


List of Detectors
------------------

.. autoclass:: scatman.Detectors.MSFT
    :exclude-members: acquire_dataset
    :show-inheritance:
    :members:
    :inherited-members:

.. autoclass:: scatman.Detectors.Ideal
    :exclude-members: acquire_dataset
    :show-inheritance:
    :members:
    :inherited-members:
    
.. autoclass:: scatman.Detectors.MCP
    :exclude-members: acquire_dataset
    :show-inheritance:
    :members:
    :inherited-members:

    
..     .. automethod:: scatman.Detectors.Detector.acquire
    

..     .. automethod:: acquire(self, shape)

.. .. automodule:: scatman.Detectors
..     :undoc-members:
..     :noindex:
.. .. autoclass:: scatman.Detectors.Detector
..     :members:
    
.. raw:: latex

    \newpage