Dark Current Varation anpassen

reflectance_main.py

@@ -24,9 +24,12 @@ import numpy as np
 # processor = ReflectanceProcessor("24", "Aug", "140824", "kla", 31, 31, 83, 81, 85)
 # processor = ReflectanceProcessor("24", "Aug", "200824", "kla", 48, 48, 49, 47, 50)
 # processor = ReflectanceProcessor("24", "Aug", "210824", "kla", 16, 18, 17, 15, 19)
-processor = ReflectanceProcessor("24", "Aug", "140824", "kla", 31, 31, 30, 28, 32)
+# processor = ReflectanceProcessor("24", "Aug", "140824", "kla", 31, 31, 30, 28, 32)
+# processor = ReflectanceProcessor("25", "Jun", "230625", "kla", 29, 31, 30, 28, 32)
+processor = ReflectanceProcessor("25", "Jun", "230625", "kla", 29, 31, 30, 28, 32, [25, 27, 29, 31, 33, 35])
+
 wavelen, theta, refl_diode, ref_before, ref_after, uncert_rel, uncert_abs = processor.compute_reflectance()
 
-print('refl_diode', refl_diode)
-print('wavelen, uncert_rel', wavelen, uncert_rel, np.amax(uncert_rel))
+# print('refl_diode', refl_diode)
+# print('wavelen, uncert_rel', wavelen, uncert_rel, np.amax(uncert_rel))
 

reflectance_processor.py

@@ -4,7 +4,7 @@ import os
 from EveHDF import EveHDF
 import scipy
 import sys
-from typing import List, Tuple  # Import necessary types
+from typing import Optional, List, Tuple  # Import necessary types
 
 if sys.version_info >= (3, 9):
     ColumnListType = list[str]  # Native list[str] for Python 3.9+
@@ -95,35 +95,6 @@ def get_var_data():
 
     return DC_diode_var_s, DC_diode_var_p, DC_detector_var_s, DC_detector_var_p
 
-# def calculate_uncertainty(I0_drift: np.ndarray, RT: np.ndarray, Reflectance: np.ndarray,
-#                           cov_factor: int = 2, *extra_terms: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
-#     """
-#     Berechnet die relative und absolute Unsicherheit der Reflektanzmessung.
-#
-#     :param I0_drift: Drift der Intensität I0, als Array
-#     :param RT: Unsicherheitskomponente durch Messungenauigkeit, als Array
-#     :param Reflectance: Berechnete Reflektanzwerte, als Array
-#     :param cov_factor: Ganzzahliger Coverage-Faktor für Unsicherheitsberechnung (1, 2 oder 3) (Standard: 2 für 95% Konfidenz)
-#     :param extra_terms: Optionale weitere Unsicherheitsbeiträge, die quadratisch addiert werden
-#     :return: Tuple (relative Unsicherheit, absolute Unsicherheit)
-#     """
-#
-#     if cov_factor not in {1, 2, 3}:
-#         raise ValueError("cov_factor muss 1, 2 oder 3 sein.")
-#
-#     # hier stattdessen das Package uncertainties nutzen?
-#     # Quadratische Addition aller Unsicherheitsbeiträge
-#     total_uncert = (100 * I0_drift) ** 2 + (100 * RT) ** 2
-#
-#     # Falls zusätzliche Unsicherheitsquellen übergeben wurden, addiere sie quadratisch
-#     for term in extra_terms:
-#         total_uncert += (100 * term) ** 2
-#
-#     uncert_rel = np.sqrt(total_uncert)
-#     uncert_abs = (cov_factor * uncert_rel / 100) * abs(Reflectance)
-#
-#     return uncert_rel, uncert_abs
-
 def get_rel_error_diode_homogenity(wavelength):
     """
     Gibt den relativen Fehler durch Diodenhomogenität für eine gegebene Wellenlänge.
@@ -194,7 +165,8 @@ def calculate_uncertainty(wavelength, I0_drift: np.ndarray, RT: np.ndarray, Refl
 class ReflectanceProcessor:
     def __init__(self, year: str, month: str, date: str, acronym: str,
                  dark_current_1: int, dark_current_2: int, reflect_raw: int,
-                 reference_before: int, reference_after: int):
+                 reference_before: int, reference_after: int,
+                 dark_current_var_indices: Optional[List[int]] = None):
         """
         Initialisiert einen ReflectanceProcessor mit den erforderlichen Parametern.
 
@@ -218,6 +190,7 @@ class ReflectanceProcessor:
         self.reflect_raw_data = self.generate_filename(reflect_raw)
         self.reference_before = self.generate_filename(reference_before)
         self.reference_after = self.generate_filename(reference_after)
+        self.dark_current_var_indices = dark_current_var_indices or []
 
         # Spaltennamen für verschiedene Datentypen
         self.usecols_reflec_names = ['Wavlen', 'Theta', 'Keysi_03', 'Keysi_04', 'mlsRingMaschine', 'mlsRing_1']
@@ -354,6 +327,32 @@ class ReflectanceProcessor:
 
         return wavelen, theta
 
+    def compute_dark_current_variation(self, indizes: List[int]) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Berechnet die Variation (Betrag Mittelwert + Standardabweichung) der Dunkelstrommessungen
+        für Detektor und Diode.
+
+        :param indizes: Liste von Indizes für die Dunkelstromdateien
+        :return: Tuple (variation_diode, variation_detector)
+        """
+        if len(indizes) < 6:
+            print(
+                f'Dark current length of {len(indizes)} is too small to significantly calculate a standard deviation.')
+
+        dfs = []
+        for index in indizes:
+            file_path = self.generate_filename(index)
+            df = self._load_h5(file_path, self.usecols_dark_names).dropna()
+            dfs.append(df)
+
+        def stats(key: str) -> np.ndarray:
+            combined = pd.concat([df[key].reset_index(drop=True) for df in dfs], axis=1)
+            return (combined.abs().mean(axis=1) + combined.std(axis=1)).to_numpy()
+
+        var_diode = stats("Measur_Diod")
+        var_detector = stats("Measur_Detec")
+        return var_diode, var_detector
+
     def compute_reflectance(self) -> TupleType[np.ndarray, np.ndarray, np.ndarray]:
         """
         Berechnet die Reflexion basierend auf den Messdaten und Referenzdaten.
@@ -362,20 +361,31 @@ class ReflectanceProcessor:
         """
         ######
         # hier klarstellen, wo diese Daten herkommen
-        # ist es klug die hie rhard reinzucoden?
+        # ist es klug die hier hard reinzucoden?
         # lieber woanders hin?
         # oder irgendwie dynamisch laden?
         ######
 
-        DC_diode_var_s, DC_diode_var_p, DC_detector_var_s, DC_detector_var_p = get_var_data()
-
-        # def compute_reference_diode(reference_data):
-        #     # nutze self.mean_dc_diode[-1], um den Dunkelstromwert bei 2Theta=180° auszulesen
-        #     denominator = -reference_data['Measur_Diod'].values - self.mean_dc_diode[-1]
-        #     numerator = -reference_data['Measur_Detec'].values - self.mean_dc_detect[-1]
-        #     if np.any(denominator == 0):
-        #         raise ValueError("Division durch Null in compute_reference_diode!")
-        #     return numerator / denominator
+        # if not self.dark_current_var_indices:
+        #     raise ValueError("Keine Indizes für Dunkelstrom-Variation gesetzt!")
+        #
+        # # DC_diode_var_s, DC_diode_var_p, DC_detector_var_s, DC_detector_var_p = get_var_data()
+        # # indizes_dark_var = [25, 27, 29, 31, 33, 35]
+        #
+        # DC_diode_var_p, DC_detector_var_p = self.compute_dark_current_variation(self.dark_current_var_indices)
+
+        if not self.dark_current_var_indices:
+            # Fallback: vorab berechnete, hardcodierte Varianz-Daten laden
+            DC_diode_var_s, DC_diode_var_p, DC_detector_var_s, DC_detector_var_p = get_var_data()
+            DC_diode_var = DC_diode_var_p
+            DC_detector_var = DC_detector_var_p
+            print('Use hard coded values of TiN for dark current variation')
+        else:
+            # Variationen basierend auf übergebenen Indizes neu berechnen
+            DC_diode_var, DC_detector_var = self.compute_dark_current_variation(self.dark_current_var_indices)
+            # Optional: Werte für "s" auf "p" setzen, falls keine getrennte Berechnung erfolgt
+            # DC_diode_var_s = DC_diode_var_p
+            # DC_detector_var_s = DC_detector_var_p
 
         # 🔹 Reshape der Dunkelstrom-Durchschnitte (auf Anzahl der Thetawerte)
         self.mean_dc_detect = self.mean_dc_detect[:len(self.theta)]
@@ -414,8 +424,8 @@ class ReflectanceProcessor:
         #####
 
         # 🔹 Reflektanzdaten berechnen mit korrekt skalierten Dunkelstromwerten
-        # TODO Warum hier nur p?
-        reflectance_data_var = (Ref_detect - DC_detector_var_p) / (Ref_diode - DC_diode_var_p)
+        # reflectance_data_var = (Ref_detect - DC_detector_var_p) / (Ref_diode - DC_diode_var_p)
+        reflectance_data_var = (Ref_detect - DC_detector_var) / (Ref_diode - DC_diode_var)
         # TODO
         reflectance_data = (Ref_detect - mean_dc_detect_expanded) / \
                            (Ref_diode - mean_dc_diode_expanded)
@@ -440,11 +450,14 @@ class ReflectanceProcessor:
         # könnte man sicher umstrukturieren
         #####
 
-        Reference_diode_before_var = compute_reference_diode(self.data_ref_before, DC_diode_var_p,
-                                                        DC_detector_var_p)
-        Reference_diode_after_var = compute_reference_diode(self.data_ref_after, DC_diode_var_p,
-                                                        DC_detector_var_p)
-
+        # Reference_diode_before_var = compute_reference_diode(self.data_ref_before, DC_diode_var_p,
+        #                                                 DC_detector_var_p)
+        # Reference_diode_after_var = compute_reference_diode(self.data_ref_after, DC_diode_var_p,
+        #                                                 DC_detector_var_p)
+        Reference_diode_before_var = compute_reference_diode(self.data_ref_before, DC_diode_var,
+                                                             DC_detector_var)
+        Reference_diode_after_var = compute_reference_diode(self.data_ref_after, DC_diode_var,
+                                                            DC_detector_var)
         #####
         # I0 NUR nötig für die Unsicherheiten
         # zusammen mit den ..._var packen?
@@ -482,7 +495,7 @@ class ReflectanceProcessor:
 
         # uncert_rel, uncert_abs = calculate_uncertainty(I0_drift, RT, Reflectance_Diode)
         uncert_rel, uncert_abs = calculate_uncertainty(self.data_meas['Wavlen'].round().astype(int).unique(), I0_drift, RT, Reflectance_Diode)
-        print('wl, unc_rel',self.data_meas['Wavlen'].round().astype(int).unique()[17], uncert_rel[17])
+        # print('wl, unc_rel',self.data_meas['Wavlen'].round().astype(int).unique()[17], uncert_rel[17])
 
         return self.wavelen, self.theta, Reflectance_Diode, Reference_diode_before, Reference_diode_after, uncert_rel, uncert_abs