# uplim/views.py # from astropy_healpix import HEALPix does not have an option to # search for pixels non-inclusively import healpy as hp import astropy.units as u from astropy.coordinates import SkyCoord import scipy.special as sp from astropy.stats import poisson_conf_interval from django.db.models import Sum from rest_framework.views import APIView from rest_framework.response import Response from rest_framework import status from django.shortcuts import get_object_or_404 from uplim.models import Pixel # SURVEY PARAMETER PARSER # ************************************************************** def parse_survey_param(raw): surveys = set() for part in raw.split(","): if "-" in part: start, end = part.split("-", 1) surveys.update(range(int(start), int(end) + 1)) else: surveys.add(int(part)) return sorted(surveys) # PIXEL VIEW (MOSTLY FOR TESTING) # ************************************************************** class PixelAggregateView(APIView): def get(self, request): # GET PARAMETERS FROM THE QUERY # ************************************************************** raw_pixel = request.query_params.get("pixel") raw_survey = request.query_params.get("survey") # 400 BADREQUEST # ************************************************************** if raw_pixel is None or raw_survey is None: return Response( {"detail": "Both `pixel` and `survey` parameters are required."}, status=status.HTTP_400_BAD_REQUEST ) # FILTER THE INPUTS # ************************************************************** try: hpid = int(raw_pixel) except ValueError: return Response( {"detail": "`pixel` must be an integer."}, status=status.HTTP_400_BAD_REQUEST ) try: survey_numbers = parse_survey_param(raw_survey) except ValueError: return Response( {"detail": "Malformed `survey`; use N, N,M or N-M formats."}, status=status.HTTP_400_BAD_REQUEST ) # FILTER AND AGGREGATE # ************************************************************** qs = Pixel.objects.filter( hpid=hpid, survey__number__in=survey_numbers ) if not qs.exists(): # no matching pixel(s) → 404 get_object_or_404(Pixel, hpid=hpid, survey__number__in=survey_numbers) result = qs.aggregate( total_counts=Sum("counts"), total_exposure=Sum("exposure") ) # RETURN THE SUMS # ************************************************************** return Response(result, status=status.HTTP_200_OK) # class PixelDetailView(APIView): # """ # API endpoint that returns the pixel data (counts, exposure, rate) # for a given hpid. # """ # def get(self, request, hpid): # # Get the survey using the survey_number field. # # survey = get_object_or_404(Survey, number=survey_number) # # Get the pixel corresponding to the survey and hpid. # pixel = get_object_or_404(Pixel, hpid=hpid) # # Serialize the pixel data to JSON. # serializer = PixelSerializer(pixel) # return Response(serializer.data, status=status.HTTP_200_OK) # UPPER LIMIT COMPUTATION VIEW # ************************************************************** class UpperLimitView(APIView): """ Calculate confidence bounds based on aperture photometry using classic and bayesian methods """ def get(self, request): # GET PARAMETERS FROM THE REQUEST # ************************************************************** try: ra = float(request.query_params.get('ra')) dec = float(request.query_params.get('dec')) confidence_level = float(request.query_params.get('cl')) except (TypeError, ValueError): return Response( {"error": "Invalud parameters, provide RA, DEC, and CL"}, status = status.HTTP_400_BAD_REQUEST ) # ── NEW: pull & parse survey selection ── raw_survey = request.query_params.get('survey') if raw_survey is None: return Response( {"error": "Missing required `survey` parameter (e.g. ?survey=1,3-5)"}, status=status.HTTP_400_BAD_REQUEST ) try: survey_numbers = parse_survey_param(raw_survey) except ValueError: return Response( {"error": "Malformed `survey`; use formats like 1, 2-5, or 1,3-4"}, status=status.HTTP_400_BAD_REQUEST ) # hp = HEALPix( # nside = 4096, # order = 'ring', # frame = 'icrs' # ) # CREATE SKYCOORD, CONVERT TO GALACTIC BECAUSE THE HEALPIX # MAP ITSELF WAS MADE IN GALACTIC COORDINATES # ************************************************************** src_coord = SkyCoord( ra, dec, unit = 'deg', frame = 'icrs' ) gal = src_coord.galactic src_vec = hp.ang2vec(gal.l.deg, gal.b.deg, lonlat = True) # DEFINE APERTURE AND ANNULUS RADII # ************************************************************** aperture_radius = 71 # radius of the aperture in arc seconds # HPD ~48 arcseconds # 90% ~100 arcseconds annulus_inner = 142 # 2 * aperture_radius annulus_outer = 284 # 4 * aperture_radius # FETCH PIXEL DATA DEFINED VIA HP.QUERY_DISC (INCLUSIVE=FALSE) # ************************************************************** source_pixel_list = hp.query_disc( nside = 4096, vec = src_vec, inclusive = False, nest = False, radius = (aperture_radius * u.arcsecond).to(u.radian).value ) inner_pixel_list = hp.query_disc( nside = 4096, vec = src_vec, inclusive = False, nest = False, radius = (annulus_inner * u.arcsecond).to(u.radian).value ) outer_pixel_list = hp.query_disc( nside = 4096, vec = src_vec, inclusive = False, nest = False, radius = (annulus_outer * u.arcsecond).to(u.radian).value ) annulus_pixel_list = [ item for item in outer_pixel_list if item not in inner_pixel_list ] source_pixels = Pixel.objects.filter( hpid__in = source_pixel_list, survey__number__in = survey_numbers ) annulus_pixels = Pixel.objects.filter( hpid__in = annulus_pixel_list, survey__number__in = survey_numbers ) if not source_pixels.exists() and not annulus_pixels.exists(): return Response( {"detail": "No pixel data for the given survey selection."}, status=status.HTTP_404_NOT_FOUND ) # COMPUTE COUNTS, BACKGROUND ESTIMATE, EXPOSURE # ************************************************************** N = sum(obj.counts for obj in source_pixels) Nnpix = len(source_pixels) Bcounts = sum(obj.counts for obj in annulus_pixels) Bnpix = len(annulus_pixels) B = Bcounts / Bnpix * Nnpix tsum = sum(obj.exposure for obj in source_pixels) t = tsum / Nnpix # CONSTANTS # ************************************************************** EEF = .9 # eclosed energy fraction, .5 for hpd, .9 for w90 ECF = 4e-11 # energy conversion factor # BAYESIAN IMPLEMENTATION VIA POISSON_CONF_INTERVAL # ************************************************************** low, high = poisson_conf_interval( n = N, background = B, interval = 'kraft-burrows-nousek', confidence_level=confidence_level ) bayesian_count_ul = high bayesian_count_ll = low bayesian_rate_ul = bayesian_count_ul / t / EEF # count rate limits bayesian_rate_ll = bayesian_count_ll / t / EEF bayesian_flux_ul = bayesian_rate_ul * ECF # flux limits bayesian_flux_ll = bayesian_rate_ll * ECF # CLASSICAL IMPLEMENTATION VIA GAMMAINCCINV # **************************************************************** classic_count_ul = sp.gammainccinv(N+1, 1 - confidence_level) - B classic_count_ll = sp.gammainccinv(N, confidence_level) - B classic_count_ll = max(classic_count_ll, 0) classic_rate_ul = classic_count_ul / t / EEF # count rate limits classic_rate_ll = classic_count_ll / t / EEF classic_flux_ul = classic_rate_ul * ECF # flux limits classic_flux_ll = classic_rate_ll * ECF # FLUX ESTIMATION # **************************************************************** S = N - B # counts as simply counts within aperture # with the background estimate subtracted CR = S / t / EEF # count rate FL = CR * ECF # conversion to flux Flux = max(FL, 0) # flux cannot be lower than zero return Response({ 'ClassicUpperLimit' : classic_count_ul, 'ClassicLowerLimit' : classic_count_ll, 'ClassicCountRateUpperLimit' : classic_rate_ul, 'ClassicCountRateLowerLimit' : classic_rate_ll, 'ClassicFluxUpperLimit' : classic_flux_ul, 'ClassicFluxLowerLimit' : classic_flux_ll, 'BayesianUpperLimit' : bayesian_count_ul, 'BayesianLowerLimit' : bayesian_count_ll, 'BayesianCountRateUpperLimit' : bayesian_rate_ul, 'BayesianCountRateLowerLimit' : bayesian_rate_ll, 'BayesianFluxUpperLimit' : bayesian_flux_ul, 'BayesianFluxLowerLimit' : bayesian_flux_ll, 'FluxEstimate' : Flux, # 'N' : N, # 'Nnpix' : Nnpix, # 'Bcounts' : Bcounts, # 'Bnpix' : Bnpix, # 'B' : B, # 'tsum' : tsum, # 't' : t }, status=status.HTTP_200_OK)