diff --git a/Tracking/CMakeLists.txt b/Tracking/CMakeLists.txt index 3a6654d50..f6099c849 100644 --- a/Tracking/CMakeLists.txt +++ b/Tracking/CMakeLists.txt @@ -44,3 +44,6 @@ setup_library(module Tracking sources ${SRC_FILES}) setup_python(package_name LDMX/Tracking) + +# Install data files (e.g. DAQ maps) to install/data/Tracking. +setup_data(module Tracking) diff --git a/Tracking/data/daqmap_esa25_slice_test.json b/Tracking/data/daqmap_esa25_slice_test.json new file mode 100644 index 000000000..66577e919 --- /dev/null +++ b/Tracking/data/daqmap_esa25_slice_test.json @@ -0,0 +1,57 @@ +{ + "_comment": [ + "Tracker DAQ map for the ESA slice test (Detectors/data/ldmx-reduced-v3).", + "Maps (feb_id, hybrid_id) -> Acts surface id + strip transform. All four", + "hybrids are serviced by a single FEB (0). Layer ordering, per hardware:", + " station 1 = hybrid 0 (stereo) + hybrid 1 (axial)", + " station 2 = hybrid 3 (stereo) + hybrid 2 (axial)", + "In each station the stereo sensor sits upstream (more negative z) and the", + "axial sensor downstream.", + "layer_id is the Acts surface id, keyed as vol*1000 + layer*100 + sensor", + "(TrackingGeometry.cxx:327); recoil volume = 3. The sensor suffix (00/01)", + "follows z-order, so the upstream stereo sensor is xx00 and the downstream", + "axial one is xx01." + ], + "sensors": [ + { + "feb": 0, + "hybrid": 0, + "station": 1, + "orientation": "stereo", + "layer_id": 3100, + "n_strips": 640, + "first_strip": 0, + "reversed": false + }, + { + "feb": 0, + "hybrid": 1, + "station": 1, + "orientation": "axial", + "layer_id": 3101, + "n_strips": 640, + "first_strip": 0, + "reversed": false + }, + { + "feb": 0, + "hybrid": 2, + "station": 2, + "orientation": "axial", + "layer_id": 3201, + "n_strips": 640, + "first_strip": 0, + "reversed": false + }, + { + "feb": 0, + "hybrid": 3, + "station": 2, + "orientation": "stereo", + "layer_id": 3200, + "n_strips": 640, + "first_strip": 0, + "reversed": false + } + ] +} diff --git a/Tracking/exampleConfigs/README.md b/Tracking/exampleConfigs/README.md index bc8a59adb..0d0e1e905 100644 --- a/Tracking/exampleConfigs/README.md +++ b/Tracking/exampleConfigs/README.md @@ -78,10 +78,11 @@ Step 2 — decode → subtract → build waveforms (per physics run) -------------------------------------------------------------- ```bash -just fire Tracking/exampleConfigs/decode_to_waveforms.py -- \ +just fire Tracking/exampleConfigs/raw_to_measurements.py -- \ --dat /sdf/data/hps/users/mgignac/hardware/data/LDMX/Run_182_20251213_174649.dat \ --pedestal-file pedestals.json \ --max-events 5 \ + --stop-at fit \ --output tracker_waveforms_run182.root # → [TrackerPedestalProvider] Loaded 2560 channel pedestals from 'pedestals.json' # → [SiStripWaveformBuilder] Built N waveforms (>=4 samples @5s, streak>=5 @3s) from 25800 hits diff --git a/Tracking/exampleConfigs/compute_pedestals.py b/Tracking/exampleConfigs/compute_pedestals.py index 421da5c52..91f2ab46c 100644 --- a/Tracking/exampleConfigs/compute_pedestals.py +++ b/Tracking/exampleConfigs/compute_pedestals.py @@ -2,13 +2,13 @@ This is step 1 of the real-data tracker waveform chain: - 1. compute_pedestals.py (this script) -> pedestals.json - 2. decode_to_waveforms.py -> .root (TrackerWaveforms) - 3. plot_waveforms.py -> PNG figures + 1. compute_pedestals.py (this script) -> pedestals.json + 2. raw_to_measurements.py --stop-at fit -> .root (TrackerWaveforms) + 3. plot_waveforms.py -> PNG figures Runs SingleSubsystemUnpacker -> RawTrackerDecoder -> PedestalCalculator on a baseline (pedestal) run, computing per-channel, per-sample mean and RMS noise, -and writes them to a JSON file. Pass that JSON to decode_to_waveforms.py for +and writes them to a JSON file. Pass that JSON to raw_to_measurements.py for physics runs. Usage diff --git a/Tracking/exampleConfigs/decode_to_waveforms.py b/Tracking/exampleConfigs/decode_to_waveforms.py deleted file mode 100644 index 64b53d121..000000000 --- a/Tracking/exampleConfigs/decode_to_waveforms.py +++ /dev/null @@ -1,143 +0,0 @@ -"""Decode tracker Rogue data, subtract pedestals, and build per-channel waveforms. - -This is step 2 of the real-data tracker waveform chain: - - 1. compute_pedestals.py -> pedestals.json - 2. decode_to_waveforms.py (this script) -> .root (TrackerWaveforms) - 3. plot_waveforms.py -> PNG figures - -Pipeline --------- - SingleSubsystemUnpacker -> TrackerRawData (bytes) - RawTrackerDecoder -> RawSiStripHits - PedestalSubtractor -> TrackerHits (RawSiStripHit, pedestal-subtracted) - SiStripWaveformBuilder -> TrackerWaveforms (SiStripWaveform) - -Pedestals/noise are supplied as a conditions object (TrackerPedestals) by -TrackerPedestalProvider, not stored in the hits. Run compute_pedestals.py first -to produce the pedestal JSON, then pass it here via --pedestal-file (it is handed -to the provider). The output ROOT file contains a SiStripWaveform collection -('TrackerWaveforms') ready for waveform analysis / plotting. - -Usage ------ - ldmx fire Tracking/exampleConfigs/decode_to_waveforms.py \ - [-- --dat /path/to/physics.dat] \ - [-- --pedestal-file pedestals.json] \ - [-- --max-events N] \ - [-- --output tracker_waveforms.root] -""" - -import argparse -import sys - -# just fire passes '--' through to the script; strip it so argparse sees only flags. -sys.argv = [a for a in sys.argv if a != "--"] - -parser = argparse.ArgumentParser(f"ldmx fire {sys.argv[0]}") -parser.add_argument( - "--dat", - default="/sdf/data/hps/users/mgignac/hardware/data/LDMX/Run_037_20251210_145218.dat", - help="Path to the physics Rogue .dat file", -) -parser.add_argument( - "--pedestal-file", - default="pedestals.json", - help="Path to the pedestal JSON file from compute_pedestals.py", -) -parser.add_argument( - "--max-events", - type=int, - default=10, - help="Maximum number of physics events to process (default: 10)", -) -parser.add_argument( - "--frame-offset", - type=int, - default=0, - help="Skip this many tracker frames at the start of the file (default: 0)", -) -parser.add_argument( - "--output", - default="tracker_waveforms.root", - help="Output ROOT file (default: tracker_waveforms.root)", -) -parser.add_argument( - "--verbose-waveforms", - action="store_true", - help="Drop the waveform builder to the trace logging level, printing the " - "per-channel fit results and full ASCII waveform traces", -) -parser.add_argument( - "--high-threshold", - type=float, - default=5.0, - help="Per-sample significance (ADC/noise) for the high-threshold count cut", -) -parser.add_argument( - "--min-high-samples", - type=int, - default=4, - help="Min samples that must exceed high_threshold", -) -parser.add_argument( - "--low-threshold", - type=float, - default=3.0, - help="Per-sample significance (ADC/noise) for the consecutive-streak cut", -) -parser.add_argument( - "--min-consecutive-low", - type=int, - default=5, - help="Min consecutive samples that must exceed low_threshold", -) -parser.add_argument( - "--n-triggers", type=int, default=10, help="Expected number of APV triggers per RoR" -) -arg = parser.parse_args() - -from LDMX.Framework import ldmxcfg -from LDMX.Packing import rawio -from LDMX.Tracking import rawdecoder - -p = ldmxcfg.Process("trackerReco") -p.log_frequency = 1 -p.max_events = arg.max_events -p.output_files = [arg.output] - -# Stage 1: unpack raw Rogue frames. -unpacker = rawio.SingleSubsystemUnpacker( - dat_file=arg.dat, - output_name="TrackerRawData", - subsystem_name="tracker", - frame_offset=arg.frame_offset, -) - -# Stage 2: decode into RawSiStripHit. -decoder = rawdecoder.RawTrackerDecoder() -decoder.output_collection = "RawSiStripHits" - -# Pedestals are delivered as a conditions object (a service): the provider reads -# the JSON written by compute_pedestals.py and serves it to the processors below. -peds = rawdecoder.TrackerPedestalProvider(pedestal_file=arg.pedestal_file) - -# Stage 3: subtract pedestals -> subtracted RawSiStripHit collection. -subtractor = rawdecoder.PedestalSubtractor() -subtractor.input_collection = decoder.output_collection -subtractor.output_collection = "TrackerHits" - -# Stage 4: assemble per-channel waveforms -> SiStripWaveform collection. -builder = rawdecoder.SiStripWaveformBuilder() -builder.input_collection = subtractor.output_collection -builder.output_collection = "TrackerWaveforms" -builder.high_threshold = arg.high_threshold -builder.min_high_samples = arg.min_high_samples -builder.low_threshold = arg.low_threshold -builder.min_consecutive_low = arg.min_consecutive_low -builder.n_triggers = arg.n_triggers - -if arg.verbose_waveforms: - p.logger.trace(builder) - -p.sequence = [unpacker, decoder, subtractor, builder] diff --git a/Tracking/exampleConfigs/plot_measurements.py b/Tracking/exampleConfigs/plot_measurements.py new file mode 100644 index 000000000..26cc2e56f --- /dev/null +++ b/Tracking/exampleConfigs/plot_measurements.py @@ -0,0 +1,188 @@ +"""Basic diagnostic plots for a StripMeasurements collection. + +Reads the ROOT file written by raw_to_measurements.py and produces a small +set of PNGs: per-layer beam-spot occupancy, local-U distributions, cluster +amplitude and cluster size, and the global-position spread that shows the two +recoil stations. + +Coordinate note: the Measurement's global position is stored in the Acts frame +ordering, so component 0 is the beam axis (the two stations sit at ~ -243/-237 +and -143/-137 mm) and components 1,2 are the in-plane coordinates. + +Usage +----- + denv_workspace="$PWD" denv fire Tracking/exampleConfigs/plot_measurements.py \ + -- [--input meas.root] [--outdir plots_meas] + +--input accepts a single file or a glob over per-chunk outputs, e.g. + --input '/path/to/prod/meas_chunk_*.root' +so the batch chunks can be plotted directly without an hadd merge. +""" + +import argparse +import os +import sys + +sys.argv = [a for a in sys.argv if a != "--"] + +parser = argparse.ArgumentParser(f"ldmx fire {sys.argv[0]}") +parser.add_argument( + "--input", + nargs="+", + default=["meas.root"], + help="one file, or many (a shell glob over per-chunk outputs)", +) +parser.add_argument("--outdir", default="plots_meas") +parser.add_argument("--tree", default="LDMX_Events") +parser.add_argument("--branch", default="StripMeasurements_trackerReco") +arg = parser.parse_args() + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import ROOT + +ROOT.gSystem.Load("libTracking_Event") + +os.makedirs(arg.outdir, exist_ok=True) + +# --input is a list of files (a single file, or a shell glob over per-chunk +# outputs, e.g. '.../meas_chunk_*.root'). TChain reads them as one dataset, so +# there is no need to hadd the per-chunk outputs together first. +t = ROOT.TChain(arg.tree) +for pat in arg.input: + t.Add(pat) +input_label = arg.input[0] if len(arg.input) == 1 else f"{len(arg.input)} files" +print(f"chained {len(arg.input)} input path(s)") +meas = ROOT.std.vector("ldmx::Measurement")() +t.SetBranchAddress(arg.branch, meas) + +# Collect per-layer arrays. +layers = {} # layer_id -> dict of lists +n_per_event = [] +for i in range(t.GetEntries()): + t.GetEntry(i) + n_per_event.append(meas.size()) + for m in meas: + lid = m.getLayerID() + d = layers.setdefault( + lid, + {"gx": [], "gy": [], "gz": [], "u": [], "amp": [], "nstrips": [], "t": []}, + ) + g = m.getGlobalPosition() + lp = m.getLocalPosition() + d["gx"].append(g[0]) + d["gy"].append(g[1]) + d["gz"].append(g[2]) + d["u"].append(lp[0]) + d["amp"].append(m.getClusterAmplitude()) + d["nstrips"].append(m.getNStrips()) + d["t"].append(m.getTime()) + +layer_ids = sorted(layers) +total = sum(len(layers[lyr]["u"]) for lyr in layer_ids) +print( + f"{input_label}: {t.GetEntries()} events, {total} measurements, layers {layer_ids}" +) +for lid in layer_ids: + d = layers[lid] + print(f" layer {lid}: {len(d['u']):5d} meas, beam-pos {np.mean(d['gx']):.2f} mm") + +colors = dict( + zip(layer_ids, plt.cm.viridis(np.linspace(0.1, 0.85, len(layer_ids))), strict=True) +) + + +def save(fig, name): + path = os.path.join(arg.outdir, name) + fig.tight_layout() + fig.savefig(path, dpi=110) + plt.close(fig) + print("wrote", path) + + +# 1. Global in-plane occupancy (beam spot) per layer. +fig, axes = plt.subplots( + 1, len(layer_ids), figsize=(4 * len(layer_ids), 4), squeeze=False +) +for ax, lid in zip(axes[0], layer_ids, strict=True): + d = layers[lid] + ax.scatter(d["gy"], d["gz"], s=6, alpha=0.4, color=colors[lid]) + ax.set_title(f"layer {lid} (beam {np.mean(d['gx']):.1f} mm)") + ax.set_xlabel("global y [mm]") + ax.set_ylabel("global z [mm]") + ax.set_aspect("equal", "datalim") +save(fig, "occupancy_global.png") + +# 2. Local-U distributions. +fig, ax = plt.subplots(figsize=(7, 4.5)) +for lid in layer_ids: + ax.hist( + layers[lid]["u"], + bins=60, + histtype="step", + label=f"layer {lid}", + color=colors[lid], + ) +ax.set_xlabel("local U [mm]") +ax.set_ylabel("measurements") +ax.set_title("Cluster local-U position per layer") +ax.legend() +save(fig, "local_u.png") + +# 3. Cluster amplitude. +fig, ax = plt.subplots(figsize=(7, 4.5)) +amp_max = max(max(layers[lyr]["amp"]) for lyr in layer_ids) +bins = np.linspace(0, amp_max, 60) +for lid in layer_ids: + ax.hist( + layers[lid]["amp"], + bins=bins, + histtype="step", + label=f"layer {lid}", + color=colors[lid], + ) +ax.set_xlabel("cluster amplitude [ADC]") +ax.set_ylabel("measurements") +ax.set_title("Cluster amplitude per layer") +ax.legend() +save(fig, "cluster_amplitude.png") + +# 4. Cluster size. +fig, ax = plt.subplots(figsize=(7, 4.5)) +smax = max(max(layers[lyr]["nstrips"]) for lyr in layer_ids) +bins = np.arange(0.5, smax + 1.5, 1) +for lid in layer_ids: + ax.hist( + layers[lid]["nstrips"], + bins=bins, + histtype="step", + label=f"layer {lid}", + color=colors[lid], + ) +ax.set_xlabel("strips per cluster") +ax.set_ylabel("measurements") +ax.set_title("Cluster size per layer") +ax.legend() +save(fig, "cluster_size.png") + +# 5. Global beam-axis positions (shows the two stations + axial/stereo split). +fig, ax = plt.subplots(figsize=(7, 4.5)) +all_gx = [gx for lid in layer_ids for gx in layers[lid]["gx"]] +ax.hist(all_gx, bins=120, color="0.3") +ax.set_xlabel("global beam-axis position [mm]") +ax.set_ylabel("measurements") +ax.set_title("Measurement beam-axis position (two recoil stations)") +save(fig, "beam_axis_positions.png") + +# 6. Measurements per event. +fig, ax = plt.subplots(figsize=(7, 4.5)) +ax.hist(n_per_event, bins=np.arange(-0.5, max(n_per_event) + 1.5, 1), color="0.3") +ax.set_xlabel("measurements per event") +ax.set_ylabel("events") +ax.set_title(f"Measurements per event (mean {np.mean(n_per_event):.2f})") +save(fig, "meas_per_event.png") + +print("done") diff --git a/Tracking/exampleConfigs/plot_waveforms.py b/Tracking/exampleConfigs/plot_waveforms.py index 591a0ce62..9219d3991 100644 --- a/Tracking/exampleConfigs/plot_waveforms.py +++ b/Tracking/exampleConfigs/plot_waveforms.py @@ -1,10 +1,10 @@ -"""Plot SiStripWaveforms straight from a decode_to_waveforms.py ROOT file. +"""Plot SiStripWaveforms straight from a raw_to_measurements.py ROOT file. This is step 3 of the real-data tracker waveform chain: - 1. compute_pedestals.py -> pedestals.json - 2. decode_to_waveforms.py -> .root (TrackerWaveforms) - 3. plot_waveforms.py (this script) -> PNG figures + 1. compute_pedestals.py -> pedestals.json + 2. raw_to_measurements.py --stop-at fit -> .root (TrackerWaveforms) + 3. plot_waveforms.py (this script) -> PNG figures Unlike the in-framework approach, this is a STANDALONE script: it reads the ldmx::SiStripWaveform collection directly out of the ROOT event tree with @@ -50,7 +50,11 @@ parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter ) -parser.add_argument("root_file", help="ROOT file from decode_to_waveforms.py") +parser.add_argument( + "root_file", + help="ROOT file from raw_to_measurements.py (--stop-at fit " + "or measurements, so it carries the fit collection)", +) parser.add_argument( "--tree", default="LDMX_Events", help="Event tree name (default: LDMX_Events)" ) @@ -59,6 +63,21 @@ default="TrackerWaveforms", help="SiStripWaveform collection name (default: TrackerWaveforms)", ) +parser.add_argument( + "--fitted-collection", + default="FittedSiStripHits", + help="FittedSiStripHit collection holding the pulse fits " + "(default: FittedSiStripHits). The fit is no longer " + "stored on the waveform; it is joined back on via the " + "DAQ map. If the collection or the map is missing, the " + "fit overlay is simply omitted.", +) +parser.add_argument( + "--daq-map", + default=None, + help="DAQ map JSON used to join fits back onto waveforms " + "(default: search the install tree and ./Tracking/data)", +) parser.add_argument( "--n-examples", type=int, @@ -171,6 +190,18 @@ print(f"Reading '{branch_name}' from '{args.tree}' in {args.root_file}") + +def find_branch(collection): + """First branch named '' or '_', else None.""" + for b in tree.GetListOfBranches(): + name = b.GetName() + if name == collection or name.startswith(collection + "_"): + return name + return None + + +fitted_branch = find_branch(args.fitted_collection) + # --------------------------------------------------------------------------- # Scan the tree: collect top-N waveforms by peak significance + occupancy # --------------------------------------------------------------------------- @@ -233,10 +264,97 @@ def load_pedestal_noise(path): have_noise = bool(ped_noise) +# --------------------------------------------------------------------------- +# DAQ map: needed to join the fit results back onto the waveforms +# --------------------------------------------------------------------------- +# The pulse fit is not stored on the SiStripWaveform; SiStripWaveformFitProcessor +# writes it to a FittedSiStripHit keyed by (layer_id, strip_id). Applying the +# same DAQ-map strip transform the C++ uses (channelmap::stripId) turns a +# waveform's (feb, hybrid, pchannel) into that key, so the join is exact. +DAQ_MAP_CANDIDATES = [ + os.path.join( + os.environ.get("LDMX_INSTALL_PREFIX", ""), + "data/Tracking/daqmap_esa25_slice_test.json", + ), + "Tracking/data/daqmap_esa25_slice_test.json", + os.path.join( + os.path.dirname(os.path.abspath(__file__)), + "../data/daqmap_esa25_slice_test.json", + ), +] + + +def load_daq_map(path): + """Map (feb, hybrid) -> sensor record from the DAQ map JSON.""" + with open(path) as fh: + doc = json.load(fh) + sensors = {} + for s in doc["sensors"]: + sensors[(int(s["feb"]), int(s["hybrid"]))] = { + "layer_id": int(s["layer_id"]), + "n_strips": int(s["n_strips"]), + "first_strip": int(s["first_strip"]), + "reversed": bool(s["reversed"]), + } + return sensors + + +def strip_id(pch, sensor): + """Sensor strip index for a physical channel, matching channelmap::stripId.""" + if sensor["reversed"]: + return sensor["first_strip"] + sensor["n_strips"] - 1 - pch + return sensor["first_strip"] + pch + + +daq_map = {} +if fitted_branch is not None: + candidates = [args.daq_map] if args.daq_map else DAQ_MAP_CANDIDATES + for cand in candidates: + if cand and os.path.exists(cand): + daq_map = load_daq_map(cand) + print(f"Loaded DAQ map for {len(daq_map)} sensors from '{cand}'") + break + if not daq_map: + print( + "WARNING: no DAQ map found; cannot join fit results onto " + "waveforms, so the fit overlay will be omitted. Pass --daq-map." + ) + +have_fits = bool(fitted_branch) and bool(daq_map) +if fitted_branch is None: + print( + f"WARNING: no '{args.fitted_collection}' branch; fit overlay omitted. " + "Re-run raw_to_measurements.py with --stop-at fit (or measurements) " + "to produce it." + ) +elif have_fits: + print(f"Joining fits from '{fitted_branch}' via the DAQ map") + +NO_FIT = { + "fit_converged": False, + "fit_amplitude": 0.0, + "fit_t0": 0.0, + "fit_chi2": 0.0, + "fit_ndf": 0, +} + + top = [] # list of waveform record dicts occupancy = {} # (feb, hybrid) -> np.array(N_STRIPS) of counts for ievt, entry in enumerate(tree): + # Index this event's fits by the address the fit processor assigned them. + fits = {} + if have_fits: + for h in getattr(entry, fitted_branch): + fits[(int(h.getLayerID()), int(h.getStripID()))] = { + "fit_converged": True, + "fit_amplitude": float(h.getAmplitude()), + "fit_t0": float(h.getT0()), + "fit_chi2": float(h.getChi2()), + "fit_ndf": int(h.getNDF()), + } + collection = getattr(entry, branch_name) for wf in collection: samples = list(wf.getSamples()) @@ -253,25 +371,28 @@ def load_pedestal_noise(path): peak_amp = max(samples) if samples else 0 peak_sigma = (peak_amp / noise) if noise > 0 else 0.0 - top.append( - { - "event": ievt, - "feb": feb, - "hybrid": hybrid, - "pchannel": pch, - "n_triggers": int(wf.getNTriggers()), - "noise": noise, - "peak_amp": peak_amp, - "peak_sigma": peak_sigma, - "samples": samples, - "label": hybrid_label(feb, hybrid), - "fit_converged": bool(wf.isFitConverged()), - "fit_amplitude": float(wf.getFitAmplitude()), - "fit_t0": float(wf.getFitT0()), - "fit_chi2": float(wf.getFitChi2()), - "fit_ndf": int(wf.getFitNDF()), - } - ) + # Join the fit back on. A waveform with no entry either failed the fit + # or came from a hybrid the DAQ map does not cover; either way it is + # plotted without an overlay. + fit = NO_FIT + sensor = daq_map.get(key) + if sensor is not None: + fit = fits.get((sensor["layer_id"], strip_id(pch, sensor)), NO_FIT) + + record = { + "event": ievt, + "feb": feb, + "hybrid": hybrid, + "pchannel": pch, + "n_triggers": int(wf.getNTriggers()), + "noise": noise, + "peak_amp": peak_amp, + "peak_sigma": peak_sigma, + "samples": samples, + "label": hybrid_label(feb, hybrid), + } + record.update(fit) + top.append(record) tfile.Close() diff --git a/Tracking/exampleConfigs/raw_to_measurements.py b/Tracking/exampleConfigs/raw_to_measurements.py new file mode 100644 index 000000000..a60769821 --- /dev/null +++ b/Tracking/exampleConfigs/raw_to_measurements.py @@ -0,0 +1,174 @@ +"""Real-data tracker chain: raw bytes -> global-position Measurements. + +This is the step that connects the electronics-addressed waveform chain to the +geometry-aware reconstruction, using a DAQ map to turn (feb, hybrid, pchannel) +into (layer_id, strip_id). + +The chain can be truncated with --stop-at, which subsumes the old +decode_to_waveforms.py workflows: + + --stop-at waveforms -> TrackerWaveforms (no fit) + --stop-at fit -> TrackerWaveforms + FittedSiStripHits + --stop-at measurements -> ... + StripMeasurements (default) + +The 'fit' output is exactly what plot_waveforms.py reads -- same collection +names, same 'trackerReco' pass -- so waveform QA runs straight off this script +(point plot_waveforms.py at the output ROOT file). + +Pipeline +-------- + SingleSubsystemUnpacker -> TrackerRawData (bytes) + RawTrackerDecoder -> RawSiStripHits + PedestalSubtractor -> TrackerHits (pedestal-subtracted) + SiStripWaveformBuilder -> TrackerWaveforms (SiStripWaveform) + SiStripWaveformFitProcessor -> FittedSiStripHits (geometry-addressed) + StripClusterProcessor -> StripMeasurements (global x/y/z) + +Pedestals/noise are supplied as a conditions object (TrackerPedestals) by +TrackerPedestalProvider; the DAQ map is a plain JSON file. Run +compute_pedestals.py first to produce the pedestal JSON. + +Usage +----- + denv_workspace="$PWD" denv fire Tracking/exampleConfigs/raw_to_measurements.py -- \ + [--dat /path/to/physics.dat] \ + [--pedestal-file pedestals.json] \ + [--daq-map /path/to/daqmap.json] \ + [--detector ldmx-reduced-v3] \ + [--max-events N] \ + [--stop-at waveforms|fit|measurements] \ + [--output measurements.root] +""" + +import argparse +import sys + +# fire passes '--' through to the script; strip it so argparse sees only flags. +sys.argv = [a for a in sys.argv if a != "--"] + +parser = argparse.ArgumentParser(f"ldmx fire {sys.argv[0]}") +parser.add_argument( + "--dat", + default="/sdf/data/hps/users/mgignac/hardware/data/LDMX/Run_037_20251210_145218.dat", + help="Path to the physics Rogue .dat file", +) +parser.add_argument( + "--pedestal-file", + default="pedestals.json", + help="Pedestal JSON from compute_pedestals.py", +) +parser.add_argument( + "--daq-map", + default=None, + help="DAQ map JSON (default: the installed ESA slice-test map)", +) +parser.add_argument( + "--detector", + default="ldmx-reduced-v3", + help="Detector name for the tracking geometry (default: ldmx-reduced-v3)", +) +parser.add_argument("--max-events", type=int, default=300) +parser.add_argument("--frame-offset", type=int, default=0) +parser.add_argument("--output", default="measurements.root") +parser.add_argument("--n-triggers", type=int, default=10) +# Waveform-builder significance cuts (both conditions must pass). +parser.add_argument("--high-threshold", type=float, default=5.0) +parser.add_argument("--min-high-samples", type=int, default=4) +parser.add_argument("--low-threshold", type=float, default=3.0) +parser.add_argument("--min-consecutive-low", type=int, default=5) +# Clustering significance cuts. +parser.add_argument("--seed-threshold", type=float, default=4.0) +parser.add_argument("--neighbor-threshold", type=float, default=3.0) +parser.add_argument("--cluster-threshold", type=float, default=4.0) +# Chain truncation (subsumes the old decode_to_waveforms.py --no-fit / QA modes). +parser.add_argument( + "--stop-at", + choices=("waveforms", "fit", "measurements"), + default="measurements", + help="Truncate the chain: 'waveforms' (only TrackerWaveforms, no fit / no " + "DAQ map), 'fit' (adds FittedSiStripHits -- what plot_waveforms.py " + "wants), or 'measurements' (full chain to StripMeasurements, default).", +) +parser.add_argument( + "--verbose-waveforms", + action="store_true", + help="Set the waveform builder (and fitter, if run) to trace logging, " + "printing per-channel fit results and full ASCII waveform traces.", +) +arg = parser.parse_args() + +from LDMX.Framework import ldmxcfg +from LDMX.Packing import rawio +from LDMX.Tracking import rawdecoder, tracking +from LDMX.Tracking.geo import TrackersTrackingGeometryProvider as TrackGeo + +daq_map = arg.daq_map if arg.daq_map is not None else rawdecoder.daq_map_path() + +p = ldmxcfg.Process("trackerReco") +p.log_frequency = 1 +p.max_events = arg.max_events +p.output_files = [arg.output] + +# Tracking geometry for the local -> global transform in StripClusterProcessor. +TrackGeo.get_instance().set_detector(arg.detector) + +# Stage 1: unpack raw Rogue frames. +unpacker = rawio.SingleSubsystemUnpacker( + dat_file=arg.dat, + output_name="TrackerRawData", + subsystem_name="tracker", + frame_offset=arg.frame_offset, +) + +# Stage 2: decode into RawSiStripHit. +decoder = rawdecoder.RawTrackerDecoder() +decoder.output_collection = "RawSiStripHits" + +# Pedestals delivered as a conditions object (service). +peds = rawdecoder.TrackerPedestalProvider(pedestal_file=arg.pedestal_file) + +# Stage 3: subtract pedestals. +subtractor = rawdecoder.PedestalSubtractor() +subtractor.input_collection = decoder.output_collection +subtractor.output_collection = "TrackerHits" + +# Stage 4: assemble per-channel waveforms. +builder = rawdecoder.SiStripWaveformBuilder() +builder.input_collection = subtractor.output_collection +builder.output_collection = "TrackerWaveforms" +builder.n_triggers = arg.n_triggers +builder.high_threshold = arg.high_threshold +builder.min_high_samples = arg.min_high_samples +builder.low_threshold = arg.low_threshold +builder.min_consecutive_low = arg.min_consecutive_low + +# Stage 5: fit the pulse shape and map electronics addresses to +# (layer_id, strip_id) via the DAQ map. +fitter = rawdecoder.SiStripWaveformFitProcessor() +fitter.input_collection = builder.output_collection +fitter.output_collection = "FittedSiStripHits" +fitter.daq_map_file = daq_map + +# Stage 6: cluster into global-position Measurements. +clusterer = tracking.StripClusterProcessor() +clusterer.in_collection = fitter.output_collection +clusterer.out_collection = "StripMeasurements" +clusterer.daq_map_file = daq_map +clusterer.seed_threshold = arg.seed_threshold +clusterer.neighbor_threshold = arg.neighbor_threshold +clusterer.cluster_threshold = arg.cluster_threshold + +# Assemble the sequence, truncating at --stop-at. The builder always runs; the +# fitter is added for 'fit'/'measurements', the clusterer only for the full +# 'measurements' chain. (In 'waveforms' mode the fitter/clusterer objects above +# are simply left out of the sequence, and the DAQ map is never opened.) +p.sequence = [unpacker, decoder, subtractor, builder] +if arg.stop_at in ("fit", "measurements"): + p.sequence.append(fitter) +if arg.stop_at == "measurements": + p.sequence.append(clusterer) + +if arg.verbose_waveforms: + p.logger.trace(builder) + if arg.stop_at in ("fit", "measurements"): + p.logger.trace(fitter) diff --git a/Tracking/include/Tracking/Digitization/StripClusterer.h b/Tracking/include/Tracking/Digitization/StripClusterer.h index 352fdd887..2d1a7839d 100644 --- a/Tracking/include/Tracking/Digitization/StripClusterer.h +++ b/Tracking/include/Tracking/Digitization/StripClusterer.h @@ -99,6 +99,15 @@ class StripClusterer { double cluster_weighted_t, double cluster_total_amp) const; + /// Per-strip noise RMS to use for @p h: the hit's own measured noise if set + /// (> 0), otherwise the uniform noise passed to the constructor. This is how + /// real-data hits (with per-channel pedestal noise) and MC hits (uniform) + /// share one clustering path. + double hitNoise(const ldmx::FittedSiStripHit& h) const { + const double n = h.getNoise(); + return n > 0.0 ? n : noise_sigma_adc_; + } + double seed_threshold_; double neighbor_threshold_; double cluster_threshold_; diff --git a/Tracking/include/Tracking/Event/FittedSiStripHit.h b/Tracking/include/Tracking/Event/FittedSiStripHit.h index 59e85c942..ebccfe7e4 100644 --- a/Tracking/include/Tracking/Event/FittedSiStripHit.h +++ b/Tracking/include/Tracking/Event/FittedSiStripHit.h @@ -26,7 +26,7 @@ class FittedSiStripHit { FittedSiStripHit(int layer_id, int strip_id, float amplitude, float t0, float chi2, int ndf, int track_id = -1, int pdg_id = 0, - int sim_hit_id = -1, float edep = 0.f) + int sim_hit_id = -1, float edep = 0.f, float noise = 0.f) : layer_id_(layer_id), strip_id_(strip_id), amplitude_(amplitude), @@ -36,7 +36,8 @@ class FittedSiStripHit { track_id_(track_id), pdg_id_(pdg_id), sim_hit_id_(sim_hit_id), - edep_(edep) {} + edep_(edep), + noise_(noise) {} virtual ~FittedSiStripHit() = default; @@ -51,6 +52,7 @@ class FittedSiStripHit { pdg_id_ = 0; sim_hit_id_ = -1; edep_ = 0.f; + noise_ = 0.f; } // --- Getters --- @@ -71,6 +73,9 @@ class FittedSiStripHit { int getSimHitID() const { return sim_hit_id_; } /// Energy deposited by the parent SimTrackerHit [MeV] (0 if unknown). float getEdep() const { return edep_; } + /// Per-strip noise RMS [ADC counts] (0 if unknown; clustering then falls back + /// to the global StripClusterer noise). + float getNoise() const { return noise_; } // --- Setters --- void setLayerID(int v) { layer_id_ = v; } @@ -83,13 +88,15 @@ class FittedSiStripHit { void setPdgID(int v) { pdg_id_ = v; } void setSimHitID(int v) { sim_hit_id_ = v; } void setEdep(float v) { edep_ = v; } + void setNoise(float v) { noise_ = v; } friend std::ostream& operator<<(std::ostream& o, const FittedSiStripHit& h) { o << "[ FittedSiStripHit ]: layer=" << h.layer_id_ << " strip=" << h.strip_id_ << " amp=" << h.amplitude_ << " ADC" << " t0=" << h.t0_ << " ns" << " chi2/ndf=" << h.chi2_ << "/" << h.ndf_ << " track_id=" << h.track_id_ << " pdg_id=" << h.pdg_id_ - << " sim_hit_id=" << h.sim_hit_id_ << " edep=" << h.edep_ << " MeV"; + << " sim_hit_id=" << h.sim_hit_id_ << " edep=" << h.edep_ << " MeV" + << " noise=" << h.noise_ << " ADC"; return o; } @@ -111,8 +118,10 @@ class FittedSiStripHit { int sim_hit_id_{-1}; /// Energy deposited by the parent SimTrackerHit [MeV]. float edep_{0.f}; + /// Per-strip noise RMS [ADC counts]; 0 means unset. + float noise_{0.f}; - ClassDef(FittedSiStripHit, 1); + ClassDef(FittedSiStripHit, 2); }; } // namespace ldmx diff --git a/Tracking/include/Tracking/Event/SiStripWaveform.h b/Tracking/include/Tracking/Event/SiStripWaveform.h index 68bd27291..719e6d7fc 100644 --- a/Tracking/include/Tracking/Event/SiStripWaveform.h +++ b/Tracking/include/Tracking/Event/SiStripWaveform.h @@ -38,24 +38,6 @@ class SiStripWaveform { uint8_t getFebId() const { return feb_id_; } uint8_t getNTriggers() const { return n_triggers_; } - /// Store the result of a pulse-shape fit to this waveform. - void setFitResult(float amplitude, float t0, float chi2, int ndf, - bool converged) { - fit_amplitude_ = amplitude; - fit_t0_ = t0; - fit_chi2_ = chi2; - fit_ndf_ = ndf; - fit_converged_ = converged; - } - - /// Fitted pulse amplitude [ADC counts] (peak of the fitted pulse shape). - float getFitAmplitude() const { return fit_amplitude_; } - /// Fitted hit arrival time T [ns] in the sample-window frame (t_i = i*dt). - float getFitT0() const { return fit_t0_; } - float getFitChi2() const { return fit_chi2_; } - int getFitNDF() const { return fit_ndf_; } - bool isFitConverged() const { return fit_converged_; } - /// Sample at trigger index t (0-based), APV sample s (0-2). short getSample(uint8_t t, uint8_t s) const { return samples_[t * 3 + s]; } @@ -84,14 +66,7 @@ class SiStripWaveform { uint8_t feb_id_{0}; uint8_t n_triggers_{0}; ///< number of APV triggers assembled - // Pulse-shape fit result (filled by SiStripWaveformBuilder). - float fit_amplitude_{0}; ///< fitted amplitude [ADC] - float fit_t0_{0}; ///< fitted hit arrival time T [ns] - float fit_chi2_{0}; ///< chi-squared at the minimum - int fit_ndf_{0}; ///< degrees of freedom = n_samples - 2 - bool fit_converged_{false}; ///< true if the fit succeeded - - ClassDef(SiStripWaveform, 2); + ClassDef(SiStripWaveform, 3); }; } // namespace ldmx diff --git a/Tracking/include/Tracking/Reco/SiStripChannelMap.h b/Tracking/include/Tracking/Reco/SiStripChannelMap.h index fb1c03e74..75bfa9ba4 100644 --- a/Tracking/include/Tracking/Reco/SiStripChannelMap.h +++ b/Tracking/include/Tracking/Reco/SiStripChannelMap.h @@ -41,6 +41,38 @@ inline constexpr int16_t pchannel(uint8_t apv_id, uint8_t channel) { (apv_id * K_CHANNELS_PER_APV + (K_CHANNELS_PER_APV - 1) - channel)); } +/** + * Inverse of pchannel(): recover the (APV id, APV channel) pair from a physical + * strip number. Needed because SiStripWaveform stores only the pchannel, but + * the pedestal table is keyed by the full electronics address. + * + * Round-trips exactly with pchannel() for all valid (apv, channel): + * apvChannelFromPchannel(pchannel(a, c), a', c') gives a' == a, c' == c. + */ +inline constexpr void apvChannelFromPchannel(int16_t pchannel, uint8_t& apv_id, + uint8_t& channel) { + const int x = (K_CHANNELS_PER_HYBRID - 1) - pchannel; + apv_id = static_cast(x / K_CHANNELS_PER_APV); + channel = + static_cast((K_CHANNELS_PER_APV - 1) - (x % K_CHANNELS_PER_APV)); +} + +/** + * Map a physical strip number (pchannel) onto a sensor strip index using a + * DAQ-map transform: an optional readout reversal plus a per-sensor offset. + * + * @param pchannel physical strip within the hybrid, [0, kChannelsPerHybrid). + * @param n_strips number of bonded strips on the sensor. + * @param first_strip sensor strip index that pchannel 0 (or the reversed end) + * maps to. + * @param reversed true if the hybrid reads the sensor in descending order. + */ +inline constexpr int stripId(int16_t pchannel, int n_strips, int first_strip, + bool reversed) { + return reversed ? first_strip + n_strips - 1 - pchannel + : first_strip + pchannel; +} + /// Build the per-channel pedestal-map key "feb:hybrid:apv:channel". inline std::string channelKey(uint8_t feb, uint8_t hybrid, uint8_t apv, uint8_t channel) { diff --git a/Tracking/include/Tracking/Reco/SiStripWaveformBuilder.h b/Tracking/include/Tracking/Reco/SiStripWaveformBuilder.h index af7713947..57250f1f2 100644 --- a/Tracking/include/Tracking/Reco/SiStripWaveformBuilder.h +++ b/Tracking/include/Tracking/Reco/SiStripWaveformBuilder.h @@ -1,11 +1,9 @@ #ifndef TRACKING_RECO_SISTRIPWAVEFORMBUILDER_H_ #define TRACKING_RECO_SISTRIPWAVEFORMBUILDER_H_ -#include #include #include "Framework/EventProcessor.h" -#include "Tracking/Digitization/PulseShape.h" #include "Tracking/Event/SiStripWaveform.h" namespace tracking::reco { @@ -13,25 +11,6 @@ namespace tracking::reco { /** * Assemble per-trigger pedestal-subtracted hits into full per-channel * waveforms. - * - * Each RoR issues n_triggers APV triggers, each producing a 3-sample - * (pedestal-subtracted) RawSiStripHit per channel. This producer groups all - * trigger hits for the same (feb, hybrid, pchannel), sorts by apv_trigger, and - * concatenates their samples into a SiStripWaveform with n_triggers*3 samples - * ordered as [s0_t0, s1_t0, s2_t0, s0_t1, ...]. - * - * The per-channel noise used for the significance cuts is obtained from the - * TrackerPedestals conditions object (a service); the physical strip number - * (pchannel) is derived from the hit's apv_id/channel via channelmap. - * - * Two-stage threshold applied before writing: - * 1. High-threshold count: at least min_high_samples samples with - * ADC/noise > high_threshold (default: 4 samples > 5σ). - * 2. Consecutive low-threshold streak: at least min_consecutive_low - * consecutive samples with ADC/noise > low_threshold - * (default: 5 consecutive samples > 3σ). - * Both conditions must be satisfied; this strongly suppresses noise while - * retaining real APV25 signal pulses. */ class SiStripWaveformBuilder : public framework::Producer { public: @@ -40,7 +19,6 @@ class SiStripWaveformBuilder : public framework::Producer { void configure(framework::config::Parameters& ps) override; void produce(framework::Event& event) override; - void onProcessEnd() override; private: std::string input_collection_{"TrackerHits"}; @@ -52,16 +30,8 @@ class SiStripWaveformBuilder : public framework::Producer { double low_threshold_{ 3.0}; ///< per-sample significance for consecutive-streak cut int min_consecutive_low_{ - 5}; ///< min consecutive samples exceeding low_threshold + 5}; ///< min consecutive samples exceeding low_threshold int n_triggers_{10}; ///< expected APV triggers per RoR - - /// Pulse shape used for the per-waveform fit test (built lazily in produce). - std::unique_ptr pulse_shape_; - - // Fit monitoring counters (summed over the whole job, reported in - // onProcessEnd). - long n_fit_attempted_{0}; ///< waveforms passed to the fitter - long n_fit_failed_{0}; ///< fits that did not converge }; } // namespace tracking::reco diff --git a/Tracking/include/Tracking/Reco/SiStripWaveformFitProcessor.h b/Tracking/include/Tracking/Reco/SiStripWaveformFitProcessor.h new file mode 100644 index 000000000..b1d7879ef --- /dev/null +++ b/Tracking/include/Tracking/Reco/SiStripWaveformFitProcessor.h @@ -0,0 +1,81 @@ +#ifndef TRACKING_RECO_SISTRIPWAVEFORMFITPROCESSOR_H_ +#define TRACKING_RECO_SISTRIPWAVEFORMFITPROCESSOR_H_ + +#include +#include +#include + +#include "Framework/EventProcessor.h" +#include "Tracking/Digitization/PulseShape.h" +#include "Tracking/Reco/TrackerDaqMap.h" + +namespace tracking::reco { + +/** + * Fit a pulse shape to each SiStripWaveform and produce FittedSiStripHits. + * This is the real-data counterpart of StripFitProcessor, and the bridge into + * the geometry-aware reconstruction. + * + * Input : collection of ldmx::SiStripWaveform (pedestal-subtracted samples) + * Output : collection of ldmx::FittedSiStripHit + * + * Configuration parameters + * ------------------------ + * input_collection SiStripWaveform input collection. + * input_pass_name Pass name for the input collection. + * output_collection FittedSiStripHit output collection. + * daq_map_file Path to the DAQ map JSON (required). + * t_scan_min_ns Lower bound of the hit-time scan [ns] (default -50). + * t_scan_max_ns Upper bound of the hit-time scan [ns]; <= 0 means auto, + * i.e. n_samples * sampling interval (default -1). + * t_scan_step_ns Step size of the coarse scan [ns] (default 1). + * max_chi2_ndf If > 0, discard fits with chi2/ndf above this value + * (default -1 = off). + */ +class SiStripWaveformFitProcessor : public framework::Producer { + public: + SiStripWaveformFitProcessor(const std::string& name, + framework::Process& process) + : framework::Producer(name, process) {} + + void configure(framework::config::Parameters& ps) override; + void onProcessStart() override; + void produce(framework::Event& event) override; + void onProcessEnd() override; + + private: + std::string input_collection_{"TrackerWaveforms"}; + std::string input_pass_name_{""}; + std::string output_collection_{"FittedSiStripHits"}; + std::string daq_map_file_{""}; + + // Scan range for the hit-time search. A non-positive t_scan_max_ns_ sizes + // the scan to each waveform, since T may peak anywhere from before sample 0 + // to the last sample and waveforms differ in length. + double t_scan_min_ns_{-50.0}; + double t_scan_max_ns_{-1.0}; + double t_scan_step_ns_{1.0}; + + // Quality cut (<= 0 means disabled) + double max_chi2_ndf_{-1.0}; + + TrackerDaqMap daq_map_; + + /// Pulse shape shared by every per-channel fitter (built in onProcessStart). + std::unique_ptr pulse_shape_; + + // Diagnostics accumulated over the job. + long n_waveforms_{0}; ///< waveforms seen + long n_unmapped_{0}; ///< skipped: (feb, hybrid) absent from the DAQ map + long n_out_of_range_{0}; ///< skipped: strip_id outside the sensor + long n_fit_attempted_{0}; ///< waveforms passed to the fitter + long n_unconverged_{0}; ///< skipped because the pulse fit did not converge + long n_bad_chi2_{0}; ///< skipped by the chi2/ndf cut + long n_hits_{0}; ///< emitted FittedSiStripHits + /// (feb, hybrid) pairs already warned about, so each is reported only once. + std::map unmapped_sensors_; +}; + +} // namespace tracking::reco + +#endif // TRACKING_RECO_SISTRIPWAVEFORMFITPROCESSOR_H_ diff --git a/Tracking/include/Tracking/Reco/StripClusterProcessor.h b/Tracking/include/Tracking/Reco/StripClusterProcessor.h index 872404a1f..0edfff901 100644 --- a/Tracking/include/Tracking/Reco/StripClusterProcessor.h +++ b/Tracking/include/Tracking/Reco/StripClusterProcessor.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -7,6 +8,7 @@ #include "Framework/Event.h" #include "Tracking/Digitization/SiStripConstants.h" #include "Tracking/Digitization/StripClusterer.h" +#include "Tracking/Reco/TrackerDaqMap.h" #include "Tracking/Reco/TrackingGeometryUser.h" namespace tracking::reco { @@ -59,6 +61,13 @@ class StripClusterProcessor : public TrackingGeometryUser { std::string in_pass_{""}; std::string out_collection_{"StripMeasurements"}; + // Optional DAQ map: when set, the local-U centre offset for a layer uses that + // sensor's real strip count instead of the fixed N_READOUT_STRIPS constant. + // Empty (the default, MC) keeps the constant and leaves the MC path + // unchanged. + std::string daq_map_file_{""}; + std::map layer_n_strips_; ///< layer_id -> n_strips (from DAQ map) + // Clustering parameters (forwarded to StripClusterer) double seed_threshold_{4.0}; double neighbor_threshold_{3.0}; diff --git a/Tracking/include/Tracking/Reco/TrackerDaqMap.h b/Tracking/include/Tracking/Reco/TrackerDaqMap.h new file mode 100644 index 000000000..554893fd1 --- /dev/null +++ b/Tracking/include/Tracking/Reco/TrackerDaqMap.h @@ -0,0 +1,75 @@ +#ifndef TRACKING_RECO_TRACKERDAQMAP_H_ +#define TRACKING_RECO_TRACKERDAQMAP_H_ + +#include +#include +#include + +namespace tracking::reco { + +/** + * Tracker DAQ map: the correspondence between an electronics sensor address + * (feb, hybrid) and the detector layer it reads, plus the strip transform that + * turns a physical channel (pchannel) into a sensor strip index. + */ +class TrackerDaqMap { + public: + /// Everything the mapper needs to know about one physical sensor. + struct SensorInfo { + int layer_id_{-1}; ///< Acts surface id (volume*1000 + layer*100 + sensor). + int n_strips_{0}; ///< Number of bonded strips on the sensor. + int first_strip_{0}; ///< Sensor strip index that pchannel 0 maps to. + bool reversed_{false}; ///< True if the hybrid reads the sensor descending. + }; + + TrackerDaqMap() = default; + + /** + * Load a DAQ map from a JSON file. + * + * Expected schema: + * @code + * { "sensors": [ { "feb": 0, "hybrid": 0, "layer_id": 3101, + * "n_strips": 640, "first_strip": 0, "reversed": false }, + * ... ] } + * @endcode + * Extra keys (e.g. "station", "orientation", "_comment") are ignored. + * + * @throws framework::exception::Exception if the file cannot be opened, is + * not valid JSON, lacks a non-empty "sensors" array, has a malformed entry, + * or defines the same (feb, hybrid) twice. There is no silent empty-map + * failure mode. + */ + static TrackerDaqMap fromJsonFile(const std::string& path); + + /// True if the map has an entry for this electronics sensor. + bool has(uint8_t feb, uint8_t hybrid) const { + return sensors_.find(key(feb, hybrid)) != sensors_.end(); + } + + /** + * Look up the sensor read by (feb, hybrid). + * @throws framework::exception::Exception if there is no such entry; callers + * that want to skip unmapped sensors must guard with has() first. + */ + const SensorInfo& at(uint8_t feb, uint8_t hybrid) const; + + /// Number of sensors in the map. + std::size_t size() const { return sensors_.size(); } + + /// All sensors, keyed by (feb << 8) | hybrid. For consumers that need to + /// index by something other than the electronics address (e.g. build a + /// layer_id -> n_strips lookup). + const std::map& sensors() const { return sensors_; } + + private: + static uint16_t key(uint8_t feb, uint8_t hybrid) { + return static_cast((static_cast(feb) << 8) | hybrid); + } + + std::map sensors_; ///< key = (feb << 8) | hybrid +}; + +} // namespace tracking::reco + +#endif // TRACKING_RECO_TRACKERDAQMAP_H_ diff --git a/Tracking/python/rawdecoder.py b/Tracking/python/rawdecoder.py index 08f376038..3acbcdf7f 100644 --- a/Tracking/python/rawdecoder.py +++ b/Tracking/python/rawdecoder.py @@ -6,6 +6,16 @@ ) +def daq_map_path(name="daqmap_esa25_slice_test.json"): + """Absolute path to a DAQ-map JSON installed under data/Tracking. + + The install prefix is substituted by cmake's configure_file when this module + is installed (the same mechanism LDMX.Detectors.make_path uses), so this + resolves to the real installed location at runtime. + """ + return "@CMAKE_INSTALL_PREFIX@/data/Tracking/" + name + + @processor("tracking::reco::RawTrackerDecoder", "Tracking") class RawTrackerDecoder(Processor): """Decode raw Rogue frame bytes into a RawSiStripHit collection. @@ -130,8 +140,9 @@ class SiStripWaveformBuilder(Processor): n_triggers : int Expected number of APV triggers per RoR (default: 10). - Per-fit results and full waveform traces are emitted at the 'trace' logging - level; set the processor's logging level to see them. + No pulse fitting happens here; that is SiStripWaveformFitProcessor's job. + Full waveform traces are emitted at the 'trace' logging level; set the + processor's logging level to see them. """ input_collection: str = "TrackerHits" @@ -142,3 +153,49 @@ class SiStripWaveformBuilder(Processor): low_threshold: float = 3.0 min_consecutive_low: int = 5 n_triggers: int = 10 + + +@processor("tracking::reco::SiStripWaveformFitProcessor", "Tracking") +class SiStripWaveformFitProcessor(Processor): + """Fit SiStripWaveforms and produce geometry-addressed FittedSiStripHits. + + The real-data counterpart of tracking.StripFitProcessor: it fits the same + pulse shape with the same fitter, differing only in that the sensor address + comes from a DAQ map (JSON) instead of the hit itself, and the per-sample + noise comes from the TrackerPedestals conditions object instead of a fixed + constant. Each waveform's (feb, hybrid, pchannel) becomes a + (layer_id, strip_id); the resulting FittedSiStripHit collection feeds the + standard, geometry-aware StripClusterProcessor exactly as the MC chain does. + + Attributes + ---------- + input_collection : str + SiStripWaveform collection to read (default: 'TrackerWaveforms'). + input_pass_name : str + Pass name of the upstream producer (empty = any pass). + output_collection : str + Name for the output FittedSiStripHit collection (default: + 'FittedSiStripHits'). + daq_map_file : str + Path to the DAQ map JSON file; required (loaded at onProcessStart). + t_scan_min_ns : float + Lower bound of the hit-time scan [ns] (default: -50). + t_scan_max_ns : float + Upper bound of the hit-time scan [ns]. Non-positive means auto, i.e. + sized to each waveform as n_samples * 25 ns (default: -1). + t_scan_step_ns : float + Step size of the coarse hit-time scan [ns] (default: 1). + max_chi2_ndf : float + If > 0, discard fits with chi2/ndf above this value (default: -1, off). + + Per-fit results are emitted at the 'trace' logging level. + """ + + input_collection: str = "TrackerWaveforms" + input_pass_name: str = "" + output_collection: str = "FittedSiStripHits" + daq_map_file: str = "" + t_scan_min_ns: float = -50.0 + t_scan_max_ns: float = -1.0 + t_scan_step_ns: float = 1.0 + max_chi2_ndf: float = -1.0 diff --git a/Tracking/python/tracking.py b/Tracking/python/tracking.py index 0afce8608..e708cc586 100644 --- a/Tracking/python/tracking.py +++ b/Tracking/python/tracking.py @@ -58,6 +58,9 @@ class DigitizationProcessor(Processor): Minimum number of charge deposition segments (Mode 1). out_raw_collection : str Output RawSiStripHit collection name; empty disables (Mode 1). + dump_geo_csv : str + If non-empty, write every ACTS surface (id, centre, U/V/W axes) to this + CSV at onProcessStart and continue; empty disables. """ merge_hits: bool = True @@ -81,6 +84,7 @@ class DigitizationProcessor(Processor): deposition_granularity: float = 0.10 n_segments_min: int = 5 out_raw_collection: str = "" + dump_geo_csv: str = "" @processor("tracking::reco::SeedFinderProcessor", "Tracking") @@ -593,6 +597,10 @@ class StripClusterProcessor(Processor): <= 0 disables (default -1). max_chi2_ndf : float Max chi2/ndf for a fitted hit to be used; <= 0 disables (default -1). + daq_map_file : str + Optional DAQ map JSON. When set, the local-U centre offset for each + layer uses that sensor's real strip count instead of the fixed + N_READOUT_STRIPS constant. Empty (default, MC) keeps the constant. """ in_collection: str = "FittedSiStripHits" @@ -605,3 +613,4 @@ class StripClusterProcessor(Processor): time_window_ns: float = -1.0 neighbor_delta_t_ns: float = -1.0 max_chi2_ndf: float = -1.0 + daq_map_file: str = "" diff --git a/Tracking/src/Tracking/Digitization/StripClusterer.cxx b/Tracking/src/Tracking/Digitization/StripClusterer.cxx index 34e241038..c29257082 100644 --- a/Tracking/src/Tracking/Digitization/StripClusterer.cxx +++ b/Tracking/src/Tracking/Digitization/StripClusterer.cxx @@ -71,19 +71,22 @@ std::vector StripClusterer::findClusters( // ------------------------------------------------------------------------- // Determine which strips are clusterable (≥ neighbor threshold) and which // can seed a cluster (≥ seed threshold + timing/chi2 cuts). + // + // Thresholds are in units of the per-strip noise RMS. Each hit may carry its + // own measured noise (getNoise() > 0, from the real-data pedestal table); + // when it does not (MC), we fall back to the uniform ctor noise so the MC + // path is unchanged. // ------------------------------------------------------------------------- - const double snr_neighbor = neighbor_threshold_ * noise_sigma_adc_; - const double snr_seed = seed_threshold_ * noise_sigma_adc_; - std::set clusterable_set; std::vector seed_channels; for (const auto& [ch, hp] : channel_map) { const double amp = hp->getAmplitude(); - if (amp >= snr_neighbor) { + const double noise = hitNoise(*hp); + if (amp >= neighbor_threshold_ * noise) { clusterable_set.insert(ch); } - if (amp >= snr_seed && passesSeedCuts(*hp)) { + if (amp >= seed_threshold_ * noise && passesSeedCuts(*hp)) { seed_channels.push_back(ch); } } @@ -120,10 +123,11 @@ std::vector StripClusterer::findClusters( const double amp = hit.getAmplitude(); // Accumulate cluster quantities. + const double noise = hitNoise(hit); cand.strip_ids.push_back(cur_ch); cluster_total_amp += amp; cluster_weighted_t += amp * hit.getT0(); - cluster_noise_sq += noise_sigma_adc_ * noise_sigma_adc_; + cluster_noise_sq += noise * noise; // Check nearest neighbours (strip ± 1). for (int delta : {-1, +1}) { diff --git a/Tracking/src/Tracking/Event/SiStripWaveform.cxx b/Tracking/src/Tracking/Event/SiStripWaveform.cxx index e84c4c5cf..74e6e67e4 100644 --- a/Tracking/src/Tracking/Event/SiStripWaveform.cxx +++ b/Tracking/src/Tracking/Event/SiStripWaveform.cxx @@ -21,11 +21,6 @@ void SiStripWaveform::clear() { hybrid_id_ = 0; feb_id_ = 0; n_triggers_ = 0; - fit_amplitude_ = 0; - fit_t0_ = 0; - fit_chi2_ = 0; - fit_ndf_ = 0; - fit_converged_ = false; } std::ostream& operator<<(std::ostream& output, const SiStripWaveform& w) { @@ -34,12 +29,7 @@ std::ostream& operator<<(std::ostream& output, const SiStripWaveform& w) { << " Hybrid=" << static_cast(w.hybrid_id_) << " PCh=" << w.pchannel_ << " NTrig=" << static_cast(w.n_triggers_) - << " PeakAmp=" << w.peakAmplitude(); - if (w.fit_converged_) { - output << " | fit amp=" << w.fit_amplitude_ << " t0=" << w.fit_t0_ << "ns" - << " chi2/ndf=" << w.fit_chi2_ << "/" << w.fit_ndf_; - } - output << "\n"; + << " PeakAmp=" << w.peakAmplitude() << "\n"; if (w.samples_.empty()) return output; diff --git a/Tracking/src/Tracking/Reco/SiStripWaveformBuilder.cxx b/Tracking/src/Tracking/Reco/SiStripWaveformBuilder.cxx index 0901dca06..487a2ffa2 100644 --- a/Tracking/src/Tracking/Reco/SiStripWaveformBuilder.cxx +++ b/Tracking/src/Tracking/Reco/SiStripWaveformBuilder.cxx @@ -3,8 +3,6 @@ #include #include -#include "Tracking/Digitization/SiStripConstants.h" -#include "Tracking/Digitization/StripPulseFitter.h" #include "Tracking/Event/RawSiStripHit.h" #include "Tracking/Reco/SiStripChannelMap.h" #include "Tracking/Reco/TrackerPedestals.h" @@ -32,14 +30,6 @@ void SiStripWaveformBuilder::produce(framework::Event& event) { const auto& hits = event.getCollection(input_collection_, input_pass_name_); - // Build the pulse shape once (reused for every channel fit this job). - if (!pulse_shape_) { - pulse_shape_ = tracking::digitization::PulseShape::make( - std::string(tracking::digitization::PULSE_SHAPE_NAME), - tracking::digitization::PEAKING_TIME_NS, - tracking::digitization::SECOND_TIME_CONST_NS); - } - // Key: encodes (feb, hybrid, pchannel) uniquely. // Value: vector of (apv_trigger, samples) pairs. struct TriggerSamples { @@ -114,38 +104,8 @@ void SiStripWaveformBuilder::produce(framework::Event& event) { uint8_t n_trig = static_cast( std::min(static_cast(ch.triggers_.size()), 255)); - // --- TEST: fit a CR-RC pulse shape to the full assembled waveform. --- - // Samples are pedestal-subtracted (ped = 0) and lie on a uniform 25 ns - // grid, so the scan range is sized to the waveform: T may peak anywhere - // from before sample 0 to the last sample. - const int n_samples = static_cast(samples.size()); - const double t_scan_max = - n_samples * tracking::digitization::SAMPLING_INTERVAL_NS; - tracking::digitization::StripPulseFitter fitter( - *pulse_shape_, - /*t0_offset_ns=*/0.0, - tracking::digitization::SAMPLING_INTERVAL_NS, - /*pedestal_adc=*/0.0, - /*noise_sigma_adc=*/ch.noise_, - /*t_scan_min_ns=*/-50.0, t_scan_max, /*t_scan_step_ns=*/1.0); - const auto fit = fitter.fit(samples); - ++n_fit_attempted_; - if (!fit.converged) ++n_fit_failed_; - - ldmx_log(trace) << "fit feb=" << static_cast(ch.feb_id_) - << " hyb=" << static_cast(ch.hybrid_id_) - << " pch=" << pchannel << " nsamp=" << n_samples - << " noise=" << ch.noise_ - << " -> converged=" << (fit.converged ? "yes" : "no") - << " amp=" << fit.amplitude << " t0=" << fit.t0 << "ns" - << " chi2/ndf=" << fit.chi2 << "/" << fit.ndf << " (" - << (fit.ndf > 0 ? fit.chi2 / fit.ndf : 0.0) << ")"; - waveforms.emplace_back(std::move(samples), pchannel, ch.hybrid_id_, ch.feb_id_, n_trig); - waveforms.back().setFitResult( - static_cast(fit.amplitude), static_cast(fit.t0), - static_cast(fit.chi2), fit.ndf, fit.converged); } ldmx_log(debug) << "Built " << waveforms.size() @@ -161,17 +121,6 @@ void SiStripWaveformBuilder::produce(framework::Event& event) { event.add(output_collection_, waveforms); } -void SiStripWaveformBuilder::onProcessEnd() { - const long n_ok = n_fit_attempted_ - n_fit_failed_; - const double fail_pct = - n_fit_attempted_ > 0 - ? 100.0 * static_cast(n_fit_failed_) / n_fit_attempted_ - : 0.0; - ldmx_log(info) << "Fit summary: " << n_fit_attempted_ << " attempted, " - << n_ok << " converged, " << n_fit_failed_ << " failed (" - << fail_pct << "%)"; -} - } // namespace tracking::reco DECLARE_PRODUCER(tracking::reco::SiStripWaveformBuilder) diff --git a/Tracking/src/Tracking/Reco/SiStripWaveformFitProcessor.cxx b/Tracking/src/Tracking/Reco/SiStripWaveformFitProcessor.cxx new file mode 100644 index 000000000..fa072be78 --- /dev/null +++ b/Tracking/src/Tracking/Reco/SiStripWaveformFitProcessor.cxx @@ -0,0 +1,179 @@ +#include "Tracking/Reco/SiStripWaveformFitProcessor.h" + +#include "Framework/Exception/Exception.h" +#include "Tracking/Digitization/SiStripConstants.h" +#include "Tracking/Digitization/StripPulseFitter.h" +#include "Tracking/Event/FittedSiStripHit.h" +#include "Tracking/Event/SiStripWaveform.h" +#include "Tracking/Reco/SiStripChannelMap.h" +#include "Tracking/Reco/TrackerPedestals.h" + +namespace tracking::reco { + +void SiStripWaveformFitProcessor::configure(framework::config::Parameters& ps) { + input_collection_ = + ps.get("input_collection", input_collection_); + input_pass_name_ = ps.get("input_pass_name", input_pass_name_); + output_collection_ = + ps.get("output_collection", output_collection_); + daq_map_file_ = ps.get("daq_map_file", daq_map_file_); + + t_scan_min_ns_ = ps.get("t_scan_min_ns", t_scan_min_ns_); + t_scan_max_ns_ = ps.get("t_scan_max_ns", t_scan_max_ns_); + t_scan_step_ns_ = ps.get("t_scan_step_ns", t_scan_step_ns_); + max_chi2_ndf_ = ps.get("max_chi2_ndf", max_chi2_ndf_); +} + +void SiStripWaveformFitProcessor::onProcessStart() { + using namespace tracking::digitization; + + if (daq_map_file_.empty()) { + EXCEPTION_RAISE("BadConfig", + "SiStripWaveformFitProcessor requires a daq_map_file."); + } + // Loads eagerly so a missing/malformed map fails here, at start-up, with the + // path in the message -- never as a silently empty event stream. + daq_map_ = TrackerDaqMap::fromJsonFile(daq_map_file_); + + pulse_shape_ = PulseShape::make(std::string(PULSE_SHAPE_NAME), + PEAKING_TIME_NS, SECOND_TIME_CONST_NS); + + ldmx_log(info) << "SiStripWaveformFitProcessor configured:" + << " daq_map='" << daq_map_file_ << "' (" << daq_map_.size() + << " sensors)" << " shape=" << PULSE_SHAPE_NAME + << " tp=" << PEAKING_TIME_NS << " ns" << " T scan [" + << t_scan_min_ns_ << ", " + << (t_scan_max_ns_ > 0.0 ? std::to_string(t_scan_max_ns_) + : std::string("auto")) + << "] ns" << " step=" << t_scan_step_ns_ << " ns"; +} + +void SiStripWaveformFitProcessor::produce(framework::Event& event) { + using namespace tracking::digitization; + + const auto& peds = + getCondition(TrackerPedestals::CONDITIONS_NAME); + + const auto& waveforms = event.getCollection( + input_collection_, input_pass_name_); + + std::vector hits; + hits.reserve(waveforms.size()); + + for (const auto& wf : waveforms) { + ++n_waveforms_; + + // ----------------------------------------------------------------------- + // Address first: an unmapped hybrid or an unbonded channel is dropped + // before paying for the fit. + // ----------------------------------------------------------------------- + const uint8_t feb = wf.getFebId(); + const uint8_t hybrid = wf.getHybridId(); + + if (!daq_map_.has(feb, hybrid)) { + ++n_unmapped_; + const uint16_t k = + static_cast((static_cast(feb) << 8) | hybrid); + if (unmapped_sensors_[k]++ == 0) { + ldmx_log(warn) << "No DAQ-map entry for feb=" << static_cast(feb) + << " hybrid=" << static_cast(hybrid) + << " -- dropping its waveforms (reported once)"; + } + continue; + } + + const auto& sensor = daq_map_.at(feb, hybrid); + const int16_t pchannel = wf.getPchannel(); + const int strip_id = channelmap::stripId( + pchannel, sensor.n_strips_, sensor.first_strip_, sensor.reversed_); + + // Reject channels that fall outside the bonded strip range (e.g. a read-out + // but unbonded APV). Counted so an unexpected layout shows up loudly. + if (strip_id < sensor.first_strip_ || + strip_id >= sensor.first_strip_ + sensor.n_strips_) { + ++n_out_of_range_; + continue; + } + + // ----------------------------------------------------------------------- + // Fit. The noise enters the chi2, so the fitter is per channel; it is a + // cheap value type over a shared pulse shape. + // ----------------------------------------------------------------------- + uint8_t apv_id, channel; + channelmap::apvChannelFromPchannel(pchannel, apv_id, channel); + const float noise = peds.noise(feb, hybrid, apv_id, channel); + + const auto& samples = wf.getSamples(); + const int n_samples = static_cast(samples.size()); + // Samples are pedestal-subtracted and lie on a uniform grid measured from + // sample 0, so T may peak anywhere from before sample 0 to the last sample. + const double t_scan_max = (t_scan_max_ns_ > 0.0) + ? t_scan_max_ns_ + : n_samples * SAMPLING_INTERVAL_NS; + + StripPulseFitter fitter(*pulse_shape_, + /*t0_offset_ns=*/0.0, SAMPLING_INTERVAL_NS, + /*pedestal_adc=*/0.0, + /*noise_sigma_adc=*/noise, t_scan_min_ns_, + t_scan_max, t_scan_step_ns_); + const auto fit = fitter.fit(samples); + ++n_fit_attempted_; + + ldmx_log(trace) << "fit feb=" << static_cast(feb) + << " hyb=" << static_cast(hybrid) + << " pch=" << pchannel << " nsamp=" << n_samples + << " noise=" << noise + << " -> converged=" << (fit.converged ? "yes" : "no") + << " amp=" << fit.amplitude << " t0=" << fit.t0 << "ns" + << " chi2/ndf=" << fit.chi2 << "/" << fit.ndf << " (" + << (fit.ndf > 0 ? fit.chi2 / fit.ndf : 0.0) << ")"; + + if (!fit.converged) { + ++n_unconverged_; + continue; + } + + if (max_chi2_ndf_ > 0.0 && fit.ndf > 0) { + if (fit.chi2 / fit.ndf > max_chi2_ndf_) { + ++n_bad_chi2_; + continue; + } + } + + hits.emplace_back( + sensor.layer_id_, strip_id, static_cast(fit.amplitude), + static_cast(fit.t0), static_cast(fit.chi2), fit.ndf, + /*track_id=*/-1, /*pdg_id=*/0, /*sim_hit_id=*/-1, + /*edep=*/0.f, noise); + ++n_hits_; + } + + ldmx_log(debug) << "Produced " << hits.size() << " FittedSiStripHits from " + << waveforms.size() << " waveforms"; + + event.add(output_collection_, hits); +} + +void SiStripWaveformFitProcessor::onProcessEnd() { + const long n_converged = n_fit_attempted_ - n_unconverged_; + const double fail_pct = + n_fit_attempted_ > 0 + ? 100.0 * static_cast(n_unconverged_) / n_fit_attempted_ + : 0.0; + ldmx_log(info) << "Fit summary: " << n_fit_attempted_ << " attempted, " + << n_converged << " converged, " << n_unconverged_ + << " failed (" << fail_pct << "%)"; + ldmx_log(info) << "SiStripWaveformFitProcessor summary: " << n_waveforms_ + << " waveforms -> " << n_hits_ << " hits (" << n_unmapped_ + << " unmapped, " << n_out_of_range_ << " out-of-range, " + << n_unconverged_ << " unconverged, " << n_bad_chi2_ + << " bad chi2/ndf)"; + for (const auto& [k, n] : unmapped_sensors_) { + ldmx_log(warn) << " unmapped feb=" << (k >> 8) << " hybrid=" << (k & 0xFF) + << ": " << n << " waveforms dropped"; + } +} + +} // namespace tracking::reco + +DECLARE_PRODUCER(tracking::reco::SiStripWaveformFitProcessor) diff --git a/Tracking/src/Tracking/Reco/StripClusterProcessor.cxx b/Tracking/src/Tracking/Reco/StripClusterProcessor.cxx index 43685e74f..884b31184 100644 --- a/Tracking/src/Tracking/Reco/StripClusterProcessor.cxx +++ b/Tracking/src/Tracking/Reco/StripClusterProcessor.cxx @@ -33,6 +33,7 @@ void StripClusterProcessor::configure( time_window_ns_ = parameters.get("time_window_ns", -1.0); neighbor_delta_t_ns_ = parameters.get("neighbor_delta_t_ns", -1.0); max_chi2_ndf_ = parameters.get("max_chi2_ndf", -1.0); + daq_map_file_ = parameters.get("daq_map_file", ""); } // --------------------------------------------------------------------------- @@ -44,6 +45,20 @@ void StripClusterProcessor::onProcessStart() { seed_threshold_, neighbor_threshold_, cluster_threshold_, NOISE_SIGMA_ADC, mean_time_ns_, time_window_ns_, neighbor_delta_t_ns_, max_chi2_ndf_); + // Optional DAQ map: build a layer_id -> n_strips lookup so the local-U centre + // offset can use the real per-sensor strip count for real data. Left empty + // for MC, in which case the fixed N_READOUT_STRIPS constant is used below. + layer_n_strips_.clear(); + if (!daq_map_file_.empty()) { + const auto map = TrackerDaqMap::fromJsonFile(daq_map_file_); + for (const auto& [key, sensor] : map.sensors()) { + layer_n_strips_[sensor.layer_id_] = sensor.n_strips_; + } + ldmx_log(info) << "StripClusterProcessor loaded DAQ map from '" + << daq_map_file_ << "' (" << layer_n_strips_.size() + << " layers) for centre-strip offsets"; + } + ldmx_log(info) << "StripClusterProcessor configured:" << " seed_thr=" << seed_threshold_ << " nbr_thr=" << neighbor_threshold_ << " cls_thr=" << cluster_threshold_ @@ -102,7 +117,14 @@ void StripClusterProcessor::produce(framework::Event& event) { // µm, etc. // ------------------------------------------------------------------- using namespace tracking::digitization; - const int n_int = N_READOUT_STRIPS / 2; // integer division = 383 + // Centre-strip offset: N/2 (integer division). For MC this is the fixed + // N_READOUT_STRIPS constant; for real data, if a DAQ map was supplied, + // use that sensor's real strip count so the local origin sits at its + // centre. + int n_strips = N_READOUT_STRIPS; + auto it_ns = layer_n_strips_.find(layer_id); + if (it_ns != layer_n_strips_.end()) n_strips = it_ns->second; + const int n_int = n_strips / 2; const double offset = static_cast(n_int); const double local_u = (cl.centroid_strip - offset) * READOUT_PITCH_MM; diff --git a/Tracking/src/Tracking/Reco/StripFitProcessor.cxx b/Tracking/src/Tracking/Reco/StripFitProcessor.cxx index 230529a51..211f454b6 100644 --- a/Tracking/src/Tracking/Reco/StripFitProcessor.cxx +++ b/Tracking/src/Tracking/Reco/StripFitProcessor.cxx @@ -75,7 +75,8 @@ void StripFitProcessor::produce(framework::Event& event) { raw.getLayerID(), raw.getStripID(), static_cast(result.amplitude), static_cast(result.t0), static_cast(result.chi2), result.ndf, raw.getTrackID(), - raw.getPdgID(), raw.getSimHitID(), raw.getEdep()); + raw.getPdgID(), raw.getSimHitID(), raw.getEdep(), + tracking::digitization::NOISE_SIGMA_ADC); ldmx_log(trace) << "Fitted: layer=" << raw.getLayerID() << " strip=" << raw.getStripID() diff --git a/Tracking/src/Tracking/Reco/TrackerDaqMap.cxx b/Tracking/src/Tracking/Reco/TrackerDaqMap.cxx new file mode 100644 index 000000000..3955f15fa --- /dev/null +++ b/Tracking/src/Tracking/Reco/TrackerDaqMap.cxx @@ -0,0 +1,77 @@ +#include "Tracking/Reco/TrackerDaqMap.h" + +#include +#include + +#include "Framework/Exception/Exception.h" + +namespace tracking::reco { + +TrackerDaqMap TrackerDaqMap::fromJsonFile(const std::string& path) { + std::ifstream in(path); + if (!in) { + EXCEPTION_RAISE( + "FileNotFound", + "TrackerDaqMap could not open DAQ map file '" + path + "'."); + } + + nlohmann::json doc; + try { + in >> doc; + } catch (const nlohmann::json::parse_error& e) { + EXCEPTION_RAISE("BadFormat", + "TrackerDaqMap failed to parse DAQ map file '" + path + + "': " + e.what()); + } + + if (!doc.contains("sensors") || !doc.at("sensors").is_array() || + doc.at("sensors").empty()) { + EXCEPTION_RAISE("BadFormat", "TrackerDaqMap file '" + path + + "' has no non-empty 'sensors' array."); + } + + TrackerDaqMap map; + for (const auto& s : doc.at("sensors")) { + for (const char* required : + {"feb", "hybrid", "layer_id", "n_strips", "first_strip", "reversed"}) { + if (!s.contains(required)) { + EXCEPTION_RAISE("BadFormat", "TrackerDaqMap file '" + path + + "' has a sensor entry missing '" + + required + "': " + s.dump()); + } + } + + const auto feb = s.at("feb").get(); + const auto hybrid = s.at("hybrid").get(); + const uint16_t k = key(feb, hybrid); + if (map.sensors_.count(k) != 0u) { + EXCEPTION_RAISE("BadFormat", "TrackerDaqMap file '" + path + + "' defines (feb=" + std::to_string(feb) + + ", hybrid=" + std::to_string(hybrid) + + ") more than once."); + } + + SensorInfo info; + info.layer_id_ = s.at("layer_id").get(); + info.n_strips_ = s.at("n_strips").get(); + info.first_strip_ = s.at("first_strip").get(); + info.reversed_ = s.at("reversed").get(); + map.sensors_[k] = info; + } + + return map; +} + +const TrackerDaqMap::SensorInfo& TrackerDaqMap::at(uint8_t feb, + uint8_t hybrid) const { + auto it = sensors_.find(key(feb, hybrid)); + if (it == sensors_.end()) { + EXCEPTION_RAISE("NotFound", "TrackerDaqMap has no entry for (feb=" + + std::to_string(feb) + + ", hybrid=" + std::to_string(hybrid) + + "); guard with has() before calling at()."); + } + return it->second; +} + +} // namespace tracking::reco