From 7e6ff243037b074a65e144ff02877d5ac2533a0e Mon Sep 17 00:00:00 2001
From: Laurent Mercadier <laurent.mercadier@xfel.eu>
Date: Mon, 14 Mar 2022 21:52:49 +0100
Subject: [PATCH] add fluence calibration routine

---
 doc/changelog.rst                   |  2 +-
 src/toolbox_scs/misc/laser_utils.py | 90 +++++++++++++++++++++++++++--
 2 files changed, 87 insertions(+), 5 deletions(-)

diff --git a/doc/changelog.rst b/doc/changelog.rst
index ee15b95..ed49140 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -19,7 +19,7 @@ unreleased
     - improve SLURM scripts with named arguments :mr:`176`
   
 - **New Features**
-
+    - add routine for fluence calibration :mr:`180`
 
 1.5.0
 -----
diff --git a/src/toolbox_scs/misc/laser_utils.py b/src/toolbox_scs/misc/laser_utils.py
index b284999..99698e1 100644
--- a/src/toolbox_scs/misc/laser_utils.py
+++ b/src/toolbox_scs/misc/laser_utils.py
@@ -1,10 +1,14 @@
 __all__ = [
     'degToRelPower',
     'positionToDelay',
+    'fluenceCalibration',
 ]
 
+import numpy as np
+import matplotlib.pyplot as plt
 
-def positionToDelay(pos, origin=0, invert = False, reflections=1):
+
+def positionToDelay(pos, origin=0, invert=True, reflections=1):
     ''' converts a motor position in mm into optical delay in picosecond
         Inputs:
             pos: array-like delay stage motor position
@@ -15,18 +19,96 @@ def positionToDelay(pos, origin=0, invert = False, reflections=1):
         Output:
             delay in picosecond
     '''
-    c_ = 299792458 *1e-9 # speed of light in mm/ps
+    c_ = 299792458 * 1e-9  # speed of light in mm/ps
     x = -1 if invert else 1
     return 2*reflections*(pos-origin)*x/c_
-    
+
+
 def degToRelPower(x, theta0=0):
     ''' converts a half-wave plate position in degrees into relative power
         between 0 and 1.
         Inputs:
             x: array-like positions of half-wave plate, in degrees
             theta0: position for which relative power is zero
-            
+
         Output:
             array-like relative power
     '''
     return np.sin(2*(x-theta0)*np.pi/180)**2
+
+
+def fluenceCalibration(hwp, power_mW, npulses, w0x, w0y=None,
+                       train_rep_rate=10, fit_order=1,
+                       plot=True, xlabel='HWP [%]'):
+    """
+    Given a measurement of relative powers or half wave plate angles
+    and averaged powers in mW, this routine calculates the corresponding
+    fluence and fits a polynomial to the data.
+
+    Parameters
+    ----------
+    hwp: array-like (N)
+        angle or relative power from the half wave plate
+    power_mW: array-like (N)
+        measured power in mW by powermeter
+    npulses: int
+        number of pulses per train during power measurement
+    w0x: float
+        radius at 1/e^2 in x-axis in meter
+    w0y: float, optional
+        radius at 1/e^2 in y-axis in meter. If None, w0y=w0x is assumed.
+    train_rep_rate: float
+        repetition rate of the FEL, by default equals to 10 Hz.
+    fit_order: int
+        order of the polynomial fit
+    plot: bool
+        Plot the results if True
+    xlabel: str
+        xlabel for the plot
+
+    Output
+    ------
+    F: ndarray (N)
+        fluence in mJ/cm^2
+    fit_F: ndarray
+        coefficients of the fluence polynomial fit
+    E: ndarray (N)
+        pulse energy in microJ
+    fit_E: ndarray
+        coefficients of the fluence polynomial fit
+    """
+    power = np.array(power_mW)
+    hwp = np.array(hwp)
+    E = power/(train_rep_rate*npulses)*1e-3  # pulse energy in J
+    if w0y is None:
+        w0y = w0x
+    F = 2*E/(np.pi*w0x*w0y)  # fluence in J/m^2
+
+    fit_E = np.polyfit(hwp, E*1e6, fit_order)
+    fit_F = np.polyfit(hwp, F*1e-1, fit_order)
+    x = np.linspace(hwp.min(), hwp.max(), 100)
+    if plot:
+        fig, ax = plt.subplots(figsize=(6, 4))
+        ax.set_title(f'w0x = {w0x*1e6:.0f} $\mu$m, w0y = {w0y*1e6:.0f} $\mu$m')
+        ax.plot(hwp, F*1e-1, 'o', label='data')
+        fit_label = 'F = '
+        for i in range(len(fit_F)-1, 1, -1):
+            fit_label += f'{fit_F[i]:.2g}x$^{i}$ + '
+            if i % 2 == 0:
+                fit_label += '\n'
+        fit_label += f'{fit_F[-2]:.2g}x + {fit_F[-1]:.2g}'
+        ax.plot(x, np.poly1d(fit_F)(x), label=fit_label)
+        ax.set_ylabel('Fluence [mJ/cm$^2$]')
+        ax.set_xlabel(xlabel)
+        ax.legend()
+        ax.grid()
+
+        def eTf(x):
+            return 1e-7*2*x/(np.pi*w0x*w0y)
+
+        def fTe(x):
+            return 1e7*x*np.pi*w0x*w0y/2
+        ax2 = ax.secondary_yaxis('right', functions=(fTe, eTf))
+        ax2.set_ylabel(r'Pulse energy [$\mu$J]')
+        fig.tight_layout()
+    return F*1e-1, fit_F, E*1e6, fit_E
-- 
GitLab