gondola_core/calibration/
tracker.rs

1//! The following file is part of gaps-online-software and published 
2//! under the GPLv3 license
3//!
4//! Calibration routines for the GAPS tracker system
5
6use crate::prelude::*;
7
8#[cfg_attr(feature="pybindings", pyclass)] 
9pub struct TrackerOfflineCalibration {
10  pub mask_map      : HashMap<u32,TrackerStripMask>, 
11  pub tf_map        : HashMap<u32,TrackerStripTransferFunction>, 
12  pub ped_map       : HashMap<u32,TrackerStripPedestal>,
13  pub gain_map      : HashMap<u32,TrackerStripGain>,
14  pub pulse_map     : HashMap<u32,TrackerStripPulse>,
15  pub adc_sig_cut   : HashMap<u32,f32>,
16  pub remove_cmn    : bool,
17  pub ped_sigma_cut : f32,
18}
19
20impl TrackerOfflineCalibration {
21
22  pub fn new() -> Self {
23    Self {
24      mask_map      : HashMap::<u32, TrackerStripMask>::new(),
25      tf_map        : HashMap::<u32, TrackerStripTransferFunction>::new(), 
26      ped_map       : HashMap::<u32, TrackerStripPedestal>::new(),
27      adc_sig_cut   : HashMap::<u32,f32>::new(),
28      gain_map      : HashMap::<u32,TrackerStripGain>::new(),
29      pulse_map     : HashMap::<u32,TrackerStripPulse>::new(),
30      remove_cmn    : false,
31      ped_sigma_cut : 0.0,
32    }
33  }
34
35  /// Calculate the common noise level, according to the research of the GAPS tracker group
36  ///
37  /// General scheme is to find the pulsed channel on the strip, get the pulse value, 
38  /// subtract the pedestal and subtract it from the hit adc (taking gain into account)
39  pub fn get_common_noise(&self, hit : &TrackerHit, event_hits : &Vec<TrackerHit>) -> (f32, bool) {
40    let mut cmn_level     = 0.0f32;
41    // per default gains are 1 if they can not be looked up
42    let strip_gain     : f32;
43    let mut pulse_gain = 1.0f32; 
44    let mut pulse_avg  = 0.0f32;
45    let mut pulse_chn  = -1i32;
46    let mut pulse_adc  = 0.0f32;
47    let is_pulser      : bool;
48    //let mut pulse_is_mean = false;
49    if let Some(cmn_strip) = self.gain_map.get(&hit.get_stripid()) {
50      strip_gain = cmn_strip.gain;
51    } else {
52      debug!("There is no entry for the strip gain for {}, default is 1.0", &hit.get_stripid());
53      strip_gain = 1.0;
54    }
55    if let Some(pulse_chn_) = self.pulse_map.get(&hit.get_stripid()) { 
56      pulse_chn  = pulse_chn_.pulse_chn;  
57      pulse_avg  = pulse_chn_.pulse_avg;
58      //pulse_is_mean = pulse_chn_.pulse_is_mean;
59    }
60    is_pulser = pulse_chn == hit.channel as i32;
61    if pulse_chn < 0 || pulse_avg > 400.0 {
62      return (0.0, is_pulser) 
63    }
64    let pulse_id = TrackerStrip::create_stripid(hit.layer, hit.row, hit.module, pulse_chn as u8);
65    // now we need to find the hit in this event which is caused by the pulser 
66    let mut p_hit  = TrackerHit::new();
67    p_hit.layer    = hit.layer;
68    p_hit.row      = hit.row;
69    p_hit.module   = hit.module;
70    p_hit.channel  = pulse_chn as u8;
71    let mut found  = false;
72    for h in event_hits {
73      if h.get_stripid() == p_hit.get_stripid() {
74        pulse_adc = h.adc as f32;
75        found = true;
76        break; // FIXME - there might be multiple hits
77      }
78    }
79    if let Some(pulse_gain_) = self.gain_map.get(&pulse_id) {
80      pulse_gain = pulse_gain_.gain; 
81    } 
82    if found {
83      let adc_pulc_diff = pulse_adc - pulse_avg;
84      if adc_pulc_diff < 250.0 && pulse_avg > 0.0 {
85        cmn_level = adc_pulc_diff / pulse_gain;   
86      }
87    } 
88    debug!("Calculated common level of {} gain {}", cmn_level, strip_gain);
89    if cmn_level < 0.0 {
90      //error!("CMN is negative! {}", cmn_level);
91      //cmn_level = 0.0;
92    }
93    return (cmn_level * strip_gain, is_pulser);  
94  }
95
96  pub fn mask_hits(&self, event_hits : &mut Vec<TrackerHit>) {
97    let mut active_strips = Vec::<u32>::new();
98    for h in event_hits.iter() {
99      if let Some(mask) = self.mask_map.get(&h.get_stripid()) {
100        if mask.active {
101          active_strips.push(h.get_stripid());
102        }
103      } else {
104        warn!("No mask information for strip {}", h.get_stripid());
105      }
106    }
107    event_hits.retain(|x| active_strips.contains(&x.get_stripid()));
108  }
109
110  pub fn calibrate_event(&self, event : &mut TelemetryEvent, remove_pulsed : bool) -> Result<(),CalibrationError> {
111    self.calibrate(&mut event.tracker_hits, remove_pulsed)?;
112    Ok(())
113  }
114
115  /// The actual energy calibration of the individual tracker strips 
116  ///
117  /// This is a multistep process which includes 
118  /// 1) Pedestal subtraction
119  /// 2) Common noise reduction (if so desired) 
120  /// 3) ADC -> Energy (Transfer functions) 
121  /// 3) Deal with strips which had the pulser active 
122  ///
123  /// While this might not be the best way to do things, this 
124  /// function matches the implementation in SimpleDet very 
125  /// closely and produces the same results. 
126  ///
127  /// See also the compontents in the `database` related part 
128  /// of this code.
129  ///
130  /// # Arguments:
131  ///   * remove_pulsed: If true, remove pulser hits with adc > 400
132  pub fn calibrate(&self, event_hits : &mut Vec<TrackerHit>, remove_pulsed : bool) -> Result<(),CalibrationError> {
133    let mut calibrated_hits = Vec::<TrackerHit>::with_capacity(event_hits.len());
134    let mut c_hit : TrackerHit; //= TrackerHit::new();
135    let mv_2_kev = 0.841f32;// mV to keV
136    for hit in event_hits.iter() {
137      let mut hit_ped = 0.0f32;
138      let mut energy  = hit.adc as f32;
139      let hit_tf  : Option<&TrackerStripTransferFunction>;
140      if let Some(ped) = self.ped_map.get(&hit.get_stripid()) {
141        hit_ped = ped.pedestal_mean;
142        energy -= hit_ped;
143        //if energy < self.ped_sigma_cut * ped.pedestal_sigma {
144        //  // this basically means that the energy is within the given 
145        //  // range of the pedestal. In this case, we set it basiccally to 0 
146        //  //
147        //  // Note - we don't raise the error, because this simply means this 
148        //  // is a noisy strip!
149        //  //return Ok(());i
150        //  continue;
151        //}
152      } else {
153        error!("No entry for {} in pedestal map", hit);
154        //continue;
155        //return Err(CalibrationError::NoPedestalAvailable); 
156      }
157      let mut hit_is_pulser = false;
158      if self.remove_cmn {
159        let cmn = self.get_common_noise(hit, event_hits); 
160        energy -= cmn.0;
161        hit_is_pulser = cmn.1;
162      }
163      // apply transfer functions 
164      hit_tf = self.tf_map.get(&hit.get_stripid());
165      if let Some(tf) = hit_tf {
166        //println!("energy before : {}", energy);
167        energy = tf.evaluate(energy);
168        //println!("energy after : {}", energy);
169      } else {
170        error!("Trying to calibrate {}, but we don't have a transfer function for that!", hit.get_stripid());
171        energy = 0.0;
172        //return Err(CalibrationError::NoTransferFnAvailable);
173      }
174      // now we need to check if this is a pulsed channel 
175      if hit_is_pulser {
176        if let Some(pulse) = self.pulse_map.get(&hit.get_stripid()) {
177          let mut p_avg = pulse.pulse_avg; 
178          if p_avg > 0.0 {
179            if p_avg > 400.0 && remove_pulsed {
180              continue;
181            }
182            p_avg -= hit_ped;
183            if let Some(tf) = hit_tf {
184              p_avg = tf.evaluate(p_avg); 
185              energy -= p_avg;
186            }
187          }
188        } else {
189          warn!("Hit is from the pulser, but we don't have any information about that strip in the cmn map!"); 
190        }
191      }
192      c_hit = hit.clone();
193      energy       *= mv_2_kev/1000.0; 
194      c_hit.energy = energy;
195      //println!("c_hit {}", c_hit);
196      calibrated_hits.push(c_hit);
197      //if energy < -0.13 {
198      //  error!("Small energy!");
199      //  error!("{}", c_hit);
200      //}
201    }
202    event_hits.clear(); 
203    *event_hits = calibrated_hits;
204    Ok(())
205  }
206}
207
208impl fmt::Display for TrackerOfflineCalibration {
209  fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
210    let mut repr = String::from("<TrackerOfflineCalibration:");
211    repr += &(format!("\n  N masks  : {}", self.mask_map.len()));
212    repr += &(format!("\n  N pedest : {}", self.ped_map.len()));
213    repr += &(format!("\n  N tf_fns : {}", self.tf_map.len()));
214    repr += &(format!("\n  N gains  : {}", self.gain_map.len()));
215    repr += &(format!("\n  N pulses : {}", self.pulse_map.len()));
216    repr += ">";
217    write!(f, "{}", repr)
218  }
219}
220
221#[cfg(feature="pybindings")] 
222#[pymethods]
223impl TrackerOfflineCalibration {
224 
225  #[pyo3(name="calibrate_hits")]
226  fn calibrate_hits_py(&self, mut hits : Vec<TrackerHit>, remove_pulsed : bool) -> Vec<TrackerHit> {
227    let _ = self.calibrate(&mut hits, remove_pulsed);
228    hits 
229  }
230   
231  // all the getters/setters are here because I did not get
232  // #[cfg_attr(feature="pybindings", pyo3(set,get)] 
233  // to work. If we find out how to use that, get rid 
234  // of all the getters/setters here
235
236  #[getter]
237  #[pyo3(name="remove_cmn")] 
238  fn get_remove_cmn_py(&self) -> bool {
239    self.remove_cmn
240  }
241  #[getter]
242  #[pyo3(name="ped_sigma_cut")] 
243  fn get_ped_sigma_cut_py(&self) -> f32 {
244    self.ped_sigma_cut
245  }
246 
247  #[setter]
248  #[pyo3(name="remove_cmn")] 
249  fn set_remove_cmn_py(&mut self, value : bool) {
250    self.remove_cmn = value;
251  }
252
253  #[setter]
254  #[pyo3(name="ped_sigma_cut")] 
255  fn set_ped_sigma_cut_py(&mut self, value : f32) {
256    self.ped_sigma_cut = value;
257  }
258
259  #[getter]
260  #[pyo3(name="mask_map")]
261  fn get_mask_map_py(&self) -> HashMap<u32,TrackerStripMask> {
262    self.mask_map.clone()
263  }
264  
265  #[setter]
266  #[pyo3(name="mask_map")]
267  fn set_mask_map_py(&mut self, value : HashMap<u32,TrackerStripMask>) {
268    self.mask_map = value;
269  }
270  
271  #[getter]
272  #[pyo3(name="tf_map")]
273  fn get_tf_map_py(&self) -> HashMap<u32,TrackerStripTransferFunction> {
274    self.tf_map.clone()
275  }
276  #[setter]
277  #[pyo3(name="tf_map")]
278  fn set_tf_map_py(&mut self, value : HashMap<u32, TrackerStripTransferFunction>) {
279    self.tf_map = value;
280  }
281
282  #[getter]
283  #[pyo3(name="ped_map")]
284  fn get_ped_map_py(&self) -> HashMap<u32,TrackerStripPedestal> {
285    self.ped_map.clone()
286  }
287  #[setter] 
288  #[pyo3(name="ped_map")]
289  fn set_ped_map_py(&mut self, value : HashMap<u32,TrackerStripPedestal>) {
290    self.ped_map = value;
291  }
292
293  #[getter]
294  #[pyo3(name="pulse_map")]
295  fn get_pulse_map_py(&self) -> HashMap<u32, TrackerStripPulse> {
296    self.pulse_map.clone()
297  }
298  
299  #[setter]
300  #[pyo3(name="pulse_map")]
301  fn set_pulse_map_py(&mut self, value : HashMap<u32,TrackerStripPulse>) {
302    self.pulse_map = value;
303  }
304  
305  #[getter]
306  #[pyo3(name="gain_map")]
307  fn get_gain_map_py(&self) -> HashMap<u32, TrackerStripGain> {
308    self.gain_map.clone()
309  }
310  #[setter]
311  #[pyo3(name="gain_map")]
312  fn set_gain_map_py(&mut self, value : HashMap<u32,TrackerStripGain>) {
313    self.gain_map = value;
314  }
315
316  #[getter]
317  #[pyo3(name="adc_sig_cut")]
318  fn get_adc_sig_cut_py(&self) -> HashMap<u32,f32> {
319    self.adc_sig_cut.clone()
320  }
321  #[setter]
322  #[pyo3(name="adc_sig_cut")]
323  fn set_adc_sig_cut_py(&mut self, value : HashMap<u32,f32>) {
324    self.adc_sig_cut = value;
325  }
326}
327
328#[cfg(feature="pybindings")]
329pythonize!(TrackerOfflineCalibration);
330
gondola_core/calibration/tracker.rs

gondola_core/calibration/
tracker.rs
