MUSIC unknown
Analysis for the MUSIC active-target ionization chamber
Loading...
Searching...
No Matches
CalibrateBeam.cpp File Reference
#include "CalibrateBeam.hpp"
#include <Rtypes.h>
#include <TH2D.h>
#include <utility>
Include dependency graph for CalibrateBeam.cpp:

Go to the source code of this file.

Classes

struct  StripPairSamples
 Paired (L, R) raw-ADC samples for one strip, collected UNGATED from events where both ends fire, plus the uncapped "shoulder" slice: short-side values from events where the LONG side reads low (charge went mostly to the short end). More...
struct  RidgeFit
 Everything SaveRidgeFitPlots needs to redraw a fit: the slice medians it was fitted to, the window they were taken from, and the resulting line. More...
struct  GainMatchState

Functions

Double_t StripTotalAdc (const UShort_t *l, const UShort_t *r, UShort_t s0, UShort_t s17, Int_t s)
Int_t GatePartner (Int_t s)
 The strip whose total is paired with strip s to gate it: the one before it, except strips 0 and 1 which share the (0, 1) pair because there is nothing before strip 0.
void RobustPeakSeed (const std::vector< Float_t > &v, Double_t &mode, Double_t &sigma)
 Robust peak/width estimate used to seed the beam-peak Gaussian fit, and as the fallback anchor when the fit fails.
Char_t LongSide (Int_t strip)
Bool_t IsCalibrated (const ChannelCal &c)
Double_t Gain (const ChannelCal &c)
BeamFit2D FindBeamGateStrips (const FileSpec &spec, Int_t sx, Int_t sy, const TString &run_label, const TString &plot_subdir, Bool_t save_plot=kTRUE)
 The beam gate for one strip: a 2D Gaussian on the (strip sx, strip sy) raw totals.
template<class T>
Double_t Median (std::vector< T > &v)
void CollectAnchorSamplesOneSubfile (const FileSpec &spec, const std::vector< ChannelCal > &chans, const BeamFit2D gate[18], std::vector< std::vector< Float_t > > &samples, StripPairSamples pairs[18])
 Every strip carries its own beam gate, on its own total and its neighbour's (GatePartner).
Bool_t TheilSenLine (const std::vector< Double_t > &x, const std::vector< Double_t > &y, Double_t &slope, Double_t &intercept)
 Theil-Sen robust line: slope = median of pairwise slopes, intercept = median residual.
Double_t RidgeShortAnchor (const std::vector< Float_t > &v_short, const std::vector< Float_t > &v_long, Double_t c_long, Double_t &slope_out, Double_t &intercept_out, RidgeFit &dbg)
 Robust line through the single-particle charge-sharing ridge, returning the fitted line's short-axis crossing, -intercept/slope (nonzero means the ridge was measurable), and the line itself through slope_out/intercept_out.
void SaveRidgeFitPlots (const StripPairSamples pairs[18], const RidgeFit dbg[18], const TString &plot_subdir)
 One plot per strip under <plot_subdir>/ridge, named ridge_s<NN>: the beam-gated long-vs-short plane, the slice medians the fit was actually given (black), and the fitted line (violet).
Double_t GmFindPeak (const std::vector< Float_t > &v, Double_t lo, Double_t hi, Double_t skip_frac)
 Histogram-mode peak finder, mirroring the notebook's find_peak(): histogram v over [lo, hi] with kGmBins bins, skip the first skip_frac of bins (to avoid the threshold pile), return the max-bin centre.
void ComputeLRGainMatch (std::vector< ChannelCal > &chans, const StripPairSamples pairs[18], const TString &plot_subdir)
void ReduceToAnchors (std::vector< ChannelCal > &chans, std::vector< std::vector< Float_t > > &samples, const StripPairSamples pairs[18], const TString &plot_subdir)
 Cathode uses median + IQR/1.349 (asymmetric tail, no clean peak).
void WriteCalibrationTree (TFile *dst, const std::vector< ChannelCal > &chans, const StripAlignmentResult *align)
 Writes a one-row calibration tree into the open file dst (typically a per-subfile .cal.root).
void WriteCalibrationToEvents (const FileSpec &spec, const std::vector< ChannelCal > &chans, const StripAlignmentResult *align)
 Writes the per-channel gain table (tree "calibration") into the subfile's own events file.
void SaveDynamicRangeOverlay (const FileSpec &spec, const TString &plot_subdir, const TString &file_label)
Double_t SmoothedPeakIn (TH1D *proj, Double_t lo, Double_t hi)
 Peak of a smoothed projection in [lo, hi]: the highest bin of the window, refined to sub-bin precision by the parabola through it and its two neighbours (kept within half a bin).
StripAlignmentResult FindStripCentroidAlignment (const FileSpec &spec, const TString &plot_subdir, const TString &file_label)
 Per-strip two-point alignment of the strip total.

Variables

const Int_t kMaxChannels = 36
const Long64_t kMinSamples = 200
const Long64_t kSampleCap = 20000
const Long64_t kPairCap = 100000
 (L, R) pair cap for the L/R gain-matching passes.
const Int_t kGmBins = 512
 L/R gain matching, after the check_LR notebook's recipe (37Cl_an_check_LR.ipynb, "Save 2-pass calibration" cell) minus the eta/position correction, with a short-end offset read off the ridge line:
const Double_t kRidgeShortMinFrac
 Ridge-slope short anchor.
const Double_t kRidgeShortMaxFrac
const Double_t kRidgeBandLo = 0.30
const Double_t kRidgeBandHi = 1.45
const Int_t kRidgeSlices = 60
const Int_t kRidgeMinSlices = 6
const Long64_t kRidgeMinPerSlice = 200
const Long64_t kRidgeNoiseFloor = 20
const Int_t kRidgeMinPts = 6
const Double_t kRidgeRatioLo = 0.50
const Double_t kRidgeRatioHi = 2.00
const Double_t kGmEsumLo = 0.8
const Double_t kGmEsumHi = 2.5

Function Documentation

◆ CollectAnchorSamplesOneSubfile()

void CollectAnchorSamplesOneSubfile ( const FileSpec & spec,
const std::vector< ChannelCal > & chans,
const BeamFit2D gate[18],
std::vector< std::vector< Float_t > > & samples,
StripPairSamples pairs[18] )

Every strip carries its own beam gate, on its own total and its neighbour's (GatePartner).

A strip is therefore anchored on events that were beam-like THERE, rather than on events that were beam-like at strips 1 and 2 and whatever they happened to be doing further down the chamber.

Definition at line 292 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ ComputeLRGainMatch()

void ComputeLRGainMatch ( std::vector< ChannelCal > & chans,
const StripPairSamples pairs[18],
const TString & plot_subdir )

Definition at line 1005 of file CalibrateBeam.cpp.

Referenced by ReduceToAnchors().

◆ FindBeamGateStrips()

BeamFit2D FindBeamGateStrips ( const FileSpec & spec,
Int_t sx,
Int_t sy,
const TString & run_label,
const TString & plot_subdir,
Bool_t save_plot = kTRUE )

The beam gate for one strip: a 2D Gaussian on the (strip sx, strip sy) raw totals.

Gating a strip on itself and its neighbour selects single-particle events in THAT strip; the gate is on totals, so it says nothing about how the charge divides between the two ends, which is what the anchor fit reads.

Use 1024 bins (previously 256) so each bin is narrower (16 ADC for 37Cl). The sigma floor in ComputeMoments is 2 * bin_width; at the old 64 ADC/bin it clamped sigma to 128 ADC, hiding the true beam width and washing out correlation (rho ≈ 0.09 even though strip1/strip2 track the same beam).

With 1024 bins instead of 256, scale from 10→40 to keep the ADC seed window roughly 640 ADC (10 bins × 16384/256 = 640; 40 bins × 16384/1024 = 640).

Definition at line 80 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ FindStripCentroidAlignment()

StripAlignmentResult FindStripCentroidAlignment ( const FileSpec & spec,
const TString & plot_subdir,
const TString & file_label )

Per-strip two-point alignment of the strip total.

Decodes every event with the per-channel gains already on disk, projects each strip's total, and locates two peaks in the smoothed projection: the beam (B, the highest bin in [0.6, 1.5]) and the two-particle pile-up (P, the highest bin in [1.7 B, 2.3 B], a window placed from the beam found so a strip whose beam sits at 0.9 is searched around 1.8 and the 3-particle shoulder at 3 B stays out). Two beam particles in one event deposit exactly twice what one does, so the total is put on the scale that has B at 1.0 and P at 2.0:

total' = factor * total + offset,
factor = 1 / (P - B),  offset = 1 - factor * B.

A beam-only anchor cannot separate a gain from a pedestal; the pile-up peak is the second point that fixes both, and the summed spectra then show the beam at 1 and the pile-up at 2 on every strip instead of pile-up peaks spread over 1.8-2.2 by per-channel baseline offsets. Where no pile-up peak is found the strip keeps a factor-only alignment, 1 / B, with no offset. The offset is a per-strip quantity; EnergyView carries it on the long end (the unsegmented value on strips 0 and 17), so a strip's total is still the plain sum of its ends.

Uses ALL events (not beam-gated): the beam dominates every strip's histogram by a wide margin, and the pile-up peak is the next feature. The peaks come from 0.005 a.u. bins with parabolic sub-bin refinement (SmoothedPeakIn); a percent-level alignment needs peaks located well inside a percent, which the earlier 0.05 a.u. bins could not do.

Definition at line 1341 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ Gain()

Double_t Gain ( const ChannelCal & c)

◆ GatePartner()

Int_t GatePartner ( Int_t s)
inline

The strip whose total is paired with strip s to gate it: the one before it, except strips 0 and 1 which share the (0, 1) pair because there is nothing before strip 0.

Definition at line 21 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile(), and CollectAnchorSamplesOneSubfile().

◆ GmFindPeak()

Double_t GmFindPeak ( const std::vector< Float_t > & v,
Double_t lo,
Double_t hi,
Double_t skip_frac )

Histogram-mode peak finder, mirroring the notebook's find_peak(): histogram v over [lo, hi] with kGmBins bins, skip the first skip_frac of bins (to avoid the threshold pile), return the max-bin centre.

Returns 0 when empty.

Definition at line 733 of file CalibrateBeam.cpp.

◆ IsCalibrated()

Bool_t IsCalibrated ( const ChannelCal & c)

◆ LongSide()

◆ Median()

template<class T>
Double_t Median ( std::vector< T > & v)
inline

◆ ReduceToAnchors()

void ReduceToAnchors ( std::vector< ChannelCal > & chans,
std::vector< std::vector< Float_t > > & samples,
const StripPairSamples pairs[18],
const TString & plot_subdir )

Cathode uses median + IQR/1.349 (asymmetric tail, no clean peak).

Every other channel is anchored on the modal bin of its beam-peak histogram (RobustPeakSeed: 100 bins over the 5th-95th percentile of the samples), with no fit on top and no width: nothing downstream uses one.

Definition at line 1027 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ RidgeShortAnchor()

Double_t RidgeShortAnchor ( const std::vector< Float_t > & v_short,
const std::vector< Float_t > & v_long,
Double_t c_long,
Double_t & slope_out,
Double_t & intercept_out,
RidgeFit & dbg )

Robust line through the single-particle charge-sharing ridge, returning the fitted line's short-axis crossing, -intercept/slope (nonzero means the ridge was measurable), and the line itself through slope_out/intercept_out.

The caller turns slope and intercept into the short anchor and the short-end offset against the long end's modal peak (see ComputeLRGainMatch). c_long is only the seed that places the short window and the long band; the line's intercept must land near it. Returns 0 when the ridge is not measurable.

Definition at line 553 of file CalibrateBeam.cpp.

◆ RobustPeakSeed()

void RobustPeakSeed ( const std::vector< Float_t > & v,
Double_t & mode,
Double_t & sigma )

Robust peak/width estimate used to seed the beam-peak Gaussian fit, and as the fallback anchor when the fit fails.

Histograms only the 5th-95th-percentile core so outlier ADC values can't stretch the binning, takes the modal bin centre as the peak and IQR/1.349 as the width. Peak-like: unlike the raw sample mean it is not pulled up by the straggling beam-dE tail.

Definition at line 234 of file CalibrateBeam.cpp.

Referenced by ReduceToAnchors().

◆ SaveDynamicRangeOverlay()

void SaveDynamicRangeOverlay ( const FileSpec & spec,
const TString & plot_subdir,
const TString & file_label )

Definition at line 1180 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ SaveRidgeFitPlots()

void SaveRidgeFitPlots ( const StripPairSamples pairs[18],
const RidgeFit dbg[18],
const TString & plot_subdir )

One plot per strip under <plot_subdir>/ridge, named ridge_s<NN>: the beam-gated long-vs-short plane, the slice medians the fit was actually given (black), and the fitted line (violet).

Drawn for every strip, including the ones whose fit was rejected, so a bad ridge can be seen rather than inferred from the slope in the log.

Definition at line 665 of file CalibrateBeam.cpp.

Referenced by ComputeLRGainMatch().

◆ SmoothedPeakIn()

Double_t SmoothedPeakIn ( TH1D * proj,
Double_t lo,
Double_t hi )

Peak of a smoothed projection in [lo, hi]: the highest bin of the window, refined to sub-bin precision by the parabola through it and its two neighbours (kept within half a bin).

No fit. Returns 0 when the window is empty or its highest bin sits on the window edge, which means the peak itself lies outside the window. The window is what selects the feature: picking the local maximum nearest a target instead let a ripple of the smoothed tail at 2.1-2.2 a.u. win over the real pile-up peak at 1.8 on the strips whose beam sits low, which threw their alignment by 20%.

Definition at line 1285 of file CalibrateBeam.cpp.

Referenced by FindStripCentroidAlignment().

◆ StripTotalAdc()

Double_t StripTotalAdc ( const UShort_t * l,
const UShort_t * r,
UShort_t s0,
UShort_t s17,
Int_t s )
inline

Definition at line 9 of file CalibrateBeam.cpp.

Referenced by CollectAnchorSamplesOneSubfile(), and FindBeamGateStrips().

◆ TheilSenLine()

Bool_t TheilSenLine ( const std::vector< Double_t > & x,
const std::vector< Double_t > & y,
Double_t & slope,
Double_t & intercept )

Theil-Sen robust line: slope = median of pairwise slopes, intercept = median residual.

Unlike OLS ("pol1" via a chi-square fit) a single outlying slice median cannot bend the slope: it corrects the intercept, complementing the slice median's own robustness. Returns kFALSE when the slice set is too small or degenerate (all-same x).

Definition at line 519 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ WriteCalibrationToEvents()

void WriteCalibrationToEvents ( const FileSpec & spec,
const std::vector< ChannelCal > & chans,
const StripAlignmentResult * align )

Writes the per-channel gain table (tree "calibration") into the subfile's own events file.

No per-event calibrated tree is produced: downstream readers recover a.u. on the fly via gain x raw ADC (EnergyView), so the raw events tree plus this one-row gain table fully determine every calibrated value.

Definition at line 1159 of file CalibrateBeam.cpp.

Referenced by CalibrateBeam::CalibrateBeamOneSubfile().

◆ WriteCalibrationTree()

void WriteCalibrationTree ( TFile * dst,
const std::vector< ChannelCal > & chans,
const StripAlignmentResult * align )

Writes a one-row calibration tree into the open file dst (typically a per-subfile .cal.root).

align carries the per-strip alignment factors; pass nullptr when they have not been computed yet (identity is written).

Definition at line 1064 of file CalibrateBeam.cpp.

Referenced by WriteCalibrationToEvents().

Variable Documentation

◆ kGmBins

const Int_t kGmBins = 512

L/R gain matching, after the check_LR notebook's recipe (37Cl_an_check_LR.ipynb, "Save 2-pass calibration" cell) minus the eta/position correction, with a short-end offset read off the ridge line:

Pass 1 — per-side anchors and the short-end offset:

  • a fixed deposit divided between the two ends traces the line long/C_long + (short - S0)/C_short = 1 in the (short, long) plane. S0 is an additive offset the short end's energy carries whenever that end fired: on 37Cl the R-chain short ends read 100-230 ADC more than the charge they collected (R5 231, R3 137, R13 127 ADC in run 97; R9 about 0), and the L-chain short ends a few tens of ADC, which is what a trigger threshold looks like from here. An end that did not fire reads 0 and carries no offset. Without S0 the strip total is bimodal: events without a short hit sum to C_long, events with one to C_long + S0 x C_long/C_short, and the second peak is the high-side tail the odd strips showed.
  • LONG anchor C_long = the long end's own modal beam peak (ReduceToAnchors): most beam events leave the short end silent, so that peak is the full deposit on the long end to within the charge a silent short end can hide, which is below its trigger threshold.
  • the robust line through the beam-gated pairs (RidgeShortAnchor) has slope -k = -C_long/C_short and intercept a = C_long + k S0, so C_short = C_long / k, S0 = (a - C_long) / k. Both come from the same line; S0 is the sideways shift of the ridge relative to the long-only peak.
  • where the ridge is not measurable the short anchor is the parity's median C_short/C_long (a preamp gain ratio, so it transfers across strips) times C_long, or the global median, and S0 is the parity's median offset in units of C_long, or 0.
  • gain = TARGET / anchor per side, with TARGET = 1.0 a.u. (the notebook uses 1000 ADC; only the overall scale differs); the decode applies gain x (ADC - S0) to an end that fired (EnergyView::Decode).

Pass 2 — check, do not correct: with the offset removed the summed beam peak of the pairs on the ridge (short inside the fit window) equals the long-only peak, 1.0 a.u., by construction; a departure exposes a bad ridge fit or a poor fallback. Censored events (short end 0) sum to the long end alone and are left out so they cannot mask a bad fit.

  • the old pass 2 rescaled both gains by 1/peak, which left the L/R ratio untouched, so a wrong short anchor misallocated charge within a strip but no longer shifts the strip total.
  • when IGNORE_SHORT_STRIPS is set the decode keeps only the long end, so the long end alone is normalised to 1.0 and pass 1 already provides it.

Pairs are collected UNGATED (all events with both ends firing) because the ridge needs the off-centre crossings that a beam gate removes. Strips whose ridge cannot be fitted fall back to the MEDIAN anchor of the strips that could; that fallback assumes a common electronics scale across strips, which the measured ridge slopes can be used to check rather than assume.

Definition at line 453 of file CalibrateBeam.cpp.

Referenced by GmFindPeak().

◆ kGmEsumHi

const Double_t kGmEsumHi = 2.5

Definition at line 493 of file CalibrateBeam.cpp.

◆ kGmEsumLo

const Double_t kGmEsumLo = 0.8

Definition at line 492 of file CalibrateBeam.cpp.

◆ kMaxChannels

const Int_t kMaxChannels = 36

Definition at line 6 of file CalibrateBeam.cpp.

Referenced by WriteCalibrationTree().

◆ kMinSamples

const Long64_t kMinSamples = 200

◆ kPairCap

const Long64_t kPairCap = 100000

(L, R) pair cap for the L/R gain-matching passes.

Larger than kSampleCap because the short-side "shoulder" anchor is found in a narrow slice of the pairs (long side ≈ 300-350 ADC) that only a small fraction of events populate.

Definition at line 28 of file CalibrateBeam.cpp.

Referenced by CollectAnchorSamplesOneSubfile().

◆ kRidgeBandHi

const Double_t kRidgeBandHi = 1.45

Definition at line 480 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeBandLo

const Double_t kRidgeBandLo = 0.30

Definition at line 479 of file CalibrateBeam.cpp.

◆ kRidgeMinPerSlice

const Long64_t kRidgeMinPerSlice = 200

Definition at line 483 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeMinPts

const Int_t kRidgeMinPts = 6

Definition at line 487 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeMinSlices

const Int_t kRidgeMinSlices = 6

Definition at line 482 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeNoiseFloor

const Long64_t kRidgeNoiseFloor = 20

Definition at line 486 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeRatioHi

const Double_t kRidgeRatioHi = 2.00

Definition at line 491 of file CalibrateBeam.cpp.

◆ kRidgeRatioLo

const Double_t kRidgeRatioLo = 0.50

Definition at line 490 of file CalibrateBeam.cpp.

◆ kRidgeShortMaxFrac

const Double_t kRidgeShortMaxFrac
Initial value:
=
0.5

Definition at line 477 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor(), and SaveRidgeFitPlots().

◆ kRidgeShortMinFrac

const Double_t kRidgeShortMinFrac
Initial value:
=
0.15

Ridge-slope short anchor.

In the long-vs-short plane a fixed deposit divided between the two ends of a strip traces ADC_long/C_long + ADC_short/C_short = 1, a line of slope -C_long/C_short. The slope is the same for every deposit energy, so the ridge direction alone gives the anchor ratio and C_short need never be observed directly – which matters because the beam is collimated onto the long end and rarely deposits its full charge on the short one. The plane also contains pile-up bands at 2x, 3x the single-particle deposit; the 2-particle band falls into the single-particle window once short > ~0.55*C_short, so the fit is capped well below that.

The fit window also starts well above the short end's threshold. Just above threshold the short reading carries no position information: on the ends that are censored for most beam events (the L end of the even strips on 37Cl) the pairs that survive pile up in a flat blob at 0.04-0.15 C_long with the long end at its full deposit, and on the other ends the slice medians bend upward there too. A slice fit that sees that blob comes out nearly flat (slopes of -0.1 to -0.5 on strips 2 and 4, ratios of 2-10) and puts the normed band far below -1. Only the falling ridge beyond the blob is fitted; it is sparse on the censored ends, so the slice count adapts to the events available.

Definition at line 475 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kRidgeSlices

const Int_t kRidgeSlices = 60

Definition at line 481 of file CalibrateBeam.cpp.

Referenced by RidgeShortAnchor().

◆ kSampleCap

const Long64_t kSampleCap = 20000

Definition at line 23 of file CalibrateBeam.cpp.

Referenced by CollectAnchorSamplesOneSubfile().