Files
2026-04-26 22:52:37 -04:00

1513 lines
132 KiB
HTML
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.6.41">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<meta name="author" content="John Murphy and Steven Carlson">
<meta name="dcterms.date" content="2025-11-11">
<title>Lab 3: Acquiring and Analyzing Simulation Data – Geant4 Radiation Labs</title>
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1em; /* quarto-specific, see https://github.com/quarto-dev/quarto-cli/issues/4556 */
vertical-align: middle;
}
/* CSS for syntax highlighting */
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
}
pre.numberSource { margin-left: 3em; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
</style>
<script src="../site_libs/quarto-nav/quarto-nav.js"></script>
<script src="../site_libs/quarto-nav/headroom.min.js"></script>
<script src="../site_libs/clipboard/clipboard.min.js"></script>
<script src="../site_libs/quarto-search/autocomplete.umd.js"></script>
<script src="../site_libs/quarto-search/fuse.min.js"></script>
<script src="../site_libs/quarto-search/quarto-search.js"></script>
<meta name="quarto:offset" content="../">
<link href="../lab4_exampleB3/lab4.html" rel="next">
<link href="../lab2_exampleB1/lab2.html" rel="prev">
<script src="../site_libs/quarto-html/quarto.js"></script>
<script src="../site_libs/quarto-html/popper.min.js"></script>
<script src="../site_libs/quarto-html/tippy.umd.min.js"></script>
<script src="../site_libs/quarto-html/anchor.min.js"></script>
<link href="../site_libs/quarto-html/tippy.css" rel="stylesheet">
<link href="../site_libs/quarto-html/quarto-syntax-highlighting-48ffa3e5b9d089919c6712c39e5b00f2.css" rel="stylesheet" id="quarto-text-highlighting-styles">
<script src="../site_libs/bootstrap/bootstrap.min.js"></script>
<link href="../site_libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="../site_libs/bootstrap/bootstrap-c85aa45615bc724621e561d3dca1e930.min.css" rel="stylesheet" append-hash="true" id="quarto-bootstrap" data-mode="light">
<script id="quarto-search-options" type="application/json">{
"location": "navbar",
"copy-button": false,
"collapse-after": 3,
"panel-placement": "end",
"type": "overlay",
"limit": 50,
"keyboard-shortcut": [
"f",
"/",
"s"
],
"show-item-context": false,
"language": {
"search-no-results-text": "No results",
"search-matching-documents-text": "matching documents",
"search-copy-link-title": "Copy link to search",
"search-hide-matches-text": "Hide additional matches",
"search-more-match-text": "more match in this document",
"search-more-matches-text": "more matches in this document",
"search-clear-button-title": "Clear",
"search-text-placeholder": "",
"search-detached-cancel-button-title": "Cancel",
"search-submit-button-title": "Submit",
"search-label": "Search"
}
}</script>
</head>
<body class="nav-sidebar docked nav-fixed">
<div id="quarto-search-results"></div>
<header id="quarto-header" class="headroom fixed-top">
<nav class="navbar navbar-expand-lg " data-bs-theme="dark">
<div class="navbar-container container-fluid">
<div class="navbar-brand-container mx-auto">
<a class="navbar-brand" href="../index.html">
<span class="navbar-title">Geant4 Radiation Labs</span>
</a>
</div>
<div id="quarto-search" class="" title="Search"></div>
<button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarCollapse" aria-controls="navbarCollapse" role="menu" aria-expanded="false" aria-label="Toggle navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
<span class="navbar-toggler-icon"></span>
</button>
<div class="collapse navbar-collapse" id="navbarCollapse">
<ul class="navbar-nav navbar-nav-scroll me-auto">
<li class="nav-item">
<a class="nav-link active" href="../index.html" aria-current="page">
<span class="menu-text">Home</span></a>
</li>
<li class="nav-item">
<a class="nav-link" href="../lab1_setup/lab1.html">
<span class="menu-text">Lab 1: Setup</span></a>
</li>
<li class="nav-item">
<a class="nav-link" href="../lab2_exampleB1/lab2.html">
<span class="menu-text">Lab 2: B1</span></a>
</li>
<li class="nav-item">
<a class="nav-link" href="../lab3_modifyB1/lab3.html">
<span class="menu-text">Lab 3: Modify B1</span></a>
</li>
<li class="nav-item">
<a class="nav-link" href="../lab4_exampleB3/lab4.html">
<span class="menu-text">Lab 4: PET (B3a)</span></a>
</li>
</ul>
<ul class="navbar-nav navbar-nav-scroll ms-auto">
<li class="nav-item compact">
<a class="nav-link" href="https://github.com/jwinterm/geant4-radiation-labs"> <i class="bi bi-github" role="img" aria-label="GitHub repo">
</i>
<span class="menu-text"></span></a>
</li>
</ul>
</div> <!-- /navcollapse -->
<div class="quarto-navbar-tools">
</div>
</div> <!-- /container-fluid -->
</nav>
<nav class="quarto-secondary-nav">
<div class="container-fluid d-flex">
<button type="button" class="quarto-btn-toggle btn" data-bs-toggle="collapse" role="button" data-bs-target=".quarto-sidebar-collapse-item" aria-controls="quarto-sidebar" aria-expanded="false" aria-label="Toggle sidebar navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
<i class="bi bi-layout-text-sidebar-reverse"></i>
</button>
<nav class="quarto-page-breadcrumbs" aria-label="breadcrumb"><ol class="breadcrumb"><li class="breadcrumb-item"><a href="../lab1_setup/lab1.html">Labs</a></li><li class="breadcrumb-item"><a href="../lab3_modifyB1/lab3.html">Lab 3: Acquiring and Analyzing Simulation Data</a></li></ol></nav>
<a class="flex-grow-1" role="navigation" data-bs-toggle="collapse" data-bs-target=".quarto-sidebar-collapse-item" aria-controls="quarto-sidebar" aria-expanded="false" aria-label="Toggle sidebar navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
</a>
<button type="button" class="btn quarto-search-button" aria-label="Search" onclick="window.quartoOpenSearch();">
<i class="bi bi-search"></i>
</button>
</div>
</nav>
</header>
<!-- content -->
<div id="quarto-content" class="quarto-container page-columns page-rows-contents page-layout-article page-navbar">
<!-- sidebar -->
<nav id="quarto-sidebar" class="sidebar collapse collapse-horizontal quarto-sidebar-collapse-item sidebar-navigation docked overflow-auto">
<div class="mt-2 flex-shrink-0 align-items-center">
<div class="sidebar-search">
<div id="quarto-search" class="" title="Search"></div>
</div>
</div>
<div class="sidebar-menu-container">
<ul class="list-unstyled mt-1">
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="../index.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Introduction</span></a>
</div>
</li>
<li class="sidebar-item sidebar-item-section">
<div class="sidebar-item-container">
<a class="sidebar-item-text sidebar-link text-start" data-bs-toggle="collapse" data-bs-target="#quarto-sidebar-section-1" role="navigation" aria-expanded="true">
<span class="menu-text">Labs</span></a>
<a class="sidebar-item-toggle text-start" data-bs-toggle="collapse" data-bs-target="#quarto-sidebar-section-1" role="navigation" aria-expanded="true" aria-label="Toggle section">
<i class="bi bi-chevron-right ms-2"></i>
</a>
</div>
<ul id="quarto-sidebar-section-1" class="collapse list-unstyled sidebar-section depth1 show">
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="../lab1_setup/lab1.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Lab 1: Setting Up Linux and Building Geant4</span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="../lab2_exampleB1/lab2.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Lab 2: First Simulation - Building and Running Geant4 Example B1</span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="../lab3_modifyB1/lab3.html" class="sidebar-item-text sidebar-link active">
<span class="menu-text">Lab 3: Acquiring and Analyzing Simulation Data</span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="../lab4_exampleB3/lab4.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Lab 4: Medical Physics Application - A PET Scanner</span></a>
</div>
</li>
</ul>
</li>
</ul>
</div>
</nav>
<div id="quarto-sidebar-glass" class="quarto-sidebar-collapse-item" data-bs-toggle="collapse" data-bs-target=".quarto-sidebar-collapse-item"></div>
<!-- margin-sidebar -->
<div id="quarto-margin-sidebar" class="sidebar margin-sidebar">
<nav id="TOC" role="doc-toc" class="toc-active">
<h2 id="toc-title">On this page</h2>
<ul>
<li><a href="#learning-objectives" id="toc-learning-objectives" class="nav-link active" data-scroll-target="#learning-objectives">Learning Objectives</a></li>
<li><a href="#introduction" id="toc-introduction" class="nav-link" data-scroll-target="#introduction">Introduction</a></li>
<li><a href="#prerequisites" id="toc-prerequisites" class="nav-link" data-scroll-target="#prerequisites">Prerequisites</a></li>
<li><a href="#pre-lab-questions" id="toc-pre-lab-questions" class="nav-link" data-scroll-target="#pre-lab-questions">Pre-Lab Questions</a></li>
<li><a href="#step-1-set-up-the-lab-and-review-the-code-structure" id="toc-step-1-set-up-the-lab-and-review-the-code-structure" class="nav-link" data-scroll-target="#step-1-set-up-the-lab-and-review-the-code-structure">Step 1: Set Up the Lab and Review the Code Structure</a></li>
<li><a href="#step-2-modify-c-code-to-save-data" id="toc-step-2-modify-c-code-to-save-data" class="nav-link" data-scroll-target="#step-2-modify-c-code-to-save-data">Step 2: Modify C++ Code to Save Data</a>
<ul class="collapse">
<li><a href="#what-the-example-already-does" id="toc-what-the-example-already-does" class="nav-link" data-scroll-target="#what-the-example-already-does">What the example already does</a></li>
<li><a href="#what-well-change" id="toc-what-well-change" class="nav-link" data-scroll-target="#what-well-change">What we’ll change</a></li>
<li><a href="#modification-1-runaction.cc" id="toc-modification-1-runaction.cc" class="nav-link" data-scroll-target="#modification-1-runaction.cc">Modification 1: RunAction.cc</a></li>
<li><a href="#modification-2-eventaction.cc" id="toc-modification-2-eventaction.cc" class="nav-link" data-scroll-target="#modification-2-eventaction.cc">Modification 2: EventAction.cc</a></li>
<li><a href="#summary-of-modifications" id="toc-summary-of-modifications" class="nav-link" data-scroll-target="#summary-of-modifications">Summary of Modifications</a></li>
</ul></li>
<li><a href="#step-3-compile-and-run-the-modified-application" id="toc-step-3-compile-and-run-the-modified-application" class="nav-link" data-scroll-target="#step-3-compile-and-run-the-modified-application">Step 3: Compile and Run the Modified Application</a></li>
<li><a href="#step-4-install-python-analysis-libraries" id="toc-step-4-install-python-analysis-libraries" class="nav-link" data-scroll-target="#step-4-install-python-analysis-libraries">Step 4: Install Python Analysis Libraries</a></li>
<li><a href="#step-5-analyze-the-data-with-python---basic-analysis" id="toc-step-5-analyze-the-data-with-python---basic-analysis" class="nav-link" data-scroll-target="#step-5-analyze-the-data-with-python---basic-analysis">Step 5: Analyze the Data with Python - Basic Analysis</a></li>
<li><a href="#step-6-advanced-analysis---data-filtering-and-comparison" id="toc-step-6-advanced-analysis---data-filtering-and-comparison" class="nav-link" data-scroll-target="#step-6-advanced-analysis---data-filtering-and-comparison">Step 6: Advanced Analysis - Data Filtering and Comparison</a></li>
<li><a href="#step-7-comparing-different-particle-types" id="toc-step-7-comparing-different-particle-types" class="nav-link" data-scroll-target="#step-7-comparing-different-particle-types">Step 7: Comparing Different Particle Types</a></li>
<li><a href="#analysis-questions" id="toc-analysis-questions" class="nav-link" data-scroll-target="#analysis-questions">Analysis Questions</a>
<ul class="collapse">
<li><a href="#part-1-understanding-the-distribution" id="toc-part-1-understanding-the-distribution" class="nav-link" data-scroll-target="#part-1-understanding-the-distribution">Part 1: Understanding the Distribution</a></li>
<li><a href="#part-2-understanding-zero-energy-events" id="toc-part-2-understanding-zero-energy-events" class="nav-link" data-scroll-target="#part-2-understanding-zero-energy-events">Part 2: Understanding Zero Energy Events</a></li>
<li><a href="#part-3-comparing-particle-types" id="toc-part-3-comparing-particle-types" class="nav-link" data-scroll-target="#part-3-comparing-particle-types">Part 3: Comparing Particle Types</a></li>
<li><a href="#part-4-pandas-data-analysis" id="toc-part-4-pandas-data-analysis" class="nav-link" data-scroll-target="#part-4-pandas-data-analysis">Part 4: Pandas Data Analysis</a></li>
<li><a href="#part-5-scientific-workflow" id="toc-part-5-scientific-workflow" class="nav-link" data-scroll-target="#part-5-scientific-workflow">Part 5: Scientific Workflow</a></li>
</ul></li>
<li><a href="#troubleshooting" id="toc-troubleshooting" class="nav-link" data-scroll-target="#troubleshooting">Troubleshooting</a></li>
<li><a href="#verification-checklist" id="toc-verification-checklist" class="nav-link" data-scroll-target="#verification-checklist">Verification Checklist</a></li>
<li><a href="#conclusion" id="toc-conclusion" class="nav-link" data-scroll-target="#conclusion">Conclusion</a></li>
<li><a href="#references" id="toc-references" class="nav-link" data-scroll-target="#references">References</a></li>
<li><a href="#additional-resources" id="toc-additional-resources" class="nav-link" data-scroll-target="#additional-resources">Additional Resources</a></li>
<li><a href="#optional-extensions" id="toc-optional-extensions" class="nav-link" data-scroll-target="#optional-extensions">Optional Extensions</a></li>
</ul>
</nav>
</div>
<!-- main -->
<main class="content" id="quarto-document-content">
<header id="title-block-header" class="quarto-title-block default"><nav class="quarto-page-breadcrumbs quarto-title-breadcrumbs d-none d-lg-block" aria-label="breadcrumb"><ol class="breadcrumb"><li class="breadcrumb-item"><a href="../lab1_setup/lab1.html">Labs</a></li><li class="breadcrumb-item"><a href="../lab3_modifyB1/lab3.html">Lab 3: Acquiring and Analyzing Simulation Data</a></li></ol></nav>
<div class="quarto-title">
<h1 class="title">Lab 3: Acquiring and Analyzing Simulation Data</h1>
</div>
<div class="quarto-title-meta">
<div>
<div class="quarto-title-meta-heading">Author</div>
<div class="quarto-title-meta-contents">
<p>John Murphy and Steven Carlson </p>
</div>
</div>
<div>
<div class="quarto-title-meta-heading">Published</div>
<div class="quarto-title-meta-contents">
<p class="date">November 11, 2025</p>
</div>
</div>
</div>
</header>
<section id="learning-objectives" class="level2">
<h2 class="anchored" data-anchor-id="learning-objectives">Learning Objectives</h2>
<p>In this lab, you will bridge the gap between simulation and analysis. You will modify the C++ source code of a Geant4 application to save data and then use Python to analyze and visualize that data. By the end of this session, you will be able to:</p>
<ul>
<li>Understand the basic structure of Geant4 user action classes</li>
<li>Modify C++ code to implement the <code>G4AnalysisManager</code> for data output</li>
<li>Save simulation results to CSV file format for universal compatibility</li>
<li>Write Python scripts using the <code>pandas</code> and <code>matplotlib</code> libraries</li>
<li>Load simulation data into pandas DataFrames</li>
<li>Create and customize histograms to visualize distributions</li>
<li>Calculate descriptive statistics (mean, standard deviation, percentiles) from data</li>
<li>Filter and manipulate data using pandas operations</li>
<li>Compare multiple datasets and create multi-plot figures</li>
</ul>
</section>
<section id="introduction" class="level2">
<h2 class="anchored" data-anchor-id="introduction">Introduction</h2>
<p>A simulation that produces no output is not very useful. To perform scientific analysis, we need to “score” quantities of interest—like the energy deposited in a detector—and save this data to a file. In this lab, you will learn the modern, standard way to do this in Geant4 using a special class called <code>G4AnalysisManager</code>. We will configure it to save the energy deposited in each event into a simple text-based CSV (Comma-Separated Values) file.</p>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/workflow_diagram.png" class="img-fluid figure-img"></p>
<figcaption>Data Analysis Workflow</figcaption>
</figure>
</div>
<p>Once we have the data, we will switch hats from a simulation user to a data analyst. We will use the powerful Python programming language along with two essential libraries: <code>pandas</code> for data manipulation and <code>matplotlib</code> for plotting. You will write a script to read the CSV file and generate a histogram, which is a graph that shows the frequency distribution of the energy deposition values. This process—simulate, save, analyze, visualize—is a fundamental cycle in computational science.</p>
</section>
<section id="prerequisites" class="level2">
<h2 class="anchored" data-anchor-id="prerequisites">Prerequisites</h2>
<p>Before you begin, ensure you have the following:</p>
<ul>
<li>Completed Lab 2 (successfully built and ran Example B1)</li>
<li>Text editor for C++ code (gedit, nano, or VS Code)</li>
<li>Basic understanding of C++ syntax (we’ll guide you through the modifications)</li>
<li>Python 3 installed (should be pre-installed on Ubuntu)</li>
</ul>
</section>
<section id="pre-lab-questions" class="level2">
<h2 class="anchored" data-anchor-id="pre-lab-questions">Pre-Lab Questions</h2>
<ol type="1">
<li>What is a CSV file? Why is it a good choice for this lab compared to a more complex format like ROOT?</li>
<li>What is a histogram and what does it show?</li>
<li>Briefly describe the roles of the <code>pandas</code> and <code>matplotlib</code> Python libraries.</li>
<li>In Geant4, what is the difference between a “Run”, an “Event”, and a “Step”?</li>
<li>What does it mean to “score” a quantity in a simulation?</li>
</ol>
</section>
<section id="step-1-set-up-the-lab-and-review-the-code-structure" class="level2">
<h2 class="anchored" data-anchor-id="step-1-set-up-the-lab-and-review-the-code-structure">Step 1: Set Up the Lab and Review the Code Structure</h2>
<ol type="1">
<li><p>Create a new directory for this lab and copy the B1 example into it, just as you did in Lab 2:</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="fu">mkdir</span> <span class="at">-p</span> ~/g4-labs/lab3</span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> ~/g4-labs/lab3</span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a><span class="fu">cp</span> <span class="at">-r</span> ~/geant4-v11.3.2/examples/basic/B1 .</span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> B1</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Before we start modifying, let’s understand the key files and the Geant4 action hierarchy:</p>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/action_classes.png" class="img-fluid figure-img"></p>
<figcaption>Geant4 User Action Classes</figcaption>
</figure>
</div>
<p>In the modern B1 example, all of these classes live inside <code>namespace B1 { ... }</code>. The header and source files do <strong>not</strong> carry a <code>B1</code> filename prefix — for example, the run-action class is declared as <code>class RunAction</code> inside <code>namespace B1</code>, in <code>include/RunAction.hh</code>. Every code addition you make in this lab will go inside an existing <code>namespace B1</code> block in the corresponding <code>.cc</code> file.</p>
<p>Familiarize yourself with these key files:</p>
<ul>
<li><strong><code>include/RunAction.hh</code></strong> &amp; <strong><code>src/RunAction.cc</code></strong>:
<ul>
<li>Controls actions at the beginning and end of a <strong>Run</strong> (a collection of events)</li>
<li>We will edit the source file to open, write, and close the CSV output</li>
</ul></li>
<li><strong><code>include/EventAction.hh</code></strong> &amp; <strong><code>src/EventAction.cc</code></strong>:
<ul>
<li>Controls actions at the beginning and end of each <strong>Event</strong> (simulation of one primary particle)</li>
<li>We will edit the source file to record one row of CSV data per event</li>
</ul></li>
<li><strong><code>src/SteppingAction.cc</code></strong>:
<ul>
<li>Executed at every single <strong>Step</strong> a particle takes through the geometry</li>
<li>The example already collects per-step energy here. We will <strong>not</strong> modify this file.</li>
</ul></li>
</ul></li>
<li><p>Open <code>src/SteppingAction.cc</code> in a text editor and read through <code>UserSteppingAction</code>:</p>
<div class="sourceCode" id="cb2"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="ex">gedit</span> src/SteppingAction.cc <span class="kw">&amp;</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Notice that the example already extracts each step’s energy deposit and adds it to the event total via <code>fEventAction-&gt;AddEdep(edepStep)</code>. The data lives in memory but never reaches disk — fixing that is what this lab is about.</p></li>
</ol>
</section>
<section id="step-2-modify-c-code-to-save-data" class="level2">
<h2 class="anchored" data-anchor-id="step-2-modify-c-code-to-save-data">Step 2: Modify C++ Code to Save Data</h2>
<section id="what-the-example-already-does" class="level3">
<h3 class="anchored" data-anchor-id="what-the-example-already-does">What the example already does</h3>
<p>Before we add anything, look at how energy already flows through the unmodified B1 example. Open these three files in your editor and find the lines below — you do <strong>not</strong> need to change them:</p>
<ul>
<li><p><code>include/EventAction.hh</code> already declares the helper and the per-event accumulator:</p>
<div class="sourceCode" id="cb3"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="dt">void</span> AddEdep<span class="op">(</span>G4double edep<span class="op">)</span> <span class="op">{</span> fEdep <span class="op">+=</span> edep<span class="op">;</span> <span class="op">}</span></span>
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true" tabindex="-1"></a>G4double fEdep <span class="op">=</span> <span class="fl">0.</span><span class="op">;</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p><code>src/EventAction.cc</code> already resets <code>fEdep</code> at the start of each event and forwards it to the run action at the end:</p>
<div class="sourceCode" id="cb4"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a><span class="dt">void</span> EventAction<span class="op">::</span>BeginOfEventAction<span class="op">(</span><span class="at">const</span> G4Event<span class="op">*)</span> <span class="op">{</span> fEdep <span class="op">=</span> <span class="fl">0.</span><span class="op">;</span> <span class="op">}</span></span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a><span class="dt">void</span> EventAction<span class="op">::</span>EndOfEventAction<span class="op">(</span><span class="at">const</span> G4Event<span class="op">*)</span> <span class="op">{</span> fRunAction<span class="op">-&gt;</span>AddEdep<span class="op">(</span>fEdep<span class="op">);</span> <span class="op">}</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p><code>src/SteppingAction.cc</code> already accumulates each step’s deposit into the event total:</p>
<div class="sourceCode" id="cb5"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>G4double edepStep <span class="op">=</span> step<span class="op">-&gt;</span>GetTotalEnergyDeposit<span class="op">();</span></span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a>fEventAction<span class="op">-&gt;</span>AddEdep<span class="op">(</span>edepStep<span class="op">);</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
</ul>
<p>So the per-step → per-event → per-run pipeline is <strong>already wired up</strong>. The energy lands in <code>fEdep</code> at the end of each event. Our only job in this lab is to also write that value to a CSV file.</p>
</section>
<section id="what-well-change" class="level3">
<h3 class="anchored" data-anchor-id="what-well-change">What we’ll change</h3>
<p>There are exactly two files to edit:</p>
<ol type="1">
<li><strong><code>src/RunAction.cc</code></strong> — set up the CSV file at the start of the run and write/close it at the end.</li>
<li><strong><code>src/EventAction.cc</code></strong> — at the end of each event, write <code>fEdep</code> as one new row of the CSV.</li>
</ol>
<p>That’s it. Two <code>#include</code> additions and ~12 new lines of code in total.</p>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/analysis_manager_flow.png" class="img-fluid figure-img"></p>
<figcaption>G4AnalysisManager Workflow</figcaption>
</figure>
</div>
</section>
<section id="modification-1-runaction.cc" class="level3">
<h3 class="anchored" data-anchor-id="modification-1-runaction.cc">Modification 1: RunAction.cc</h3>
<p>We use <code>G4AnalysisManager</code>, Geant4’s built-in writer that supports CSV, ROOT, XML, and HDF5. We’ll use CSV.</p>
<ol type="1">
<li><p>Open <code>src/RunAction.cc</code>:</p>
<div class="sourceCode" id="cb6"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a><span class="ex">gedit</span> src/RunAction.cc <span class="kw">&amp;</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Add the analysis-manager header near the other <code>#include</code> lines at the top:</p>
<div class="sourceCode" id="cb7"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"G4AnalysisManager.hh"</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Inside <code>RunAction::BeginOfRunAction(const G4Run*)</code>, add the following block (anywhere inside the function body is fine — putting it at the top is conventional):</p>
<div class="sourceCode" id="cb8"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a><span class="co">// --- Lab 3 addition: open CSV and define the per-event ntuple ---</span></span>
<span id="cb8-2"><a href="#cb8-2" aria-hidden="true" tabindex="-1"></a><span class="kw">auto</span> analysisManager <span class="op">=</span> G4AnalysisManager<span class="op">::</span>Instance<span class="op">();</span></span>
<span id="cb8-3"><a href="#cb8-3" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>SetDefaultFileType<span class="op">(</span><span class="st">"csv"</span><span class="op">);</span></span>
<span id="cb8-4"><a href="#cb8-4" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>OpenFile<span class="op">(</span><span class="st">"B1_output"</span><span class="op">);</span></span>
<span id="cb8-5"><a href="#cb8-5" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb8-6"><a href="#cb8-6" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>CreateNtuple<span class="op">(</span><span class="st">"B1"</span><span class="op">,</span> <span class="st">"Energy Deposition Data"</span><span class="op">);</span></span>
<span id="cb8-7"><a href="#cb8-7" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>CreateNtupleDColumn<span class="op">(</span><span class="st">"Edep"</span><span class="op">);</span> <span class="co">// column 0</span></span>
<span id="cb8-8"><a href="#cb8-8" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>FinishNtuple<span class="op">();</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p><strong>What this does:</strong> Tells the analysis manager to write CSV, gives the file a name (<code>B1_output</code>), and declares one ntuple column called <code>Edep</code>. Each row will hold one event’s <code>fEdep</code>.</p></li>
<li><p>Inside <code>RunAction::EndOfRunAction(const G4Run* run)</code>, add the following at the <strong>very end</strong> of the function (after the existing <code>G4cout</code> print block):</p>
<div class="sourceCode" id="cb9"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="co">// --- Lab 3 addition: flush and close the CSV ---</span></span>
<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a><span class="kw">auto</span> analysisManager <span class="op">=</span> G4AnalysisManager<span class="op">::</span>Instance<span class="op">();</span></span>
<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>Write<span class="op">();</span></span>
<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>CloseFile<span class="op">();</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
</ol>
</section>
<section id="modification-2-eventaction.cc" class="level3">
<h3 class="anchored" data-anchor-id="modification-2-eventaction.cc">Modification 2: EventAction.cc</h3>
<ol type="1">
<li><p>Open <code>src/EventAction.cc</code>:</p>
<div class="sourceCode" id="cb10"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a><span class="ex">gedit</span> src/EventAction.cc <span class="kw">&amp;</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Add the analysis-manager header near the top:</p>
<div class="sourceCode" id="cb11"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"G4AnalysisManager.hh"</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Find the existing <code>EventAction::EndOfEventAction</code> function. It currently looks like this — <strong>do not retype it</strong>, you are only locating it:</p>
<div class="sourceCode" id="cb12"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a><span class="dt">void</span> EventAction<span class="op">::</span>EndOfEventAction<span class="op">(</span><span class="at">const</span> G4Event<span class="op">*)</span></span>
<span id="cb12-2"><a href="#cb12-2" aria-hidden="true" tabindex="-1"></a><span class="op">{</span></span>
<span id="cb12-3"><a href="#cb12-3" aria-hidden="true" tabindex="-1"></a> <span class="co">// accumulate statistics in run action</span></span>
<span id="cb12-4"><a href="#cb12-4" aria-hidden="true" tabindex="-1"></a> fRunAction<span class="op">-&gt;</span>AddEdep<span class="op">(</span>fEdep<span class="op">);</span></span>
<span id="cb12-5"><a href="#cb12-5" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p><strong>Inside</strong> that function, <em>after</em> the existing <code>fRunAction-&gt;AddEdep(fEdep);</code> line and <em>before</em> the closing <code>}</code>, paste these four lines (and only these four lines):</p>
<div class="sourceCode" id="cb13"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a><span class="co">// --- Lab 3 addition: also write fEdep to the CSV ---</span></span>
<span id="cb13-2"><a href="#cb13-2" aria-hidden="true" tabindex="-1"></a><span class="kw">auto</span> analysisManager <span class="op">=</span> G4AnalysisManager<span class="op">::</span>Instance<span class="op">();</span></span>
<span id="cb13-3"><a href="#cb13-3" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>FillNtupleDColumn<span class="op">(</span><span class="dv">0</span><span class="op">,</span> fEdep<span class="op">);</span></span>
<span id="cb13-4"><a href="#cb13-4" aria-hidden="true" tabindex="-1"></a>analysisManager<span class="op">-&gt;</span>AddNtupleRow<span class="op">();</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>After your edit, the complete function should read:</p>
<div class="sourceCode" id="cb14"><pre class="sourceCode cpp code-with-copy"><code class="sourceCode cpp"><span id="cb14-1"><a href="#cb14-1" aria-hidden="true" tabindex="-1"></a><span class="dt">void</span> EventAction<span class="op">::</span>EndOfEventAction<span class="op">(</span><span class="at">const</span> G4Event<span class="op">*)</span></span>
<span id="cb14-2"><a href="#cb14-2" aria-hidden="true" tabindex="-1"></a><span class="op">{</span></span>
<span id="cb14-3"><a href="#cb14-3" aria-hidden="true" tabindex="-1"></a> <span class="co">// accumulate statistics in run action</span></span>
<span id="cb14-4"><a href="#cb14-4" aria-hidden="true" tabindex="-1"></a> fRunAction<span class="op">-&gt;</span>AddEdep<span class="op">(</span>fEdep<span class="op">);</span></span>
<span id="cb14-5"><a href="#cb14-5" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb14-6"><a href="#cb14-6" aria-hidden="true" tabindex="-1"></a> <span class="co">// --- Lab 3 addition: also write fEdep to the CSV ---</span></span>
<span id="cb14-7"><a href="#cb14-7" aria-hidden="true" tabindex="-1"></a> <span class="kw">auto</span> analysisManager <span class="op">=</span> G4AnalysisManager<span class="op">::</span>Instance<span class="op">();</span></span>
<span id="cb14-8"><a href="#cb14-8" aria-hidden="true" tabindex="-1"></a> analysisManager<span class="op">-&gt;</span>FillNtupleDColumn<span class="op">(</span><span class="dv">0</span><span class="op">,</span> fEdep<span class="op">);</span></span>
<span id="cb14-9"><a href="#cb14-9" aria-hidden="true" tabindex="-1"></a> analysisManager<span class="op">-&gt;</span>AddNtupleRow<span class="op">();</span></span>
<span id="cb14-10"><a href="#cb14-10" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p><strong>What this does:</strong> Each event ends with one new row in the CSV containing the total energy deposited in the scoring volume during that event.</p></li>
</ol>
</section>
<section id="summary-of-modifications" class="level3">
<h3 class="anchored" data-anchor-id="summary-of-modifications">Summary of Modifications</h3>
<p>You’ve now modified two files:</p>
<table class="caption-top table">
<colgroup>
<col style="width: 42%">
<col style="width: 57%">
</colgroup>
<thead>
<tr class="header">
<th>File</th>
<th>Change</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td><code>src/RunAction.cc</code></td>
<td>Add <code>G4AnalysisManager.hh</code>; open the CSV in <code>BeginOfRunAction</code>, write/close it in <code>EndOfRunAction</code></td>
</tr>
<tr class="even">
<td><code>src/EventAction.cc</code></td>
<td>Add <code>G4AnalysisManager.hh</code>; in <code>EndOfEventAction</code>, fill one ntuple row per event</td>
</tr>
</tbody>
</table>
<p>Note that <code>EventAction.hh</code>, <code>SteppingAction.cc</code>, and the energy-accumulation logic itself were left untouched — the example already had everything we needed up to the point of writing to disk.</p>
</section>
</section>
<section id="step-3-compile-and-run-the-modified-application" class="level2">
<h2 class="anchored" data-anchor-id="step-3-compile-and-run-the-modified-application">Step 3: Compile and Run the Modified Application</h2>
<p>Now let’s compile your modified code and generate some data!</p>
<ol type="1">
<li><p>Navigate to the build directory (create one if it doesn’t exist):</p>
<div class="sourceCode" id="cb15"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb15-1"><a href="#cb15-1" aria-hidden="true" tabindex="-1"></a><span class="fu">mkdir</span> <span class="at">-p</span> build</span>
<span id="cb15-2"><a href="#cb15-2" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> build</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Configure the build with CMake:</p>
<div class="sourceCode" id="cb16"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb16-1"><a href="#cb16-1" aria-hidden="true" tabindex="-1"></a><span class="fu">cmake</span> ..</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Compile the modified application:</p>
<div class="sourceCode" id="cb17"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb17-1"><a href="#cb17-1" aria-hidden="true" tabindex="-1"></a><span class="fu">make</span> <span class="at">-j</span><span class="va">$(</span><span class="fu">nproc</span><span class="va">)</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>If you encounter errors, carefully check your code modifications against the instructions. Common issues:</p>
<ul>
<li>Missing semicolons</li>
<li>Incorrect placement of code</li>
<li>Typos in variable names</li>
</ul></li>
</ol>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/compile_success.png" class="img-fluid figure-img"></p>
<figcaption>Successful Compilation</figcaption>
</figure>
</div>
<ol start="4" type="1">
<li><p>Create a macro file for data collection. In the <code>B1</code> directory (not <code>build</code>), create <code>lab3.mac</code>:</p>
<div class="sourceCode" id="cb18"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb18-1"><a href="#cb18-1" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> ..</span>
<span id="cb18-2"><a href="#cb18-2" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> lab3.mac</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Add these contents:</p>
<div class="sourceCode" id="cb19"><pre class="sourceCode default code-with-copy"><code class="sourceCode default"><span id="cb19-1"><a href="#cb19-1" aria-hidden="true" tabindex="-1"></a># Lab 3 Data Collection Run</span>
<span id="cb19-2"><a href="#cb19-2" aria-hidden="true" tabindex="-1"></a># Simulate 10,000 events with default 6 MeV gammas</span>
<span id="cb19-3"><a href="#cb19-3" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb19-4"><a href="#cb19-4" aria-hidden="true" tabindex="-1"></a># Force single-threaded mode so we get ONE CSV file, not one per worker.</span>
<span id="cb19-5"><a href="#cb19-5" aria-hidden="true" tabindex="-1"></a># Must come BEFORE /run/initialize.</span>
<span id="cb19-6"><a href="#cb19-6" aria-hidden="true" tabindex="-1"></a>/run/numberOfThreads 1</span>
<span id="cb19-7"><a href="#cb19-7" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb19-8"><a href="#cb19-8" aria-hidden="true" tabindex="-1"></a>/run/initialize</span>
<span id="cb19-9"><a href="#cb19-9" aria-hidden="true" tabindex="-1"></a>/run/printProgress 1000</span>
<span id="cb19-10"><a href="#cb19-10" aria-hidden="true" tabindex="-1"></a>/run/beamOn 10000</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>This will run 10,000 events and print progress every 1,000 events. <code>/run/numberOfThreads 1</code> keeps Geant4 from running multi-threaded, which would otherwise produce one CSV file per worker thread.</p></li>
<li><p>Run the simulation:</p>
<div class="sourceCode" id="cb20"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb20-1"><a href="#cb20-1" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> build</span>
<span id="cb20-2"><a href="#cb20-2" aria-hidden="true" tabindex="-1"></a><span class="ex">./exampleB1</span> ../lab3.mac</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>After completion, list the CSV files:</p>
<div class="sourceCode" id="cb21"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb21-1"><a href="#cb21-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ls</span> <span class="at">-lh</span> <span class="pp">*</span>.csv</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>You should see one file named <code>B1_output_nt_B1_t0.csv</code>. Geant4’s CSV writer always appends a <code>_tN</code> worker-thread suffix; with one thread, that’s <code>_t0</code>. If you see multiple <code>_tN.csv</code> files, you forgot <code>/run/numberOfThreads 1</code> — delete them and re-run.</p></li>
<li><p>View the first few lines:</p>
<div class="sourceCode" id="cb22"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb22-1"><a href="#cb22-1" aria-hidden="true" tabindex="-1"></a><span class="fu">head</span> <span class="at">-20</span> B1_output_nt_B1_t0.csv</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>You should see a header row followed by numerical energy values in MeV.</p></li>
</ol>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/csv_preview.png" class="img-fluid figure-img"></p>
<figcaption>CSV Output Preview</figcaption>
</figure>
</div>
</section>
<section id="step-4-install-python-analysis-libraries" class="level2">
<h2 class="anchored" data-anchor-id="step-4-install-python-analysis-libraries">Step 4: Install Python Analysis Libraries</h2>
<p>Before we can analyze our data, we need to install the required Python packages.</p>
<ol type="1">
<li><p>First, check if pip is installed:</p>
<div class="sourceCode" id="cb23"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb23-1"><a href="#cb23-1" aria-hidden="true" tabindex="-1"></a><span class="ex">pip3</span> <span class="at">--version</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>If needed, install pip:</p>
<div class="sourceCode" id="cb24"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb24-1"><a href="#cb24-1" aria-hidden="true" tabindex="-1"></a><span class="fu">sudo</span> apt install python3-pip</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Install pandas and matplotlib:</p>
<div class="sourceCode" id="cb25"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb25-1"><a href="#cb25-1" aria-hidden="true" tabindex="-1"></a><span class="ex">pip3</span> install <span class="at">--break-system-packages</span> pandas matplotlib numpy</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p><strong>Why the <code>--break-system-packages</code> flag?</strong> Recent versions of Ubuntu (24.04+) follow <a href="https://peps.python.org/pep-0668/">PEP 668</a> and mark the system Python environment as “externally managed”, so a plain <code>pip3 install pandas</code> will refuse to run with an <code>error: externally-managed-environment</code> message. The flag tells pip to install into the user site-packages directory anyway, which is fine for this VM-based course where we are not maintaining the host system. (In a production environment you would create a virtual environment with <code>python3 -m venv</code> instead.)</p>
<p>These packages provide:</p>
<ul>
<li><strong>pandas</strong>: Data manipulation and analysis</li>
<li><strong>matplotlib</strong>: Creating visualizations and plots</li>
<li><strong>numpy</strong>: Numerical computing (dependency for pandas)</li>
</ul></li>
</ol>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/pip_install.png" class="img-fluid figure-img"></p>
<figcaption>Python Package Installation</figcaption>
</figure>
</div>
<ol start="4" type="1">
<li><p>Verify the installation:</p>
<div class="sourceCode" id="cb26"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb26-1"><a href="#cb26-1" aria-hidden="true" tabindex="-1"></a><span class="ex">python3</span> <span class="at">-c</span> <span class="st">"import pandas; import matplotlib; print('Success!')"</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
</ol>
</section>
<section id="step-5-analyze-the-data-with-python---basic-analysis" class="level2">
<h2 class="anchored" data-anchor-id="step-5-analyze-the-data-with-python---basic-analysis">Step 5: Analyze the Data with Python - Basic Analysis</h2>
<p>Now let’s create a Python script to analyze our simulation data.</p>
<p>In your <code>B1/build</code> directory, create <code>analyze_basic.py</code>:</p>
<div class="sourceCode" id="cb27"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb27-1"><a href="#cb27-1" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> analyze_basic.py</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Paste the following content. <strong>Important</strong>: the code must start at column 0 (no leading whitespace), or Python will reject it with <code>IndentationError: unexpected indent</code>.</p>
<div class="sourceCode" id="cb28"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb28-1"><a href="#cb28-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> pandas <span class="im">as</span> pd</span>
<span id="cb28-2"><a href="#cb28-2" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> matplotlib.pyplot <span class="im">as</span> plt</span>
<span id="cb28-3"><a href="#cb28-3" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-4"><a href="#cb28-4" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Load the Data ---</span></span>
<span id="cb28-5"><a href="#cb28-5" aria-hidden="true" tabindex="-1"></a><span class="co"># Geant4's CSV output begins with metadata lines starting with '#'</span></span>
<span id="cb28-6"><a href="#cb28-6" aria-hidden="true" tabindex="-1"></a><span class="co"># (e.g. "#class tools::wcsv::ntuple", "#column double Edep").</span></span>
<span id="cb28-7"><a href="#cb28-7" aria-hidden="true" tabindex="-1"></a><span class="co"># We tell pandas to skip those and supply the column name ourselves.</span></span>
<span id="cb28-8"><a href="#cb28-8" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb28-9"><a href="#cb28-9" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"GEANT4 SIMULATION DATA ANALYSIS"</span>)</span>
<span id="cb28-10"><a href="#cb28-10" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb28-11"><a href="#cb28-11" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">Reading data from B1_output_nt_B1_t0.csv..."</span>)</span>
<span id="cb28-12"><a href="#cb28-12" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-13"><a href="#cb28-13" aria-hidden="true" tabindex="-1"></a>data <span class="op">=</span> pd.read_csv(<span class="st">'B1_output_nt_B1_t0.csv'</span>,</span>
<span id="cb28-14"><a href="#cb28-14" aria-hidden="true" tabindex="-1"></a> comment<span class="op">=</span><span class="st">'#'</span>, header<span class="op">=</span><span class="va">None</span>, names<span class="op">=</span>[<span class="st">'Edep'</span>])</span>
<span id="cb28-15"><a href="#cb28-15" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-16"><a href="#cb28-16" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Data loaded: </span><span class="sc">{</span><span class="bu">len</span>(data)<span class="sc">}</span><span class="ss"> events, shape </span><span class="sc">{</span>data<span class="sc">.</span>shape<span class="sc">}</span><span class="ss">"</span>)</span>
<span id="cb28-17"><a href="#cb28-17" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">First 10 energy deposition values (MeV):"</span>)</span>
<span id="cb28-18"><a href="#cb28-18" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(data.head(<span class="dv">10</span>))</span>
<span id="cb28-19"><a href="#cb28-19" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-20"><a href="#cb28-20" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Basic Statistics ---</span></span>
<span id="cb28-21"><a href="#cb28-21" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span> <span class="op">+</span> <span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb28-22"><a href="#cb28-22" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"BASIC STATISTICS"</span>)</span>
<span id="cb28-23"><a href="#cb28-23" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb28-24"><a href="#cb28-24" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-25"><a href="#cb28-25" aria-hidden="true" tabindex="-1"></a>mean_edep <span class="op">=</span> data[<span class="st">'Edep'</span>].mean()</span>
<span id="cb28-26"><a href="#cb28-26" aria-hidden="true" tabindex="-1"></a>std_edep <span class="op">=</span> data[<span class="st">'Edep'</span>].std()</span>
<span id="cb28-27"><a href="#cb28-27" aria-hidden="true" tabindex="-1"></a>median_edep <span class="op">=</span> data[<span class="st">'Edep'</span>].median()</span>
<span id="cb28-28"><a href="#cb28-28" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-29"><a href="#cb28-29" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Mean Energy Deposition: </span><span class="sc">{</span>mean_edep<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-30"><a href="#cb28-30" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Standard Deviation: </span><span class="sc">{</span>std_edep<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-31"><a href="#cb28-31" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Median: </span><span class="sc">{</span>median_edep<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-32"><a href="#cb28-32" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Minimum: </span><span class="sc">{</span>data[<span class="st">'Edep'</span>]<span class="sc">.</span><span class="bu">min</span>()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-33"><a href="#cb28-33" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Maximum: </span><span class="sc">{</span>data[<span class="st">'Edep'</span>]<span class="sc">.</span><span class="bu">max</span>()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-34"><a href="#cb28-34" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-35"><a href="#cb28-35" aria-hidden="true" tabindex="-1"></a>percentiles <span class="op">=</span> data[<span class="st">'Edep'</span>].quantile([<span class="fl">0.25</span>, <span class="fl">0.5</span>, <span class="fl">0.75</span>, <span class="fl">0.90</span>, <span class="fl">0.95</span>])</span>
<span id="cb28-36"><a href="#cb28-36" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"</span><span class="ch">\n</span><span class="ss">Percentiles:"</span>)</span>
<span id="cb28-37"><a href="#cb28-37" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" 25th percentile (Q1): </span><span class="sc">{</span>percentiles[<span class="fl">0.25</span>]<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-38"><a href="#cb28-38" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" 50th percentile (median): </span><span class="sc">{</span>percentiles[<span class="fl">0.50</span>]<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-39"><a href="#cb28-39" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" 75th percentile (Q3): </span><span class="sc">{</span>percentiles[<span class="fl">0.75</span>]<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-40"><a href="#cb28-40" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" 90th percentile: </span><span class="sc">{</span>percentiles[<span class="fl">0.90</span>]<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-41"><a href="#cb28-41" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" 95th percentile: </span><span class="sc">{</span>percentiles[<span class="fl">0.95</span>]<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb28-42"><a href="#cb28-42" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-43"><a href="#cb28-43" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Histogram ---</span></span>
<span id="cb28-44"><a href="#cb28-44" aria-hidden="true" tabindex="-1"></a>plt.figure(figsize<span class="op">=</span>(<span class="dv">12</span>, <span class="dv">8</span>))</span>
<span id="cb28-45"><a href="#cb28-45" aria-hidden="true" tabindex="-1"></a>plt.hist(data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">100</span>, alpha<span class="op">=</span><span class="fl">0.75</span>,</span>
<span id="cb28-46"><a href="#cb28-46" aria-hidden="true" tabindex="-1"></a> color<span class="op">=</span><span class="st">'steelblue'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb28-47"><a href="#cb28-47" aria-hidden="true" tabindex="-1"></a>plt.axvline(mean_edep, color<span class="op">=</span><span class="st">'red'</span>, linestyle<span class="op">=</span><span class="st">'--'</span>,</span>
<span id="cb28-48"><a href="#cb28-48" aria-hidden="true" tabindex="-1"></a> linewidth<span class="op">=</span><span class="dv">2</span>, label<span class="op">=</span><span class="ss">f'Mean: </span><span class="sc">{</span>mean_edep<span class="sc">:.3f}</span><span class="ss"> MeV'</span>)</span>
<span id="cb28-49"><a href="#cb28-49" aria-hidden="true" tabindex="-1"></a>plt.axvline(median_edep, color<span class="op">=</span><span class="st">'green'</span>, linestyle<span class="op">=</span><span class="st">'--'</span>,</span>
<span id="cb28-50"><a href="#cb28-50" aria-hidden="true" tabindex="-1"></a> linewidth<span class="op">=</span><span class="dv">2</span>, label<span class="op">=</span><span class="ss">f'Median: </span><span class="sc">{</span>median_edep<span class="sc">:.3f}</span><span class="ss"> MeV'</span>)</span>
<span id="cb28-51"><a href="#cb28-51" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-52"><a href="#cb28-52" aria-hidden="true" tabindex="-1"></a>plt.title(<span class="ss">f'Energy Deposition Distribution for 6 MeV Gammas</span><span class="ch">\n</span><span class="ss">(</span><span class="sc">{</span><span class="bu">len</span>(data)<span class="sc">}</span><span class="ss"> events)'</span>,</span>
<span id="cb28-53"><a href="#cb28-53" aria-hidden="true" tabindex="-1"></a> fontsize<span class="op">=</span><span class="dv">14</span>, fontweight<span class="op">=</span><span class="st">'bold'</span>)</span>
<span id="cb28-54"><a href="#cb28-54" aria-hidden="true" tabindex="-1"></a>plt.xlabel(<span class="st">'Energy Deposited (MeV)'</span>, fontsize<span class="op">=</span><span class="dv">12</span>)</span>
<span id="cb28-55"><a href="#cb28-55" aria-hidden="true" tabindex="-1"></a>plt.ylabel(<span class="st">'Frequency (log scale)'</span>, fontsize<span class="op">=</span><span class="dv">12</span>)</span>
<span id="cb28-56"><a href="#cb28-56" aria-hidden="true" tabindex="-1"></a>plt.yscale(<span class="st">'log'</span>)</span>
<span id="cb28-57"><a href="#cb28-57" aria-hidden="true" tabindex="-1"></a>plt.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb28-58"><a href="#cb28-58" aria-hidden="true" tabindex="-1"></a>plt.legend(fontsize<span class="op">=</span><span class="dv">10</span>)</span>
<span id="cb28-59"><a href="#cb28-59" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-60"><a href="#cb28-60" aria-hidden="true" tabindex="-1"></a>stats_text <span class="op">=</span> <span class="ss">f'mu = </span><span class="sc">{</span>mean_edep<span class="sc">:.3f}</span><span class="ss"> MeV</span><span class="ch">\n</span><span class="ss">sigma = </span><span class="sc">{</span>std_edep<span class="sc">:.3f}</span><span class="ss"> MeV</span><span class="ch">\n</span><span class="ss">Events = </span><span class="sc">{</span><span class="bu">len</span>(data)<span class="sc">}</span><span class="ss">'</span></span>
<span id="cb28-61"><a href="#cb28-61" aria-hidden="true" tabindex="-1"></a>plt.text(<span class="fl">0.98</span>, <span class="fl">0.97</span>, stats_text, transform<span class="op">=</span>plt.gca().transAxes,</span>
<span id="cb28-62"><a href="#cb28-62" aria-hidden="true" tabindex="-1"></a> fontsize<span class="op">=</span><span class="dv">11</span>, verticalalignment<span class="op">=</span><span class="st">'top'</span>, horizontalalignment<span class="op">=</span><span class="st">'right'</span>,</span>
<span id="cb28-63"><a href="#cb28-63" aria-hidden="true" tabindex="-1"></a> bbox<span class="op">=</span><span class="bu">dict</span>(boxstyle<span class="op">=</span><span class="st">'round'</span>, facecolor<span class="op">=</span><span class="st">'wheat'</span>, alpha<span class="op">=</span><span class="fl">0.8</span>))</span>
<span id="cb28-64"><a href="#cb28-64" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb28-65"><a href="#cb28-65" aria-hidden="true" tabindex="-1"></a>plt.tight_layout()</span>
<span id="cb28-66"><a href="#cb28-66" aria-hidden="true" tabindex="-1"></a>plt.savefig(<span class="st">'edep_histogram.png'</span>, dpi<span class="op">=</span><span class="dv">300</span>)</span>
<span id="cb28-67"><a href="#cb28-67" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">Histogram saved to edep_histogram.png"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Save and close, then run:</p>
<div class="sourceCode" id="cb29"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb29-1"><a href="#cb29-1" aria-hidden="true" tabindex="-1"></a><span class="ex">python3</span> analyze_basic.py</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>View the generated histogram:</p>
<div class="sourceCode" id="cb30"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb30-1"><a href="#cb30-1" aria-hidden="true" tabindex="-1"></a><span class="ex">eog</span> edep_histogram.png</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>(Or use your preferred image viewer.)</p>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/histogram_basic.png" class="img-fluid figure-img"></p>
<figcaption>Energy Deposition Histogram</figcaption>
</figure>
</div>
</section>
<section id="step-6-advanced-analysis---data-filtering-and-comparison" class="level2">
<h2 class="anchored" data-anchor-id="step-6-advanced-analysis---data-filtering-and-comparison">Step 6: Advanced Analysis - Data Filtering and Comparison</h2>
<p>Let’s explore pandas’ powerful data manipulation capabilities by filtering the data and creating comparative visualizations.</p>
<p>Create <code>analyze_advanced.py</code>:</p>
<div class="sourceCode" id="cb31"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb31-1"><a href="#cb31-1" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> analyze_advanced.py</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Paste the following content (must start at column 0 — no leading whitespace):</p>
<div class="sourceCode" id="cb32"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb32-1"><a href="#cb32-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> pandas <span class="im">as</span> pd</span>
<span id="cb32-2"><a href="#cb32-2" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> matplotlib.pyplot <span class="im">as</span> plt</span>
<span id="cb32-3"><a href="#cb32-3" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> numpy <span class="im">as</span> np</span>
<span id="cb32-4"><a href="#cb32-4" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-5"><a href="#cb32-5" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"Loading simulation data..."</span>)</span>
<span id="cb32-6"><a href="#cb32-6" aria-hidden="true" tabindex="-1"></a>data <span class="op">=</span> pd.read_csv(<span class="st">'B1_output_nt_B1_t0.csv'</span>,</span>
<span id="cb32-7"><a href="#cb32-7" aria-hidden="true" tabindex="-1"></a> comment<span class="op">=</span><span class="st">'#'</span>, header<span class="op">=</span><span class="va">None</span>, names<span class="op">=</span>[<span class="st">'Edep'</span>])</span>
<span id="cb32-8"><a href="#cb32-8" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-9"><a href="#cb32-9" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Data Filtering Examples ---</span></span>
<span id="cb32-10"><a href="#cb32-10" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span> <span class="op">+</span> <span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-11"><a href="#cb32-11" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"DATA FILTERING EXAMPLES"</span>)</span>
<span id="cb32-12"><a href="#cb32-12" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-13"><a href="#cb32-13" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-14"><a href="#cb32-14" aria-hidden="true" tabindex="-1"></a>high_energy <span class="op">=</span> data[data[<span class="st">'Edep'</span>] <span class="op">&gt;</span> <span class="fl">4.0</span>]</span>
<span id="cb32-15"><a href="#cb32-15" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"</span><span class="ch">\n</span><span class="ss">Events with Edep &gt; 4.0 MeV: </span><span class="sc">{</span><span class="bu">len</span>(high_energy)<span class="sc">}</span><span class="ss"> "</span></span>
<span id="cb32-16"><a href="#cb32-16" aria-hidden="true" tabindex="-1"></a> <span class="ss">f"(</span><span class="sc">{</span><span class="bu">len</span>(high_energy)<span class="op">/</span><span class="bu">len</span>(data)<span class="op">*</span><span class="dv">100</span><span class="sc">:.1f}</span><span class="ss">%)"</span>)</span>
<span id="cb32-17"><a href="#cb32-17" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-18"><a href="#cb32-18" aria-hidden="true" tabindex="-1"></a>low_energy <span class="op">=</span> data[data[<span class="st">'Edep'</span>] <span class="op">&lt;</span> <span class="fl">1.0</span>]</span>
<span id="cb32-19"><a href="#cb32-19" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Events with Edep &lt; 1.0 MeV: </span><span class="sc">{</span><span class="bu">len</span>(low_energy)<span class="sc">}</span><span class="ss"> "</span></span>
<span id="cb32-20"><a href="#cb32-20" aria-hidden="true" tabindex="-1"></a> <span class="ss">f"(</span><span class="sc">{</span><span class="bu">len</span>(low_energy)<span class="op">/</span><span class="bu">len</span>(data)<span class="op">*</span><span class="dv">100</span><span class="sc">:.1f}</span><span class="ss">%)"</span>)</span>
<span id="cb32-21"><a href="#cb32-21" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-22"><a href="#cb32-22" aria-hidden="true" tabindex="-1"></a>zero_energy <span class="op">=</span> data[data[<span class="st">'Edep'</span>] <span class="op">==</span> <span class="fl">0.0</span>]</span>
<span id="cb32-23"><a href="#cb32-23" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Events with Edep = 0.0 MeV: </span><span class="sc">{</span><span class="bu">len</span>(zero_energy)<span class="sc">}</span><span class="ss"> "</span></span>
<span id="cb32-24"><a href="#cb32-24" aria-hidden="true" tabindex="-1"></a> <span class="ss">f"(</span><span class="sc">{</span><span class="bu">len</span>(zero_energy)<span class="op">/</span><span class="bu">len</span>(data)<span class="op">*</span><span class="dv">100</span><span class="sc">:.1f}</span><span class="ss">%)"</span>)</span>
<span id="cb32-25"><a href="#cb32-25" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">" -&gt; These are events where the particle passed through "</span></span>
<span id="cb32-26"><a href="#cb32-26" aria-hidden="true" tabindex="-1"></a> <span class="st">"without interacting"</span>)</span>
<span id="cb32-27"><a href="#cb32-27" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-28"><a href="#cb32-28" aria-hidden="true" tabindex="-1"></a>mid_range <span class="op">=</span> data[(data[<span class="st">'Edep'</span>] <span class="op">&gt;=</span> <span class="fl">2.0</span>) <span class="op">&amp;</span> (data[<span class="st">'Edep'</span>] <span class="op">&lt;=</span> <span class="fl">4.0</span>)]</span>
<span id="cb32-29"><a href="#cb32-29" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f"Events with 2.0 &lt;= Edep &lt;= 4.0 MeV: </span><span class="sc">{</span><span class="bu">len</span>(mid_range)<span class="sc">}</span><span class="ss"> "</span></span>
<span id="cb32-30"><a href="#cb32-30" aria-hidden="true" tabindex="-1"></a> <span class="ss">f"(</span><span class="sc">{</span><span class="bu">len</span>(mid_range)<span class="op">/</span><span class="bu">len</span>(data)<span class="op">*</span><span class="dv">100</span><span class="sc">:.1f}</span><span class="ss">%)"</span>)</span>
<span id="cb32-31"><a href="#cb32-31" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-32"><a href="#cb32-32" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Multi-Panel Figure ---</span></span>
<span id="cb32-33"><a href="#cb32-33" aria-hidden="true" tabindex="-1"></a>fig, axes <span class="op">=</span> plt.subplots(<span class="dv">2</span>, <span class="dv">2</span>, figsize<span class="op">=</span>(<span class="dv">14</span>, <span class="dv">10</span>))</span>
<span id="cb32-34"><a href="#cb32-34" aria-hidden="true" tabindex="-1"></a>fig.suptitle(<span class="st">'Advanced Energy Deposition Analysis'</span>,</span>
<span id="cb32-35"><a href="#cb32-35" aria-hidden="true" tabindex="-1"></a> fontsize<span class="op">=</span><span class="dv">16</span>, fontweight<span class="op">=</span><span class="st">'bold'</span>)</span>
<span id="cb32-36"><a href="#cb32-36" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-37"><a href="#cb32-37" aria-hidden="true" tabindex="-1"></a><span class="co"># Panel 1: Full distribution (log y to see the tail past the zero spike)</span></span>
<span id="cb32-38"><a href="#cb32-38" aria-hidden="true" tabindex="-1"></a>ax1 <span class="op">=</span> axes[<span class="dv">0</span>, <span class="dv">0</span>]</span>
<span id="cb32-39"><a href="#cb32-39" aria-hidden="true" tabindex="-1"></a>ax1.hist(data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">100</span>, color<span class="op">=</span><span class="st">'steelblue'</span>,</span>
<span id="cb32-40"><a href="#cb32-40" aria-hidden="true" tabindex="-1"></a> alpha<span class="op">=</span><span class="fl">0.7</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb32-41"><a href="#cb32-41" aria-hidden="true" tabindex="-1"></a>ax1.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb32-42"><a href="#cb32-42" aria-hidden="true" tabindex="-1"></a>ax1.set_ylabel(<span class="st">'Frequency (log scale)'</span>)</span>
<span id="cb32-43"><a href="#cb32-43" aria-hidden="true" tabindex="-1"></a>ax1.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb32-44"><a href="#cb32-44" aria-hidden="true" tabindex="-1"></a>ax1.set_title(<span class="ss">f'Full Distribution (n=</span><span class="sc">{</span><span class="bu">len</span>(data)<span class="sc">}</span><span class="ss">)'</span>)</span>
<span id="cb32-45"><a href="#cb32-45" aria-hidden="true" tabindex="-1"></a>ax1.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb32-46"><a href="#cb32-46" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-47"><a href="#cb32-47" aria-hidden="true" tabindex="-1"></a><span class="co"># Panel 2: Excluding zero-energy events (log y for the long tail)</span></span>
<span id="cb32-48"><a href="#cb32-48" aria-hidden="true" tabindex="-1"></a>ax2 <span class="op">=</span> axes[<span class="dv">0</span>, <span class="dv">1</span>]</span>
<span id="cb32-49"><a href="#cb32-49" aria-hidden="true" tabindex="-1"></a>non_zero <span class="op">=</span> data[data[<span class="st">'Edep'</span>] <span class="op">&gt;</span> <span class="dv">0</span>]</span>
<span id="cb32-50"><a href="#cb32-50" aria-hidden="true" tabindex="-1"></a>ax2.hist(non_zero[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">80</span>, color<span class="op">=</span><span class="st">'coral'</span>,</span>
<span id="cb32-51"><a href="#cb32-51" aria-hidden="true" tabindex="-1"></a> alpha<span class="op">=</span><span class="fl">0.7</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb32-52"><a href="#cb32-52" aria-hidden="true" tabindex="-1"></a>ax2.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb32-53"><a href="#cb32-53" aria-hidden="true" tabindex="-1"></a>ax2.set_ylabel(<span class="st">'Frequency (log scale)'</span>)</span>
<span id="cb32-54"><a href="#cb32-54" aria-hidden="true" tabindex="-1"></a>ax2.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb32-55"><a href="#cb32-55" aria-hidden="true" tabindex="-1"></a>ax2.set_title(<span class="ss">f'Non-Zero Events Only (n=</span><span class="sc">{</span><span class="bu">len</span>(non_zero)<span class="sc">}</span><span class="ss">)'</span>)</span>
<span id="cb32-56"><a href="#cb32-56" aria-hidden="true" tabindex="-1"></a>ax2.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb32-57"><a href="#cb32-57" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-58"><a href="#cb32-58" aria-hidden="true" tabindex="-1"></a><span class="co"># Panel 3: Cumulative distribution</span></span>
<span id="cb32-59"><a href="#cb32-59" aria-hidden="true" tabindex="-1"></a>ax3 <span class="op">=</span> axes[<span class="dv">1</span>, <span class="dv">0</span>]</span>
<span id="cb32-60"><a href="#cb32-60" aria-hidden="true" tabindex="-1"></a>sorted_data <span class="op">=</span> np.sort(data[<span class="st">'Edep'</span>])</span>
<span id="cb32-61"><a href="#cb32-61" aria-hidden="true" tabindex="-1"></a>cumulative <span class="op">=</span> np.arange(<span class="dv">1</span>, <span class="bu">len</span>(sorted_data) <span class="op">+</span> <span class="dv">1</span>) <span class="op">/</span> <span class="bu">len</span>(sorted_data) <span class="op">*</span> <span class="dv">100</span></span>
<span id="cb32-62"><a href="#cb32-62" aria-hidden="true" tabindex="-1"></a>ax3.plot(sorted_data, cumulative, linewidth<span class="op">=</span><span class="dv">2</span>, color<span class="op">=</span><span class="st">'darkgreen'</span>)</span>
<span id="cb32-63"><a href="#cb32-63" aria-hidden="true" tabindex="-1"></a>ax3.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb32-64"><a href="#cb32-64" aria-hidden="true" tabindex="-1"></a>ax3.set_ylabel(<span class="st">'Cumulative Percentage (%)'</span>)</span>
<span id="cb32-65"><a href="#cb32-65" aria-hidden="true" tabindex="-1"></a>ax3.set_title(<span class="st">'Cumulative Distribution Function'</span>)</span>
<span id="cb32-66"><a href="#cb32-66" aria-hidden="true" tabindex="-1"></a>ax3.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>)</span>
<span id="cb32-67"><a href="#cb32-67" aria-hidden="true" tabindex="-1"></a>ax3.axhline(<span class="dv">50</span>, color<span class="op">=</span><span class="st">'red'</span>, linestyle<span class="op">=</span><span class="st">'--'</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, label<span class="op">=</span><span class="st">'50%'</span>)</span>
<span id="cb32-68"><a href="#cb32-68" aria-hidden="true" tabindex="-1"></a>ax3.axhline(<span class="dv">95</span>, color<span class="op">=</span><span class="st">'orange'</span>, linestyle<span class="op">=</span><span class="st">'--'</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, label<span class="op">=</span><span class="st">'95%'</span>)</span>
<span id="cb32-69"><a href="#cb32-69" aria-hidden="true" tabindex="-1"></a>ax3.legend()</span>
<span id="cb32-70"><a href="#cb32-70" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-71"><a href="#cb32-71" aria-hidden="true" tabindex="-1"></a><span class="co"># Panel 4: Box plot</span></span>
<span id="cb32-72"><a href="#cb32-72" aria-hidden="true" tabindex="-1"></a>ax4 <span class="op">=</span> axes[<span class="dv">1</span>, <span class="dv">1</span>]</span>
<span id="cb32-73"><a href="#cb32-73" aria-hidden="true" tabindex="-1"></a>box_data <span class="op">=</span> [data[<span class="st">'Edep'</span>], non_zero[<span class="st">'Edep'</span>]]</span>
<span id="cb32-74"><a href="#cb32-74" aria-hidden="true" tabindex="-1"></a>bp <span class="op">=</span> ax4.boxplot(box_data, labels<span class="op">=</span>[<span class="st">'All Events'</span>, <span class="st">'Non-Zero Only'</span>],</span>
<span id="cb32-75"><a href="#cb32-75" aria-hidden="true" tabindex="-1"></a> patch_artist<span class="op">=</span><span class="va">True</span>, showmeans<span class="op">=</span><span class="va">True</span>)</span>
<span id="cb32-76"><a href="#cb32-76" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> patch, color <span class="kw">in</span> <span class="bu">zip</span>(bp[<span class="st">'boxes'</span>], [<span class="st">'lightblue'</span>, <span class="st">'lightcoral'</span>]):</span>
<span id="cb32-77"><a href="#cb32-77" aria-hidden="true" tabindex="-1"></a> patch.set_facecolor(color)</span>
<span id="cb32-78"><a href="#cb32-78" aria-hidden="true" tabindex="-1"></a>ax4.set_ylabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb32-79"><a href="#cb32-79" aria-hidden="true" tabindex="-1"></a>ax4.set_title(<span class="st">'Box Plot Comparison'</span>)</span>
<span id="cb32-80"><a href="#cb32-80" aria-hidden="true" tabindex="-1"></a>ax4.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, axis<span class="op">=</span><span class="st">'y'</span>)</span>
<span id="cb32-81"><a href="#cb32-81" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-82"><a href="#cb32-82" aria-hidden="true" tabindex="-1"></a>plt.tight_layout()</span>
<span id="cb32-83"><a href="#cb32-83" aria-hidden="true" tabindex="-1"></a>plt.savefig(<span class="st">'edep_advanced_analysis.png'</span>, dpi<span class="op">=</span><span class="dv">300</span>)</span>
<span id="cb32-84"><a href="#cb32-84" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">Advanced analysis figure saved to edep_advanced_analysis.png"</span>)</span>
<span id="cb32-85"><a href="#cb32-85" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-86"><a href="#cb32-86" aria-hidden="true" tabindex="-1"></a><span class="co"># --- Statistical Comparison ---</span></span>
<span id="cb32-87"><a href="#cb32-87" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span> <span class="op">+</span> <span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-88"><a href="#cb32-88" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"STATISTICAL COMPARISON"</span>)</span>
<span id="cb32-89"><a href="#cb32-89" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-90"><a href="#cb32-90" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">All Events:"</span>)</span>
<span id="cb32-91"><a href="#cb32-91" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Mean: </span><span class="sc">{</span>data[<span class="st">'Edep'</span>]<span class="sc">.</span>mean()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-92"><a href="#cb32-92" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Median: </span><span class="sc">{</span>data[<span class="st">'Edep'</span>]<span class="sc">.</span>median()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-93"><a href="#cb32-93" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Std: </span><span class="sc">{</span>data[<span class="st">'Edep'</span>]<span class="sc">.</span>std()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-94"><a href="#cb32-94" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">Non-Zero Events Only:"</span>)</span>
<span id="cb32-95"><a href="#cb32-95" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Mean: </span><span class="sc">{</span>non_zero[<span class="st">'Edep'</span>]<span class="sc">.</span>mean()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-96"><a href="#cb32-96" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Median: </span><span class="sc">{</span>non_zero[<span class="st">'Edep'</span>]<span class="sc">.</span>median()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-97"><a href="#cb32-97" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="ss">f" Std: </span><span class="sc">{</span>non_zero[<span class="st">'Edep'</span>]<span class="sc">.</span>std()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb32-98"><a href="#cb32-98" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb32-99"><a href="#cb32-99" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span> <span class="op">+</span> <span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-100"><a href="#cb32-100" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"DETAILED SUMMARY STATISTICS"</span>)</span>
<span id="cb32-101"><a href="#cb32-101" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">50</span>)</span>
<span id="cb32-102"><a href="#cb32-102" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(data[<span class="st">'Edep'</span>].describe())</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Run it and view the figure:</p>
<div class="sourceCode" id="cb33"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb33-1"><a href="#cb33-1" aria-hidden="true" tabindex="-1"></a><span class="ex">python3</span> analyze_advanced.py</span>
<span id="cb33-2"><a href="#cb33-2" aria-hidden="true" tabindex="-1"></a><span class="ex">eog</span> edep_advanced_analysis.png</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/advanced_analysis.png" class="img-fluid figure-img"></p>
<figcaption>Advanced Analysis Multi-Panel</figcaption>
</figure>
</div>
</section>
<section id="step-7-comparing-different-particle-types" class="level2">
<h2 class="anchored" data-anchor-id="step-7-comparing-different-particle-types">Step 7: Comparing Different Particle Types</h2>
<p>Let’s generate data for different particles and compare their energy deposition patterns.</p>
<ol type="1">
<li><p>Create macro files for different particles. First, protons:</p>
<div class="sourceCode" id="cb34"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb34-1"><a href="#cb34-1" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> ~/g4-labs/lab3/B1</span>
<span id="cb34-2"><a href="#cb34-2" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> proton_beam.mac</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Add:</p>
<div class="sourceCode" id="cb35"><pre class="sourceCode default code-with-copy"><code class="sourceCode default"><span id="cb35-1"><a href="#cb35-1" aria-hidden="true" tabindex="-1"></a>/run/numberOfThreads 1</span>
<span id="cb35-2"><a href="#cb35-2" aria-hidden="true" tabindex="-1"></a>/run/initialize</span>
<span id="cb35-3"><a href="#cb35-3" aria-hidden="true" tabindex="-1"></a>/gun/particle proton</span>
<span id="cb35-4"><a href="#cb35-4" aria-hidden="true" tabindex="-1"></a>/gun/energy 150 MeV</span>
<span id="cb35-5"><a href="#cb35-5" aria-hidden="true" tabindex="-1"></a>/run/printProgress 1000</span>
<span id="cb35-6"><a href="#cb35-6" aria-hidden="true" tabindex="-1"></a>/run/beamOn 10000</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Create electron beam macro:</p>
<div class="sourceCode" id="cb36"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb36-1"><a href="#cb36-1" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> electron_beam.mac</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Add:</p>
<div class="sourceCode" id="cb37"><pre class="sourceCode default code-with-copy"><code class="sourceCode default"><span id="cb37-1"><a href="#cb37-1" aria-hidden="true" tabindex="-1"></a>/run/numberOfThreads 1</span>
<span id="cb37-2"><a href="#cb37-2" aria-hidden="true" tabindex="-1"></a>/run/initialize</span>
<span id="cb37-3"><a href="#cb37-3" aria-hidden="true" tabindex="-1"></a>/gun/particle e-</span>
<span id="cb37-4"><a href="#cb37-4" aria-hidden="true" tabindex="-1"></a>/gun/energy 50 MeV</span>
<span id="cb37-5"><a href="#cb37-5" aria-hidden="true" tabindex="-1"></a>/run/printProgress 1000</span>
<span id="cb37-6"><a href="#cb37-6" aria-hidden="true" tabindex="-1"></a>/run/beamOn 10000</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
<li><p>Run simulations and rename output files. Remember: each run produces <code>B1_output_nt_B1_t0.csv</code>, so we rename it after each run before it gets overwritten.</p>
<div class="sourceCode" id="cb38"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb38-1"><a href="#cb38-1" aria-hidden="true" tabindex="-1"></a><span class="bu">cd</span> build</span>
<span id="cb38-2"><a href="#cb38-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb38-3"><a href="#cb38-3" aria-hidden="true" tabindex="-1"></a><span class="co"># Run gamma simulation (already done, but let's be explicit)</span></span>
<span id="cb38-4"><a href="#cb38-4" aria-hidden="true" tabindex="-1"></a><span class="ex">./exampleB1</span> ../lab3.mac</span>
<span id="cb38-5"><a href="#cb38-5" aria-hidden="true" tabindex="-1"></a><span class="fu">mv</span> B1_output_nt_B1_t0.csv gamma_6MeV.csv</span>
<span id="cb38-6"><a href="#cb38-6" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb38-7"><a href="#cb38-7" aria-hidden="true" tabindex="-1"></a><span class="co"># Run proton simulation</span></span>
<span id="cb38-8"><a href="#cb38-8" aria-hidden="true" tabindex="-1"></a><span class="ex">./exampleB1</span> ../proton_beam.mac</span>
<span id="cb38-9"><a href="#cb38-9" aria-hidden="true" tabindex="-1"></a><span class="fu">mv</span> B1_output_nt_B1_t0.csv proton_150MeV.csv</span>
<span id="cb38-10"><a href="#cb38-10" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb38-11"><a href="#cb38-11" aria-hidden="true" tabindex="-1"></a><span class="co"># Run electron simulation</span></span>
<span id="cb38-12"><a href="#cb38-12" aria-hidden="true" tabindex="-1"></a><span class="ex">./exampleB1</span> ../electron_beam.mac</span>
<span id="cb38-13"><a href="#cb38-13" aria-hidden="true" tabindex="-1"></a><span class="fu">mv</span> B1_output_nt_B1_t0.csv electron_50MeV.csv</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div></li>
</ol>
<p>Create the comparison script:</p>
<div class="sourceCode" id="cb39"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb39-1"><a href="#cb39-1" aria-hidden="true" tabindex="-1"></a><span class="fu">nano</span> compare_particles.py</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Paste the following content (must start at column 0 — no leading whitespace):</p>
<div class="sourceCode" id="cb40"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb40-1"><a href="#cb40-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> pandas <span class="im">as</span> pd</span>
<span id="cb40-2"><a href="#cb40-2" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> matplotlib.pyplot <span class="im">as</span> plt</span>
<span id="cb40-3"><a href="#cb40-3" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-4"><a href="#cb40-4" aria-hidden="true" tabindex="-1"></a><span class="co"># All three CSVs were produced by the modified B1, so they share the</span></span>
<span id="cb40-5"><a href="#cb40-5" aria-hidden="true" tabindex="-1"></a><span class="co"># same Geant4 header format. We supply the column name ourselves.</span></span>
<span id="cb40-6"><a href="#cb40-6" aria-hidden="true" tabindex="-1"></a>read_kwargs <span class="op">=</span> <span class="bu">dict</span>(comment<span class="op">=</span><span class="st">'#'</span>, header<span class="op">=</span><span class="va">None</span>, names<span class="op">=</span>[<span class="st">'Edep'</span>])</span>
<span id="cb40-7"><a href="#cb40-7" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-8"><a href="#cb40-8" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"Loading particle data..."</span>)</span>
<span id="cb40-9"><a href="#cb40-9" aria-hidden="true" tabindex="-1"></a>gamma_data <span class="op">=</span> pd.read_csv(<span class="st">'gamma_6MeV.csv'</span>, <span class="op">**</span>read_kwargs)</span>
<span id="cb40-10"><a href="#cb40-10" aria-hidden="true" tabindex="-1"></a>proton_data <span class="op">=</span> pd.read_csv(<span class="st">'proton_150MeV.csv'</span>, <span class="op">**</span>read_kwargs)</span>
<span id="cb40-11"><a href="#cb40-11" aria-hidden="true" tabindex="-1"></a>electron_data <span class="op">=</span> pd.read_csv(<span class="st">'electron_50MeV.csv'</span>, <span class="op">**</span>read_kwargs)</span>
<span id="cb40-12"><a href="#cb40-12" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-13"><a href="#cb40-13" aria-hidden="true" tabindex="-1"></a>fig, axes <span class="op">=</span> plt.subplots(<span class="dv">2</span>, <span class="dv">2</span>, figsize<span class="op">=</span>(<span class="dv">14</span>, <span class="dv">10</span>))</span>
<span id="cb40-14"><a href="#cb40-14" aria-hidden="true" tabindex="-1"></a>fig.suptitle(<span class="st">'Particle Type Comparison: Energy Deposition'</span>,</span>
<span id="cb40-15"><a href="#cb40-15" aria-hidden="true" tabindex="-1"></a> fontsize<span class="op">=</span><span class="dv">16</span>, fontweight<span class="op">=</span><span class="st">'bold'</span>)</span>
<span id="cb40-16"><a href="#cb40-16" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-17"><a href="#cb40-17" aria-hidden="true" tabindex="-1"></a>ax1 <span class="op">=</span> axes[<span class="dv">0</span>, <span class="dv">0</span>]</span>
<span id="cb40-18"><a href="#cb40-18" aria-hidden="true" tabindex="-1"></a>ax1.hist(gamma_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">80</span>, alpha<span class="op">=</span><span class="fl">0.7</span>,</span>
<span id="cb40-19"><a href="#cb40-19" aria-hidden="true" tabindex="-1"></a> color<span class="op">=</span><span class="st">'blue'</span>, label<span class="op">=</span><span class="st">'6 MeV gamma'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb40-20"><a href="#cb40-20" aria-hidden="true" tabindex="-1"></a>ax1.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb40-21"><a href="#cb40-21" aria-hidden="true" tabindex="-1"></a>ax1.set_ylabel(<span class="st">'Frequency (log scale)'</span>)</span>
<span id="cb40-22"><a href="#cb40-22" aria-hidden="true" tabindex="-1"></a>ax1.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb40-23"><a href="#cb40-23" aria-hidden="true" tabindex="-1"></a>ax1.set_title(<span class="st">'Gamma Rays (6 MeV)'</span>)</span>
<span id="cb40-24"><a href="#cb40-24" aria-hidden="true" tabindex="-1"></a>ax1.legend()</span>
<span id="cb40-25"><a href="#cb40-25" aria-hidden="true" tabindex="-1"></a>ax1.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb40-26"><a href="#cb40-26" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-27"><a href="#cb40-27" aria-hidden="true" tabindex="-1"></a>ax2 <span class="op">=</span> axes[<span class="dv">0</span>, <span class="dv">1</span>]</span>
<span id="cb40-28"><a href="#cb40-28" aria-hidden="true" tabindex="-1"></a>ax2.hist(proton_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">80</span>, alpha<span class="op">=</span><span class="fl">0.7</span>,</span>
<span id="cb40-29"><a href="#cb40-29" aria-hidden="true" tabindex="-1"></a> color<span class="op">=</span><span class="st">'red'</span>, label<span class="op">=</span><span class="st">'150 MeV p+'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb40-30"><a href="#cb40-30" aria-hidden="true" tabindex="-1"></a>ax2.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb40-31"><a href="#cb40-31" aria-hidden="true" tabindex="-1"></a>ax2.set_ylabel(<span class="st">'Frequency (log scale)'</span>)</span>
<span id="cb40-32"><a href="#cb40-32" aria-hidden="true" tabindex="-1"></a>ax2.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb40-33"><a href="#cb40-33" aria-hidden="true" tabindex="-1"></a>ax2.set_title(<span class="st">'Protons (150 MeV)'</span>)</span>
<span id="cb40-34"><a href="#cb40-34" aria-hidden="true" tabindex="-1"></a>ax2.legend()</span>
<span id="cb40-35"><a href="#cb40-35" aria-hidden="true" tabindex="-1"></a>ax2.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb40-36"><a href="#cb40-36" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-37"><a href="#cb40-37" aria-hidden="true" tabindex="-1"></a>ax3 <span class="op">=</span> axes[<span class="dv">1</span>, <span class="dv">0</span>]</span>
<span id="cb40-38"><a href="#cb40-38" aria-hidden="true" tabindex="-1"></a>ax3.hist(electron_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">80</span>, alpha<span class="op">=</span><span class="fl">0.7</span>,</span>
<span id="cb40-39"><a href="#cb40-39" aria-hidden="true" tabindex="-1"></a> color<span class="op">=</span><span class="st">'green'</span>, label<span class="op">=</span><span class="st">'50 MeV e-'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>)</span>
<span id="cb40-40"><a href="#cb40-40" aria-hidden="true" tabindex="-1"></a>ax3.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb40-41"><a href="#cb40-41" aria-hidden="true" tabindex="-1"></a>ax3.set_ylabel(<span class="st">'Frequency (log scale)'</span>)</span>
<span id="cb40-42"><a href="#cb40-42" aria-hidden="true" tabindex="-1"></a>ax3.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb40-43"><a href="#cb40-43" aria-hidden="true" tabindex="-1"></a>ax3.set_title(<span class="st">'Electrons (50 MeV)'</span>)</span>
<span id="cb40-44"><a href="#cb40-44" aria-hidden="true" tabindex="-1"></a>ax3.legend()</span>
<span id="cb40-45"><a href="#cb40-45" aria-hidden="true" tabindex="-1"></a>ax3.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb40-46"><a href="#cb40-46" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-47"><a href="#cb40-47" aria-hidden="true" tabindex="-1"></a><span class="co"># Normalized overlay (log y so the three shapes are comparable across orders)</span></span>
<span id="cb40-48"><a href="#cb40-48" aria-hidden="true" tabindex="-1"></a>ax4 <span class="op">=</span> axes[<span class="dv">1</span>, <span class="dv">1</span>]</span>
<span id="cb40-49"><a href="#cb40-49" aria-hidden="true" tabindex="-1"></a>ax4.hist(gamma_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">60</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, color<span class="op">=</span><span class="st">'blue'</span>,</span>
<span id="cb40-50"><a href="#cb40-50" aria-hidden="true" tabindex="-1"></a> label<span class="op">=</span><span class="ss">f'gamma (mu=</span><span class="sc">{</span>gamma_data[<span class="st">"Edep"</span>]<span class="sc">.</span>mean()<span class="sc">:.2f}</span><span class="ss">)'</span>, density<span class="op">=</span><span class="va">True</span>)</span>
<span id="cb40-51"><a href="#cb40-51" aria-hidden="true" tabindex="-1"></a>ax4.hist(proton_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">60</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, color<span class="op">=</span><span class="st">'red'</span>,</span>
<span id="cb40-52"><a href="#cb40-52" aria-hidden="true" tabindex="-1"></a> label<span class="op">=</span><span class="ss">f'p+ (mu=</span><span class="sc">{</span>proton_data[<span class="st">"Edep"</span>]<span class="sc">.</span>mean()<span class="sc">:.2f}</span><span class="ss">)'</span>, density<span class="op">=</span><span class="va">True</span>)</span>
<span id="cb40-53"><a href="#cb40-53" aria-hidden="true" tabindex="-1"></a>ax4.hist(electron_data[<span class="st">'Edep'</span>], bins<span class="op">=</span><span class="dv">60</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, color<span class="op">=</span><span class="st">'green'</span>,</span>
<span id="cb40-54"><a href="#cb40-54" aria-hidden="true" tabindex="-1"></a> label<span class="op">=</span><span class="ss">f'e- (mu=</span><span class="sc">{</span>electron_data[<span class="st">"Edep"</span>]<span class="sc">.</span>mean()<span class="sc">:.2f}</span><span class="ss">)'</span>, density<span class="op">=</span><span class="va">True</span>)</span>
<span id="cb40-55"><a href="#cb40-55" aria-hidden="true" tabindex="-1"></a>ax4.set_xlabel(<span class="st">'Energy Deposited (MeV)'</span>)</span>
<span id="cb40-56"><a href="#cb40-56" aria-hidden="true" tabindex="-1"></a>ax4.set_ylabel(<span class="st">'Probability Density (log scale)'</span>)</span>
<span id="cb40-57"><a href="#cb40-57" aria-hidden="true" tabindex="-1"></a>ax4.set_yscale(<span class="st">'log'</span>)</span>
<span id="cb40-58"><a href="#cb40-58" aria-hidden="true" tabindex="-1"></a>ax4.set_title(<span class="st">'Normalized Overlay Comparison'</span>)</span>
<span id="cb40-59"><a href="#cb40-59" aria-hidden="true" tabindex="-1"></a>ax4.legend()</span>
<span id="cb40-60"><a href="#cb40-60" aria-hidden="true" tabindex="-1"></a>ax4.grid(<span class="va">True</span>, alpha<span class="op">=</span><span class="fl">0.3</span>, which<span class="op">=</span><span class="st">'both'</span>)</span>
<span id="cb40-61"><a href="#cb40-61" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-62"><a href="#cb40-62" aria-hidden="true" tabindex="-1"></a>plt.tight_layout()</span>
<span id="cb40-63"><a href="#cb40-63" aria-hidden="true" tabindex="-1"></a>plt.savefig(<span class="st">'particle_comparison.png'</span>, dpi<span class="op">=</span><span class="dv">300</span>)</span>
<span id="cb40-64"><a href="#cb40-64" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"Comparison figure saved to particle_comparison.png"</span>)</span>
<span id="cb40-65"><a href="#cb40-65" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-66"><a href="#cb40-66" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span> <span class="op">+</span> <span class="st">"="</span> <span class="op">*</span> <span class="dv">60</span>)</span>
<span id="cb40-67"><a href="#cb40-67" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"STATISTICAL COMPARISON BY PARTICLE TYPE"</span>)</span>
<span id="cb40-68"><a href="#cb40-68" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(<span class="st">"="</span> <span class="op">*</span> <span class="dv">60</span>)</span>
<span id="cb40-69"><a href="#cb40-69" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb40-70"><a href="#cb40-70" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> name, df <span class="kw">in</span> [(<span class="st">'Gamma (6 MeV)'</span>, gamma_data),</span>
<span id="cb40-71"><a href="#cb40-71" aria-hidden="true" tabindex="-1"></a> (<span class="st">'Proton (150 MeV)'</span>, proton_data),</span>
<span id="cb40-72"><a href="#cb40-72" aria-hidden="true" tabindex="-1"></a> (<span class="st">'Electron (50 MeV)'</span>, electron_data)]:</span>
<span id="cb40-73"><a href="#cb40-73" aria-hidden="true" tabindex="-1"></a> n_zero <span class="op">=</span> (df[<span class="st">'Edep'</span>] <span class="op">==</span> <span class="dv">0</span>).<span class="bu">sum</span>()</span>
<span id="cb40-74"><a href="#cb40-74" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f"</span><span class="ch">\n</span><span class="sc">{</span>name<span class="sc">}</span><span class="ss">:"</span>)</span>
<span id="cb40-75"><a href="#cb40-75" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f" Mean: </span><span class="sc">{</span>df[<span class="st">'Edep'</span>]<span class="sc">.</span>mean()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb40-76"><a href="#cb40-76" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f" Median: </span><span class="sc">{</span>df[<span class="st">'Edep'</span>]<span class="sc">.</span>median()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb40-77"><a href="#cb40-77" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f" Std Dev: </span><span class="sc">{</span>df[<span class="st">'Edep'</span>]<span class="sc">.</span>std()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb40-78"><a href="#cb40-78" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f" Max: </span><span class="sc">{</span>df[<span class="st">'Edep'</span>]<span class="sc">.</span><span class="bu">max</span>()<span class="sc">:.4f}</span><span class="ss"> MeV"</span>)</span>
<span id="cb40-79"><a href="#cb40-79" aria-hidden="true" tabindex="-1"></a> <span class="bu">print</span>(<span class="ss">f" Zero events: </span><span class="sc">{</span>n_zero<span class="sc">}</span><span class="ss"> (</span><span class="sc">{</span>n_zero<span class="op">/</span><span class="bu">len</span>(df)<span class="op">*</span><span class="dv">100</span><span class="sc">:.1f}</span><span class="ss">%)"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<p>Run it:</p>
<div class="sourceCode" id="cb41"><pre class="sourceCode bash code-with-copy"><code class="sourceCode bash"><span id="cb41-1"><a href="#cb41-1" aria-hidden="true" tabindex="-1"></a><span class="ex">python3</span> compare_particles.py</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="quarto-figure quarto-figure-center">
<figure class="figure">
<p><img src="img/particle_comparison.png" class="img-fluid figure-img"></p>
<figcaption>Particle Type Comparison</figcaption>
</figure>
</div>
</section>
<section id="analysis-questions" class="level2">
<h2 class="anchored" data-anchor-id="analysis-questions">Analysis Questions</h2>
<p>Now that you have generated and analyzed the data, answer these questions:</p>
<section id="part-1-understanding-the-distribution" class="level3">
<h3 class="anchored" data-anchor-id="part-1-understanding-the-distribution">Part 1: Understanding the Distribution</h3>
<ol type="1">
<li><p>Describe the shape of the energy deposition histogram for gamma rays. Is the energy deposition the same for every event?</p></li>
<li><p>Why is there a distribution of values instead of a single value? Consider:</p>
<ul>
<li>Stochastic nature of particle interactions</li>
<li>Different interaction mechanisms (photoelectric, Compton, pair production)</li>
<li>Variation in secondary particle production</li>
</ul></li>
<li><p>The peak of the distribution represents the most probable energy deposition. What is this value approximately? How does it compare to the initial energy of the gamma particle (6 MeV)?</p></li>
<li><p>Explain why the deposited energy is typically lower than the incident energy. (Hint: Does the particle always deposit all of its energy in the detector?)</p></li>
</ol>
</section>
<section id="part-2-understanding-zero-energy-events" class="level3">
<h3 class="anchored" data-anchor-id="part-2-understanding-zero-energy-events">Part 2: Understanding Zero Energy Events</h3>
<ol start="5" type="1">
<li><p>What percentage of events have zero energy deposition? What does this represent physically?</p></li>
<li><p>Why would a 6 MeV gamma ray pass through the detector without depositing any energy?</p></li>
</ol>
</section>
<section id="part-3-comparing-particle-types" class="level3">
<h3 class="anchored" data-anchor-id="part-3-comparing-particle-types">Part 3: Comparing Particle Types</h3>
<ol start="7" type="1">
<li><p>Compare the mean energy deposition for:</p>
<ul>
<li>6 MeV gammas</li>
<li>150 MeV protons<br>
</li>
<li>50 MeV electrons</li>
</ul></li>
<li><p>Which particle type deposits the most energy on average? Why?</p></li>
<li><p>Which particle type has the most consistent energy deposition (smallest standard deviation)? Explain this in terms of the physics of how each particle type interacts with matter.</p></li>
<li><p>Notice the shape differences in the distributions:</p>
<ul>
<li>Gammas: Broad, exponential-like distribution</li>
<li>Protons: Sharp, peaked distribution</li>
<li>Electrons: Intermediate characteristics</li>
</ul>
<p>Explain these shape differences based on the interaction mechanisms of each particle type.</p></li>
</ol>
</section>
<section id="part-4-pandas-data-analysis" class="level3">
<h3 class="anchored" data-anchor-id="part-4-pandas-data-analysis">Part 4: Pandas Data Analysis</h3>
<ol start="11" type="1">
<li><p>In the advanced analysis script, we filtered events with <code>data[data['Edep'] &gt; 4.0]</code>. Explain in plain English what this pandas operation does.</p></li>
<li><p>What is the advantage of using pandas DataFrames versus manipulating the CSV file directly with basic Python?</p></li>
<li><p>The cumulative distribution plot shows what percentage of events have energy deposition below a certain value. What is the energy value at which 50% of events fall below? What is this called statistically?</p></li>
</ol>
</section>
<section id="part-5-scientific-workflow" class="level3">
<h3 class="anchored" data-anchor-id="part-5-scientific-workflow">Part 5: Scientific Workflow</h3>
<ol start="14" type="1">
<li><p>Why is CSV a good intermediate format for this analysis workflow (Geant4 → CSV → Python)?</p></li>
<li><p>How would you modify the code to save additional information, such as:</p>
<ul>
<li>The number of steps per event</li>
<li>Energy deposited in different detector volumes</li>
<li>Particle type information for secondary particles</li>
</ul></li>
</ol>
</section>
</section>
<section id="troubleshooting" class="level2">
<h2 class="anchored" data-anchor-id="troubleshooting">Troubleshooting</h2>
<p>Common issues and solutions:</p>
<ul>
<li><strong>“No module named ‘pandas’”</strong>: Install with <code>pip3 install --break-system-packages pandas matplotlib numpy</code>.</li>
<li><strong><code>error: externally-managed-environment</code> from pip</strong>: Ubuntu 24.04+ marks the system Python as managed (PEP 668). Re-run with <code>--break-system-packages</code>, as shown in Step 4.</li>
<li><strong>Many <code>B1_output_nt_B1_tN.csv</code> files instead of one <code>_t0.csv</code></strong>: Geant4 ran multi-threaded and each worker thread wrote its own CSV. Add <code>/run/numberOfThreads 1</code> to your macro <strong>before</strong> <code>/run/initialize</code>, delete the per-thread files (<code>rm B1_output_nt_B1_t*.csv</code>), and re-run.</li>
<li><strong>CSV file not found</strong>: Make sure you’re running Python from the <code>build</code> directory and that the filename in <code>pd.read_csv(...)</code> matches what <code>ls *.csv</code> shows.</li>
<li><strong>Empty histogram</strong>: Check that your CSV file has data and proper column headers</li>
<li><strong>Compilation errors</strong>: Verify all code modifications carefully match the instructions</li>
<li><strong>G4AnalysisManager errors</strong>: Ensure you’ve included the header file in all necessary source files</li>
</ul>
</section>
<section id="verification-checklist" class="level2">
<h2 class="anchored" data-anchor-id="verification-checklist">Verification Checklist</h2>
<p>Before completing this lab, verify you can:</p>
<ul class="task-list">
<li><label><input type="checkbox">Modify Geant4 C++ source code to add data output</label></li>
<li><label><input type="checkbox">Use G4AnalysisManager to create and write CSV files</label></li>
<li><label><input type="checkbox">Compile modified Geant4 applications</label></li>
<li><label><input type="checkbox">Generate simulation data with multiple particle types</label></li>
<li><label><input type="checkbox">Load CSV data into pandas DataFrames</label></li>
<li><label><input type="checkbox">Calculate statistical measures (mean, std, median, percentiles)</label></li>
<li><label><input type="checkbox">Create histograms with matplotlib</label></li>
<li><label><input type="checkbox">Filter and manipulate data using pandas operations</label></li>
<li><label><input type="checkbox">Create multi-panel figures for comparative analysis</label></li>
<li><label><input type="checkbox">Interpret energy deposition distributions physically</label></li>
</ul>
</section>
<section id="conclusion" class="level2">
<h2 class="anchored" data-anchor-id="conclusion">Conclusion</h2>
<p>You have successfully modified a Geant4 application to output data, run simulations with different particle types, and performed sophisticated data analysis using Python’s pandas and matplotlib libraries. You’ve learned:</p>
<ul>
<li>How Geant4’s action classes work together to record simulation data</li>
<li>The modern approach to data output using G4AnalysisManager</li>
<li>Essential pandas operations for data manipulation and analysis</li>
<li>Creating publication-quality visualizations with matplotlib</li>
<li>The scientific workflow from simulation to analysis to visualization</li>
</ul>
<p>These skills are fundamental for computational physics and will be used throughout the remainder of this course and in professional scientific computing work.</p>
</section>
<section id="references" class="level2">
<h2 class="anchored" data-anchor-id="references">References</h2>
<ul>
<li>Geant4 Collaboration. (n.d.). <em>Geant4 Analysis Tools Documentation</em></li>
<li>Hrivnacova, I. (2018). <em>Using Geant4 Analysis Manager</em></li>
<li>Hrivnacova, I. (2023). <em>Geant4 Analysis Tools Tutorial</em></li>
<li>McKinney, W. (2010). <em>Data Structures for Statistical Computing in Python</em></li>
<li>Hunter, J. D. (2007). <em>Matplotlib: A 2D Graphics Environment</em></li>
<li>The Python Graph Gallery. (n.d.). <em>Histogram Examples</em></li>
</ul>
</section>
<section id="additional-resources" class="level2">
<h2 class="anchored" data-anchor-id="additional-resources">Additional Resources</h2>
<ul>
<li><a href="https://pandas.pydata.org/docs/">Pandas Documentation</a></li>
<li><a href="https://matplotlib.org/stable/gallery/index.html">Matplotlib Gallery</a></li>
<li><a href="https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/analysis.html">Geant4 Analysis Manager Examples</a></li>
<li><a href="https://jakevdp.github.io/PythonDataScienceHandbook/">Python Data Science Handbook</a></li>
</ul>
</section>
<section id="optional-extensions" class="level2">
<h2 class="anchored" data-anchor-id="optional-extensions">Optional Extensions</h2>
<p>If you finish early or want to explore further:</p>
<ol type="1">
<li><p><strong>Add More Columns</strong>: Modify the code to save additional data:</p>
<ul>
<li>Track length</li>
<li>Number of steps</li>
<li>Particle position at end of event</li>
</ul></li>
<li><p><strong>Different Geometries</strong>: Modify the detector geometry in <code>DetectorConstruction.cc</code> and observe how it affects energy deposition</p></li>
<li><p><strong>Energy Scan</strong>: Create a script that runs simulations at multiple energies (e.g., 1, 2, 4, 6, 8, 10 MeV) and plots how mean energy deposition varies with incident energy</p></li>
<li><p><strong>Interactive Plotting</strong>: Modify the Python scripts to use <code>plt.show()</code> for interactive exploration of the data</p></li>
<li><p><strong>Statistical Tests</strong>: Use scipy to perform statistical tests comparing different datasets (t-tests, KS tests, etc.)</p></li>
</ol>
</section>
</main> <!-- /main -->
<script id="quarto-html-after-body" type="application/javascript">
window.document.addEventListener("DOMContentLoaded", function (event) {
const toggleBodyColorMode = (bsSheetEl) => {
const mode = bsSheetEl.getAttribute("data-mode");
const bodyEl = window.document.querySelector("body");
if (mode === "dark") {
bodyEl.classList.add("quarto-dark");
bodyEl.classList.remove("quarto-light");
} else {
bodyEl.classList.add("quarto-light");
bodyEl.classList.remove("quarto-dark");
}
}
const toggleBodyColorPrimary = () => {
const bsSheetEl = window.document.querySelector("link#quarto-bootstrap");
if (bsSheetEl) {
toggleBodyColorMode(bsSheetEl);
}
}
toggleBodyColorPrimary();
const icon = "";
const anchorJS = new window.AnchorJS();
anchorJS.options = {
placement: 'right',
icon: icon
};
anchorJS.add('.anchored');
const isCodeAnnotation = (el) => {
for (const clz of el.classList) {
if (clz.startsWith('code-annotation-')) {
return true;
}
}
return false;
}
const onCopySuccess = function(e) {
// button target
const button = e.trigger;
// don't keep focus
button.blur();
// flash "checked"
button.classList.add('code-copy-button-checked');
var currentTitle = button.getAttribute("title");
button.setAttribute("title", "Copied!");
let tooltip;
if (window.bootstrap) {
button.setAttribute("data-bs-toggle", "tooltip");
button.setAttribute("data-bs-placement", "left");
button.setAttribute("data-bs-title", "Copied!");
tooltip = new bootstrap.Tooltip(button,
{ trigger: "manual",
customClass: "code-copy-button-tooltip",
offset: [0, -8]});
tooltip.show();
}
setTimeout(function() {
if (tooltip) {
tooltip.hide();
button.removeAttribute("data-bs-title");
button.removeAttribute("data-bs-toggle");
button.removeAttribute("data-bs-placement");
}
button.setAttribute("title", currentTitle);
button.classList.remove('code-copy-button-checked');
}, 1000);
// clear code selection
e.clearSelection();
}
const getTextToCopy = function(trigger) {
const codeEl = trigger.previousElementSibling.cloneNode(true);
for (const childEl of codeEl.children) {
if (isCodeAnnotation(childEl)) {
childEl.remove();
}
}
return codeEl.innerText;
}
const clipboard = new window.ClipboardJS('.code-copy-button:not([data-in-quarto-modal])', {
text: getTextToCopy
});
clipboard.on('success', onCopySuccess);
if (window.document.getElementById('quarto-embedded-source-code-modal')) {
const clipboardModal = new window.ClipboardJS('.code-copy-button[data-in-quarto-modal]', {
text: getTextToCopy,
container: window.document.getElementById('quarto-embedded-source-code-modal')
});
clipboardModal.on('success', onCopySuccess);
}
var localhostRegex = new RegExp(/^(?:http|https):\/\/localhost\:?[0-9]*\//);
var mailtoRegex = new RegExp(/^mailto:/);
var filterRegex = new RegExp("https:\/\/jwinterm\.github\.io\/geant4-radiation-labs\/");
var isInternal = (href) => {
return filterRegex.test(href) || localhostRegex.test(href) || mailtoRegex.test(href);
}
// Inspect non-navigation links and adorn them if external
var links = window.document.querySelectorAll('a[href]:not(.nav-link):not(.navbar-brand):not(.toc-action):not(.sidebar-link):not(.sidebar-item-toggle):not(.pagination-link):not(.no-external):not([aria-hidden]):not(.dropdown-item):not(.quarto-navigation-tool):not(.about-link)');
for (var i=0; i<links.length; i++) {
const link = links[i];
if (!isInternal(link.href)) {
// undo the damage that might have been done by quarto-nav.js in the case of
// links that we want to consider external
if (link.dataset.originalHref !== undefined) {
link.href = link.dataset.originalHref;
}
// target, if specified
link.setAttribute("target", "_blank");
if (link.getAttribute("rel") === null) {
link.setAttribute("rel", "noopener");
}
}
}
function tippyHover(el, contentFn, onTriggerFn, onUntriggerFn) {
const config = {
allowHTML: true,
maxWidth: 500,
delay: 100,
arrow: false,
appendTo: function(el) {
return el.parentElement;
},
interactive: true,
interactiveBorder: 10,
theme: 'quarto',
placement: 'bottom-start',
};
if (contentFn) {
config.content = contentFn;
}
if (onTriggerFn) {
config.onTrigger = onTriggerFn;
}
if (onUntriggerFn) {
config.onUntrigger = onUntriggerFn;
}
window.tippy(el, config);
}
const noterefs = window.document.querySelectorAll('a[role="doc-noteref"]');
for (var i=0; i<noterefs.length; i++) {
const ref = noterefs[i];
tippyHover(ref, function() {
// use id or data attribute instead here
let href = ref.getAttribute('data-footnote-href') || ref.getAttribute('href');
try { href = new URL(href).hash; } catch {}
const id = href.replace(/^#\/?/, "");
const note = window.document.getElementById(id);
if (note) {
return note.innerHTML;
} else {
return "";
}
});
}
const xrefs = window.document.querySelectorAll('a.quarto-xref');
const processXRef = (id, note) => {
// Strip column container classes
const stripColumnClz = (el) => {
el.classList.remove("page-full", "page-columns");
if (el.children) {
for (const child of el.children) {
stripColumnClz(child);
}
}
}
stripColumnClz(note)
if (id === null || id.startsWith('sec-')) {
// Special case sections, only their first couple elements
const container = document.createElement("div");
if (note.children && note.children.length > 2) {
container.appendChild(note.children[0].cloneNode(true));
for (let i = 1; i < note.children.length; i++) {
const child = note.children[i];
if (child.tagName === "P" && child.innerText === "") {
continue;
} else {
container.appendChild(child.cloneNode(true));
break;
}
}
if (window.Quarto?.typesetMath) {
window.Quarto.typesetMath(container);
}
return container.innerHTML
} else {
if (window.Quarto?.typesetMath) {
window.Quarto.typesetMath(note);
}
return note.innerHTML;
}
} else {
// Remove any anchor links if they are present
const anchorLink = note.querySelector('a.anchorjs-link');
if (anchorLink) {
anchorLink.remove();
}
if (window.Quarto?.typesetMath) {
window.Quarto.typesetMath(note);
}
if (note.classList.contains("callout")) {
return note.outerHTML;
} else {
return note.innerHTML;
}
}
}
for (var i=0; i<xrefs.length; i++) {
const xref = xrefs[i];
tippyHover(xref, undefined, function(instance) {
instance.disable();
let url = xref.getAttribute('href');
let hash = undefined;
if (url.startsWith('#')) {
hash = url;
} else {
try { hash = new URL(url).hash; } catch {}
}
if (hash) {
const id = hash.replace(/^#\/?/, "");
const note = window.document.getElementById(id);
if (note !== null) {
try {
const html = processXRef(id, note.cloneNode(true));
instance.setContent(html);
} finally {
instance.enable();
instance.show();
}
} else {
// See if we can fetch this
fetch(url.split('#')[0])
.then(res => res.text())
.then(html => {
const parser = new DOMParser();
const htmlDoc = parser.parseFromString(html, "text/html");
const note = htmlDoc.getElementById(id);
if (note !== null) {
const html = processXRef(id, note);
instance.setContent(html);
}
}).finally(() => {
instance.enable();
instance.show();
});
}
} else {
// See if we can fetch a full url (with no hash to target)
// This is a special case and we should probably do some content thinning / targeting
fetch(url)
.then(res => res.text())
.then(html => {
const parser = new DOMParser();
const htmlDoc = parser.parseFromString(html, "text/html");
const note = htmlDoc.querySelector('main.content');
if (note !== null) {
// This should only happen for chapter cross references
// (since there is no id in the URL)
// remove the first header
if (note.children.length > 0 && note.children[0].tagName === "HEADER") {
note.children[0].remove();
}
const html = processXRef(null, note);
instance.setContent(html);
}
}).finally(() => {
instance.enable();
instance.show();
});
}
}, function(instance) {
});
}
let selectedAnnoteEl;
const selectorForAnnotation = ( cell, annotation) => {
let cellAttr = 'data-code-cell="' + cell + '"';
let lineAttr = 'data-code-annotation="' + annotation + '"';
const selector = 'span[' + cellAttr + '][' + lineAttr + ']';
return selector;
}
const selectCodeLines = (annoteEl) => {
const doc = window.document;
const targetCell = annoteEl.getAttribute("data-target-cell");
const targetAnnotation = annoteEl.getAttribute("data-target-annotation");
const annoteSpan = window.document.querySelector(selectorForAnnotation(targetCell, targetAnnotation));
const lines = annoteSpan.getAttribute("data-code-lines").split(",");
const lineIds = lines.map((line) => {
return targetCell + "-" + line;
})
let top = null;
let height = null;
let parent = null;
if (lineIds.length > 0) {
//compute the position of the single el (top and bottom and make a div)
const el = window.document.getElementById(lineIds[0]);
top = el.offsetTop;
height = el.offsetHeight;
parent = el.parentElement.parentElement;
if (lineIds.length > 1) {
const lastEl = window.document.getElementById(lineIds[lineIds.length - 1]);
const bottom = lastEl.offsetTop + lastEl.offsetHeight;
height = bottom - top;
}
if (top !== null && height !== null && parent !== null) {
// cook up a div (if necessary) and position it
let div = window.document.getElementById("code-annotation-line-highlight");
if (div === null) {
div = window.document.createElement("div");
div.setAttribute("id", "code-annotation-line-highlight");
div.style.position = 'absolute';
parent.appendChild(div);
}
div.style.top = top - 2 + "px";
div.style.height = height + 4 + "px";
div.style.left = 0;
let gutterDiv = window.document.getElementById("code-annotation-line-highlight-gutter");
if (gutterDiv === null) {
gutterDiv = window.document.createElement("div");
gutterDiv.setAttribute("id", "code-annotation-line-highlight-gutter");
gutterDiv.style.position = 'absolute';
const codeCell = window.document.getElementById(targetCell);
const gutter = codeCell.querySelector('.code-annotation-gutter');
gutter.appendChild(gutterDiv);
}
gutterDiv.style.top = top - 2 + "px";
gutterDiv.style.height = height + 4 + "px";
}
selectedAnnoteEl = annoteEl;
}
};
const unselectCodeLines = () => {
const elementsIds = ["code-annotation-line-highlight", "code-annotation-line-highlight-gutter"];
elementsIds.forEach((elId) => {
const div = window.document.getElementById(elId);
if (div) {
div.remove();
}
});
selectedAnnoteEl = undefined;
};
// Handle positioning of the toggle
window.addEventListener(
"resize",
throttle(() => {
elRect = undefined;
if (selectedAnnoteEl) {
selectCodeLines(selectedAnnoteEl);
}
}, 10)
);
function throttle(fn, ms) {
let throttle = false;
let timer;
return (...args) => {
if(!throttle) { // first call gets through
fn.apply(this, args);
throttle = true;
} else { // all the others get throttled
if(timer) clearTimeout(timer); // cancel #2
timer = setTimeout(() => {
fn.apply(this, args);
timer = throttle = false;
}, ms);
}
};
}
// Attach click handler to the DT
const annoteDls = window.document.querySelectorAll('dt[data-target-cell]');
for (const annoteDlNode of annoteDls) {
annoteDlNode.addEventListener('click', (event) => {
const clickedEl = event.target;
if (clickedEl !== selectedAnnoteEl) {
unselectCodeLines();
const activeEl = window.document.querySelector('dt[data-target-cell].code-annotation-active');
if (activeEl) {
activeEl.classList.remove('code-annotation-active');
}
selectCodeLines(clickedEl);
clickedEl.classList.add('code-annotation-active');
} else {
// Unselect the line
unselectCodeLines();
clickedEl.classList.remove('code-annotation-active');
}
});
}
const findCites = (el) => {
const parentEl = el.parentElement;
if (parentEl) {
const cites = parentEl.dataset.cites;
if (cites) {
return {
el,
cites: cites.split(' ')
};
} else {
return findCites(el.parentElement)
}
} else {
return undefined;
}
};
var bibliorefs = window.document.querySelectorAll('a[role="doc-biblioref"]');
for (var i=0; i<bibliorefs.length; i++) {
const ref = bibliorefs[i];
const citeInfo = findCites(ref);
if (citeInfo) {
tippyHover(citeInfo.el, function() {
var popup = window.document.createElement('div');
citeInfo.cites.forEach(function(cite) {
var citeDiv = window.document.createElement('div');
citeDiv.classList.add('hanging-indent');
citeDiv.classList.add('csl-entry');
var biblioDiv = window.document.getElementById('ref-' + cite);
if (biblioDiv) {
citeDiv.innerHTML = biblioDiv.innerHTML;
}
popup.appendChild(citeDiv);
});
return popup.innerHTML;
});
}
}
});
</script>
<nav class="page-navigation">
<div class="nav-page nav-page-previous">
<a href="../lab2_exampleB1/lab2.html" class="pagination-link" aria-label="Lab 2: First Simulation - Building and Running Geant4 Example B1">
<i class="bi bi-arrow-left-short"></i> <span class="nav-page-text">Lab 2: First Simulation - Building and Running Geant4 Example B1</span>
</a>
</div>
<div class="nav-page nav-page-next">
<a href="../lab4_exampleB3/lab4.html" class="pagination-link" aria-label="Lab 4: Medical Physics Application - A PET Scanner">
<span class="nav-page-text">Lab 4: Medical Physics Application - A PET Scanner</span> <i class="bi bi-arrow-right-short"></i>
</a>
</div>
</nav>
</div> <!-- /content -->
</body></html>