Merge pull request #263 from FrancaCassol/update_muon_script

Update muon script and calibration scripts

lstchain/calib/camera/flatfield.py

@@ -175,7 +175,6 @@ def store_results(self, event):
         # mask the part of the array not filled
         self.sample_masked_pixels[self.num_events_seen:] = 1
 
-
         relative_gain_results = self.calculate_relative_gain_results(
             self.charge_medians,
             self.charges,
@@ -202,8 +201,6 @@ def store_results(self, event):
         event.mon.tel[self.tel_id].pixel_status.flatfield_failing_pixels = \
             np.logical_or(ff_charge_failing_pixels, container.time_median_outliers)
 
-
-
     def setup_sample_buffers(self, waveform, sample_size):
         """Initialize sample buffers"""
 

lstchain/calib/camera/pedestals.py

@@ -166,7 +166,6 @@ def store_results(self, event):
         if self.num_events_seen == 0:
             raise ValueError("No pedestal events in statistics, zero results")
 
-
         container = event.mon.tel[self.tel_id].pedestal
 
         # mask the part of the array not filled

lstchain/calib/camera/r0.py

@@ -144,14 +144,15 @@ def calibrate(self, event):
             self.time_lapse_corr(event, tel_id)
             self.interpolate_spikes(event, tel_id)
 
-            event.r1.tel[self.tel_id].trigger_type = event.r0.tel[self.tel_id].trigger_type
+            event.r1.tel[tel_id].trigger_type = event.r0.tel[tel_id].trigger_type
 
-            event.r1.tel[self.tel_id].trigger_time = event.r0.tel[self.tel_id].trigger_time
+            event.r1.tel[tel_id].trigger_time = event.r0.tel[tel_id].trigger_time
 
             samples = event.r1.tel[tel_id].waveform[:, :, self.r1_sample_start:self.r1_sample_end]
 
             event.r1.tel[tel_id].waveform = samples.astype('int16') - self.offset
 
+
     def subtract_pedestal(self, event, tel_id):
         """
         Subtract cell offset using pedestal file.

lstchain/scripts/lstchain_data_muon_analysis.py

@@ -40,22 +40,21 @@
 parser.add_argument("--calibration_file", help = "Path to the file containing the calibration constants",
                     type = str)
 
-parser.add_argument("--time_calibration_file", help = "Path to the calibraiton file with time corrections ",
-                    type = str)
+parser.add_argument("--time_calibration_file", help = "Path to the calibration file with time corrections ",
+                    type = str, default='')
 
 # Optional argument
 parser.add_argument("--max_events", help = "Maximum numbers of events to read. "
                                          "Default = 100",
                     type = int, default = 100)
 parser.add_argument("--gain_threshold", help = "Gain threshold in ADC. "
-                                               "Default = 4094", type = int, default=4094)
+                                               "Default = 3500", type = int, default=3500)
 
 parser.add_argument("--plot_rings", help = "Plot figures of the stored rings", 
                     default = False, action='store_true')
 
 parser.add_argument("--plots_path", help = "Path to the plots",
                     default = None, type = str)
-
 parser.add_argument("--tel_id", help = "telescope id. "
                                        "Default = 1",type = int, default = 1)
 
@@ -69,10 +68,11 @@ def main():
     print("output file: {}".format(args.output_file))
     print("pedestal file: {}".format(args.pedestal_file))
     print("calibration file: {}".format(args.calibration_file))
+    print("time calibration file: {}".format(args.time_calibration_file))
     print("max events: {}".format(args.max_events))
 
     # Camera geometry
-    geom = CameraGeometry.from_name("LSTCam-003")
+    #geom = CameraGeometry.from_name("LSTCam-003")
 
     # Definition of the output parameters for the table
     output_parameters = {'event_id': [],
@@ -111,6 +111,7 @@ def main():
                                             gain_threshold=args.gain_threshold,
                                             allowed_tels=[tel_id])
 
+
     # Maximum number of events
     if args.max_events:
         max_events = args.max_events
@@ -121,9 +122,16 @@ def main():
     num_muons = 0
     source = event_source(input_url = args.input_file, max_events = max_events)
 
+    # geometry
+    subarray = source.subarray(tel_id)
+    telescope_description = subarray.tel[tel_id]
+    equivalent_focal_length = telescope_description.optics.equivalent_focal_length
+    mirror_area = telescope_description.optics.mirror_area.to("m2")
+    geom = telescope_description.camera
+
+
     for event in source:
         event_id = event.lst.tel[tel_id].evt.event_id
-        telescope_description = event.inst.subarray.tel[tel_id]
 
         # drs4 calibration
         r0calib.calibrate(event)
@@ -135,14 +143,13 @@ def main():
         if not tag_pix_thr(image): #default skipps pedestal and calibration events
             continue
 
-        if not muon_filter(image): #default values apply no filtering
+        if not muon_filter(image,800,5000): #default values apply no filtering
             continue
 
+
         print("--> Event {}. Number of pixels above 10 phe: {}".format(event_id,
                                                                    np.size(image[image > 10.])))
 
-        equivalent_focal_length = telescope_description.optics.equivalent_focal_length
-        mirror_area = telescope_description.optics.mirror_area.to("m2")
 
         muonintensityparam, size_outside_ring, muonringparam, good_ring = \
             analyze_muon_event(event_id, image, geom, equivalent_focal_length,

lstchain/scripts/onsite/onsite_create_calibration_file.py

@@ -29,13 +29,21 @@
 optional.add_argument('-s', '--statistics', help="Number of events for the flat-field and pedestal statistics",
                       type=int, default=10000)
 optional.add_argument('-b','--base_dir', help="Root dir for the output directory tree",type=str, default='/fefs/aswg/data/real')
+optional.add_argument('--default_time_run', help="If 0 time calibration is calculated otherwise create a link to the give run time calibration",type=int, default='1625')
+optional.add_argument('--ff_calibration', help="Perform the charge calibration (yes/no)",type=str, default='yes')
+optional.add_argument('--tel_id', help="telescope id. Default = 1", type=int, default=1)
+
 
 args = parser.parse_args()
 run = args.run_number
 ped_run = args.pedestal_run
 prod_id = 'v%02d'%args.version
 stat_events = args.statistics
 base_dir = args.base_dir
+default_time_run = args.default_time_run
+ff_calibration = args.ff_calibration
+tel_id = args.tel_id
+
 max_events = 1000000
 
 
@@ -69,18 +77,21 @@ def main():
             print(f">>> Error: The pedestal file {pedestal_file} do not exist.\n Exit")
             exit(0)
 
+
         #
         # produce ff calibration file
         #
 
+
         # define charge file names
         output_file = f"{output_dir}/calibration.Run{run}.0000.hdf5"
         log_file = f"{output_dir}/log/calibration.Run{run}.0000.log"
         print(f"\n--> Output file {output_file}")
-        if os.path.exists(output_file):
+        if os.path.exists(output_file) and ff_calibration is 'yes':
             if query_yes_no(">>> Output file exists already. Do you want to remove it?"):
                 os.remove(output_file)
             else:
+                print(f"\n--> Stop")
                 exit(1)
 
         print(f"\n--> Log file {log_file}")
@@ -92,31 +103,49 @@ def main():
             exit(1)
         print(f"\n--> Config file {config_file}")
 
-        # run lstchain script
-        cmd = f"lstchain_data_create_calibration_file " \
-              f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \
-              f"--FlatFieldCalculator.sample_size={stat_events} --PedestalCalculator.sample_size={stat_events}  " \
-              f"--EventSource.max_events={max_events} --config={config_file}  >  {log_file} 2>&1"
+        if ff_calibration is 'yes':
+            # run lstchain script
+            cmd = f"lstchain_data_create_calibration_file " \
+                  f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \
+                  f"--FlatFieldCalculator.sample_size={stat_events} --PedestalCalculator.sample_size={stat_events}  " \
+                  f"--EventSource.max_events={max_events} --config={config_file}  >  {log_file} 2>&1"
 
-        print("\n--> RUNNING...")
-        os.system(cmd)
+            print("\n--> RUNNING...")
+            os.system(cmd)
 
-        # plot and save some results
-        plot_file=f"{output_dir}/log/calibration.Run{run}.0000.pdf"
-        print(f"\n--> PRODUCING PLOTS in {plot_file} ...")
-        calib.read_file(output_file)
-        calib.plot_all(calib.ped_data, calib.ff_data, calib.calib_data, run, plot_file)
+            # plot and save some results
+            plot_file=f"{output_dir}/log/calibration.Run{run}.0000.pedestal.Run{ped_run}.0000.pdf"
+            print(f"\n--> PRODUCING PLOTS in {plot_file} ...")
+            calib.read_file(output_file,tel_id)
+            calib.plot_all(calib.ped_data, calib.ff_data, calib.calib_data, run, plot_file)
 
         #
         # produce drs4 time calibration file
         #
-
         time_file = f"{output_dir}/time_calibration.Run{run}.0000.hdf5"
-        print(f"\n--> PRODUCING TIME CALIBRATION in {time_file} ...")
-        cmd = f"lstchain_data_create_time_calibration_file  --input_file {input_file} " \
-              f"--output_file {time_file} -conf {config_file} -ped {pedestal_file} 2>&1"
-        print("\n--> RUNNING...")
-        os.system(cmd)
+
+        if default_time_run is 0:
+            print(f"\n--> PRODUCING TIME CALIBRATION in {time_file} ...")
+            cmd = f"lstchain_data_create_time_calibration_file  --input_file {input_file} " \
+                  f"--output_file {time_file} -conf {config_file} -ped {pedestal_file} 2>&1"
+            print("\n--> RUNNING...")
+            os.system(cmd)
+        else:
+            # otherwise perform a link to the default time calibration file
+            print(f"\n--> PRODUCING LINK TO DEFAULT TIME CALIBRATION (run {default_time_run})")
+            file_list = sorted(Path(f"{base_dir}/calibration/").rglob(f'*/{prod_id}/time_calibration.Run{default_time_run}*'))
+
+
+            if len(file_list) == 0:
+                print(f">>> Error: time calibration file for run {default_time_run} not found\n")
+                raise NameError()
+            else:
+                time_calibration_file = file_list[0]
+                input_dir, name = os.path.split(os.path.abspath(time_calibration_file ))
+                cmd=f"ln -sf {time_calibration_file} {time_file}"
+                os.system(cmd)
+
+        print(f"\n--> Time calibration file: {time_file}")
 
         print("\n--> END")
 

lstchain/scripts/onsite/onsite_create_drs4_pedestal_file.py

@@ -26,12 +26,16 @@
                       type=int, default=8000)
 optional.add_argument('-b','--base_dir', help="Base dir for the output directory tree",
                       type=str, default='/fefs/aswg/data/real')
+optional.add_argument('--tel_id', help="telescope id. Default = 1",
+                      type=int, default=1)
+
 
 args = parser.parse_args()
 run = args.run_number
 prod_id = 'v%02d'%args.prod_version
 max_events = args.max_events
 base_dir = args.base_dir
+tel_id = args.tel_id
 
 
 def main():
@@ -81,12 +85,10 @@ def main():
         # plot and save some results
         plot_file=f"{output_dir}/log/drs4_pedestal.Run{run}.0000.pdf"
         print(f"\n--> PRODUCING PLOTS in {plot_file} ...")
-        drs4.plot_pedestals(input_file, output_file, run, plot_file, tel_id=1, offset_value=300)
+        drs4.plot_pedestals(input_file, output_file, run, plot_file, tel_id=tel_id, offset_value=300)
 
         print("\n--> END")
 
-
-
     except Exception as e:
         print(f"\n >>> Exception: {e}")
 

lstchain/tools/lstchain_data_create_calibration_file.py

@@ -195,6 +195,7 @@ def start(self):
                 if event.r1.tel[self.tel_id].trigger_type == 32:
                     new_ped = self.pedestal.calculate_pedestals(event)
 
+
                 # if flat-field event: no calibration  TIB for the moment,
                 # use a cut on the charge for ff events and on std for rejecting Magic Lidar events
                 elif event.r1.tel[self.tel_id].trigger_type == 4 or (

lstchain/visualization/plot_calib.py

@@ -24,6 +24,7 @@
 
 plot_dir = "none"
 
+
 def read_file(file_name, tel_id=1):
     """
      read camera calibration quantities
@@ -62,14 +63,15 @@ def plot_all(ped_data, ff_data, calib_data, run=0, plot_file="none"):
      calib_data: calibration container WaveformCalibrationContainer()
 
      """
-    camera = CameraGeometry.from_name("LSTCam", 2)
+    camera = CameraGeometry.from_name("LSTCam", 4)
 
     # plot open pdf
     if plot_file != "none":
         pp = PdfPages(plot_file)
 
     plt.rc('font', size=15)
 
+
     ### first figure
     fig = plt.figure(1, figsize=(12, 24))
     plt.tight_layout()
@@ -263,7 +265,7 @@ def plot_all(ped_data, ff_data, calib_data, run=0, plot_file="none"):
         # select good pixels
         select = np.logical_not(mask[chan])
         fig = plt.figure(chan+10, figsize=(12, 18))
-        fig.tight_layout(rect=[0, 0.03, 1, 0.95])
+        fig.tight_layout(rect=[0, 0.0, 1, 0.95])
 
         fig.suptitle(f"Run {run} channel: {channel[chan]}", fontsize=25)
 

lstchain/visualization/plot_drs4.py

@@ -5,7 +5,6 @@
 from matplotlib.backends.backend_pdf import PdfPages
 from ctapipe.io import event_source
 from lstchain.calib.camera.r0 import LSTR0Corrections
-from ctapipe.io.containers import PedestalContainer
 from ctapipe.instrument import  CameraGeometry
 
 
@@ -31,6 +30,10 @@ def plot_pedestals(data_file, pedestal_file, run=0 , plot_file="none", tel_id=1,
      run: run number of data to be corrected
 
      plot_file:  name of output pdf file
+
+     tel_id: id of the telescope
+
+     offset_value: baseline off_set
      """
 
     # plot open pdf
@@ -65,6 +68,8 @@ def plot_pedestals(data_file, pedestal_file, run=0 , plot_file="none", tel_id=1,
                                   config=charge_config)
 
     for i, event in enumerate(reader):
+        if tel_id != event.r0.tels_with_data[0]:
+            raise Exception(f"Given wrong telescope id {tel_id}, files has id {event.r0.tels_with_data[0]}")
 
         # move from R0 to R1
         r0_calib.calibrate(event)
@@ -189,7 +194,6 @@ def plot_pedestals(data_file, pedestal_file, run=0 , plot_file="none", tel_id=1,
             plt.plot([0, 40], [offset_value, offset_value], 'k--',  label="offset")
             plt.xlabel("time sample [ns]")
             plt.ylabel("counts [ADC]")
-
             plt.legend()
             plt.ylim([-50, 500])