add reading of exposure info from coadds

py/picca/delta_extraction/astronomical_objects/desi_forest.py

@@ -37,7 +37,7 @@ class DesiForest(Forest):
     targetid: int
     Targetid of the object
 
-    tile: list of int
+    tileid: list of int
     Identifier of the tile used in the observation. None for no info
     """
     def __init__(self, **kwargs):
@@ -68,10 +68,20 @@ def __init__(self, **kwargs):
                                           "Missing variable 'targetid'")
         del kwargs["targetid"]
 
-        self.tile = []
-        if kwargs.get("tile") is not None:
-            self.tile.append(kwargs.get("tile"))
-            del kwargs["tile"]
+        self.tileid = []
+        if kwargs.get("tileid") is not None:
+            self.tileid.append(kwargs.get("tileid"))
+            del kwargs["tileid"]
+
+        self.fiber = []
+        if kwargs.get("fiber") is not None:
+            self.fiber.append(kwargs.get("fiber"))
+            del kwargs["fiber"]
+
+        self.expid = []
+        if kwargs.get("expid") is not None:
+            self.expid.append(kwargs.get("expid"))
+            del kwargs["expid"]
 
         # call parent constructor
         kwargs["los_id"] = self.targetid
@@ -98,7 +108,9 @@ def coadd(self, other):
                                           f"{type(other).__name__}")
         self.night += other.night
         self.petal += other.petal
-        self.tile += other.tile
+        self.tileid += other.tileid
+        self.expid += other.expid
+        self.fiber += other.fiber
         super().coadd(other)
 
     def get_header(self):
@@ -129,10 +141,20 @@ def get_header(self):
                 'comment': 'Observation petal(s)'
             },
             {
-                'name': 'TILE',
-                'value': "-".join(str(tile) for tile in self.tile),
+                'name': 'TILEID',
+                'value': "-".join(str(tileid) for tileid in self.tileid),
                 'comment': 'Observation tile(s)'
             },
+            {
+                'name': 'EXPID',
+                'value': "-".join(str(expid) for expid in self.expid),
+                'comment': 'Observation expid(s)'
+            },
+            {
+                'name': 'FIBER',
+                'value': "-".join(str(fiber) for fiber in self.fiber),
+                'comment': 'Observation fiber(s)'
+            },
         ]
 
         return header
@@ -150,9 +172,11 @@ def get_metadata(self):
         metadata = super().get_metadata()
         metadata += [
             self.targetid,
-            "-".join(str(night) for night in self.night),
-            "-".join(str(petal) for petal in self.petal),
-            "-".join(str(tile) for tile in self.tile),
+            "-".join(str(n) for night in self.night for n in night),
+            "-".join(str(p) for petal in self.petal for p in petal),
+            "-".join(str(t) for tileid in self.tileid for t in tileid),
+            "-".join(str(e) for expid in self.expid for e in expid),
+            "-".join(str(f) for fiber in self.fiber for f in fiber),
         ]
         return metadata
 
@@ -168,7 +192,8 @@ def get_metadata_dtype(cls):
         data
         """
         dtype = super().get_metadata_dtype()
-        dtype += [('TARGETID', int), ('NIGHT', 'S12'), ('PETAL', 'S12'), ('TILE', 'S12')]
+        dtype += [('TARGETID', int), ('NIGHT', 'S150'), ('PETAL', 'S150'),
+                  ('TILEID', 'S150'), ('EXPID', 'S150'),('FIBER', 'S150')]
         return dtype
 
     @classmethod
@@ -182,5 +207,5 @@ def get_metadata_units(cls):
         A list with the units of the line-of-sight data
         """
         units = super().get_metadata_units()
-        units += ["", "", "", ""]
+        units += ["", "", "", "", "", ""]
         return units

py/picca/delta_extraction/data_catalogues/desi_data.py

@@ -352,7 +352,8 @@ def format_data(self,
                     catalogue,
                     spectrographs_data,
                     targetid_spec,
-                    reso_from_truth=False):
+                    reso_from_truth=False,
+                    metadata_dict=None):
         """After data has been read, format it into DesiForest instances
 
         Instances will be DesiForest or DesiPk1dForest depending on analysis_type
@@ -399,6 +400,51 @@ def format_data(self,
                         f"for {targetid}")
             else:
                 w_t = w_t[0]
+            if metadata_dict is not None and not self.use_non_coadded_spectra:
+                exp_w_t = np.where(metadata_dict["EXP_TARGETID"] == targetid)[0]
+                expid = metadata_dict["EXP_EXPID"][exp_w_t]
+                night = metadata_dict["EXP_NIGHT"][exp_w_t]
+                petal = metadata_dict["EXP_PETAL"][exp_w_t]
+                tileid = metadata_dict["EXP_TILEID"][exp_w_t]
+                fiber = metadata_dict["EXP_FIBER"][exp_w_t]
+                metadata_dict_targetid = {'expid': expid,
+                                          'night': night,
+                                          'petal': petal,
+                                          'fiber': fiber,
+                                          'tileid': tileid}
+            elif metadata_dict is not None and self.use_non_coadded_spectra:
+                try:
+                    len(metadata_dict["EXPID"][w_t])
+                    expid = metadata_dict["EXPID"][w_t]
+                except TypeError:
+                    expid = [metadata_dict["EXPID"][w_t]]
+                try:
+                    len(metadata_dict["NIGHT"][w_t])
+                    night = metadata_dict["NIGHT"][w_t]
+                except TypeError:
+                    night = [metadata_dict["NIGHT"][w_t]]
+                try:
+                    len(metadata_dict["PETAL"][w_t])
+                    petal = metadata_dict["PETAL"][w_t]
+                except TypeError:
+                    petal = [metadata_dict["PETAL"][w_t]]
+                try:
+                    len(metadata_dict["FIBER"][w_t])
+                    fiber = metadata_dict["FIBER"][w_t]
+                except TypeError:
+                    fiber = [metadata_dict["FIBER"][w_t]]
+                try:
+                    len(metadata_dict["TILEID"][w_t])
+                    tileid = metadata_dict["TILEID"][w_t]
+                except TypeError:
+                    tileid = [metadata_dict["TILEID"][w_t]]
+                metadata_dict_targetid = {'expid': expid,
+                                          'night': night,
+                                          'petal': petal,
+                                          'fiber': fiber,
+                                          'tileid': tileid}
+            else:
+                metadata_dict_targetid = None
             # Construct DesiForest instance
             # Fluxes from the different spectrographs will be coadded
             for spec in spectrographs_data.values():
@@ -431,7 +477,8 @@ def format_data(self,
                                 ivar_i,
                                 w_t,
                                 reso_from_truth,
-                                num_data)
+                                num_data,
+                                metadata_dict=metadata_dict_targetid)
                 else:
                     forests_by_targetid, num_data  = self.update_forest_dictionary(
                             forests_by_targetid,
@@ -443,7 +490,8 @@ def format_data(self,
                             ivar,
                             w_t,
                             reso_from_truth,
-                            num_data)
+                            num_data,
+                            metadata_dict=metadata_dict_targetid)
         return forests_by_targetid, num_data
 
     def update_forest_dictionary(self,
@@ -456,7 +504,8 @@ def update_forest_dictionary(self,
                                  ivar,
                                  w_t,
                                  reso_from_truth,
-                                 num_data):
+                                 num_data,
+                                 metadata_dict=None):
         """Add new forests to the current forest dictonary
 
         Arguments
@@ -508,6 +557,8 @@ def update_forest_dictionary(self,
             "dec": row['DEC'],
             "z": row['Z'],
         }
+        if metadata_dict is not None:
+            args.update(metadata_dict)
         args["log_lambda"] = np.log10(spec['WAVELENGTH'])
 
 

py/picca/delta_extraction/data_catalogues/desi_healpix.py

@@ -208,6 +208,8 @@ def read_file(self, filename, catalogue):
             return {}, 0
         # Read targetid from fibermap to match to catalogue later
         fibermap = hdul['FIBERMAP'].read()
+        if not self.use_non_coadded_spectra:
+            exp_fibermap = hdul['EXP_FIBERMAP'].read()
 
         index_unique = np.full(fibermap.shape,True)
         if self.uniquify_night_targetid:
@@ -275,10 +277,38 @@ def _read_resolution(color, indices):
         if hdul_truth is not None:
             hdul_truth.close()
 
+        if not self.use_non_coadded_spectra:
+            exp_targetid = exp_fibermap['TARGETID']
+            exp_expid = exp_fibermap['EXPID']
+            exp_petal = exp_fibermap['PETAL_LOC']
+            exp_fiber = exp_fibermap['FIBER']
+            exp_night = exp_fibermap['NIGHT']
+            exp_tileid = exp_fibermap['TILEID']
+            metadata_dict = {'EXP_PETAL': exp_petal,
+                            'EXP_TILEID': exp_tileid,
+                            'EXP_NIGHT': exp_night,
+                            'EXP_EXPID': exp_expid,
+                            'EXP_FIBER': exp_fiber,
+                            'EXP_TARGETID': exp_targetid}
+        else:
+            #the indexing in this case is because the targetid is indexed the same way below
+            expid = fibermap['EXPID'][index_unique]
+            petal = fibermap['PETAL_LOC'][index_unique]
+            fiber = fibermap['FIBER'][index_unique]
+            night = fibermap['NIGHT'][index_unique]
+            tileid = fibermap['TILEID'][index_unique]
+
+            metadata_dict = {'PETAL': petal,
+                            'TILEID': tileid,
+                            'NIGHT': night,
+                            'EXPID': expid,
+                            'FIBER': fiber}
+
         forests_by_targetid, num_data = self.format_data(
             catalogue,
             spectrographs_data,
             fibermap["TARGETID"][index_unique],
-            reso_from_truth=reso_from_truth)
+            reso_from_truth=reso_from_truth,
+            metadata_dict=metadata_dict)
 
         return forests_by_targetid, num_data

py/picca/delta_extraction/data_catalogues/desi_healpix_fast.py

@@ -122,6 +122,9 @@ def read_data(self, is_mock=False):
                 self.forests = combine_results(imap_it)
                 t1 = time.time()
                 self.logger.progress(f"Time spent meerging threads: {t1-t0}")
+        else:
+            raise NotImplementedError('fast healpix reading is not implemented'
+                                      'for analyses with "num processors=1"')
 
         if len(self.forests) == 0:
             raise DataError("No quasars found, stopping here")
@@ -198,6 +201,7 @@ def read_file(self, filename, catalogue):
             return {}, 0
         # Read targetid from fibermap to match to catalogue later
         fibermap = hdul['FIBERMAP'].read()
+        exp_fibermap = hdul['EXP_FIBERMAP'].read()
 
         colors = ["B", "R"]
         if "Z_FLUX" in hdul:
@@ -222,11 +226,18 @@ def read_file(self, filename, catalogue):
             # It should be there by construction
             targetid = row["TARGETID"]
             w_t = np.where(fibermap["TARGETID"] == targetid)[0]
+            w_t_exp = np.where(exp_fibermap["TARGETID"] == targetid)[0]
             if len(w_t) == 0:
                 self.logger.warning(
                     f"Error reading {targetid}. Ignoring object")
                 continue
 
+            nights=exp_fibermap["NIGHT"][w_t_exp]
+            petals=exp_fibermap["PETAL_LOC"][w_t_exp]
+            fibers=exp_fibermap["FIBER"][w_t_exp]
+            tileids=exp_fibermap["TILEID"][w_t_exp]
+            expids=exp_fibermap["EXPID"][w_t_exp]
+
             w_t = w_t[0]
             # Construct DesiForest instance
             # Fluxes from the different spectrographs will be coadded
@@ -241,6 +252,11 @@ def read_file(self, filename, catalogue):
                 "dec": row['DEC'],
                 "z": row['Z'],
                 "log_lambda": log_lambda,
+                "night": nights,
+                "petal": petals,
+                "fiber": fibers,
+                "tileid": tileids,
+                "expid": expids,
             }
             forest = DesiForest(**args)
             forest.rebin()

py/picca/delta_extraction/data_catalogues/desi_tile.py

@@ -232,6 +232,9 @@ def read_file(self, filename, catalogue):
             return {}, 0
 
         fibermap = hdul['FIBERMAP'].read()
+        if not self.use_non_coadded_spectra:
+            exp_fibermap = hdul['EXP_FIBERMAP'].read()
+
 
         ra = fibermap['TARGET_RA']
         dec = fibermap['TARGET_DEC']
@@ -251,6 +254,33 @@ def read_file(self, filename, catalogue):
 
         petal_spec = fibermap['PETAL_LOC'][0]
 
+        if not self.use_non_coadded_spectra:
+            exp_targetid = exp_fibermap['TARGETID']
+            exp_expid = exp_fibermap['EXPID']
+            exp_petal = exp_fibermap['PETAL_LOC']
+            exp_fiber = exp_fibermap['FIBER']
+            exp_night = exp_fibermap['NIGHT']
+            exp_tileid = exp_fibermap['TILEID']
+            metadata_dict = {'EXP_PETAL': exp_petal,
+                            'EXP_TILEID': exp_tileid,
+                            'EXP_NIGHT': exp_night,
+                            'EXP_EXPID': exp_expid,
+                            'EXP_FIBER': exp_fiber,
+                            'EXP_TARGETID': exp_targetid}
+        else:
+            expid = fibermap['EXPID']
+            petal = fibermap['PETAL_LOC']
+            fiber = fibermap['FIBER']
+            night = fibermap['NIGHT']
+            tileid = fibermap['TILEID']
+
+            metadata_dict = {'PETAL': petal,
+                            'TILEID': tileid,
+                            'NIGHT': night,
+                            'EXPID': expid,
+                            'FIBER': fiber}
+
+
         spectrographs_data = {}
         for color in colors:
             try:
@@ -294,6 +324,7 @@ def read_file(self, filename, catalogue):
             catalogue[select],
             spectrographs_data,
             fibermap["TARGETID"],
+            metadata_dict=metadata_dict
         )
 
         return forests_by_targetid, num_data