diff --git a/pyproject.toml b/pyproject.toml index b5440dc..fd2ded1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "srcsim" -version = "1.0.0" +version = "1.1.0" authors = [ { name="Ievgen Vovk", email="vovk@icrr.u-tokyo.ac.jp" }, { name="Marcel Strzys", email="strzys@icrr.u-tokyo.ac.jp" }, @@ -20,12 +20,12 @@ classifiers = [ "Operating System :: OS Independent", ] dependencies = [ - "astropy==5.1", + "astropy>=5.1", "numpy>=1.21.0", - "pandas==1.4.1", + "pandas>=1.4.1", "scipy>=1.5.0", - "tables==3.7.0", - "PyYAML==5.3.1" + "tables>=3.7.0", + "PyYAML>=5.3.1" ] [project.urls] @@ -35,3 +35,4 @@ dependencies = [ [project.scripts] getruns = "srcsim.scripts.get_runs:main" simrun = "srcsim.scripts.sim_run:main" +magic_to_dl2 = "srcsim.scripts.magic_to_dl2:main" diff --git a/src/srcsim/magic/mc/__init__.py b/src/srcsim/magic/mc/__init__.py new file mode 100644 index 0000000..147d797 --- /dev/null +++ b/src/srcsim/magic/mc/__init__.py @@ -0,0 +1,3 @@ +from .info import MagicMcInfo +from .events import MagicStereoEvents, MagicMcOrigEvents +from .file import MagicMcFile \ No newline at end of file diff --git a/src/srcsim/magic/mc/events.py b/src/srcsim/magic/mc/events.py new file mode 100644 index 0000000..1ca03b8 --- /dev/null +++ b/src/srcsim/magic/mc/events.py @@ -0,0 +1,384 @@ +import numpy +import uproot +import pandas +import astropy.units as u +from dataclasses import dataclass + + +@dataclass(frozen=True) +class EventSample: + tel_id: u.Quantity + obs_id: u.Quantity + event_id: u.Quantity + reco_src_x: u.Quantity + reco_src_y: u.Quantity + reco_alt: u.Quantity + reco_az: u.Quantity + reco_ra: u.Quantity + reco_dec: u.Quantity + reco_energy: u.Quantity + src_x: u.Quantity + src_y: u.Quantity + ra_tel: u.Quantity + dec_tel: u.Quantity + az_tel: u.Quantity + alt_tel: u.Quantity + delta_t: u.Quantity + gammaness: u.Quantity + mjd: u.Quantity = u.Quantity([], unit='d') + mc_alt: u.Quantity = u.Quantity([], unit='deg') + mc_az: u.Quantity = u.Quantity([], unit='deg') + mc_energy: u.Quantity = u.Quantity([], unit='TeV') + mc_alt_tel: u.Quantity = u.Quantity([], unit='deg') + mc_az_tel: u.Quantity = u.Quantity([], unit='deg') + + +@dataclass(frozen=True) +class McOrigEventSample: + az_tel: u.Quantity + alt_tel: u.Quantity + energy: u.Quantity + src_x: u.Quantity = u.Quantity([], unit='m') + src_y: u.Quantity = u.Quantity([], unit='m') + alt: u.Quantity = u.Quantity([], unit='deg') + az: u.Quantity = u.Quantity([], unit='deg') + + +@dataclass(frozen=True) +class MagicStereoEvents(EventSample): + file_name: str = '' + + def __str__(self) -> str: + summary = \ +f"""{type(self).__name__} instance + {'File name':.<20s}: {self.file_name} + {'MC':.<20s}: {self.is_mc} + {'Events':.<20s}: {self.n_events} +""" + return summary + + def __repr__(self): + print(str(self)) + + return super().__repr__() + + @property + def n_events(self): + return len(self.reco_alt) + + @property + def is_mc(self): + return len(self.mc_alt) > 0 + + @classmethod + def from_file(cls, file_name): + cam2angle = 1.5 / 450 *u.Unit('deg/mm') + north_offset = 7 + magic_tel_id = 5 * u.one # tel_id expected from the MAGIC+LST simulations + + event_data = dict() + + data_array_list = [ + 'MTriggerPattern_1.fPrescaled', + 'MRawEvtHeader_1.fStereoEvtNumber', + 'MRawEvtHeader_1.fDAQEvtNumber', + 'MRawEvtHeader_1.fTimeDiff', + 'MPointingPos_1.fZd', + 'MPointingPos_1.fAz', + 'MPointingPos_1.fRa', + 'MPointingPos_1.fDec', + 'MEnergyEst.fEnergy', + 'MHadronness.fHadronness', + 'MRawEvtHeader_1.fDAQEvtNumber', + 'MSrcPosCam_1.fX', + 'MSrcPosCam_1.fY', + 'MStereoParDisp.fDirectionX', + 'MStereoParDisp.fDirectionY', + 'MStereoParDisp.fDirectionRA', + 'MStereoParDisp.fDirectionDec', + 'MStereoParDisp.fDirectionAz', + 'MStereoParDisp.fDirectionZd', + ] + + time_array_list = ['MTime_1.fMjd', 'MTime_1.fTime.fMilliSec', 'MTime_1.fNanoSec'] + + mc_array_list = [ + 'MMcEvt_1.fEnergy', + 'MMcEvt_1.fTheta', + 'MMcEvt_1.fPhi', + ] + + names_mapping = { + 'MTriggerPattern_1.fPrescaled': 'trigger_pattern', + 'MRawEvtHeader_1.fStereoEvtNumber': 'stereo_event_number', + 'MRawEvtHeader_1.fDAQEvtNumber': 'daq_event_number', + 'MRawEvtHeader_1.fTimeDiff': 'delta_t', + 'MPointingPos_1.fZd': 'zd_tel', + 'MPointingPos_1.fAz': 'az_tel', + 'MPointingPos_1.fRa': 'ra_tel', + 'MPointingPos_1.fDec': 'dec_tel', + 'MSrcPosCam_1.fX': 'src_x', + 'MSrcPosCam_1.fY': 'src_y', + 'MStereoParDisp.fDirectionX': 'reco_src_x', + 'MStereoParDisp.fDirectionY': 'reco_src_y', + 'MStereoParDisp.fDirectionRA': 'reco_ra', + 'MStereoParDisp.fDirectionDec': 'reco_dec', + 'MStereoParDisp.fDirectionAz': 'reco_az', + 'MStereoParDisp.fDirectionZd': 'reco_zd', + 'MHadronness.fHadronness': 'hadronness', + 'MEnergyEst.fEnergy': 'reco_energy', + 'MMcEvt_1.fEnergy': 'mc_energy', + 'MMcEvt_1.fTheta': 'mc_zd', + 'MMcEvt_1.fPhi': 'mc_az' + } + + data_units = { + 'run_number': u.one, + 'daq_event_number': u.one, + 'delta_t': u.s, + 'src_x': u.mm, + 'src_y': u.mm, + 'reco_src_x': u.deg, + 'reco_src_y': u.deg, + 'reco_alt': u.deg, + 'reco_az': u.deg, + 'reco_zd': u.deg, + 'reco_ra': u.deg, + 'reco_dec': u.deg, + 'reco_energy': u.GeV, + 'gammaness': u.one, + 'hadronness': u.one, + 'mjd': u.d, + 'alt_tel': u.deg, + 'az_tel': u.deg, + 'zd_tel': u.deg, + 'ra_tel': u.deg, + 'dec_tel':u.deg, + 'mc_alt':u.deg, + 'mc_az':u.deg, + 'mc_zd':u.deg, + 'mc_energy': u.GeV, + } + + with uproot.open(file_name) as input_file: + if 'Events' in input_file: + data = input_file['Events'].arrays(data_array_list) + + # Mapping the read structure to the alternative names + for key in data.fields: + name = names_mapping[key] + event_data[name] = data[key] + + is_mc = 'MMcEvt_1.' in input_file['Events'] + if is_mc: + data = input_file['Events'].arrays(mc_array_list) + + # Mapping the read structure to the alternative names + for key in data.fields: + name = names_mapping[key] + event_data[name] = data[key] + + # Post processing + event_data['mc_zd'] = numpy.degrees(event_data['mc_zd']) + event_data['mc_az'] = numpy.degrees(event_data['mc_az']) + # Transformation from Monte Carlo to usual azimuth + event_data['mc_az'] = -1 * (event_data['mc_az'] - 180 + north_offset) + else: + # Reading the event arrival time information + data = input_file['Events'].arrays(time_array_list) + + # Computing the event arrival time + mjd = data['MTime_1.fMjd'] + millisec = data['MTime_1.fTime.fMilliSec'] + nanosec = data['MTime_1.fNanoSec'] + + event_data['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400.0 + + run_number = input_file['RunHeaders']['MRawRunHeader_1./MRawRunHeader_1.fRunNumber'].array()[0] + + else: + # The file is likely corrupted, so return empty arrays + for key in names_mapping: + name = names_mapping[key] + event_data[name] = numpy.zeros(0) + run_number = 0 + + event_data['mc_alt'] = 90 - event_data['mc_zd'] + event_data['alt_tel'] = 90 - event_data['zd_tel'] + event_data['reco_alt'] = 90 - event_data['reco_zd'] + event_data['gammaness'] = 1 - event_data['hadronness'] + + for key in event_data: + event_data[key] = event_data[key].to_numpy() + + event_data['run_number'] = numpy.repeat(run_number, len(event_data['az_tel'])) + + for key in event_data: + if key in data_units: + event_data[key] = event_data[key] * data_units[key] + + event_data['reco_src_x'] /= cam2angle + event_data['reco_src_y'] /= cam2angle + + is_mc = 'mc_energy' in event_data + + if is_mc: + events = MagicStereoEvents( + tel_id = numpy.repeat(magic_tel_id, len(event_data['run_number'])), + obs_id = event_data['run_number'], + event_id = event_data['daq_event_number'], + reco_src_x = event_data['reco_src_x'], + reco_src_y = event_data['reco_src_y'], + reco_alt = event_data['reco_alt'], + reco_az = event_data['reco_az'], + reco_ra = event_data['reco_ra'], + reco_dec = event_data['reco_dec'], + reco_energy = event_data['reco_energy'], + src_x = event_data['src_x'], + src_y = event_data['src_y'], + ra_tel = event_data['ra_tel'], + dec_tel = event_data['dec_tel'], + az_tel = event_data['az_tel'], + alt_tel = event_data['alt_tel'], + delta_t = event_data['delta_t'], + gammaness = event_data['gammaness'], + mc_alt = event_data['mc_alt'], + mc_az = event_data['mc_az'], + mc_energy = event_data['mc_energy'], + mc_alt_tel = event_data['alt_tel'], + mc_az_tel = event_data['az_tel'], + file_name = file_name + ) + else: + events = MagicStereoEvents( + tel_id = numpy.repeat(magic_tel_id, len(event_data['run_number'])), + obs_id = event_data['run_number'], + event_id = event_data['daq_event_number'], + reco_src_x = event_data['reco_src_x'], + reco_src_y = event_data['reco_src_y'], + reco_alt = event_data['reco_alt'], + reco_az = event_data['reco_az'], + reco_ra = event_data['reco_ra'], + reco_dec = event_data['reco_dec'], + reco_energy = event_data['reco_energy'], + src_x = event_data['src_x'], + src_y = event_data['src_y'], + ra_tel = event_data['ra_tel'], + dec_tel = event_data['dec_tel'], + az_tel = event_data['az_tel'], + alt_tel = event_data['alt_tel'], + delta_t = event_data['delta_t'], + mjd = event_data['mjd'], + gammaness = event_data['gammaness'], + file_name = file_name + ) + + return events + + def to_df(self): + units = dict( + # delta_t = u.s, + tel_id = u.one, + obs_id = u.one, + event_id = u.one, + src_x = u.m, + src_y = u.m, + reco_src_x = u.m, + reco_src_y = u.m, + reco_alt = u.rad, + reco_az = u.rad, + reco_ra = u.rad, + reco_dec = u.rad, + reco_energy = u.TeV, + gammaness = u.one, + mjd = u.d, + alt_tel = u.rad, + az_tel = u.rad, + ra_tel = u.rad, + dec_tel =u.rad, + mc_alt =u.rad, + mc_az = u.rad, + mc_energy = u.TeV, + mc_alt_tel = u.rad, + mc_az_tel = u.rad, + ) + data = { + key: self.__getattribute__(key).to(units[key]).value + for key in units + if len(self.__getattribute__(key)) + } + return pandas.DataFrame(data=data) + + +@dataclass(frozen=True) +class MagicMcOrigEvents(McOrigEventSample): + file_name: str = '' + + def __str__(self) -> str: + summary = \ +f"""{type(self).__name__} instance + {'File name':.<20s}: {self.file_name} + {'Events':.<20s}: {self.n_events} +""" + return summary + + def __repr__(self): + print(str(self)) + + return super().__repr__() + + @property + def n_events(self): + return len(self.az_tel) + + @classmethod + def from_file(cls, file_name): + north_offset = 7 * u.deg + + data_array_list = [ + 'MMcEvtBasic_1.fEnergy', + 'MMcEvtBasic_1.fTelescopePhi', + 'MMcEvtBasic_1.fTelescopeTheta', + 'MSrcPosCam_1.fX', + 'MSrcPosCam_1.fY', + ] + + data_units = { + 'MMcEvtBasic_1.fEnergy': u.GeV, + 'MSrcPosCam_1.fX': u.mm, + 'MSrcPosCam_1.fY': u.mm, + 'MMcEvtBasic_1.fTelescopePhi': u.rad, + 'MMcEvtBasic_1.fTelescopeTheta': u.rad, + } + + with uproot.open(file_name) as input_file: + if 'OriginalMC' in input_file: + data = input_file['OriginalMC'].arrays(data_array_list) + else: + # The file is likely corrupted, so return empty arrays + data = { + key: numpy.zeros(0) + for key in data_array_list + } + + event_data = { + key: data[key].to_numpy() * data_units[key] + for key in data.fields + } + + event_data['MMcEvtBasic_1.fTelescopeAlt'] = numpy.pi/2 * u.rad - event_data['MMcEvtBasic_1.fTelescopeTheta'] + # Transformation from Monte Carlo to usual azimuth + data['MMcEvtBasic_1.fTelescopePhi'] = -1 * ( + data['MMcEvtBasic_1.fTelescopePhi'] - (180 * u.deg + north_offset).to(data_units['MMcEvtBasic_1.fTelescopePhi']) + ) + + events = MagicMcOrigEvents( + az_tel = event_data['MMcEvtBasic_1.fTelescopePhi'], + alt_tel = event_data['MMcEvtBasic_1.fTelescopeAlt'], + energy = event_data['MMcEvtBasic_1.fEnergy'], + src_x = event_data['MSrcPosCam_1.fX'], + src_y = event_data['MSrcPosCam_1.fY'], + file_name = file_name + ) + + return events \ No newline at end of file diff --git a/src/srcsim/magic/mc/file.py b/src/srcsim/magic/mc/file.py new file mode 100644 index 0000000..b28c19e --- /dev/null +++ b/src/srcsim/magic/mc/file.py @@ -0,0 +1,176 @@ +import numpy +import healpy +import astropy.units as u +from dataclasses import dataclass + +from tables import open_file +from tables import IsDescription +from tables import UInt32Col, UInt64Col, Float32Col + +from .info import MagicMcInfo +from .events import MagicStereoEvents, MagicMcOrigEvents + + +class TablesMcInfo(IsDescription): + obs_id = UInt64Col() + num_showers = UInt64Col() + shower_reuse = UInt32Col() + energy_range_min = Float32Col() + energy_range_max = Float32Col() + spectral_index = Float32Col() + max_scatter_range = Float32Col() + min_scatter_range = Float32Col() + max_viewcone_radius = Float32Col() + min_viewcone_radius = Float32Col() + + +@dataclass(frozen=True) +class MagicMcFile: + file_name: str + meta: MagicMcInfo + events_triggered: MagicStereoEvents + events_simulated: MagicMcOrigEvents + + def __str__(self) -> str: + summary = \ +f"""{type(self).__name__} instance + {'File name':.<20s}: {self.file_name} + {'Simulated events':.<20s}: {self.events_simulated.n_events} + {'Triggered events':.<20s}: {self.events_triggered.n_events} + {'Energy range':.<20s}: {self.meta.energy_range_min} - {self.meta.energy_range_max} + {'Spectral index':.<20s}: {self.meta.spectral_index} + {'Scatter range':.<20s}: {self.meta.min_scatter_range} - {self.meta.max_scatter_range} + {'Viewcone':.<20s}: {self.meta.min_viewcone_radius} - {self.meta.max_viewcone_radius} + {'Telescope azimuth':.<20s}: {self.events_simulated.az_tel.to('deg').min()} - {self.events_simulated.az_tel.to('deg').max()} + {'Telescope altitude':.<20s}: {self.events_simulated.alt_tel.to('deg').min()} - {self.events_simulated.alt_tel.to('deg').max()} +""" + return summary + + @classmethod + def from_file(cls, file_name): + self = MagicMcFile( + file_name = file_name, + meta = MagicMcInfo.from_file(file_name), + events_triggered = MagicStereoEvents.from_file(file_name), + events_simulated = MagicMcOrigEvents.from_file(file_name) + ) + return self + + def write(self, file_name): + units = dict( + energy = u.TeV, + angle = u.rad, + distance = u.m, + viewcone = u.deg + ) + + # Special treatment for run config as it's structure may be too complex for Pandas + with open_file(file_name, mode="w", title="MAGIC MC file") as output: + group = output.create_group("/", 'simulation', 'Simulation information') + table = output.create_table(group, 'run_config', TablesMcInfo, 'Simulation run configuration') + + entry = table.row + entry['obs_id'] = self.meta.obs_id.value + entry['num_showers'] = self.meta.num_showers.value + entry['shower_reuse'] = 1 + entry['spectral_index'] = self.meta.spectral_index.value + entry['energy_range_min'] = self.meta.energy_range_min.to(units['energy']).value + entry['energy_range_max'] = self.meta.energy_range_max.to(units['energy']).value + entry['max_scatter_range'] = self.meta.max_scatter_range.to(units['distance']).value + entry['min_scatter_range'] = self.meta.min_scatter_range.to(units['distance']).value + entry['max_viewcone_radius'] = self.meta.max_viewcone_radius.to(units['viewcone']).value + entry['min_viewcone_radius'] = self.meta.min_viewcone_radius.to(units['viewcone']).value + entry.append() + table.flush() + + df = self.events_triggered.to_df() + # Overwriting obs_id as it is not properly stored for events in the MC files + df['obs_id'] = self.meta.obs_id.value + # TODO: check if one can use "MAGIC_MAGICCam" instead + df.to_hdf( + file_name, + key='/dl2/event/telescope/parameters/LST_LSTCam', + mode='a' + ) + + def healpy_split(self, nside): + npix = healpy.nside2npix(nside) + + pix_ids_triggered = healpy.ang2pix( + nside, + numpy.pi/2 - self.events_triggered.alt_tel.to('rad').value, + self.events_triggered.az_tel.to('rad').value, + ) + + pix_ids_simulated = healpy.ang2pix( + nside, + numpy.pi/2 - self.events_simulated.alt_tel.to('rad').value, + self.events_simulated.az_tel.to('rad').value, + ) + + files = [] + + for pix_id in range(npix): + n_simulated = u.one * numpy.sum(pix_ids_simulated == pix_id) + meta = MagicMcInfo( + obs_id = self.meta.obs_id, + num_showers = n_simulated, + energy_range_min = self.meta.energy_range_min, + energy_range_max = self.meta.energy_range_max, + spectral_index = self.meta.spectral_index, + max_scatter_range = self.meta.max_scatter_range, + min_scatter_range = self.meta.min_scatter_range, + max_viewcone_radius = self.meta.max_viewcone_radius, + min_viewcone_radius = self.meta.min_viewcone_radius + ) + + data = self.events_simulated + selection = pix_ids_simulated == pix_id + events_simulated = MagicMcOrigEvents( + az_tel = data.az_tel[selection], + alt_tel = data.alt_tel[selection], + energy = data.energy[selection], + src_x = data.src_x[selection], + src_y = data.src_y[selection], + file_name = data.file_name + ) + + data = self.events_triggered + selection = pix_ids_triggered == pix_id + events_triggered = MagicStereoEvents( + obs_id = data.obs_id[selection], + event_id = data.event_id[selection], + reco_src_x = data.reco_src_x[selection], + reco_src_y = data.reco_src_y[selection], + reco_alt = data.reco_alt[selection], + reco_az = data.reco_az[selection], + reco_ra = data.reco_ra[selection], + reco_dec = data.reco_dec[selection], + reco_energy = data.reco_energy[selection], + src_x = data.src_x[selection], + src_y = data.src_y[selection], + ra_tel = data.ra_tel[selection], + dec_tel = data.dec_tel[selection], + az_tel = data.az_tel[selection], + alt_tel = data.alt_tel[selection], + delta_t = data.delta_t[selection], + gammaness = data.gammaness[selection], + mjd = data.mjd[selection] if len(data.mjd) else data.mjd, + mc_alt = data.mc_alt[selection] if len(data.mc_alt) else data.mc_alt, + mc_az = data.mc_az[selection] if len(data.mc_az) else data.mc_az, + mc_energy = data.mc_energy[selection] if len(data.mc_energy) else data.mc_energy, + mc_alt_tel = data.mc_alt_tel[selection] if len(data.mc_alt_tel) else data.mc_alt_tel, + mc_az_tel = data.mc_az_tel[selection] if len(data.mc_az_tel) else data.mc_az_tel, + file_name = self.events_triggered.file_name + ) + + files.append( + MagicMcFile( + file_name = self.file_name, + meta = meta, + events_triggered = events_triggered, + events_simulated = events_simulated, + ) + ) + + return files \ No newline at end of file diff --git a/src/srcsim/magic/mc/info.py b/src/srcsim/magic/mc/info.py new file mode 100644 index 0000000..34b0b93 --- /dev/null +++ b/src/srcsim/magic/mc/info.py @@ -0,0 +1,90 @@ +import uproot +import pandas +import astropy.units as u +from dataclasses import dataclass + + +@dataclass(frozen=True) +class McInfo: + obs_id: u.Quantity + num_showers: u.Quantity + energy_range_min: u.Quantity + energy_range_max: u.Quantity + spectral_index: u.Quantity + max_scatter_range: u.Quantity + min_scatter_range: u.Quantity + max_viewcone_radius: u.Quantity + min_viewcone_radius: u.Quantity + + +@dataclass(frozen=True) +class MagicMcInfo(McInfo): + file_name: str = '' + + def __repr__(self): + print( +f"""{type(self).__name__} instance + {'File name':.<20s}: {self.file_name} + {'Events':.<20s}: {self.num_showers} + {'Energy min':.<20s}: {self.energy_range_min} + {'Energy max':.<20s}: {self.energy_range_max} + {'Spectral index':.<20s}: {self.spectral_index} + {'Scatter range min':.<20s}: {self.min_scatter_range} + {'Scatter range max':.<20s}: {self.max_scatter_range} + {'Viewcone min':.<20s}: {self.min_viewcone_radius} + {'Viewcone max':.<20s}: {self.max_viewcone_radius} +""" + ) + + return super().__repr__() + + @classmethod + def from_file(cls, file_name): + meta = dict() + with uproot.open(file_name) as input_file: + meta['obs_id'] = u.one * int( + input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fMcRunNumber'].array()[0] + ) + meta['num_showers'] = u.one * int( + input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fNumSimulatedShowers'].array()[0] + ) + meta['energy_range_min'] = u.GeV * input_file['RunHeaders']['MMcCorsikaRunHeader./MMcCorsikaRunHeader.fELowLim'].array()[0] + meta['energy_range_max'] = u.GeV * input_file['RunHeaders']['MMcCorsikaRunHeader./MMcCorsikaRunHeader.fEUppLim'].array()[0] + meta['spectral_index'] = u.one * input_file['RunHeaders']['MMcRunHeader_1.fSlopeSpec'].array()[0] + meta['max_scatter_range'] = u.cm * input_file['RunHeaders']['MMcRunHeader_1.fImpactMax'].array()[0] + meta['min_scatter_range'] = u.cm * 0 + meta['max_viewcone_radius'] = u.deg * input_file['RunHeaders']['MMcRunHeader_1./MMcRunHeader_1.fRandomPointingConeSemiAngle'].array()[0] + meta['min_viewcone_radius'] = u.deg * 0 + + info = MagicMcInfo( + obs_id = meta['obs_id'], + num_showers = meta['num_showers'], + energy_range_min = meta['energy_range_min'], + energy_range_max = meta['energy_range_max'], + spectral_index = meta['spectral_index'], + max_scatter_range = meta['max_scatter_range'], + min_scatter_range = meta['min_scatter_range'], + max_viewcone_radius = meta['max_viewcone_radius'], + min_viewcone_radius = meta['min_viewcone_radius'] + ) + + return info + + def to_df(self): + units = dict( + num_showers = u.one, + energy_range_min = u.TeV, + energy_range_max = u.TeV, + spectral_index = u.one, + max_scatter_range = u.m, + min_scatter_range = u.m, + max_viewcone_radius = u.deg, + min_viewcone_radius = u.deg, + ) + data = { + key: [ + self.__getattribute__(key).to(units[key]).value + ] + for key in units + } + return pandas.DataFrame(data=data) diff --git a/src/srcsim/mc.py b/src/srcsim/mc.py index c74eca3..9b57038 100644 --- a/src/srcsim/mc.py +++ b/src/srcsim/mc.py @@ -11,6 +11,12 @@ def power_law(e, e0, norm, index): class MCSample: + # Effective focal lengths + camscale = { + 1: 1 / 29.30565 * u.Unit('rad/m'), + 5: 1 / 18.2121 * u.Unit('rad/m'), + } + def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=None): self.units = dict( energy = u.TeV, @@ -18,10 +24,6 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No distance = u.m, viewcone = u.deg ) - - # TODO: refine this value - lst_focal_length = 28.01 * u.m - self.cam2angle = 1 * u.rad / lst_focal_length # TODO: refine the logic below / implement nicer if data_table is not None and config_table is not None: @@ -40,7 +42,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No self.tel_pos = SkyCoord(pointing_data['mc_az_tel'], pointing_data['mc_alt_tel'], unit=self.units['angle'], frame='altaz') # Working out the simulation spectrum - rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0] * self.units['distance'] + rmin, rmax = self.config_table[['min_scatter_range', 'max_scatter_range']].iloc[0].tolist() * self.units['distance'] ground_area = np.pi * (rmax**2 - rmin**2) nevents = self.config_table['num_showers'].iloc[0] * self.config_table['shower_reuse'].iloc[0] emin = self.config_table['energy_range_min'].iloc[0] * self.units['energy'] @@ -49,13 +51,13 @@ def __init__(self, file_name=None, obs_id=None, data_table=None, config_table=No self.spec_data = self.get_spec_data(nevents, emin, emax, index) self.spec_data['norm'] /= ground_area - cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle + cam_x, cam_y = self.data_table[['src_x', 'src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist()) self.evt_coord = SkyCoord(cam_x, cam_y, frame=self.tel_pos.skyoffset_frame()) self.evt_energy = self.data_table['mc_energy'].to_numpy() * self.units['energy'] # Filtering out events with excessive offsets (e.g. due to the simulation numerical accuracy) - offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0] * self.units['viewcone'] + offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone'] evt_offset = self.evt_coord.separation(self.tel_pos) in_fov = (evt_offset >= offset_min) & (evt_offset <= offset_max) @@ -106,6 +108,14 @@ def read_config(cls, file_name, obs_id): return pd.DataFrame(data=data) + def cam2angle(self, tel_id): + cam2angle = u.Quantity( + list( + map(lambda idx: self.camscale[idx], tel_id) + ) + ) + return cam2angle + def get_spec_data(self, n_events, emin, emax, index=-1): e0 = (emin * emax)**0.5 @@ -128,7 +138,7 @@ def dnde(self, energy): return power_law(energy, **self.spec_data) def dndo(self, coord): - offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0] * self.units['viewcone'] + offset_min, offset_max = self.config_table[['min_viewcone_radius', 'max_viewcone_radius']].iloc[0].tolist() * self.units['viewcone'] sky_area = 2 * np.pi * (np.cos(offset_min) - np.cos(offset_max)) * u.sr norm = 1 / sky_area diff --git a/src/srcsim/off.py b/src/srcsim/off.py index 6ab4533..cc5d357 100644 --- a/src/srcsim/off.py +++ b/src/srcsim/off.py @@ -7,6 +7,12 @@ class OffSample: + # Effective focal lengths + camscale = { + 1: 1 / 29.30565 * u.Unit('rad/m'), + 5: 1 / 18.2121 * u.Unit('rad/m'), + } + def __init__(self, file_name=None, obs_id=None, data_table=None): self.units = dict( energy = u.TeV, @@ -43,7 +49,7 @@ def __init__(self, file_name=None, obs_id=None, data_table=None): unit=self.units['angle'], frame='altaz' ) - cam_x, cam_y = self.data_table[['reco_src_x', 'reco_src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle + cam_x, cam_y = self.data_table[['reco_src_x', 'reco_src_y']].to_numpy().transpose() * self.units['distance'] * self.cam2angle(self.data_table['tel_id'].tolist()) self.evt_coord = SkyCoord(cam_x, cam_y, frame=tel_pos.skyoffset_frame()) self.evt_energy = self.data_table['reco_energy'].to_numpy() * self.units['energy'] @@ -72,6 +78,14 @@ def calc_obs_duration(self, data_table): return t_elapsed + def cam2angle(self, tel_id): + cam2angle = u.Quantity( + list( + map(lambda idx: self.camscale[idx], tel_id) + ) + ) + return cam2angle + def dndedo(self, energy, coord): dummy = 1 + 0 * energy.value val = dummy * 1 / (1 / self.obs_duration * u.Unit('1/(s TeV sr)')) diff --git a/src/srcsim/run.py b/src/srcsim/run.py index 46b0abe..45416f2 100644 --- a/src/srcsim/run.py +++ b/src/srcsim/run.py @@ -209,8 +209,8 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m # Reconstructed events coordinates reco_coords = SkyCoord( - evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle, - evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle, + evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle(evt['tel_id'].to_numpy()), + evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle(evt['tel_id'].to_numpy()), frame=offset_frame ) evt = evt.assign( diff --git a/src/srcsim/scripts/magic_to_dl2.py b/src/srcsim/scripts/magic_to_dl2.py new file mode 100644 index 0000000..5db1f6f --- /dev/null +++ b/src/srcsim/scripts/magic_to_dl2.py @@ -0,0 +1,95 @@ +import argparse +import logging +import healpy +import numpy + +from srcsim.magic.mc import MagicMcFile + + +def main(): + placeholder = "[HEALPY]" + + arg_parser = argparse.ArgumentParser( + description=""" + Convert MAGIC Melibea/Superstar MC files to the DL2 (HDF5) format. + """ + ) + + arg_parser.add_argument( + "-i", + "--input-file", + help='MAGIC MC file in Root format' + ) + arg_parser.add_argument( + "-n", + "--nside", + type=int, + default=0, + help='NSIDE to use in HEALPix split of the input file.' \ + ' If nside=0, no spliting is applied. Defaults to 0.' + ) + arg_parser.add_argument( + "-o", + "--output-file", + help='Output HDF5 file name.' \ + f' If nside > 0 is given, must containt a "{placeholder}" placeholder,' \ + ' that will be replaced with "healpy-nside-pixid" when writing the file.' + ) + arg_parser.add_argument( + '-v', + "--verbose", + action='store_true', + help='extra verbosity' + ) + args = arg_parser.parse_args() + + logging.basicConfig( + format='%(asctime)s %(name)-10s : %(levelname)-8s %(message)s', + datefmt='%Y-%m-%dT%H:%M:%S', + ) + + log = logging.getLogger(__name__) + + if args.verbose: + log.setLevel(logging.DEBUG) + else: + log.setLevel(logging.INFO) + + log.info(f'reading MC file') + mcfile = MagicMcFile.from_file(args.input_file) + log.debug(str(mcfile)) + + if args.nside > 0: + if not placeholder in args.output_file: + raise ValueError( + f'nside > 0 specified, but no "{placeholder}" given in output name "{args.output_file}"' + ) + + log.info('running HEALPix spliting') + log.info( + f'approximate pixel resolution: {numpy.degrees(healpy.nside2resol(args.nside)):.1f} deg' \ + f' (nside={args.nside})' + ) + files = mcfile.healpy_split(args.nside) + + npix = healpy.nside2npix(args.nside) + ids = range(npix) + with_events = tuple( + filter(lambda x: x[1].events_simulated.n_events > 0, zip(ids, files)) + ) + log.info(f'split complete: {len(with_events)} non-empty files (npix={npix})') + + log.info(f'writing the DL2 files') + for pix_id, mfile in with_events: + output_name = args.output_file.replace(placeholder, f'healpy-{args.nside}-{pix_id}') + log.debug(f'writing {output_name}') + mfile.write(output_name) + else: + log.info(f'writing the DL2 file') + mcfile.write(args.output_file) + + log.info('conversion complete') + + +if __name__ == '__main__': + main()