diff --git a/analyze_qwenapp_trace.py b/analyze_qwenapp_trace.py
new file mode 100644
index 0000000..6942e3c
--- /dev/null
+++ b/analyze_qwenapp_trace.py
@@ -0,0 +1,339 @@
+#!/usr/bin/env python3
+"""Generate all characterization CSV/JSON data for the qwenapp trace.
+
+Reads the formatted qwenapp trace (release format with hash_ids) and produces:
+  1. token_length_cdf_qwenapp.csv + .stats.json
+  2. f2a_block_lifetime_cdf_qwenapp.csv + .stats.json
+  3. f2a_reuse_gap_cdf_qwenapp.csv
+  4. f2a_reuse_timeline_qwenapp.csv
+  5. f2b_session_skew_cdf_qwenapp.json
+
+Run on dash1:
+    python paper/data/analyze_qwenapp_trace.py \
+        ~/ali-trace/trace-qwenapp-formatted/0529-1000-1200-ttp.jsonl
+"""
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+import sys
+from collections import Counter, defaultdict
+from pathlib import Path
+
+import numpy as np
+
+
+def _percentile(sorted_arr, p):
+    idx = int(len(sorted_arr) * p / 100.0)
+    return sorted_arr[min(idx, len(sorted_arr) - 1)]
+
+
+def _log_bins(lo, hi, n=200):
+    return np.logspace(np.log10(max(lo, 1e-4)), np.log10(hi), n)
+
+
+# ---------------------------------------------------------------------------
+# 1. Token length CDF
+# ---------------------------------------------------------------------------
+def build_token_length(rows, sessions, outdir):
+    input_toks = sorted(r["input_length"] for r in rows)
+    output_toks = sorted(r["output_length"] for r in rows)
+
+    append_toks = []
+    for sid, sess_rows in sessions.items():
+        if len(sess_rows) < 2:
+            continue
+        for i in range(1, len(sess_rows)):
+            delta = sess_rows[i]["input_length"] - sess_rows[i - 1]["input_length"]
+            append_toks.append(delta)
+    append_toks.sort()
+
+    n = max(len(input_toks), len(output_toks), len(append_toks))
+    cdf_points = np.linspace(0, 1, min(n, 2000))
+
+    with (outdir / "token_length_cdf_qwenapp.csv").open("w", newline="") as f:
+        w = csv.writer(f)
+        w.writerow(["cdf", "input_tokens", "append_tokens", "output_tokens"])
+        for c in cdf_points:
+            ii = min(int(c * len(input_toks)), len(input_toks) - 1)
+            ai = min(int(c * len(append_toks)), len(append_toks) - 1) if append_toks else 0
+            oi = min(int(c * len(output_toks)), len(output_toks) - 1)
+            w.writerow([
+                f"{c:.3f}",
+                input_toks[ii],
+                append_toks[ai] if append_toks else 0,
+                output_toks[oi],
+            ])
+
+    n_root = sum(1 for sid, sr in sessions.items() for i, r in enumerate(sr) if i == 0)
+    n_neg = sum(1 for a in append_toks if a < 0)
+    stats = {
+        "trace": "trace-qwenapp-formatted/0529-1000-1200.jsonl",
+        "n_records": len(rows),
+        "input": {
+            "n": len(input_toks), "p50": _percentile(input_toks, 50),
+            "p90": _percentile(input_toks, 90), "p99": _percentile(input_toks, 99),
+            "mean": sum(input_toks) / len(input_toks),
+            "min": input_toks[0], "max": input_toks[-1],
+        },
+        "append": {
+            "n": len(append_toks),
+            "p50": _percentile(append_toks, 50) if append_toks else 0,
+            "p90": _percentile(append_toks, 90) if append_toks else 0,
+            "p99": _percentile(append_toks, 99) if append_toks else 0,
+            "mean": sum(append_toks) / len(append_toks) if append_toks else 0,
+            "min": append_toks[0] if append_toks else 0,
+            "max": append_toks[-1] if append_toks else 0,
+            "n_root": n_root, "n_negative": n_neg,
+        },
+        "output": {
+            "n": len(output_toks), "p50": _percentile(output_toks, 50),
+            "p90": _percentile(output_toks, 90), "p99": _percentile(output_toks, 99),
+            "mean": sum(output_toks) / len(output_toks),
+            "min": output_toks[0], "max": output_toks[-1],
+        },
+    }
+    (outdir / "token_length_cdf_qwenapp.stats.json").write_text(
+        json.dumps(stats, indent=4, ensure_ascii=False))
+    print(f"  token_length: {len(rows)} rows, append={len(append_toks)}")
+
+
+# ---------------------------------------------------------------------------
+# 2. Block lifetime CDF
+# ---------------------------------------------------------------------------
+def build_block_lifetime(rows, outdir):
+    block_first: dict[int, float] = {}
+    block_last: dict[int, float] = {}
+    for r in rows:
+        ts = r["timestamp"]
+        for h in r["hash_ids"]:
+            if h not in block_first:
+                block_first[h] = ts
+            block_last[h] = ts
+
+    lifetimes = []
+    for h in block_first:
+        lt = block_last[h] - block_first[h]
+        if lt > 0:
+            lifetimes.append(lt)
+    lifetimes.sort()
+
+    if not lifetimes:
+        print("  block_lifetime: no reused blocks!")
+        return
+
+    bins = np.concatenate([
+        _log_bins(0.01, lifetimes[-1] * 1.01, 300),
+        [float("inf")],
+    ])
+    counts, _ = np.histogram(lifetimes, bins=bins)
+
+    with (outdir / "f2a_block_lifetime_cdf_qwenapp.csv").open("w", newline="") as f:
+        w = csv.writer(f)
+        w.writerow(["bucket_upper_s", "count"])
+        for edge, cnt in zip(bins[1:], counts):
+            w.writerow([edge, int(cnt)])
+
+    stats = {
+        "n_records": len(rows),
+        "n_distinct_blocks": len(block_first),
+        "n_reused_blocks": len(lifetimes),
+        "reused_fraction": len(lifetimes) / len(block_first) if block_first else 0,
+        "p50_s": _percentile(lifetimes, 50),
+        "p90_s": _percentile(lifetimes, 90),
+        "p99_s": _percentile(lifetimes, 99),
+        "mean_s": sum(lifetimes) / len(lifetimes),
+        "max_s": lifetimes[-1],
+    }
+    (outdir / "f2a_block_lifetime_cdf_qwenapp.stats.json").write_text(
+        json.dumps(stats, indent=4, ensure_ascii=False))
+    print(f"  block_lifetime: {len(lifetimes)} reused blocks / {len(block_first)} total")
+
+
+# ---------------------------------------------------------------------------
+# 3. Reuse gap CDF
+# ---------------------------------------------------------------------------
+def build_reuse_gap(rows, sessions, outdir):
+    block_last_touch: dict[int, float] = {}
+    block_session: dict[int, str] = {}
+
+    # For qwenapp: no system/content split via tools. Instead:
+    #   intra-session = reuse within the same session
+    #   cross-session = reuse from a different session
+    # We don't have a system-prompt block identifier, so approximate:
+    #   system-prompt blocks = the first N blocks of each request that match
+    #   across sessions (shared system prompt prefix)
+    gaps_intra = []
+    gaps_cross = []
+
+    for r in rows:
+        ts = r["timestamp"]
+        sid = r["session_id"]
+        for h in r["hash_ids"]:
+            if h in block_last_touch:
+                gap = ts - block_last_touch[h]
+                if gap > 0:
+                    prev_sid = block_session[h]
+                    if prev_sid == sid:
+                        gaps_intra.append(gap)
+                    else:
+                        gaps_cross.append(gap)
+            block_last_touch[h] = ts
+            block_session[h] = sid
+
+    max_gap = max(
+        max(gaps_intra) if gaps_intra else 1,
+        max(gaps_cross) if gaps_cross else 1,
+    )
+    bins = np.concatenate([
+        _log_bins(0.01, max_gap * 1.01, 300),
+        [float("inf")],
+    ])
+    h_intra, _ = np.histogram(gaps_intra, bins=bins) if gaps_intra else (np.zeros(len(bins) - 1), None)
+    h_cross, _ = np.histogram(gaps_cross, bins=bins) if gaps_cross else (np.zeros(len(bins) - 1), None)
+
+    with (outdir / "f2a_reuse_gap_cdf_qwenapp.csv").open("w", newline="") as f:
+        w = csv.writer(f)
+        w.writerow(["bucket_upper_s", "intra", "cross_session"])
+        for edge, ci, cc in zip(bins[1:], h_intra, h_cross):
+            w.writerow([edge, int(ci), int(cc)])
+
+    print(f"  reuse_gap: intra={len(gaps_intra)}, cross={len(gaps_cross)}")
+
+
+# ---------------------------------------------------------------------------
+# 4. Reuse timeline
+# ---------------------------------------------------------------------------
+def build_reuse_timeline(rows, sessions, outdir, bin_s=30.0):
+    max_ts = max(r["timestamp"] for r in rows)
+    n_bins = int(math.ceil(max_ts / bin_s))
+
+    bin_data = [{"n_req": 0, "intra_tok": 0, "cross_tok": 0,
+                 "new_tok": 0, "total_tok": 0} for _ in range(n_bins)]
+
+    block_last_touch: dict[int, float] = {}
+    block_session: dict[int, str] = {}
+
+    for r in rows:
+        ts = r["timestamp"]
+        b = min(int(ts / bin_s), n_bins - 1)
+        bin_data[b]["n_req"] += 1
+
+        n_blocks = len(r["hash_ids"])
+        reuse_intra = 0
+        reuse_cross = 0
+        new_blocks = 0
+        sid = r["session_id"]
+
+        for h in r["hash_ids"]:
+            if h in block_last_touch:
+                prev_sid = block_session[h]
+                if prev_sid == sid:
+                    reuse_intra += 1
+                else:
+                    reuse_cross += 1
+            else:
+                new_blocks += 1
+            block_last_touch[h] = ts
+            block_session[h] = sid
+
+        tok_per_block = r["input_length"] / n_blocks if n_blocks > 0 else 0
+        bin_data[b]["intra_tok"] += reuse_intra * tok_per_block
+        bin_data[b]["cross_tok"] += reuse_cross * tok_per_block
+        bin_data[b]["new_tok"] += new_blocks * tok_per_block
+        bin_data[b]["total_tok"] += r["input_length"]
+
+    with (outdir / "f2a_reuse_timeline_qwenapp.csv").open("w", newline="") as f:
+        w = csv.writer(f)
+        w.writerow(["bin_start_s", "n_req", "qps", "intra_tok_s",
+                     "cross_tok_s", "new_tok_s", "reuse_tok_s", "total_tok_s", "reuse_frac"])
+        for i, d in enumerate(bin_data):
+            t0 = i * bin_s
+            qps = d["n_req"] / bin_s
+            reuse = d["intra_tok"] + d["cross_tok"]
+            total = d["total_tok"]
+            frac = reuse / total if total > 0 else 0
+            w.writerow([
+                f"{t0:.1f}", d["n_req"], f"{qps:.3f}",
+                f"{d['intra_tok']:.1f}", f"{d['cross_tok']:.1f}",
+                f"{d['new_tok']:.1f}", f"{reuse:.1f}", f"{total:.1f}",
+                f"{frac:.4f}",
+            ])
+
+    print(f"  reuse_timeline: {n_bins} bins")
+
+
+# ---------------------------------------------------------------------------
+# 5. Session skew
+# ---------------------------------------------------------------------------
+def build_session_skew(rows, sessions, outdir):
+    session_tokens = {}
+    for sid, sess_rows in sessions.items():
+        session_tokens[sid] = sum(r["input_length"] for r in sess_rows)
+
+    total = sum(session_tokens.values())
+    ranked = sorted(session_tokens.values(), reverse=True)
+    cumsum = np.cumsum(ranked) / total * 100.0
+    n = len(ranked)
+    rank_pct = np.arange(1, n + 1) / n * 100.0
+
+    step = max(1, n // 500)
+    samples = [{"rank_pct": float(rank_pct[i]), "cum_pct": float(cumsum[i])}
+               for i in range(0, n, step)]
+    if samples[-1]["rank_pct"] < 100.0:
+        samples.append({"rank_pct": 100.0, "cum_pct": 100.0})
+
+    anchors = {}
+    for pct in [1, 5, 10, 25, 50]:
+        idx = max(0, int(n * pct / 100.0) - 1)
+        anchors[f"top{pct}pct"] = float(cumsum[idx])
+
+    data = {
+        "n_sessions": n,
+        "total_tokens": int(total),
+        "anchors_check": anchors,
+        "source_trace": "trace-qwenapp-formatted/0529-1000-1200.jsonl",
+        "samples": samples,
+    }
+    (outdir / "f2b_session_skew_cdf_qwenapp.json").write_text(
+        json.dumps(data, indent=2, ensure_ascii=False))
+    print(f"  session_skew: {n} sessions, top1%={anchors['top1pct']:.1f}%")
+
+
+# ---------------------------------------------------------------------------
+def main():
+    p = argparse.ArgumentParser()
+    p.add_argument("trace", type=Path)
+    p.add_argument("--outdir", type=Path, default=None)
+    args = p.parse_args()
+
+    outdir = args.outdir or Path(__file__).resolve().parent
+    outdir.mkdir(parents=True, exist_ok=True)
+
+    print(f"Reading {args.trace} ...")
+    rows = []
+    sessions: dict[str, list[dict]] = defaultdict(list)
+    with args.trace.open() as f:
+        for line in f:
+            r = json.loads(line)
+            rows.append(r)
+            sessions[r["session_id"]].append(r)
+
+    for sess_rows in sessions.values():
+        sess_rows.sort(key=lambda r: r["turn"])
+
+    print(f"Loaded {len(rows)} rows, {len(sessions)} sessions")
+
+    build_token_length(rows, sessions, outdir)
+    build_block_lifetime(rows, outdir)
+    build_reuse_gap(rows, sessions, outdir)
+    build_reuse_timeline(rows, sessions, outdir)
+    build_session_skew(rows, sessions, outdir)
+
+    print("Done.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/format_qwenapp_trace.py b/format_qwenapp_trace.py
new file mode 100644
index 0000000..a33b00c
--- /dev/null
+++ b/format_qwenapp_trace.py
@@ -0,0 +1,269 @@
+"""Format raw qwenapp Aquila access-log trace into the standard release format.
+
+qwenapp traces differ from coding-agent traces in two ways:
+  1. Prompts are pre-rendered strings (not a `messages` list).
+  2. No tool calls — this is a general chat + RAG workload.
+
+The script reads the raw JSONL, groups rows by session_id, assigns
+chat_id / parent_chat_id / turn, tokenizes the prompt to build hash_ids
+for prefix-cache replay, and writes the standard format understood by
+the replayer and sample_trace.py.
+
+Usage (run on dash1 where the trace and tokenizer live):
+    python scripts/format_qwenapp_trace.py \
+        --input ~/ali-trace/trace-qwenapp/0529-1000-1200.jsonl \
+        --output ~/ali-trace/trace-qwenapp-formatted/0529-1000-1200.jsonl \
+        --tokenizer ~/ali-trace/trace_model_meta/Qwen/Qwen3-Coder-480B-A35B-Instruct/tokenizer.json
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import sys
+from array import array
+from collections import defaultdict
+from pathlib import Path
+
+from tokenizers import Tokenizer
+from tqdm.auto import tqdm
+
+
+def _block_digest(block: list[int]) -> bytes:
+    digest = hashlib.blake2b(digest_size=16)
+    digest.update(len(block).to_bytes(4, "little", signed=False))
+    digest.update(array("I", block).tobytes())
+    return digest.digest()
+
+
+def _parse_row(raw: dict) -> dict | None:
+    """Extract the fields we need from one raw Aquila access-log row."""
+    request_params = raw.get("request_params")
+    if isinstance(request_params, str):
+        request_params = json.loads(request_params)
+    response_params = raw.get("response_params")
+    if isinstance(response_params, str):
+        response_params = json.loads(response_params)
+
+    if not isinstance(request_params, dict) or not isinstance(response_params, dict):
+        return None
+
+    inp = request_params.get("payload", {}).get("input", {})
+    resp_payload = response_params.get("payload", {})
+    usage = resp_payload.get("usage", {})
+    resp_header = response_params.get("header", {})
+
+    prompt = inp.get("prompt", "")
+    if not prompt:
+        return None
+
+    input_tokens = int(usage.get("input_tokens", 0))
+    output_tokens = int(usage.get("output_tokens", 0))
+    if input_tokens <= 0:
+        return None
+
+    status_code = str(resp_header.get("status_code", raw.get("status_code", "")))
+    if status_code and status_code != "200" and status_code != "1000":
+        return None
+
+    # Timestamp: prefer the `time` field (ms precision), fall back to __time__
+    time_str = raw.get("time", "")
+    wall_clock_s = _parse_time_to_epoch(time_str) if time_str else None
+    if wall_clock_s is None:
+        epoch_s = raw.get("__time__")
+        wall_clock_s = float(epoch_s) if epoch_s else 0.0
+
+    session_id = str(inp.get("session_id", ""))
+    user_id = str((inp.get("source") or {}).get("user_id", ""))
+    request_id = str(raw.get("request_id", inp.get("request_id", "")))
+
+    return {
+        "session_id": session_id,
+        "user_id": user_id,
+        "request_id": request_id,
+        "wall_clock_s": wall_clock_s,
+        "input_tokens": input_tokens,
+        "output_tokens": output_tokens,
+        "prompt": prompt,
+        "total_cost_time_ms": int(raw.get("total_cost_time", 0)),
+    }
+
+
+def _parse_time_to_epoch(time_str: str) -> float | None:
+    """Parse '2026-05-29 09:59:59.505' (UTC+8) to epoch seconds."""
+    import datetime
+    try:
+        dt = datetime.datetime.strptime(time_str.strip(), "%Y-%m-%d %H:%M:%S.%f")
+    except ValueError:
+        try:
+            dt = datetime.datetime.strptime(time_str.strip(), "%Y-%m-%d %H:%M:%S")
+        except ValueError:
+            return None
+    tz_utc8 = datetime.timezone(datetime.timedelta(hours=8))
+    dt = dt.replace(tzinfo=tz_utc8)
+    return dt.timestamp()
+
+
+def main() -> None:
+    p = argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawDescriptionHelpFormatter)
+    p.add_argument("--input", type=Path, required=True)
+    p.add_argument("--output", type=Path, required=True)
+    p.add_argument("--tokenizer", type=Path, required=True,
+                   help="Path to tokenizer.json (Qwen3 family)")
+    p.add_argument("--block-size", type=int, default=512)
+    p.add_argument("--ttp", action="store_true",
+                   help="Also emit time_to_parent_chat field")
+    args = p.parse_args()
+
+    print(f"Loading tokenizer from {args.tokenizer} ...")
+    tokenizer = Tokenizer.from_file(str(args.tokenizer))
+
+    # --- Pass 1: read all rows, group by session ---
+    print(f"Reading {args.input} ...")
+    sessions: dict[str, list[dict]] = defaultdict(list)
+    skipped = 0
+    total = 0
+
+    with args.input.open("r", encoding="utf-8") as fh:
+        for line in tqdm(fh, desc="Read", unit="row"):
+            total += 1
+            raw = json.loads(line)
+            parsed = _parse_row(raw)
+            if parsed is None:
+                skipped += 1
+                continue
+            sessions[parsed["session_id"]].append(parsed)
+
+    print(f"  Read {total} rows, kept {total - skipped}, skipped {skipped}")
+    print(f"  Sessions: {len(sessions)}")
+
+    # --- Sort each session by timestamp and assign turn/chat_id ---
+    for rows in sessions.values():
+        rows.sort(key=lambda r: r["wall_clock_s"])
+
+    # Assign globally-unique chat_ids
+    next_chat_id = 1
+    chat_id_map: dict[str, int] = {}  # request_id -> chat_id
+    all_rows: list[dict] = []
+
+    for sid, rows in sessions.items():
+        prev_chat_id = -1
+        prev_end_s = None
+        for turn_idx, row in enumerate(rows, start=1):
+            cid = next_chat_id
+            next_chat_id += 1
+            chat_id_map[row["request_id"]] = cid
+
+            ttp = None
+            if prev_end_s is not None:
+                ttp = row["wall_clock_s"] - prev_end_s
+                if ttp < 0:
+                    ttp = 0.0
+
+            all_rows.append({
+                "session_id": sid,
+                "chat_id": cid,
+                "parent_chat_id": prev_chat_id,
+                "turn": turn_idx,
+                "wall_clock_s": row["wall_clock_s"],
+                "input_tokens": row["input_tokens"],
+                "output_tokens": row["output_tokens"],
+                "prompt": row["prompt"],
+                "time_to_parent_chat": ttp,
+                "total_cost_time_ms": row["total_cost_time_ms"],
+            })
+
+            prev_chat_id = cid
+            prev_end_s = row["wall_clock_s"] + row["total_cost_time_ms"] / 1000.0
+
+    # Sort globally by timestamp
+    all_rows.sort(key=lambda r: r["wall_clock_s"])
+
+    # Re-zero timestamps
+    t0 = all_rows[0]["wall_clock_s"] if all_rows else 0.0
+
+    # --- Pass 2: tokenize and hash ---
+    print(f"Tokenizing {len(all_rows)} rows (block_size={args.block_size}) ...")
+    block_ids_by_digest: dict[str, int] = {}
+    next_block_id = 0
+
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    with args.output.open("w", encoding="utf-8") as out:
+        for row in tqdm(all_rows, desc="Tokenize+write", unit="row"):
+            token_ids = tokenizer.encode(row["prompt"]).ids
+
+            # Build block hashes
+            hash_ids = []
+            for i in range(0, len(token_ids), args.block_size):
+                block = token_ids[i:i + args.block_size]
+                digest_hex = _block_digest(block).hex()
+                bid = block_ids_by_digest.get(digest_hex)
+                if bid is None:
+                    bid = next_block_id
+                    next_block_id += 1
+                    block_ids_by_digest[digest_hex] = bid
+                hash_ids.append(bid)
+
+            # Pad hash_ids to match provider input_tokens if needed
+            target_blocks = (row["input_tokens"] + args.block_size - 1) // args.block_size
+            while len(hash_ids) < target_blocks:
+                pad_key = f"pad:{row['chat_id']}:{len(hash_ids)}"
+                pad_hex = hashlib.blake2b(pad_key.encode(), digest_size=16).hexdigest()
+                bid = block_ids_by_digest.get(pad_hex)
+                if bid is None:
+                    bid = next_block_id
+                    next_block_id += 1
+                    block_ids_by_digest[pad_hex] = bid
+                hash_ids.append(bid)
+            hash_ids = hash_ids[:target_blocks]
+
+            out_row = {
+                "chat_id": row["chat_id"],
+                "parent_chat_id": row["parent_chat_id"],
+                "timestamp": round(row["wall_clock_s"] - t0, 3),
+                "input_length": row["input_tokens"],
+                "output_length": row["output_tokens"],
+                "type": "app",
+                "turn": row["turn"],
+                "hash_ids": hash_ids,
+                "session_id": row["session_id"],
+            }
+            if args.ttp and row["time_to_parent_chat"] is not None:
+                out_row["time_to_parent_chat"] = round(row["time_to_parent_chat"], 3)
+
+            out.write(json.dumps(out_row, ensure_ascii=False) + "\n")
+
+    # --- Summary ---
+    n_sessions = len(sessions)
+    n_requests = len(all_rows)
+    turns_per = [len(v) for v in sessions.values()]
+    multi = sum(1 for t in turns_per if t > 1)
+    span_s = all_rows[-1]["wall_clock_s"] - t0 if all_rows else 0
+    qps = n_requests / span_s if span_s > 0 else 0
+    inp_lens = [r["input_tokens"] for r in all_rows]
+    out_lens = [r["output_tokens"] for r in all_rows]
+
+    unique_blocks = next_block_id
+    shared_blocks = sum(1 for d, bid in block_ids_by_digest.items()
+                        if not d.startswith("pad:"))  # rough
+    # Count actually shared (appear in >1 request)
+    block_ref_count: dict[int, int] = defaultdict(int)
+    # Re-scan is expensive; use the hash_ids we just wrote
+    # Instead, just report unique count
+    print(f"\n=== Summary ===")
+    print(f"  Requests: {n_requests}")
+    print(f"  Sessions: {n_sessions}")
+    print(f"  Multi-turn: {multi} ({multi*100/n_sessions:.1f}%)")
+    print(f"  Turns/session: min={min(turns_per)} max={max(turns_per)} avg={sum(turns_per)/len(turns_per):.1f}")
+    print(f"  Input tokens: min={min(inp_lens)} max={max(inp_lens)} avg={sum(inp_lens)/len(inp_lens):.0f}")
+    print(f"  Output tokens: min={min(out_lens)} max={max(out_lens)} avg={sum(out_lens)/len(out_lens):.0f}")
+    print(f"  Trace span: {span_s:.1f}s ({span_s/60:.1f} min)")
+    print(f"  QPS: {qps:.2f}")
+    print(f"  Unique hash blocks: {unique_blocks}")
+    print(f"\nWrote {n_requests} rows to {args.output}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/session_sticky_imbalance.py b/session_sticky_imbalance.py
new file mode 100644
index 0000000..017d093
--- /dev/null
+++ b/session_sticky_imbalance.py
@@ -0,0 +1,144 @@
+#!/usr/bin/env python3
+"""Session-sticky routing imbalance simulation over a formatted ali-trace.
+
+Routing policy:
+  - turn==1 (parent_chat_id == -1): dispatch to the instance with the FEWEST
+    total requests received in the trailing LOAD_WIN seconds (least-loaded).
+  - turn>1: pinned to the same instance as the request's parent turn (sticky).
+
+For every tumbling WINDOW-second time window we count the total requests landed
+on each of the N instances and report imbalance = max / mean over all N
+instances (idle instances counted in the mean). A CDF of the per-window
+imbalance values is plotted.
+
+Run on dash0 (trace lives there). Trace is timestamp-sorted and every parent
+precedes its child, so a single in-order pass suffices.
+"""
+import argparse
+import json
+from collections import deque
+
+import numpy as np
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+
+def simulate(trace, n_inst, window, load_win, weight="req", route_by="req"):
+    # `route_by` is the trailing-load_win quantity the turn-1 least-loaded
+    # routing balances on; `weight` is what each window accumulates for the
+    # imbalance statistic. Both take values:
+    #   "req"    -> 1 per request (request count)
+    #   "tokens" -> input_length + output_length per request (total tokens)
+    assign = {}                       # chat_id -> instance
+    cnt = np.zeros(n_inst, dtype=np.int64)   # trailing-load_win route load per inst
+    dq = deque()                      # (ts, inst, rw) pending expiry from the window
+    wincounts = {}                    # window_idx -> np.array(n_inst)
+
+    with open(trace) as fh:
+        for line in fh:
+            d = json.loads(line)
+            t = d["timestamp"]
+            tok = d["input_length"] + d["output_length"]
+
+            # expire dispatches older than the trailing load window
+            thr = t - load_win
+            while dq and dq[0][0] < thr:
+                _, old, ow = dq.popleft()
+                cnt[old] -= ow
+
+            if d["parent_chat_id"] == -1:          # turn 1: least-loaded
+                inst = int(cnt.argmin())
+            else:                                  # turn>1: sticky to parent
+                inst = assign[d["parent_chat_id"]]
+            assign[d["chat_id"]] = inst            # chain so deeper turns follow
+
+            rw = (tok if route_by == "tokens"
+                  else d["input_length"] if route_by == "input"
+                  else 1)
+            dq.append((t, inst, rw))
+            cnt[inst] += rw
+
+            w = int(t // window)
+            arr = wincounts.get(w)
+            if arr is None:
+                arr = np.zeros(n_inst, dtype=np.int64)
+                wincounts[w] = arr
+            arr[inst] += tok if weight == "tokens" else 1
+
+    ws = sorted(wincounts)
+    imb = np.empty(len(ws))
+    for i, w in enumerate(ws):
+        arr = wincounts[w]
+        mean = arr.sum() / n_inst
+        imb[i] = arr.max() / mean if mean > 0 else 0.0
+    return imb, [wincounts[w] for w in ws]
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--trace", required=True)
+    ap.add_argument("--n-inst", type=int, default=380)
+    ap.add_argument("--window", type=float, default=60.0)
+    ap.add_argument("--load-win", type=float, default=60.0)
+    ap.add_argument("--weight", choices=["req", "tokens"], default="req",
+                    help="per-window quantity: request count or total tokens")
+    ap.add_argument("--route-by", choices=["req", "tokens", "input"],
+                    default="req",
+                    help="turn-1 least-loaded routing balances on this quantity "
+                         "(input = arrival-known input_length only)")
+    ap.add_argument("--out-prefix", required=True)
+    args = ap.parse_args()
+
+    imb, _ = simulate(args.trace, args.n_inst, args.window, args.load_win,
+                      args.weight, args.route_by)
+
+    pcts = {p: float(np.percentile(imb, p)) for p in (50, 90, 95, 99)}
+    stats = {
+        "trace": args.trace,
+        "n_inst": args.n_inst,
+        "window_s": args.window,
+        "load_win_s": args.load_win,
+        "weight": args.weight,
+        "route_by": args.route_by,
+        "n_windows": len(imb),
+        "imbalance_min": float(imb.min()),
+        "imbalance_mean": float(imb.mean()),
+        "imbalance_p50": pcts[50],
+        "imbalance_p90": pcts[90],
+        "imbalance_p95": pcts[95],
+        "imbalance_p99": pcts[99],
+        "imbalance_max": float(imb.max()),
+    }
+    with open(args.out_prefix + "_stats.json", "w") as f:
+        json.dump(stats, f, indent=2)
+    np.save(args.out_prefix + "_imbalance.npy", imb)
+
+    # CDF
+    x = np.sort(imb)
+    y = np.arange(1, len(x) + 1) / len(x)
+    fig, ax = plt.subplots(figsize=(6, 4))
+    ax.plot(x, y, marker=".", ms=4, lw=1.4, color="#1f77b4")
+    ax.axvline(1.0, color="gray", ls="--", lw=1, label="perfect balance (1.0)")
+    ax.axvline(pcts[50], color="#2ca02c", ls=":", lw=1,
+               label=f"p50={pcts[50]:.2f}")
+    ax.axvline(pcts[90], color="#d62728", ls=":", lw=1,
+               label=f"p90={pcts[90]:.2f}")
+    qty = "total tokens" if args.weight == "tokens" else "num_req"
+    ax.set_xlabel("per-window %s imbalance (max / mean over %d instances)"
+                  % (qty, args.n_inst))
+    ax.set_ylabel("CDF over %d one-minute windows" % len(imb))
+    routeq = {"tokens": "tokens", "input": "input-only"}.get(
+        args.route_by, "num_req")
+    ax.set_title("Session-sticky imbalance: measure=%s, route-by=%s"
+                 % (qty, routeq))
+    ax.grid(True, alpha=0.3)
+    ax.legend(fontsize=8)
+    fig.tight_layout()
+    fig.savefig(args.out_prefix + "_cdf.png", dpi=140)
+    print(json.dumps(stats, indent=2))
+
+
+if __name__ == "__main__":
+    main()