Improve fit

• Estimate peak using a simple Gaussian filter and np.argmax to speed up fitting.
• Calculate outliers with IsolationForest and per channel.
• Add possibility of a polynomial fit on pes and sqrt(pes).
• Add possibility of fitting with weights and function to get weights based on the XGM intensity to get a flatter response a as a function of intensity.

pes_to_spec/test/offline_analysis.py

@@ -8,7 +8,7 @@ import os
 import argparse
 
 import numpy as np
-from extra_data import open_run, by_id
+from extra_data import open_run, by_id, RunDirectory
 from pes_to_spec.model import Model, matching_two_ids
 
 from itertools import product
@@ -126,19 +126,21 @@ def main():
     Main entry point. Reads some data, trains and predicts.
     """
     parser = argparse.ArgumentParser(prog="offline_analysis", description="Test pes2spec doing an offline analysis of the data.")
-    parser.add_argument('-p', '--proposal', type=int, meta='INT', help='Proposal number', default=2828)
-    parser.add_argument('-r', '--run', type=int, meta='INT', help='Run number', default=206)
-    parser.add_argument('-m', '--model', type=str, meta='FILENAME', default="", help='Model to load. If given, do not train a model and just do inference with this one.')
-    parser.add_argument('-d', '--directory', type=str, meta='DIRECTORY', default=".", help='Where to save the results.')
-    parser.add_argument('-S', '--spec', type=str, meta='NAME', default="SA3_XTD10_SPECT/MDL/SPECTROMETER_SQS_NAVITAR:output", help='SPEC name')
-    parser.add_argument('-P', '--pes', type=str, meta='NAME', default="SA3_XTD10_PES/ADC/1:network", help='PES name')
-    parser.add_argument('-X', '--xgm', type=str, meta='NAME', default="SA3_XTD10_XGM/XGM/DOOCS:output", help='XGM name')
-    parser.add_argument('-o', '--offset', type=int, meta='INT', default=0, help='Train ID offset')
+    parser.add_argument('-p', '--proposal', type=int, metavar='INT', help='Proposal number', default=2828)
+    parser.add_argument('-r', '--run', type=int, metavar='INT', help='Run number', default=206)
+    parser.add_argument('-m', '--model', type=str, metavar='FILENAME', default="", help='Model to load. If given, do not train a model and just do inference with this one.')
+    parser.add_argument('-d', '--directory', type=str, metavar='DIRECTORY', default=".", help='Where to save the results.')
+    parser.add_argument('-S', '--spec', type=str, metavar='NAME', default="SA3_XTD10_SPECT/MDL/SPECTROMETER_SQS_NAVITAR:output", help='SPEC name')
+    parser.add_argument('-P', '--pes', type=str, metavar='NAME', default="SA3_XTD10_PES/ADC/1:network", help='PES name')
+    parser.add_argument('-X', '--xgm', type=str, metavar='NAME', default="SA3_XTD10_XGM/XGM/DOOCS:output", help='XGM name')
+    parser.add_argument('-o', '--offset', type=int, metavar='INT', default=0, help='Train ID offset')
 
     args = parser.parse_args()
 
+    print("Opening run ...")
     # get run
     run = open_run(proposal=args.proposal, run=args.run)
+    #run = RunDirectory("/gpfs/exfel/data/scratch/tmichela/data/r0206")
 
     # ----------------Used in the first tests-------------------------
     # get train IDs and match them, so we are sure to have information from all needed sources
@@ -184,6 +186,8 @@ def main():
                for ch in channels}
     pes_raw_t = {ch: run[pes_name, f"digitizers.{ch}.raw.samples"].select_trains(by_id[test_tids]).ndarray()
                for ch in channels}
+    
+    print("Data in memory.")
 
     # read the XGM information
     #xgm_pressure = run['SA3_XTD10_XGM/XGM/DOOCS', "pressure.pressureFiltered.value"].select_trains(by_id[tids]).ndarray()
@@ -198,7 +202,7 @@ def main():
 
     train_idx = np.isin(tids, train_tids)
 
-    model.debug_peak_finding(pes_raw, os.path.join(args.directory, "test_peak_finding.png"))
+    #model.debug_peak_finding(pes_raw, os.path.join(args.directory, "test_peak_finding.png"))
     if len(args.model) == 0:
         print("Fitting")
         start = time_ns()
@@ -230,7 +234,7 @@ def main():
     rmse = model.check_compatibility(pes_raw_t)
     print("Consistency check RMSE ratios:", rmse)
     rmse = model.check_compatibility_per_channel(pes_raw_t)
-    print("Consistency per channel check RMSE ratios:", rmse)
+    print("Consistency per channel check (chi2 - icdf(p=0.05))/ndof:", rmse)
     t += [time_ns() - start]
     t_names += ["Consistency"]
 
@@ -253,7 +257,7 @@ def main():
         spec_smooth = model.preprocess_high_res(spec_raw_int)
     first, last = model.get_low_resolution_range()
     first += 10
-    last -= 100
+    last -= 10
     pes_to_show = 'channel_1_D'
     # plot
     for tid in test_tids: