Stats PESTO

Pourquoi un histogramme, pas une simple moyenne ?

Why a histogram, not a simple average?

Une EMCCD (Electron-Multiplying CCD) ne lit pas directement des « photons ». Chaque pixel, à chaque trame, subit une chaîne physique en trois étapes stochastiques : le nombre réel d'électrons cette trame-là suit une loi de Poisson de moyenne égale au flux réel ; ces électrons sont ensuite amplifiés par un registre à multiplication d'électrons (EM) selon une loi Gamma(n, gain) — un vrai processus stochastique, pas juste « n × gain » ; enfin, l'électronique de lecture ajoute un bruit gaussien (RON) et un décalage fixe (bias) pour que le convertisseur ne voie jamais de nombre négatif.

An EMCCD (Electron-Multiplying CCD) does not read out "photons" directly. Every pixel, every frame, goes through a three-stage stochastic physical chain: the actual number of electrons that frame follows a Poisson law with mean equal to the true flux; those electrons are then amplified by an electron-multiplying (EM) register following a Gamma(n, gain) law — a genuinely stochastic process, not just "n × gain"; finally the readout electronics add Gaussian noise (RON) and a fixed offset (bias) so the ADC never has to represent a negative number.

Aux flux typiques d'une caméra d'acquisition à haute cadence (ce projet : de l'ordre de 0,01–0,1 e⁻/trame pour le ciel, parfois quelques e⁻/trame au cœur d'une étoile brillante), une seule trame est presque entièrement du bruit de lecture — impossible de distinguer « 0 photon » de « 1 photon peu amplifié » à l'œil. C'est la forme de la distribution des valeurs ADU accumulée sur des centaines ou milliers de trames qui porte l'information : un pic gaussien net au niveau du bias (trames à zéro électron) et une petite traîne en forme de Gamma (trames à un ou plusieurs électrons amplifiés). Ajuster cette forme au maximum de vraisemblance — pas une simple moyenne ni un écrêtage sigma — est la seule façon de retrouver le flux exactement, même très en-dessous de 1 e⁻/trame, et cela reste tout aussi exact au-dessus de 1 e⁻/trame.

At the faint fluxes typical of a high-cadence acquisition camera (this project: roughly 0.01–0.1 e⁻/frame for the sky, occasionally a few e⁻/frame at a bright star's core), a single frame is almost pure read noise — you cannot tell "0 photons" from "1 photon that wasn't amplified much" by eye. What carries the information is the shape of the ADU distribution built up over hundreds or thousands of frames: a sharp Gaussian spike at the bias level (zero-electron frames) plus a small Gamma-shaped tail (amplified one- or multi-electron frames). Fitting that whole shape by maximum likelihood — not a mean, not a sigma-clipped average — is the only way to recover the flux exactly, even well below 1 e⁻/frame, and it stays just as exact above 1 e⁻/frame too.

La dérivation complète du modèle (Harpsøe et al. 2012, A&A 537, A50 ; Basden et al. 2003, MNRAS 345, 985) est utilisée telle quelle par ce projet via la bibliothèque sœur emccd_histo_fit.

The full model derivation (Harpsøe et al. 2012, A&A 537, A50; Basden et al. 2003, MNRAS 345, 985) is used as-is by this project via the sibling library emccd_histo_fit.

Les 7 étapes du pipeline

The pipeline's 7 stages

1. Configuration — pipeline_config.yaml, un seul fichier qui documente chaque paramètre.
2. Régions DS9 — définir à la souris les ouvertures photométriques (toi1452.reg), puis check_apertures.py impose un rayon commun.
3. Suivi de dérive — register_frames.py, position de l'étoile la plus brillante au fil du temps.
4. Enregistrement & empilement — register_and_stack.py (sous-pixel) et build_registered_cube.py (pixel entier).
5. Calibration du détecteur — fit_emccd_real_data.py et l'étape 3 de pixel_flux_map.py : bias, RON, gain EM.
6. Statistiques par pixel — pixel_flux_map.py : carte de flux moyen et carte d'incertitude, pixel par pixel.
7. Photométrie résolue — aperture_lightcurve.py : somme d'ouverture, une mesure par trame brute, sans soustraction de ciel.

1. Configuration — pipeline_config.yaml, one file documenting every parameter.
2. DS9 regions — point-and-click photometric apertures (toi1452.reg), then check_apertures.py enforces a common radius.
3. Drift tracking — register_frames.py, brightest-star position over time.
4. Registration & stacking — register_and_stack.py (sub-pixel) and build_registered_cube.py (whole pixel).
5. Detector calibration — fit_emccd_real_data.py and step 3 of pixel_flux_map.py: bias, RON, EM gain.
6. Per-pixel statistics — pixel_flux_map.py: mean-flux map and uncertainty map, pixel by pixel.
7. Resolved photometry — aperture_lightcurve.py: aperture sums, one measurement per raw frame, no sky subtraction.

Configuration — pipeline_config.yaml

Configuration — pipeline_config.yaml

Un seul fichier YAML rassemble et explique tous les paramètres du pipeline : où se trouvent les trames brutes sur votre machine, le seuil de détection, les réglages d'enregistrement, de calibration, de statistiques par pixel et de photométrie. Extrait :

One YAML file gathers and explains every pipeline parameter: where the raw frames live on your machine, the detection threshold, and the registration, calibration, per-pixel-statistics and photometry settings. Excerpt:

data:
  # Chemin local vers vos propres trames -- JAMAIS commité sur GitHub.
  raw_dir: "data/TOI1452Science"          # <-- À MODIFIER
  glob_pattern: "*.fits"
  date_header_keyword: "DATE"              # tri chronologique, pas par nom de fichier

detection:
  threshold_adu: 350   # ~10 sigma de bruit de lecture au-dessus du bias

registration:
  batch_size: 20        # trames par estimation de position
  cutout_half: 12        # demi-taille (pix) de la boîte pour l'ajustement gaussien 2D

calibration:
  sky_segmentation: {n_sky_sample: 20000, sky_sigma: 5.0}

pixel_stats:
  nbins: 20
  nmax_electrons: 40
  mu_max: 10.0           # e-/trame, borne haute de la grille de flux

La ligne raw_dir est la seule à modifier pour utiliser ce pipeline avec votre propre jeu de données : le dossier data/ est explicitement exclu du dépôt (.gitignore) — ce site publie le code et une poignée d'images de résultats, jamais les trames EMCCD elles-mêmes (souvent plusieurs Go, et des données d'observation propriétaires).

The raw_dir line is the only one you need to change to run this pipeline on your own dataset: the data/ folder is explicitly excluded from the repository (.gitignore) — this site publishes the code and a handful of result images, never the EMCCD frames themselves (often several GB, and proprietary observing data).

Régions DS9 — définir les ouvertures photométriques

DS9 regions — defining photometric apertures

aperture_lightcurve.py lit ses ouvertures circulaires depuis un fichier de région DS9 (.reg), la façon standard et robuste de marquer « où est l'étoile » sur une image sans coder les coordonnées en dur dans Python.

aperture_lightcurve.py reads its circular apertures from a DS9 region file (.reg), the standard, GUI-friendly way to mark "where the star is" on an image without hard-coding pixel coordinates in Python.

1. Ouvrir une trame ou, mieux, une image empilée dans SAOImageDS9 : ds9 figures/registered_binary_stack.fits
2. Region → Shape → Circle, puis cliquer-glisser un cercle sur l'étoile. Double-cliquer pour affiner le centre (x, y) et le rayon en pixels.
3. Répéter pour chaque ouverture voulue (une étoile cible, une étoile de comparaison, etc.).
4. Region → Save Regions…, format « DS9/Funtools », système de coordonnées « physical » (pas « wcs » — ce pipeline n'a pas de solution astrométrique, il travaille en pixels détecteur purs).

1. Open a frame, or better, a stacked image, in SAOImageDS9: ds9 figures/registered_binary_stack.fits
2. Region → Shape → Circle, then click-drag a circle over the star. Double-click to fine-tune the center (x, y) and radius in pixels.
3. Repeat for every aperture you need (a target star, a comparison star, etc.).
4. Region → Save Regions…, format "DS9/Funtools", coordinate system "physical" (not "wcs" — this pipeline has no astrometric solution, it works purely in detector pixel coordinates).

# Region file format: DS9 version 4.1
global color=green dashlist=8 3 width=1 ...
physical
circle(651.07706,146.99909,3)
circle(650.99917,140.37874,3)

Piège d'indexation : les coordonnées « physical » de DS9 sont indexées à partir de 1 (convention FITS), alors que tous les tableaux numpy de ce pipeline sont indexés à partir de 0. read_circles() soustrait 1 en x et en y à la lecture — sans ce correctif, chaque ouverture serait décalée d'un pixel.

Indexing gotcha: DS9 "physical" coordinates are 1-indexed (the FITS convention), while every numpy array in this pipeline is 0-indexed. read_circles() subtracts 1 from both x and y on load — without that correction every aperture would be off by one pixel.

Ouverture commune obligatoire — check_apertures.py

Enforced common aperture radius — check_apertures.py

Une ouverture circulaire ne capture jamais 100 % de la lumière d'une étoile — ce qu'elle rate dans les ailes de la PSF dépend de son rayon (la « correction d'ouverture »). Comparer deux sommes d'ouverture n'a de sens que si les deux ouvertures ratent la même fraction de lumière. C'est pourquoi aperture_lightcurve.py refuse de s'exécuter si les cercles d'un fichier .reg n'ont pas exactement le même rayon.

A circular aperture never captures 100% of a star's light — how much it misses in the PSF's wings depends on its radius (the "aperture correction"). Comparing two aperture sums only makes sense if both apertures miss the same fraction of light. This is why aperture_lightcurve.py refuses to run if the circles in a .reg file don't share the exact same radius.

check_apertures.py mesure, pour chaque ouverture, un profil de croissance (flux cumulatif soustrait du fond, en fonction du rayon) centré sur cette étoile, à partir de l'empilement enregistré. Le rayon suggéré est celui où ce profil atteint 90 % de sa valeur asymptotique — plafonné à 45 % de la distance à l'ouverture voisine la plus proche, pour ne jamais suggérer un rayon qui ferait se chevaucher deux ouvertures. Un prompt affiche ensuite la suggestion et demande confirmation avant de réécrire le fichier .reg (les centres ne sont jamais modifiés, seul le rayon change).

check_apertures.py measures, for every aperture, a curve of growth (background-subtracted cumulative flux vs. radius) centered on that star, using the registered stack. The suggested radius is where that curve reaches 90% of its asymptotic value — capped at 45% of the distance to the nearest other aperture, so it never suggests a radius that would make two apertures overlap. A prompt then shows the suggestion and asks for confirmation before rewriting the .reg file (centers are never touched, only the radius changes).

$ python check_apertures.py
Read 2 apertures from toi1452.reg:
  aperture 1: x=651.077  y=146.999  r=3.000 px  (DS9 physical)
  aperture 2: x=650.999  y=140.379  r=3.000 px  (DS9 physical)

All apertures already share a common radius: 3.0 px.

Nearest-neighbour aperture separation: 6.62 px -> capping any suggested
radius at 2.98 px to avoid overlapping apertures.
  aperture 1: curve-of-growth suggests r = 2.50 px (background = 113.10)
  aperture 2: curve-of-growth suggests r = 2.00 px (background = 123.24)

Suggested COMMON radius for every aperture: 2.5 px
Apply this common radius (2.5 px) to all apertures in toi1452.reg? [y/N]

Cet exemple montre le cas « déjà en accord » : les deux rayons de toi1452.reg étaient déjà égaux (3 px chacun) avant même de lancer l'outil. S'ils avaient différé, le même prompt aurait affiché un avertissement explicite avant la suggestion.

This example shows the "already agreeing" case: the two radii in toi1452.reg were already equal (3 px each) before ever running the tool. Had they differed, the same prompt would have printed an explicit warning before the suggestion.

Suivi de la dérive du pointage

Pointing drift tracking

register_frames.py coadditionne des lots de 20 trames (seuillées en binaire à 350 ADU, puis sommées — un « photon-counting » simple et robuste au bruit stochastique du gain EM), repère l'étoile la plus brillante et ajuste une gaussienne 2D pour un centroïde sous-pixellaire. Un spline cubique interpole ensuite ces ~50 mesures pour obtenir une position à l'instant exact de chaque trame individuelle.

register_frames.py coadds batches of 20 frames (binary-thresholded at 350 ADU, then summed — a simple photon-counting scheme robust to the EM gain's stochastic noise), locates the brightest star, and fits a 2D Gaussian for a sub-pixel centroid. A cubic spline then interpolates those ~50 measurements to get a position at the exact time of every individual frame.

Deux enregistrements distincts sont ensuite construits à partir de ce même spline, pour deux usages différents : un enregistrement sous-pixel (register_and_stack.py, par interpolation bilinéaire) pour l'imagerie et la photométrie d'ouverture, où l'on veut le meilleur alignement spatial possible ; et un enregistrement au pixel entier (build_registered_cube.py, simple découpage de tableau, aucune interpolation) pour les statistiques physiques par pixel, où chaque valeur ADU doit rester un échantillon réel et non un mélange de pixels voisins.

Two different registration outputs are then built from that same spline, for two different purposes: a sub-pixel registration (register_and_stack.py, via bilinear interpolation) for imaging and aperture photometry, where the best possible spatial alignment is wanted; and a whole-pixel registration (build_registered_cube.py, plain array slicing, no interpolation) for per-pixel physical statistics, where every ADU value must stay a real sample, not a blend of neighbouring pixels.

Chiffres clés

Key numbers

Cube enregistré & empilement seuillé

Registered cube & thresholded stack

Une fois la dérive corrigée, sommer les masques binaires (seuillés, enregistrés) de toutes les trames donne une image de champ propre, robuste au bruit de gain EM — utile à la fois comme aperçu visuel et comme cible pour dessiner les régions DS9 (section 3 ci-dessus). build_registered_cube.py construit en parallèle un cube 3D (trame, y, x) en ADU brut, au pixel entier, qui sert d'entrée à toutes les statistiques par pixel qui suivent.

Once the drift is corrected, summing the (thresholded, registered) binary masks of every frame gives a clean field image, robust to EM-gain noise — useful both as a visual preview and as the target image for drawing the DS9 regions (section 3 above). build_registered_cube.py builds, in parallel, a 3-D (frame, y, x) raw-ADU cube at whole-pixel registration, which feeds every per-pixel statistic that follows.

Cube brut : 1001 × 426 × 1024, uint16, 873 Mo — décalages entiers dx ∈ [−4, 4] px, dy ∈ [−1, 3] px, seulement 0,32 % de pixels de remplissage en bordure.

Raw cube: 1001 × 426 × 1024, uint16, 873 MB — integer shifts dx ∈ [−4, 4] px, dy ∈ [−1, 3] px, only 0.32% edge-fill pixels.

Pourquoi deux enregistrements ?

Why two registrations?

Un décalage sous-pixel (interpolation bilinéaire) donne le meilleur alignement spatial, mais chaque valeur devient un mélange pondéré de pixels voisins — d'où les valeurs fractionnaires visibles dans le stack ci-dessus. Un décalage entier garde chaque échantillon ADU intact : c'est ce qu'exige le modèle physique par pixel de la section suivante.

A sub-pixel shift (bilinear interpolation) gives the best spatial alignment, but every value becomes a weighted blend of neighbouring pixels — hence the fractional values visible in the stack above. A whole-pixel shift keeps every ADU sample intact, which is exactly what the per-pixel physical model in the next section requires.

Segmentation du ciel

Sky segmentation

Avant de calibrer le détecteur, il faut savoir quels pixels sont du ciel pur (aucune source persistante) et lesquels contiennent une étoile. pixel_flux_map.py seuille et somme le cube enregistré, puis utilise des statistiques robustes (écrêtage sigma itératif, astropy.stats.sigma_clipped_stats) sur cette carte de détection pour identifier le niveau de fond « normal » sans être biaisé par la poignée de pixels d'étoiles brillantes.

Before calibrating the detector, we need to know which pixels are pure sky (no persistent source) and which contain a star. pixel_flux_map.py thresholds and sums the registered cube, then uses robust statistics (iterative sigma-clipping, astropy.stats.sigma_clipped_stats) on that detection map to identify the "normal" background level without being biased by the handful of bright-star pixels.

Calibration du détecteur : bias, RON, gain EM

Detector calibration: bias, RON, EM gain

Le bias, le bruit de lecture (RON) et le gain EM sont des propriétés physiques de l'électronique, identiques pour tous les pixels. fit_emccd_real_data.py ajuste les quatre paramètres (bias, RON, gain, flux) simultanément par MCMC (paquet emcee) sur l'histogramme ADU d'une région fixe supposée être du ciel pur, et affiche le posterior conjoint sous forme de « triangle plot ». pixel_flux_map.py refait la même calibration physique, mais sur un masque de ciel détecté automatiquement (section précédente) et avec un simple point de maximum de vraisemblance (Nelder-Mead) — les deux méthodes indépendantes s'accordent étroitement, un bon test de cohérence.

Bias, read noise (RON), and EM gain are physical properties of the electronics, identical for every pixel. fit_emccd_real_data.py fits all four parameters (bias, RON, gain, flux) simultaneously via MCMC (the emcee package) on the ADU histogram of a fixed region assumed to be pure sky, and shows the joint posterior as a "triangle plot". pixel_flux_map.py reruns the same physical calibration, but on an automatically detected sky mask (previous section) and with a simple maximum-likelihood point estimate (Nelder-Mead) — the two independent methods agree closely, a useful consistency check.

Avant l'ajustement, tout échantillon de ciel inférieur à 50 % de la médiane de l'échantillon est écarté : le modèle physique ne peut produire aucune valeur aussi loin sous le bias (seule la traîne Gamma s'étend au-dessus), donc ces valeurs sont des artéfacts de lecture, pas de vrais échantillons du détecteur. Sur ce jeu de données : 63 894 échantillons (0,3 %) écartés sur 20 020 000.

Before fitting, any sky sample below 50% of the sample median is dropped: the physical model cannot produce a value that far below bias (only the Gamma tail extends above it), so such values are readout artifacts, not real detector samples. On this dataset: 63,894 samples (0.3%) dropped out of 20,020,000.

Carte de flux par pixel & carte d'incertitude

Per-pixel flux map & uncertainty map

Avec bias/RON/gain désormais fixés, le seul paramètre libre restant — le flux moyen μ en e⁻/trame — est ajusté indépendamment pour chaque pixel du champ, à partir de son propre histogramme sur 1001 trames. Ceci réutilise le solveur MLE vectorisé de la bibliothèque sœur emccd_histo_fit : un histogramme non uniforme (fin près du pic à zéro électron, plus grossier dans la traîne), un produit matriciel sur une grille dense de μ plutôt qu'un optimiseur par pixel — ce qui rend l'ajustement de ~436 000 pixels praticable en quelques secondes.

With bias/RON/gain now fixed, the one remaining free parameter — mean flux μ in e⁻/frame — is fit independently for every pixel in the field, from its own 1001-frame histogram. This reuses the sibling emccd_histo_fit library's vectorised MLE solver: a non-uniform histogram (fine near the zero-electron peak, coarser in the tail), a matrix product over a dense μ grid rather than a per-pixel optimizer — which is what makes fitting ~436,000 pixels tractable in a few seconds.

L'incertitude statistique de chaque pixel vient directement du critère du rapport de vraisemblance : les bornes μ_lo/μ_hi où la log-vraisemblance chute de 0,5 par rapport à son maximum (l'équivalent standard « 1σ », souvent asymétrique). La carte d'incertitude publiée est σ = (μ_hi − μ_lo) / 2, avec les bornes asymétriques exactes conservées à part.

Each pixel's statistical uncertainty comes directly from the likelihood-ratio criterion: the μ_lo/μ_hi bounds where the log-likelihood drops by 0.5 from its peak (the standard "1-sigma-equivalent", often asymmetric). The published uncertainty map is σ = (μ_hi − μ_lo) / 2, with the exact asymmetric bounds kept separately.

Attention : certains pixels dépassent bel et bien 1 e⁻/trame (le cœur d'une étoile brillante) — le modèle Poisson/Gamma complet et le bin de saturation restent exacts dans ce régime, contrairement à une approximation linéaire « flux bien en-dessous de 1 e⁻/trame ».

Note: some pixels genuinely exceed 1 e⁻/frame (a bright star's core) — the full Poisson/Gamma model and its saturation bin remain exact in that regime, unlike a "well below 1 e⁻/frame" linear approximation would be.

Photométrie résolue (ouverture) : seuillée vs somme simple

Resolved (aperture) photometry: thresholded vs. simple sum

aperture_lightcurve.py calcule deux photométries indépendantes par trame brute (la cadence la plus fine disponible, ~49 ms ici), pour les deux cercles DS9 (section 3), sans soustraction de ciel : la somme seuillée habituelle (comptage de photons, panneau du haut) et une somme simple, sans seuillage : chaque valeur ADU brute enregistrée est sommée (recouvrement fractionnaire exact, photutils), le bias est soustrait (un bias par pixel couvert, y compris fractionnaire) puis le résultat est divisé par le gain EM pour obtenir des électrons (panneau du bas).

aperture_lightcurve.py computes two independent photometries per raw frame (the finest cadence available, ~49 ms here), for the two DS9 circles (section 3), with no sky subtraction: the usual thresholded sum (photon counting, top panel) and a simple sum with no thresholding at all: every raw ADU value is summed (exact fractional-pixel overlap, photutils), the bias is subtracted (one bias per covered pixel, including fractional coverage), and the result divided by the EM gain to get electrons (bottom panel).

Pourquoi comparer les deux : le gain EM lui-même est un processus stochastique (Gamma, pas une simple multiplication — voir fit_emccd_real_data.py), donc sommer les amplitudes brutes importe ce « bruit d'excès » dans la mesure, alors que le seuillage l'évite en perdant l'information d'amplitude. On s'attend à ce que la méthode seuillée ait une dispersion relative (RMS/moyenne) plus faible.

Why compare the two: the EM gain itself is a stochastic process (Gamma, not a simple multiplication — see fit_emccd_real_data.py), so summing raw amplitudes imports that "excess noise" into the measurement, while thresholding avoids it at the cost of amplitude information. The thresholded method is expected to have the lower relative scatter (RMS/mean).

Sur ce jeu de données, la méthode seuillée est nettement moins bruitée (relativement) pour les deux ouvertures. Attention à l'explication : ce n'est PAS le « bruit d'excès » classique (facteur √2, Robbins & Hadwen), qui décrit la variance ajoutée par des événements multi-électrons — négligeables ici (~98 % des trames non nulles n'ont qu'un seul électron). L'écart vient plutôt de l'amplitude d'un seul électron amplifié : Gamma(1, gain) = Exponentielle(gain), dont le coefficient de variation vaut exactement 1 (aussi bruyant que son intensité moyenne). La somme simple utilise cette amplitude brute, donc en hérite ; le seuillage ne garde que « détecté ou non » et l'évite complètement.

On this dataset, the thresholded method is clearly less noisy (relatively) for both apertures. A note on why: this is NOT the classic "excess noise" factor (√2, Robbins & Hadwen), which describes variance added by multi-electron events — negligible here (~98% of non-empty frames have exactly one electron). The gap instead comes from a single amplified electron's own amplitude: Gamma(1, gain) = Exponential(gain), whose coefficient of variation is exactly 1 (as noisy as its own mean). The simple sum uses that raw amplitude and inherits this; thresholding keeps only "detected or not" and sidesteps it entirely.

Les deux cercles constituent la mesure complète, telle que demandée — aucune couronne de fond n'est estimée ni soustraite (le bias, une propriété de l'électronique, n'est pas un fond de ciel). Voir pipeline_config.yaml, section 7, pour ajouter une soustraction de ciel plus tard (une troisième région, plus grande, ou une couronne).

The two circles ARE the full measurement, exactly as requested — no background annulus is estimated or subtracted (the bias, an electronics property, is not a sky background). See pipeline_config.yaml, section 7, for how to add sky subtraction later (a third, larger region, or an annulus).

Le code

The code

Chaque script est commenté de façon très explicite — pas seulement le « quoi », mais le « pourquoi » physique et algorithmique à chaque étape — pour qu'on puisse le lire du début à la fin comme un tutoriel, pas seulement l'exécuter.

Every script is commented very explicitly — not just the "what", but the physical and algorithmic "why" at each step — so it can be read start to finish like a tutorial, not just executed.