feat: updatad analysis script, thanks to Bhawani

Author: c-dilks

qa-physics/charge_analysis/analyze.py

@@ -115,11 +115,12 @@ def main():
     chunk_size = 2000
     num_chunks = len(timestamps) // chunk_size
 
-    chunk_caseA, chunk_caseB, chunk_caseC, chunk_default = [], [], [], []
-    cum_caseA, cum_caseB, cum_caseC, cum_default = [], [], [], []
+
+    chunk_caseA, chunk_caseB, chunk_caseC, chunk_default, chunk_default_ungated = [], [], [], [], []
+    cum_caseA, cum_caseB, cum_caseC, cum_default, cum_default_ungated = [], [], [], [], []
     chunk_indices, skipped_counts = [], []
 
-    runA, runB, runC, runDef = 0, 0, 0, 0
+    runA, runB, runC, runDef, runDefUng = 0, 0, 0, 0, 0
     total_skipped = 0
 
     corrected_livetimes_A = []
@@ -137,7 +138,7 @@ def main():
         fcupgated_diff = np.diff(fcupgateds[start:end])
         live_sub = live_times[start+1:end]
 
-        sumA, sumB, sumC, sumDef = 0, 0, 0, 0
+        sumA, sumB, sumC, sumDef, sumDefUng = 0, 0, 0, 0, 0
         skipped_in_chunk = 0
 
         for j, lt in enumerate(live_sub):
@@ -153,6 +154,8 @@ def main():
                 corrected_livetimes_C.append(lt)
                 # Default
                 sumDef += fcupgated_diff[j]
+                # Default ungated
+                sumDefUng += fcup_diff[j]
             else:
                 # ----- Case A/B nearest-neighbor substitution -----
                 idx_candidates = []
@@ -178,6 +181,9 @@ def main():
 
                     # Default
                     sumDef += fcupgated_diff[j]
+
+                    # Default ungated
+                    sumDefUng += fcup_diff[j]
                 else:
                     skipped_in_chunk += 1
                     total_skipped += 1
@@ -199,25 +205,28 @@ def main():
         runB += sumB
         runC += sumC
         runDef += sumDef
+        runDefUng += sumDefUng
 
         chunk_caseA.append(sumA)
         chunk_caseB.append(sumB)
         chunk_caseC.append(sumC)
         chunk_default.append(sumDef)
+        chunk_default_ungated.append(sumDefUng)
         cum_caseA.append(runA)
         cum_caseB.append(runB)
         cum_caseC.append(runC)
         cum_default.append(runDef)
+        cum_default_ungated.append(runDefUng)
         chunk_indices.append(i)
         skipped_counts.append(skipped_in_chunk)
 
     logger.info(f"Computed chunked FCUP Gated values with neighbor handling (Cases A, B, C).")
     logger.info(f"Total skipped events (no valid LT neighbor): {total_skipped}")
 
     # ---------- Plot 2: Chunked FCUP Gated + Ratios + Skips + LT Distribution ----------
-    fig2, (ax_top, ax_mid, ax_bottom, ax_ltdist) = plt.subplots(
-        4, 1, figsize=(12, 14), sharex=False,
-        gridspec_kw={'height_ratios': [3, 1, 1, 2]}
+    fig2, (ax_top, ax_mid, ax_gatedrat, ax_bottom, ax_ltdist) = plt.subplots(
+        5, 1, figsize=(12, 17), sharex=False,
+        gridspec_kw={'height_ratios': [3, 1, 1, 1, 2]}
     )
     fig2.suptitle(f'Run {run_number} - Chunked FCUP Gated (Neighbor Handling)', fontsize=16)
 
@@ -226,6 +235,7 @@ def main():
     #ax_top.plot(chunk_indices, cum_caseB, label='Cumulative Case B (LT_nn × FCUPungated_nn)', color='darkgreen', marker='s')
     ax_top.plot(chunk_indices, cum_caseC, label='Cumulative Case C (20-NN mean × FCUPungated)', color='darkorange', marker='d')
     ax_top.plot(chunk_indices, cum_default, label='Cumulative Default (FCUPgated)', color='blue', marker='^')
+    ax_top.plot(chunk_indices, cum_default_ungated, label='Cumulative Default Ungated (FCUPungated)', color='teal', marker='x',linestyle='--')
     ax_top.set_ylabel('Cumulative Σ', fontsize=11)
     ax_top.grid(True, linestyle='--', alpha=0.6)
     ax_top.legend(fontsize=10)
@@ -235,16 +245,27 @@ def main():
     ratioA = np.divide(cum_caseA, cum_default, out=np.full_like(cum_caseA, np.nan, dtype=float), where=np.array(cum_default) != 0)
     #ratioB = np.divide(cum_caseB, cum_default, out=np.full_like(cum_caseB, np.nan, dtype=float), where=np.array(cum_default) != 0)
     ratioC = np.divide(cum_caseC, cum_default, out=np.full_like(cum_caseC, np.nan, dtype=float), where=np.array(cum_default) != 0)
+    ratioDefUng = np.divide(cum_default_ungated, cum_default, out=np.full_like(cum_default_ungated, np.nan, dtype=float), where=np.array(cum_default) != 0)
 
     ax_mid.plot(chunk_indices, ratioA, label='Case A / Default', color='darkred', marker='o')
     #ax_mid.plot(chunk_indices, ratioB, label='Case B / Default', color='darkgreen', marker='s')
     ax_mid.plot(chunk_indices, ratioC, label='Case C / Default', color='darkorange', marker='d')
+    ax_mid.plot(chunk_indices, ratioDefUng, label='Default FCUP Ungated / Default FCUP gated', color='teal', marker='x',linestyle='--')
     ax_mid.axhline(1.0, color='black', linestyle='--', linewidth=1)
     ax_mid.set_ylabel('Ratio', fontsize=11)
     ax_mid.grid(True, linestyle='--', alpha=0.6)
     ax_mid.legend(fontsize=10)
     ax_mid.tick_params(axis='both', labelsize=10)
 
+    # Gated / Ungated ratio panel
+    ratio_gated_ung = np.divide(cum_default, cum_default_ungated, out=np.full_like(cum_default, np.nan, dtype=float), where=np.array(cum_default_ungated) != 0)
+    ax_gatedrat.plot(chunk_indices, ratio_gated_ung, label='FCUPgated / FCUPungated', color='navy', marker='o')
+    ax_gatedrat.axhline(1.0, color='black', linestyle='--', linewidth=1)
+    ax_gatedrat.set_ylabel('Gated / Ungated', fontsize=11)
+    ax_gatedrat.grid(True, linestyle='--', alpha=0.6)
+    ax_gatedrat.legend(fontsize=10)
+    ax_gatedrat.tick_params(axis='both', labelsize=10)
+
     # Bottom-1: skipped events count
     ax_bottom.bar(chunk_indices, skipped_counts, color='gray', alpha=0.7)
     ax_bottom.set_xlabel(f'Chunk Index (Each = {chunk_size} events)', fontsize=11)