@@ -0,0 +1,2388 @@
+#!/usr/bin/env python%0A%22%22%22%0AStacks boss spectra.%0A%22%22%22%0A%0Aimport argparse%0Aimport numpy as np%0Aimport h5py%0A%0Aimport qusp%0A%0Adef main():%0A    # parse command-line arguments%0A    parser = argparse.ArgumentParser(%0A        formatter_class=argparse.ArgumentDefaultsHelpFormatter)%0A    parser.add_argument(%22-i%22, %22--input%22, type=str, default=None,%0A        help=%22target list%22)%0A    parser.add_argument(%22-n%22, %22--ntargets%22, type=int, default=0,%0A        help=%22number of targets to use, 0 for all%22)%0A    parser.add_argument(%22-o%22, %22--output%22, type=str, default=None,%0A        help=%22output file name%22)%0A    parser.add_argument(%22--verbose%22, action=%22store_true%22,%0A        help=%22more verbose output%22)%0A    qusp.Paths.add_args(parser)%0A    args = parser.parse_args()%0A%0A    # setup boss data directory path%0A    paths = qusp.Paths(**qusp.Paths.from_args(args))%0A    # read target list%0A    targets = qusp.target.load_target_list(args.input)%0A%0A    # trim target list if requested%0A    ntargets = args.ntargets if args.ntargets %3E 0 else len(targets)%0A    targets = targets%5B:ntargets%5D%0A%0A    # initialize stack arrays%0A    npixels_fiducial = 4800%0A    wavelength = qusp.wavelength.get_fiducial_wavelength(np.arange(npixels_fiducial))%0A%0A    flux_wsum = np.zeros(npixels_fiducial)%0A    weight_sum = np.zeros_like(flux_wsum)%0A%0A    # loop over targets%0A    target_plate_generator = qusp.target.read_target_plates(paths.boss_path, targets, verbose=args.verbose)%0A    for target, spplate in target_plate_generator:%0A        # read this target's combined spectrum%0A        combined = qusp.read_combined_spectrum(spplate, target)%0A%0A        offset = qusp.wavelength.get_fiducial_pixel_index_offset(np.log10(combined.wavelength%5B0%5D))%0A%0A        indices = slice(offset,offset+combined.npixels)%0A%0A        flux_wsum%5Bindices%5D += combined.ivar*combined.flux%0A        weight_sum%5Bindices%5D += combined.ivar%0A%0A    flux_wmean = np.empty_like(flux_wsum)%0A    nonzero_weights = np.nonzero(weight_sum)%0A    flux_wmean%5Bnonzero_weights%5D = flux_wsum%5Bnonzero_weights%5D/weight_sum%5Bnonzero_weights%5D%0A%0A    if args.output:%0A        # save target list with sn column%0A        outfile = h5py.File(args.output+'.hdf5', 'w')%0A%0A        outfile.create_dataset('flux_wmean', data=flux_wmean)%0A        outfile.create_dataset('weight_sum', data=weight_sum)%0A        outfile.create_dataset('wavelength', data=wavelength)%0A        outfile.attrs%5B'ntargets'%5D = ntargets%0A%0A%0Aif __name__ == '__main__':%0A    main()%0A

add example program to stack spectra