__init__.py 80.1 KB
Newer Older
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
1
"""
2
3
# Event source for MAGIC calibrated data files.
# Requires uproot package (https://github.com/scikit-hep/uproot).
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
4
5
"""

Moritz Huetten's avatar
Moritz Huetten committed
6
import logging
7
8
9

import glob
import re
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
10
import os.path
11
from pathlib import Path
12
from enum import Enum, auto
13
14

import numpy as np
15
16

import scipy
17
18
19
20
import scipy.interpolate

from astropy import units as u
from astropy.time import Time
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
21

22
from ctapipe.io.eventsource import EventSource
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
23
24
25
26
27
from ctapipe.io.datalevels import DataLevel

from ctapipe.core import Container
from ctapipe.core import Field

Alessio Berti's avatar
Alessio Berti committed
28
29
30
from ctapipe.containers import (
    ArrayEventContainer,
    SimulatedShowerContainer,
31
    SimulationConfigContainer,
Alessio Berti's avatar
Alessio Berti committed
32
33
34
35
36
37
38
39
40
41
42
43
44
    PointingContainer,
    TelescopePointingContainer,
    TelescopeTriggerContainer,
    MonitoringCameraContainer,
    PedestalContainer,
)

from ctapipe.instrument import (
    TelescopeDescription,
    SubarrayDescription,
    OpticsDescription,
    CameraDescription,
)
45
46
47

__all__ = ['MAGICEventSource']

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
48
LOGGER = logging.getLogger(__name__)
49

Ievgen Vovk's avatar
Ievgen Vovk committed
50

51
# MAGIC telescope positions in m wrt. to the center of CTA simulations
52
53
54
55
56
57
#MAGIC_TEL_POSITIONS = {
#    1: [-27.24, -146.66, 50.00] * u.m,
#    2: [-96.44, -96.77, 51.00] * u.m
#}

# MAGIC telescope positions in m wrt. to the center of MAGIC simulations, from reflector input card
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
58
MAGIC_TEL_POSITIONS = {
59
60
    1: [31.80, -28.10, 0.00] * u.m,
    2: [-31.80, 28.10, 0.00] * u.m
61
}
62

63
# MAGIC telescope description
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
64
OPTICS = OpticsDescription.from_name('MAGIC')
65
66
MAGICCAM = CameraDescription.from_name("MAGICCam")
GEOM = MAGICCAM.geometry
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
67
MAGIC_TEL_DESCRIPTION = TelescopeDescription(
68
    name='MAGIC', tel_type='MAGIC', optics=OPTICS, camera=MAGICCAM)
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
69
MAGIC_TEL_DESCRIPTIONS = {1: MAGIC_TEL_DESCRIPTION, 2: MAGIC_TEL_DESCRIPTION}
70

71
# trigger patterns:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
72
73
74
75
MC_TRIGGER_PATTERN = 1
PEDESTAL_TRIGGER_PATTERN = 8
DATA_TRIGGER_PATTERN = 128

76
77
78
79
80
81
82
83
class MARSDataLevel(Enum):
    """Enum of the different MARS Data Levels"""

    CALIBRATED = auto()  # Calibrated images in charge and time (no waveforms)
    STAR       = auto()  # Cleaned images, with Hillas parametrization
    SUPERSTAR  = auto()  # Stereo parameters reconstructed
    MELIBEA    = auto()  # Reconstruction of hadronness, event direction and energy

84
85
86
87
88
89
90
class L3JumpError(Exception):
    """
    Exception raised when L3 trigger number jumps backward.
    """

    def __init__(self, message):
        self.message = message
91

92
93
94
95
96
97
98
99
class MissingDriveReportError(Exception):
    """
    Exception raised when a subrun does not have drive reports.
    """

    def __init__(self, message):
        self.message = message

100
101
102
103
104
105
106
107
108
class MAGICEventSource(EventSource):
    """
    EventSource for MAGIC calibrated data.

    This class operates with the MAGIC data run-wise. This means that the files
    corresponding to the same data run are loaded and processed together.
    """
    _count = 0

Alessio Berti's avatar
Alessio Berti committed
109
    def __init__(self, input_url=None, config=None, parent=None, **kwargs):
110
111
112
113
114
115
116
117
118
        """
        Constructor

        Parameters
        ----------
        config: traitlets.loader.Config
            Configuration specified by config file or cmdline arguments.
            Used to set traitlet values.
            Set to None if no configuration to pass.
119
        parent : ctapipe.core.Tool
120
121
122
123
124
125
126
127
128
            Tool executable that is calling this component.
            Passes the correct logger to the component.
            Set to None if no Tool to pass.
        kwargs: dict
            Additional parameters to be passed.
            NOTE: The file mask of the data to read can be passed with
            the 'input_url' parameter.
        """

129
        super().__init__(input_url=input_url, config=config, parent=parent, **kwargs)
130
131

        # Retrieving the list of run numbers corresponding to the data files
Alessio Berti's avatar
Alessio Berti committed
132
        run_info = self.get_run_info_from_name(str(input_url))
133
134
        run_numbers = run_info[0]
        is_mc_runs = run_info[1]
135
        telescope = run_info[2]
Moritz Huetten's avatar
Moritz Huetten committed
136

Alessio Berti's avatar
Alessio Berti committed
137
138
        self.run_numbers = run_numbers
        self.is_mc = is_mc_runs
139
        self.telescope = telescope
140

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
141
142
        # Retrieving the data level (so far HARDCODED Sorcerer)
        self.datalevel = DataLevel.DL1_IMAGES
143
        self.mars_datalevel = run_info[3]
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
144

145
146
147
        if self.is_mc:
            self.simulation_config = self.parse_simulation_header()

148
149
150
151
        # # Setting up the current run with the first run present in the data
        # self.current_run = self._set_active_run(run_number=0)
        self.current_run = None

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
152
153
        self._subarray_info = SubarrayDescription(
            'MAGIC', MAGIC_TEL_POSITIONS, MAGIC_TEL_DESCRIPTIONS)
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175

    @staticmethod
    def is_compatible(file_mask):
        """
        This method checks if the specified file mask corresponds
        to MAGIC data files. The result will be True only if all
        the files are of ROOT format and contain an 'Events' tree.

        Parameters
        ----------
        file_mask: str
            A file mask to check

        Returns
        -------
        bool:
            True if the masked files are MAGIC data runs, False otherwise.

        """

        is_magic_root_file = True

Alessio Berti's avatar
Alessio Berti committed
176
        file_list = glob.glob(str(file_mask))
177
178
179

        for file_path in file_list:
            try:
Alessio Berti's avatar
Alessio Berti committed
180
                import uproot
181
182
183
184
185
186
187
188
189
190

                try:
                    with uproot.open(file_path) as input_data:
                        if 'Events' not in input_data:
                            is_magic_root_file = False
                except ValueError:
                    # uproot raises ValueError if the file is not a ROOT file
                    is_magic_root_file = False

            except ImportError:
Ievgen Vovk's avatar
Ievgen Vovk committed
191
                if re.match(r'.+_m\d_.+root', file_path.lower()) is None:
192
193
194
195
196
                    is_magic_root_file = False

        return is_magic_root_file

    @staticmethod
197
    def get_run_info_from_name(file_name):
198
        """
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
199
        This internal method extracts the run number and
Moritz Huetten's avatar
Moritz Huetten committed
200
        type (data/MC) from the specified file name.
201
202
203

        Parameters
        ----------
204
        file_name : str
205
206
207
208
            A file name to process.

        Returns
        -------
209
210
211
212
213
214
215
216
217
218
219
220
221
        run_number: int
            The run number of the file.
        is_mc: Bool
            Flag to tag MC files
        telescope: int
            Number of the telescope
        datalevel: MARSDataLevel
            Data level according to MARS

        Raises
        ------
        IndexError
            Description
222
223
        """

224
225
226
227
228
229
230
231
        mask_data_calibrated = r"\d{6}_M(\d+)_(\d+)\.\d+_Y_.*"
        mask_data_star       = r"\d{6}_M(\d+)_(\d+)\.\d+_I_.*"
        mask_data_superstar  = r"\d{6}_(\d+)_S_.*"
        mask_data_melibea    = r"\d{6}_(\d+)_Q_.*"
        mask_mc_calibrated   = r"GA_M(\d)_za\d+to\d+_\d_(\d+)_Y_.*"
        mask_mc_star         = r"GA_M(\d)_za\d+to\d+_\d_(\d+)_I_.*"
        mask_mc_superstar    = r"GA_za\d+to\d+_\d_S_.*"
        mask_mc_melibea      = r"GA_za\d+to\d+_\d_Q_.*"
232
        mask_mc_alt = r".*_M(\d)_\d_(\d+)_.*"
233
234
235
236
237
        if re.findall(mask_data_calibrated, file_name):
            parsed_info = re.findall(mask_data_calibrated, file_name)
            telescope  = int(parsed_info[0][0])
            run_number = int(parsed_info[0][1])
            datalevel  = MARSDataLevel.CALIBRATED
Moritz Huetten's avatar
Moritz Huetten committed
238
            is_mc = False
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
        elif re.findall(mask_data_star, file_name):
            parsed_info = re.findall(mask_data_star, file_name)
            telescope  = int(parsed_info[0][0])
            run_number = int(parsed_info[0][1])
            datalevel  = MARSDataLevel.STAR
            is_mc = False
        elif re.findall(mask_data_superstar, file_name):
            parsed_info = re.findall(mask_data_superstar, file_name)
            telescope  = None
            run_number = int(parsed_info[0])
            datalevel  = MARSDataLevel.SUPERSTAR
            is_mc = False
        elif re.findall(mask_data_melibea, file_name):
            parsed_info = re.findall(mask_data_melibea, file_name)
            telescope  = None
            run_number = int(parsed_info[0])
            datalevel  = MARSDataLevel.MELIBEA
            is_mc = False
        elif re.findall(mask_mc_calibrated, file_name):
            parsed_info = re.findall(mask_mc_calibrated, file_name)
            telescope  = int(parsed_info[0][0])
            run_number = int(parsed_info[0][1])
            datalevel  = MARSDataLevel.CALIBRATED
Moritz Huetten's avatar
Moritz Huetten committed
262
            is_mc = True
263
264
        elif re.findall(mask_mc_star, file_name):
            parsed_info = re.findall(mask_mc_star, file_name)
265
266
            telescope  = int(parsed_info[0][0])
            run_number = int(parsed_info[0][1])
267
268
269
270
271
            datalevel  = MARSDataLevel.STAR
            is_mc = True
        elif re.findall(mask_mc_superstar, file_name):
            parsed_info = re.findall(mask_mc_superstar, file_name)
            telescope  = None
272
            run_number = None
273
274
275
276
277
            datalevel  = MARSDataLevel.SUPERSTAR
            is_mc = True
        elif re.findall(mask_mc_melibea, file_name):
            parsed_info = re.findall(mask_mc_melibea, file_name)
            telescope  = None
278
            run_number = None
279
280
281
            datalevel  = MARSDataLevel.MELIBEA
            is_mc = True
        else:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
282
283
284
            raise IndexError(
                'Can not identify the run number and type (data/MC) of the file'
                '{:s}'.format(file_name))
285

286
        return run_number, is_mc, telescope, datalevel
287

288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
    def parse_simulation_header(self):

        try:
            import uproot
            with uproot.open(self.input_url) as input_file:
                run_header_tree = input_file['RunHeaders']
                spectral_index  = run_header_tree['MMcCorsikaRunHeader.fSlopeSpec'].array(library="np")[0]
                view_cone       = run_header_tree['MMcRunHeader.fRandomPointingConeSemiAngle'].array(library="np")[0]
                e_low           = run_header_tree['MMcCorsikaRunHeader.fELowLim'].array(library="np")[0]
                e_high          = run_header_tree['MMcCorsikaRunHeader.fEUppLim'].array(library="np")[0]
                max_impact      = run_header_tree['MMcRunHeader.fImpactMax'].array(library="np")[0]
                n_showers       = np.sum(run_header_tree['MMcRunHeader.fNumSimulatedShowers'].array(library="np"))
                corsika_version = run_header_tree['MMcCorsikaRunHeader.fCorsikaVersion'].array(library="np")[0]
                site_height     = run_header_tree['MMcCorsikaRunHeader.fHeightLev[10]'].array(library="np")[0][0]
                atm_model       = run_header_tree['MMcCorsikaRunHeader.fAtmosphericModel'].array(library="np")[0]
                max_zd          = run_header_tree['MMcRunHeader.fShowerThetaMax'].array(library="np")[0]
                min_zd          = run_header_tree['MMcRunHeader.fShowerThetaMin'].array(library="np")[0]
                max_az          = run_header_tree['MMcRunHeader.fShowerPhiMax'].array(library="np")[0]
                min_az          = run_header_tree['MMcRunHeader.fShowerPhiMin'].array(library="np")[0]
                max_wavelength  = run_header_tree['MMcRunHeader.fCWaveUpper'].array(library="np")[0]
                min_wavelength  = run_header_tree['MMcRunHeader.fCWaveLower'].array(library="np")[0]

                return SimulationConfigContainer(
                    corsika_version=corsika_version,
                    energy_range_min=e_low * u.GeV,
                    energy_range_max=e_high * u.GeV,
                    prod_site_alt=site_height * u.cm,
                    spectral_index=spectral_index,
                    num_showers=n_showers,
                    shower_reuse=1, #not written in the magic root file, but since the sim_events already include shower reuse we artificially set it to 1 (actually every shower reused 5 times for std MAGIC MC)
318
319
                    max_alt=(90. - min_zd) * u.deg,
                    min_alt=(90. - max_zd) * u.deg,
320
321
322
323
324
325
326
327
328
329
330
331
332
333
                    max_az=max_az * u.deg,
                    min_az=min_az * u.deg,
                    max_viewcone_radius=view_cone * u.deg,
                    min_viewcone_radius=0.0 * u.deg,
                    max_scatter_range=max_impact * u.m,
                    min_scatter_range=0.0 * u.m,
                    atmosphere=atm_model,
                    corsika_wlen_min=min_wavelength * u.nm,
                    corsika_wlen_max=max_wavelength * u.nm,
                )

        except ImportError:
            print("Error while importing uproot.")

Moritz Huetten's avatar
Moritz Huetten committed
334
    def _set_active_run(self, run_number):
335
336
337
338
339
340
341
342
343
344
        """
        This internal method sets the run that will be used for data loading.

        Parameters
        ----------
        run_number: int
            The run number to use.

        Returns
        -------
345
        run: MarsRun
Moritz Huetten's avatar
Moritz Huetten committed
346
            The run to use
347
        """
Moritz Huetten's avatar
Moritz Huetten committed
348

349
350
351
        run = dict()
        run['number'] = run_number
        run['read_events'] = 0
Alessio Berti's avatar
Alessio Berti committed
352
        run['data'] = MarsRun(run_file_mask=str(self.input_url))
353
354
355

        return run

356
357
358
359
    @property
    def subarray(self):
        return self._subarray_info

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
360
361
362
363
364
365
366
367
368
    @property
    def is_simulation(self):
        return self.is_mc

    @property
    def datalevels(self):
        return (self.datalevel, )

    @property
369
    def obs_ids(self):
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
370
371
        return self.run_numbers

372
373
374
375
376
377
378
379
380
381
    def _generator(self):
        """
        The default event generator. Return the stereo event
        generator instance.

        Returns
        -------

        """

382
        return self._mono_event_generator(telescope=f"M{self.telescope}")
383
384
385
386
387
388
389
390
391
392
393
394
395
396

    def _stereo_event_generator(self):
        """
        Stereo event generator. Yields DataContainer instances, filled
        with the read event data.

        Returns
        -------

        """

        counter = 0

        # Data container - is initialized once, and data is replaced within it after each yield
397
        data = ArrayEventContainer()
398
399
400
401
402
403

        # Telescopes with data:
        tels_in_file = ["m1", "m2"]
        tels_with_data = {1, 2}

        # Loop over the available data runs
404
        for run_number in self.run_numbers:
405
406
407
408
409
410

            # Removing the previously read data run from memory
            if self.current_run is not None:
                if 'data' in self.current_run:
                    del self.current_run['data']

Moritz Huetten's avatar
Moritz Huetten committed
411
            # Setting the new active run (class MarsRun object)
Moritz Huetten's avatar
Moritz Huetten committed
412
            self.current_run = self._set_active_run(run_number)
413
            
414
415
            # Set monitoring data:
            if not self.is_mc:
Moritz Huetten's avatar
Moritz Huetten committed
416

417
                monitoring_data = self.current_run['data'].monitoring_data
Moritz Huetten's avatar
Moritz Huetten committed
418

419
420
421
                for tel_i, tel_id in enumerate(tels_in_file):
                    monitoring_camera = MonitoringCameraContainer()
                    pedestal_info = PedestalContainer()
422
                    badpixel_info = PixelStatusContainer()
Moritz Huetten's avatar
Moritz Huetten committed
423

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
424
425
426
427
428
429
                    time_tmp = Time(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
                    pedestal_info.sample_time = Time(
                        time_tmp, format='unix', scale='utc', precision=9)
                    # hardcoded number of pedestal events averaged over:
                    pedestal_info.n_events = 500
430
                    pedestal_info.charge_mean = []
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
431
432
433
434
435
436
                    pedestal_info.charge_mean.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
                    pedestal_info.charge_mean.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
                    pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
437
                    pedestal_info.charge_std = []
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
438
439
440
441
442
443
444
445
446
447
448
449
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
                    pedestal_info.charge_std.append(
                        monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

                    t_range = Time(monitoring_data['M{:d}'.format(
                        tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

                    badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                        tel_i + 1)]['badpixelinfo']
Moritz Huetten's avatar
Moritz Huetten committed
450
                    badpixel_info.sample_time_range = t_range
451

452
                    monitoring_camera.pedestal = pedestal_info
453
                    monitoring_camera.pixel_status = badpixel_info
454
455

                    data.mon.tels_with_data = {1, 2}
456
                    data.mon.tel[tel_i + 1] = monitoring_camera
457
            else:
458
459
                assert self.current_run['data'].mcheader_data['M1'] == self.current_run['data'].mcheader_data['M2'], "Simulation configurations are different for M1 and M2 !!!"
                data.mcheader.num_showers = self.current_run['data'].mcheader_data['M1']['sim_nevents']
460
                data.mcheader.shower_reuse = self.current_run['data'].mcheader_data['M1']['sim_reuse']
461
462
463
464
465
466
467
                data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data['M1']['sim_emin']).to(u.TeV) # GeV->TeV
                data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data['M1']['sim_emax']).to(u.TeV) # GeV->TeV
                data.mcheader.spectral_index = self.current_run['data'].mcheader_data['M1']['sim_eslope']
                data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data['M1']['sim_max_impact']).to(u.m) # cm->m
                data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data['M1']['sim_conesemiangle']).to(u.deg)# deg->deg
                if data.mcheader.max_viewcone_radius != 0.:
                    data.mcheader.diffuse = True
468
                else:
469
470
                    data.mcheader.diffuse = False
                    
471

472
473
474
475
476
477
478
479
            # Loop over the events
            for event_i in range(self.current_run['data'].n_stereo_events):
                # Event and run ids
                event_order_number = self.current_run['data'].stereo_ids[event_i][0]
                event_id = self.current_run['data'].event_data['M1']['stereo_event_number'][event_order_number]
                obs_id = self.current_run['number']

                # Reading event data
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
480
481
482
                event_data = self.current_run['data'].get_stereo_event_data(
                    event_i)

483
                data.meta = event_data['mars_meta']
484
485
486

                # Event counter
                data.count = counter
Moritz Huetten's avatar
Moritz Huetten committed
487
488
                data.index.obs_id = obs_id
                data.index.event_id = event_id
489
490
491
492
493
494
495
496
497
498

                # Setting up the R0 container
                data.r0.tel.clear()

                # Setting up the R1 container
                data.r1.tel.clear()

                # Setting up the DL0 container
                data.dl0.tel.clear()

Alessio Berti's avatar
Alessio Berti committed
499
500
501
                # Setting up the DL1 container
                data.dl1.tel.clear()

502
                pointing = PointingContainer()
503
504
505
                # Filling the DL1 container with the event data
                for tel_i, tel_id in enumerate(tels_in_file):
                    # Creating the telescope pointing container
506
507
508
                    pointing_tel = TelescopePointingContainer()
                    
                    pointing_tel.azimuth = np.deg2rad(
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
509
                        event_data['{:s}_pointing_az'.format(tel_id)]) * u.rad
510
511
                    
                    pointing_tel.altitude = np.deg2rad(
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
512
                        90 - event_data['{:s}_pointing_zd'.format(tel_id)]) * u.rad
513
                    
514
                    # pointing.ra = np.deg2rad(
Moritz Huetten's avatar
Moritz Huetten committed
515
                    #    event_data['{:s}_pointing_ra'.format(tel_id)]) * u.rad
516
                    # pointing.dec = np.deg2rad(
Moritz Huetten's avatar
Moritz Huetten committed
517
                    #    event_data['{:s}_pointing_dec'.format(tel_id)]) * u.rad
518

519
                    pointing.tel[tel_i + 1] = pointing_tel
520

Moritz Huetten's avatar
Moritz Huetten committed
521
522
523
524
525
                    # Adding trigger id (MAGIC nomenclature)
                    data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]
                    data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]
                    data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data['M1']['trigger_pattern'][event_order_number]

526
                    # Adding event charge and peak positions per pixel
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
527
528
529
                    data.dl1.tel[tel_i +
                                 1].image = event_data['{:s}_image'.format(tel_id)]
                    data.dl1.tel[tel_i +
Moritz Huetten's avatar
Moritz Huetten committed
530
                                 1].peak_time = event_data['{:s}_pulse_time'.format(tel_id)]
531
532
533
534
535
536
                
                pointing.array_azimuth = np.deg2rad(event_data['m1_pointing_az']) * u.rad
                pointing.array_altitude = np.deg2rad(90 - event_data['m1_pointing_zd']) * u.rad
                pointing.array_ra = np.deg2rad(event_data['m1_pointing_ra']) * u.rad
                pointing.array_dec = np.deg2rad(90 - event_data['m1_pointing_dec']) * u.rad
                data.pointing = pointing
Moritz Huetten's avatar
Moritz Huetten committed
537

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
538
                if not self.is_mc:
Moritz Huetten's avatar
Moritz Huetten committed
539
                    # Adding the event arrival time
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
540
541
                    time_tmp = Time(
                        event_data['mjd'], scale='utc', format='mjd')
Moritz Huetten's avatar
Moritz Huetten committed
542
                    data.trigger.time = Time(
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
543
                        time_tmp, format='unix', scale='utc', precision=9)
Moritz Huetten's avatar
Moritz Huetten committed
544
545
                else:
                    data.mc.energy = event_data['true_energy'] * u.GeV
546
                    data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
547
548
549
550
551
                    # check meaning of 7deg transformation (I.Vovk)
                    data.mc.az = -1 * \
                        (event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
                    data.mc.shower_primary_id = 1 - \
                        event_data['true_shower_primary_id']
Moritz Huetten's avatar
Moritz Huetten committed
552
                    data.mc.h_first_int = event_data['true_h_first_int'] * u.cm
553
554
555
556
557
558
                    
                    # adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
                    # magnetic north is 7 deg westward w.r.t. geographical north
                    rot_corsika = 7 *u.deg
                    data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
                    data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
559
560
561
562
563

                # Setting the telescopes with data
                data.r0.tels_with_data = tels_with_data
                data.r1.tels_with_data = tels_with_data
                data.dl0.tels_with_data = tels_with_data
Moritz Huetten's avatar
Moritz Huetten committed
564
                data.trigger.tels_with_trigger = tels_with_data
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589

                yield data
                counter += 1

        return

    def _mono_event_generator(self, telescope):
        """
        Mono event generator. Yields DataContainer instances, filled
        with the read event data.

        Parameters
        ----------
        telescope: str
            The telescope for which to return events. Can be either "M1" or "M2".

        Returns
        -------

        """

        counter = 0
        telescope = telescope.upper()

        # Data container - is initialized once, and data is replaced within it after each yield
590
        data = ArrayEventContainer()
591
592
593
594
595

        # Telescopes with data:
        tels_in_file = ["M1", "M2"]

        if telescope not in tels_in_file:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
596
597
            raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(
                telescope, tels_in_file))
598
599
600
601

        tel_i = tels_in_file.index(telescope)
        tels_with_data = {tel_i + 1, }

Alessio Berti's avatar
Alessio Berti committed
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
        # Removing the previously read data run from memory
        if self.current_run is not None:
            if 'data' in self.current_run:
                del self.current_run['data']

        # Setting the new active run
        self.current_run = self._set_active_run(self.run_numbers)

        # Set monitoring data:
        if not self.is_mc:

            monitoring_data = self.current_run['data'].monitoring_data

            monitoring_camera = MonitoringCameraContainer()
            pedestal_info = PedestalContainer()
            badpixel_info = PixelStatusContainer()

            time_tmp = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
            pedestal_info.sample_time = Time(
                time_tmp, format='unix', scale='utc', precision=9)
            pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
            pedestal_info.charge_mean = []
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
            pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
            pedestal_info.charge_std = []
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

            t_range = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

            badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfo']
            badpixel_info.sample_time_range = t_range

            monitoring_camera.pedestal = pedestal_info
            monitoring_camera.pixel_status = badpixel_info

            data.mon.tel[tel_i + 1] = monitoring_camera

        if telescope == 'M1':
            n_events = self.current_run['data'].n_mono_events_m1
        else:
            n_events = self.current_run['data'].n_mono_events_m2

        # Loop over the events
        for event_i in range(n_events):
            # Event and run ids
            event_order_number = self.current_run['data'].mono_ids[telescope][event_i]
            event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
            obs_id = self.current_run['number']

            # Reading event data
            event_data = self.current_run['data'].get_mono_event_data(
                event_i, telescope=telescope)

            data.meta = event_data['mars_meta']
Alessio Berti's avatar
Alessio Berti committed
668
669
            data.meta["max_events"] = self.max_events

Alessio Berti's avatar
Alessio Berti committed
670
            data.trigger.event_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
Alessio Berti's avatar
Alessio Berti committed
671
672
673
674
675
676
677
678
679
680
681
682
            data.trigger.tels_with_trigger = tels_with_data
            if self.allowed_tels:
                data.trigger.tels_with_trigger = np.intersect1d(
                    data.trigger.tels_with_trigger,
                    self.subarray.tel_ids,
                    assume_unique=True,)
            if not self.is_mc:
                # Adding the event arrival time
                time_tmp = Time(
                    event_data['mjd'], scale='utc', format='mjd')
                data.trigger.tel[tel_i + 1] = TelescopeTriggerContainer(time=Time(
                    time_tmp, format='unix', scale='utc', precision=9))
Alessio Berti's avatar
Alessio Berti committed
683
684
685
686
687
688
689
690
691
692
693

            # Event counter
            data.count = counter
            data.index.obs_id = obs_id
            data.index.event_id = event_id

            # Setting up the R0 container
            data.r0.tel.clear()
            data.r1.tel.clear()
            data.dl0.tel.clear()
            data.dl1.tel.clear()
Alessio Berti's avatar
Alessio Berti committed
694
            data.pointing.tel.clear()
Alessio Berti's avatar
Alessio Berti committed
695
696
697

            # Creating the telescope pointing container
            pointing = PointingContainer()
Alessio Berti's avatar
Alessio Berti committed
698
699
700
701
702
703
            pointing_tel = TelescopePointingContainer(
                azimuth = np.deg2rad(
                event_data['pointing_az']) * u.rad,
                altitude = np.deg2rad(
                90 - event_data['pointing_zd']) * u.rad,)

Alessio Berti's avatar
Alessio Berti committed
704
            pointing.tel[tel_i + 1] = pointing_tel
Alessio Berti's avatar
Alessio Berti committed
705
706

            pointing.array_azimuth  = np.deg2rad(event_data['pointing_az']) * u.rad
Alessio Berti's avatar
Alessio Berti committed
707
            pointing.array_altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
Alessio Berti's avatar
Alessio Berti committed
708
709
710
            pointing.array_ra       = np.deg2rad(event_data['pointing_ra']) * u.rad
            pointing.array_dec      = np.deg2rad(90 - event_data['pointing_dec']) * u.rad

Alessio Berti's avatar
Alessio Berti committed
711
            data.pointing = pointing
Moritz Huetten's avatar
Moritz Huetten committed
712

Alessio Berti's avatar
Alessio Berti committed
713
714
715
            # Adding event charge and peak positions per pixel
            data.dl1.tel[tel_i + 1].image = event_data['image']
            data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']
Moritz Huetten's avatar
Moritz Huetten committed
716

Alessio Berti's avatar
Alessio Berti committed
717
718
            if self.is_mc:
                rot_corsika = 7 *u.deg
Alessio Berti's avatar
Alessio Berti committed
719
720
721
                # check meaning of 7deg transformation (I.Vovk)
                # adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
                # magnetic north is 7 deg westward w.r.t. geographical north
Alessio Berti's avatar
Alessio Berti committed
722
723
724
                data.simulation.shower = SimulatedShowerContainer(
                    energy = u.Quantity(event_data['true_energy'], u.GeV),
                    alt = Angle((np.pi/2 - event_data['true_zd']), u.rad),
Alessio Berti's avatar
Alessio Berti committed
725
                    az = Angle(-1 * (event_data['true_az'] - np.deg2rad(180 - rot_corsika.value)), u.rad),
Alessio Berti's avatar
Alessio Berti committed
726
727
728
729
730
                    shower_primary_id = 1 - event_data['true_shower_primary_id'],
                    h_first_int = u.Quantity(event_data['true_h_first_int'], u.cm),
                    core_x = u.Quantity((event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika)), u.cm),
                    core_y = u.Quantity((event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika)), u.cm),
                )
Alessio Berti's avatar
Alessio Berti committed
731
732
733

            yield data
            counter += 1
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755

        return

    def _pedestal_event_generator(self, telescope):
        """
        Pedestal event generator. Yields DataContainer instances, filled
        with the read event data.

        Parameters
        ----------
        telescope: str
            The telescope for which to return events. Can be either "M1" or "M2".

        Returns
        -------

        """

        counter = 0
        telescope = telescope.upper()

        # Data container - is initialized once, and data is replaced within it after each yield
Moritz Huetten's avatar
Moritz Huetten committed
756
        data = EventAndMonDataContainer()
757
758
759
760
761

        # Telescopes with data:
        tels_in_file = ["M1", "M2"]

        if telescope not in tels_in_file:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
762
763
            raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(
                telescope, tels_in_file))
764
765
766
767
768

        tel_i = tels_in_file.index(telescope)
        tels_with_data = {tel_i + 1, }

        # Loop over the available data runs
769
        for run_number in self.run_numbers:
770
771
772
773
774
775
776

            # Removing the previously read data run from memory
            if self.current_run is not None:
                if 'data' in self.current_run:
                    del self.current_run['data']

            # Setting the new active run
Moritz Huetten's avatar
Moritz Huetten committed
777
            self.current_run = self._set_active_run(run_number)
778

Moritz Huetten's avatar
Moritz Huetten committed
779
            monitoring_data = self.current_run['data'].monitoring_data
Moritz Huetten's avatar
Moritz Huetten committed
780

Moritz Huetten's avatar
Moritz Huetten committed
781
782
783
            monitoring_camera = MonitoringCameraContainer()
            pedestal_info = PedestalContainer()
            badpixel_info = PixelStatusContainer()
Moritz Huetten's avatar
Moritz Huetten committed
784

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
785
786
787
788
            time_tmp = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
            pedestal_info.sample_time = Time(
                time_tmp, format='unix', scale='utc', precision=9)
789
            pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
Moritz Huetten's avatar
Moritz Huetten committed
790
            pedestal_info.charge_mean = []
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
791
792
793
794
795
796
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
            pedestal_info.charge_mean.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
            pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
                tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
Moritz Huetten's avatar
Moritz Huetten committed
797
            pedestal_info.charge_std = []
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
798
799
800
801
802
803
804
805
806
807
808
809
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
            pedestal_info.charge_std.append(
                monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])

            t_range = Time(monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')

            badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
                tel_i + 1)]['badpixelinfo']
Moritz Huetten's avatar
Moritz Huetten committed
810
            badpixel_info.sample_time_range = t_range
Moritz Huetten's avatar
Moritz Huetten committed
811

Moritz Huetten's avatar
Moritz Huetten committed
812
813
814
815
816
817
            monitoring_camera.pedestal = pedestal_info
            monitoring_camera.pixel_status = badpixel_info

            data.mon.tels_with_data = tels_with_data
            data.mon.tel[tel_i + 1] = monitoring_camera

818
819
820
821
822
823
824
825
826
827
828
829
830
            if telescope == 'M1':
                n_events = self.current_run['data'].n_pedestal_events_m1
            else:
                n_events = self.current_run['data'].n_pedestal_events_m2

            # Loop over the events
            for event_i in range(n_events):
                # Event and run ids
                event_order_number = self.current_run['data'].pedestal_ids[telescope][event_i]
                event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
                obs_id = self.current_run['number']

                # Reading event data
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
831
832
833
                event_data = self.current_run['data'].get_pedestal_event_data(
                    event_i, telescope=telescope)

834
                data.meta = event_data['mars_meta']
835
836
837

                # Event counter
                data.count = counter
Moritz Huetten's avatar
Moritz Huetten committed
838
839
                data.index.obs_id = obs_id
                data.index.event_id = event_id
840

841
842
                # Setting up the R0 container
                data.r0.tel.clear()
843
                data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
844
845
846

                # Setting up the R1 container
                data.r1.tel.clear()
847
                data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
848
849
850

                # Setting up the DL0 container
                data.dl0.tel.clear()
851
                data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
852
853
854

                # Creating the telescope pointing container
                pointing = TelescopePointingContainer()
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
855
856
857
858
                pointing.azimuth = np.deg2rad(
                    event_data['pointing_az']) * u.rad
                pointing.altitude = np.deg2rad(
                    90 - event_data['pointing_zd']) * u.rad
Moritz Huetten's avatar
Moritz Huetten committed
859
860
                #pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
                #pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad
861
862

                # Adding the pointing container to the event data
Moritz Huetten's avatar
Moritz Huetten committed
863
                data.pointing.tel[tel_i + 1] = pointing
864
865
866

                # Adding event charge and peak positions per pixel
                data.dl1.tel[tel_i + 1].image = event_data['image']
Moritz Huetten's avatar
Moritz Huetten committed
867
                data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']
868
869
870

                # Adding the event arrival time
                time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
Moritz Huetten's avatar
Moritz Huetten committed
871
                data.trigger.time = Time(
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
872
                    time_tmp, format='unix', scale='utc', precision=9)
873
874
875
876
877

                # Setting the telescopes with data
                data.r0.tels_with_data = tels_with_data
                data.r1.tels_with_data = tels_with_data
                data.dl0.tels_with_data = tels_with_data
Moritz Huetten's avatar
Moritz Huetten committed
878
                data.trigger.tels_with_trigger = tels_with_data
879
880
881
882
883
884
885

                yield data
                counter += 1

        return


Moritz Huetten's avatar
Moritz Huetten committed
886
class MarsRun:
887
888
889
890
    """
    This class implements reading of the event data from a single MAGIC data run.
    """

Moritz Huetten's avatar
Moritz Huetten committed
891
    def __init__(self, run_file_mask, filter_list=None):
892
893
894
895
896
897
898
899
900
        """
        Constructor of the class. Defines the run to use and the camera pixel arrangement.

        Parameters
        ----------
        run_file_mask: str
            A path mask for files belonging to the run. Must correspond to a single run
            or an exception will be raised. Must correspond to calibrated ("sorcerer"-level)
            data.
901
902
903
        filter_list: list, optional
            A list of files, to which the run_file_mask should be applied. If None, all the
            files satisfying run_file_mask will be used. Defaults to None.
904
905
        """

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
906
        self.n_camera_pixels = GEOM.n_pixels
907

908
909
910
911
        self.run_file_mask = run_file_mask

        # Preparing the lists of M1/2 data files
        file_list = glob.glob(run_file_mask)
912
913
914
915
916

        # Filtering out extra files if necessary
        if filter_list is not None:
            file_list = list(set(file_list) & set(filter_list))

Moritz Huetten's avatar
astyled    
Moritz Huetten committed
917
918
919
920
        self.m1_file_list = list(
            filter(lambda name: '_M1_' in name, file_list))
        self.m2_file_list = list(
            filter(lambda name: '_M2_' in name, file_list))
921
        self.m1_file_list.sort()
922
923
924
        self.m2_file_list.sort()

        # Retrieving the list of run numbers corresponding to the data files
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
925
        run_info = list(
926
            map(MAGICEventSource.get_run_info_from_name, file_list))
Moritz Huetten's avatar
Moritz Huetten committed
927
        run_numbers = [i[0] for i in run_info]
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
928
        is_mc_runs = [i[1] for i in run_info]
929
930
931

        run_numbers = np.unique(run_numbers)
        is_mc_runs = np.unique(is_mc_runs)
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
932
        # Checking if run type (data/MC) is consistent:
933
        if len(is_mc_runs) > 1:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
934
935
936
            raise ValueError(
                "Run type is not consistently data or MC: {}".format(is_mc_runs))

937
        self.is_mc = is_mc_runs[0]
938
939
940

        # Checking if a single run is going to be read
        if len(run_numbers) > 1:
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
941
942
            raise ValueError(
                "Run mask corresponds to more than one run: {}".format(run_numbers))
943

944
        # Reading the data
Moritz Huetten's avatar
astyled    
Moritz Huetten committed
945
946
947
948
949
        m1_data = self.load_events(
            self.m1_file_list, self.is_mc, self.n_camera_pixels)
        m2_data = self.load_events(
            self.m2_file_list, self.is_mc, self.n_camera_pixels)

950
        # Getting the event data
951
        self.event_data = dict()
952
953
954
955
956
957
958
        self.event_data['M1'] = m1_data[0]
        self.event_data['M2'] = m2_data[0]

        # Getting the monitoring data
        self.monitoring_data = dict()
        self.monitoring_data['M1'] = m1_data[1]
        self.monitoring_data['M2'] = m2_data[1]
959

960
        # Getting the run-wise MC header data
961
        if self.is_mc:
962
963
964
            self.mcheader_data = dict()
            self.mcheader_data['M1'] = m1_data[2]
            self.mcheader_data['M2'] = m2_data[2]
965

966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
        # Detecting pedestal events
        self.pedestal_ids = self._find_pedestal_events()
        # Detecting stereo events
        self.stereo_ids = self._find_stereo_events()
        # Detecting mono events
        self.mono_ids = self._find_mono_events()

    @property
    def n_events_m1(self):
        return len(self.event_data['M1']['MJD'])

    @property
    def n_events_m2(self):
        return len(self.event_data['M2']['MJD'])

    @property
    def n_stereo_events(self):
        return len(self.stereo_ids)

    @property
    def n_mono_events_m1(self):
        return len(self.mono_ids['M1'])

    @property
    def n_mono_events_m2(self):
        return len(self.mono_ids['M2'])

    @property
    def n_pedestal_events_m1(self):
        return len(self.pedestal_ids['M1'])

    @property
    def n_pedestal_events_m2(self):
        return len(self.pedestal_ids['M2'])