gondola_core/calibration/
tof.rs

1//! Calibration routines for the GAPS TOF system
2//!
3// This file is part of gaps-online-software and published 
4// under the GPLv3 license
5
6use crate::prelude::*;
7
8/// Roll over the entries from the end of a vector 
9/// to the beginning by a given offset.
10/// 
11/// This is similar to 
12/// https://numpy.org/doc/2.2/reference/generated/numpy.roll.html
13///
14/// # Arguments:
15///   * `vec`   : The vector to be rolled over. It will 
16///               be changed in place 
17///   * `offset`: The signed number to shift elements by (can be to 
18///               the left or right)
19pub fn roll<T: Clone>(vec: &mut Vec<T>, offset: isize) {
20  let len = vec.len() as isize;
21  if len <= 1 {
22      return;
23  }
24  let offset = offset % len;
25  if offset == 0 {
26      return;
27  }
28  let split_point = if offset > 0 {
29      len - offset
30  } else {
31      -offset
32  } as usize;
33
34  let mut temp = Vec::with_capacity(len as usize);
35  temp.extend_from_slice(&vec[split_point..]);
36  temp.extend_from_slice(&vec[..split_point]);
37
38  vec.clear();
39  vec.extend_from_slice(&temp);
40}
41
42//-----------------------------------------------
43
44/// Simplified version of spike cleaning 
45///
46/// Taken over from Jamie's python code
47pub fn clean_spikes(trace : &mut Vec<f32>, vcaldone : bool) {
48  //# TODO: make robust (symmetric, doubles, fixed/estimated spike height)
49  let mut thresh : f32 = 360.0;
50  if vcaldone {
51    thresh = 16.0;
52  }
53
54  let mut spf_allch = vec![0usize;1023];
55  let mut spf_sum   = vec![0f32;1024];
56  let tracelen      = trace.len();
57  let spikefilter0 = &trace[0..tracelen-3];
58  let spikefilter1 = &trace[1..tracelen-2];
59  let spikefilter2 = &trace[2..tracelen-1];
60  let spikefilter3 = &trace[3..tracelen];
61  let spf_len      = spikefilter0.len();
62  for k in 0..spf_len {
63    spf_sum[k] += spikefilter1[k] - spikefilter0[k] + spikefilter2[k] - spikefilter3[k];
64  }
65  for k in 0..spf_len {
66    if spf_sum[k] > thresh {
67      spf_allch[k] += 1;
68    }
69  }
70  let mut spikes = Vec::<usize>::new();
71  for k in 0..spf_allch.len() {
72    if spf_allch[k] >= 2 {
73      spikes.push(k);
74    }
75  }
76  for spike in spikes.iter() {
77    let d_v : f32 = (trace[spike+3] - trace[*spike])/3.0;
78    trace[spike+1] = trace[*spike] + d_v;
79    trace[spike+2] = trace[*spike] + 2.0*d_v;
80  }
81}
82
83//-----------------------------------------------
84
85//#[cfg(feature="pybindings")]
86//#[pyfunction]
87//#[pyo3(name="clean_spikes")]
88//pub fn clean_spikes_pyx<'_py>(value : Bound<'_py,PyArray1<f32>>, vcal_done : bool) {
89//  unsafe {
90//    match clean_spikes(value.as_slice().unwrap(), vcal_done) {
91//      Err(err) => {
92//        return Err(PyValueError::new_err(err.to_string()));
93//      }
94//      Ok(max_val) => {
95//        return Ok(max_val);
96//      }
97//    }
98//  }
99//}
100
101
102//-----------------------------------------------
103
104#[derive(Debug, Copy, Clone)]
105#[repr(u8)]
106pub enum Edge {
107  Unknown = 0u8,
108  Rising  = 10u8,
109  Falling = 20u8,
110  Average = 30u8,
111  None    = 40u8
112}
113
114impl fmt::Display for Edge {
115  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
116    let r : &str;
117    match self {
118      Edge::Unknown => {r = "Unknown"} 
119      Edge::Rising  => {r = "Rising"},
120      Edge::Falling => {r = "Falling"},
121      Edge::Average => {r = "Average"},
122      Edge::None    => {r = "None"}
123    }
124    write!(f, "<Edge: {}>", r)
125  }
126}
127
128impl From<u8> for Edge {
129  fn from(value: u8) -> Self {
130    match value {
131      0u8  => Edge::Unknown,
132      10u8 => Edge::Rising,
133      20u8 => Edge::Falling,
134      30u8 => Edge::Average,
135      40u8 => Edge::None,
136      _    => Edge::Unknown
137    }
138  }
139}
140
141#[cfg(feature = "random")]
142impl FromRandom for Edge {
143  
144  fn from_random() -> Self {
145    let choices = [
146      Edge::Rising,
147      Edge::Falling,
148      Edge::Average,
149      Edge::None,
150    ];
151    let mut rng  = rand::rng();
152    let idx = rng.random_range(0..choices.len());
153    choices[idx]
154  }
155}
156
157
158/// smaller packet to bring through the gcu
159#[derive(Debug, Clone, PartialEq)]
160#[cfg_attr(feature="pybindings", pyclass)]
161pub struct RBCalibrationFlightT {
162  pub rb_id      : u8,
163  pub timestamp  : u32, // time the calibration was taken
164  pub tbins      : [[f16;NWORDS];NCHN], // cell width (ns)
165}
166
167//-----------------------------------------------
168
169/// Get the periods of a sine wave
170///
171/// # Arguments:
172///   * dts : some version of tcal data? FIXME
173pub fn get_periods(trace   : &Vec<f32>,
174                   dts     : &Vec<f32>,
175                   nperiod : f32,
176                   nskip   : f32,
177                   edge    : &Edge) -> (Vec<usize>, Vec<f32>) {
178  let mut trace_c = trace.clone();
179  let mut periods = Vec::<f32>::new();
180  if trace_c.len() == 0 {
181    let foo = Vec::<usize>::new();
182    return (foo, periods);
183  }
184  let firstbin : usize = 20;
185  let lastbin = firstbin + (nperiod * (900.0/nperiod).floor()).floor() as usize;
186  //info!("firstbin {} lastbin {}", firstbin, lastbin);
187  let mut vec_for_median = Vec::<f32>::new();
188  for bin in firstbin..lastbin {
189    if trace[bin].is_nan() {
190      continue;
191    }
192    vec_for_median.push(trace[bin]);
193  }
194  let median_val = median(&vec_for_median);
195  debug!("median val {median_val}");
196  for k in 0..trace_c.len() {
197    trace_c[k] -= median_val;
198  }
199  let mut zcs = find_zero_crossings(&trace_c);
200  //trace!("Found {} zero crossings!", zcs.len());
201  let mut zcs_nskip = Vec::<usize>::with_capacity(zcs.len());
202  for k in 0..zcs.len() {
203    if zcs[k] > nskip as usize {
204      zcs_nskip.push(zcs[k]);
205    }
206  }
207  zcs = zcs_nskip;
208  let mut zcs_filter = Vec::<usize>::with_capacity(zcs.len());
209  for k in 0..zcs.len() {
210    match edge {
211      Edge::Rising => {
212        if trace_c[zcs[k]] < 0.0 {
213          zcs_filter.push(zcs[k]);
214        }
215      }
216      Edge::Falling => {
217        // What about the equal case?
218        if trace_c[zcs[k]] > 0.0 {
219          zcs_filter.push(zcs[k]);
220        }
221      },
222      Edge::None => {
223        warn!("Unsure what to do for Edge::None");
224      },
225      Edge::Average => {
226        warn!("Unsure what to do for Edge::Average");
227      },
228      Edge::Unknown => {
229        warn!("Unsure what to do for Edge::Unknown");
230      }
231    }
232  }
233  debug!("Found {} zero crossings!", zcs_filter.len()); 
234  zcs = zcs_filter;
235  if zcs.len() < 3 {
236    return (zcs, periods);
237  }
238 
239  for k in 0..zcs.len() -1 {
240    let zcs_a  = &zcs[k];
241    let zcs_b  = &zcs[k+1];
242    // FIXME - there is an issue with the last
243    // zero crossings
244    if zcs_a + 2 > trace_c.len() || zcs_b + 2 > trace_c.len() {
245      continue;
246    }
247    let tr_a   = &trace_c[*zcs_a..*zcs_a+2];
248    let tr_b   = &trace_c[*zcs_b..*zcs_b+2];
249    let mut period : f32 = 0.0;
250    for n in zcs_a+1..*zcs_b {
251      period += dts[n];
252    }
253    period += dts[*zcs_a]*f32::abs(tr_a[1]/(tr_a[1] - tr_a[0])); // first semi bin
254    period += dts[*zcs_b]*f32::abs(tr_b[0]/(tr_b[1] - tr_b[0])); // first semi bin
255    if period.is_nan() {
256      error!("NAN in period found!");
257      continue;
258    }
259    
260    if f32::abs(*zcs_b as f32 - *zcs_a as f32 - nperiod) > 5.0 {
261      let mut zcs_tmp = Vec::<usize>::new();
262      zcs_tmp.extend_from_slice(&zcs[0..k+1]);
263      zcs = zcs_tmp;
264      break;
265    }
266    trace!("zcs_a, zcs_b, period, nperiod {} {} {} {}", zcs_a, zcs_b, period, nperiod);
267    periods.push(period);
268  }
269  debug!("Calculated {} zero-crossings and {} periods!", zcs.len(), periods.len());
270  (zcs, periods)
271}
272
273//-----------------------------------------------
274
275/// Designed to match np.where(np.diff(np.signbit(trace)))\[0\] 
276/// FIXME -> This needs to be moved to analysis!
277pub fn find_zero_crossings(trace: &Vec<f32>) -> Vec<usize> {
278  let mut zero_crossings = Vec::new();
279  for i in 1..trace.len() {
280    // acccount for the fact that sometimes the second/first point can be 0
281    if (trace[i - 1] > 0.0 && trace[i] < 0.0) || (trace[i - 1] < 0.0 && trace[i] > 0.0) {
282      zero_crossings.push(i - 1);
283    }
284    if i < trace.len() - 1 {
285      if trace[i - 1] > 0.0 && trace[i] == 0.0 && trace[i+1] < 0.0 {
286        zero_crossings.push(1);
287      }
288      if trace[i - 1] < 0.0 && trace[i] == 0.0 && trace[i+1] > 0.0 {
289        zero_crossings.push(i);
290      }
291    }
292  }
293  zero_crossings
294}
295
296//-----------------------------------------------
297
298impl RBCalibrationFlightT {
299  pub fn new() -> Self {
300    Self {
301      rb_id     : 0,
302      timestamp : 0,
303      tbins     : [[f16::from_f32(0.0);NWORDS];NCHN],
304    }
305  }
306}
307
308impl Serialization for RBCalibrationFlightT {
309  const SIZE            : usize = NCHN*NWORDS*2 + 4 + 1; 
310  const HEAD            : u16   = 0xAAAA; // 43690 
311  const TAIL            : u16   = 0x5555; // 21845 
312  
313  fn from_bytestream(bytestream : &Vec<u8>, 
314                     pos        : &mut usize)
315    -> Result<Self, SerializationError> { 
316    let mut rb_cal = Self::new();
317    if parse_u16(bytestream, pos) != Self::HEAD {
318      return Err(SerializationError::HeadInvalid {});
319    }
320    rb_cal.rb_id                = parse_u8(bytestream, pos);
321    rb_cal.timestamp            = parse_u32(bytestream, pos);
322    for ch in 0..NCHN {
323      for k in 0..NWORDS {
324        let value         = parse_f16(bytestream, pos);
325        rb_cal.tbins[ch][k]      = value;
326      }
327    }
328    *pos += 2;
329    Ok(rb_cal)
330  }
331
332  fn to_bytestream(&self) -> Vec<u8> {
333    
334    let mut bs = Vec::<u8>::with_capacity(Self::SIZE);
335    bs.extend_from_slice(&Self::HEAD.to_le_bytes());
336    bs.extend_from_slice(&self.rb_id.to_le_bytes());
337    bs.extend_from_slice(&self.timestamp.to_le_bytes());
338    for ch in 0..NCHN {
339      for k in 0..NWORDS {
340        bs.extend_from_slice(&self.tbins[ch][k]     .to_le_bytes());
341      }
342    }
343    bs.extend_from_slice(&Self::TAIL.to_le_bytes());
344    bs
345  }
346}
347
348impl TofPackable for RBCalibrationFlightT {
349  const TOF_PACKET_TYPE : TofPacketType = TofPacketType::RBCalibrationFlightT;
350}
351
352impl fmt::Display for RBCalibrationFlightT {
353  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
354    let mut timestamp_str = String::from("?");
355    match Utc.timestamp_opt(self.timestamp.into(), 0) {
356      LocalResult::Single(datetime_utc) => {
357        timestamp_str = datetime_utc.format("%Y/%m/%d %H:%M:%S").to_string();
358      },
359      LocalResult::Ambiguous(_, _) => {
360        error!("The given timestamp is ambiguous.");
361      },
362      LocalResult::None => {
363        error!("The given timestamp is not valid.");
364      },
365    }
366
367    //let datetime_utc: DateTime<Utc> = Utc.timestamp_opt(self.timestamp as i64, 0).earliest().unwrap_or(DateTime::<Utc>::from_timestamp_millis(0).unwrap());
368    write!(f, 
369  "<RBCalibrationFlightT [{} UTC]:
370      RB             : {}
371      T Bins    [ch0]: .. {:?} {:?} ..>",
372      timestamp_str,
373      self.rb_id,
374      self.tbins[0][98],
375      self.tbins[0][99]
376    )
377  } 
378}
379
380#[cfg(feature = "random")]
381impl FromRandom for RBCalibrationFlightT {
382    
383  fn from_random() -> Self {
384    let mut cali   = Self::new();
385    let mut rng    = rand::rng();
386    cali.rb_id     = rng.random::<u8>();
387    cali.timestamp = rng.random::<u32>();
388    for ch in 0..NCHN {
389      for n in 0..NWORDS { 
390        cali.tbins     [ch][n] = f16::from_f32(rng.random::<f32>());
391      }
392    }
393    cali
394  }
395}
396
397//-----------------------------------------------
398
399#[cfg(feature="pybindings")]
400#[pymethods]
401impl RBCalibrationFlightT {
402  #[getter]
403  fn get_rb_id(&self) -> u8 {
404    self.rb_id
405  }
406  
407  #[getter]
408  fn get_timestamp(&self) -> u32 {
409    self.timestamp
410  }
411
412  fn get_tbins<'_py>(&self, py: Python<'_py>, channel : u8 ) -> PyResult<Bound<'_py, PyArray1<f32>>> {
413    if channel < 1 || channel > 9 {
414      return Err(PyValueError::new_err("Channel must be > 0 and < 9"));
415    }
416    let mut py_tbins = [0.0f32;NWORDS];
417    for k in 0..NWORDS {
418      py_tbins[k] = self.tbins[channel as usize][k].to_f32();
419    };
420    let py_array = py_tbins.to_pyarray(py);
421    Ok(py_array)
422  }
423}
424
425#[cfg(feature="pybindings")]
426pythonize_packable!(RBCalibrationFlightT);
427
428//-----------------------------------------------
429
430#[derive(Debug, Clone, PartialEq)]
431#[cfg_attr(feature="pybindings", pyclass)]
432pub struct RBCalibrationFlightV {
433  pub rb_id      : u8,
434  pub d_v        : f32, // input voltage difference between 
435                        // vcal_data and noi data
436  pub timestamp  : u32, // time the calibration was taken
437  pub v_offsets : [[f16;NWORDS];NCHN], // voltage offset
438  pub v_dips    : [[f16;NWORDS];NCHN], // voltage "dip" (time-dependent correction)
439  pub v_inc     : [[f16;NWORDS];NCHN], // voltage increment (mV/ADC unit)
440}
441
442impl RBCalibrationFlightV {
443  pub fn new() -> Self {
444    Self {
445      rb_id     : 0,
446      d_v       : 0.9,
447      timestamp : 0,
448      v_offsets : [[f16::from_f32(0.0);NWORDS];NCHN],
449      v_dips    : [[f16::from_f32(0.0);NWORDS];NCHN],
450      v_inc     : [[f16::from_f32(0.0);NWORDS];NCHN],
451    }
452  }
453}
454
455impl Serialization for RBCalibrationFlightV {
456  const SIZE            : usize = NCHN*NWORDS*2*3 + 4 + 4 + 1; 
457  const HEAD            : u16   = 0xAAAA; // 43690 
458  const TAIL            : u16   = 0x5555; // 21845 
459
460  fn from_bytestream(bytestream : &Vec<u8>, 
461                     pos        : &mut usize)
462    -> Result<Self, SerializationError> { 
463    let mut rb_cal = Self::new();
464    if parse_u16(bytestream, pos) != Self::HEAD {
465      return Err(SerializationError::HeadInvalid {});
466    }
467    rb_cal.rb_id                = parse_u8(bytestream, pos);
468    rb_cal.d_v                  = parse_f32(bytestream, pos);
469    rb_cal.timestamp            = parse_u32(bytestream, pos);
470    for ch in 0..NCHN {
471      for k in 0..NWORDS {
472        let mut value = parse_f16(bytestream, pos);
473        rb_cal.v_offsets[ch][k] = value;
474        value         = parse_f16(bytestream, pos);
475        rb_cal.v_dips[ch][k]    = value;
476        value         = parse_f16(bytestream, pos);
477        rb_cal.v_inc[ch][k]     = value;
478      }
479    }
480    *pos += 2;
481    Ok(rb_cal)
482  }
483  
484  fn to_bytestream(&self) -> Vec<u8> {
485    let mut bs = Vec::<u8>::with_capacity(Self::SIZE);
486    bs.extend_from_slice(&Self::HEAD.to_le_bytes());
487    bs.extend_from_slice(&self.rb_id.to_le_bytes());
488    bs.extend_from_slice(&self.d_v.to_le_bytes());
489    bs.extend_from_slice(&self.timestamp.to_le_bytes());
490    for ch in 0..NCHN {
491      for k in 0..NWORDS {
492        bs.extend_from_slice(&self.v_offsets[ch][k].to_le_bytes());
493        bs.extend_from_slice(&self.v_dips[ch][k]   .to_le_bytes());
494        bs.extend_from_slice(&self.v_inc[ch][k]    .to_le_bytes());
495      }
496    }
497    bs.extend_from_slice(&Self::TAIL.to_le_bytes());
498    bs
499  }
500}
501
502impl TofPackable for RBCalibrationFlightV {
503  const TOF_PACKET_TYPE : TofPacketType = TofPacketType::RBCalibrationFlightV;
504}
505
506impl fmt::Display for RBCalibrationFlightV {
507  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
508    let mut timestamp_str = String::from("?");
509    match Utc.timestamp_opt(self.timestamp.into(), 0) {
510      LocalResult::Single(datetime_utc) => {
511        timestamp_str = datetime_utc.format("%Y/%m/%d %H:%M:%S").to_string();
512      },
513      LocalResult::Ambiguous(_, _) => {
514        error!("The given timestamp is ambiguous.");
515      },
516      LocalResult::None => {
517        error!("The given timestamp is not valid.");
518      },
519    }
520    write!(f, 
521  "<RBCalibrationFlightV [{} UTC]:
522      RB             : {}
523      V Offsets [ch0]: .. {:?} {:?} ..
524      V Incrmts [ch0]: .. {:?} {:?} ..
525      V Dips    [ch0]: .. {:?} {:?} ..>",
526      timestamp_str,
527      self.rb_id,
528      self.v_offsets[0][98],
529      self.v_offsets[0][99],
530      self.v_inc[0][98],
531      self.v_inc[0][99],
532      self.v_dips[0][98],
533      self.v_dips[0][99]
534    )
535  } 
536}
537
538#[cfg(feature = "random")]
539impl FromRandom for RBCalibrationFlightV {
540    
541  fn from_random() -> Self {
542    let mut cali     = Self::new();
543    let mut rng      = rand::rng();
544    cali.rb_id       = rng.random::<u8>();
545    cali.d_v         = rng.random::<f32>();
546    cali.timestamp   = rng.random::<u32>();
547    for ch in 0..NCHN {
548      for n in 0..NWORDS { 
549        cali.v_offsets  [ch][n] = f16::from_f32(rng.random::<f32>());
550        cali.v_dips     [ch][n] = f16::from_f32(rng.random::<f32>());
551        cali.v_inc      [ch][n] = f16::from_f32(rng.random::<f32>());
552      }
553    }
554    cali
555  }
556}
557
558//-----------------------------------------------
559
560#[cfg(feature="pybindings")]
561#[pymethods]
562impl RBCalibrationFlightV {
563  #[getter]
564  fn get_rb_id(&self) -> u8 {
565    self.rb_id
566  }
567 
568  #[getter]
569  fn get_d_v(&self) -> f32 {
570    self.d_v
571  }
572
573  #[getter]
574  fn get_timestamp(&self) -> u32 {
575    self.timestamp
576  }
577
578  fn get_v_offsets<'_py>(&self, py: Python<'_py>, channel : u8 ) -> PyResult<Bound<'_py, PyArray1<f32>>> {
579    if channel < 1 || channel > 9 {
580      return Err(PyValueError::new_err("Channel must be > 0 and < 9"));
581    }
582    let mut py_voffsets = [0.0f32;NWORDS];
583    for k in 0..NWORDS {
584      py_voffsets[k] = self.v_offsets[channel as usize][k].to_f32();
585    };
586    let py_array = py_voffsets.to_pyarray(py);
587    Ok(py_array)
588  }
589
590  fn get_v_dips<'_py>(&self, py: Python<'_py>, channel : u8 ) -> PyResult<Bound<'_py, PyArray1<f32>>> {
591    if channel < 1 || channel > 9 {
592      return Err(PyValueError::new_err("Channel must be > 0 and < 9"));
593    }
594    let mut py_vdips = [0.0f32;NWORDS];
595    for k in 0..NWORDS {
596      py_vdips[k] = self.v_dips[channel as usize][k].to_f32();
597    };
598    let py_array = py_vdips.to_pyarray(py);
599    Ok(py_array)
600  }
601
602  fn get_v_inc<'_py>(&self, py: Python<'_py>, channel : u8 ) -> PyResult<Bound<'_py, PyArray1<f32>>> {
603    if channel < 1 || channel > 9 {
604      return Err(PyValueError::new_err("Channel must be > 0 and < 9"));
605    }
606    let mut py_vinc = [0.0f32;NWORDS];
607    for k in 0..NWORDS {
608      py_vinc[k] = self.v_inc[channel as usize][k].to_f32();
609    };
610    let py_array = py_vinc.to_pyarray(py);
611    Ok(py_array)
612  }
613}
614
615#[cfg(feature="pybindings")]
616pythonize_packable!(RBCalibrationFlightV);
617
618//-----------------------------------------------
619
620// new way of property testing using quicktest!
621#[test]
622fn prop_roll_then_unroll_gives_original() {
623  fn prop(mut vec: Vec<u8>, offset: i8) -> bool {
624    let original = vec.clone();
625    let offset = offset as isize;
626
627    roll(&mut vec, offset);
628    roll(&mut vec, -offset);
629
630    vec == original
631  }
632  quickcheck::QuickCheck::new().tests(100).quickcheck(prop as fn(Vec<u8>, i8) -> bool);
633}
634
635//-----------------------------------------------
636
637#[derive(Debug, Clone, PartialEq)]
638#[cfg_attr(feature="pybindings", pyclass)]
639pub struct RBCalibrations {
640  pub rb_id      : u8,
641  pub d_v        : f32, // input voltage difference between 
642                        // vcal_data and noi data
643  pub timestamp  : u32, // time the calibration was taken
644  pub serialize_event_data : bool,
645  // calibration constants
646  pub v_offsets : [[f32;NWORDS];NCHN], // voltage offset
647  pub v_dips    : [[f32;NWORDS];NCHN], // voltage "dip" (time-dependent correction)
648  pub v_inc     : [[f32;NWORDS];NCHN], // voltage increment (mV/ADC unit)
649  pub tbin      : [[f32;NWORDS];NCHN], // cell width (ns)
650
651  // calibration data
652  pub vcal_data : Vec::<RBEvent>,
653  pub tcal_data : Vec::<RBEvent>,
654  pub noi_data  : Vec::<RBEvent>,
655}
656
657impl RBCalibrations {
658  // skip the first n cells for the 
659  // voltage calibration. Historically,
660  // this had been 2.
661  pub const NSKIP       : usize = 2;
662  pub const SINMAX      : usize = 60; // ~1000 ADC units
663  pub const DVCUT       : f32   = 15.0; // ns away that should be considered
664  pub const NOMINALFREQ : f32   = 2.0; // nominal sampling frequency,
665                                       // GHz
666  pub const CALFREQ     : f32   = 0.025; // calibration sine wave frequency (25MHz)
667 
668  /// Re-assemble a RBCalibration from chopped up parts
669  pub fn assemble_from_flightcal(fcal_t : RBCalibrationFlightT,
670                                 fcal_v : RBCalibrationFlightV)
671    -> Result<Self, CalibrationError> {
672    if (fcal_t.timestamp != fcal_v.timestamp) || (fcal_t.rb_id != fcal_v.rb_id) {
673      error!("These calibrations do not match! {} , {}", fcal_t, fcal_v);
674      return Err(CalibrationError::IncompatibleFlightCalibrations);
675    }
676    let mut cal              = Self::new(fcal_t.rb_id);
677    cal.rb_id                = fcal_t.rb_id;
678    cal.timestamp            = fcal_t.timestamp;
679    cal.d_v                  = fcal_v.d_v;
680    cal.serialize_event_data = false;
681    for ch in 0..NCHN {
682      for k in 0..NWORDS {
683        cal.tbin     [ch][k] = f16::to_f32(fcal_t.tbins[ch][k]);
684        cal.v_offsets[ch][k] = f16::to_f32(fcal_v.v_offsets[ch][k]); 
685        cal.v_dips   [ch][k] = f16::to_f32(fcal_v.v_dips[ch][k]);
686        cal.v_inc    [ch][k] = f16::to_f32(fcal_v.v_inc[ch][k]); 
687      }
688    }
689    Ok(cal)
690  }
691
692  /// Return the timing part of the calibration in a 
693  /// package digestable by the flight computer.
694  ///
695  /// Additonal compression by using f16
696  pub fn emit_flighttcal(&self) -> RBCalibrationFlightT {
697    let mut cal = RBCalibrationFlightT::new();
698    cal.rb_id = self.rb_id;
699    cal.timestamp = self.timestamp;
700    for ch in 0..NCHN {
701      for n in 0..NWORDS { 
702        cal.tbins  [ch][n] = f16::from_f32(self.tbin[ch][n]);
703      }
704    }
705    cal
706  }
707  
708  /// Return the voltage part of the calibration in a 
709  /// package digestable by the flight computer.
710  ///
711  /// Additional compression by using f16
712  pub fn emit_flightvcal(&self) -> RBCalibrationFlightV {
713    let mut cal = RBCalibrationFlightV::new();
714    cal.rb_id = self.rb_id;
715    cal.timestamp = self.timestamp;
716    cal.d_v   = self.d_v;
717    for ch in 0..NCHN {
718      for n in 0..NWORDS { 
719        cal.v_offsets  [ch][n] = f16::from_f32(self.v_offsets[ch][n]);
720        cal.v_dips     [ch][n] = f16::from_f32(self.v_dips[ch][n]);
721        cal.v_inc      [ch][n] = f16::from_f32(self.v_inc[ch][n]);
722      }
723    }
724    cal
725  }
726
727  /// Remove events with invalid traces or event fragment bits set
728  pub fn clean_input_data(&mut self) {
729    self.vcal_data.retain(|x|  !x.header.drs_lost_trigger()
730        && !x.header.is_event_fragment()
731        && x.trace_check()); 
732    self.tcal_data.retain(|x|  !x.header.drs_lost_trigger()
733        && !x.header.is_event_fragment()
734        && x.trace_check()); 
735    self.noi_data.retain(|x|   !x.header.drs_lost_trigger() 
736        && !x.header.is_event_fragment()
737        && x.trace_check()); 
738    self.noi_data.sort_by(|a, b| a.header.event_id.cmp(&b.header.event_id));
739    self.vcal_data.sort_by(|a, b| a.header.event_id.cmp(&b.header.event_id));
740    self.tcal_data.sort_by(|a, b| a.header.event_id.cmp(&b.header.event_id));
741  }
742
743  // apply the vcal to a dataset of the calibration
744  // (e.g. timing calibration)
745  fn apply_vcal(&self, 
746                data      : &Vec<RBEvent>)
747      -> (Vec<Vec<Vec<f32>>>,Vec<isize>) {
748    let nevents          = data.len();
749    let mut traces       = Vec::<Vec::<Vec::<f32>>>::new();
750    let mut trace        = Vec::<f32>::with_capacity(NWORDS);
751    let mut stop_cells   = Vec::<isize>::new();
752    let mut empty_events = Vec::<Vec::<f32>>::new();
753    for _ in 0..nevents {
754        empty_events.push(trace.clone());
755    }
756    for ch in 0..NCHN {
757      traces.push(empty_events.clone());
758      for (k,ev) in data.iter().enumerate() {
759        trace.clear();
760        stop_cells.push(ev.header.stop_cell as isize);
761        for k in 0..NWORDS {
762          trace.push(ev.adc[ch][k] as f32);
763        }
764        self.voltages(ch + 1, ev.header.stop_cell as usize,
765                      &ev.adc[ch], &mut trace);
766        traces[ch][k] = trace.clone();
767      }
768    }
769    (traces, stop_cells)
770  }
771
772  // channel is from 0-8
773  pub fn apply_vcal_constants(&self,
774                              adc       : &Vec<f32>,
775                              channel   : usize,  
776                              stop_cell : usize) -> Vec<f32> {
777    let mut waveform = Vec::<f32>::with_capacity(adc.len());
778    let mut value : f32;
779    for k in 0..adc.len() {
780      value  = adc[k] as f32;
781      value -= self.v_offsets[channel][(k + (stop_cell)) %NWORDS];
782      value -= self.v_dips   [channel][k];
783      value *= self.v_inc    [channel][(k + (stop_cell)) %NWORDS];
784      waveform.push(value);
785    } 
786    waveform
787  }
788
789  /// Calculate the offset and dips calibration constants 
790  /// for input data. 
791  ///
792  /// # Return:
793  ///
794  /// offsets, dips
795  fn voltage_offset_and_dips(input_vcal_data : &Vec<RBEvent>) 
796  -> Result<(Vec<Vec<f32>>, Vec<Vec<f32>>), CalibrationError> {
797    if input_vcal_data.len() == 0 {
798      return Err(CalibrationError::EmptyInputData);
799    }
800    let mut all_v_offsets = Vec::<Vec::<f32>>::new();
801    let mut all_v_dips    = Vec::<Vec::<f32>>::new();
802    let nchn = input_vcal_data[0].header.get_nchan();
803
804    debug!("Found {nchn} channels!");
805    for _ in 0..nchn {
806      let empty_vec_off : Vec<f32> = vec![0.0;NWORDS];
807      all_v_offsets.push(empty_vec_off);  
808      let empty_vec_dip : Vec<f32> = vec![0.0;NWORDS];
809      all_v_dips.push(empty_vec_dip);  
810    }
811
812    // we temporarily get the adc traces
813    // traces are [channel][event][adc_cell]
814    let mut traces        = unpack_traces(&input_vcal_data);
815    let mut rolled_traces = traces.clone();
816    for ch in 0..nchn {
817      for n in 0..input_vcal_data.len() {
818        for k in 0..Self::NSKIP {
819          traces[ch][n][k] = f32::NAN;
820          rolled_traces[ch][n][k] = f32::NAN;
821        }
822        // the traces are filled and the first 2 bins
823        // marked with nan, now need to get "rolled over",
824        // so that they start with the stop cell
825        roll(&mut rolled_traces[ch][n],
826             input_vcal_data[n].header.stop_cell as isize); 
827      }// first loop over events done
828      all_v_offsets[ch] = calculate_column_stat(&rolled_traces[ch], mean);
829      //let v_offsets = calculate_column_medians(&rolled_traces[ch]);
830      debug!("We calculated {} voltage offset values for ch {}", all_v_offsets[ch].len(), ch);
831      // fill these in the prepared array structure
832      //for k in 0..v_offsets.len() {
833      //  all_v_offsets[ch][k] = v_offsets[k];
834      //}
835      for (n, ev) in input_vcal_data.iter().enumerate() {
836        // now we roll the v_offsets back
837        let mut v_offsets_rolled = all_v_offsets[ch].clone();
838        roll(&mut v_offsets_rolled, -1*ev.header.stop_cell as isize);
839        for k in 0..traces[ch][n].len() {
840          traces[ch][n][k] -= v_offsets_rolled[k];
841        }
842      }
843      let v_dips = calculate_column_stat(&traces[ch], median);
844      for k in 0..v_dips.len() {
845        if k < Self::NSKIP {
846          all_v_dips[ch][k] = 0.0;
847        } else {
848          all_v_dips[ch][k] = v_dips[k];
849        }
850      }
851    }
852    Ok((all_v_offsets, all_v_dips))
853  }
854
855
856  /// Voltage calibration has to be applied
857  /// 
858  /// # Returns
859  ///
860  ///   vec\[ch\[9\], tbin\[1024\]\]
861  ///
862  pub fn timing_calibration(&self,
863                            edge       : &Edge,
864                            apply_vcal : bool) 
865  -> Result<Vec<Vec<f32>>, CalibrationError> {
866    if self.tcal_data.len() == 0 {
867      error!("Input data for timing calibration is empty!");
868      return Err(CalibrationError::EmptyInputData);
869    }
870    // tcal values are [channel][adc_cell] 
871    let mut all_tcal = Vec::<Vec::<f32>>::new();
872    for _ in 0..NCHN {
873      all_tcal.push(Vec::<f32>::new());
874    }
875    // traces are [channel][event][adc_cell]
876    let mut traces       : Vec<Vec<Vec<f32>>>;
877    let mut stop_cells   = Vec::<isize>::new();
878    if apply_vcal {
879      let result = self.apply_vcal(&self.tcal_data);
880      traces     = result.0;
881      stop_cells = result.1;
882    } else {
883      warn!("Not applying voltage calibration to tcal data. This most likely makes no sense!");
884      traces  = unpack_traces(&self.tcal_data);
885      for ev in self.tcal_data.iter() {
886        stop_cells.push(ev.header.event_id as isize);
887      }
888    }
889    let do_spike_cleaning = true;
890    if do_spike_cleaning {
891      for k_ch in 0..traces.len() {
892        for k_ev in 0..traces[k_ch].len() {
893          clean_spikes(&mut traces[k_ch][k_ev], true);
894        }
895      }
896    }
897    let nwords = traces[0][0].len();
898    for ch in 0..NCHN {
899      for ev in 0..traces[ch].len() {
900        for k in 0..nwords {
901          if k < Self::NSKIP {
902            traces[ch][ev][k]  = f32::NAN;
903          }
904          if f32::abs(traces[ch][ev][k]) > Self::SINMAX as f32 { 
905            traces[ch][ev][k]  = f32::NAN;
906          }// the traces are filled and the first 2 bins
907        }
908      }  
909    }
910
911    let mut rolled_traces = traces.clone();
912    let mut drolled_traces = traces.clone();
913    for ch in 0..NCHN {
914      for ev in 0..traces[ch].len() {
915        roll(&mut rolled_traces[ch][ev],
916             stop_cells[ev]); 
917      }
918    }
919    for ch in 0..NCHN {
920      for ev in 0..traces[ch].len() {
921        for n in 0..traces[ch][ev].len() {
922          let mut dval : f32;
923          if n == traces[ch][ev].len() - 1 {
924            dval = rolled_traces[ch][ev][0] - rolled_traces[ch][ev][traces[ch][ev].len() -1];
925          } else {
926            dval = rolled_traces[ch][ev][n + 1] - rolled_traces[ch][ev][n];      
927          }
928          match edge {
929            Edge::Rising | Edge::Average => {
930              if dval < 0.0 {
931                dval = f32::NAN;
932              }
933            },
934            Edge::Falling => {
935              if dval > 0.0 {
936                dval = f32::NAN;
937              }
938            },
939            _ => {
940              // FIXME - better error handling
941              error!("Only average, rising or falling edge supported!");
942            }
943          } // end match
944          dval = f32::abs(dval); 
945          if f32::abs(dval - 15.0) > Self::DVCUT {
946            dval = f32::NAN;
947          }
948          drolled_traces[ch][ev][n] = dval;
949        } // end loop over adc bins
950      } // end loop over events
951      let col_means = calculate_column_stat(&drolled_traces[ch], mean);
952      let nfreq_vec : Vec<f32> = vec![1.0/Self::NOMINALFREQ;NWORDS];
953      all_tcal[ch]  = nfreq_vec;
954      let ch_mean   = mean(&col_means);
955      for n in 0..all_tcal[ch].len() {
956        all_tcal[ch][n] *= col_means[n]/ch_mean;  
957      }
958    } // end loop over channels
959    Ok(all_tcal)
960  }
961
962  /// Call to the calibration routine, using
963  /// the set input data
964  pub fn calibrate(&mut self) -> Result<(), CalibrationError>{
965    if self.vcal_data.len() == 0
966    || self.tcal_data.len() == 0 
967    || self.noi_data.len() == 0 {
968      return Err(CalibrationError::EmptyInputData);
969    }
970    info!("Starting calculating voltage calibration constants!");
971    let (v_offsets_high, _v_dips_high) 
972        = Self::voltage_offset_and_dips(&self.vcal_data)?;
973    let (v_offsets_low, v_dips_low) 
974        = Self::voltage_offset_and_dips(&self.noi_data)?;
975    // which of the v_offsets do we actually use?
976    for ch in 0..NCHN {
977      for k in 0..NWORDS {
978        self.v_offsets[ch][k] = v_offsets_low[ch][k];
979        self.v_dips[ch][k]    = v_dips_low[ch][k];
980        self.v_inc[ch][k]     = self.d_v/(v_offsets_high[ch][k] - v_offsets_low[ch][k]);
981      }
982    }
983    // at this point, the voltage calibration is complete
984    info!("Filnished calculating voltage calibration constants!");
985    info!("Starting calculating timing calibration constants!");
986    warn!("Calibration only supported for Edge::Average!");
987    // this just suppresses a warning
988    // We have to think if edge will be
989    // a parameter or a constant.
990    let mut edge    = Edge::None;
991    if matches!(edge, Edge::None) {
992      edge = Edge::Average;
993    }
994
995    let mut tcal_av = self.timing_calibration( &edge, true)?;
996    if matches!(edge, Edge::Average) {
997      edge = Edge::Falling;
998      let tcal_falling = self.timing_calibration(&edge, true)?;
999      for ch in 0..NCHN {
1000        for k in 0..tcal_av.len() {
1001          tcal_av[ch][k] += tcal_falling[ch][k];
1002          tcal_av[ch][k] /= 2.0;
1003        }
1004      } 
1005      // for further calibration procedure
1006      edge = Edge::Rising;
1007    } else {
1008      error!("This is not implemented for any other case yet!");
1009      todo!();
1010    }
1011    
1012    // another set of constants
1013    //nevents,nchan,tracelen = gbf.traces.shape
1014    let mut damping   : f32 = 0.1;
1015    let corr_limit    : f32 = 0.05;
1016    //let n_iter_period : f32 = 1000; //#500 or nevents #
1017
1018    let nperiod = Self::NOMINALFREQ/Self::CALFREQ; 
1019    let global = true;
1020    if global {
1021      
1022      //let mut tcal_traces   = Vec::<Vec::<Vec::<f32>>>::new();
1023      //let mut stop_cells    = Vec::<isize>::new();
1024      
1025      let result  = self.apply_vcal(&self.tcal_data);
1026      let tcal_traces = result.0;
1027      let stop_cells  = result.1;
1028
1029      //for n in 0..1000 {
1030      for ch in 0..NCHN {
1031        for ev in 0..tcal_traces[ch].len() {
1032          let tracelen = tcal_traces[ch][ev].len();
1033          let stop_cell = stop_cells[ev];
1034          let mut tcal_av_cp = tcal_av.clone();
1035          roll(&mut tcal_av_cp[ch], -1* (stop_cell as isize));
1036          
1037          let (zcs, periods) = get_periods(&tcal_traces[ch][ev],
1038                                           &tcal_av_cp[ch],
1039                                           nperiod,
1040                                           Self::NSKIP as f32,
1041                                           &edge);
1042          debug!("Will iterate over {} periods!", periods.len());
1043          for (n_p,period) in periods.iter().enumerate() {
1044            if *period == 0.0 {
1045              warn!("period is 0 {:?}", periods);
1046            }
1047            if period.is_nan() {
1048              warn!("period is nan! {:?}", periods);
1049            }
1050            let zcs_a = zcs[n_p]     + stop_cell as usize;
1051            let zcs_b = zcs[n_p + 1] + stop_cell as usize;
1052            let mut corr = (1.0/Self::CALFREQ)/period;
1053            if matches!(edge, Edge::None) {
1054              corr *= 0.5;
1055            }
1056            if f32::abs(corr - 1.0) > corr_limit {
1057              continue;
1058            }
1059            corr = (corr-1.0)*damping + 1.0;
1060            let zcs_a_tl = zcs_a%tracelen;
1061            let zcs_b_tl = zcs_b%tracelen;
1062            if zcs_a < tracelen && zcs_b > tracelen {
1063              for m in zcs_a..tcal_av[ch].len() {
1064                tcal_av[ch][m] *= corr;
1065              }
1066              for m in 0..zcs_b_tl {
1067                tcal_av[ch][m] *= corr;
1068              }
1069            } else {
1070              for m in zcs_a_tl..zcs_b_tl {
1071                tcal_av[ch][m] *= corr;
1072              }
1073            }
1074          }
1075          //n_correct[ch] += 1.0;
1076        } // end over nchannel
1077        damping *= 0.99;
1078      } // end loop over n_iter_period
1079   
1080    } // end global
1081    for ch in 0..NCHN {
1082      for k in 0..NWORDS {
1083        self.tbin[ch][k] = tcal_av[ch][k];
1084      }
1085    }
1086    Ok(())
1087  }
1088
1089  /// Apply the spike cleaning to all channels
1090  pub fn spike_cleaning(voltages  : &mut Vec<Vec<f32>>,
1091                        stop_cell : u16) 
1092    -> Result<(), WaveformError> {
1093
1094    let mut spikes      : [i32;10] = [0;10];
1095    let mut filter      : f32;
1096    let mut dfilter     : f32;
1097    let mut n_neighbor  : usize;
1098    let mut n_rsp       = 0usize;
1099    let mut rsp : [i32;10]    = [-1;10];
1100    // to me, this seems that should be u32
1101    // the 10 is for a maximum of 10 spikes (Jeff)
1102    let mut sp   : [[usize;10];NCHN] = [[0;10];NCHN];
1103    let mut n_sp : [usize;10]        = [0;10];
1104
1105    for j in 0..NWORDS as usize {
1106      for i in 0..NCHN as usize {
1107        filter = -voltages[i][j] + voltages[i][(j + 1) % NWORDS] + voltages[i][(j + 2) % NWORDS] - voltages[i][(j + 3) % NWORDS];
1108        dfilter = filter + 2.0 * voltages[i][(j + 3) % NWORDS] + voltages[i][(j + 4) % NWORDS] - voltages[i][(j + 5) % NWORDS];
1109        if filter > 20.0  && filter < 100.0 {
1110          if n_sp[i] < 10 {   // record maximum of 10 spikes
1111            sp[i][n_sp[i] as usize] = (j + 1) % NWORDS ;
1112            n_sp[i] += 1;
1113          } else {
1114            return Err(WaveformError::TooSpiky);
1115          }            // too many spikes -> something wrong
1116        }// end of if
1117        else if dfilter > 40.0 && dfilter < 100.0 && filter > 10.0 {
1118          if n_sp[i] < 9 {  // record maximum of 10 spikes
1119            sp[i][n_sp[i] as usize] = (j + 1) % NWORDS ;
1120            sp[i][(n_sp[i] + 1) as usize] = (j + 3) % NWORDS ;
1121            n_sp[i] += 2;
1122          } else {
1123            return Err(WaveformError::TooSpiky);
1124          } // too many spikes -> something wrong
1125        } // end of else if
1126
1127      }// end loop over NCHN
1128    } // end loop over NWORDS
1129
1130    // go through all spikes and look for neighbors */
1131    for i in 0..NCHN {
1132      for j in 0..n_sp[i] as usize {
1133        //n_symmetric = 0;
1134        n_neighbor = 0;
1135        for k in 0..NCHN {
1136          for l in 0..n_sp[k] as usize {
1137          //check if this spike has a symmetric partner in any channel
1138            if (sp[i][j] as i32 + sp[k][l] as i32 - 2 * stop_cell as i32) as i32 % NWORDS as i32 == 1022 {
1139              //n_symmetric += 1;
1140              break;
1141            }
1142          }
1143        } // end loop over k
1144        // check if this spike has same spike is in any other channels */
1145        //for (k = 0; k < nChn; k++) {
1146        for k in 0..NCHN {
1147          if i != k {
1148            for l in 0..n_sp[k] {
1149              if sp[i][j] == sp[k][l] {
1150              n_neighbor += 1;
1151              break;
1152              }
1153            } // end loop over l   
1154          } // end if
1155        } // end loop over k
1156
1157        if n_neighbor >= 2 {
1158          for k in 0..n_rsp {
1159            if rsp[k] == sp[i][j] as i32 {break;} // ignore repeats
1160            if n_rsp < 10 && k == n_rsp {
1161              rsp[n_rsp] = sp[i][j] as i32;
1162              n_rsp += 1;
1163            }
1164          }  
1165        }
1166
1167      } // end loop over j
1168    } // end loop over i
1169
1170    // recognize spikes if at least one channel has it */
1171    //for (k = 0; k < n_rsp; k++)
1172    let magic_value : f32 = 14.8;
1173    let mut x       : f32;
1174    let mut y       : f32;
1175
1176    let mut skip_next : bool = false;
1177    for k in 0..n_rsp {
1178      if skip_next {
1179        skip_next = false;
1180        continue;
1181      }
1182      spikes[k] = rsp[k];
1183      //for (i = 0; i < nChn; i++)
1184      for i in 0..NCHN {
1185        if k < n_rsp && i32::abs(rsp[k] as i32 - rsp[k + 1] as i32 % NWORDS as i32) == 2
1186        {
1187          // remove double spike 
1188          let j = if rsp[k] > rsp[k + 1] {rsp[k + 1] as usize}  else {rsp[k] as usize};
1189          x = voltages[i][(j - 1) % NWORDS];
1190          y = voltages[i][(j + 4) % NWORDS];
1191          if f32::abs(x - y) < 15.0 {
1192            voltages[i][j % NWORDS] = x + 1.0 * (y - x) / 5.0;
1193            voltages[i][(j + 1) % NWORDS] = x + 2.0 * (y - x) / 5.0;
1194            voltages[i][(j + 2) % NWORDS] = x + 3.0 * (y - x) / 5.0;
1195            voltages[i][(j + 3) % NWORDS] = x + 4.0 * (y - x) / 5.0;
1196          } else {
1197            voltages[i][j % NWORDS] -= magic_value;
1198            voltages[i][(j + 1) % NWORDS] -= magic_value;
1199            voltages[i][(j + 2) % NWORDS] -= magic_value;
1200            voltages[i][(j + 3) % NWORDS] -= magic_value;
1201          }
1202        } else {
1203          // remove single spike 
1204          x = voltages[i][((rsp[k] - 1) % NWORDS as i32) as usize];
1205          y = voltages[i][(rsp[k] + 2) as usize % NWORDS];
1206          if f32::abs(x - y) < 15.0 {
1207            voltages[i][rsp[k] as usize] = x + 1.0 * (y - x) / 3.0;
1208            voltages[i][(rsp[k] + 1) as usize % NWORDS] = x + 2.0 * (y - x) / 3.0;
1209          } else {
1210            voltages[i][rsp[k] as usize] -= magic_value;
1211            voltages[i][(rsp[k] + 1) as usize % NWORDS] -= magic_value;
1212          }
1213        } // end loop over nchn
1214      } // end loop over n_rsp
1215      if k < n_rsp && i32::abs(rsp[k] - rsp[k + 1] % NWORDS as i32) == 2
1216        {skip_next = true;} // skip second half of double spike
1217    } // end loop over k
1218  Ok(())
1219  }
1220
1221  /// Apply the voltage calibration to a single channel 
1222  /// FIXME - mixing of naming conventions for the channels
1223  ///
1224  /// FIXME - make it return Result<(), CalibrationError>
1225  ///
1226  /// # Arguments
1227  ///
1228  /// * channel   : Channel id 1-9
1229  /// * stop_cell : This channels stop cell 
1230  /// * adc       : Uncalibrated channel data
1231  /// * waveform  : Pre-allocated array to hold 
1232  ///               calibrated waveform data.
1233  pub fn voltages(&self,
1234                  channel   : usize,
1235                  stop_cell : usize,
1236                  adc       : &Vec<u16>,
1237                  waveform  : &mut Vec<f32>) {
1238                  //waveform  : &mut [f32;NWORDS]) {
1239    if channel > 9 || channel == 0 {
1240      error!("There is no channel larger than 9 and no channel 0! Channel {channel} was requested. Can not perform voltage calibration!");
1241      return;
1242    }
1243    if adc.len() != waveform.len() {
1244      error!("Ch{} has {} adc values, however we are expecting {}!", channel,  adc.len(), waveform.len());
1245      return;
1246    }
1247
1248    let mut value : f32; 
1249    for k in 0..NWORDS {
1250      value  = adc[k] as f32;
1251      value -= self.v_offsets[channel -1][(k + (stop_cell)) %NWORDS];
1252      value -= self.v_dips   [channel -1][k];
1253      value *= self.v_inc    [channel -1][(k + (stop_cell)) %NWORDS];
1254      waveform[k] = value;
1255    }
1256  }
1257  
1258  /// Apply the timing calibration to a single channel 
1259  /// 
1260  /// This will allocate the array for the waveform 
1261  /// time bins (unit is ns)
1262  ///
1263  /// # Arguments
1264  ///
1265  /// * channel   : Channel id 1-9
1266  /// * stop_cell : This channels stop cell 
1267  pub fn nanoseconds(&self,
1268                     channel   : usize,
1269                     stop_cell : usize,
1270                     times     : &mut Vec<f32>) {
1271    if channel > 9 || channel == 0 {
1272      error!("There is no channel larger than 9 and no channel 0! Channel {channel} was requested. Can not perform timing calibration!");
1273    }
1274    for k in 1..NWORDS { 
1275      times[k] = times[k-1] + self.tbin[channel -1][(k-1+(stop_cell))%NWORDS];
1276    }
1277  }
1278
1279  pub fn new(rb_id : u8) -> Self {
1280    let timestamp = Utc::now().timestamp() as u32;
1281    Self {
1282      rb_id     : rb_id,
1283      d_v       : 182.0, // FIXME - this needs to be a constant
1284      timestamp : timestamp,
1285      serialize_event_data : false, // per default, don't serialize the data 
1286      v_offsets : [[0.0;NWORDS];NCHN], 
1287      v_dips    : [[0.0;NWORDS];NCHN], 
1288      v_inc     : [[0.0;NWORDS];NCHN], 
1289      tbin      : [[0.0;NWORDS];NCHN],
1290      vcal_data : Vec::<RBEvent>::new(),
1291      tcal_data : Vec::<RBEvent>::new(),
1292      noi_data  : Vec::<RBEvent>::new()
1293    }
1294  }
1295
1296  /// Discard the data to reduce the memory footprint
1297  pub fn discard_data(&mut self) {
1298    self.vcal_data = Vec::<RBEvent>::new();
1299    self.tcal_data = Vec::<RBEvent>::new();
1300    self.noi_data  = Vec::<RBEvent>::new();
1301  }
1302
1303  /// Gets the calibration from a file which 
1304  /// has the RBCalibration stored in a 
1305  /// TofPacket
1306  ///
1307  /// E.g. if it was written with TofPacketWriter
1308  pub fn from_file(filename : String, discard_data : bool) -> Result<Self, SerializationError> {
1309    let mut reader = TofPacketReader::new(&filename);
1310    loop {
1311      match reader.next() {
1312        None => {
1313          error!("Can't load calibration!");
1314          break;
1315        },
1316        Some(pack) => {
1317          if pack.packet_type == TofPacketType::RBCalibration { 
1318            let mut cali = RBCalibrations::from_bytestream(&pack.payload, &mut 0)?;
1319            if discard_data {
1320              cali.discard_data();
1321            }
1322            return Ok(cali);
1323          } else {
1324            continue;
1325          }
1326        }
1327      }
1328    }
1329    Err(SerializationError::StreamTooShort)
1330  }
1331
1332
1333  /// Infer the readoutboard id from the filename
1334  ///
1335  /// Assuming a certain naming scheme for the filename "rbXX_cal.txt"
1336  /// we extract the readoutboard id
1337  pub fn get_id_from_filename(&mut self, filename : &Path) -> u8 {
1338    let rb_id : u8;
1339    match filename.file_name() {
1340      None   => {
1341        error!("Path {} seems non-sensical!", filename.display());
1342        self.rb_id = 0;
1343        return 0;
1344      }
1345      Some(name) => {
1346        // TODO This might panic! Is it ok?
1347        let fname = name.to_os_string().into_string().unwrap();
1348        let id    = &fname[2..4];
1349        // TODO This might panic! Is it ok?
1350        rb_id     = id.parse::<u8>().unwrap();
1351        debug!("Extracted RB ID {} from filename {}", rb_id, fname);
1352      }
1353    }
1354  self.rb_id = rb_id;
1355  rb_id
1356  }
1357  
1358  /// Self check if the timing constants are sane 
1359  pub fn passes_timing_checks(&self) -> bool {
1360    for ch in 0..9 {
1361      let mut mean = 0.0;
1362      for k in 0..NWORDS {
1363        mean += self.tbin[ch][k];
1364      }
1365      mean /= NWORDS as f32;
1366      if mean < 0.48824 || mean > 0.48834 {
1367        error!("Timing calibration for ch {} / RB {} failed. Got mean of {}", ch + 1, self.rb_id, mean);
1368        return false;
1369      }
1370      let tbin_ch = self.tbin[ch].to_vec();
1371      let var = standard_deviation(&tbin_ch).unwrap_or(0.0);
1372      if var < 0.08 || var > 0.15 {
1373        error!("Timing calibration for ch {}/ RB {} failed. Got st dev of {}", ch + 1, self.rb_id, var);
1374        return false;
1375      }
1376    }
1377    debug!("Passed timing calibration sanity checks for RB {}!", self.rb_id);
1378    true
1379  }
1380
1381  /// Self check if the voltage constants are sane
1382  pub fn passes_voltage_checks(&self) -> bool {
1383    for ch in 0..9 {
1384      let mut mean = 0.0;
1385      for k in 0..NWORDS {
1386        mean += self.v_offsets[ch][k];
1387      }
1388      mean /= NWORDS as f32;
1389      if mean < 4200.0 || mean > 5200.0 {
1390        error!("Voltage offset calibration for ch {} / RB {} failed. Got mean of {}", ch + 1, self.rb_id, mean);
1391        return false;
1392      }
1393      let v_off = self.tbin[ch].to_vec();
1394      let var = standard_deviation(&v_off).unwrap_or(0.0);
1395      if var > 150.0 {
1396        error!("Voltage offset calibration for ch {} / RB {} failed. Got st dev of {}", ch + 1, self.rb_id, var);
1397        return false;
1398      }
1399      mean = 0.0;
1400      for k in 0..NWORDS {
1401        mean += self.v_dips[ch][k];
1402      }
1403      mean /= NWORDS as f32;
1404      if mean < -0.5 || mean > 0.5 {
1405        error!("Voltage droop calibration for ch {} / RB {} failed. Got mean of {}", ch + 1, self.rb_id, mean);
1406        return false;
1407      }
1408      mean = 0.0;
1409      for k in 0..NWORDS {
1410        mean += self.v_inc[ch][k];
1411      }
1412      mean /= NWORDS as f32;
1413      if mean < 0.06 || mean > 0.07 {
1414        error!("Voltage gain calibration for ch {} / RB {} failed. Got mean of {}", ch + 1, self.rb_id, mean);
1415        return false;
1416      }
1417      let v_inc = self.v_inc[ch].to_vec();
1418      let var = standard_deviation(&v_inc).unwrap_or(0.0);
1419      if var > 0.00025 {
1420        error!("Voltage gain calibration for ch {} / RB {} failed. Got st dev of {}", ch + 1, self.rb_id, var);
1421        //return false;
1422      }
1423    }
1424    debug!("Passed voltage calibration sanity checks for RB {}!", self.rb_id);
1425    true 
1426  }
1427
1428
1429  /// Check voltage and timing constants for sanity 
1430  ///
1431  /// for each RB channel
1432  /// take distribution of 1024 calibration constants and confirm
1433  /// 
1434  /// Tcal:
1435  ///     average value between 0.48824 and 0.48834
1436  ///     standard deviation between 0.08 and 0.15.
1437  ///     Most common problem: avg of tcal distribution is < 0.48824
1438  /// 
1439  /// Vcal1 (offsets):
1440  ///     average value between 4200 and 5200
1441  ///     standard deviation less than 150
1442  /// 
1443  /// Vcal2 (droop):
1444  ///     average value between -0.5 and 0.5
1445  ///     print out a warning if the standard deviation for any channel is > 5.0
1446  ///     Note that ch1 on any particular board has a more dramatic droop so if ch1 droop is the only channel on the board with this warning, it can be ignored
1447  /// 
1448  /// Vcal3 (gain):
1449  ///     average value between 0.06 and 0.07
1450  ///     standard deviation less than 0.00025
1451  pub fn check(&self) -> bool {
1452    self.passes_timing_checks() && self.passes_voltage_checks() 
1453  }
1454}
1455
1456impl TofPackable for RBCalibrations {
1457  const TOF_PACKET_TYPE : TofPacketType = TofPacketType::RBCalibration;
1458}
1459
1460impl Serialization for RBCalibrations {
1461  const SIZE            : usize = NCHN*NWORDS*4*8 + 4 + 1 + 1; 
1462  const HEAD            : u16   = 0xAAAA; // 43690 
1463  const TAIL            : u16   = 0x5555; // 21845 
1464  
1465  fn from_bytestream(bytestream : &Vec<u8>, 
1466                     pos        : &mut usize)
1467    -> Result<Self, SerializationError> { 
1468    let mut rb_cal = Self::new(0);
1469    if parse_u16(bytestream, pos) != Self::HEAD {
1470      return Err(SerializationError::HeadInvalid {});
1471    }
1472    rb_cal.rb_id                = parse_u8(bytestream, pos);
1473    rb_cal.d_v                  = parse_f32(bytestream, pos);
1474    rb_cal.timestamp            = parse_u32(bytestream, pos);
1475    rb_cal.serialize_event_data = parse_bool(bytestream, pos);
1476    for ch in 0..NCHN {
1477      for k in 0..NWORDS {
1478        let mut value = parse_f32(bytestream, pos);
1479        rb_cal.v_offsets[ch][k] = value;
1480        value         = parse_f32(bytestream, pos);
1481        rb_cal.v_dips[ch][k]    = value;
1482        value         = parse_f32(bytestream, pos);
1483        rb_cal.v_inc[ch][k]     = value;
1484        value         = parse_f32(bytestream, pos);
1485        rb_cal.tbin[ch][k]      = value;
1486      }
1487    }
1488    if rb_cal.serialize_event_data {
1489      let broken_event = RBEvent::new();
1490      let n_noi  = parse_u16(bytestream, pos);
1491      debug!("Found {n_noi} no input data events!");
1492      for _ in 0..n_noi {
1493        match RBEvent::from_bytestream(bytestream, pos) {
1494          Ok(ev) => {
1495            rb_cal.noi_data.push(ev);            
1496          }
1497          Err(err) => {
1498            error!("Unable to read RBEvent! {err}");
1499          }
1500        }
1501        // FIXME - broken event won't advance the pos marker
1502      }
1503      let n_vcal = parse_u16(bytestream, pos); 
1504      debug!("Found {n_vcal} VCal data events!");
1505      for _ in 0..n_vcal {
1506        match RBEvent::from_bytestream(bytestream, pos) {
1507          Err(err) => {
1508            error!("Found broken event {err}");
1509          },
1510          Ok(good_ev) => {
1511            rb_cal.vcal_data.push(good_ev);
1512          }
1513        }
1514      }
1515      let n_tcal = parse_u16(bytestream, pos); 
1516      debug!("Found {n_tcal} TCal data events!");
1517      for _ in 0..n_tcal {
1518        rb_cal.tcal_data.push(RBEvent::from_bytestream(bytestream, pos).unwrap_or(broken_event.clone()));
1519      }
1520    } else {
1521      // we can skip the next 6 bytes, 
1522      // which just contain 0
1523      *pos += 6;
1524    }
1525    if parse_u16(bytestream, pos) != Self::TAIL {
1526      return Err(SerializationError::TailInvalid {});
1527    }
1528    Ok(rb_cal)
1529  }
1530
1531  fn to_bytestream(&self) -> Vec<u8> {
1532    let mut bs = Vec::<u8>::with_capacity(Self::SIZE);
1533    bs.extend_from_slice(&Self::HEAD.to_le_bytes());
1534    bs.extend_from_slice(&self.rb_id.to_le_bytes());
1535    bs.extend_from_slice(&self.d_v.to_le_bytes());
1536    bs.extend_from_slice(&self.timestamp.to_le_bytes());
1537    let serialize_event_data = self.serialize_event_data as u8;
1538    bs.push(serialize_event_data);
1539    for ch in 0..NCHN {
1540      for k in 0..NWORDS {
1541        bs.extend_from_slice(&self.v_offsets[ch][k].to_le_bytes());
1542        bs.extend_from_slice(&self.v_dips[ch][k]   .to_le_bytes());
1543        bs.extend_from_slice(&self.v_inc[ch][k]    .to_le_bytes());
1544        bs.extend_from_slice(&self.tbin[ch][k]     .to_le_bytes());
1545      }
1546    }
1547    if self.serialize_event_data {
1548      info!("Serializing calibration event data!");
1549      let n_noi  = self.noi_data.len()  as u16;
1550      let n_vcal = self.vcal_data.len() as u16;
1551      let n_tcal = self.tcal_data.len() as u16;
1552      bs.extend_from_slice(&n_noi.to_le_bytes());
1553      for ev in &self.noi_data {
1554        bs.extend_from_slice(&ev.to_bytestream());
1555      }
1556      bs.extend_from_slice(&n_vcal.to_le_bytes());
1557      for ev in &self.vcal_data {
1558        bs.extend_from_slice(&ev.to_bytestream());
1559      }
1560      bs.extend_from_slice(&n_tcal.to_le_bytes());
1561      for ev in &self.tcal_data {
1562        bs.extend_from_slice(&ev.to_bytestream());
1563      }
1564    } else { // if we don't serialize event data, write 0 
1565             // for the empty data
1566      // (3 16bit 0s) for noi, vcal, tcal
1567      for _ in 0..6 {
1568        bs.push(0);
1569      }
1570      //bs.push(0); // noi data
1571      //bs.push(0); // vcal data
1572      //bs.push(0); // tcal data
1573    }
1574    bs.extend_from_slice(&Self::TAIL.to_le_bytes());
1575    bs
1576  }
1577}
1578
1579impl Default for RBCalibrations {
1580  fn default() -> Self {
1581    Self::new(0)
1582  }
1583}
1584
1585impl fmt::Display for RBCalibrations {
1586  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
1587    let mut timestamp_str = String::from("?");
1588    match Utc.timestamp_opt(self.timestamp.into(), 0) {
1589      LocalResult::Single(datetime_utc) => {
1590        timestamp_str = datetime_utc.format("%Y/%m/%d %H:%M:%S").to_string();
1591      },
1592      LocalResult::Ambiguous(_, _) => {
1593        println!("The given timestamp is ambiguous.");
1594      },
1595      LocalResult::None => {
1596        println!("The given timestamp is not valid.");
1597      },
1598    }
1599
1600    //let datetime_utc: DateTime<Utc> = Utc.timestamp_opt(self.timestamp as i64, 0).earliest().unwrap_or(DateTime::<Utc>::from_timestamp_millis(0).unwrap());
1601    if !self.check() {
1602      return write!(f, "<RBCalibrations [{} UTC] for board {}:
1603   -- fails self checks!>", timestamp_str, self.rb_id);
1604    }  
1605    write!(f, 
1606  "<RBCalibrations [{} UTC]:
1607      ** all self checks passed! ** 
1608      RB             : {}
1609      VCalData       : {} (events)
1610      TCalData       : {} (events)
1611      NoInputData    : {} (events)
1612      V Offsets [ch0]: .. {:?} {:?} ..
1613      V Incrmts [ch0]: .. {:?} {:?} ..
1614      V Dips    [ch0]: .. {:?} {:?} ..
1615      T Bins    [ch0]: .. {:?} {:?} ..>",
1616      timestamp_str,
1617      self.rb_id,
1618      self.vcal_data.len(),
1619      self.tcal_data.len(),
1620      self.noi_data.len(),
1621      self.v_offsets[0][98],
1622      self.v_offsets[0][99],
1623      self.v_inc[0][98],
1624      self.v_inc[0][99],
1625      self.v_dips[0][98],
1626      self.v_dips[0][99],
1627      self.tbin[0][98],
1628      self.tbin[0][99])
1629  } 
1630}
1631
1632#[cfg(feature = "pybindings")] 
1633#[pymethods]
1634impl RBCalibrations {
1635  #[getter]
1636  fn rb_id(&self) -> u8 {
1637    self.rb_id
1638  }
1639
1640  #[getter]
1641  fn d_v(&self) -> f32 {
1642    self.d_v
1643  }
1644
1645  #[getter]
1646  fn vcal_data(&self) -> Vec<RBEvent> {
1647    self.vcal_data.clone()
1648  }
1649  
1650  #[getter]
1651  fn tcal_data(&self) -> Vec<RBEvent> {
1652    self.tcal_data.clone()
1653  }
1654  
1655  #[getter]
1656  fn noi_data(&self) -> Vec<RBEvent> {
1657    self.noi_data.clone()
1658  }
1659 
1660  #[getter]
1661  fn v_offsets<'_py>(&self, py: Python<'_py>) -> PyResult<Bound<'_py, PyArray2<f32>>> {  
1662    let mut data = Vec::<Vec<f32>>::with_capacity(9);
1663    for ch in 0..9 {
1664      data.push(self.v_offsets[ch].to_vec());
1665    }
1666    let pyarray = PyArray2::from_vec2(py, &data).unwrap();
1667    Ok(pyarray)
1668  }
1669  
1670  #[getter]
1671  fn v_dips<'_py>(&self, py: Python<'_py>) -> PyResult<Bound<'_py, PyArray2<f32>>> {  
1672    let mut data = Vec::<Vec<f32>>::with_capacity(9);
1673    for ch in 0..9 {
1674      data.push(self.v_dips[ch].to_vec());
1675    }
1676    let pyarray = PyArray2::from_vec2(py, &data).unwrap();
1677    Ok(pyarray)
1678  }
1679  
1680  #[getter]
1681  fn v_inc<'_py>(&self, py: Python<'_py>) -> PyResult<Bound<'_py, PyArray2<f32>>> {  
1682    let mut data = Vec::<Vec<f32>>::with_capacity(9);
1683    for ch in 0..9 {
1684      data.push(self.v_inc[ch].to_vec());
1685    }
1686    let pyarray = PyArray2::from_vec2(py, &data).unwrap();
1687    Ok(pyarray)
1688  }
1689  
1690  #[getter]
1691  fn tbin<'_py>(&self, py: Python<'_py>) -> PyResult<Bound<'_py, PyArray2<f32>>> {  
1692    let mut data = Vec::<Vec<f32>>::with_capacity(9);
1693    for ch in 0..9 {
1694      data.push(self.tbin[ch].to_vec());
1695    }
1696    let pyarray = PyArray2::from_vec2(py, &data).unwrap();
1697    Ok(pyarray)
1698  }
1699  
1700  /// Apply the voltage calibration to a single channel 
1701  /// FIXME - mixing of naming conventions for the channels
1702  ///
1703  /// FIXME - make it return Result<(), CalibrationError>
1704  ///
1705  /// # Arguments
1706  ///
1707  /// * channel   : Channel id 1-9
1708  /// * stop_cell : This channels stop cell 
1709  /// * adc       : Uncalibrated channel data
1710  /// * waveform  : Pre-allocated array to hold 
1711  ///               calibrated waveform data.
1712  #[pyo3(name="voltages")]
1713  pub fn voltages_py<'_py>(&self,
1714                           py        : Python<'_py>,
1715                           channel   : usize,
1716                           stop_cell : usize,
1717                           adc       : Bound<'_py, PyArray1<u16>>)
1718      -> PyResult<Bound<'_py, PyArray1<f32>>>{
1719                  //waveform  : &mut [f32;NWORDS]) {
1720    let mut voltages = vec![0.0f32; 1024];
1721    let adc_data = adc.to_vec().unwrap();
1722    self.voltages(channel, stop_cell, &adc_data, &mut voltages); 
1723    let pyarray = PyArray1::from_vec(py, voltages);
1724    Ok(pyarray)
1725  }
1726  
1727  /// Apply the timing calibration to a single channel 
1728  /// 
1729  /// This will allocate the array for the waveform 
1730  /// time bins (unit is ns)
1731  ///
1732  /// # Arguments
1733  ///
1734  /// * channel   : Channel id 1-9
1735  /// * stop_cell : This channels stop cell 
1736  #[pyo3(name="nanoseconds")]
1737  pub fn nanoseconds_py<'_py>(&self,
1738                              py        : Python<'_py>,                           
1739                              channel   : usize,
1740                              stop_cell : usize) 
1741      -> PyResult<Bound<'_py, PyArray1<f32>>> {
1742    let mut times = vec![0.0f32; 1024];
1743    self.nanoseconds(channel, stop_cell, &mut times); 
1744    let pyarray = PyArray1::from_vec(py, times);
1745    Ok(pyarray)
1746  }
1747
1748  /// Load the calibration from a file with a 
1749  /// TofPacket of type RBCalibration in it
1750  ///
1751  /// # Arguments:
1752  ///
1753  /// * filename     : File with a TofPacket of type RBCalibration in it
1754  /// * discard_data : Throw away event data after loading
1755  #[pyo3(name = "from_file", signature = (filename, discard_data = true))]
1756  #[staticmethod]
1757  fn from_file_py(filename : String, discard_data : bool) -> PyResult<Self> {
1758    let cali = RBCalibrations::from_file(filename, discard_data);
1759    match cali {
1760      Ok(c) => {
1761        return Ok(c);
1762      },
1763      Err(err) => {
1764        return Err(PyValueError::new_err(err.to_string()));
1765      }
1766    }
1767  }
1768
1769  #[pyo3(name="passes_voltage_checks")]
1770  fn passes_voltage_checks_py(&self) -> bool {
1771    self.passes_voltage_checks() 
1772  }
1773
1774  #[pyo3(name="passes_timing_checks")]
1775  fn passes_timing_checks_py(&self) -> bool {
1776    self.passes_timing_checks()
1777  }
1778
1779  #[pyo3(name="check")]
1780  fn check_py(&self) -> bool {
1781    self.check()
1782  }
1783}
1784
1785#[cfg(feature = "pybindings")]
1786pythonize_packable_no_new!(RBCalibrations);
1787
1788#[cfg(feature = "random")]
1789impl FromRandom for RBCalibrations {
1790    
1791  fn from_random() -> Self {
1792    let mut cali   = Self::new(0);
1793    let mut rng    = rand::rng();
1794    cali.rb_id     = rng.random::<u8>();
1795    cali.d_v       = rng.random::<f32>(); 
1796    cali.serialize_event_data = rng.random::<bool>();
1797    for ch in 0..NCHN {
1798      for n in 0..NWORDS { 
1799        cali.v_offsets[ch][n] = rng.random::<f32>();
1800        cali.v_dips   [ch][n] = rng.random::<f32>(); 
1801        cali.v_inc    [ch][n] = rng.random::<f32>(); 
1802        cali.tbin     [ch][n] = rng.random::<f32>();
1803      }
1804    }
1805    if cali.serialize_event_data {
1806      for _ in 0..1000 {
1807        let mut ev = RBEvent::from_random();
1808        cali.vcal_data.push(ev);
1809        ev = RBEvent::from_random();
1810        cali.noi_data.push(ev);
1811        ev = RBEvent::from_random();
1812        cali.tcal_data.push(ev);
1813      }
1814    }
1815    cali
1816  }
1817}
1818
1819//-----------------------------------------------
1820
1821#[cfg(feature = "random")]
1822#[test]
1823fn serialization_rbcalibration_noeventpayload() {
1824  let mut calis = Vec::<RBCalibrations>::new();
1825  for _ in 0..100 {
1826    let cali = RBCalibrations::from_random();
1827    if cali.serialize_event_data {
1828      continue;
1829    }
1830    calis.push(cali);
1831    break;
1832  }
1833  let test = RBCalibrations::from_bytestream(&calis[0].to_bytestream(), &mut 0).unwrap();
1834  assert_eq!(calis[0], test);
1835}
1836
1837#[cfg(feature = "random")]
1838#[test]
1839fn serialization_rbcalibration_witheventpayload() {
1840  let mut calis = Vec::<RBCalibrations>::new();
1841  for _ in 0..100 {
1842    let cali = RBCalibrations::from_random();
1843    if !cali.serialize_event_data {
1844      continue;
1845    }
1846    calis.push(cali);
1847    break;
1848  }
1849  let test = RBCalibrations::from_bytestream(&calis[0].to_bytestream(), &mut 0).unwrap();
1850  assert_eq!(calis[0], test);
1851}
1852
1853#[cfg(feature = "random")]
1854#[test]
1855fn pack_rbcalibrationfilghtt() {
1856  for _ in 0..100 {
1857    let cfg  = RBCalibrationFlightT::from_random();
1858    let test : RBCalibrationFlightT = cfg.pack().unpack().unwrap();
1859    assert_eq!(cfg, test);
1860  }
1861}
1862
1863#[cfg(feature = "random")]
1864#[test]
1865fn pack_rbcalibfilghtv() {
1866  for _ in 0..100 {
1867    let cfg  = RBCalibrationFlightV::from_random();
1868    let test : RBCalibrationFlightV = cfg.pack().unpack().unwrap();
1869    assert_eq!(cfg, test);
1870  }
1871}
1872
1873#[cfg(feature = "random")]
1874#[test]
1875fn assemble_flightcal() {
1876  for _ in 0..10 {
1877    let cal  = RBCalibrations::from_random();
1878    let fct  = cal.emit_flighttcal();
1879    let fcv  = cal.emit_flightvcal();
1880    let test = RBCalibrations::assemble_from_flightcal(fct, fcv).unwrap();
1881    assert_eq!(cal.rb_id, test.rb_id);
1882    assert_eq!(cal.d_v  , test.d_v);
1883    assert_eq!(cal.timestamp, test.timestamp);
1884    assert_eq!(test.serialize_event_data, false);
1885    for ch in 0..NCHN {
1886      for k in 0..NWORDS {
1887        assert_eq!(f16::from_f32(cal.tbin[ch][k]),     f16::from_f32(test.tbin[ch][k])); 
1888        assert_eq!(f16::from_f32(cal.v_offsets[ch][k]),f16::from_f32(test.v_offsets[ch][k])); 
1889        assert_eq!(f16::from_f32(cal.v_dips[ch][k]),   f16::from_f32(test.v_dips[ch][k])); 
1890        assert_eq!(f16::from_f32(cal.v_inc[ch][k]),    f16::from_f32(test.v_inc[ch][k]));
1891      }
1892    }
1893  }
1894}
1895
1896
gondola_core/calibration/tof.rs

gondola_core/calibration/
tof.rs
