Merge branch 'master' into update_dl3

environment.yml

@@ -16,6 +16,7 @@ dependencies:
   - joblib
   - toml
   - protozfits=2.0
+  - pyparsing~=2.4
   - pip:
       - pytest_runner
       - pytest-ordering

lstchain/calib/camera/calibration_calculator.py

@@ -3,9 +3,10 @@
 """
 
 import numpy as np
+import h5py
 from ctapipe.core import Component
 from ctapipe.core import traits
-from ctapipe.core.traits import  Float, List
+from ctapipe.core.traits import  Float, List, Path
 from lstchain.calib.camera.flatfield import FlatFieldCalculator
 from lstchain.calib.camera.pedestals import PedestalCalculator
 from ctapipe.containers import EventType
@@ -35,6 +36,12 @@ class CalibrationCalculator(Component):
     kwargs
 
     """
+
+    systematic_correction_path = Path(
+        exists=True, directory_ok=False,
+        help='Path to systematic correction file ',
+    ).tag(config=True)
+
     squared_excess_noise_factor = Float(
         1.222,
         help='Excess noise factor squared: 1+ Var(gain)/Mean(Gain)**2'
@@ -86,6 +93,10 @@ def __init__(
 
         super().__init__(parent=parent, config=config,**kwargs)
 
+        if self.squared_excess_noise_factor<=0:
+            msg="Argument squared_excess_noise_factor must have a positive value"
+            raise ValueError(msg)
+
         self.flatfield = FlatFieldCalculator.from_name(
             self.flatfield_product,
             parent=self,
@@ -103,6 +114,16 @@ def __init__(
 
         self.tel_id = self.flatfield.tel_id
 
+        # load systematic correction term B
+        self.quadratic_term = 0
+        if self.systematic_correction_path is not None:
+            try:
+                with h5py.File(self.systematic_correction_path, 'r') as hf:
+                    self.quadratic_term = np.array(hf['B_term'])
+
+            except:
+                raise IOError(f"Problem in reading quadratic term file {self.systematic_correction_path}")
+
         self.log.debug(f"{self.pedestal}")
         self.log.debug(f"{self.flatfield}")
 
@@ -138,12 +159,22 @@ def calculate_calibration_coefficients(self, event):
         calib_data.unusable_pixels = np.logical_or(monitoring_unusable_pixels,
                                                    status_data.hardware_failing_pixels)
 
+        signal = ff_data.charge_median - ped_data.charge_median
+
         # Extract calibration coefficients with F-factor method
         # Assume fixed excess noise factor must be known from elsewhere
+        numerator = ff_data.charge_std ** 2 - ped_data.charge_std ** 2
+        denominator = self.squared_excess_noise_factor  * signal
+        gain = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator != 0)
+
+        # correct for the quadratic term (which is zero if not given)
+        systematic_correction = self.quadratic_term**2 * signal / self.squared_excess_noise_factor
+        gain -= systematic_correction
 
         # calculate photon-electrons
-        numerator = self.squared_excess_noise_factor  * (ff_data.charge_median - ped_data.charge_median) ** 2
-        denominator = ff_data.charge_std ** 2 - ped_data.charge_std ** 2
+        numerator = signal
+        denominator = gain
+
         n_pe = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator != 0)
 
         # fill WaveformCalibrationContainer
@@ -152,21 +183,14 @@ def calculate_calibration_coefficients(self, event):
         calib_data.time_max = ff_data.sample_time_max
         calib_data.n_pe = n_pe
 
-        # find signal median of good pixels
+        # find signal median of good pixels over the camera (FF factor=<npe>/npe)
         masked_npe = np.ma.array(n_pe, mask=calib_data.unusable_pixels)
         npe_signal_median = np.ma.median(masked_npe, axis=1)
 
-        # Flat field factor
-        numerator = npe_signal_median[:, np.newaxis]
-        denominator = n_pe
-        ff = np.divide(numerator, denominator, out=np.zeros_like(denominator), where=denominator != 0)
-
-        # calibration coefficients
-        numerator = n_pe * ff
-
-        # correct the signal for the integration window
-        denominator = (ff_data.charge_median - ped_data.charge_median) 
-        calib_data.dc_to_pe = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator != 0)
+        # flat-fielded calibration coefficients
+        numerator = npe_signal_median[:,np.newaxis]
+        denominator = signal
+        calib_data.dc_to_pe = np.divide(numerator, denominator, out=np.zeros_like(denominator), where=denominator != 0)
 
         # flat-field time corrections
         calib_data.time_correction = -ff_data.relative_time_median

lstchain/calib/camera/time_sampling_correction.py

@@ -3,7 +3,7 @@
 
 
 from ctapipe.core import Component
-from ctapipe.core.traits import Unicode
+from ctapipe.core.traits import Path
 
 __all__ = [
     'TimeSamplingCorrection',
@@ -17,13 +17,11 @@ class TimeSamplingCorrection(Component):
         using Fourier series expansion.
     """
 
-    time_sampling_correction_path = Unicode(
-        '',
-        help='Path to the waveform sampling correction file',
-        allow_none = True,
+    time_sampling_correction_path = Path(
+        exists=True, directory_ok=False,
+        help='Path to time sampling correction file'
     ).tag(config=True)
 
-
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 

lstchain/data/filters_transmission.dat

@@ -0,0 +1,37 @@
+# filters, transmission estimated from camera data (June 2021, Y. Kobayashi)
+11 1.00000000
+12 0.30000000
+13 0.14000000
+14 0.07300000
+15 0.01800000
+16 0.00300000
+21 0.35000000
+22 0.10500000
+23 0.04900000
+24 0.02555000
+25 0.00630000
+26 0.00105000
+31 0.11000000
+32 0.03300000
+33 0.01540000
+34 0.00803000
+35 0.00198000
+36 0.00033000
+41 0.03000000
+42 0.00900000
+43 0.00420000
+44 0.00219000
+45 0.00054000
+46 0.00009000
+51 0.00300000
+52 0.00090000
+53 0.00042000
+54 0.00021900
+55 0.00005400
+56 0.00000900
+61 0.00050000
+62 0.00015000
+63 0.00007000
+64 0.00003650
+65 0.00000900
+66 0.00000150

lstchain/data/lstchain_fit_intensity_scan_config_example.json

@@ -0,0 +1,21 @@
+{
+  "FitIntensityScan":{
+
+      "log_level":"DEBUG",  
+      "signal_range": [[1500,14000],[200,14000]],				     
+      "run_list":[4753,4754,4755,4756,4757,4758,4759,4760,4761,4762,4763,4764,4765,4766,4766,4767,4769,4770,4771,4772],                 
+
+      "log_config": {
+       "loggers": {          
+
+          "ctapipe.io.hdf5tableio.HDF5TableReader":{
+             "level": "INFO",
+	     "handlers": ["console"]  
+          }	                
+        } 
+
+      }    
+    } 
+
+
+} 

lstchain/data/lstchain_src_dep_config.json

@@ -29,9 +29,9 @@
     "intensity": [0, Infinity],
     "width": [0, Infinity],
     "length": [0, Infinity],
-    "wl": [0, 1],
-    "r": [0, 1],
-    "leakage_intensity_width_2": [0, 1]
+    "wl": [0, Infinity],
+    "r": [0, Infinity],
+    "leakage_intensity_width_2": [0, Infinity]
   },
 
   "tailcuts_clean_with_pedestal_threshold": {
@@ -43,8 +43,32 @@
     "use_only_main_island":false,
     "delta_time": 2
   },
-
-    "random_forest_regressor_args": {
+
+  "dynamic_cleaning": {
+    "apply": true,
+    "threshold": 267,
+    "fraction_cleaning_intensity": 0.03
+  },
+
+  "random_forest_energy_regressor_args": {
+    "max_depth": 50,
+    "min_samples_leaf": 2,
+    "n_jobs": 4,
+    "n_estimators": 150,
+    "bootstrap": true,
+    "criterion": "mse",
+    "max_features": "auto",
+    "max_leaf_nodes": null,
+    "min_impurity_decrease": 0.0,
+    "min_samples_split": 2,
+    "min_weight_fraction_leaf": 0.0,
+    "oob_score": false,
+    "random_state": 42,
+    "verbose": 0,
+    "warm_start": false
+  },
+
+  "random_forest_disp_regressor_args": {
     "max_depth": 50,
     "min_samples_leaf": 2,
     "n_jobs": 4,
@@ -62,7 +86,7 @@
     "warm_start": false
   },
 
-  "random_forest_classifier_args": {
+  "random_forest_particle_classifier_args": {
     "max_depth": 100,
     "min_samples_leaf": 2,
     "n_jobs": 4,
@@ -81,7 +105,26 @@
     "class_weight": null
   },
 
-  "regression_features": [
+  "random_forest_disp_classifier_args": {
+    "max_depth": 100,
+    "min_samples_leaf": 2,
+    "n_jobs": 4,
+    "n_estimators": 100,
+    "criterion": "gini",
+    "min_samples_split": 2,
+    "min_weight_fraction_leaf": 0.0,
+    "max_features": "auto",
+    "max_leaf_nodes": null,
+    "min_impurity_decrease": 0.0,
+    "bootstrap": true,
+    "oob_score": false,
+    "random_state": 42,
+    "verbose": 0.0,
+    "warm_start": false,
+    "class_weight": null
+  },
+
+  "energy_regression_features": [
     "log_intensity",
     "width",
     "length",
@@ -93,7 +136,19 @@
     "dist"
   ],
 
-  "classification_features": [
+  "disp_regression_features": [
+    "log_intensity",
+    "width",
+    "length",
+    "wl",
+    "skewness_from_source",
+    "kurtosis",
+    "time_gradient_from_source",
+    "leakage_intensity_width_2",
+    "dist"
+  ],
+
+  "particle_classification_features": [
     "log_intensity",
     "width",
     "length",
@@ -106,6 +161,19 @@
     "dist"
   ],
 
+  "disp_classification_features": [
+    "log_intensity",
+    "width",
+    "length",
+    "wl",
+    "skewness_from_source",
+    "kurtosis",
+    "time_gradient_from_source",
+    "leakage_intensity_width_2",
+    "log_reco_energy",
+    "dist"
+  ],
+
   "allowed_tels": [1, 2, 3, 4],
   "max_events": null,
   "custom_calibration": false,
@@ -131,6 +199,8 @@
     },
   "timestamps_pointing":"ucts",
 
+  "train_gamma_src_r_deg": [0, Infinity],
+
   "source_dependent": true,
   "mc_nominal_source_x_deg": 0.4,
   "mc_nominal_source_y_deg": 0.0,