gondola_core/tof/analysis_engine.rs
1// This file is part of gaps-online-software and published
2// under the GPLv3 license
3
4use crate::prelude::*;
5
6/// Waveform analysis engine - identify waveform variables
7///
8/// This will populate the TofHits in an RBEvent
9///
10/// TofHits contain information about peak location,
11/// charge, timing.
12///
13/// FIXME - I think this should take a HashMap with
14/// algorithm settings, which we can load from a
15/// json file
16///
17/// # Arguments
18///
19/// * event : current RBEvent with waveforms to
20/// work on
21/// * rb : ReadoutBoard as loaded from the DB,
22/// with latest calibration attached
23/// * settings : Parameters to configure the waveform
24/// analysis & peak finding
25#[cfg(feature="database")]
26#[cfg_attr(feature="pybindings", pyfunction)]
27pub fn waveform_analysis(event : &mut RBEvent,
28 rb : &ReadoutBoard,
29 settings : AnalysisEngineSettings)
30-> Result<(), AnalysisError> {
31 // Don't do analysis for mangled events!
32 if event.has_any_mangling_flag() {
33 warn!("Event for RB {} has data mangling! Not doing analysis!", rb.rb_id);
34 return Err(AnalysisError::DataMangling);
35 }
36 match event.self_check() {
37 Err(_err) => {
38 // Phlip want to ahve all hits even if they are broken
39 },
40 Ok(_) => ()
41 }
42 let active_channels = event.header.get_channels();
43 // will become a parameter
44 let fit_sinus = true;
45 // allocate memory for the calbration results
46 let mut voltages : Vec<f32>= vec![0.0; NWORDS];
47 let mut times : Vec<f32>= vec![0.0; NWORDS];
48
49 // Step 0 : If desired, fit sine
50 let mut fit_result = (0.0f32, 0.0f32, 0.0f32);
51 if fit_sinus {
52 if !active_channels.contains(&8) {
53 warn!("RB {} does not have ch9 data!", rb.rb_id);
54 //println!("{}", event.header);
55 return Err(AnalysisError::NoChannel9);
56 }
57 rb.calibration.voltages(9,
58 event.header.stop_cell as usize,
59 &event.adc[8],
60 &mut voltages);
61 //warn!("We have to rework the spike cleaning!");
62 //match RBCalibrations::spike_cleaning(&mut ch_voltages,
63 // event.header.stop_cell) {
64 // Err(err) => {
65 // error!("Spike cleaning failed! {err}");
66 // }
67 // Ok(_) => ()
68 //}
69 rb.calibration.nanoseconds(9,
70 event.header.stop_cell as usize,
71 &mut times);
72 fit_result = fit_sine_simple(&voltages, ×);
73
74 //println!("FIT RESULT = {:?}", fit_result);
75 //event.header.set_sine_fit(fit_result);
76 }
77
78 // structure to store final result
79 // extend with Vec<TofHit> in case
80 // we want to have multiple hits
81 let mut paddles = HashMap::<u8, TofHit>::new();
82 //println!("RBID {}, Paddles {:?}", rb.rb_id ,rb.get_paddle_ids());
83 for pid in rb.get_paddle_ids() {
84 // cant' fail by constructon of pid
85 let ch_a = rb.get_pid_rbchA(pid).unwrap() as usize;
86 let ch_b = rb.get_pid_rbchB(pid).unwrap() as usize;
87 let mut hit = TofHit::new();
88 hit.paddle_id = pid;
89 //println!("{ch_a}, {ch_b}, active_channels {:?}", active_channels);
90 for (k, ch) in [ch_a, ch_b].iter().enumerate() {
91 // Step 1: Calibration
92 //println!("Ch {}, event {}", ch, event);
93 //println!("---------------------------");
94 //println!("pid {}, active channels : {:?}, ch {}",pid, active_channels, ch);
95 if !active_channels.contains(&(*ch as u8 -1)) {
96 trace!("Skipping channel {} because it is not marked to be readout in the event header channel mask!", ch);
97 continue;
98 }
99 //println!("Will do waveform analysis for ch {}", ch);
100 rb.calibration.voltages(*ch,
101 event.header.stop_cell as usize,
102 &event.adc[*ch as usize -1],
103 &mut voltages);
104 //FIXME - spike cleaning!
105 //match RBCalibrations::spike_cleaning(&mut ch_voltages,
106 // event.header.stop_cell) {
107 // Err(err) => {
108 // error!("Spike cleaning failed! {err}");
109 // }
110 // Ok(_) => ()
111 //}
112 rb.calibration.nanoseconds(*ch,
113 event.header.stop_cell as usize,
114 &mut times);
115 // Step 2: Pedestal subtraction
116 let (ped, ped_err) = calculate_pedestal(&voltages,
117 settings.pedestal_thresh,
118 settings.pedestal_begin_bin,
119 settings.pedestal_win_bins);
120 trace!("Calculated pedestal of {} +- {}", ped, ped_err);
121 for n in 0..voltages.len() {
122 voltages[n] -= ped;
123 }
124 let mut charge : f32 = 0.0;
125 //let peaks : Vec::<(usize, usize)>;
126 let mut cfd_times = Vec::<f32>::new();
127 let mut max_volts = 0.0f32;
128 // Step 4 : Find peaks
129 // FIXME - what do we do for multiple peaks?
130 // Currently we basically throw them away
131 match find_peaks(&voltages ,
132 × ,
133 settings.find_pks_t_start ,
134 settings.find_pks_t_window,
135 settings.min_peak_size ,
136 settings.find_pks_thresh ,
137 settings.max_peaks ) {
138 Err(err) => {
139 // FIXME - if this happens, most likely the channel is dead.
140 debug!("Unable to find peaks for RB{:02} ch {ch}! Ignoring this channel!", rb.rb_id);
141 debug!("We won't be able to calculate timing information for this channel! Err {err}");
142 },
143 Ok(peaks) => {
144 //peaks = pks;
145 // Step 5 : Find tdcs
146 //println!("Found {} peaks for ch {}! {:?}", peaks.len(), raw_ch, peaks);
147 for pk in peaks.iter() {
148 match cfd_simple(&voltages,
149 ×,
150 settings.cfd_fraction,
151 pk.0, pk.1) {
152 Err(err) => {
153 debug!("Unable to calculate cfd for peak {} {}! {}", pk.0, pk.1, err);
154 }
155 Ok(cfd) => {
156 cfd_times.push(cfd);
157 }
158 }
159 let pk_height = voltages[pk.0..pk.1].iter().max_by(|a,b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Less)).unwrap();
160 max_volts = *pk_height;
161 let max_index = voltages.iter().position(|element| *element == max_volts).unwrap();
162
163 let (start_q_int, stop_q_int) = if max_index - 40 < 10 {
164 (10, 210)
165 } else {
166 (max_index - 40, max_index + 160)
167 };
168
169
170 //debug!("Check impedance value! Just using 50 [Ohm]");
171 // Step 3 : charge integration
172 // FIXME - make impedance a settings parameter
173 match integrate(&voltages,
174 ×,
175 //settings.integration_start,
176 //settings.integration_window,
177 //pk.0,
178 //pk.1,
179 start_q_int,
180 stop_q_int,
181 50.0) {
182 Err(err) => {
183 error!("Integration failed! Err {err}");
184 }
185 Ok(chrg) => {
186 charge = chrg;
187 }
188 }
189 // // just do the first peak for now
190 // let pk_height = voltages[pk.0..pk.1].iter().max_by(|a,b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Less)).unwrap();
191 // max_volts = *pk_height;
192 // //debug!("Check impedance value! Just using 50 [Ohm]");
193 // // Step 3 : charge integration
194 // // FIXME - make impedance a settings parameter
195 // match integrate(&voltages,
196 // ×,
197 // //settings.integration_start,
198 // //settings.integration_window,
199 // pk.0,
200 // pk.1,
201 // 50.0) {
202 // Err(err) => {
203 // error!("Integration failed! Err {err}");
204 // }
205
206 break;
207 }
208 }// end OK
209 } // end match find_peaks
210 let mut tdc : f32 = 0.0;
211 if cfd_times.len() > 0 {
212 tdc = cfd_times[0];
213 }
214 //println!("Calucalated tdc {}, charge {}, max {} for ch {}!", tdc, charge, max_volts, ch);
215 if k == 0 {
216 hit.set_time_a(tdc);
217 hit.set_charge_a(charge);
218 hit.set_peak_a(max_volts);
219 hit.baseline_a = f16::from_f32(ped);
220 hit.baseline_a_rms = f16::from_f32(ped_err);
221 // calculate time over threshold
222 match settings.tot_threshold_low {
223 Some(thr) => {
224 let th_low_sl = time_over_threshold(&voltages, ×, thr);
225 hit.tot_low_a = f16::from_f32(th_low_sl.0);
226 hit.tot_slp_low_a = f16::from_f32(th_low_sl.1);
227 }
228 None => ()
229 }
230 match settings.tot_threshold_high {
231 Some(thr) => {
232 let th_high_sl = time_over_threshold(&voltages, ×, thr);
233 hit.tot_high_a = f16::from_f32(th_high_sl.0);
234 hit.tot_slp_high_a = f16::from_f32(th_high_sl.1);
235 }
236 None => ()
237 }
238 } else {
239
240 hit.set_time_b(tdc);
241 hit.set_charge_b(charge);
242 hit.set_peak_b(max_volts);
243 hit.baseline_b = f16::from_f32(ped);
244 hit.baseline_b_rms = f16::from_f32(ped_err);
245 // this is the seoond iteration,
246 // we are done!
247 hit.phase = f16::from_f32(fit_result.2);
248 // calculate time over threshold
249 match settings.tot_threshold_low {
250 Some(thr) => {
251 let th_low_sl = time_over_threshold(&voltages, ×, thr);
252 hit.tot_low_b = f16::from_f32(th_low_sl.0);
253 hit.tot_slp_low_b = f16::from_f32(th_low_sl.1);
254 }
255 None => ()
256 }
257 match settings.tot_threshold_high {
258 Some(thr) => {
259 let th_high_sl = time_over_threshold(&voltages, ×, thr);
260 hit.tot_high_b = f16::from_f32(th_high_sl.0);
261 hit.tot_slp_high_b = f16::from_f32(th_high_sl.1);
262 }
263 None => ()
264 }
265 paddles.insert(pid, hit);
266 }
267 }
268 }
269 let result = paddles.into_values().collect();
270 event.hits = result;
271 //print ("EVENT {}", event);
272 Ok(())
273}
274