325 lines
11 KiB
Python
325 lines
11 KiB
Python
# axc_ul/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 axc_ul.models import Pixel
|
|
#from axc_ul_server.axc_ul.serializers import PixelSerializer
|
|
|
|
|
|
# SURVEY PARAMETER PARSER
|
|
# **************************************************************
|
|
|
|
def parse_survey_param(raw: str) -> list[int]:
|
|
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) |