diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 7c2031968c..9756b92634 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/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""" diff --git a/lstchain/calib/camera/pedestals.py b/lstchain/calib/camera/pedestals.py index 52b25234ae..7e8aac5542 100644 --- a/lstchain/calib/camera/pedestals.py +++ b/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 diff --git a/lstchain/calib/camera/r0.py b/lstchain/calib/camera/r0.py index 1cd91a6610..3c3ba7792c 100644 --- a/lstchain/calib/camera/r0.py +++ b/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. diff --git a/lstchain/scripts/lstchain_data_muon_analysis.py b/lstchain/scripts/lstchain_data_muon_analysis.py index a3584bd6bc..dd123dbd8b 100644 --- a/lstchain/scripts/lstchain_data_muon_analysis.py +++ b/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, diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index ec9d064383..eec17d5481 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -29,6 +29,10 @@ 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 @@ -36,6 +40,10 @@ 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") diff --git a/lstchain/scripts/onsite/onsite_create_drs4_pedestal_file.py b/lstchain/scripts/onsite/onsite_create_drs4_pedestal_file.py index dcc98c0345..b6a50224ce 100755 --- a/lstchain/scripts/onsite/onsite_create_drs4_pedestal_file.py +++ b/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}") diff --git a/lstchain/tools/lstchain_data_create_calibration_file.py b/lstchain/tools/lstchain_data_create_calibration_file.py index c01cb66ac9..2633de195a 100644 --- a/lstchain/tools/lstchain_data_create_calibration_file.py +++ b/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 ( diff --git a/lstchain/visualization/plot_calib.py b/lstchain/visualization/plot_calib.py index dab88250f2..09bff9a260 100644 --- a/lstchain/visualization/plot_calib.py +++ b/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,7 +63,7 @@ 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": @@ -70,6 +71,7 @@ def plot_all(ped_data, ff_data, calib_data, run=0, plot_file="none"): 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) diff --git a/lstchain/visualization/plot_drs4.py b/lstchain/visualization/plot_drs4.py index 08991d0e11..c7a660d653 100644 --- a/lstchain/visualization/plot_drs4.py +++ b/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]) diff --git a/notebooks/Analyze_real_muon_data.ipynb b/notebooks/Analyze_real_muon_data.ipynb index 9e524403b9..37850fbddf 100644 --- a/notebooks/Analyze_real_muon_data.ipynb +++ b/notebooks/Analyze_real_muon_data.ipynb @@ -983,7 +983,6 @@ { "ename": "KeyboardInterrupt", "evalue": "", - "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", @@ -991,7 +990,8 @@ "\u001b[0;32m/Users/rubenlopez/anaconda/envs/cta-dev/lib/python3.6/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36mamax\u001b[0;34m(a, axis, out, keepdims, initial)\u001b[0m\n\u001b[1;32m 2332\u001b[0m \"\"\"\n\u001b[1;32m 2333\u001b[0m return _wrapreduction(a, np.maximum, 'max', axis, None, out, keepdims=keepdims,\n\u001b[0;32m-> 2334\u001b[0;31m initial=initial)\n\u001b[0m\u001b[1;32m 2335\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2336\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/Users/rubenlopez/anaconda/envs/cta-dev/lib/python3.6/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36m_wrapreduction\u001b[0;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_wrapreduction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mufunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 65\u001b[0m \u001b[0mpasskwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 66\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 67\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_NoValue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0mpasskwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mKeyboardInterrupt\u001b[0m: " - ] + ], + "output_type": "error" } ], "source": [ @@ -1402,10 +1402,10 @@ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 28)", - "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m28\u001b[0m\n\u001b[0;31m print(f\"Event {event.lst.tel[0].evt.event_id}, Max: {np.max(std_signal)} counts\")\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" - ] + ], + "output_type": "error" } ], "source": [ @@ -1471,13 +1471,13 @@ { "ename": "NameError", "evalue": "name 'np' is not defined", - "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;31m# If you want to make a movie with all the slices\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mstd_signal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1855\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mpixel\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mN_modules\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mstd_signal\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mpixel\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mevent\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtel\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwaveform\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpixel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m38\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" - ] + ], + "output_type": "error" } ], "source": [