Line data Source code
1 : /**
2 : * @file TAMakerMichelElectronAlgorithm.cpp
3 : *
4 : * This is part of the DUNE DAQ Application Framework, copyright 2021.
5 : * Licensing/copyright details are in the COPYING file that you should have
6 : * received with this code.
7 : */
8 :
9 : #include "triggeralgs/MichelElectron/TAMakerMichelElectronAlgorithm.hpp"
10 : #include "TRACE/trace.h"
11 : #define TRACE_NAME "TAMakerMichelElectronAlgorithm"
12 : #include <vector>
13 : #include <algorithm>
14 :
15 : using namespace triggeralgs;
16 :
17 : using Logging::TLVL_DEBUG_MEDIUM;
18 :
19 : void
20 0 : TAMakerMichelElectronAlgorithm::process(const TriggerPrimitive& input_tp,
21 : std::vector<TriggerActivity>& output_ta)
22 : {
23 :
24 : // The first time operator() is called, reset the window object.
25 0 : if (m_current_window.is_empty()) {
26 0 : m_current_window.reset(input_tp);
27 0 : m_primitive_count++;
28 0 : return;
29 : }
30 :
31 : // If the difference between the current TP's start time and the start of the window
32 : // is less than the specified window size, add the TP to the window.
33 0 : if ((input_tp.time_start - m_current_window.time_start) < m_window_length) {
34 0 : m_current_window.add(input_tp);
35 : }
36 :
37 : // Check Michel Candidate ========================================================
38 : // We've filled the window, now require a sufficient length track AND that the track
39 : // has a potential Bragg P, and then a kink.
40 0 : else if (longest_activity().size() > m_adjacency_threshold) {
41 :
42 :
43 : // We have a good length acitivity, now search for Bragg peak and kinks
44 0 : std::vector<TriggerPrimitive> trackHits = longest_activity();
45 :
46 0 : if (check_bragg_peak(trackHits)){
47 0 : if (check_kinks(trackHits)){
48 0 : TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TAM:ME] Emitting a trigger for candidate Michel event.";
49 0 : output_ta.push_back(construct_ta());
50 0 : m_current_window.reset(input_tp);
51 : } // Kinks
52 : } // Bragg peak
53 :
54 0 : }
55 :
56 : // Otherwise, slide the window along using the current TP.
57 : else {
58 0 : m_current_window.move(input_tp, m_window_length);
59 : }
60 :
61 0 : m_primitive_count++;
62 0 : return;
63 : }
64 :
65 : void
66 0 : TAMakerMichelElectronAlgorithm::configure(const nlohmann::json& config)
67 : {
68 0 : TriggerActivityMaker::configure(config);
69 :
70 : // FIX ME: Use some schema here. Also can't work out how to pass booleans.
71 0 : if (config.is_object()) {
72 0 : if (config.contains("window_length"))
73 0 : m_window_length = config["window_length"];
74 0 : if (config.contains("adjacency_tolerance"))
75 0 : m_adj_tolerance = config["adjacency_tolerance"];
76 0 : if (config.contains("adjacency_threshold"))
77 0 : m_adjacency_threshold = config["adjacency_threshold"];
78 : }
79 :
80 0 : }
81 :
82 : TriggerActivity
83 0 : TAMakerMichelElectronAlgorithm::construct_ta() const
84 : {
85 :
86 0 : TriggerPrimitive latest_tp_in_window = m_current_window.inputs.back();
87 :
88 0 : TriggerActivity ta;
89 0 : ta.time_start = m_current_window.time_start;
90 0 : ta.time_end = latest_tp_in_window.time_start + latest_tp_in_window.samples_over_threshold * 32; // FIXME: Replace the hard-coded SOT to TOT scaling.
91 0 : ta.time_peak = latest_tp_in_window.samples_to_peak * 32 + latest_tp_in_window.time_start; // FIXME: Replace STP to `time_peak` conversion.
92 0 : ta.time_activity = ta.time_peak;
93 0 : ta.channel_start = latest_tp_in_window.channel;
94 0 : ta.channel_end = latest_tp_in_window.channel;
95 0 : ta.channel_peak = latest_tp_in_window.channel;
96 0 : ta.adc_integral = m_current_window.adc_integral;
97 0 : ta.adc_peak = latest_tp_in_window.adc_peak;
98 0 : ta.detid = latest_tp_in_window.detid;
99 0 : ta.type = TriggerActivity::Type::kTPC;
100 0 : ta.algorithm = TriggerActivity::Algorithm::kMichelElectron;
101 0 : ta.inputs = m_current_window.inputs;
102 :
103 0 : return ta;
104 0 : }
105 :
106 : std::vector<TriggerPrimitive>
107 0 : TAMakerMichelElectronAlgorithm::longest_activity() const
108 : {
109 : // This function attempts to return a vector of hits that correspond to the longest
110 : // piece of activity in the current window. The logic follows that from the HMA
111 : // check_adjacency() function and further details can be found there.
112 0 : std::vector<TriggerPrimitive> trackHits;
113 0 : std::vector<TriggerPrimitive> finalHits; // The vector of track hits, which we return
114 :
115 0 : uint16_t adj = 1; // Initialise adjacency, 1 for the first wire.
116 0 : uint16_t max = 0;
117 0 : unsigned int channel = 0; // Current channel ID
118 0 : unsigned int next_channel = 0; // Next channel ID
119 0 : unsigned int next = 0; // The next position in the hit channels vector
120 0 : unsigned int tol_count = 0; // Tolerance count, should not pass adj_tolerance
121 :
122 : // Generate a channelID ordered list of hit channels for this window
123 0 : std::vector<TriggerPrimitive> hitList;
124 0 : for (auto tp : m_current_window.inputs) {
125 0 : hitList.push_back(tp);
126 : }
127 0 : std::sort(hitList.begin(), hitList.end(), [](TriggerPrimitive a, TriggerPrimitive b)
128 0 : { return a.channel < b.channel; });
129 :
130 : // ADAJACENCY LOGIC ====================================================================
131 : // =====================================================================================
132 : // Adjcancency Tolerance = Number of times prepared to skip missed hits before resetting
133 : // the adjacency count. This accounts for things like dead channels / missed TPs. The
134 : // maximum gap is 4 which comes from tuning on December 2021 coldbox data, and June 2022
135 : // coldbox runs.
136 0 : for (int i = 0; i < hitList.size(); ++i) {
137 :
138 0 : next = (i + 1) % hitList.size(); // Loops back when outside of channel list range
139 0 : channel = hitList.at(i).channel;
140 0 : next_channel = hitList.at(next).channel; // Next channel with a hit
141 :
142 0 : if (trackHits.size() == 0 ){ trackHits.push_back(hitList.at(i)); }
143 :
144 : // End of vector condition.
145 0 : if (next_channel == 0) { next_channel = channel - 1; }
146 :
147 : // Skip same channel hits for adjacency counting, but add to the track!
148 0 : if (next_channel == channel) {
149 0 : trackHits.push_back(hitList.at(next));
150 0 : continue; }
151 :
152 : // If next hit is on next channel, increment the adjacency count.
153 0 : else if (next_channel == channel + 1){
154 0 : trackHits.push_back(hitList.at(next));
155 0 : ++adj; }
156 :
157 : // If next channel is not on the next hit, but the 'second next', increase adjacency
158 : // but also tally up with the tolerance counter.
159 0 : else if (((next_channel == channel + 2) || (next_channel == channel + 3) ||
160 0 : (next_channel == channel + 4) || (next_channel == channel + 5))
161 0 : && (tol_count < m_adj_tolerance)) {
162 0 : trackHits.push_back(hitList.at(next));
163 0 : ++adj;
164 0 : for (int i = 0 ; i < next_channel-channel ; ++i){ ++tol_count; }
165 : }
166 :
167 : // If next hit isn't within reach, end the adjacency count and check for a new max.
168 : // Reset variables for next iteration.
169 : else {
170 0 : if (adj > max) {
171 0 : max = adj;
172 0 : finalHits.clear(); // Clear previous track
173 0 : for (auto h : trackHits){
174 0 : finalHits.push_back(h);
175 : }
176 : }
177 0 : adj = 1;
178 0 : tol_count = 0;
179 0 : trackHits.clear();
180 : }
181 : }
182 :
183 0 : return finalHits;
184 0 : }
185 :
186 :
187 : // Function that tries to identify a Bragg peak via a running
188 : // mean average of the ADC values. We use the running mean as it's less susceptible to
189 : // spikes of activity that might trick the algorithm. We establish a baseline, then
190 : // count up clusters of charge deposition above that baseline. If the largest is at
191 : // one of the ends of that collection, signal a potential Bragg peak.
192 : bool
193 0 : TAMakerMichelElectronAlgorithm::check_bragg_peak(std::vector<TriggerPrimitive> trackHits)
194 : {
195 0 : bool bragg = false;
196 0 : std::vector<float> adc_means_list;
197 0 : uint16_t convolve_value = 6;
198 :
199 : // Loop over hits that correspond to high adjacency activity
200 0 : for (uint16_t i = 0; i < trackHits.size(); ++i){
201 0 : float adc_sum = 0;
202 : float adc_mean = 0;
203 :
204 : // Calculate running ADC mean of this track
205 0 : for (uint16_t j = i; j < i+convolve_value; ++j){
206 0 : int hit = (j) % trackHits.size();
207 0 : adc_sum += trackHits.at(hit).adc_integral;
208 : }
209 :
210 0 : adc_mean = adc_sum / convolve_value;
211 0 : adc_means_list.push_back(adc_mean);
212 0 : adc_sum = 0;
213 : }
214 :
215 : // We now have a list of convolved adc means.
216 0 : float ped = std::accumulate(adc_means_list.begin(), adc_means_list.end(), 0.0) / adc_means_list.size();
217 0 : float charge = 0;
218 0 : std::vector<float> charge_dumps;
219 :
220 : // Now go through the list, picking up clusters of charge above the baseline/ped
221 0 : for (auto a : adc_means_list){
222 0 : if (a > ped){
223 0 : charge += a;
224 : }
225 0 : else if( a < ped && charge !=0 ){
226 0 : charge_dumps.push_back(charge);
227 0 : charge = 0;
228 : }
229 : }
230 :
231 : // If the maximum of that list of charge dumps is near(at?) either end of it
232 0 : float max_charge = *max_element(charge_dumps.begin(), charge_dumps.end());
233 0 : if(max_charge == charge_dumps.front() || max_charge == charge_dumps.back()){ bragg=true; }
234 :
235 0 : return bragg;
236 0 : }
237 :
238 : bool
239 0 : TAMakerMichelElectronAlgorithm::check_kinks(std::vector<TriggerPrimitive> finalHits)
240 : {
241 0 : bool kinks = false; // We actually required two kinks in the coldbox, the michel kink and the wes kink
242 0 : std::vector<float> runningGradient;
243 0 : std::vector<float> runningMeanGradient;
244 :
245 : // Choice to be made here. Do we want to scane in collection (z) or time (x) direction when calculating gradient between hits. I
246 : // would say if we have already made the request to pass a track of length specific threshold which is longer than the drift
247 : // direction for the coldbox, it makes sense to scan across channels a little more.
248 0 : std::sort(finalHits.begin(), finalHits.end(), [](TriggerPrimitive a, TriggerPrimitive b) { return a.channel < b.channel; });
249 :
250 : // Populate the runningGradient with the track hits. Do this between ith and i+kth TPs, to small scale fluctuations of the track
251 : // Yet k should be kept small, so that enough gradient information is preserved at the end of the track to identify kinks
252 0 : for (int i=0 ; i < finalHits.size()-2; i++){
253 :
254 : // Skip same channel hits or if the start times are the same - no div by zero!
255 0 : if (finalHits.at(i+2).channel == finalHits.at(i).channel || (finalHits.at(i+2).time_start == finalHits.at(i).time_start) ) { continue; }
256 :
257 : // Check that the next TP is closeby; enough in space and time directions so as to avoid obtaining a gradient value from
258 : // same channel hits at large time difference or vice versa due to kink topology or showers. Clearly we shouldn't be very far in
259 : // channel number, but since we might later try to do this check in the time direction, leave the condition in.
260 0 : int diff = finalHits.at(i+2).time_start - finalHits.at(i).time_start;
261 0 : if((std::abs(diff) > 1000) || ((std::abs(channel_diff_t(finalHits.at(i+2).channel) - channel_diff_t(finalHits.at(i).channel)) > 6))) { continue; }
262 :
263 : // Gradient is just change in z (collection) over change in x (drift). x is admitedly roughly converted from
264 : // hit start time, but I don't think diffusion effects are a huge concern over 20cm. Using mm for readability/visualisation
265 0 : float dz = (finalHits.at(i+2).channel - finalHits.at(i).channel)*4.67; // Change in collection wire z to separation in mm
266 0 : long long int dt = finalHits.at(i+2).time_start - finalHits.at(i).time_start;
267 0 : float dx = dt*0.028; // Change time to separation in x mm
268 0 : float g = dz/dx;
269 :
270 0 : runningGradient.push_back(g);
271 : }
272 :
273 : // Require a decent length of the gradients vector. Otherwise some adjacent events are showers and the conditions above mean
274 : // we don't get enough entries. In essence, this provides some confidence that it's track-like rather than shower-like. Which
275 : // is what we want for a Michel event.
276 0 : if ( runningGradient.size() > 10 ){
277 :
278 : // Now lets take a running mean of the gradients between TPs, less susceptible to wild changes due to deltas/etc
279 0 : for(int g=0 ; g < runningGradient.size()-1 ; g++){
280 0 : float gsum = runningGradient.at(g) + runningGradient.at(g+1);
281 0 : runningMeanGradient.push_back(gsum/2);
282 : }
283 :
284 : // We have a list of gradients, now just demand that the two ends have gradients that differ significantly
285 : // from the mean gradient of the activity at both ends. This aims to pick out wesKinks and michelKinks
286 0 : if(runningMeanGradient.size() > 10 ){
287 :
288 0 : float mean = (std::abs(std::accumulate(runningMeanGradient.begin(), runningMeanGradient.end(), 0.0)))/(runningMeanGradient.size());
289 :
290 : // If you're testing on simulation or december data, you won't see the wes kink, so use an || instead of &&
291 0 : if((std::abs(runningMeanGradient.front()) + mean > 2.5*mean) || ((std::abs(runningMeanGradient.back() + mean)) > 2.5*mean)){ kinks=true; }
292 : }
293 : }
294 :
295 0 : return kinks;
296 0 : }
297 :
298 : // ===============================================================================================
299 : // ===============================================================================================
300 : // Functions below this line are for debugging purposes.
301 : // ===============================================================================================
302 :
303 : void
304 0 : TAMakerMichelElectronAlgorithm::add_window_to_record(Window window)
305 : {
306 0 : m_window_record.push_back(window);
307 0 : return;
308 : }
309 :
310 :
311 : // Function to dump the details of the TA window currently on record
312 : void
313 0 : TAMakerMichelElectronAlgorithm::dump_window_record()
314 : {
315 : // FIX ME: Need to index this outfile in the name by detid or something similar.
316 0 : std::ofstream outfile;
317 0 : outfile.open("window_record_tam.csv", std::ios_base::app);
318 :
319 0 : for (auto window : m_window_record) {
320 0 : outfile << window.time_start << ",";
321 0 : outfile << window.inputs.back().time_start << ",";
322 0 : outfile << window.inputs.back().time_start - window.time_start << ","; // window_length - from TP start times
323 0 : outfile << window.adc_integral << ",";
324 0 : outfile << window.n_channels_hit() << ","; // Number of unique channels with hits
325 0 : outfile << window.inputs.size() << ","; // Number of TPs in Window
326 0 : outfile << window.inputs.back().channel << ","; // Last TP Channel ID
327 0 : outfile << window.inputs.front().channel << ","; // First TP Channel ID
328 0 : outfile << longest_activity().size() << std::endl; // New adjacency value for the window
329 0 : }
330 :
331 0 : outfile.close();
332 :
333 0 : m_window_record.clear();
334 :
335 0 : return;
336 0 : }
337 :
338 : // Function to add current TP details to a text file for testing and debugging.
339 : void
340 0 : TAMakerMichelElectronAlgorithm::dump_tp(TriggerPrimitive const& input_tp)
341 : {
342 0 : std::ofstream outfile;
343 0 : outfile.open("coldbox_tps.txt", std::ios_base::app);
344 :
345 : // Output relevant TP information to file
346 0 : outfile << input_tp.time_start << " "; // Start time of TP
347 0 : outfile << input_tp.samples_over_threshold << " "; // in multiples of 25
348 0 : outfile << input_tp.samples_to_peak << " "; //
349 0 : outfile << input_tp.channel << " "; // Offline channel ID
350 0 : outfile << input_tp.adc_integral << " "; // ADC Sum
351 0 : outfile << input_tp.adc_peak << " "; // ADC Peak Value
352 0 : outfile << input_tp.detid << " "; // Det ID - Identifies detector element, APA or PDS part etc...
353 0 : outfile.close();
354 :
355 0 : return;
356 0 : }
357 :
358 : /*
359 : void
360 : TAMakerMichelElectronAlgorithm::flush(timestamp_t, std::vector<TriggerActivity>& output_ta)
361 : {
362 : // Check the status of the current window, construct TA if conditions are met. Regardless
363 : // of whether the conditions are met, reset the window.
364 : if(m_current_window.adc_integral > m_adc_threshold && m_trigger_on_adc){
365 : //else if(m_current_window.adc_integral > m_conf.adc_threshold && m_conf.trigger_on_adc){
366 : //TLOG_DEBUG(TRACE_NAME) << "ADC integral in window is greater than specified threshold.";
367 : output_ta.push_back(construct_ta());
368 : }
369 : else if(m_current_window.n_channels_hit() > m_n_channels_threshold && m_trigger_on_n_channels){
370 : //else if(m_current_window.n_channels_hit() > m_conf.n_channels_threshold && m_conf.trigger_on_n_channels){
371 : //TLOG_DEBUG(TRACE_NAME) << "Number of channels hit in the window is greater than specified threshold.";
372 : output_ta.push_back(construct_ta());
373 : }
374 :
375 : //TLOG_DEBUG(TRACE_NAME) << "Clearing the current window, on the arrival of the next input_tp, the window will be
376 : reset."; m_current_window.clear();
377 :
378 : return;
379 : }*/
380 :
381 : // Register algo in TA Factory
382 12 : REGISTER_TRIGGER_ACTIVITY_MAKER(TRACE_NAME, TAMakerMichelElectronAlgorithm)
|
