Processing many files using SLURM jobs

Author: Philipp Schreiner
Created: Aug. 29, 2025
Last updated: Jul. 2, 2026 by Philipp Schreiner


Defining the elements of an analysis on a single file and scaling it up to many files can be done in many different ways (mostly depending on personal preference).

Here, we show a version where you start from a notebook, send out the tasks to the cluster, and collect the results in the same notebook again, for you to continue your analysis. This approach easily scales and you can add more files by just extending a list of file names.

We assume that at this point, you already know how to create SEVs and OFs, apply quality cuts and calculate spectra, and that you just want to trigger many files and perform those processing steps on them. Of course, the principle can be applied to all sorts of different tasks, so feel free to adapt the scripts to your needs.

Tip

Even though the procedure might look intimidating at first, we recommend to start following it as early as possible because it automatically results in a clean notebook/script structure even for many analysis files.

To get started, make sure to copy the Files below to your working directory, together with the analysis notebook my_analysis.ipynb that we will be working in:

- analysis_info.py
- my_analysis.ipynb
- submit.py
- workflow.py

We also recommend to create thre folders files, job_files, and hdf5, to keep things organized.

Check the files that you copied for user specific variables (e.g. paths of streams, user names, container names, etc.) and update them.

The most important one is analysis_info.py. Here, all analysis information like file paths, channel names, record length, trigger configurations, etc. are saved in one place.

Warning

This tutorial serves as an illustration. You cannot just run it like the previous ones.

Tip

Here, we work in a single notebook, my_analysis.ipynb. What the author of this tutorial prefers is to have two notebooks. One for template/filter creation (e.g. 0_templates_and_filters.ipynb) and one for the main analysis (e.g. 1_analysis.ipynb). They share and use the same analysis_info.py file, which grows as the analysis develops. See below.

Notebook

import os
import numpy as np
import scipy as sp

import cait as ai
import cait.versatile as vai
# The analysis_info file stores things like the record length that we want to 
# use, and the SEVs/NPSs/OFs if we created them already. We load them here.
from analysis_info import fdir, fdirh5, ch_name, record_length, OF, SEV, NPS

# Have a look at analysis_info.py! You'll find a list of file names for the
# template creation and the analysis. There, you can easily extend the list
# even after you started this notebook to scale your analysis.
from analysis_info import fnames_analysis as fnames

n_channels = 2

# This is just to keep things organised.
# I usually format my HDF5 file names in that way.
# For the appendix, I use things like '_templates' for when I'm creating SEVs, etc.
# Feel free to adjust it.
detector_name = "myDet"
run = "042"
hdf5_appendix = ""
hdf5_name = f"{detector_name}_run{run}{hdf5_appendix}"
# Here we create a dictionary of all the stream files to organize them naicely.
# The exact way you do it depends on the stream hardware. For VDAQ2 files, the one
# below works quite well
streams = {fname.split('.')[0]: vai.Stream("vdaq2", os.path.join(fdir, fname)) for fname in fnames}
streams
# We also create a dictionary of DataHandlers where each entry corresponds to one stream file
dhs = {
    fname: ai.DataHandler(
        record_length=record_length, 
        nmbr_channels=len(channel_names)+len(passive_names), 
        sample_frequency=s.sample_frequency,
    )
    for fname, stream in streams.items()
}

# Initialize the empty DataHandlers
for k, v in dhs.items():
    v.set_filepath(path_h5=fdirh5, fname=f"{k}_{detector_name}{hdf5_appendix}", appendix=False)
    v.init_empty()

Preprocessing

We now have (Stream, DataHandler)-pairs and we want to perform a number of operations on each pair (and save the results in the DataHandlers). What operations to perform is defined in the STAGES_CONFIG dictionary in the analysis_info.py file (have a look below!). We see triggering stages, one for main parameter calculation (CMP) and optimum filter/template fit amplitude reconstruction.

Submitting to the cluster

Let’s say we want to run stages ["OF_TRIGGER", "CMP", "REAPPLY_OF", "TEMPLATE_FIT"]. For that, we can use the define_job function from the workflow.py script. After calling it, a .json file containing information for submitting SLURM jobs is created.

The fname will be the name of the .json file. We recommend to store it in the job_files folder. The stages must be valid keys in the STAGES_CONFIG dictionary in the analysis_info.py file. If you want to run the job only on a subset of Streams and DataHandlers, you can also filter them here. E.g. if you have bck and ncal files, you can define a job for ncal files only by doing something like streams=[s for name, s in streams.items() if 'ncal' in name].

from workflow import define_job

define_job(
    fname="job_files/of_trigger_cmp_of_tf", 
    stages=["OF_TRIGGER", "CMP", "REAPPLY_OF", "TEMPLATE_FIT"],
    streams=list(streams.values()),
    dhs=list(dhs.values())
)

You can then log into your computing cluster, go to

!pwd

From here, we run

python3 submit.py job_files/of_trigger_cmp_of_tf.json

We can check if the jobs are running using squeue -u $USER or slurm q and what they are doing using sattach <jobid>.0.

Additionally, job logs will be written to job_files/logs and you must check them for problems. If all jobs completed (i.e. they disappear from slurm q) and they didn’t stop with an error (as seen from the logs), you can move on to the next step.

See also

Saving the DataHandler and Stream in .json files is done through a process called serialization. cait implements its own serialization mechanism which you can read more about in the documentation of cait.serialize.

Running in notebook (not on cluster):

Alternatively (if you don’t want to use SLURM jobs), you can just run all calculations on the files sequentially in the notebook (otherwise this loop would just be run in parallel on the cluster, the result is identical):

from workflow import job

for dh, s in zip(dhs.values(), streams.values()):
    print(f"Working on file '{dh.fname}' ...")
    job(stages, s, dh)

Saving file information

Sometimes, you have different ‘kinds of files’. E.g., background data files, 57Co calibration files, neutron calibration files, etc. You can either create separate DataHandlers for them, or (as preferred by the author of this tutorial) you can just have a single DataHandler and save flags that let you select the respective files easily.

for i, (name, dh) in enumerate(dhs.items()):
    for g in [
        "events", 
        "testpulses", 
        "controlpulses",
        "noise",
    ]:
        bck_flag = np.zeros(dh[f'{g}/hours'].shape[-1], dtype=bool)
        ncal_flag = np.zeros(dh[f'{g}/hours'].shape[-1], dtype=bool)
        file_index = i*np.ones(dh[f'{g}/hours'].shape[-1], dtype=np.int32)
        
        if name.startswith("bck"):
            bck_flag[:] = True
        elif name.startswith("ncal"):
            ncal_flag[:] = True
        else:
            raise NotImplementedError

        dh.set(g, flag_bck=bck_flag, flag_ncal=ncal_flag, dtype=bool, overwrite_existing=True)
        dh.set(g, file_index=file_index, dtype=np.int32, overwrite_existing=True)

Combining the processed HDF5 files

Running the job produced an HDF5 file for every (Stream, DataHandler)-pair. We now combine them to do analysis on a single DataHandler.

# Combine HDF5 files into a single one for analysis
ai.data.combine_h5(
    fname=hdf5_name, 
    files=[dh.fname for dh in dhs.values()], 
    src_dir=fdirh5, 
    out_dir=fdirh5, 
    groups_combine=[
        "events", 
        "testpulses", 
        "noise", 
        "controlpulses",
    ],
) 

Start Analysis

Now that you have a single HDF5 file, you can continue your analysis as usual by creating a DataHandler instance:

dh = ai.DataHandler(
    record_length=record_length, 
    nmbr_channels=len(channel_names) + len(passive_names), 
    sample_frequency=list(streams.values())[0].sample_frequency,
)
dh.set_filepath(path_h5=fdirh5, fname=hdf5_name, appendix=False)
print(dh)

From there, you can define cuts, plot spectra, etc.

Note

We recommend to create two such notebooks. One for template/filter creation (e.g. 0_templates_and_filters.ipynb) and one for the main analysis (e.g. 1_analysis.ipynb). In the first one, the goal is just to save SEV, NPS, OF to .xy files and to estimate the baseline resolution, respectively threshold to be used for the of-trigger. Then, one starts with the of-trigger in the second notebook for the main analysis.

The working directory would then look something like this:

files
    - OF.xy
    - NPS.xy
    - SEV.xy
hdf5
    - filter_creation_file0.h5
    - filter_creation_file1.h5
    - filter_creation.h5
    - main_analysis_file0.h5
    - main_analysis_file1.h5
    - main_analysis_file2.h5
    - main_analysis_file3.h5
    - main_analysis.h5
job_files
    - trigger_templates.json
    - trigger_of.json
    - template_fit_and_reapply_of.json

0_templates_and_filters.ipynb
1_analysis.ipynb
analysis_info.py
submit.py
workflow.py

Note

Note that you don’t have to define everything in the analysis_info.py script from the beginning. Usually, one starts with the z-score trigger and only runs stages=["Z-SCORE_TRIGGER", "CMP"]. Then one creates SEV, NPS, OF and fills in the of-trigger information, then runs it again for stages=["OF_TRIGGER", "CMP"]. After validating which parameters to use for dh.apply_ofilter and dh.apply_template_fit (using the preview feature), one fills it into analysis_info.py and runs stages=["REAPPLY_OF", "TEMPLATE_FIT"], and so on. In the end, all the relevant information is nicely collected in analysis_info.py. And if someone asks you which trigger threshold you used in your analysis, you know exactly where to find it.

Tip

Often, the analysis on a few files is already done while more files are being written. The approach presented above is very convenient for that situation because you just have to extend the fnames list in analysis_info.py and run all stages again. Since the stages check whether the information has been computed before, only the new files are processed. Once the new files are processed, you combine the HDF5 files again and work with the same DataHandler as before.

For this to be as convenient as possible, we recommend to do all costly calculations on the single DataHandlers individually, not on the combined DataHandler (since those results would be lost once you combine HDF5 files again). An example that loops through all single files is given in Files.

Files

Below, you find the script files that collect the analyis configuration and the helper functions for cluster job submission.

analysis_info.py:

import numpy as np
import cait.versatile as vai

# BASIC FILE INFO (added immediately)
fdir = "/where/the/stream/files/are/located/"
fdirh5 = "/your/working/directory/"

# Files used for template/filter creation
fnames_creation = [
    "bck_001",
    "bck_002",
]

# Files that are included in the analysis (in temporally increasing order)
fnames_analysis = [
    "bck_001",
    "bck_002",
    "bck_003",
    "ncal_001",
    "ncal_002",
    "ncal_003",
]

# BASIC TRIGGER INFO (added immediately)
ch_name = {
    "Phonon": "ADC1",
    "Phonon_TP": "DAC1",
    "Light": "ADC2",
    "Light_TP": "DAC2",
}
# Record length and random noise sampling frequency for triggering
record_length = 2**14
f_noise = 1000 # per hour

# CONTROLPULSE INFO (added after inspecting the available TPAs)
CP_above = {
    "Phonon": 10, # V
    "Light": 10, # V
}

# ENERGY RESOLUTIONS (added after calculating them in the template creation notebook)
bl_res = {
    "Phonon": 0.123e-3, # V
    "Light": 1.234e-3, # V
}

# TRUNCATION LIMITS 
# (added after determining it in the template creation notebook)
truncation_limit = {
    "Phonon": 0.8, # V
    "Light": 1.5, # V
}

# SEVs (added after building and saving them in the template creation notebook)
SEV = {
    "Phonon": vai.SEV.from_file("files/SEV_phonon"),
    "Phonon_TP": vai.SEV.from_file("files/SEV_phonon_TP"),
    "Light": vai.SEV.from_file("files/SEV_light"),
    "Light_TP": vai.SEV.from_file("files/SEV_light_TP"),
}


# NPSs (added after building and saving them in the template creation notebook)
NPS = {
    "Phonon": vai.NPS.from_file("files/NPS_phonon"),
    "Light": vai.NPS.from_file("files/NPS_light"),
}

# OFs (added after building and saving them in the template creation notebook)
OF = {
    "Phonon": vai.OF.from_file("files/OF_phonon"),
    "Phonon_TP": vai.OF.from_file("files/OF_phonon_TP"),
    "Light": vai.OF.from_file("files/OF_light"),
    "Light_TP": vai.OF.from_file("files/OF_light_TP"),
}

# JOB STAGES (added and extended as the analysis progresses)
STAGES_CONFIG = {
    # ALWAYS DONE AFTER TRIGGERING (added immediately)
    # Needs arguments for dh.cmp()
    "CMP": [
        { "group": "events" },
        { "group": "testpulses" },
        { "group": "noise" },
        { "group": "controlpulses" },
    ],
    # FIRST ROUND OF TRIGGERING 
    # (added after 'controlpulses_above' is known)
    # Needs arguments for dh.trigger_zscore()
    "Z-SCORE_TRIGGER": {
        "trigger_channels": [ch_name["Phonon"]],
        "passive_channels": [ch_name["Light"]],
        "thresholds": [5], # sigma
        "testpulse_channels": [
            ch_name["Phonon_TP"], 
            ch_name["Light_TP"],
        ],
        "controlpulses_above": [
            CP_above["Phonon"], 
            CP_above["Light"],
        ],
        "copy_events": False,
        "f_noise": f_noise,
        "reuse_triggers": False,
    },
    # SECOND ROUND OF TRIGGERING 
    # (added after building templates/filters in the template creation notebook)
    # Needs arguments for dh.trigger_of()
    "OF_TRIGGER": {
        "trigger_channels": [ch_name["Phonon"]],
        "passive_channels": [ch_name["Light"]],
        "of": OF["Phonon"], 
        "thresholds": [bl_res["Phonon"]*5],
        "testpulse_channels": [
            ch_name["Phonon_TP"], 
            ch_name["Light_TP"],
        ],
        "controlpulses_above": [
            CP_above["Phonon"], 
            CP_above["Light"],
        ],
        "copy_events": False,
        "f_noise": f_noise,
        "reuse_triggers": False
    },
    # REAPPLYING THE FILTER
    # (added after building templates/filters in the template creation notebook)
    # Needs arguments for dh.apply_ofilter()
    "REAPPLY_OF": [
        {
            "group": "events",
            "of": OF["Phonon"],
            "sev": SEV["Light"],
            "on_stream": True,
            "tag": "uncorrelated",
        },
        {
            "group": "events",
            "of": [OF["Phonon"], OF["Light"]],
            "sev": [SEV["Phonon"], SEV["Light"]],
            "max_search": [(0.2, 0.4), (-38, -28)],
            "relative_to": [None, 0],
            "on_stream": True,
            "tag": "correlated",
        },
        # ... add more configurations as needed
    ],
    # TEMPLATE FIT
    # (added after building templates/filters in the template creation notebook)
    # Needs arguments for dh.apply_template_fit()
    "TEMPLATE_FIT": [
        {
            "group": "events",
            "sev": [SEV["Phonon"], SEV["Light"]],
            "bl_poly_order": [3, 1], 
            "truncation_limit": [truncation_limit["Phonon"], truncation_limit["Light"]],
            "correlated": True,
            "fit_onset": [True, False],
        },
        # ... add more configurations as needed
    ]
}

workflow.py:

from typing import List
import json

import cait as ai
import cait.versatile as vai

from analysis_info import STAGES_CONFIG

# Function that checks if all requested stages are present in the STAGES_CONFIG
def _check_stages(stages: List[str]):
    if not all([stage in STAGES_CONFIG.keys() for stage in stages]):
        raise KeyError(f"All selected stages must be present in the configured STAGES dictionary. Available: {list(STAGES.keys())}. Got: {stages}.")

def define_job(
    fname: str, 
    stages: List[str], 
    streams: List[vai.stream.streambase.StreamBaseClass], 
    dhs: List[ai.DataHandler],
):
    """
    Save a json-file with 'fname' (omit .json extension) that defines a job consisting of stage names 'stages' and is applied to each (Stream, DataHandler)-pair in 'streams'/'dhs'.
    
    Jobs can be submitted via 'python3 submit.py fname.json'
    """
    _check_stages(stages)
    d = {
        "stages": stages,
        "streams": [s.to_dict() for s in streams],
        "dhs": [d.to_dict() for d in dhs]
    }
    
    with open(f"{fname}.json", "w") as f:
        json.dump(d, f, indent=4)

def job(
    stages: List[str], 
    stream: vai.stream.streambase.StreamBaseClass, 
    dh: ai.DataHandler,
):
    """Perform each analysis step defined in 'stages' to the (Stream, DataHandler)-pair."""
    _check_stages(stages)
        
    # In each stage, we check, whether the dataset is available already
    # and only calculate it, if it's not
    for stage in stages:
        print(ai.styles.txt_fmt(f"Entering stage {stage} ...", color="yellow", style="bold"))
        
        if stage.startswith("Z-SCORE_TRIGGER"):
            if not dh.exists(f"event_building-z-score"):
                dh.trigger_zscore(stream=stream, **STAGES_CONFIG[stage])
                
        if stage.startswith("OF_TRIGGER"):
            if not dh.exists(f"event_building-of"):
                dh.trigger_of(stream=stream, **STAGES_CONFIG[stage])
                
        if stage.startswith("CMP"):
            for substage in STAGES_CONFIG[stage]:
                has_events = True
                try: 
                    # Note that the check 'dh.exists("group/events")' is unreliable as the events
                    # could be virtually stored. The easiest check is therefore to see if
                    # get_event_iterator() fails
                    dh.get_event_iterator(substage['group'])
                except: 
                    has_events = False
                if has_events:
                    if not dh.exists(f"{substage['group']}/pulse_height"):
                        dh.cmp(**substage)
                    
        if stage.startswith("REAPPLY_OF"):
            for substage in STAGES_CONFIG[stage]:
                appendix = substage.get("tag", "")
                appendix = "-" + appendix if appendix else ""
                print(f"{substage['group']=}, {appendix=}")
                if dh.exists(f"{substage['group']}") and not dh.exists(f"{substage['group']}/of_ph{appendix}"):
                    dh.apply_ofilter(**substage)
                    
        if stage.startswith("TEMPLATE_FIT"):
            for substage in STAGES_CONFIG[stage]:
                appendix = substage.get("tag", "")
                appendix = "-" + appendix if appendix else ""
                print(f"{substage['group']=}, {appendix=}")
                if dh.exists(f"{substage['group']}") and not dh.exists(f"{substage['group']}/templatefit_pars{appendix}"):
                    dh.apply_template_fit(**substage)
                    
        # YOU CAN ADD CUSTOM STAGE HANDLING BELOW. 
        # Add the respective entry in the STAGES_CONFIG dictionary in analysis_info.py and handle
        # it as demonstrated above for 'CMP', 'REAPPLY_OF', ...
        
        # ...
        
        print(ai.styles.txt_fmt(f"DONE with stage {stage} ..", color="green", style="bold"))
        
if __name__ == "__main__":
    import os
    import sys
    import json
    
    print(f"Received script arguments: {sys.argv[1:]}")
    
    json_path = sys.argv[1]
    
    try:
        with open(json_path, "r") as f:
            job_dict = json.load(f)

        print(f"Running job defined by dict: {job_dict}")

        job(
            job_dict["stages"], 
            ai.serialize.load(job_dict["stream"]),
            ai.serialize.load(job_dict["dh"]),
        )
        
    finally:
        print("Cleaning up temporary file ...")
        try:
            os.remove(json_path)
            print(f"Deleted temp file {json_path}")
        except Exception as e:
            print(f"Could not delete temp file {json_path}: {e}")

submit.py:

import os
import subprocess

# The analysis directory (where the workflow.py, analysis_info.py and submit.py files are).
DIR = os.path.dirname(os.path.abspath(__file__))

# We will write temporary files and job logs into this directory
JOB_FILES_DIR = f"{DIR}/job_files"

# Path of the container
CONTAINER = "/path/to/cait_develop-slim.sif"

# Script to be executed in the jobs
SCRIPT = "workflow.py"

# Job name in the 'squeue -u $USER' list
JOB_NAME = "cait-analysis"

# RESOURCES
TIME, CORES, MEMORY = "02:00:00", 4, "16G"

# Validation function for a dictionary that defines a job
def _check_valid_json(json_dict: dict):
    for k in ["stages", "streams", "dhs"]:
        if k not in json_dict.keys():
            raise KeyError(f"JSON files defining a valid job need to have a '{k}' key.")

        if not isinstance(json_dict[k], list):
            raise TypeError(f"JSON files defining a valid job need to have a '{k}' key which contains a list.")
            
    if len(json_dict["streams"]) != len(json_dict["dhs"]):
        raise ValueError(f"Lists of streams and DataHandlers have to have same length. Got {len(json_dict['streams'])} and {len(json_dict['dhs'])}.")
        
# Function that submits the job to the cluster's queue (for each (Stream, DataHandler)-pair separately; temp_json_files are created below)
def submit_job(job_tag: str, temp_json_file: str):
    command = f"srun singularity run {CONTAINER} python {SCRIPT} {temp_json_file}"
    
    sbatch_command = [
        "sbatch",
        f"--job-name={JOB_NAME}-{job_tag}",
        f"--chdir={DIR}",
        f"--output={JOB_FILES_DIR}/logs/{job_tag}.out",
        f"--error={JOB_FILES_DIR}/logs/{job_tag}.err",
        f"--time={TIME}",
        f"--mem={MEMORY}",
        f"--cpus-per-task={CORES}",
        "--partition=c",
        "--qos=c_medium",
        "--parsable",
        "--wrap", command
    ]
    
    print(f"Job command: {sbatch_command}")
    
    subprocess.run(sbatch_command, check=True)

if __name__ == "__main__":
    import sys
    import json
    import tempfile
    
    # Load json information
    json_path = sys.argv[1]
    # Name of json file is used to make the job more
    # recognizable (job name and log folder)
    job_tag = os.path.basename(json_path).split(".")[0]
    
    os.makedirs(JOB_FILES_DIR, exist_ok=True)
    
    # Load the .json file that defines all job stages for all (Stream, DataHandler)-pairs
    with open(json_path, "r") as f:
        job_dict = json.load(f)
    
    # Validate it
    _check_valid_json(job_dict)
    stages = job_dict["stages"]
    
    # Submit a job for each (Stream, DataHandler)-pair in the json file.
    # For that purpose, a smaller json file (per job) is created to
    # transfer them to the job script more easily.
    for i, (s, dh) in enumerate(zip(job_dict["streams"], job_dict["dhs"])):
        with tempfile.NamedTemporaryFile(mode="w", delete=False, suffix=".json", dir=JOB_FILES_DIR) as tmpfile:
            json.dump({"stages": stages, "stream": s, "dh": dh}, tmpfile)
            tmpfile_path = tmpfile.name
            
            # Ensure the file is readable (some clusters might restrict tmp file access)
            os.chmod(tmpfile_path, 0o644)
 
        submit_job(f"{job_tag}-{i}", tmpfile_path)