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

gondola_core/calibration/
tracker.rs
