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

gondola_core/events/
rb_waveform.rs
