use a background image for fluorescence imaging

fluorescence_imaging.py

@@ -12,7 +12,7 @@ Fluorescence imaging
 # CAMERA OPTIONS
 camera_config = {
     "exposure_time": 12549.1,
-    "gain": 18.0,
+    "gain": 10.0,
     "gamma": 1.25,
     "black_level": 5.0,
     "bin_size": 1,
@@ -23,19 +23,19 @@ t0 = 3.0 # take the picture
 
 SEQ = SuperSequence(outdir, "sequence", {
     "Sequence": (Sequence
-        .digital_hilo(*C.dummy, 0.0, t0 + 4000e-3)
+        .digital_hilo(*C.dummy, 0.0, t0 + 5500e-3)
         .with_color("k").with_stack_idx(0)
     ),
     "Camera 1": (Sequence
         .digital_pulse(*C.flir_trig, t0, 12.5491e-3)
         .with_color("C2").with_stack_idx(4)
     ),
-    #"Camera 2": (Sequence
-    #    .digital_pulse(*C.flir_trig, t0 + 3800e-3, 12.5491e-3)
-    #    .with_color("C2").with_stack_idx(4)
-    #),
+    "Camera 2": (Sequence
+        .digital_pulse(*C.flir_trig, t0 + 5000e-3, 12.5491e-3)
+        .with_color("C2").with_stack_idx(4)
+    ),
     "Push beam": (Sequence
-        .digital_lohi(*C.push_sh, t0 - 10e-3, t0 + 4000e-3)
+        .digital_lohi(*C.push_sh, t0 - 10e-3, t0 + 5500e-3)
         .with_color("C3").with_stack_idx(3)
     ),
     "Scope": (Sequence
@@ -70,10 +70,10 @@ class FluorescenceImaging(Controller):
         SEQ.save()
 
     def run_camera(self, *args):
-        frames = self.cam.acquire_frames(1)
+        frames = self.cam.acquire_frames(2)
         data = FluorescenceData(
             outdir=outdir,
-            arrays={"image": frames[0]},
+            arrays={"image": frames[0], "background": frames[1]},
             config=camera_config,
             comments=comments
         )

lib/imaging.py

@@ -165,15 +165,19 @@ class FluorescenceData(ImagingData):
         _dA = self.config.get("bin_size", 1)**2 \
                 * (DEF_DX**2 if dA is None else dA)
         if printflag: print("[imaging] Compute results")
-        #N_bkgd = compute_mot_number(
-        #    self.arrays["background"].astype(np.float64),
-        #    0.40,
-        #    self.config["gain"], # dB
-        #    self.config["exposure_time"] * 1e-6, # s
-        #    (0.0127**2 * np.pi) / (4 * np.pi * 0.125**2),
-        #    50e6 * 2 * np.pi,
-        #    1.0
-        #) if "background" in self.arrays.keys() else 0.0
+        compute_mot_number_params = {
+            "QE": 0.40,
+            "gain": self.config["gain"],
+            "exposure_time": self.config["exposure_time"] * 1e-6,
+            "solid_angle": (0.0127**2 * np.pi) / (4 * np.pi * 0.125**2),
+            "detuning": 50e6 * 2 * np.pi,
+            "intensity_parameter": 1.0,
+        }
+        N_bkgd = compute_mot_number(
+            self.arrays["background"].astype(np.float64),
+            **compute_mot_number_params,
+            k=0
+        ) if "background" in self.arrays.keys() else 0.0
         for label, array in self.arrays.items():
             if label == "background":
                 continue
@@ -182,18 +186,14 @@ class FluorescenceData(ImagingData):
                 print(f"[imaging] WARNING: image '{label}' may contain clipping")
             N = compute_mot_number(
                 array.astype(np.float64),
-                0.40,
-                self.config["gain"], # dB
-                self.config["exposure_time"] * 1e-6, # s
-                (0.0127**2 * np.pi) / (4 * np.pi * 0.125**2),
-                50e6 * 2 * np.pi,
-                1.0
+                **compute_mot_number_params,
+                k=1
             )
             sx, sy = lls_fit_gaussian(array, _dA)
             self.results = dict() if self.results is None else self.results
             self.results.update({
                 label: {
-                    "N": float(N),
+                    "N": float(N) - float(N_bkgd),
                     "sx": float(sx),
                     "sy": float(sy)
                 }
@@ -201,38 +201,30 @@ class FluorescenceData(ImagingData):
         return self
 
 def compute_mot_number(image, QE, gain, exposure_time, solid_angle,
-        detuning, intensity_parameter):
+        detuning, intensity_parameter, k=0):
     H, W = image.shape
-    K = 4
-    im = image[H//K:(K - 1)*H//K, W//K:(K - 1)*W//K]
-    im_bkgd = image.copy()
-    im_bkgd[H//K:(K - 1)*H//K, W//K:(K - 1)*W//K] = 0
-    #pd.Plotter().imshow(image).savefig("image.png").close()
-    #pd.Plotter().imshow(im).savefig("img.png").close()
-    #pd.Plotter().imshow(im_bkgd).savefig("img_bkgd.png").close()
+    K_h = 3.5
+    slice_h = pd.S[int(H / K_h) : int((K_h - 1) * H / K_h)]
+    K_w = 3.0
+    slice_w = pd.S[int(W / K_w) : int((K_w - 1) * W / K_w)]
+    im = image[slice_h, slice_w]
+    #pd.Plotter().imshow(image).savefig(f"image_{k}.png").close()
+    #pd.Plotter().imshow(im).savefig(f"img_{k}.png").close()
     electron_rate = im.astype(np.float64).sum()
-    electron_rate_bkgd = im_bkgd.astype(np.float64).sum()
-    photon_rate \
-            = electron_rate \
-            / 10**(gain / 10) \
-            / QE \
-            / solid_angle \
-            / exposure_time
-    photon_rate_bkgd \
-            = electron_rate_bkgd \
-            / 10**(gain / 10) \
-            / QE \
-            / solid_angle \
-            / exposure_time
+    photon_rate = (
+        electron_rate
+        / 10**(gain / 10)
+        / QE
+        / solid_angle
+        / exposure_time
+    )
     Y = 29.1e6 * 2 * np.pi # transition linewidth; s^-1
     l = 399e-9 # transition wavelength; m
     I_sat = (Y * H * C * np.pi) / (3 * l**3)
     scatter_rate = lambda D, I: \
             (Y / 2) * I / (1 + 4 * (D / Y)**2 + I)
     N = photon_rate / scatter_rate(detuning, intensity_parameter)
-    N_bkgd = (K - 2)**2 / (4 * K - 4) \
-            * photon_rate_bkgd / scatter_rate(detuning, intensity_parameter)
-    return N - N_bkgd
+    return N
 
 def lls_fit_gaussian(A, dA):
     """