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")]
26pub fn waveform_analysis(event : &mut RBEvent,
27 rb : &ReadoutBoard,
28 settings : AnalysisEngineSettings)
29-> Result<(), AnalysisError> {
30 // Don't do analysis for mangled events!
31 if event.has_any_mangling_flag() {
32 warn!("Event for RB {} has data mangling! Not doing analysis!", rb.rb_id);
33 return Err(AnalysisError::DataMangling);
34 }
35 match event.self_check() {
36 Err(_err) => {
37 // Phlip want to ahve all hits even if they are broken
38 },
39 Ok(_) => ()
40 }
41 let active_channels = event.header.get_channels();
42 // will become a parameter
43 let fit_sinus = true;
44 // allocate memory for the calbration results
45 let mut voltages : Vec<f32>= vec![0.0; NWORDS];
46 let mut times : Vec<f32>= vec![0.0; NWORDS];
47
48 // Step 0 : If desired, fit sine
49 let mut fit_result = (0.0f32, 0.0f32, 0.0f32);
50 if fit_sinus {
51 if !active_channels.contains(&8) {
52 warn!("RB {} does not have ch9 data!", rb.rb_id);
53 //println!("{}", event.header);
54 return Err(AnalysisError::NoChannel9);
55 }
56 rb.calibration.voltages(9,
57 event.header.stop_cell as usize,
58 &event.adc[8],
59 &mut voltages);
60 //warn!("We have to rework the spike cleaning!");
61 //match RBCalibrations::spike_cleaning(&mut ch_voltages,
62 // event.header.stop_cell) {
63 // Err(err) => {
64 // error!("Spike cleaning failed! {err}");
65 // }
66 // Ok(_) => ()
67 //}
68 rb.calibration.nanoseconds(9,
69 event.header.stop_cell as usize,
70 &mut times);
71 fit_result = fit_sine_simple(&voltages, ×);
72
73 //println!("FIT RESULT = {:?}", fit_result);
74 //event.header.set_sine_fit(fit_result);
75 }
76
77 // structure to store final result
78 // extend with Vec<TofHit> in case
79 // we want to have multiple hits
80 let mut paddles = HashMap::<u8, TofHit>::new();
81 //println!("RBID {}, Paddles {:?}", rb.rb_id ,rb.get_paddle_ids());
82 for pid in rb.get_paddle_ids() {
83 // cant' fail by constructon of pid
84 let ch_a = rb.get_pid_rbchA(pid).unwrap() as usize;
85 let ch_b = rb.get_pid_rbchB(pid).unwrap() as usize;
86 let mut hit = TofHit::new();
87 hit.paddle_id = pid;
88 //println!("{ch_a}, {ch_b}, active_channels {:?}", active_channels);
89 for (k, ch) in [ch_a, ch_b].iter().enumerate() {
90 // Step 1: Calibration
91 //println!("Ch {}, event {}", ch, event);
92 //println!("---------------------------");
93 //println!("pid {}, active channels : {:?}, ch {}",pid, active_channels, ch);
94 if !active_channels.contains(&(*ch as u8 -1)) {
95 trace!("Skipping channel {} because it is not marked to be readout in the event header channel mask!", ch);
96 continue;
97 }
98 //println!("Will do waveform analysis for ch {}", ch);
99 rb.calibration.voltages(*ch,
100 event.header.stop_cell as usize,
101 &event.adc[*ch as usize -1],
102 &mut voltages);
103 //FIXME - spike cleaning!
104 //match RBCalibrations::spike_cleaning(&mut ch_voltages,
105 // event.header.stop_cell) {
106 // Err(err) => {
107 // error!("Spike cleaning failed! {err}");
108 // }
109 // Ok(_) => ()
110 //}
111 rb.calibration.nanoseconds(*ch,
112 event.header.stop_cell as usize,
113 &mut times);
114 // Step 2: Pedestal subtraction
115 let (ped, ped_err) = calculate_pedestal(&voltages,
116 settings.pedestal_thresh,
117 settings.pedestal_begin_bin,
118 settings.pedestal_win_bins);
119 trace!("Calculated pedestal of {} +- {}", ped, ped_err);
120 for n in 0..voltages.len() {
121 voltages[n] -= ped;
122 }
123 let mut charge : f32 = 0.0;
124 //let peaks : Vec::<(usize, usize)>;
125 let mut cfd_times = Vec::<f32>::new();
126 let mut max_volts = 0.0f32;
127 // Step 4 : Find peaks
128 // FIXME - what do we do for multiple peaks?
129 // Currently we basically throw them away
130 match find_peaks(&voltages ,
131 × ,
132 settings.find_pks_t_start ,
133 settings.find_pks_t_window,
134 settings.min_peak_size ,
135 settings.find_pks_thresh ,
136 settings.max_peaks ) {
137 Err(err) => {
138 // FIXME - if this happens, most likely the channel is dead.
139 debug!("Unable to find peaks for RB{:02} ch {ch}! Ignoring this channel!", rb.rb_id);
140 debug!("We won't be able to calculate timing information for this channel! Err {err}");
141 },
142 Ok(peaks) => {
143 //peaks = pks;
144 // Step 5 : Find tdcs
145 //println!("Found {} peaks for ch {}! {:?}", peaks.len(), raw_ch, peaks);
146 for pk in peaks.iter() {
147 match cfd_simple(&voltages,
148 ×,
149 settings.cfd_fraction,
150 pk.0, pk.1) {
151 Err(err) => {
152 debug!("Unable to calculate cfd for peak {} {}! {}", pk.0, pk.1, err);
153 }
154 Ok(cfd) => {
155 cfd_times.push(cfd);
156 }
157 }
158 let pk_height = voltages[pk.0..pk.1].iter().max_by(|a,b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Less)).unwrap();
159 max_volts = *pk_height;
160 let max_index = voltages.iter().position(|element| *element == max_volts).unwrap();
161
162 let (start_q_int, stop_q_int) = if max_index - 40 < 10 {
163 (10, 210)
164 } else {
165 (max_index - 40, max_index + 160)
166 };
167
168
169 //debug!("Check impedance value! Just using 50 [Ohm]");
170 // Step 3 : charge integration
171 // FIXME - make impedance a settings parameter
172 match integrate(&voltages,
173 ×,
174 //settings.integration_start,
175 //settings.integration_window,
176 //pk.0,
177 //pk.1,
178 start_q_int,
179 stop_q_int,
180 50.0) {
181 Err(err) => {
182 error!("Integration failed! Err {err}");
183 }
184 Ok(chrg) => {
185 charge = chrg;
186 }
187 }
188 // // just do the first peak for now
189 // let pk_height = voltages[pk.0..pk.1].iter().max_by(|a,b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Less)).unwrap();
190 // max_volts = *pk_height;
191 // //debug!("Check impedance value! Just using 50 [Ohm]");
192 // // Step 3 : charge integration
193 // // FIXME - make impedance a settings parameter
194 // match integrate(&voltages,
195 // ×,
196 // //settings.integration_start,
197 // //settings.integration_window,
198 // pk.0,
199 // pk.1,
200 // 50.0) {
201 // Err(err) => {
202 // error!("Integration failed! Err {err}");
203 // }
204
205 break;
206 }
207 }// end OK
208 } // end match find_peaks
209 let mut tdc : f32 = 0.0;
210 if cfd_times.len() > 0 {
211 tdc = cfd_times[0];
212 }
213 //println!("Calucalated tdc {}, charge {}, max {} for ch {}!", tdc, charge, max_volts, ch);
214 //if rb.channel_to_paddle_end_id[*raw_ch as usize] > 2000 {
215 if k == 0 {
216 //hit.ftime_a = tdc;
217 //hit.fpeak_a = max_volts;
218 hit.set_time_a(tdc);
219 hit.set_charge_a(charge);
220 hit.set_peak_a(max_volts);
221 hit.baseline_a = f16::from_f32(ped);
222 hit.baseline_a_rms = f16::from_f32(ped_err);
223 } else {
224 //hit.ftime_b = tdc;
225 //hit.fpeak_b = max_volts;
226 hit.set_time_b(tdc);
227 hit.set_charge_b(charge);
228 hit.set_peak_b(max_volts);
229 hit.baseline_b = f16::from_f32(ped);
230 hit.baseline_b_rms = f16::from_f32(ped_err);
231 // this is the seoond iteration,
232 // we are done!
233 hit.phase = f16::from_f32(fit_result.2);
234 paddles.insert(pid, hit);
235 }
236 }
237 }
238 let result = paddles.into_values().collect();
239 event.hits = result;
240 //print ("EVENT {}", event);
241 Ok(())
242}
243