gondola_core/events/
rb_waveform.rs

1// This file is part of gaps-online-software and published 
2// under the GPLv3 license
3
4use crate::prelude::*;
5
6/// Waveform container for Tof waveforms
7/// This holds the waveforms for both 
8/// paddle ends. Fields are available to 
9/// hold calibrated waveforms, however,
10/// only adc will be saved to disk.
11#[derive(Debug, Clone, PartialEq)]
12#[cfg_attr(feature="pybindings", pyclass)]
13pub struct RBWaveform {
14  pub event_id      : u32,
15  pub rb_id         : u8,
16  ///// FIXME - this is form 0-8, but should it be from 1-9?
17  pub rb_channel_a  : u8,
18  pub rb_channel_b  : u8,
19  /// DRS4 stop cell
20  pub stop_cell     : u16,
21  pub adc_a         : Vec<u16>,
22  pub adc_b         : Vec<u16>,
23  // FIXME - to be compatible with Antarctica data from 
24  // 2024/25, we do not serialize the paddle id in this version
25  // HOWEVER, let's change that in v0.12
26  pub paddle_id     : u8,
27  pub voltages_a    : Vec<f32>,
28  pub nanoseconds_a : Vec<f32>,
29  pub voltages_b    : Vec<f32>,
30  pub nanoseconds_b : Vec<f32>
31}
32
33impl RBWaveform {
34  
35  pub fn new() -> Self {
36    Self {
37      event_id       : 0,
38      rb_id          : 0,
39      rb_channel_a   : 0,
40      rb_channel_b   : 0,
41      stop_cell      : 0,
42      paddle_id      : 0,
43      adc_a          : Vec::<u16>::new(),
44      voltages_a     : Vec::<f32>::new(),
45      nanoseconds_a  : Vec::<f32>::new(),
46      adc_b          : Vec::<u16>::new(),
47      voltages_b     : Vec::<f32>::new(),
48      nanoseconds_b  : Vec::<f32>::new()
49    }
50  }
51
52  /// Calculate the time in ns for which the waveform is 
53  /// above a certain threshold for paddle end A
54  pub fn time_over_threshold_a(&self, threshold : f32) -> f32 {
55    let mut tot : f32 = 0.0;
56    for k in 1..self.voltages_a.len() {
57      if self.voltages_a[k] > threshold {
58        tot += self.nanoseconds_a[k] - self.nanoseconds_a[k-1];
59      }
60    }
61    return tot;
62  }
63
64  /// Calculate the time in ns for which the waveform is 
65  /// above a certain threshold for paddle end B
66  pub fn time_over_threshold_b(&self, threshold : f32) -> f32 {
67    let mut tot : f32 = 0.0;
68    for k in 1..self.voltages_b.len() {
69      if self.voltages_b[k] > threshold {
70        tot += self.nanoseconds_b[k] - self.nanoseconds_b[k-1];
71      }
72    }
73    return tot;
74  }
75  
76  pub fn charge_a_below_500(&self) -> f32 {
77
78    let mut total_area = 0.0f32;
79    for i in 0..(self.nanoseconds_a.len() - 1) {
80        let h = self.nanoseconds_a[i + 1] - self.nanoseconds_a[i]; // The width of the trapezoid.
81        let mut v1 = self.voltages_a[i];
82        let mut v2 = self.voltages_a[i + 1];
83        if v1 > 500.0 {
84          v1 = 500.0;
85        }
86        if v2 > 500.0 { 
87          v2 = 500.0; 
88        }
89        let area = h * (v1 + v2) / 2.0;
90        total_area += area;
91    } 
92    total_area
93  }
94  
95  pub fn charge_b_below_500(&self) -> f32 {
96
97    let mut total_area = 0.0f32;
98    for i in 0..(self.nanoseconds_b.len() - 1) {
99        let h = self.nanoseconds_b[i + 1] - self.nanoseconds_b[i]; // The width of the trapezoid.
100        let mut v1 = self.voltages_b[i];
101        let mut v2 = self.voltages_b[i + 1];
102        if v1 > 500.0 {
103          v1 = 500.0;
104        }
105        if v2 > 500.0 { 
106          v2 = 500.0; 
107        }
108        let area = h * (v1 + v2) / 2.0;
109        total_area += area;
110    } 
111    total_area
112  }
113
114  //pub fn charge_b_trap_bottom_fwhm(&self) -> f32 {
115  //  let mut total_area = 0.0f32;
116  //  for i in 0..(self.nanoseconds_a.len() - 1) {
117  //      let h = self.nanoseconds_a[i + 1] - self.nanoseconds_a[i]; // The width of the trapezoid.
118  //      let area = h * (self.voltages_a[i] + self.voltages_a[i + 1]) / 2.0;
119  //      total_area += area;
120  //  } 
121  //  total_area
122  //}
123  
124  pub fn charge_a_trap(&self) -> f32 {
125    let mut total_area = 0.0f32;
126    for i in 0..(self.nanoseconds_a.len() - 1) {
127        let h = self.nanoseconds_a[i + 1] - self.nanoseconds_a[i]; // The width of the trapezoid.
128        let area = h * (self.voltages_a[i] + self.voltages_a[i + 1]) / 2.0;
129        total_area += area;
130    } 
131    total_area
132  }
133  
134  pub fn charge_b_trap(&self) -> f32 {
135    let mut total_area = 0.0f32;
136    for i in 0..(self.nanoseconds_b.len() - 1) {
137        let h = self.nanoseconds_b[i + 1] - self.nanoseconds_b[i]; // The width of the trapezoid.
138        let area = h * (self.voltages_b[i] + self.voltages_b[i + 1]) / 2.0;
139        total_area += area;
140    } 
141    total_area
142  }
143
144  pub fn guess_max_peak_a(&self) -> f32 {
145    let mut current_max = -1000.0f32;
146    for k in 20..self.voltages_a.len() {
147      if self.voltages_a[k] > current_max {
148        current_max = self.voltages_a[k];
149      }
150    }
151    current_max
152  }
153  
154  pub fn guess_max_peak_b(&self) -> f32 {
155    let mut current_max = -1000.0f32;
156    for k in 20..self.voltages_b.len() {
157      if self.voltages_b[k] > current_max {
158        current_max = self.voltages_b[k];
159      }
160    }
161    current_max
162  }
163
164  pub fn subtract_pedestals(&mut self) {
165    let (ped_a, _ped_a_err) = calculate_pedestal(&self.voltages_a,
166                                                 10.0,
167                                                 700,
168                                                 200);
169    let (ped_b, _ped_b_err) = calculate_pedestal(&self.voltages_b,
170                                                 10.0,
171                                                 700,
172                                                 200);
173    for k in 0..self.voltages_a.len() {
174      self.voltages_a[k] = self.voltages_a[k] - ped_a;
175      self.voltages_b[k] = self.voltages_b[k] - ped_b;
176    }
177  }
178
179  /// Apply a RB calibration to the waveform, filling the voltages and 
180  /// nanoseconds fields
181  pub fn calibrate(&mut self, cali : &RBCalibrations) -> Result<(), CalibrationError>  {
182    if cali.rb_id != self.rb_id {
183      error!("Calibration is for board {}, but wf is for {}", cali.rb_id, self.rb_id);
184      return Err(CalibrationError::WrongBoardId);
185    }
186    let mut voltages = vec![0.0f32;1024];
187    let mut nanosecs = vec![0.0f32;1024];
188    cali.voltages(self.rb_channel_a as usize + 1,
189                  self.stop_cell as usize,
190                  &self.adc_a,
191                  &mut voltages);
192    self.voltages_a = voltages.clone();
193    cali.nanoseconds(self.rb_channel_a as usize + 1,
194                     self.stop_cell as usize,
195                     &mut nanosecs);
196    self.nanoseconds_a = nanosecs.clone();
197    cali.voltages(self.rb_channel_b as usize + 1,
198                  self.stop_cell as usize,
199                  &self.adc_b,
200                  &mut voltages);
201    self.voltages_b = voltages;
202    cali.nanoseconds(self.rb_channel_b as usize + 1,
203                     self.stop_cell as usize,
204                     &mut nanosecs);
205    self.nanoseconds_b = nanosecs;
206    Ok(())
207  }
208
209  /// Apply Jamie's simple spike filter to the calibrated voltages
210  pub fn apply_spike_filter(&mut self) {
211    clean_spikes(&mut self.voltages_a, true);
212    clean_spikes(&mut self.voltages_b, true);
213  }
214}
215
216impl TofPackable for RBWaveform {
217  const TOF_PACKET_TYPE : TofPacketType = TofPacketType::RBWaveform;
218}
219
220impl Serialization for RBWaveform {
221  const HEAD               : u16    = 43690; //0xAAAA
222  const TAIL               : u16    = 21845; //0x5555
223  
224  fn from_bytestream(stream : &Vec<u8>, pos : &mut usize)
225    -> Result<Self, SerializationError> {
226    let mut wf           = RBWaveform::new();
227    if parse_u16(stream, pos) != Self::HEAD {
228      error!("The given position {} does not point to a valid header signature of {}", pos, Self::HEAD);
229      return Err(SerializationError::HeadInvalid {});
230    }
231    wf.event_id          = parse_u32(stream, pos);
232    wf.rb_id             = parse_u8 (stream, pos);
233    wf.rb_channel_a      = parse_u8 (stream, pos);
234    wf.rb_channel_b      = parse_u8 (stream, pos);
235    wf.stop_cell         = parse_u16(stream, pos);
236    wf.paddle_id         = parse_u8 (stream, pos);
237    if stream.len() < *pos+2*NWORDS {
238      return Err(SerializationError::StreamTooShort);
239    }
240    let data_a           = &stream[*pos..*pos+2*NWORDS];
241    wf.adc_a             = u8_to_u16(data_a);
242    *pos += 2*NWORDS;
243    if stream.len() < *pos+2*NWORDS {
244      return Err(SerializationError::StreamTooShort);
245    }
246    let data_b           = &stream[*pos..*pos+2*NWORDS];
247    wf.adc_b             = u8_to_u16(data_b);
248    *pos += 2*NWORDS;
249    if parse_u16(stream, pos) != Self::TAIL {
250      error!("The given position {} does not point to a tail signature of {}", pos, Self::TAIL);
251      return Err(SerializationError::TailInvalid);
252    }
253    Ok(wf)
254  }
255
256  fn to_bytestream(&self) -> Vec<u8> {
257    let mut stream = Vec::<u8>::new();
258    stream.extend_from_slice(&Self::HEAD.to_le_bytes());
259    stream.extend_from_slice(&self.event_id.to_le_bytes());
260    stream.extend_from_slice(&self.rb_id.to_le_bytes());
261    stream.extend_from_slice(&self.rb_channel_a.to_le_bytes());
262    stream.extend_from_slice(&self.rb_channel_b.to_le_bytes());
263    stream.extend_from_slice(&self.stop_cell.to_le_bytes());
264    stream.push(self.paddle_id);
265    if self.adc_a.len() != 0 {
266      for k in 0..NWORDS {
267        stream.extend_from_slice(&self.adc_a[k].to_le_bytes());  
268      }
269    }
270    if self.adc_b.len() != 0 {
271      for k in 0..NWORDS {
272        stream.extend_from_slice(&self.adc_b[k].to_le_bytes());  
273      }
274    }
275    stream.extend_from_slice(&Self::TAIL.to_le_bytes());
276    stream
277  }
278}
279
280impl fmt::Display for RBWaveform {
281  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
282    let mut repr = String::from("<RBWaveform:");
283    repr += &(format!("\n  Event ID  : {}", self.event_id));
284    repr += &(format!("\n  RB        : {}", self.rb_id));
285    repr += &(format!("\n  ChannelA  : {}", self.rb_channel_a));
286    repr += &(format!("\n  ChannelB  : {}", self.rb_channel_b));
287    repr += &(format!("\n  Paddle ID : {}", self.paddle_id));
288    repr += &(format!("\n  Stop cell : {}", self.stop_cell));
289    if self.adc_a.len() >= 273 {
290      repr += &(format!("\n  adc [A] [{}]      : .. {} {} {} ..",self.adc_a.len(), self.adc_a[270], self.adc_a[271], self.adc_a[272]));
291    } else {
292      repr += &(String::from("\n  adc [A] [EMPTY]"));
293    }
294    if self.adc_b.len() >= 273 {
295      repr += &(format!("\n  adc [B] [{}]      : .. {} {} {} ..",self.adc_b.len(), self.adc_b[270], self.adc_b[271], self.adc_b[272]));
296    } else {
297      repr += &(String::from("\n  adc [B] [EMPTY]"));
298    }
299    write!(f, "{}", repr)
300  }
301}
302
303//---------------------------------------------------
304
305#[cfg(feature="pybindings")]
306#[pymethods]
307impl RBWaveform {
308  
309  /// Paddle ID of this wveform (1-160)
310  #[getter]
311  fn get_paddle_id(&self) -> u8 {
312    self.paddle_id
313  }
314
315  #[getter]
316  fn max_peak_a_guess(&self) -> f32 {
317    self.guess_max_peak_a()
318  }
319  
320  #[getter]
321  fn max_peak_b_guess(&self) -> f32 {
322    self.guess_max_peak_b()
323  }
324
325  #[getter]
326  fn get_rb_id(&self) -> u8 {
327    self.rb_id
328  }
329  
330  #[getter]
331  fn get_event_id(&self) -> u32 {
332    self.event_id
333  }
334  
335  #[getter]
336  fn get_rb_channel_a(&self) -> u8 {
337    self.rb_channel_a
338  }
339  
340  #[getter]
341  fn get_rb_channel_b(&self) -> u8 {
342    self.rb_channel_b
343  }
344  
345  #[getter]
346  fn get_stop_cell(&self) -> u16 {
347    self.stop_cell
348  }
349 
350  // FIXME - make this consistent 
351  #[getter]
352  fn get_adc_a<'_py>(&self, py: Python<'_py>) ->  Bound<'_py, PyArray1<u16>> {
353    //let arr = PyArray1::<u16>::from_slice(py, self.adc_a.as_slice());
354    self.adc_a.to_pyarray(py)
355    //Ok(arr)
356  }
357  
358  #[getter]
359  fn get_adc_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<u16>>> {
360    let arr = PyArray1::<u16>::from_slice(py, self.adc_b.as_slice());
361    Ok(arr)
362  }
363  
364  #[getter]
365  fn get_voltages_a<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
366    let arr = PyArray1::<f32>::from_slice(py, self.voltages_a.as_slice());
367    Ok(arr)
368  }
369
370  #[getter]
371  fn get_times_a<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
372    let arr    = PyArray1::<f32>::from_slice(py, self.nanoseconds_a.as_slice());
373    Ok(arr)
374  }
375
376  #[getter]
377  fn get_voltages_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
378    let arr = PyArray1::<f32>::from_slice(py, self.voltages_b.as_slice());
379    Ok(arr)
380  }
381
382  /// Time over threshold - waveform needs to be 
383  /// calibrated. 
384  /// Paddle end A
385  ///
386  /// # Arguments:
387  ///   * threshold : value in mV
388  fn get_tot_a(&self, threshold : f32) -> f32 {
389    self.time_over_threshold_a(threshold)
390  }
391  
392  /// Time over threshold - waveform needs to be 
393  /// calibrated. 
394  /// Paddle end B
395  ///
396  /// # Arguments:
397  ///   * threshold : value in mV
398  fn get_tot_b(&self, threshold : f32) -> f32 {
399    self.time_over_threshold_b(threshold)
400  }
401
402  #[pyo3(name="subtract_pedestals")]
403  fn subtract_pedestals_py(&mut self) {
404    self.subtract_pedestals()
405  }
406
407  #[getter]
408  fn get_charge_a_trap(&self) -> f32 {
409    self.charge_a_trap()
410  }
411  
412  #[getter]
413  fn get_charge_a_below_500_trap(&self) -> f32 {
414    self.charge_a_below_500() 
415  } 
416  
417  #[getter]
418  fn get_charge_b_below_500_trap(&self) -> f32 {
419    self.charge_b_below_500() 
420  } 
421
422
423  #[getter]
424  fn get_charge_b_trap(&self) -> f32 {
425    self.charge_b_trap()
426  }
427
428  #[getter]
429  fn get_times_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
430    let arr = PyArray1::<f32>::from_slice(py, self.nanoseconds_b.as_slice());
431    Ok(arr)
432  }
433  
434  /// Apply the readoutboard calibration to convert adc/bins
435  /// to millivolts and nanoseconds
436  #[pyo3(name="calibrate")]
437  fn calibrate_py(&mut self, cali : &RBCalibrations) -> PyResult<()> {
438    match self.calibrate(&cali) {
439      Ok(_) => {
440        return Ok(());
441      }
442      Err(err) => {
443        return Err(PyValueError::new_err(err.to_string()));
444      }
445    }
446  }
447  
448  #[pyo3(name="apply_spike_filter")]
449  fn apply_spike_filter_py(&mut self) {
450    self.apply_spike_filter();
451  }
452}
453
454#[cfg(feature="pybindings")]
455pythonize_packable!(RBWaveform);
456
457// needs to fix something up
458//#[cfg(feature="pybindings")]
459//pythonize_telemetry!(RBWaveform);
460
461//---------------------------------------------------
462
463#[cfg(feature = "random")]
464impl FromRandom for RBWaveform {
465    
466  fn from_random() -> Self {
467    let mut wf      = Self::new();
468    let mut rng     = rand::rng();
469    wf.event_id     = rng.random::<u32>();
470    wf.rb_id        = rng.random::<u8>();
471    wf.rb_channel_a = rng.random::<u8>();
472    wf.rb_channel_b = rng.random::<u8>();
473    wf.stop_cell    = rng.random::<u16>();
474    wf.paddle_id    = rng.random::<u8>();
475    let random_numbers_a: Vec<u16> = (0..NWORDS).map(|_| rng.random()).collect();
476    wf.adc_a        = random_numbers_a;
477    let random_numbers_b: Vec<u16> = (0..NWORDS).map(|_| rng.random()).collect();
478    wf.adc_b        = random_numbers_b;
479    wf
480  }
481}
482
483//---------------------------------------------------
484
485#[test]
486#[cfg(feature = "random")]
487fn pack_rbwaveform() {
488  for _ in 0..100 {
489    let wf   = RBWaveform::from_random();
490    let test : RBWaveform = wf.pack().unpack().unwrap();
491    assert_eq!(wf, test);
492  }
493}
494
495
gondola_core/events/rb_waveform.rs

gondola_core/events/
rb_waveform.rs
