IceCube Master Thesis Preparation Project

Narrowing the Gap Between Simulation and Real Data in IceCube

Niels Bohr Institute. Repository: IceCube-Master-Thesis-Preparation-Project.

This project compares simulated Muon Gun atmospheric muons with real 2021 burnsample atmospheric muons in IceCube. The repository contains the project story, the project figures, and the copied analysis code needed to trace each result.

The central question is:

If an IceCube machine-learning model is trained on simulation and applied to real data, how visible is the remaining simulation-to-data mismatch?

The answer is that the mismatch is very visible. A pulse-level transformer can separate real 2021 burnsample atmospheric muons from simulated Muon Gun atmospheric muons almost perfectly. The project then applies a sequence of targeted corrections. They reduce the mismatch, especially for stopped muons, but the final samples are still clearly distinguishable.

IceCube instruments Antarctic ice with strings of DOMs. The analysis compares real and simulated atmospheric muons after both samples have been reconstructed into pulse-level detector readout.	Through-going and stopped muons leave visibly different pulse patterns, so the MC-vs-data comparison is performed separately for the two event classes.

Pulse-Level Input

The models operate on SplitInIcePulses: each event is a set of pulse rows, and each pulse row carries detector position, timing, charge, DOM efficiency, and HLC/SLC status.

Variable	Example from 2021 burnsample	Meaning
charge	1.275 PE	Reconstructed pulse charge.
dom_time	6530 ns	Reconstructed pulse time.
dom_x, dom_y, dom_z	-256.14 m, -521.08 m, 121.57 m	DOM position in IceCube detector coordinates.
hlc	0	Pulse-type flag: 1 for HLC, 0 for SLC.
rde	1	Relative DOM efficiency.
event_time	59337.887082 MJD	Event-level time label.
event_no	0	Internal event index tying pulses to the same event.
string, dom_number	1, 23	DOM identifier.

This pulse-format example is the data representation used throughout the project. The code that builds and loads the pulse-level rows lives in analysis/MC_vs_BS_analysis/scripts and mc_vs_data_parquet_dataset.py.

Read together, these two figures are the shortest version of the project. The ROC curves show that MC and data remain rank-separable after all corrections. The logit distributions show that the corrections still matter: the classifier loses much of its confidence, especially for stopped muons.

Main Result

The main diagnostic is a binary MC-vs-data classifier. An AUC of 0.5 would mean that data and simulation are indistinguishable to the benchmark model. The table shows the project in one line: each correction helps, but none closes the gap.

Stage	What changed	Stopped AUC	Through-going AUC
Baseline	No correction	0.9882	0.9960
Angular GB reweighting	MC direction distribution reweighted to data	0.9688	0.9935
Pulse merging	Low-charge satellite pulses merged on each DOM	0.9583	0.9892
HLC re-labelling	Most HLC-like simulated SLC pulses flipped to HLC	0.9281	0.9851
Low-kappa removal	Low-confidence vMF direction events removed	0.9218	0.9848

The largest single improvement comes from the data-driven HLC re-labelling. The final ROC is still far from chance, but the corresponding logit distributions show that the classifier becomes much less confident after the corrections. In other words: the gap is narrowed, not solved.

How To Read This Repository

• Start with this README for the story and the central figures.
• Open docs/project_summary.md for a shorter project walkthrough.
• Open docs/figure_index.md when you want figure to code provenance in table form.
• Open docs/project_traceability.md to audit which project sections, figures, tables, and numerical results are backed by which code or source.
• Open docs/analysis_runbook.md when you want the practical run order: stage purpose, main scripts, required inputs, and generated outputs.
• Open docs/code_map.md to navigate the analysis tree by scientific task.
• Open docs/reproduction_notes.md for what is included, what is deliberately excluded, and what would be needed to rerun the analysis.
• Browse analysis/ for the copied analysis source. Raw data, generated parquet files, SQLite databases, checkpoints, and logs are not included.

Scientific Story

IceCube analyses often train models on Monte Carlo simulation because only simulation provides truth labels. That workflow is only reliable if simulation and real detector data are sufficiently aligned in the variables seen by the model. This project stress-tests that assumption with atmospheric muons: abundant, track-like events that exercise the same detector, ice, readout, and reconstruction chain as signal-like muons.

The analysis proceeds as follows:

1. Define the detector data representation: SplitInIcePulses, where each event is a variable-size set of pulse tokens with charge, time, DOM position, HLC/SLC flag, and relative DOM efficiency.
2. Separate atmospheric muons into stopped and through-going classes using a transformer trained on MC truth, then apply the same classifier to MC and data.
3. Train a pulse-level transformer to distinguish MC from data. This becomes the benchmark for simulation-to-data mismatch.
4. Apply four corrections or diagnostics in sequence: angular GB reweighting, pulse merging, HLC re-labelling, and a vMF uncertainty cut.
5. Compare the full correction chain with ROC curves and logit distributions.

Key Takeaways

Finding	Evidence	Where to look
MC and data are strongly separable before correction.	Baseline MC-vs-data AUC is 0.9882 for stopped and 0.9960 for through-going events.	Baseline comparison
Direction matters, but is not the whole problem.	GB reweighting in reconstructed zenith/azimuth improves the benchmark but leaves high AUCs.	Angular reweighting
Low-charge DOM behavior is an important handle.	Pulse merging reduces small-pulse mismatches and modestly lowers MC-vs-data separability.	Pulse merging
HLC/SLC modelling carries major residual information.	HLC re-labelling gives the largest single AUC improvement.	HLC re-labelling
Learned direction uncertainty exposes a problematic MC population.	Low-kappa MC events show a pole-collapse signature; removing them gives the final benchmark row.	vMF uncertainty

Repository Layout

.
├── analysis/
│   ├── ThroughOrStopped_muon/
│   ├── MC_vs_BS_analysis/
│   │   ├── scripts/
│   │   ├── zenith_azimuth_inference/
│   │   └── GBreweighting/
│   └── Classifiers/
├── docs/
│   ├── project_summary.md
│   ├── figure_index.md
│   ├── project_traceability.md
│   ├── analysis_runbook.md
│   ├── code_map.md
│   └── reproduction_notes.md
└── figures/
    ├── report/
    └── report_previews/

The tracked analysis tree contains Python source, Slurm scripts, notebooks with outputs cleared, configs, model metrics, and small text summaries. It does not contain raw detector data, generated parquet tables, SQLite databases, CSV exports, NumPy arrays, pickle files, model checkpoints, logs, or the local full write-up source or build products.

Project Walkthrough

The analysis starts with the detector readout and pulse representation, then builds the stopped/through-going split, the MC-vs-data benchmark, the corrections, and the final staged comparison.

1. Introduction: why data-MC agreement matters

Most IceCube machine-learning results depend on simulation. If real detector data and simulated events differ in the pulse-level variables used by the model, downstream reconstruction and classification can inherit simulation artifacts. Atmospheric muons are used here as a high-statistics test case because they are abundant, track-like, and recorded by the same detector.

The implementation entry points are grouped by scientific task in code_map.md.

2. Detector, DOM readout, and pulse-level events

IceCube records Cherenkov light in DOMs embedded in Antarctic ice. The raw readout is reduced to pulse-level variables in SplitInIcePulses, and the later models operate on those pulse tokens.

The key point is that the analysis is not comparing abstract event labels. It is comparing distributions of pulse charge, time, DOM position, HLC/SLC status, and DOM efficiency.

Detector and readout figures

Figure	What it shows	Code/source
icecube.png	Places the analysis in the IceCube detector geometry: strings, DOMs, DeepCore, and IceTop.	External/reference detector schematic.
shrenkov.pdf	Explains why charged particles crossing the ice produce the light recorded by DOMs.	Schematic asset.
plot_run126491_event30343391_DOM83-31-0.pdf	Shows the detailed ATWD/fADC waveform behind one real HLC hit, including later small pulses.	waveform_demo
icecube_events.png	Shows the pulse-level appearance of track, cascade, and double-bang event types.	External/reference event-type figure.
event_display_through_stopped.pdf	Motivates why stopped and through-going atmospheric muons are treated as separate classes.	plot_event_display_through_stopped.py, pulse_event_display.py

3. Machine-learning tools used later

Two families of tools appear repeatedly. Gradient-boosted decision trees are used for reweighting MC in reconstructed direction space. Transformer models are used for the main pulse-level tasks: stopped/through classification, direction reconstruction, MC-vs-data benchmarking, and HLC re-labelling.

The decision-tree and BDT schematics define the reweighting method used for the angular correction. The transformer schematic defines the pulse-token model used for the stopped/through split, direction reconstruction, MC-vs-data benchmarking, and HLC re-labelling.

Machine-learning schematic figures

Figure	What it shows	Code/source
decision_tree.pdf	Introduces the threshold-cut logic behind tree models.	make_ch3_figures.py
bdt_schematic.pdf	Explains why many shallow trees can form a stronger boosted model.	make_ch3_figures.py
transformer_architecture.pdf	Shows the event model used repeatedly: pulse tokens, embeddings, transformer blocks, pooling, and prediction head.	make_ch3_figures.py

4. Stopped/through-going classifier

Stopped and through-going muons have different light patterns, so the project does not compare them as one mixed sample. A transformer is trained on MC truth labels to split the simulated events, and then the same model is applied to both MC and real burnsample data.

Every later MC-vs-data comparison is performed separately for the two classifier-defined classes.

This split is implemented by training the classifier, running it over both samples, and then generating the diagnostic plots from the stored predictions:

• train_stopped_transformer.py
• run_inference.py
• plot_stopped_transformer_documentation.py

Stopped/through classifier figures

Figure	What it shows	Code/source
training_history.pdf	Verifies the training behavior and selected best epoch for the stopped/through model.	train_stopped_transformer.py, plot_stopped_transformer_documentation.py
test_performance.pdf	Shows the held-out MC classification performance that justifies using the split downstream.	plot_stopped_transformer_documentation.py
mc_test_score_distributions.pdf	Shows how the model output separates the two MC truth classes.	plot_stopped_transformer_documentation.py

5. Baseline MC-vs-data comparison

The baseline comparison asks what differs before any correction. Pulse-level and event-level histograms show visible discrepancies: low-charge behavior, time tails, depth structure, high-charge tails, and HLC fraction. The stronger test is the MC-vs-data transformer, which confirms that the joint feature-space mismatch is large.

The baseline comparison is produced by the MC-vs-data transformer and the distribution plotting scripts:

• transformer model code
• train_mcdata_parquet.py
• make_pulse_level_a4_figure.py
• make_event_aggregate_a4_figure.py

Baseline distribution figures

Figure	What it shows	Code/source
pulse_level_variables_unmerged_full_page1.pdf	Shows baseline dom_time, charge, dom_x, and dom_y disagreements by event class.	make_pulse_level_a4_figure.py
pulse_level_variables_unmerged_full_page2.pdf	Shows baseline dom_z, rde, and hlc behavior, including the important HLC mismatch.	make_pulse_level_a4_figure.py
event_level_aggregates_unmerged_full_page1.pdf	Compares event size and charge summaries before corrections.	make_event_aggregate_a4_figure.py
event_level_aggregates_unmerged_full_page2.pdf	Compares event time and depth summaries before corrections.	make_event_aggregate_a4_figure.py
event_level_aggregates_unmerged_full_page3.pdf	Shows depth spread and HLC fraction at event level, which motivates later HLC work.	make_event_aggregate_a4_figure.py

6. Direction reconstruction and angular GB reweighting

The first correction tests whether MC and data disagree partly because they enter the detector from different directions. A transformer reconstructs zenith and azimuth. A gradient-boosted reweighter then changes the MC weights in reconstructed (zenith, azimuth) space so that the angular distribution matches data more closely.

This correction improves the benchmark and fixes the targeted angular distributions, but much of the MC-vs-data separability remains.

This correction is produced in three linked steps: train/infer the direction model, fit the angular GB reweighter, and then regenerate the weighted comparison plots:

• direction_transformer_hlc_rde_unmerged_2M
• plot_direction_transformer_documentation.py
• fit_GBreweighter_hlc_rde_unmerged_2M.py

Direction and GB-reweighting figures

Figure	What it shows	Code/source
open_angle_performance.pdf	Shows the direction reconstruction quality on held-out MC.	plot_direction_transformer_documentation.py
mc_data_zenith_azimuth_overlay.pdf	Shows the reconstructed direction mismatch before angular reweighting.	plot_direction_transformer_documentation.py
mc_data_zenith_azimuth_overlay_with_GBR.pdf	Confirms that the GB reweighter fixes the direction distribution it was trained to fix.	fit_GBreweighter_hlc_rde_unmerged_2M.py, plot_direction_transformer_documentation.py
pulse_level_variables_unmerged_gbweighted_full_page1.pdf	Tests whether angular weights also improve pulse-level variables such as charge and horizontal position.	make_pulse_level_a4_figure.py
pulse_level_variables_unmerged_gbweighted_full_page2.pdf	Tests whether angular weights improve vertical position, DOM efficiency, and HLC behavior.	make_pulse_level_a4_figure.py
event_level_aggregates_unmerged_gbweighted_full_page1.pdf	Checks event size and charge summaries after angular weighting.	make_event_aggregate_a4_figure.py
event_level_aggregates_unmerged_gbweighted_full_page2.pdf	Checks time and depth summaries after angular weighting.	make_event_aggregate_a4_figure.py
event_level_aggregates_unmerged_gbweighted_full_page3.pdf	Shows that HLC fraction remains a strong residual mismatch after angular weighting.	make_event_aggregate_a4_figure.py

7. Pulse merging

The angular correction does not remove the low-charge excess in data. The next step targets a pulse-splitting effect: small HLC pulses below 0.3 PE are merged into the nearest above-threshold pulse on the same DOM using the PulseMerger algorithm.

This step is physically motivated by waveform behavior and improves the benchmark, but only modestly.

The pulse-merging correction is implemented by the merger itself and validated with both aggregate plots and a single-DOM illustration:

• pulse_merger.py
• plot_pulse_merging.py
• make_small_pulse_merge_plot.py

Pulse-merging figures

Figure	What it shows	Code/source
small_pulses_through_run136141_event242722_string61_dom3_final_legend_default_up.pdf	Demonstrates the pulse-merging rule on one DOM signal: small pulses are folded into a nearby larger pulse while charge is preserved.	make_small_pulse_merge_plot.py, pulse_merger.py
mc_data_charge_hlc_slc.pdf	Shows how charge distributions differ for HLC and SLC pulses and why low-charge handling matters.	plot_pulse_merging.py
pulses_per_dom.pdf	Checks whether MC and data differ in multi-pulse DOM behavior.	plot_pulse_merging.py

8. Feature importance and HLC re-labelling

After pulse merging, permutation feature importance points to the hlc flag as the strongest remaining pulse-level carrier of MC-vs-data separation. The project therefore trains an HLC/SLC classifier on real data and applies it to MC. The most HLC-like simulated SLC pulses are flipped to HLC until the event-level HLC-fraction distribution best matches data.

This is the most important correction in the project.

The HLC study first measures which feature the classifier uses, then sweeps SLC-to-HLC flip rates, applies the selected flip, and replots the HLC-fraction agreement:

• eval_transformer_perm_compare_hlcflip.py
• run_hlc_flip_sweep_merged_v2_all.py
• apply_transformer_hlc_best_flip_parquets.py
• plot_hlc_flip_sweep_merged_v2_side_by_side.py

HLC feature and re-labelling figures

Figure	What it shows	Code/source
hlc_flip_rate_sweep_merged_v2_stopped_through_side_by_side_0_to_10p0_step0p5.pdf	Chooses the SLC-to-HLC flip rate by minimizing the MC/data HLC-fraction Wasserstein distance.	run_hlc_flip_sweep_merged_v2_all.py, plot_hlc_flip_sweep_merged_v2_side_by_side.py
hlc_frac_mc_vs_data_merged_v2_stopped_through_best_transformer_flip_side_by_side.pdf	Shows the HLC-fraction distribution after the selected transformer-based flip.	plot_hlc_frac_merged_v2_best_transformer_flip.py

9. Charge-time and afterpulse search

The waveform example earlier shows delayed small pulses after a main signal. The charge-time analysis checks whether an afterpulse-like structure appears in the charge vs dom_time plane after the main corrections. The search does not find a clean isolated delayed low-charge island in the pulse-level representation.

The charge-time check is produced by the afterpulse plotting scripts, with the earlier waveform example kept as the DOM-level motivation:

• plot_afterpulse_a4_transformer_hlcflip_best.py
• plot_afterpulse_master.py
• waveform_demo

Charge-time and afterpulse figures

Figure	What it shows	Code/source
afterpulse_stopped_mc_transformer_hlcflip_best.pdf	Charge-time plane for stopped MC after the main corrections.	plot_afterpulse_a4_transformer_hlcflip_best.py
afterpulse_stopped_data_transformer_hlcflip_best.pdf	Charge-time plane for stopped data after the main corrections.	plot_afterpulse_a4_transformer_hlcflip_best.py
afterpulse_stopped_mc_over_data_transformer_hlcflip_best.pdf	Residual view for stopped events, used to look for localized MC/data structures.	plot_afterpulse_a4_transformer_hlcflip_best.py
afterpulse_through_mc_transformer_hlcflip_best.pdf	Charge-time plane for through-going MC after the main corrections.	plot_afterpulse_a4_transformer_hlcflip_best.py
afterpulse_through_data_transformer_hlcflip_best.pdf	Charge-time plane for through-going data after the main corrections.	plot_afterpulse_a4_transformer_hlcflip_best.py
afterpulse_through_mc_over_data_transformer_hlcflip_best.pdf	Residual view for through-going events, used to look for localized delayed-pulse structure.	plot_afterpulse_a4_transformer_hlcflip_best.py

10. vMF uncertainty and final diagnostic

The final direction model predicts both a direction and a von Mises-Fisher concentration parameter kappa. Large kappa means the model is confident; small kappa means the event is diffuse or ambiguous. Low-kappa MC events show a pole-collapse behavior in the direction prediction, so events with kappa < 10 are removed as the final diagnostic cut.

This improves the final benchmark slightly and identifies a concrete problematic MC population.

The vMF diagnostic is produced by the uncertainty-aware direction model and the scripts that plot kappa, diagnose low-kappa MC events, and build the final cut:

• direction_transformer_vmf_final_hlcflip
• vmf_code_bundle
• plot_vmf_uncertainty_final_hlcflip.py
• plot_vmf_pole_collapse_evidence.py

vMF uncertainty figures

Figure	What it shows	Code/source
vmf_sphere.pdf	Explains the meaning of the vMF concentration parameter kappa.	Schematic asset used by the vMF section.
vmf_training_history_loss_opening_kappa.pdf	Shows vMF direction-model training, opening angle, and predicted uncertainty behavior.	train_vmf_final_hlcflip.py, plot_vmf_uncertainty_final_hlcflip.py
vmf_kappa_mc_data_stopped_through_side_by_side.pdf	Compares predicted event confidence in MC and data for stopped and through-going samples.	plot_vmf_uncertainty_final_hlcflip.py
vmf_pole_collapse_evidence.pdf	Diagnoses the low-kappa MC population where predictions collapse toward the vertical.	plot_vmf_pole_collapse_evidence.py, diagnose_low_kappa_mc.py

11. Final benchmark and interpretation

The final section collects the staged MC-vs-data benchmark. The ROC curves show that the classifier can still separate the samples after all corrections. The logit distributions show a softer but important point: the corrections remove much of the classifier confidence, especially for stopped muons.

The remaining mismatch is likely not one single effect. Plausible handles are residual time-distribution issues, horizontal coordinate or surface-entry differences, multi-pulse DOM behavior, HLC/SLC modeling, afterpulse effects that are difficult to isolate after pulse extraction, and possible imperfections in the muon selection or stopped/through split.

The final comparison is built by collecting staged logits and plotting the ROC and raw-logit views of the same benchmark:

• build_stage_test_logits.py
• plot_stage_logit_roc_overlay.py
• plot_stage_logit_catalog.py

Final benchmark figures

Figure	What it shows	Code/source
five_stage_logit_roc_overlay_combined.pdf	The central project result: staged MC-vs-data separability after each correction.	plot_stage_logit_roc_overlay.py
logit_catalog_common_xlim.pdf	Shows how the classifier confidence changes across the correction chain.	plot_stage_logit_catalog.py

Open Questions

The final stage does not mean the remaining mismatch has been explained away. It means several specific handles have been reduced and the benchmark classifier has become less confident. The remaining threads are concrete:

• dom_time remains important in the permutation study, and the late-time data tail is not directly corrected here.
• dom_x and dom_y still carry separation, possibly connected to different surface-entry or detector-entry distributions for real and simulated muons.
• Data show more multi-pulse DOM behavior than MC, but this project only uses pulse merging as a partial correction.
• The charge-time search does not isolate a clean afterpulse population, but a more targeted sample or higher-energy selection could be more sensitive.
• The real-data muon selection comes from an external GNN classifier, and the stopped/through split is itself MC-trained. Both selections are therefore natural places to audit if the residual data-MC gap is pursued further.

Figure Coverage

All 41 project figures are represented in this README. PDF figures have PNG previews in figures/report_previews/, and the original files remain linked. The compact figure-to-code mapping is in docs/figure_index.md.

What Is Not Here

The repository intentionally excludes raw data and heavy generated products:

• IceCube .i3, SQLite .db, parquet, CSV, NumPy, pickle, HDF5, and ROOT data.
• Model checkpoints and exported weights.
• Slurm logs, cache directories, and local notebook checkpoints.
• Full write-up source and build products.

The project figures are kept so the analysis can be inspected without regenerating every cluster job.

AI Assistance Note

This README and the supporting documentation were assembled with help from OpenAI Codex. The scientific work, analysis ideas, and local project files come from Holger Klevang Christiansen's IceCube preparation project; Codex was used to organize the repository, connect figures to code, and write documentation.