{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "title"
   },
   "source": "# Advanced Quantum Machine Learning (QML) Analysis with CERN CMS J/ψ → μ⁺μ⁻ Data\n\n**Project Type:** Signal–Background Separation + Monte Carlo Sensitivity Analysis  \n**Data Source:** CERN Open Data – CMS Run2011A (Jpsimumu)  \n**Goal:** Fair and statistically robust comparison of classical ML vs Quantum Kernel / VQC methods  \n**Environment:** Google Colab (Qiskit 1.4.x + qiskit-machine-learning 0.8.x)\n\n---\n\n### Project Objective\nThis notebook uses real dimuon data from the CERN CMS experiment to compare classical and quantum machine learning methods for separating the J/ψ resonance signal from background.  \nIt also simulates detector noise and systematic uncertainties via Monte Carlo methods to evaluate model robustness.\n\n### Key Contributions\n- Rich kinematic feature space (pt, η, φ, ΔR, Δφ, …)\n- Stratified Cross-Validation + statistical significance indicators\n- Systematic scan of Quantum Feature Maps and Ansatzes\n- Monte Carlo noise injection for sensitivity analysis\n- Professional visualization and reporting"
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "setup"
   },
   "source": "## 0. Installation and Environment Setup\n\n> **Important:** If this is your first run in Colab, first do **Runtime → Restart session**, then run the installation cell below."
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "install"
   },
   "source": "# Clean previous packages and install stable, mutually compatible versions\n!pip uninstall -y qiskit qiskit-machine-learning qiskit-aer 2>/dev/null || true\n!pip install -q \"qiskit==1.4.2\" \"qiskit-machine-learning==0.8.2\" \"qiskit-aer==0.15.1\" scikit-learn pandas numpy matplotlib seaborn scipy",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "imports"
   },
   "source": "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom sklearn.model_selection import (\n    train_test_split, StratifiedKFold, cross_val_score\n)\nfrom sklearn.preprocessing import RobustScaler\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier\nfrom sklearn.svm import SVC\nfrom sklearn.metrics import (\n    accuracy_score, f1_score, roc_auc_score,\n    confusion_matrix, classification_report\n)\nfrom sklearn.utils import resample\nfrom sklearn.feature_selection import mutual_info_classif\n\n# Qiskit 1.x compatible imports (no Sampler import)\nfrom qiskit.circuit.library import ZZFeatureMap, PauliFeatureMap, RealAmplitudes, EfficientSU2\nfrom qiskit_machine_learning.algorithms import QSVC, VQC\nfrom qiskit_machine_learning.kernels import FidelityQuantumKernel\nfrom qiskit_machine_learning.optimizers import COBYLA\n\n# Plot settings\nplt.style.use('seaborn-v0_8-whitegrid')\nsns.set_palette('viridis')\npd.set_option('display.max_columns', 50)\nnp.random.seed(42)\n\nimport qiskit\nimport qiskit_machine_learning\nprint('Libraries loaded successfully.')\nprint('qiskit                  :', qiskit.__version__)\nprint('qiskit-machine-learning :', qiskit_machine_learning.__version__)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "data_section"
   },
   "source": "## 1. Data Loading and Physics Preprocessing\n\nData source: CERN Open Data Portal – CMS Derived Datasets  \nFile: `Jpsimumu_Run2011A.csv`\n\n**Physics Definition:**\n- **Signal (Label = 1):** Invariant mass window $M_{\\mu\\mu} \\in [3.0,\\, 3.2]$ GeV (J/ψ resonance)\n- **Background (Label = 0):** Events outside the window"
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "data_load"
   },
   "source": "# -------------------------------------------------\n# Data loading\n# -------------------------------------------------\n# Set your own file path or upload the file to Colab\nDATA_PATH = '/content/sample_data/Jpsimumu_Run2011A.csv'\n\ntry:\n    df = pd.read_csv(DATA_PATH)\n    print(f'Data loaded successfully. Shape: {df.shape}')\nexcept FileNotFoundError:\n    print('File not found. Please upload the file to Colab or update the DATA_PATH variable.')\n    raise",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "data_explore"
   },
   "source": "print('--- First 5 rows ---')\ndisplay(df.head())\n\nprint('\\n--- Columns ---')\nprint(df.columns.tolist())\n\nprint('\\n--- Basic statistics ---')\ndisplay(df.describe().T)",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "feature_eng"
   },
   "source": "## 2. Feature Engineering\n\nRich kinematic feature set for classical and quantum models:\n- Basic: $p_{T1},\\, p_{T2},\\, \\eta_1,\\, \\eta_2,\\, \\phi_1,\\, \\phi_2$\n- Derived: $\\Delta\\eta,\\, \\Delta\\phi,\\, \\Delta R,\\, p_T^{\\rm sum},\\, |\\eta|^{\\rm sum}$"
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "features"
   },
   "source": "def engineer_features(df):\n    \"\"\"Create kinematic features.\"\"\"\n    data = df.copy()\n\n    required = ['pt1', 'pt2', 'eta1', 'eta2', 'phi1', 'phi2', 'M']\n    missing = [c for c in required if c not in data.columns]\n    if missing:\n        raise ValueError(f'Missing columns: {missing}. Please check the dataset columns.')\n\n    # Angular differences (periodic)\n    data['delta_phi'] = np.abs(data['phi1'] - data['phi2'])\n    data['delta_phi'] = np.minimum(data['delta_phi'], 2*np.pi - data['delta_phi'])\n    data['delta_eta'] = np.abs(data['eta1'] - data['eta2'])\n    data['delta_R']   = np.sqrt(data['delta_eta']**2 + data['delta_phi']**2)\n\n    data['pt_sum']  = data['pt1'] + data['pt2']\n    data['eta_sum'] = np.abs(data['eta1']) + np.abs(data['eta2'])\n\n    # Target: J/ψ mass window\n    data['Target'] = np.where((data['M'] >= 3.0) & (data['M'] <= 3.2), 1, 0)\n\n    feature_cols = [\n        'pt1', 'pt2', 'eta1', 'eta2', 'phi1', 'phi2',\n        'delta_eta', 'delta_phi', 'delta_R', 'pt_sum', 'eta_sum'\n    ]\n    return data, feature_cols\n\ndf_fe, FEATURES = engineer_features(df)\nprint(f'Number of features : {len(FEATURES)}')\nprint(f'Signal fraction    : {df_fe[\"Target\"].mean():.3f}')\nprint(f'Total events       : {len(df_fe)}')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mass_plot"
   },
   "source": "fig, ax = plt.subplots(1, 2, figsize=(14, 5))\n\nax[0].hist(df_fe['M'], bins=80, range=(2.0, 5.0), color='steelblue', alpha=0.8, edgecolor='black')\nax[0].axvspan(3.0, 3.2, color='red', alpha=0.25, label='J/ψ window')\nax[0].set_xlabel(r'$M_{\\mu\\mu}$ [GeV]')\nax[0].set_ylabel('Number of events')\nax[0].set_title('Dimuon Invariant Mass Spectrum')\nax[0].legend()\n\nsns.countplot(x='Target', data=df_fe, ax=ax[1], palette=['#4c72b0', '#dd8452'])\nax[1].set_xticklabels(['Background (0)', 'Signal (1)'])\nax[1].set_title('Class Distribution')\n\nplt.tight_layout()\nplt.show()",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mc_section"
   },
   "source": "## 3. Monte Carlo Detector Noise Simulation\n\n$$\nx_i^{\\rm noisy} = x_i + \\mathcal{N}(0,\\, \\sigma_i), \\quad \\sigma_i = \\epsilon \\cdot \\mathrm{std}(x_i)\n$$"
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "monte_carlo"
   },
   "source": "def add_detector_noise(X, noise_level=0.02, seed=42):\n    \"\"\"Add Monte Carlo-style detector noise.\"\"\"\n    rng = np.random.default_rng(seed)\n    X_noisy = X.copy().astype(float)\n    for i in range(X.shape[1]):\n        sigma = noise_level * np.std(X[:, i])\n        X_noisy[:, i] += rng.normal(0, sigma, size=X.shape[0])\n    return X_noisy\n\nSAMPLE_SIZE = 3000\ndf_sample = df_fe.sample(n=min(SAMPLE_SIZE, len(df_fe)), random_state=42)\n\nX_raw = df_sample[FEATURES].values\ny     = df_sample['Target'].values\n\nprint(f'Monte Carlo sample size: {X_raw.shape}')\nprint(f'Signal events: {y.sum()} | Background: {(y==0).sum()}')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "classical_section"
   },
   "source": "## 4. Classical Machine Learning Baselines\n\nFair comparison using Stratified 5-Fold Cross-Validation."
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "classical_models"
   },
   "source": "scaler = RobustScaler()\nX_scaled = scaler.fit_transform(X_raw)\n\nclassical_models = {\n    'Logistic Regression': LogisticRegression(max_iter=1000, random_state=42),\n    'Random Forest'      : RandomForestClassifier(n_estimators=100, max_depth=8, random_state=42, n_jobs=-1),\n    'Gradient Boosting'  : GradientBoostingClassifier(n_estimators=80, max_depth=4, random_state=42),\n    'SVM (RBF)'          : SVC(kernel='rbf', probability=True, random_state=42)\n}\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n\nclassical_results = []\nfor name, model in classical_models.items():\n    scores_acc = cross_val_score(model, X_scaled, y, cv=cv, scoring='accuracy', n_jobs=-1)\n    scores_f1  = cross_val_score(model, X_scaled, y, cv=cv, scoring='f1', n_jobs=-1)\n    scores_auc = cross_val_score(model, X_scaled, y, cv=cv, scoring='roc_auc', n_jobs=-1)\n\n    classical_results.append({\n        'Model': name,\n        'Accuracy (mean)': scores_acc.mean(),\n        'Accuracy (std)' : scores_acc.std(),\n        'F1 (mean)'      : scores_f1.mean(),\n        'F1 (std)'       : scores_f1.std(),\n        'ROC-AUC (mean)' : scores_auc.mean(),\n        'ROC-AUC (std)'  : scores_auc.std()\n    })\n    print(f'{name:22s} | Acc: {scores_acc.mean():.3f} ± {scores_acc.std():.3f} | '\n          f'F1: {scores_f1.mean():.3f} ± {scores_f1.std():.3f} | '\n          f'AUC: {scores_auc.mean():.3f} ± {scores_auc.std():.3f}')\n\nclassical_df = pd.DataFrame(classical_results)\ndisplay(classical_df.round(4))",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qml_section"
   },
   "source": "## 5. Quantum Machine Learning (QML)\n\nDue to NISQ constraints, a smaller but balanced subsample is used for QML."
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "qml_prep"
   },
   "source": "QML_TRAIN_SIZE = 120\nQML_TEST_SIZE  = 60\n\nX_train_c, X_test_c, y_train_c, y_test_c = train_test_split(\n    X_scaled, y, test_size=0.25, stratify=y, random_state=42\n)\n\ndef balanced_subsample(X, y, n_samples, seed=42):\n    rng = np.random.default_rng(seed)\n    idx_pos = np.where(y == 1)[0]\n    idx_neg = np.where(y == 0)[0]\n    n_half = n_samples // 2\n    sel_pos = rng.choice(idx_pos, size=min(n_half, len(idx_pos)), replace=False)\n    sel_neg = rng.choice(idx_neg, size=min(n_half, len(idx_neg)), replace=False)\n    sel = np.concatenate([sel_pos, sel_neg])\n    rng.shuffle(sel)\n    return X[sel], y[sel]\n\nX_train_q, y_train_q = balanced_subsample(X_train_c, y_train_c, QML_TRAIN_SIZE)\nX_test_q,  y_test_q  = balanced_subsample(X_test_c,  y_test_c,  QML_TEST_SIZE)\n\nprint(f'QML Train: {X_train_q.shape} | Signal fraction: {y_train_q.mean():.2f}')\nprint(f'QML Test : {X_test_q.shape}  | Signal fraction: {y_test_q.mean():.2f}')\n\nmi = mutual_info_classif(X_train_q, y_train_q, random_state=42)\ntop_idx = np.argsort(mi)[::-1][:4]\nprint('Selected feature indices:', top_idx)\nprint('Corresponding features  :', [FEATURES[i] for i in top_idx])\n\nX_train_q2 = X_train_q[:, top_idx]\nX_test_q2  = X_test_q[:, top_idx]\nN_QUBITS   = X_train_q2.shape[1]\nprint(f'Number of qubits used: {N_QUBITS}')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "qml_models"
   },
   "source": "feature_maps = {\n    'ZZ-linear-r2': ZZFeatureMap(feature_dimension=N_QUBITS, reps=2, entanglement='linear'),\n    'ZZ-full-r2'  : ZZFeatureMap(feature_dimension=N_QUBITS, reps=2, entanglement='full'),\n    'Pauli-r2'    : PauliFeatureMap(feature_dimension=N_QUBITS, reps=2, paulis=['Z', 'YY'])\n}\n\nansatzes = {\n    'RealAmplitudes-r1': RealAmplitudes(num_qubits=N_QUBITS, reps=1),\n    'EfficientSU2-r1'  : EfficientSU2(num_qubits=N_QUBITS, reps=1)\n}\n\noptimizer = COBYLA(maxiter=40)\nqml_results = []\n\nprint('--- Quantum Kernel (QSVC) training ---')\nfor fm_name, fm in feature_maps.items():\n    try:\n        qkernel = FidelityQuantumKernel(feature_map=fm)\n        qsvc = QSVC(quantum_kernel=qkernel)\n        qsvc.fit(X_train_q2, y_train_q)\n        preds = qsvc.predict(X_test_q2)\n        acc = accuracy_score(y_test_q, preds)\n        f1  = f1_score(y_test_q, preds, zero_division=0)\n        qml_results.append({'Model': f'QSVC ({fm_name})', 'Accuracy': acc, 'F1-Score': f1})\n        print(f'  {fm_name:18s} → Acc: {acc:.3f} | F1: {f1:.3f}')\n    except Exception as e:\n        print(f'  {fm_name} error: {e}')\n\nprint('\\n--- VQC training ---')\nfor fm_name, fm in list(feature_maps.items())[:2]:\n    for an_name, ansatz in ansatzes.items():\n        try:\n            vqc = VQC(feature_map=fm, ansatz=ansatz, optimizer=optimizer)\n            vqc.fit(X_train_q2, y_train_q)\n            preds = vqc.predict(X_test_q2)\n            preds = np.where(np.array(preds).ravel() > 0.5, 1, 0)\n            acc = accuracy_score(y_test_q, preds)\n            f1  = f1_score(y_test_q, preds, zero_division=0)\n            qml_results.append({'Model': f'VQC ({fm_name} + {an_name})', 'Accuracy': acc, 'F1-Score': f1})\n            print(f'  {fm_name} + {an_name:18s} → Acc: {acc:.3f} | F1: {f1:.3f}')\n        except Exception as e:\n            print(f'  Error: {e}')\n\nqml_df = pd.DataFrame(qml_results)\ndisplay(qml_df.sort_values('Accuracy', ascending=False).round(4))",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "noise_robust"
   },
   "source": "## 6. Monte Carlo Noise Robustness Analysis"
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "noise_analysis"
   },
   "source": "noise_levels = [0.0, 0.01, 0.03, 0.05, 0.08]\nrobustness_results = []\n\nbest_classical = RandomForestClassifier(n_estimators=100, max_depth=8, random_state=42)\nbest_fm = ZZFeatureMap(feature_dimension=N_QUBITS, reps=2, entanglement='linear')\nbest_qsvc = QSVC(quantum_kernel=FidelityQuantumKernel(feature_map=best_fm))\n\nfor eps in noise_levels:\n    X_noisy = add_detector_noise(X_scaled, noise_level=eps, seed=42)\n    X_tr, X_te, y_tr, y_te = train_test_split(X_noisy, y, test_size=0.25, stratify=y, random_state=42)\n\n    best_classical.fit(X_tr, y_tr)\n    pred_c = best_classical.predict(X_te)\n    acc_c = accuracy_score(y_te, pred_c)\n    f1_c  = f1_score(y_te, pred_c, zero_division=0)\n\n    X_tr_q, y_tr_q = balanced_subsample(X_tr[:, top_idx], y_tr, QML_TRAIN_SIZE)\n    X_te_q, y_te_q = balanced_subsample(X_te[:, top_idx], y_te, QML_TEST_SIZE)\n\n    try:\n        best_qsvc.fit(X_tr_q, y_tr_q)\n        pred_q = best_qsvc.predict(X_te_q)\n        acc_q = accuracy_score(y_te_q, pred_q)\n        f1_q  = f1_score(y_te_q, pred_q, zero_division=0)\n    except Exception as e:\n        print(f'  QSVC error (ε={eps}): {e}')\n        acc_q, f1_q = np.nan, np.nan\n\n    robustness_results.append({\n        'Noise Level (ε)': eps,\n        'Classical Acc': acc_c, 'Classical F1': f1_c,\n        'Quantum Acc': acc_q, 'Quantum F1': f1_q\n    })\n    print(f'ε={eps:.2f} | Classical Acc={acc_c:.3f} | Quantum Acc={acc_q:.3f}')\n\nrobust_df = pd.DataFrame(robustness_results)\ndisplay(robust_df.round(4))",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "noise_plot"
   },
   "source": "fig, ax = plt.subplots(1, 2, figsize=(13, 5))\n\nax[0].plot(robust_df['Noise Level (ε)'], robust_df['Classical Acc'], 'o-', label='Random Forest', linewidth=2)\nax[0].plot(robust_df['Noise Level (ε)'], robust_df['Quantum Acc'], 's--', label='QSVC', linewidth=2)\nax[0].set_xlabel('Noise Level (ε)')\nax[0].set_ylabel('Accuracy')\nax[0].set_title('Accuracy Robustness vs Noise')\nax[0].legend()\nax[0].grid(True, alpha=0.3)\n\nax[1].plot(robust_df['Noise Level (ε)'], robust_df['Classical F1'], 'o-', label='Random Forest', linewidth=2)\nax[1].plot(robust_df['Noise Level (ε)'], robust_df['Quantum F1'], 's--', label='QSVC', linewidth=2)\nax[1].set_xlabel('Noise Level (ε)')\nax[1].set_ylabel('F1-Score')\nax[1].set_title('F1 Robustness vs Noise')\nax[1].legend()\nax[1].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "stats_section"
   },
   "source": "## 7. Results Summary"
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "final_comparison"
   },
   "source": "print('='*60)\nprint('CLASSICAL MODELS (5-Fold CV)')\nprint('='*60)\ndisplay(classical_df.round(4))\n\nprint('\\n' + '='*60)\nprint('QUANTUM MODELS')\nprint('='*60)\ndisplay(qml_df.sort_values('Accuracy', ascending=False).round(4))\n\nbest_classical_name = classical_df.loc[classical_df['Accuracy (mean)'].idxmax(), 'Model']\nbest_qml_name = qml_df.loc[qml_df['Accuracy'].idxmax(), 'Model'] if len(qml_df) > 0 else 'N/A'\n\nprint(f'\\nBest classical model : {best_classical_name}')\nprint(f'Best quantum model   : {best_qml_name}')",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "final_plot"
   },
   "source": "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\nx = np.arange(len(classical_df))\nwidth = 0.35\naxes[0].bar(x - width/2, classical_df['Accuracy (mean)'], width,\n            yerr=classical_df['Accuracy (std)'], label='Accuracy', capsize=4)\naxes[0].bar(x + width/2, classical_df['F1 (mean)'], width,\n            yerr=classical_df['F1 (std)'], label='F1-Score', capsize=4)\naxes[0].set_xticks(x)\naxes[0].set_xticklabels(classical_df['Model'], rotation=20, ha='right')\naxes[0].set_ylim(0, 1.05)\naxes[0].set_title('Classical Models (5-Fold CV)')\naxes[0].legend()\naxes[0].grid(axis='y', alpha=0.3)\n\nif len(qml_df) > 0:\n    qml_sorted = qml_df.sort_values('Accuracy', ascending=True)\n    axes[1].barh(qml_sorted['Model'], qml_sorted['Accuracy'], color='teal', alpha=0.8, label='Accuracy')\n    axes[1].barh(qml_sorted['Model'], qml_sorted['F1-Score'], color='orange', alpha=0.6, label='F1')\n    axes[1].set_xlim(0, 1.05)\n    axes[1].set_title('Quantum Models')\n    axes[1].legend()\n    axes[1].grid(axis='x', alpha=0.3)\n\nplt.tight_layout()\nplt.show()",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "conclusion"
   },
   "source": "## 8. Conclusions and Future Work\n\n### Main Findings\n- Classical models (especially Random Forest / Gradient Boosting) currently outperform QML methods under NISQ simulator constraints and show more stable performance.\n- Quantum Kernel methods (QSVC) generally produce more consistent results than VQC.\n- Monte Carlo noise analysis indicates that models remain reasonably robust under low-to-moderate detector noise.\n\n### Limitations\n- QML training is limited to small sample sizes due to simulator speed.\n- The models have not yet been executed on real quantum hardware.\n\n### Future Steps\n1. Scale to the full CMS dataset and run experiments on real quantum hardware\n2. Incorporate Quantum Neural Network (QNN) approaches\n3. Use more realistic detector simulation\n4. Add bootstrap confidence intervals\n5. Prepare results for academic publication\n\n---\n\n**Note:** Update the data file path before running in Colab. QML cells may take several minutes to complete."
  }
 ]
}