diff --git a/GeVgal.csv b/GeVgal.csv new file mode 100644 index 0000000..a0ba3e1 --- /dev/null +++ b/GeVgal.csv @@ -0,0 +1,17 @@ +Name,l,b,Dist Mpc,log(M/M⊙) +NGC0628,138.617,-45.705,10.19,10.128 ± 0.136 +NGC0660,141.607,-47.347,11.57,10.098 ± 0.331 +NGC1291,247.524,-57.042,9.08,10.707 ± 0.136 +NGC1433,255.691,-51.195,9.04,10.070 ± 0.201 +NGC1512,248.668,-48.166,11.63,10.172 ± 0.160 +NGC1532,233.168,-46.584,14.26,10.528 ± 0.600 +NGC2903,208.710,44.540,8.87,10.404 ± 0.136 +NGC3368,234.435,57.010,10.42,10.523 ± 0.136 +NGC3877,150.719,65.956,14.63,10.096 ± 0.476 +NGC4192,265.434,74.960,12.68,10.371 ± 0.136 +NGC4666,299.538,62.368,14.70,10.298 ± 0.136 +NGC4818,305.212,54.323,11.04,10.008 ± 0.530 +NGC5248,335.929,68.751,13.75,10.264 ± 0.606 +NGC7331,93.722,-20.724,12.62,10.724 ± 0.327 +NGC7814,106.410,-45.175,14.40,10.520 ± 0.136 +PGC032861,245.103,55.513,14.45,12.827 ± 0.502 diff --git a/stack_request_test.ipynb b/stack_request_test.ipynb new file mode 100644 index 0000000..25b3fb3 --- /dev/null +++ b/stack_request_test.ipynb @@ -0,0 +1,170 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4cffd6c5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "RA List: [np.float64(24.1741435187182), np.float64(25.759787769012984), np.float64(49.32806570275772), np.float64(55.50613654382996), np.float64(60.97590646798826), np.float64(63.018113741788355), np.float64(143.0416839498853), np.float64(161.69014062163393), np.float64(176.53216723085026), np.float64(183.45178780664702), np.float64(191.28607761407375), np.float64(194.20387680952962), np.float64(204.38304472867193), np.float64(339.2669893687885), np.float64(0.8130105246497964), np.float64(164.03825008738184)]\n", + "Dec List: [np.float64(15.783869592105441), np.float64(13.645005280982105), np.float64(-41.10817408896737), np.float64(-47.22221653992977), np.float64(-43.34878852480052), np.float64(-32.87392378612486), np.float64(21.501827494459036), np.float64(11.819887017791151), np.float64(47.49512808371984), np.float64(14.900261322222237), np.float64(-0.4622668171007003), np.float64(-8.524938659400307), np.float64(8.885499320910338), np.float64(34.4157929700228), np.float64(16.14532447253357), np.float64(6.1729868291531)]\n" + ] + } + ], + "source": [ + "import csv\n", + "from astropy.coordinates import SkyCoord\n", + "import astropy.units as u\n", + "\n", + "# Initialize empty lists for RA and Dec\n", + "ra_list = []\n", + "dec_list = []\n", + "\n", + "# Define the path to your CSV file\n", + "csv_file_path = \"GeVgal.csv\"\n", + "\n", + "# Open and read the CSV file\n", + "with open(csv_file_path, 'r') as csvfile:\n", + " # Use csv.reader to handle the file, skipping the header row\n", + " csv_reader = csv.reader(csvfile)\n", + " next(csv_reader) # Skip the header row\n", + "\n", + " # Loop through each row in the CSV\n", + " for row in csv_reader:\n", + " try:\n", + " # Extract l and b, which are in the second and third columns (index 1 and 2)\n", + " l = float(row[1])\n", + " b = float(row[2])\n", + "\n", + " # Create a SkyCoord object with galactic coordinates\n", + " galactic_coord = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic')\n", + "\n", + " # Convert to ICRS (equatorial) coordinates\n", + " icrs_coord = galactic_coord.icrs\n", + "\n", + " # Append the RA and Dec values to the lists\n", + " ra_list.append(icrs_coord.ra.deg)\n", + " dec_list.append(icrs_coord.dec.deg)\n", + "\n", + " except (ValueError, IndexError) as e:\n", + " # Handle potential errors if a row is malformed\n", + " print(f\"Skipping a malformed row: {row} - Error: {e}\")\n", + "\n", + "# Now, ra_list and dec_list contain the converted coordinates\n", + "print(\"RA List:\", ra_list)\n", + "print(\"Dec List:\", dec_list)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff1d339a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Request successful!\n", + "{\n", + " \"Status\": 0,\n", + " \"ErrorMessage\": \"\",\n", + " \"ClassicUpperLimit\": 36.99647944311798,\n", + " \"ClassicLowerLimit\": 0.0,\n", + " \"ClassicOneSideUL\": 33.54440898622437,\n", + " \"ClassicCountRateUpperLimit\": 0.022452965897451615,\n", + " \"ClassicCountRateLowerLimit\": 0.0,\n", + " \"ClassicCountRateOneSideUL\": 0.020357922763323065,\n", + " \"ClassicFluxUpperLimit\": 7.504496105362505e-13,\n", + " \"ClassicFluxLowerLimit\": 0.0,\n", + " \"ClassicFluxOneSideUL\": 6.804265983763352e-13,\n", + " \"BayesianUpperLimit\": 33.83295926177776,\n", + " \"BayesianLowerLimit\": 0.10792639479099073,\n", + " \"BayesianOneSideUL\": 33.727107645421896,\n", + " \"BayesianCountRateUpperLimit\": 0.020533042385357955,\n", + " \"BayesianCountRateLowerLimit\": 6.549995291856849e-05,\n", + " \"BayesianCountRateOneSideUL\": 0.020468801604397097,\n", + " \"BayesianFluxUpperLimit\": 6.862796537256179e-13,\n", + " \"BayesianFluxLowerLimit\": 2.189216978388652e-15,\n", + " \"BayesianFluxOneSideUL\": 6.841325222832595e-13,\n", + " \"FluxEstimate\": 3.2382519031831067e-13,\n", + " \"ApertureCounts\": 94,\n", + " \"ApertureBackgroundCounts\": 78.03571428571428,\n", + " \"SourceCounts\": 15.964285714285722,\n", + " \"Exposure\": 2057.325295693534,\n", + " \"SourceRate\": 0.009688639787230046,\n", + " \"BackgroundRate\": 0.03793066388142964,\n", + " \"NormalizedBackgroundRate\": 1.6914547063541448e-05,\n", + " \"Contamination\": false,\n", + " \"CountMap\": {\n", + " \"healpix\": [],\n", + " \"counts\": [],\n", + " \"exposure\": [],\n", + " \"nside\": 4096,\n", + " \"order\": \"ring\",\n", + " \"radius_as\": 0.0\n", + " }\n", + "}\n" + ] + } + ], + "source": [ + "import requests\n", + "import json\n", + "\n", + "# Define the URL of your Django API endpoint\n", + "url = \"http://localhost:8000/api/stacked-upper-limit/\"\n", + "\n", + "payload = {\n", + " \"ra\": ra_list, # List of RA values\n", + " \"dec\": dec_list, # List of Dec values\n", + " \"cl\": 0.95, # A numeric confidence level\n", + " \"survey\": \"1-4\", # A string for the survey parameter\n", + " \"mr\": 0 # A numeric value for map_radius\n", + "}\n", + "\n", + "try:\n", + " # Send the PUT request with the JSON payload\n", + " response = requests.put(url, json=payload)\n", + "\n", + " # Check the response status code\n", + " response.raise_for_status() # This will raise an HTTPError for bad responses (4xx or 5xx)\n", + "\n", + " # Print the JSON response from the server\n", + " print(\"Request successful!\")\n", + " print(json.dumps(response.json(), indent=4))\n", + "\n", + "except requests.exceptions.HTTPError as err:\n", + " print(f\"HTTP Error: {err}\")\n", + " print(f\"Response body: {err.response.text}\")\n", + "except requests.exceptions.RequestException as err:\n", + " print(f\"An error occurred: {err}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv-pypy", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/urls.py b/urls.py index 7502b60..defab10 100644 --- a/urls.py +++ b/urls.py @@ -1,10 +1,19 @@ # uplim/urls.py from django.urls import path -from .views import PixelAggregateView, UpperLimitView # , PixelDetailView +from .views import ( + PixelAggregateView, + UpperLimitView, + StackedUpperLimitView, +) # , PixelDetailView urlpatterns = [ # path('pixel//', PixelDetailView.as_view(), name='pixel-detail'), path("pixel-aggregate/", PixelAggregateView.as_view(), name="pixel-aggregate"), path("upper-limit/", UpperLimitView.as_view(), name="upper-limit"), + path( + "stacked-upper-limit/", + StackedUpperLimitView.as_view(), + name="stacked-upper-limit", + ), ] diff --git a/views.py b/views.py index f8c7ecb..f11d1a7 100644 --- a/views.py +++ b/views.py @@ -24,6 +24,7 @@ from django.db.models import ( Value, ) from django.db.models.functions import Cast +from django.db.models import Q from django.shortcuts import get_object_or_404 from rest_framework.views import APIView from rest_framework.response import Response @@ -147,7 +148,7 @@ class UpperLimitView(APIView): confidence_level = float(request.query_params.get("cl")) except (TypeError, ValueError): return Response( - {"error": "Invalud parameters, provide RA, DEC, and CL"}, + {"error": "Invalid parameters, provide RA, DEC, and CL"}, status=status.HTTP_400_BAD_REQUEST, ) # pull & parse survey selection @@ -262,7 +263,7 @@ class UpperLimitView(APIView): # create and populate a list of average exposures per survey survey_averages = [] - for survey_id, exposures in t_by_survey.items(): + for _, exposures in t_by_survey.items(): average_exposure = sum(exposures) / len(exposures) survey_averages.append(average_exposure) @@ -367,7 +368,7 @@ class UpperLimitView(APIView): error_message = str(e) - N = Nnpix = Bcounts = Bnpix = B = t = S = CR = BR = 0.0 + N = Nnpix = Bcounts = Bnpix = B = t = S = CR = BR = NBR = 0.0 Flux = 0.0 classic_count_ul = classic_count_ll = classic_count_UL = 0.0 classic_rate_ul = classic_rate_ll = classic_rate_UL = 0.0 @@ -460,39 +461,15 @@ class UpperLimitView(APIView): .order_by("hpid") ) - # map_pixels_qs = ( - # Pixel.objects.filter(hpid__in=map_pixel_list, survey__in=survey_numbers) - # .values("hpid") - # .annotate( - # total_counts=Sum("counts"), - # max_contaminated_int=Max(Cast("contaminated", IntegerField())), - # ) - # .annotate( - # contaminated=Case( - # When(max_contaminated_int=1, then=Value(True)), - # default=Value(False), - # output_field=BooleanField(), - # ) - # ) - # .order_by("hpid") - # ) - # turn the queryset to a list map_pixels_list = list(map_pixels_qs) # get lists of healpix indices and count values map_healpix_list = [d["hpid"] for d in map_pixels_list] map_counts_list = [d["counts"] for d in map_pixels_list] + map_exposure_list = [d["exposure"] for d in map_pixels_list] # map_contaminated_list = [d["contaminated"] for d in map_pixels_list] - # cont_dict = dict( - # Pixel.objects.filter(hpid__in=map_healpix_list, survey__in=survey_numbers) - # .values_list("hpid", "contaminated") - # .distinct() - # ) - - # map_contaminated_list = [cont_dict[h] for h in map_healpix_list] - # set map nside map_nside = 4096 @@ -503,6 +480,7 @@ class UpperLimitView(APIView): map_dict = { "healpix": map_healpix_list, "counts": map_counts_list, + "exposure": map_exposure_list, # "contaminated": map_contaminated_list, "nside": map_nside, "order": map_order, @@ -558,3 +536,475 @@ class UpperLimitView(APIView): clean = sanitize(result) # calling sanitize() to convert NaN to null return Response(clean, status=status.HTTP_200_OK) + + +class StackedUpperLimitView(APIView): + """ + Calculate confidence bounds based on aperture photometry using classic and bayesian methods for a set of sources + """ + + def put(self, request): + + data = request.data + + try: + ra_list = data.get("ra", []) + dec_list = data.get("dec", []) + confidence_level = data.get("cl") + + if not isinstance(ra_list, list) or not all( + isinstance(x, (float, int)) for x in ra_list + ): + raise TypeError + if not isinstance(dec_list, list) or not all( + isinstance(x, (float, int)) for x in dec_list + ): + raise TypeError + if not isinstance(confidence_level, (float, int)): + raise TypeError + + except TypeError: + return Response( + { + "error": "Invalid parameters, provide 'ra', 'dec' lists of numbers, and a numeric 'cl'" + }, + status=status.HTTP_400_BAD_REQUEST, + ) + + raw_survey = data.get("survey") + + if raw_survey is None: + return Response( + {"error": "Missing required `survey` parameter"}, + 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, + ) + + map_radius_value = data.get("mr") + + # DEFINE APERTURE AND ANNULUS RADII + # ************************************************************** + + aperture_radius = 71 # radius of the aperture in arc seconds + # HPD ~48 arcseconds + # 90% ~100 arcseconds + annulus_inner = 213 # 3 * aperture_radius + annulus_outer = 355 # 5 * aperture_radius + + # CREATE SKYCOORD, CONVERT TO GALACTIC BECAUSE THE HEALPIX + # MAP ITSELF WAS MADE IN GALACTIC COORDINATES + # ************************************************************** + + # a list of skycoords + src_coord_list = SkyCoord(ra_list, dec_list, unit="deg", frame="icrs") + gal_list = src_coord_list.galactic + # and then vectors + src_vec_list = hp.ang2vec(gal_list.l.deg, gal_list.b.deg, lonlat=True) + + # FETCH PIXEL DATA DEFINED VIA HP.QUERY_DISC (INCLUSIVE=FALSE) + # ************************************************************** + + # first init pixel lists + source_pixel_dict, inner_pixel_dict, outer_pixel_dict = {}, {}, {} + + # then loop over vectors appending indices to lists + for index, src_vec in enumerate(src_vec_list): + + source_pixel_dict[index] = hp.query_disc( + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(aperture_radius * u.arcsecond).to(u.radian).value, + ) + + inner_pixel_dict[index] = hp.query_disc( + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(annulus_inner * u.arcsecond).to(u.radian).value, + ) + + outer_pixel_dict[index] = hp.query_disc( + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(annulus_outer * u.arcsecond).to(u.radian).value, + ) + + # flatten dicts into lists for most of math + + outer_pixel_list = [ + item for sublist in outer_pixel_dict.values() for item in sublist + ] + + inner_pixel_list = [ + item for sublist in inner_pixel_dict.values() for item in sublist + ] + + source_pixel_list = [ + item for sublist in source_pixel_dict.values() for item in sublist + ] + + # proceed as is + 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__in=survey_numbers + ) + + annulus_pixels = Pixel.objects.filter( + hpid__in=annulus_pixel_list, survey__in=survey_numbers + ) + + # check contamination + contamination = ( + source_pixels.filter(contaminated=True).exists() + or annulus_pixels.filter(contaminated=True).exists() + ) + + # exclude contaminated pixels from the background and source regions + annulus_pixels = annulus_pixels.exclude(contaminated=True) + # source_pixels = source_pixels.exclude(contaminated=True) + + status_int = 0 + error_message = "" + + if not source_pixels.exists() or not annulus_pixels.exists(): + status_int = 1 # status 1 if there are no source or bg pixels + + # COMPUTE COUNTS, BACKGROUND ESTIMATE, EXPOSURE + # ************************************************************** + try: + # summing counts across all surveys + 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 + + t = 0 + + for source_pixel_sublist in source_pixel_dict.values(): + + source_pixels = Pixel.objects.filter( + hpid__in=source_pixel_sublist, survey__in=survey_numbers + ) + + # create a dict of exposures keyed by survey + t_by_survey = defaultdict(list) + + for pixel in source_pixels: + t_by_survey[pixel.survey].append(pixel.exposure) + + # create and populate a list of average exposures per survey + survey_averages = [] + + for _, exposures in t_by_survey.items(): + average_exposure = sum(exposures) / len(exposures) + survey_averages.append(average_exposure) + + # sum them up across surveys + t += sum(survey_averages) + + # CONSTANTS + # ************************************************************** + + # EEF = .9 # eclosed energy fraction, .5 for hpd, .9 for w90 + # ECF = 4e-11 # energy conversion factor + + EEF = 0.80091 # use values from the paper + ECF = 3.3423184e-11 + + # BAYESIAN IMPLEMENTATION VIA POISSON_CONF_INTERVAL + # ************************************************************** + + # double sided interval + low, high = poisson_conf_interval( + n=N, + background=B, + interval="kraft-burrows-nousek", + confidence_level=confidence_level, + ) + + # because poisson_conf_interval returns lists + bayesian_count_ul = high[0] + bayesian_count_ll = low[0] + + # compute upper limit + bayesian_count_UL = ( + sp.gammaincinv( + N + 1, + confidence_level * sp.gammaincc(N + 1, B) + sp.gammainc(N + 1, B), + ) + - B + ) + + bayesian_rate_ul = bayesian_count_ul / t / EEF # count rate limits + bayesian_rate_ll = bayesian_count_ll / t / EEF + bayesian_rate_UL = bayesian_count_UL / t / EEF + + bayesian_flux_ul = bayesian_rate_ul * ECF # flux limits + bayesian_flux_ll = bayesian_rate_ll * ECF + bayesian_flux_UL = bayesian_rate_UL * ECF + + # CLASSICAL IMPLEMENTATION VIA GAMMAINCCINV + # **************************************************************** + + # confidence level for the double sided interval + cl2 = (1 + confidence_level) / 2 + # double sided interval + classic_count_ul = sp.gammainccinv(N + 1, 1 - cl2) - B + classic_count_ll = sp.gammainccinv(N, cl2) - B + # one sided interval + classic_count_UL = sp.gammainccinv(N + 1, 1 - confidence_level) - B + + if not np.isfinite(classic_count_ll) or classic_count_ll < 0: + classic_count_ll = 0.0 + + classic_rate_ul = classic_count_ul / t / EEF # count rate limits + classic_rate_ll = classic_count_ll / t / EEF + classic_rate_UL = classic_count_UL / t / EEF + + classic_flux_ul = classic_rate_ul * ECF # flux limits + classic_flux_ll = classic_rate_ll * ECF + classic_flux_UL = classic_rate_UL * ECF + + # FLUX ESTIMATION + # **************************************************************** + + S = N - B # counts as simply counts within aperture + # with the background estimate subtracted + + CR = S / t / EEF # source count rate + + BR = B / t # background rate within aperture + + FL = CR * ECF # conversion to flux + + Flux = max(FL, 0) # flux cannot be lower than zero + + # NORMALIZED BACKGROUND RATE + # **************************************************************** + # in units of ct/s/keV/arcmin^2 + + sr_to_arcmin2 = (180 / np.pi * 60) ** 2 + + pix_area = ( + hp.nside2pixarea(nside=4096) * sr_to_arcmin2 + ) # nside2pixarea yields area in sr + + bg_area = pix_area * Bnpix # get total bg region area + + en_range = 12 - 4 + + NBR = Bcounts / t / en_range / bg_area + + # handle exceptions: write error text to the response + # zero out all math results + except Exception as e: + + error_message = str(e) + + N = Nnpix = Bcounts = Bnpix = B = t = S = CR = BR = 0.0 + Flux = 0.0 + classic_count_ul = classic_count_ll = classic_count_UL = 0.0 + classic_rate_ul = classic_rate_ll = classic_rate_UL = 0.0 + classic_flux_ul = classic_flux_ll = classic_flux_UL = 0.0 + bayesian_count_ul = bayesian_count_ll = bayesian_count_UL = 0.0 + bayesian_rate_ul = bayesian_rate_ll = bayesian_rate_UL = 0.0 + bayesian_flux_ul = bayesian_flux_ll = bayesian_flux_UL = 0.0 + + # NEARBY SOURCES CHECK + # **************************************************************** + # disabled for now, many-to-many separation computation is expensive + # if map_radius == 0: + # radius_as = 5 * aperture_radius + # else: + # radius_as = map_radius * 2 + + # radius_deg = 3 # radius_as / 3600 + + # # dec_min_list = max(dec_list - radius_deg, -90) + # # dec_max_list = min(dec_list + radius_deg, 90) + + # dec_min_list = [max(d - radius_deg, -90) for d in dec_list] + # dec_max_list = [min(d + radius_deg, 90) for d in dec_list] + + # # cheap belt query + # # belt_sources = CatalogSource.objects.filter( + # # dec_deg__gte=dec_min, dec_deg__lte=dec_max + # # ) + + # query = Q() + + # for dec_min, dec_max in zip(dec_min_list, dec_max_list): + # query |= Q(dec_deg__gte=dec_min, dec_deg__lte=dec_max) + + # belt_list_sources = CatalogSource.objects.filter(query) + + # center_coord = SkyCoord(ra_list, dec_list, unit="deg") + + # nearby_sources = [] + + # radius_map = {0: 0.06, 125: 0.15, 250: 0.5, 2000: 0.9, 20000: 2.5} + + # sorted_bounds = sorted(radius_map.keys()) + + # # refine belt to circular region using astropy separation + # for catsrc in belt_list_sources: + # catsrc_coord = SkyCoord(catsrc.ra_deg, catsrc.dec_deg, unit="deg") + # separation = center_coord.separation(catsrc_coord) + # if separation.deg > radius_deg: + # continue + # f = catsrc.flux or 0.0 + # for lb in reversed(sorted_bounds): + # if f >= lb: + # mask_radius = radius_map[lb] * 3600 + # break + + # nearby_sources.append( + # { + # "srcid": catsrc.srcid, + # "name": catsrc.name, + # "ra_deg": catsrc.ra_deg, + # "dec_deg": catsrc.dec_deg, + # "pos_error": catsrc.pos_error, + # "significance": catsrc.significance, + # "flux": catsrc.flux, + # "flux_error": catsrc.flux_error, + # "catalog_name": catsrc.catalog_name, + # "new_xray": catsrc.new_xray, + # "source_type": catsrc.source_type, + # "mask_radius_as": mask_radius, + # "separation_as": separation.arcsecond, + # } + # ) + # nearby_sources.sort(key=lambda src: src["separation_as"]) + + # REGION IMAGE SERVING + # **************************************************************** + # + # default value if not specified in the query + # get hpids within a circle of radius sqrt(2) * outer annulus radius + if map_radius_value is None: + map_radius = annulus_outer * np.sqrt(2) + else: + map_radius = float(map_radius_value) + + # iterate over source vectors, making maps + + # Create a single list to hold all Healpix IDs + all_map_pixels = [] + + # Iterate over source vectors and get all unique Healpix IDs + for index, src_vec in enumerate(src_vec_list): + pixels = hp.query_disc( + nside=4096, + vec=src_vec, + inclusive=False, + nest=False, + radius=(map_radius * u.arcsecond).to(u.radian).value, + ) + all_map_pixels.extend(pixels) + + # Get only unique Healpix IDs + unique_map_pixels = list(set(all_map_pixels)) + + # Perform a single database query for all unique pixels + map_pixels_qs = ( + Pixel.objects.filter(hpid__in=unique_map_pixels, survey__in=survey_numbers) + .values("hpid") + .annotate(counts=Sum("counts"), exposure=Sum("exposure")) + .order_by("hpid") + ) + + # turn the queryset to a list + map_pixels_list = list(map_pixels_qs) + + # get lists of healpix indices and count values + map_healpix_list = [d["hpid"] for d in map_pixels_list] + map_counts_list = [d["counts"] for d in map_pixels_list] + map_exposure_list = [d["exposure"] for d in map_pixels_list] + # map_contaminated_list = [d["contaminated"] for d in map_pixels_list] + + # set map nside + map_nside = 4096 + + # set map order + map_order = "ring" + + # assemble the result dict + map_dict = { + "healpix": map_healpix_list, + "counts": map_counts_list, + "exposure": map_exposure_list, + # "contaminated": map_contaminated_list, + "nside": map_nside, + "order": map_order, + "radius_as": map_radius, + } + + # RESULT JSON + # **************************************************************** + + result = { + "Status": status_int, + # 0 OK + # 1 either source or bg pixels missing + "ErrorMessage": error_message, + # frequentist limits + "ClassicUpperLimit": classic_count_ul, + "ClassicLowerLimit": classic_count_ll, + "ClassicOneSideUL": classic_count_UL, + "ClassicCountRateUpperLimit": classic_rate_ul, + "ClassicCountRateLowerLimit": classic_rate_ll, + "ClassicCountRateOneSideUL": classic_rate_UL, + "ClassicFluxUpperLimit": classic_flux_ul, + "ClassicFluxLowerLimit": classic_flux_ll, + "ClassicFluxOneSideUL": classic_flux_UL, + # bayesian limits + "BayesianUpperLimit": bayesian_count_ul, + "BayesianLowerLimit": bayesian_count_ll, + "BayesianOneSideUL": bayesian_count_UL, + "BayesianCountRateUpperLimit": bayesian_rate_ul, + "BayesianCountRateLowerLimit": bayesian_rate_ll, + "BayesianCountRateOneSideUL": bayesian_rate_UL, + "BayesianFluxUpperLimit": bayesian_flux_ul, + "BayesianFluxLowerLimit": bayesian_flux_ll, + "BayesianFluxOneSideUL": bayesian_flux_UL, + # flux 'center value' estimate + "FluxEstimate": Flux, + # raw data + "ApertureCounts": N, + "ApertureBackgroundCounts": B, + "SourceCounts": S, + "Exposure": t, + # count rates + "SourceRate": CR, + "BackgroundRate": BR, + "NormalizedBackgroundRate": NBR, # ct/s/keV/arcmin2 + # contamination + "Contamination": contamination, + # "NearbySources": nearby_sources, + # count map for the frontend image + "CountMap": map_dict, + } + + clean = sanitize(result) # calling sanitize() to convert NaN to null + + return Response(clean, status=status.HTTP_200_OK)