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  const SIZE               : usize  = 13 + (4*NWORDS);
224
225  fn from_bytestream(stream : &Vec<u8>, pos : &mut usize)
226    -> Result<Self, SerializationError> {
227    Self::verify_fixed(stream, pos)?;
228    let mut wf           = RBWaveform::new();
229    wf.event_id          = parse_u32(stream, pos);
230    wf.rb_id             = parse_u8(stream, pos);
231    wf.rb_channel_a      = parse_u8(stream, pos);
232    wf.rb_channel_b      = parse_u8(stream, pos);
233    wf.stop_cell         = parse_u16(stream, pos);
234    wf.paddle_id         = parse_u8 (stream, pos);
235    if stream.len() < *pos+2*NWORDS {
236      return Err(SerializationError::StreamTooShort);
237    }
238    let data_a           = &stream[*pos..*pos+2*NWORDS];
239    //println!("{} data_a stack size", mem::sizeof(data_a));
240
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    *pos +=2;
254    Ok(wf)
255  }
256
257  fn to_bytestream(&self) -> Vec<u8> {
258    let mut stream = Vec::<u8>::new();
259    stream.extend_from_slice(&Self::HEAD.to_le_bytes());
260    stream.extend_from_slice(&self.event_id.to_le_bytes());
261    stream.extend_from_slice(&self.rb_id.to_le_bytes());
262    stream.extend_from_slice(&self.rb_channel_a.to_le_bytes());
263    stream.extend_from_slice(&self.rb_channel_b.to_le_bytes());
264    stream.extend_from_slice(&self.stop_cell.to_le_bytes());
265    stream.push(self.paddle_id);
266    if self.adc_a.len() != 0 {
267      for k in 0..NWORDS {
268        stream.extend_from_slice(&self.adc_a[k].to_le_bytes());  
269      }
270    }
271    if self.adc_b.len() != 0 {
272      for k in 0..NWORDS {
273        stream.extend_from_slice(&self.adc_b[k].to_le_bytes());  
274      }
275    }
276    stream.extend_from_slice(&Self::TAIL.to_le_bytes());
277    stream
278  }
279}
280
281impl fmt::Display for RBWaveform {
282  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
283    let mut repr = String::from("<RBWaveform:");
284    repr += &(format!("\n  Event ID  : {}", self.event_id));
285    repr += &(format!("\n  RB        : {}", self.rb_id));
286    repr += &(format!("\n  ChannelA  : {}", self.rb_channel_a));
287    repr += &(format!("\n  ChannelB  : {}", self.rb_channel_b));
288    repr += &(format!("\n  Paddle ID : {}", self.paddle_id));
289    repr += &(format!("\n  Stop cell : {}", self.stop_cell));
290    if self.adc_a.len() >= 273 {
291      repr += &(format!("\n  adc [A] [{}]      : .. {} {} {} ..",self.adc_a.len(), self.adc_a[270], self.adc_a[271], self.adc_a[272]));
292    } else {
293      repr += &(String::from("\n  adc [A] [EMPTY]"));
294    }
295    if self.adc_b.len() >= 273 {
296      repr += &(format!("\n  adc [B] [{}]      : .. {} {} {} ..",self.adc_b.len(), self.adc_b[270], self.adc_b[271], self.adc_b[272]));
297    } else {
298      repr += &(String::from("\n  adc [B] [EMPTY]"));
299    }
300    write!(f, "{}", repr)
301  }
302}
303
304//---------------------------------------------------
305
306#[cfg(feature="pybindings")]
307#[pymethods]
308impl RBWaveform {
309  
310  /// Paddle ID of this wveform (1-160)
311  #[getter]
312  fn get_paddle_id(&self) -> u8 {
313    self.paddle_id
314  }
315
316  #[getter]
317  fn max_peak_a_guess(&self) -> f32 {
318    self.guess_max_peak_a()
319  }
320  
321  #[getter]
322  fn max_peak_b_guess(&self) -> f32 {
323    self.guess_max_peak_b()
324  }
325
326  #[getter]
327  fn get_rb_id(&self) -> u8 {
328    self.rb_id
329  }
330  
331  #[getter]
332  fn get_event_id(&self) -> u32 {
333    self.event_id
334  }
335  
336  #[getter]
337  fn get_rb_channel_a(&self) -> u8 {
338    self.rb_channel_a
339  }
340  
341  #[getter]
342  fn get_rb_channel_b(&self) -> u8 {
343    self.rb_channel_b
344  }
345  
346  #[getter]
347  fn get_stop_cell(&self) -> u16 {
348    self.stop_cell
349  }
350 
351  // FIXME - make this consistent 
352  #[getter]
353  fn get_adc_a<'_py>(&self, py: Python<'_py>) ->  Bound<'_py, PyArray1<u16>> {
354    //let arr = PyArray1::<u16>::from_slice(py, self.adc_a.as_slice());
355    self.adc_a.to_pyarray(py)
356    //Ok(arr)
357  }
358  
359  #[getter]
360  fn get_adc_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<u16>>> {
361    let arr = PyArray1::<u16>::from_slice(py, self.adc_b.as_slice());
362    Ok(arr)
363  }
364  
365  #[getter]
366  fn get_voltages_a<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
367    let arr = PyArray1::<f32>::from_slice(py, self.voltages_a.as_slice());
368    Ok(arr)
369  }
370
371  #[getter]
372  fn get_times_a<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
373    let arr    = PyArray1::<f32>::from_slice(py, self.nanoseconds_a.as_slice());
374    Ok(arr)
375  }
376
377  #[getter]
378  fn get_voltages_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
379    let arr = PyArray1::<f32>::from_slice(py, self.voltages_b.as_slice());
380    Ok(arr)
381  }
382
383  /// Time over threshold - waveform needs to be 
384  /// calibrated. 
385  /// Paddle end A
386  ///
387  /// # Arguments:
388  ///   * threshold : value in mV
389  fn get_tot_a(&self, threshold : f32) -> f32 {
390    self.time_over_threshold_a(threshold)
391  }
392  
393  /// Time over threshold - waveform needs to be 
394  /// calibrated. 
395  /// Paddle end B
396  ///
397  /// # Arguments:
398  ///   * threshold : value in mV
399  fn get_tot_b(&self, threshold : f32) -> f32 {
400    self.time_over_threshold_b(threshold)
401  }
402
403  #[pyo3(name="subtract_pedestals")]
404  fn subtract_pedestals_py(&mut self) {
405    self.subtract_pedestals()
406  }
407
408  #[getter]
409  fn get_charge_a_trap(&self) -> f32 {
410    self.charge_a_trap()
411  }
412  
413  #[getter]
414  fn get_charge_a_below_500_trap(&self) -> f32 {
415    self.charge_a_below_500() 
416  } 
417  
418  #[getter]
419  fn get_charge_b_below_500_trap(&self) -> f32 {
420    self.charge_b_below_500() 
421  } 
422
423
424  #[getter]
425  fn get_charge_b_trap(&self) -> f32 {
426    self.charge_b_trap()
427  }
428
429  #[getter]
430  fn get_times_b<'_py>(&self, py: Python<'_py>) ->  PyResult<Bound<'_py, PyArray1<f32>>> {
431    let arr = PyArray1::<f32>::from_slice(py, self.nanoseconds_b.as_slice());
432    Ok(arr)
433  }
434  
435  /// Apply the readoutboard calibration to convert adc/bins
436  /// to millivolts and nanoseconds
437  #[pyo3(name="calibrate")]
438  fn calibrate_py(&mut self, cali : &RBCalibrations) -> PyResult<()> {
439    match self.calibrate(&cali) {
440      Ok(_) => {
441        return Ok(());
442      }
443      Err(err) => {
444        return Err(PyValueError::new_err(err.to_string()));
445      }
446    }
447  }
448  
449  #[pyo3(name="apply_spike_filter")]
450  fn apply_spike_filter_py(&mut self) {
451    self.apply_spike_filter();
452  }
453}
454
455#[cfg(feature="pybindings")]
456pythonize_packable!(RBWaveform);
457
458// needs to fix something up
459//#[cfg(feature="pybindings")]
460//pythonize_telemetry!(RBWaveform);
461
462//---------------------------------------------------
463
464#[cfg(feature = "random")]
465impl FromRandom for RBWaveform {
466    
467  fn from_random() -> Self {
468    let mut wf      = Self::new();
469    let mut rng     = rand::rng();
470    wf.event_id     = rng.random::<u32>();
471    wf.rb_id        = rng.random_range(1..50);
472    wf.rb_channel_a = rng.random_range(1..9);
473    wf.rb_channel_b = rng.random_range(1..9);
474    wf.stop_cell    = rng.random_range(0..1024);
475    //wf.paddle_id    = rng.random::<u8>();
476    let random_numbers_a: Vec<u16> = (0..NWORDS).map(|_| rng.random()).collect();
477    wf.adc_a        = random_numbers_a;
478    let random_numbers_b: Vec<u16> = (0..NWORDS).map(|_| rng.random()).collect();
479    wf.adc_b        = random_numbers_b;
480    wf
481  }
482}
483
484//---------------------------------------------------
485
486#[test]
487#[cfg(feature = "random")]
488fn pack_rbwaveform() {
489  for _ in 0..100 {
490    let wf   = RBWaveform::from_random();
491    let test : RBWaveform = wf.pack().unpack().unwrap();
492    assert_eq!(wf, test);
493  }
494}
495
496#[test]
497#[cfg(feature="random")]
498fn emit_rbwaveform() {
499  for _ in 0..100 {
500    let ev = RBEvent::from_random();
501    for wf in ev.get_rbwaveforms() {
502      assert_eq!(ev.header.rb_id, wf.rb_id);
503    }
504  }
505}
506
507
gondola_core/events/rb_waveform.rs

gondola_core/events/
rb_waveform.rs
