diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 100644
index 0000000..273276c
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1,66 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project purpose
+
+Thermal analysis of a box passing through a production oven, measured with IR sensors. The goal is to find the mathematical law that best interpolates the temperature data in `data.csv`. The user directs the analysis and Claude implements the requested models in Python.
+
+## Setup
+
+```bat
+python -m venv venv
+venv\Scripts\activate
+pip install -r requirements.txt
+```
+
+## Running scripts
+
+```bat
+python <script_name>.py
+```
+
+Each script is standalone and produces a PNG plot. No test suite exists.
+
+## Data
+
+`data.csv` columns used for fitting:
+- `time since start [ms]` — converted to seconds in every script (`/1000`)
+- `temp_obj IR [C]` — box temperature, the fitting target
+- `temp_amb IR [C]` — ambient temperature, used to compute T∞
+
+All fits start from t = 105 s. The transition zone [115.9–117.2 s] is excluded by setting `sigma = 1e10` for those points (weighted least squares, not deletion). T∞ = 22.99°C is a fixed constant computed once in `calcola_t_amb.py` and hardcoded in the other scripts.
+
+## Fitting architecture
+
+All fits use `scipy.optimize.curve_fit(..., method="trf")` (Trust Region Reflective, nonlinear least squares). The base pattern shared across every fit script:
+
+```python
+popt, pcov = curve_fit(model, t_fit, T_fit, p0=[...], bounds=([...], [...]),
+                       sigma=sigma, absolute_sigma=True, method="trf")
+residuals = T_fit - model(t_fit, *popt)
+ss_res = np.sum((residuals / sigma_clean)**2)   # ignores excluded zone
+ss_tot = np.sum((T_fit[mask] - np.mean(T_fit[mask]))**2)
+R2 = 1 - ss_res / ss_tot
+perr = np.sqrt(np.diag(pcov))  # parameter uncertainties (1σ)
+```
+
+### Models implemented
+
+| Script | Model | Fixed params |
+|--------|-------|-------------|
+| `fit_raffreddamento_intero.py` | `T∞ + A·exp(-(t-t0)/τ)` — full range, weighted | T∞, t0=115.0s |
+| `fit_raffreddamento_1tratto.py` | same — segment [115.0–115.9s] | T∞, t0=115.0s |
+| `fit_raffreddamento_2tratto.py` | same — segment [117.5s–end] | T∞, t0=117.5s |
+| `fit_doppio_esponenziale.py` | `T∞ + A1·exp(-(t-t1)/τ1) + A2·exp(-(t-t2)/τ2)` | T∞, t1=115.0s, t2=117.5s |
+
+The double exponential (`fit_doppio_esponenziale.py`) is the best model: R² = 0.9991, decomposing cooling into a fast surface component (τ₁ ≈ 1.8 s) and a slow core component (τ₂ ≈ 17.6 s).
+
+## Adding a new fit
+
+1. Copy an existing fit script as a template.
+2. Define the model function with T∞ baked in or passed as a fixed parameter.
+3. Set `p0` from physical intuition or from a preliminary simpler fit.
+4. Use `sigma = np.where(mask_exclude, 1e10, 1.0)` to down-weight bad zones.
+5. Compute R² only over the unmasked points (`mask_clean`).
+6. Save output as a new PNG; do not overwrite existing plots.