gondola_core/database/
tracker_transfer_fn.rs

1//! Database access & entities for gaps-online-software
2//!
3//! A local .sqlite database is shipped with gaps-online-software,
4//! pre-populated with relevant meta data for the GAPS experiment.
5// This file is part of gaps-online-software and published 
6// under the GPLv3 license
7
8use crate::prelude::*;
9
10use crate::database::schema;
11use diesel::prelude::*;
12
13use std::io::{
14  self,
15  BufRead,
16  BufReader
17};
18
19/// Tracker transfer functions connect the tracker adc to a measurement of energy
20#[derive(Debug,PartialEq, Clone,Queryable, Selectable, serde::Serialize, serde::Deserialize)]
21#[diesel(table_name = schema::tof_db_trackerstriptransferfunction)]
22#[diesel(primary_key(data_id))]
23#[allow(non_snake_case)]
24#[cfg_attr(feature="pybindings", pyclass)]
25pub struct TrackerStripTransferFunction {  
26    pub data_id            : i32,
27    pub strip_id           : i32,    
28    pub volume_id          : i64,    
29    pub utc_timestamp_start: i64,    
30    pub utc_timestamp_stop : i64,    
31    pub name               : Option<String>, 
32    pub pol_a2_0           : f32, 
33    pub pol_a2_1           : f32,    
34    pub pol_a2_2           : f32, 
35    pub pol_b3_0           : f32, 
36    pub pol_b3_1           : f32, 
37    pub pol_b3_2           : f32, 
38    pub pol_b3_3           : f32, 
39    pub pol_c3_0           : f32, 
40    pub pol_c3_1           : f32, 
41    pub pol_c3_2           : f32, 
42    pub pol_c3_3           : f32, 
43    pub pol_d3_0           : f32,     
44    pub pol_d3_1           : f32, 
45    pub pol_d3_2           : f32, 
46    pub pol_d3_3           : f32, 
47} 
48
49impl TrackerStripTransferFunction {
50
51  pub fn new() -> Self {
52    Self {
53      data_id             : 0,
54      strip_id            : 0,    
55      volume_id           : 0,    
56      utc_timestamp_start : 0,    
57      utc_timestamp_stop  : 0,
58      name                : None, 
59      pol_a2_0            : 0.0, 
60      pol_a2_1            : 0.0,    
61      pol_a2_2            : 0.0, 
62      pol_b3_0            : 0.0, 
63      pol_b3_1            : 0.0, 
64      pol_b3_2            : 0.0, 
65      pol_b3_3            : 0.0, 
66      pol_c3_0            : 0.0, 
67      pol_c3_1            : 0.0, 
68      pol_c3_2            : 0.0, 
69      pol_c3_3            : 0.0, 
70      pol_d3_0            : 0.0,     
71      pol_d3_1            : 0.0, 
72      pol_d3_2            : 0.0, 
73      pol_d3_3            : 0.0, 
74    }
75  }
76  
77  pub fn parse_from_file<P: AsRef<Path>>(path: P) -> io::Result<Vec<Self>> {
78    let file = File::open(path)?;
79    let reader = BufReader::new(file);
80    let mut transfer_fns = Vec::<Self>::new();
81    let hid_vid_map = get_hid_vid_maps().unwrap().1;
82    for line in reader.lines() {
83      let line = line?;
84      if line.starts_with("#") || line.starts_with("Layer") || line.starts_with("layer") {
85        continue;
86      }
87      let mut trfn = Self::new();
88      let parts: Vec<&str> = line.split(",").collect();
89      if parts.len() == 19 {
90        // Parse the first number as a standard decimal
91        let layer   = parts[0].parse::<u8>()
92          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
93        let row     = parts[1].parse::<u8>()
94          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
95        let module   = parts[2].parse::<u8>()
96          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
97        let channel  = parts[3].parse::<u8>()
98          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
99        trfn.strip_id  = TrackerStrip::create_stripid(layer, row, module, channel) as i32; 
100        trfn.volume_id = *hid_vid_map.get(&(trfn.strip_id as u32)).unwrap() as i64; // critical error is good here,
101        trfn.pol_a2_0  = parts[4].parse::<f32>()
102          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
103        trfn.pol_a2_1  = parts[5].parse::<f32>()
104          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
105        trfn.pol_a2_2  = parts[6].parse::<f32>()
106          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
107        trfn.pol_b3_0  = parts[7].parse::<f32>()
108          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
109        trfn.pol_b3_1  = parts[8].parse::<f32>()
110          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
111        trfn.pol_b3_2  = parts[9].parse::<f32>()
112          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
113        trfn.pol_b3_3  = parts[10].parse::<f32>()
114          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
115        trfn.pol_c3_0  = parts[11].parse::<f32>()
116          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
117        trfn.pol_c3_1  = parts[12].parse::<f32>()
118          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
119        trfn.pol_c3_2  = parts[13].parse::<f32>()
120          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
121        trfn.pol_c3_3  = parts[14].parse::<f32>()
122          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
123        trfn.pol_d3_0  = parts[15].parse::<f32>()
124          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
125        trfn.pol_d3_1  = parts[16].parse::<f32>()
126          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
127        trfn.pol_d3_2  = parts[17].parse::<f32>()
128          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
129        trfn.pol_d3_3  = parts[18].parse::<f32>()
130          .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))?;
131        transfer_fns.push(trfn); 
132      }
133    }
134    Ok(transfer_fns)
135  }
136  
137  /// Get Tracker strip transfer fns for a certain dataset 
138  ///
139  /// # Returns:
140  ///   * HashMap<u32 [strip id], TrackerStripTransferFn> 
141  pub fn as_dict_by_name(fname : &str) -> Result<HashMap<u32,Self>, ConnectionError> {
142    use schema::tof_db_trackerstriptransferfunction::dsl::*;
143    let mut strips = HashMap::<u32, Self>::new();
144    if fname == "" {
145      match Self::all() {
146        None => {
147          error!("Unable to retrive ANY TrackerStripTransferFunction");
148          return Ok(strips);
149        }
150        Some(_strips) => {
151          for k in _strips {
152            strips.insert(k.strip_id as u32, k);
153          }
154          return Ok(strips);
155        }
156      }
157    }
158    let mut conn = connect_to_db()?;
159    match tof_db_trackerstriptransferfunction.filter(
160      schema::tof_db_trackerstriptransferfunction::name.eq(fname)).load::<Self>(&mut conn) {
161      Err(err) => {
162        error!("We can't find any tracker strip transferfunction in the database! {err}");
163        return Ok(strips);
164      }
165      Ok(peds_) => {
166        for s in peds_ {
167          strips.insert(s.strip_id as u32, s );
168        }
169      }
170    }
171    return Ok(strips);
172  }
173  
174  pub fn all_names() -> Result<Vec<String>, ConnectionError> {
175    let mut conn = connect_to_db()?;
176    let mut names = Vec::<String>::new();
177    let unique_names =
178      schema::tof_db_trackerstriptransferfunction::table.select(
179      schema::tof_db_trackerstriptransferfunction::name)
180      .distinct()
181      .load::<Option<String>>(&mut conn).expect("Error getting names from db!");
182    for k in unique_names {
183      if let Some(n) = k {
184        names.push(n);
185      }
186    }
187    Ok(names)
188  }
189
190  /// Get all tracker strip transfer functions from the database
191  ///
192  /// # Returns:
193  ///   * HashMap<u32 [strip id], TrackeStripTransferFunction> 
194  pub fn all() -> Option<Vec<Self>> {
195    use schema::tof_db_trackerstriptransferfunction::dsl::*;
196    let mut conn = connect_to_db().ok()?;
197    match tof_db_trackerstriptransferfunction.load::<Self>(&mut conn) {
198      Err(err) => {
199        error!("Unable to load tracker transfer functions from db! {err}");
200        return None;
201      }
202      Ok(strips) => {
203        return Some(strips);
204      }
205    }
206  }
207
208  /// The actual transfer function for this 
209  /// strip. Calculate energy from adc values
210  pub fn evaluate(&self, adc : f32) -> f32 {
211    if adc < 0.0 {
212      return 0.0;
213    }
214    if adc <= 190.0 {
215      return self.pol_a2_0 + self.pol_a2_1*adc + self.pol_a2_2*(adc.powi(2));
216    }
217    if 190.0 < adc && adc <= 500.0 {
218      return self.pol_b3_0 + self.pol_b3_1*adc + self.pol_b3_2*(adc.powi(2)) + self.pol_b3_3*(adc.powi(3));
219    }
220    if 500.0 < adc && adc <= 900.0 {
221      return self.pol_c3_0 + self.pol_c3_1*adc + self.pol_c3_2*(adc.powi(2)) + self.pol_c3_3*(adc.powi(3));
222    }
223    //if 900.0 < adc && adc <= 2047.0 {
224    if 900.0 < adc && adc <= 1600.0 {
225      return self.pol_d3_0 + self.pol_d3_1*adc + self.pol_d3_2*(adc.powi(2)) + self.pol_d3_3*(adc.powi(3));
226    }
227    0.0
228  }
229 
230  /// First derivative of the transfer function
231  pub fn derivative(&self, adc : f32) -> f32 {
232    if adc < 0.0 {
233      return 0.0;
234    }
235    if adc <= 190.0 {
236      return self.pol_a2_1 + 2.0*self.pol_a2_2*adc;
237    }
238    if 190.0 < adc && adc <= 500.0 {
239      return self.pol_b3_1 + 2.0*self.pol_b3_2*adc + 3.0*self.pol_b3_3*(adc.powi(2));
240    }
241    if 500.0 < adc && adc <= 900.0 {
242      return self.pol_c3_1 + 2.0*self.pol_c3_2 + 3.0*self.pol_c3_3*(adc.powi(2));
243    }
244    //if 900.0 < adc && adc <= 2047.0 {
245    if 900.0 < adc && adc <= 1600.0 {
246      return self.pol_d3_1 + 2.0*self.pol_d3_2 + 3.0*self.pol_d3_3*(adc.powi(2));
247    }
248    0.0
249  }
250
251  /// Create the inverse of the transfer function, which then 
252  /// can be used in the simulation
253  /// FIXME - switch to f64, however, only store f32 in the db
254  pub fn invert(&self, y : f64, epsilon_opt : Option<f64>) -> Result<f64, AnalysisError> {
255    let x_min = 0.0f32;
256    let x_max = 1600.0f32;
257    let y_min = self.evaluate(x_min);
258    let y_max = self.evaluate(x_max);
259    if y < y_min as f64 || y > y_max as f64 {
260      error!("y value of {} is out of bounds for 0 <= x <= 1600", y);
261      return Err(AnalysisError::OutOfBounds);
262    }
263    let max_iterations = 50;
264    let epsilon : f64;
265    if epsilon_opt.is_some() {
266      epsilon = epsilon_opt.unwrap();
267    } else {
268      epsilon = 1e-12; // random default
269    }
270
271    let mut low  = x_min;
272    let mut high = x_max;
273    let mut x = 0.5 * (low + high);
274
275    for _ in 0..max_iterations {
276      let f_x = self.evaluate(x) as f64 - y;
277      let df_x = self.derivative(x) as f64;
278      // If we are close enough, or if the derivative is effectively zero
279      if f_x.abs() < epsilon {
280        return Ok(x as f64);
281      }
282      // Update bisection bounds
283      if f_x > 0.0 {
284        high = x;
285      } else {
286        low = x;
287      }
288      // Attempt a Newton-Raphson step
289      if df_x.abs() > 1e-9 {
290        let next_x = x as f64 - f_x / df_x;
291        // Safety check: if Newton step jumps outside our bounded bracket,
292        // fallback to a safe Bisection step.
293        if next_x > low as f64 && next_x < high as f64 {
294          x = next_x as f32;
295          continue;
296        }
297      }
298      // Fallback to Bisection if Newton fails or is too slow
299      x = 0.5 * (low + high);
300    }
301    error!("We couldn't find a solution which is epsilon {} away! Try to reduce epsilon!", epsilon);
302    Err(AnalysisError::DidNotConverge)
303  }
304}
305
306impl fmt::Display for TrackerStripTransferFunction {
307  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
308    let mut repr = format!("<TrackerStripTransferFunction [{}]:", self.strip_id);
309    repr += &(format!("\n   vid           : {}", self.volume_id));
310    repr += "\n   UTC Timestamps (Begin/End):";
311    repr += &(format!("\n   {}/{}", self.utc_timestamp_start, self.utc_timestamp_stop));    
312    if self.name.is_some() {
313      repr += &(format!("\n   name     : {}", self.name.clone().unwrap())); 
314    }
315    repr += &(format!("\n  Poly A {}*adc + {}*adc + {}*(adc**2) for adc < 190", self.pol_a2_0, self.pol_a2_1, self.pol_a2_2));
316    repr += &(format!("\n  Poly B    :{}*adc + {}*adc + {}*(adc**2) + {}*(adc**3) for 190 < adc <= 500", self.pol_b3_0, self.pol_b3_1, self.pol_b3_2, self.pol_b3_3));
317    repr += &(format!("\n  Poly C    :{}*adc + {}*adc + {}*(adc**2) + {}*(adc**3) for 500 < adc <= 900", self.pol_c3_0, self.pol_c3_1, self.pol_c3_2, self.pol_c3_3));
318    repr += &(format!("\n  Poly D    :{}*adc + {}*adc + {}*(adc**2) + {}*(adc**3) for 900 < adc <= 1600>", self.pol_d3_0, self.pol_d3_1, self.pol_d3_2, self.pol_d3_3));
319    write!(f, "{}", repr)
320  }
321}
322
323#[cfg(feature="pybindings")]
324#[pymethods]
325impl TrackerStripTransferFunction {
326  
327  #[staticmethod]
328  #[pyo3(name="all")]
329  pub fn all_py() -> Option<Vec<Self>> {
330    Self::all()
331  } 
332  
333  #[staticmethod]
334  #[pyo3(name="all_names")]
335  /// Get all names for registered datasets. These
336  /// can be used in .as_dict_by_name() to query 
337  /// the db for a set of values
338  pub fn all_names_py() -> Option<Vec<String>> {
339    match Self::all_names() {
340      Err(_) => {
341        return None;
342      }
343      Ok(names) => {
344        return Some(names);
345      }
346    }
347  } 
348  
349  #[staticmethod]
350  #[pyo3(name="as_dict_by_name")]
351  pub fn all_as_dict_py(name : &str) -> Option<HashMap<u32,Self>> {
352    match Self::as_dict_by_name(name) {
353      Err(err) => {
354        error!("Unable to retrieve tracker strip transfer fn dictionary. {err}. Did you laod the setup-env.sh shell?");
355        return None;
356      }
357      Ok(_data) => {
358        return Some(_data);
359      }
360    }
361  } 
362  
363  #[getter]
364  fn get_strip_id     (&self) -> i32 {    
365    self.strip_id
366  }
367  
368  #[getter]
369  fn get_volume_id    (&self) -> i64 {    
370    self.volume_id
371  }
372  
373  #[getter]
374  fn get_utc_timestamp_start(&self) -> i64 {
375    self.utc_timestamp_start
376  }
377  
378  #[getter]
379  fn get_utc_timestamp_stop(&self) -> i64 {
380    self.utc_timestamp_stop
381  }
382  
383  #[setter]
384  fn set_utc_timestamp_start(&mut self, value : i64) {
385    self.utc_timestamp_start = value;
386  }
387  
388  #[setter]
389  fn set_utc_timestamp_stop(&mut self, value : i64) {
390    self.utc_timestamp_stop = value;
391  }
392
393  #[getter]
394  fn get_name(&self) -> Option<String> {
395    self.name.clone()
396  }
397  
398  #[setter]
399  fn set_name(&mut self, value : String) {
400    self.name = Some(value);
401  }
402
403  #[pyo3(name="transfer_fn")]
404  fn transfer_fn_py(&self, adc : f32) -> f32 {
405    if adc > 1600.0 {
406      warn!("ADC value larger than 1600! {}. Transfer fn not defined beyond 1600.", adc);
407    }
408    return self.evaluate(adc);
409  }
410 
411  #[pyo3(name="derivative")]
412  fn derivative_py(&self, adc : f32) -> f32 {
413    if adc > 1600.0 {
414      warn!("ADC value larger than 1600! {}. Transfer fn not defined beyond 1600.", adc);
415    }
416    return self.derivative(adc);
417  }
418
419  /// The inverse transfer function as to be used 
420  /// in the simulation. 
421  /// If this returns errors for values within 
422  /// range, try to reduce the convergence 
423  /// crtierion epsilon 
424  #[pyo3(signature = (mv , epsilon = None))] 
425  #[pyo3(name="invert")]
426  fn invert_py(&self, mv : f64, epsilon : Option<f64>) -> PyResult<f64> {
427    Ok(self.invert(mv, epsilon)?)
428  }
429
430  #[staticmethod]
431  #[pyo3(name="parse_from_file")]
432  fn parse_from_file_py(fname : &str) -> Option<Vec<Self>> {
433    let result = Self::parse_from_file(fname);
434    if result.is_ok() {
435      return Some(result.unwrap());
436    } else {
437      error!("An error occured when parsing {} : '{}'", fname, result.unwrap_err());
438      return None;
439    }
440  }
441}
442
443#[cfg(feature="pybindings")]
444pythonize!(TrackerStripTransferFunction);
445
446//-------------------------------------------------
447
448#[derive(Debug,PartialEq, Clone, Insertable)]
449#[diesel(table_name = schema::tof_db_trackerstriptransferfunction)]
450#[allow(non_snake_case)]
451#[cfg_attr(feature="pybindings", pyclass)]
452struct NewTrackerStripTransferFunction {  
453  pub strip_id           : i32,    
454  pub volume_id          : i64,    
455  pub utc_timestamp_start: i64,    
456  pub utc_timestamp_stop : i64,    
457  pub name               : Option<String>, 
458  pub pol_a2_0           : f32, 
459  pub pol_a2_1           : f32,    
460  pub pol_a2_2           : f32, 
461  pub pol_b3_0           : f32, 
462  pub pol_b3_1           : f32, 
463  pub pol_b3_2           : f32, 
464  pub pol_b3_3           : f32, 
465  pub pol_c3_0           : f32, 
466  pub pol_c3_1           : f32, 
467  pub pol_c3_2           : f32, 
468  pub pol_c3_3           : f32, 
469  pub pol_d3_0           : f32,     
470  pub pol_d3_1           : f32, 
471  pub pol_d3_2           : f32, 
472  pub pol_d3_3           : f32, 
473} 
474
475impl NewTrackerStripTransferFunction {
476  pub fn from(tf : &TrackerStripTransferFunction) -> Self {
477    Self {
478      strip_id            : tf.strip_id            ,    
479      volume_id           : tf.volume_id           ,    
480      utc_timestamp_start : tf.utc_timestamp_start ,    
481      utc_timestamp_stop  : tf.utc_timestamp_stop  ,    
482      name                : tf.name.clone()        , 
483      pol_a2_0            : tf.pol_a2_0            , 
484      pol_a2_1            : tf.pol_a2_1            ,    
485      pol_a2_2            : tf.pol_a2_2            , 
486      pol_b3_0            : tf.pol_b3_0            , 
487      pol_b3_1            : tf.pol_b3_1            , 
488      pol_b3_2            : tf.pol_b3_2            , 
489      pol_b3_3            : tf.pol_b3_3            , 
490      pol_c3_0            : tf.pol_c3_0            , 
491      pol_c3_1            : tf.pol_c3_1            , 
492      pol_c3_2            : tf.pol_c3_2            , 
493      pol_c3_3            : tf.pol_c3_3            , 
494      pol_d3_0            : tf.pol_d3_0            ,     
495      pol_d3_1            : tf.pol_d3_1            , 
496      pol_d3_2            : tf.pol_d3_2            , 
497      pol_d3_3            : tf.pol_d3_3            , 
498    }
499  }
500}
501
502#[cfg_attr(feature="pybindings", pyfunction)]
503pub fn create_trk_transfer_fn_table( db_path: &str, transfer_fns: Vec<TrackerStripTransferFunction>) { 
504  use schema::tof_db_trackerstriptransferfunction::dsl::*;
505  let mut conn = SqliteConnection::establish(db_path).ok().unwrap(); 
506  let mut _query_result = diesel::sql_query("
507      CREATE TABLE IF NOT EXISTS tof_db_trackerstriptransferfunction (
508          data_id INTEGER PRIMARY KEY AUTOINCREMENT,
509          strip_id INTEGER NOT NULL,
510          volume_id BIGINT NOT NULL,
511          utc_timestamp_start BIGINT NOT NULL,
512          utc_timestamp_stop BIGINT NOT NULL,
513          name TEXT,
514          pol_a2_0 FLOAT NOT NULL, 
515          pol_a2_1 FLOAT NOT NULL,    
516          pol_a2_2 FLOAT NOT NULL, 
517          pol_b3_0 FLOAT NOT NULL, 
518          pol_b3_1 FLOAT NOT NULL, 
519          pol_b3_2 FLOAT NOT NULL, 
520          pol_b3_3 FLOAT NOT NULL, 
521          pol_c3_0 FLOAT NOT NULL, 
522          pol_c3_1 FLOAT NOT NULL, 
523          pol_c3_2 FLOAT NOT NULL, 
524          pol_c3_3 FLOAT NOT NULL, 
525          pol_d3_0 FLOAT NOT NULL,     
526          pol_d3_1 FLOAT NOT NULL, 
527          pol_d3_2 FLOAT NOT NULL, 
528          pol_d3_3 FLOAT NOT NULL 
529      )
530  ").execute(&mut conn);
531  let mut new_data = Vec::<NewTrackerStripTransferFunction>::new();
532  for p in transfer_fns {
533    let np = NewTrackerStripTransferFunction::from(&p);
534    new_data.push(np);
535  }
536  println!("Will insert {} tfns!", new_data.len());
537  _query_result = diesel::insert_into(tof_db_trackerstriptransferfunction)
538    .values(&new_data)
539    .execute(&mut conn);
540  println!("Inserted {}", _query_result.unwrap());
541}
542
gondola_core/database/tracker_transfer_fn.rs

gondola_core/database/
tracker_transfer_fn.rs
