NIST-EPA Legionella Study

Repository for NIST-EPA Legionella study Python tools. Includes scripts for aerosol particle size distribution analysis, emission-rate calculations, data validation, visualization, and statistical modeling of shower-generated aerosols under varying temperature, humidity, and exhaust-fan conditions.

Project Overview

This project analyzes aerosol characteristics from shower-generated aerosols with a focus on potential Legionella transmission risks from shower fixtures. The research investigates:

• Aerosol particle size distribution analysis
• Emission-rate calculations using numerical approaches
• CO2 decay analysis for air change rate determination
• Temperature and relative humidity monitoring
• data validation, visualization, and statistical modeling

A second, related line of analysis (the moduair_* scripts and src/moduair_loader.py) compares the study's two QuantAQ MODULAIR-PM sensors against the broader NIST IAQ&V group's fleet of co-located MODULAIR-PM units, to characterize sensor agreement and build candidate inter-sensor correction factors. See MODULAIR-PM Fleet Comparison below.

Disclaimer

Certain commercial equipment, instruments, software, or materials are identified in this repository in order to specify the experimental and analytical procedures adequately. Such identification is not intended to imply recommendation or endorsement of any product or service by NIST, nor is it intended to imply that the materials or equipment identified are necessarily the best available for the purpose.

Measurement Uncertainty

The analysis procedures in this repository involve multiple stages of data processing and mathematical modeling, each contributing to the overall uncertainty in the final results. Key sources of uncertainty include:

CO₂-derived Air Change Rate (λ) Uncertainty:

• Statistical uncertainty from linear regression of log-transformed decay data (reported as standard error of the slope)
• Systematic uncertainty from assuming a well-mixed zone and constant infiltration fractions (α, β)
• Uncertainty in source concentration characterization evaluated through three methods: weighted average (α·C_outside + β·C_entry), outside-only, and entry-only
• Data quality filters including minimum concentration difference requirements (50 ppm) and sufficient data point thresholds
• Temporal alignment uncertainty between CO₂ injection timing and decay window selection

Particle Analysis Uncertainty:

• Penetration factor (p) uncertainty from temporal variability in indoor/outdoor concentration ratios during before/after windows
• Other process rate (β) uncertainty from numerical differentiation of concentration data, mitigated by 5th–95th percentile trimming and R²-based selection criteria
• Emission rate (E) uncertainty from propagation of errors in the mass balance equation and concentration time derivatives
• Forward Euler simulation uncertainty from accumulation of integration errors over the 2-hour prediction window
• Size-bin dependent uncertainties due to varying signal-to-noise ratios across the 0.35–10.0 µm range

General Measurement Uncertainties:

• Sensor-specific uncertainties: Aranet4 CO₂ sensors (±50 ppm + 3% of reading), QuantAQ MODULAIR-PM particle counters, HOBO UX100 temperature/relative humidity sensors
• Temporal synchronization uncertainty between different sensor systems sampling at 1-minute intervals
• Spatial representativeness limitations of point measurements in heterogeneous indoor environments
• Model structural uncertainty from assuming first-order decay processes and uniform mixing in the mass balance equations

Uncertainty quantification is implemented throughout the analysis pipeline, with standard deviations, standard errors, and confidence intervals reported where applicable. The three-source-method approach for λ calculation provides explicit assessment of structural uncertainty in ventilation rate determination. Where measurements fall below detection thresholds or fail quality criteria, results are appropriately flagged or excluded to prevent overinterpretation.

Repository Structure

NIST_EPA_Legionella/
├── .vscode/                          # VS Code configuration
│   └── settings.json
├── .env                              # Environment variables (API keys for QuantAQ, gitignored)
├── .gitignore
├── CODEMETA.yaml                     # NIST Software Portal metadata
├── CODEOWNERS                        # Repository ownership for PR reviews
├── LICENSE.md                        # NIST software licensing statement
├── README.md
├── fair-software.md                  # FAIR software principles compliance checklist
├── data_config.json                  # Active configuration (gitignored)
├── data_config.template.json         # Configuration template
├── environment.yaml                  # Conda environment specification
├── run_analysis_workflow.py          # Convenience script to run the full analysis pipeline
│
├── src/                              # Core analysis modules (imported by scripts)
│   ├── __init__.py                   # Package initialization (re-exports data_paths)
│   ├── data_paths.py                 # Portable data access via data_config.json
│   ├── env_data_loader.py            # Unified environmental sensor data loader
│   ├── event_manager.py              # Event filtering, naming (MMDD_Temp_ToD_RNN), logging
│   ├── event_matching.py             # Match CO2 injections to shower events by timing
│   ├── particle_data_loader.py       # QuantAQ particle data loading and preparation
│   ├── particle_room_correction.py   # Room-concentration correction for the "inside" series (C_room fleet average / pre-fleet ratio correction)
│   ├── particle_calculations.py      # Penetration, deposition & emission calculations
│   ├── plot_co2.py                   # CO2 decay visualization functions
│   ├── plot_environmental.py         # RH, temperature, wind visualization
│   ├── plot_particle.py              # Particle decay visualization (re-exports boxplot/comparison functions)
│   ├── plot_particle_boxplots.py     # Temperature-axis and shower-head boxplot functions
│   ├── plot_comparison.py            # Categorical comparison boxplot families (door, fan, spray pattern, etc.)
│   ├── plot_style.py                 # Consistent matplotlib styling constants
│   ├── plot_utils.py                 # Central plotting re-exports (imports all plot_* modules)
│   ├── quantaq_utils.py              # QuantAQ API client utilities
│   ├── sig_figs.py                   # Significant figures rounding and formatting
│   ├── moduair_loader.py             # IAQ&V MODULAIR-PM fleet chunk loader (moduair_pm_chunks share)
│   └── deprecated/                   # Deprecated/archived code
│       └── co2_decay_analysis.py     # Previous version of CO2 analysis
│
├── scripts/                          # Executable analysis scripts (run directly)
│   ├── __init__.py
│   │
│   │ # Data Download & Processing
│   ├── download_quantaq_data.py      # Download QuantAQ sensor data from API (manual/one-off; campaign ended 2026-07-16)
│   ├── process_quantaq_data.py       # Combine weekly chunks and process raw/final QuantAQ data
│   │
│   │ # Log Processing
│   ├── process_co2_log.py            # Consolidate daily CO2 injection logs → state-change log
│   ├── process_shower_log.py         # Consolidate daily shower logs → state-change log
│   ├── fix_log_files.py              # Repair corrupted/empty log files
│   ├── analyze_log_files.py          # Diagnose log file issues (unmatched events, etc.)
│   │
│   │ # Event Management
│   ├── event_registry.py             # Unified event registry with synthetic event creation
│   │
│   │ # Analysis
│   ├── co2_decay_analysis.py         # CO2 decay & air change rate (λ) analysis
│   ├── particle_decay_analysis.py    # Particle penetration, deposition & 4-variant emission analysis (E1-E4; per-event figures)
│   ├── particle_emission_variant_figures.py  # Run after particle_decay_analysis.py: per-bin E1/E2/E3/E4 comparison figures + all summary/boxplot/comparison figures
│   ├── rh_temp_other_analysis.py     # RH, temperature & wind condition analysis
│   ├── export_event_timeseries.py    # Export per-event predicted Ct for all size bins to CSV (--event N)
│   ├── export_config_timeseries.py   # Export 1-min avg env/particle time series per config group to Excel
│   │
│   │ # Optional Figures
│   ├── co2_lambda_source_figure.py   # Per-event λ_outside vs. λ_entry comparison figure with paired t-test
│   ├── hobo_onset_decay_figures.py   # HOBO temp/RH onset & decay figures, plus per-event pre/post summary
│   │
│   │ # MODULAIR-PM Fleet Comparison (IAQ&V) — see README section below
│   ├── moduair_correction_factor.py  # Inter-sensor concentration ratios + 195/813 Deming correction fit
│   ├── moduair_event_peak_times.py   # Per-event, per-sensor delta peak time (shower ON to peak)
│   ├── moduair_bin2_timeseries.py    # Per-bin time series across co-located bedroom fleet sensors
│   ├── moduair_cave_ratio.py         # Inside sensor (195) vs. position-weighted fleet average (C_room)
│   └── export_config_timeseries_fleet.py  # Per-sensor variant of export_config_timeseries.py across the fleet
│
└── docs/                             # Documentation & instrument references
    ├── Data Analysis.docx            # Data analysis planning notes
    └── pdf_links.md                  # Public download links for instrument manuals and data sheets (not included due to copyright)

Installation

Prerequisites

• Python 3.14+
• Conda (recommended for environment management)

Setup

1. Clone the repository:

   git clone https://github.com/usnistgov/NIST_EPA_Legionella.git
   cd NIST_EPA_Legionella

2. Create the conda environment:

   conda env create -f environment.yaml
   conda activate epa_mh

3. Configure data paths:

   cp data_config.template.json data_config.json

   Edit data_config.json to set the correct data_root path for your machine.

4. Set up API credentials (for QuantAQ data): Create a .env file with your QuantAQ API key:

   QUANTAQ_API_KEY=your_api_key_here
   

Usage

Data Access

The src/data_paths.py module provides portable data access functions:

from src import get_instrument_path, get_instrument_file, get_instrument_files_for_date_range

# Get instrument data directory
path = get_instrument_path("Aranet4")

# Get data file for a specific date
file = get_instrument_file("Aranet4", "2026-01-15")

# Get files for a date range
files = get_instrument_files_for_date_range("QuantAQ", "2026-01-05", "2026-01-15")

QuantAQ Data Pipeline

The experimental campaign ended 2026-07-16; the local QuantAQ cache is complete and frozen, so step 1 below is a manual/one-off tool, not something to re-run routinely.

1. Download data from the QuantAQ API (manual, one-off — only needed to rebuild a lost/missing local cache):

   python scripts/download_quantaq_data.py

2. Process data (parse nested structures, merge raw/final data):

   python scripts/process_quantaq_data.py

Visualization Tools

Plotting modules in src/ generate publication-quality figures used by the analysis scripts:

• src/plot_co2.py - CO2 decay curves and air change rate analysis
• src/plot_particle.py - Particle concentration decay and emission rates; re-exports boxplot/comparison functions
• src/plot_particle_boxplots.py - Temperature-axis and shower-head E_total/β/E/p boxplots
• src/plot_comparison.py - Categorical comparison boxplot families (door, fan, spray pattern, mannequin, head type)
• src/plot_environmental.py - RH, temperature, and wind time series
• src/plot_style.py - Consistent matplotlib styling constants and helpers
• src/plot_utils.py - Re-exports all plot_* functions for backward compatibility

All plot modules are imported by analysis scripts through src/plot_utils.py.

Configured Instruments

Instrument	Model	Purpose	Data Location
AIO2	Met One AIO 2	Outdoor weather (10m tower)	MH DAQ/weather_station
Aranet4	Aranet4 PRO	CO2 monitoring (Bedroom, Bathroom, Entry, Outside)	CO2
HOBO_UX100	Onset HOBO UX100	RH/Temperature (6 locations: Bathroom1×2, Bathroom2, Bath/Bed, Bedroom1–3)	HOBOs
QuantAQ	MODULAIR-PM	PM monitoring (indoor: MOD-PM-00195, outdoor: MOD-PM-00785)	QuantAQ
Setra_264	Setra 264	Differential pressure	MH DAQ/indoor_daq
Vaisala_HMP155	Vaisala HMP155	RH/Temperature	MH DAQ/indoor_daq
Vaisala_HMP45A	Vaisala HMP45A	RH/Temperature	MH DAQ/indoor_daq

DAQ System

• Chassis: NI cDAQ-9178 (8-slot USB)
• Modules: 5x NI 9201 (8-channel analog input, +/-10V)

Analysis Workflow

The complete data analysis pipeline follows this sequence. To run all steps automatically:

python run_analysis_workflow.py

Each step's output is logged to output/logs/<script_name>.log. Or run steps individually:

Data collection is complete. The experimental campaign ended 2026-07-16 (PROJECT_END_DATE in download_quantaq_data.py), so download_quantaq_data.py is no longer part of the routine pipeline in run_analysis_workflow.py. Run it manually, one-off, only if the local QuantAQ chunks/ cache needs to be rebuilt:

python scripts/download_quantaq_data.py

1. Data Processing: Process and merge raw QuantAQ sensor data

   python scripts/process_quantaq_data.py

2. Log Processing: Consolidate daily 1-second log files into state-change logs

   python scripts/process_co2_log.py
   python scripts/process_shower_log.py

3. Event Management: Run the event manager to match, name, and register events

   python scripts/event_registry.py [--force] [--no-co2] [--output-dir DIR]

   Flag	Description
   --force	Regenerate registry even if it already exists
   --no-co2	Skip CO2 analysis (registry will have no lambda values or PM exclusions)
   --output-dir DIR	Output directory (default: data_root/output)

   This must be run before the analysis scripts so that events have consistent names (see Event Naming Convention below). Produces event_log.csv.

4. CO2 Analysis: Determine air change rates (λ) from CO2 decay

   python scripts/co2_decay_analysis.py [--alpha FLOAT] [--beta FLOAT] [--output-dir DIR] [--no-plot] [--no-sig-figs] [--entry-stop]

   Flag	Description
   --alpha FLOAT	Fraction of infiltration from outside (default: 0.5)
   --beta FLOAT	Fraction of infiltration from entry zone (default: 0.5)
   --output-dir DIR	Output directory (default: data_root/output)
   --no-plot	Disable plot generation
   --no-sig-figs	Disable significant figure rounding (default: 3 sig figs for files, 2 for figures)
   --entry-stop	Truncate the decay window at the first point where C_bedroom ≤ C_entry + 50 ppm (falls back to C_outside if entry data is missing) — stops the fit once the remaining decay signal above background is too small to be reliable

   run_analysis_workflow.py always runs this step with --entry-stop. Fits with R² below MIN_R_SQUARED (0.75) are rejected (λ set to NaN); this also drives the "Lambda R² less than 0.75" PM-exclusion check in event_registry.py.

5. Environmental Analysis: Characterize RH, temperature, and wind conditions

   python scripts/rh_temp_other_analysis.py [--output-dir DIR] [--no-plots] [--max-plot-events N] [--events LIST] [--no-sig-figs]

   Flag	Description
   --output-dir DIR	Output directory (default: data_root/output)
   --no-plots	Skip plot generation
   --max-plot-events N	Limit time series plots to first N events (default: all)
   --events LIST	Comma-separated event numbers to plot, e.g. 1,3,5 (overrides --max-plot-events)
   --no-sig-figs	Disable significant figure rounding (default: 3 sig figs for files, 2 for figures)

6. Particle Analysis: Calculate penetration, deposition, and emission rates (requires step 4 results)

   python scripts/particle_decay_analysis.py [--output-dir DIR] [--no-plot] [--no-sig-figs]

   Flag	Description
   --output-dir DIR	Output directory (default: data_root/output)
   --no-plot	Disable per-event pm_decay figure generation
   --no-sig-figs	Disable significant figure rounding (default: 3 sig figs for files, 2 for figures)

   Four emission-rate variants (E1-E4) are computed per bin, per the report's Table 7 / Equation 10, rather than one value bounded by λ_outside/λ_entry alone. See Particle Decay & Emission Analysis below. Summary/boxplot/comparison figures moved to step 7 below; this step only generates the per-event pm_decay figures and the Excel results.

   New dependency: the indoor ("inside") concentration is now room-corrected by src/particle_room_correction.py before analysis (see the "Room-Concentration Correction" note under Particle Decay & Emission Analysis), which requires the MODULAIR-PM fleet chunk files to be reachable on the moduair_pm_chunks share configured in data_config.json (see MODULAIR-PM Fleet Comparison below). This reuses scripts/moduair_cave_ratio.py's functions in-memory — that script does not need to be run first. If the fleet chunks are unreachable, the step prints a [WARN] and falls back to uncorrected raw MOD-PM-00195 for the whole record.

7. Particle Emission Variant Figures: Per-bin E1/E2/E3/E4 comparison figures, plus all summary/boxplot/comparison figures (requires step 6 results)

   python scripts/particle_emission_variant_figures.py [--output-dir DIR]

   Reads particle_analysis_summary.xlsx; does not recompute p, β, or E. Produces the per-bin beta_other_{pair}_bin{N}.html / emission_{pair}_bin{N}.html comparison figures for three variant pairs (entry_outside, outside_adjusted, outside_bathroom; 12 bins × 2 metrics × 3 pairs = 72 figures) under output/plots/particle/, plus the full summary/boxplot/comparison figure suite (bar charts, fixed- and metric-axis boxplots, and the 5 condition-comparison families) under output/plots/ — moved here from step 6.

8. Export Event Time Series (optional utility): Export predicted Ct for all size bins for a single event to CSV:

   python scripts/export_event_timeseries.py --event N [--output PATH] [--output-dir DIR] [--no-sig-figs]

   Flag	Description
   --event INT	Event number to export (required)
   --output PATH	Output CSV file path (default: output/event_N_timeseries.csv)
   --output-dir DIR	Analysis output directory (default: data_root/output)
   --no-sig-figs	Disable significant figure rounding

9. Export Config Time Series (optional utility): Export 1-minute averaged environmental and particle time series per unique test configuration to Excel:

   python scripts/export_config_timeseries.py [--no-sig-figs]

   Flag	Description
   --no-sig-figs	Disable significant figure rounding

   Produces output/event_config_timeseries.xlsx with one sheet per configuration group. Each sheet contains the mean, standard deviation, maximum, and minimum (across replicate events) of bathroom RH/Temp, bath/bed hallway RH/Temp, bedroom RH/Temp, and all 12 OPC-N3 particle bins (0.35–10.0 µm) at 1-minute resolution from shower-on through the 2-hour deposition window. Flow rates in the 4.1–5.6 LPM range (standard and tagged variants) are grouped together; 1.4 LPM and 2.2 LPM events remain separate. A hyperlinked index sheet lists every configuration group with event counts and event numbers.

10. CO2 Lambda Source Comparison Figure (optional utility): Build a per-event λ_outside vs. λ_entry comparison figure with a paired t-test (requires step 4 results):

    python scripts/co2_lambda_source_figure.py

    Reads co2_lambda_summary.csv; does not recompute λ. Produces output/plots/co2/co2_lambda_entry_outside.html.

11. HOBO Onset/Decay Figures (optional utility): Build interactive HOBO temp/RH figures restricted to the onset (shower-on to +60 min) and decay (shower-off +60 to +120 min) windows around each event, plus a per-event pre/post summary:

    python scripts/hobo_onset_decay_figures.py

    Produces output/plots/hobo/hobo_{temp,rh}_onset_decay.html and output/plots/hobo/hobo_{temp,rh}_pre_post.html.

Each analysis module produces CSV/Excel summaries and optional visualizations in the output/ directory.

Event Naming Convention

Events are named using a structured format generated by scripts/event_registry.py. This naming convention must be applied (step 4 in the workflow above) before running the analysis scripts.

Format: MMDD_W##[_ShowerHead[_SprayPattern]][_Mannequin]_DoorPos[_Fan]_RNN

Component	Description	Values
MMDD	Month and day of the event	e.g., 0115 for January 15
W##	Water temperature code (°C)	W11, W14, W22, W23, W24, W25, W30, W37, W38, W40, W43, W48, W49, W52, W53
ShowerHead	Shower head type; omitted for Standard head	Pepco, FilterWand, Used
SprayPattern	Spray pattern; varies by head type	Wide, Narrow, Mid, rainfall, 12Nozzle, SingleWide1, SingleWide2
Mannequin	Appended only when a mannequin was present	Mannequin
DoorPos	Bathroom-to-bedroom (bath) door position	Open, Closed, Ajar
Fan	Appended only when the bath fan ran during the test	Fan
RNN	Replicate number	R01, R02, R03, etc.

Examples:

• 0115_W48_Open_R01 - January 15, 48 °C, standard head, bath door open, replicate 1
• 0122_W11_Open_R03 - January 22, 11 °C, standard head, replicate 3
• 0224_W52_Pepco_Wide_Open_R01 - February 24, 52 °C, Pepco head, wide spray
• 0309_W40_Pepco_Narrow_Open_R01 - March 9, 40 °C, Pepco head, narrow spray
• 0311_W40_Pepco_Narrow_Mannequin_Open_R01 - March 11, mannequin present
• 0408_W38_Pepco_Wide_Closed_R01 - April 8, 38 °C, Pepco head, wide spray, bath door closed

Note on config_key: The config_key stored in event_log.csv and spreadsheet outputs is a separate, more detailed identifier used for grouping in boxplots. It includes both door positions, fan status, and measured flow rate: W##[_ShowerHead[_SprayPattern]][_Mannequin]_BathDoorXxx_BdrmDoorXxx_FanXxx[_FlowRateX.XLPM] (e.g., W48_BathDoorOpen_BdrmDoorClosed_FanOff_FlowRate4.6LPM).

Test Parameter Timeline (Standard shower head):

• W48 (48 °C): From experiment start (2026-01-14); bath door open, bedroom door closed
• W11 (11 °C): Starting 2026-01-22 14:00
• W25 (25 °C): Starting 2026-02-02 17:00
• W30 (30 °C): Starting 2026-02-05 10:00
• W37 (37 °C): Starting 2026-02-09 10:00
• W23 (23 °C): Starting 2026-02-11 08:00 (data excluded; kept for documentation)
• W22 (22 °C): Starting 2026-02-13 11:00
• W43 (43 °C): Starting 2026-02-16 11:00
• W14 (14 °C): Starting 2026-02-18 10:23
• W53 (53 °C): Starting 2026-02-20 08:00

Test Parameter Timeline (Pepco shower head, first phase, installed 2026-02-24):

• W52 (52 °C), Wide spray: Starting 2026-02-24 08:00
• W49 (49 °C), Wide spray: Starting 2026-02-26 10:23
• W40 (40 °C), Narrow spray: Starting 2026-03-02 09:23
• W40 (40 °C), Wide spray: Starting 2026-03-04 10:23
• W40 (40 °C), Mid spray: Starting 2026-03-06 08:47
• W40 (40 °C), Narrow spray: Starting 2026-03-09 08:47
• W40 (40 °C), Narrow spray + Mannequin: Starting 2026-03-11 09:47
• W38 (38 °C), Narrow spray: Starting 2026-03-12 10:47
• W38 (38 °C), Narrow spray + Mannequin: Starting 2026-03-13 08:45
• W38 (38 °C), Narrow spray (no mannequin): Starting 2026-03-17 08:25
• W38 (38 °C), Narrow spray + Mannequin: Starting 2026-03-18 09:45
• W38 (38 °C), Wide spray: Starting 2026-03-19 08:25
• W38 (38 °C), Wide spray + Mannequin: Starting 2026-03-22 09:25
• W38 (38 °C), Wide spray (no mannequin): Starting 2026-03-27 08:25

Test Parameter Timeline (FilterWand shower head, installed 2026-04-10):

• W38 (38 °C), FilterWand (no spray pattern variable): Starting 2026-04-10 08:55

Test Parameter Timeline (Used shower heads, first phase, starting 2026-04-13):

• W38 (38 °C), Used-rainfall: Starting 2026-04-13 11:35
• W38 (38 °C), Used-12Nozzle: Starting 2026-04-15 08:25
• W38 (38 °C), Used-SingleWide: Starting 2026-04-17 08:35
• W38 (38 °C), Used-SingleNarrow: Starting 2026-04-20 08:35

Test Parameter Timeline (Pepco shower head, second phase, reinstalled 2026-04-23):

• W38 (38 °C), Narrow spray + Mannequin: Starting 2026-04-23 14:15
• W38 (38 °C), Narrow spray (no mannequin): Starting 2026-04-29 08:35

Test Parameter Timeline (FilterWand shower head, second phase, reinstalled 2026-05-22):

• W38 (38 °C), FilterWand (no spray pattern variable): Starting 2026-05-22 09:30

Test Parameter Timeline (Used shower heads, second phase, reinstalled 2026-05-26):

• W38 (38 °C), Used-rainfall: Starting 2026-05-26 10:30
• W38 (38 °C), Used-12Nozzle: Starting 2026-06-01 08:00
• W38 (38 °C), Used-SingleNarrow: Starting 2026-06-03 08:30
• W38 (38 °C), Used-SingleWide: Starting 2026-06-08 13:15
• W38 (38 °C), Used-rainfall: Starting 2026-06-11 10:30
• W38 (38 °C), Used-SingleNarrow + Mannequin: Starting 2026-06-22 11:45 (mannequin present through 2026-07-02 09:45)
• W38 (38 °C), Used-SingleNarrow: Starting 2026-06-26 07:45

Water temperature changes under the Used head (no shower head/spray change):

• W49 (49 °C): Starting 2026-07-14 10:15
• W24 (24 °C): Starting 2026-07-15 10:15

Measured shower head flow rate is tracked continuously across the whole campaign (not just at head/pattern changes) and stored per event in config_key as an optional _FlowRateX.XLPM suffix; see FLOW_RATE_TRANSITIONS in src/event_manager.py for the full record.

Bath fan testing:

• Fan off: From experiment start through 2026-03-30
• Fan on (12 min per shower): Starting 2026-03-31 12:45
• Fan off: Starting ~2026-04-07 (confirmed off from 2026-04-10 08:55)

Bath door position changes:

• Bath door open (BathDoorOpen): From experiment start through 2026-04-07
• Bath door closed (BathDoorClosed): Starting 2026-04-08 09:15
• Bath door open (BathDoorOpen): Starting 2026-04-10 08:55
• Bath door ajar (BathDoorAjar): Starting 2026-05-04 08:30
• Bath door open (BathDoorOpen): Starting 2026-05-07 08:00

Bedroom door position changes:

• Bedroom door closed (BdrmDoorClosed): From experiment start (2026-01-14)
• Bedroom door ajar (BdrmDoorAjar): Starting 2026-05-01 17:10 (door blown ajar between May 1–4; intentionally kept ajar May 4–7)
• Bedroom door closed (BdrmDoorClosed): Starting 2026-05-07 08:00

Data Analysis Modules

CO2 Air Change Rate Analysis

Run the CO2 decay analysis to calculate air change rates (λ):

# Basic analysis (plots enabled by default)
python scripts/co2_decay_analysis.py

# Custom source fractions
python scripts/co2_decay_analysis.py --alpha 0.6 --beta 0.4  # default is 0.5/0.5

# Skip plots, write to a custom directory
python scripts/co2_decay_analysis.py --no-plot --output-dir /path/to/output

Methodology:

• Combines three Aranet4 sensor files: Bedroom, Entry, and Outside
• Applies 6-minute rolling average to reduce sensor noise
• Fixed 2-hour decay analysis window starting at :50 (10 minutes before the shower hour)
• Three source concentration methods for uncertainty: λ_average (mean of outside + entry), λ_outside, λ_entry
• Calculates λ via linear regression of the log-transformed decay: y = −ln[(C(t) − C_avg) / (C₀ − C_avg)], λ = slope of y vs. time
• Fits with R² below MIN_R_SQUARED (0.75) are rejected (λ set to NaN); this also drives the "Lambda R² less than 0.75" PM-analysis exclusion check in event_registry.py
• Optional --entry-stop truncates the decay window at the first point where C_bedroom ≤ C_entry + 50 ppm (falls back to C_outside if entry data is missing); run_analysis_workflow.py always runs with this flag

Output Files:

• output/co2_lambda_summary.csv - Per-event lambda results
• output/co2_lambda_overall_summary.csv - Overall statistics
• output/plots/event_figures/co2_decay/event_NN-testname_co2_decay.png - Individual event decay plots with shower ON/OFF markers
• output/plots/lambda_summary.png - Summary bar chart of λ values
• output/plots/air_change_rate_boxplot.png - Box-and-whisker of λ by water temperature

Particle Decay & Emission Analysis

Run the particle decay analysis to calculate penetration factors, deposition rates, and emission rates across 12 OPC-N3 size bins (0.35–10.0 µm):

# Run particle analysis (requires CO2 analysis results; per-event figures enabled by default)
python scripts/particle_decay_analysis.py

# Skip the per-event pm_decay figures
python scripts/particle_decay_analysis.py --no-plot

# Write to a custom directory without significant figure rounding
python scripts/particle_decay_analysis.py --output-dir /path/to/output --no-sig-figs

Methodology:

The mass balance equation for indoor particle concentration:

V·dC/dt = p·Q·C_out - Q·C - β·C·V + E
dC/dt = p·λ·C_out - λ·C - β·C + E/V

Rather than commit to one set of inputs, four emission-rate variants (E1-E4) are computed per bin, per the report's Table 7 / Equation 10 — each variant fixes a full combination of air-change-rate source, concentration series, and control volume:

Variant	Columns	λ source	Concentration series	Volume
E(t)1	bin{n}_outside_*	λ_outside	primary (raw C_bed1(t) / fleet C_room(t))	36.1086 m³ (bedroom)
E(t)2	bin{n}_entry_*	λ_entry	primary (same as E1)	36.1086 m³
E(t)3	bin{n}_adjusted_*	λ_outside	adjusted (ratio-corrected C_adjusted_room(t))	36.1086 m³
E(t)4	bin{n}_bathroom_*	λ_outside	primary (same as E1)	54.5868 m³ (bedroom + bathroom)

E1 is the value the ~45 water-temp/spray-pattern/door/etc. boxplots and comparison figures use. E2 brackets E1 against the entry-zone λ estimate. E3 is restricted to pre-cutover (single-monitor) events — post-cutover cells are NaN — since C_adjusted_room(t) corrects a single-sensor position bias that doesn't apply once the fleet is in place. E4 reuses E1's β (no volume term in the β fit) and only reruns the emission/Ct steps with the larger volume. See src/particle_emission_variants.py for the full per-bin implementation.

Room-concentration correction (src/particle_room_correction.py): the indoor concentration C above is not the raw MOD-PM-00195 reading for every variant. Two series are built, both sharing the same fleet-period replacement:

• 2026-06-03 to 2026-07-16 (the MODULAIR-PM fleet co-location window; see MODULAIR-PM Fleet Comparison), both series: C is replaced with C_room(t), the position-weighted fleet average, at every minute where more than 8 fleet sensors report; minutes at or below that threshold are left as NaN.

• Before 2026-06-03, primary series (build_raw_inside_data, feeds E1/E2/E4): C is left as the raw MOD-PM-00195 reading.

• Before 2026-06-03, adjusted series (build_corrected_inside_data, feeds E3 only): C is corrected inside each event's analysis window (shower_on − 15 min through deposition_end) as C_adjusted_room(t) = C_bed1(t) · C_room,average(minute) / C_bed1,average(minute), using the per-minute averaged curves scripts/moduair_cave_ratio.py computes from the fleet period, bucketed by the event's water temperature: < 38 °C uses the W24 ratio (n = 2 events), 38–41 °C inclusive uses W38-W41, and > 41 °C uses W49 (n = 2 events) — the only three buckets the fleet comparison has a ratio curve for. Outside any event's window, the adjusted series is left raw.

  The penetration factor (step 1) always uses the raw MOD-PM-00195/outside reading for the primary series, both before and after 2026-06-03, since its before/after windows fall outside the range the fleet ratio curves cover; E3 computes its own independent penetration factor on the adjusted series. particle_emission_variant_figures.py (workflow step 7) plots β and E per bin, per comparison pair, with a vertical line at 2026-06-03 marking this cutover.

1. Penetration factor (p): Ratio of indoor to outdoor concentration averaged over before/after windows relative to the shower event; zeros excluded; capped at 1.
   • QA: Minimum 10 valid points per window (MIN_POINTS_PENETRATION); windows where this is not met are skipped.
   • Night event windows: before = 9 pm (day before) → 2 am; after = 9 am → 2 pm
   • Day event windows: before = 9 am → 2 pm; after = 9 pm → 2 am (next day)
2. Air change rate (λ): λ_outside and λ_entry loaded from CO2 decay analysis results (h⁻¹); steps 3–6 repeat once per variant (E4 reuses E1's β; see the variant table above).
3. Other process rate (β, numerical): Step-by-step solution during the 2-hour (DEPOSITION_WINDOW_HOURS) post-shower decay window; beta solved at each time step from the discrete mass balance:

   β = 1/dt - λ - C_{t+1}/(C_t·dt) + p·λ·(C_out,t / C_t)
   

   All estimates ≤ MAX_OTHER_PROCESS_RATE (15 h⁻¹) are collected (no lower bound, to avoid upward bias). A 5th–95th percentile trim removes extreme outliers symmetrically. Beta is then selected via an R²-based four-step procedure (threshold R² = 0.80): (a) use unclamped trimmed mean — if R² ≥ 0.80 keep it; (b) clamp to ≥ 0 — if R² ≥ 0.80 keep it; (c) set β = 0 — if R² ≥ 0.80 keep it; (d) if β = 0 also fails → β = NaN (bin invalid, no Ct prediction plotted).
   • QA: Minimum 10 valid points in the decay window and after peak (MIN_POINTS_OTHER_PROCESS); MAX_OTHER_PROCESS_RATE = 15 h⁻¹ upper cap (above this, values are noise spikes for sub-3 µm particles).
4. Emission rate (E): Mass balance rearranged and solved numerically at each time step during the shower-on-to-peak window; mean and std reported from positive values.
   • QA: Minimum 3 valid consecutive steps (MIN_POINTS_EMISSION); E_total clips negative per-step contributions to 0 before trapezoid integration.
5. Predicted concentration (Ct): Single continuous forward Euler simulation from shower ON to 2 hours after shower OFF. E = E_mean before peak_time, E = 0 after. Returned as two connected segments forming one continuous predicted Ct curve on figures. Decay phase starts from the predicted peak (end of emission simulation), not the measured peak.
   • Emission R² (E_r_squared): R² of the emission-phase simulation vs. measured indoor concentration from shower ON to peak.
6. Total emission (E_total): Trapezoidal rule area under the E_t vs. time curve from shower ON to peak; negative per-step contributions clipped to 0.

Output Files (this script):

Sheets/columns marked "doubled" carry the E1/E2 variants, bin{n}_outside_... and bin{n}_entry_...; all_results additionally carries the E3 (bin{n}_adjusted_..., including its own independent bin{n}_adjusted_p_mean/p_std/peak_time) and E4 (bin{n}_bathroom_...) variants. The Excel workbook retains every event regardless of flow rate.

• output/particle_analysis_summary.xlsx - Multi-sheet Excel workbook:
  • all_results — Full results table (all metrics, all four variants, per event and bin)
  • p_penetration — Penetration factors per event and bin (includes test_name); primary series (E1/E2/E4 share this)
  • beta_deposition — Other process rates per event and bin (includes test_name); doubled
  • beta_r_squared — R² of forward Euler decay fit per event and bin (includes test_name); doubled
  • E_emission — Emission rates per event and bin (includes test_name); doubled
  • E_total_particles — Total emitted particle counts (E_total) per bin (includes test_name); doubled
  • E_r_squared — R² of forward Euler emission-phase simulation per bin (includes test_name); doubled
  • peak_comparison — Measured vs. predicted concentration at peak_time and deposition_end (wide format, one row per event); predicted columns doubled
• output/plots/event_figures/pm_decay/event_NN-testname_pm_decay.png - Individual event decay curves (four-panel), E1 (outside) variant only
  • Top panel: Measured concentration (solid) and continuous predicted Ct (dashed); shower ON/OFF dotted markers; optional outdoor PM overlay for events in OUTDOOR_PM_EVENTS
  • Emission panels (×3): Bins 0–2, 3–6, and 7–11; per-step E_t as faint lines; E_mean as dashed horizontal line; emission R² annotation; x-axes match concentration panel
• output/plots/event_figures/excluded_events/ - Figures for duration-excluded (water-temperature-testing) events

Output Files (particle_emission_variant_figures.py, workflow step 7):

Reads particle_analysis_summary.xlsx's all_results sheet; does not recompute p, β, or E. The boxplot/comparison figures below only include events with a measured flow rate in FLOW_RATE_MIN–FLOW_RATE_MAX (4.1–5.6 LPM), outside/entry (E1/E2) only.

• output/plots/particle/beta_other_{pair}_bin{N}.html, output/plots/particle/emission_{pair}_bin{N}.html - Per-bin variant comparison figures (interactive), for pair in entry_outside (E2 vs. E1), outside_adjusted (E3 vs. E1, pre-cutover only), outside_bathroom (E4 vs. E1); 12 bins × 2 metrics × 3 pairs = 72 figures. A vertical line at 2026-06-03 marks the room-correction cutover.
• output/plots/penetration_summary.png - Bar chart summary of penetration factors by size (λ-independent, one file)
• output/plots/deposition_summary_{outside,entry}.png - Bar chart summary of other process rates by size, one file per λ source
• output/plots/emission_summary_{outside,entry}.png - Bar chart summary of emission rates by size, one file per λ source
• output/plots/emission_etotal_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total by water temperature (fixed temp axis, doubled per source: 6 files)
• output/plots/emission_etotal_boxplot.md - Event membership for the E_total boxplot (lists event numbers and config keys per box)
• output/plots/other_process_rate_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - β by water temperature (fixed temp axis, doubled per source: 6 files)
• output/plots/other_process_rate_boxplot.md - Event membership for the β boxplot
• output/plots/emission_rate_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_mean by water temperature (fixed temp axis, doubled per source: 6 files)
• output/plots/emission_rate_boxplot.md - Event membership for the E_mean boxplot
• output/plots/penetration_factor_boxplot_{bin0-2,bin3-6,bin7-11}.png - p by water temperature (fixed temp axis; λ-independent, 3 files)
• output/plots/penetration_factor_boxplot.md - Event membership for the penetration factor boxplot
• output/plots/emission_etotal_by_bedroom_rh_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total vs. bedroom RH (metric axis, doubled per source)
• output/plots/emission_etotal_by_bedroom_rh_boxplot.md - Event membership for the bedroom RH boxplot
• output/plots/emission_etotal_by_bedroom_temp_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total vs. bedroom temperature (metric axis, doubled per source)
• output/plots/emission_etotal_by_bedroom_temp_boxplot.md - Event membership for the bedroom temp boxplot
• output/plots/emission_etotal_by_acr_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total vs. air change rate (metric axis, doubled per source)
• output/plots/emission_etotal_by_acr_boxplot.md - Event membership for the ACR boxplot
• output/plots/emission_etotal_by_beta_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total vs. other process rate (metric axis, doubled per source)
• output/plots/emission_etotal_by_beta_boxplot.md - Event membership for the β boxplot
• output/plots/emission_etotal_by_p_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total vs. penetration factor (metric axis, doubled per source)
• output/plots/emission_etotal_by_p_boxplot.md - Event membership for the p boxplot
• output/plots/emission_etotal_by_showerhead_boxplot_{bin0-2,bin3-6,bin7-11}_{outside,entry}.png - E_total by shower head type (doubled per source)
• output/plots/emission_etotal_by_showerhead_boxplot.md - Event membership for the shower head boxplot
• output/plots/comparison/spray_pattern_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by spray pattern (5 metric families × 3 bin groups each)
• output/plots/comparison/head_type_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by shower head type
• output/plots/comparison/mannequin_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by mannequin presence
• output/plots/comparison/bath_door_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by bathroom door position
• output/plots/comparison/bedroom_door_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by bedroom door position
• output/plots/comparison/fan_status_{metric}_{bin0-2,bin3-6,bin7-11}.png - Comparison boxplots by fan status

See the room-concentration correction note above for what the primary vs. adjusted concentration series feed.

Relative Humidity, Temperature & Wind Analysis

Run the environmental analysis to characterize conditions before and after shower events:

# Run RH/temperature/wind analysis (plots enabled by default)
python scripts/rh_temp_other_analysis.py

# Skip plots
python scripts/rh_temp_other_analysis.py --no-plots

# Plot only specific events
python scripts/rh_temp_other_analysis.py --events 1,3,5

# Limit time series plots to first 10 events
python scripts/rh_temp_other_analysis.py --max-plot-events 10

Methodology:

• Monitors RH and temperature at multiple locations:
  • Entry, Bedroom, Bathroom, Outside (Aranet4)
  • Bedroom, Bathroom, Living Room (DAQ/Vaisala)
  • Bathroom1, Bathroom2, Bath/Bed, Bedroom1, Bedroom2, Bedroom3 (HOBO UX100)
  • Note: QuantAQ met_rh / met_temp measure inside the flow cell and are excluded from RH/temp analysis
• Analyzes windspeed and direction from outdoor weather station (AIO2)
• Compares pre-shower baseline (30 min) to post-shower response (2 hours)
• Calculates mean, standard deviation, and min-max ranges
• Generates time series and box plot comparisons

Output Files:

• output/rh_temp_wind_summary.xlsx - Multi-sheet Excel workbook:
  • RH_Summary, Temp_Summary, WindSpeed_Summary, WindDir_Summary — per-sensor summary statistics
  • RH_Details, Temp_Details, WindSpeed_Details, WindDir_Details — per-event detail statistics
  • Event_Log — event metadata (test_name, config_key, shower times, door position, fan status)
  • Bedroom_Conditions — per-event bedroom RH and temperature (30-min pre-shower window); used by particle analysis for boxplot annotations; not rounded by sig figs
• output/plots/event_figures/rh_timeseries/event_NN-testname_rh_timeseries.png - RH time series per event
• output/plots/event_figures/temperature_timeseries/event_NN-testname_temperature_timeseries.png - Temperature per event
• output/plots/event_figures/wind_timeseries/event_NN-testname_wind_timeseries.png - Wind data per event
• output/plots/event_figures/excluded_events/ - Time series figures for duration-excluded events
• output/plots/rh_pre_post_boxplot.png - Pre/post RH comparison across all events
• output/plots/temperature_pre_post_boxplot.png - Pre/post temperature comparison
• output/plots/wind_speed_pre_post_boxplot.png - Pre/post wind speed comparison
• output/plots/wind_direction_pre_post_boxplot.png - Pre/post wind direction comparison

MODULAIR-PM Fleet Comparison (IAQ&V)

A separate line of analysis compares the study's two QuantAQ MODULAIR-PM sensors (inside MOD-PM-00195, outside MOD-PM-00785) against the broader NIST Indoor Air Quality and Ventilation Group fleet of co-located MODULAIR-PM units deployed in the same test facility. It characterizes sensor agreement, builds candidate inter-sensor correction factors for the historical 195 record, and evaluates a position-weighted fleet average (C_room) as an alternative bedroom concentration reference.

This is no longer a purely exploratory side analysis. moduair_cave_ratio.py's compute_room_frame / group_average_curve functions are reused in-memory by src/particle_room_correction.py to build the room-corrected "inside" concentration the main particle pipeline analyzes (see the room-concentration correction note under Particle Decay & Emission Analysis). That means particle_decay_analysis.py (workflow step 6) now also requires the moduair_pm_chunks share below to be reachable — but it does not require running moduair_cave_ratio.py as a script first.

Fleet data are produced by the separate NIST_moduair-pm repository and read from the moduair_pm_chunks share configured in data_config.json — this is independent of the QuantAQ API data downloaded by download_quantaq_data.py. src/moduair_loader.py loads and aligns fleet sensor chunks onto a shared 1-minute time base; it does not download anything itself. All five scripts below share a 2026-06-0{3,4} through 2026-07-16 default window, the period the fleet was co-located with the study sensors, and take --start/--end/ --output-dir overrides.

# Inter-sensor correction factor ratios + 195/813 Deming fit
python scripts/moduair_correction_factor.py

# Per-event, per-sensor delta peak time (shower ON to peak)
python scripts/moduair_event_peak_times.py

# Per-bin time series across the co-located bedroom fleet
python scripts/moduair_bin2_timeseries.py

# Inside sensor (195) vs. position-weighted fleet average (C_room)
python scripts/moduair_cave_ratio.py

# Per-sensor variant of export_config_timeseries.py across the fleet
python scripts/export_config_timeseries_fleet.py [--no-sig-figs] [--sensors 195 401 ...]
python scripts/export_config_timeseries_fleet.py --c-room

Script	Purpose	Key Output
moduair_correction_factor.py	Bin-wise concentration ratios (195/813, 195/quad-average, 195/others) plus a Deming (orthogonal-distance) regression fit of 813 vs. 195 for a candidate correction	moduair_correction_factor_195_813_*.csv, moduair_correction_195_813_fit_*.csv, per-bin HTML figures
moduair_event_peak_times.py	Delta peak time (shower ON to peak) per sensor per event	moduair_event_peak_times.csv, event_delta_peak_times.html
moduair_bin2_timeseries.py	One interactive figure per bin (0–10) of the co-located bedroom fleet, styled after the NIST_moduair-pm weekly figures	plots/moduair_correction/bin{N}_timeseries.html
moduair_cave_ratio.py	C_bed1 (sensor 195) vs. C_room (position-weighted fleet average) scatter and time-ratio figures, restricted to samples where >8 fleet sensors report	plots/moduair_room/c_bed1_vs_c_room_bin{N}*.html, moduair_room_ratio_fit*.csv
export_config_timeseries_fleet.py	Per-sensor and C_room aggregated + raw event-configuration workbooks, same structure as export_config_timeseries.py	output/event_config_timeseries_fleet/MOD-PM-<sn>/..., .../C_room/...

Dependencies

Key packages (see environment.yaml for complete list):

• pandas, numpy - Data manipulation
• scipy - Numerical integration (trapezoidal rule for E_total calculation)
• matplotlib - Publication-quality figure generation
• bokeh - Interactive visualization
• requests - API communication
• python-dotenv - Environment variable management (.env file loading)
• pyyaml - Configuration management
• openpyxl - Excel file I/O

References

Instrument manuals and data sheets are not included in this repository due to copyright. Public download links are listed in docs/pdf_links.md.

Aranet4 PRO (CO2 Sensor)

• Aranet4 PRO Data Sheet
• Aranet4 User Manual

Met One AIO 2 (Weather Sensor)

• AIO 2 Manual (Rev G)
• AIO 2 Data Sheet

DAQ System (NI cDAQ-9178)

• NI cDAQ-9171/9174/9178 User Manual
• NI 9201 Getting Started Guide

Setra Model 264 (Differential Pressure)

• Setra Model 264 Data Sheet

Vaisala Humidity & Temperature Probes

• HMP155 User Guide
• HMP45A/D User Guide

---

Contact

• PI: Dustin G. Poppendieck — dustin.poppendieck@nist.gov
• Maintainer: Nathan Lima — nathan.lima@nist.gov
• NIST Organizational Unit: Engineering Laboratory
• Division: Building Energy and Environment Division
• Group: Indoor Air Quality and Ventilation Group

Data Availability

The experimental data associated with this project are not included in this repository. Data can be made available upon request by contacting the PI.

Citation

If you use this software, please cite it as:

@software{nist_epa_legionella,
  author       = {Lima, Nathan M. and Poppendieck, Dustin G.},
  title        = {NIST-EPA Legionella Study Analysis Tools},
  year         = {2026},
  publisher    = {National Institute of Standards and Technology},
  doi          = {10.18434/mds2-4153},
  url          = {https://github.com/usnistgov/NIST_EPA_Legionella}
}

Related Publication

A peer-reviewed publication and/or NIST Technical Note associated with this study is in preparation. This section will be updated with a citable BibTeX entry and DOI when available.

License

This software is in the public domain and provided "as is" according to the LICENSE.md file. It was developed by employees of the National Institute of Standards and Technology (NIST), an agency of the Federal Government and is being made available as a public service. Pursuant to Title 17 United States Code Section 105, works of NIST employees are not subject to copyright protection in the United States. This software may be subject to foreign copyright. Permission in the United States and in foreign countries, to the extent that NIST may hold copyright, to use, copy, modify, create derivative works, and distribute this software and its documentation without fee is hereby granted on a non-exclusive basis, provided that this notice and disclaimer of warranty appears in all copies.

See LICENSE.md for the full NIST licensing statement.