triggeralgs::TAMakerMichelElectronAlgorithm Class Reference

#include <TAMakerMichelElectronAlgorithm.hpp>

Inheritance diagram for triggeralgs::TAMakerMichelElectronAlgorithm:

Collaboration diagram for triggeralgs::TAMakerMichelElectronAlgorithm:

Classes
class	Window

Public Member Functions
void	process (const TriggerPrimitive &input_tp, std::vector< TriggerActivity > &output_ta)
	TP processing function that creates & fills TAs.
void	configure (const nlohmann::json &config)
Public Member Functions inherited from triggeralgs::TriggerActivityMaker
virtual	~TriggerActivityMaker ()=default
void	operator() (const TriggerPrimitive &input_tp, std::vector< TriggerActivity > &output_ta)
virtual bool	preprocess (const TriggerPrimitive &input_tp)
	TP pre-processing/filtering.
virtual void	postprocess (std::vector< TriggerActivity > &output_ta)
	Post-processing/filtering of the TAs, e.g. prescale.
virtual void	flush (timestamp_t, std::vector< TriggerActivity > &)

Private Member Functions
TriggerActivity	construct_ta () const
std::vector< TriggerPrimitive >	longest_activity () const
bool	check_bragg_peak (std::vector< TriggerPrimitive > trackHits)
bool	check_kinks (std::vector< TriggerPrimitive > trackHits)
void	add_window_to_record (Window window)
void	dump_window_record ()
void	dump_tp (TriggerPrimitive const &input_tp)

Private Attributes
Window	m_current_window
uint64_t	m_primitive_count = 0
bool	debug = false
uint16_t	m_adjacency_threshold = 15
int	m_max_adjacency = 0
uint16_t	m_adj_tolerance = 3
int	index = 0
uint16_t	ta_adc = 0
uint16_t	ta_channels = 0
timestamp_t	m_window_length = 50000
std::vector< Window >	m_window_record

Additional Inherited Members
Public Attributes inherited from triggeralgs::TriggerActivityMaker
std::atomic< uint64_t >	m_data_vs_system_time = 0
std::atomic< uint64_t >	m_initial_offset = 0
uint64_t	m_prescale = 1
	Configurable prescale factor.
uint64_t	m_ta_count = 0
	TA made count for prescaling.
uint32_t	m_max_samples_over_threshold = std::numeric_limits<uint32_t>::max()
	Time-over-threshold TP filtering.

Detailed Description

Definition at line 17 of file TAMakerMichelElectronAlgorithm.hpp.

Member Function Documentation

add_window_to_record()

void TAMakerMichelElectronAlgorithm::add_window_to_record ( Window window )

private

Definition at line 304 of file TAMakerMichelElectronAlgorithm.cpp.

{
  m_window_record.push_back(window);
  return;
}

check_bragg_peak()

bool TAMakerMichelElectronAlgorithm::check_bragg_peak ( std::vector< TriggerPrimitive > trackHits )

private

Definition at line 193 of file TAMakerMichelElectronAlgorithm.cpp.

{
  bool bragg = false; 
  std::vector<float> adc_means_list;
  uint16_t convolve_value = 6;
 
  // Loop over hits that correspond to high adjacency activity
  for (uint16_t i = 0; i < trackHits.size(); ++i){
    float adc_sum = 0;
    float adc_mean = 0;
 
    // Calculate running ADC mean of this track 
    for (uint16_t j = i; j < i+convolve_value; ++j){
       int hit = (j) % trackHits.size(); 
       adc_sum += trackHits.at(hit).adc_integral;
    }
 
    adc_mean = adc_sum / convolve_value;
    adc_means_list.push_back(adc_mean);
    adc_sum = 0;
  } 
 
  // We now have a list of convolved adc means. 
  float ped = std::accumulate(adc_means_list.begin(), adc_means_list.end(), 0.0) / adc_means_list.size();
  float charge = 0;
  std::vector<float> charge_dumps;
 
  // Now go through the list, picking up clusters of charge above the baseline/ped
  for (auto a : adc_means_list){
    if (a > ped){
       charge += a;
    }
    else if( a < ped && charge !=0 ){
     charge_dumps.push_back(charge);
     charge = 0; 
    } 
  } 
 
  // If the maximum of that list of charge dumps is near(at?) either end of it
  float max_charge = *max_element(charge_dumps.begin(), charge_dumps.end()); 
  if(max_charge == charge_dumps.front() || max_charge == charge_dumps.back()){ bragg=true; }
 
  return bragg;
 }

check_kinks()

bool TAMakerMichelElectronAlgorithm::check_kinks ( std::vector< TriggerPrimitive > trackHits )

private

Definition at line 239 of file TAMakerMichelElectronAlgorithm.cpp.

{
    bool kinks = false;  // We actually required two kinks in the coldbox, the michel kink and the wes kink
    std::vector<float> runningGradient;
    std::vector<float> runningMeanGradient;
 
    // Choice to be made here. Do we want to scane in collection (z) or time (x) direction when calculating gradient between hits. I
    // would say if we have already made the request to pass a track of length specific threshold which is longer than the drift
    // direction for the coldbox, it makes sense to scan across channels a little more.
    std::sort(finalHits.begin(), finalHits.end(), [](TriggerPrimitive a, TriggerPrimitive b) { return a.channel < b.channel; });
 
    // Populate the runningGradient with the track hits. Do this between ith and i+kth TPs, to small scale fluctuations of the track
    // Yet k should be kept small, so that enough gradient information is preserved at the end of the track to identify kinks
    for (int i=0 ; i < finalHits.size()-2; i++){
   
      // Skip same channel hits or if the start times are the same - no div by zero! 
      if (finalHits.at(i+2).channel == finalHits.at(i).channel || (finalHits.at(i+2).time_start == finalHits.at(i).time_start) ) { continue; }
 
      // Check that the next TP is closeby; enough in space and time directions so as to avoid obtaining a gradient value from
      // same channel hits at large time difference or vice versa due to kink topology or showers. Clearly we shouldn't be very far in 
      // channel number, but since we might later try to do this check in the time direction, leave the condition in.
      int diff = finalHits.at(i+2).time_start - finalHits.at(i).time_start;
      if((std::abs(diff) > 1000) || ((std::abs(channel_diff_t(finalHits.at(i+2).channel) - channel_diff_t(finalHits.at(i).channel)) > 6))) { continue; } 
 
      // Gradient is just change in z (collection) over change in x (drift). x is admitedly roughly converted from
      // hit start time, but I don't think diffusion effects are a huge concern over 20cm. Using mm for readability/visualisation 
      float dz = (finalHits.at(i+2).channel - finalHits.at(i).channel)*4.67; // Change in collection wire z to separation in mm
      long long int dt = finalHits.at(i+2).time_start - finalHits.at(i).time_start;
      float dx = dt*0.028; // Change time to separation in x mm
      float g = dz/dx;
 
      runningGradient.push_back(g); 
    }
 
    // Require a decent length of the gradients vector. Otherwise some adjacent events are showers and the conditions above mean
    // we don't get enough entries. In essence, this provides some confidence that it's track-like rather than shower-like. Which
    // is what we want for a Michel event.
    if ( runningGradient.size() > 10 ){
      
      // Now lets take a running mean of the gradients between TPs, less susceptible to wild changes due to deltas/etc
      for(int g=0 ; g < runningGradient.size()-1 ; g++){
        float gsum = runningGradient.at(g) + runningGradient.at(g+1);
        runningMeanGradient.push_back(gsum/2);
      } 
   
      // We have a list of gradients, now just demand that the two ends have gradients that differ significantly
      // from the mean gradient of the activity at both ends. This aims to pick out wesKinks and michelKinks
      if(runningMeanGradient.size() > 10 ){
 
        float mean = (std::abs(std::accumulate(runningMeanGradient.begin(), runningMeanGradient.end(), 0.0)))/(runningMeanGradient.size());
 
        // If you're testing on simulation or december data, you won't see the wes kink, so use an || instead of &&
        if((std::abs(runningMeanGradient.front()) + mean > 2.5*mean)  ||  ((std::abs(runningMeanGradient.back() + mean)) > 2.5*mean)){ kinks=true; }   
      } 
    } 
 
  return kinks;
}

configure()

void TAMakerMichelElectronAlgorithm::configure ( const nlohmann::json & config )

virtual

Reimplemented from triggeralgs::TriggerActivityMaker.

Definition at line 66 of file TAMakerMichelElectronAlgorithm.cpp.

{
  TriggerActivityMaker::configure(config);
 
  // FIX ME: Use some schema here. Also can't work out how to pass booleans.
  if (config.is_object()) {
    if (config.contains("window_length"))
      m_window_length = config["window_length"];
    if (config.contains("adjacency_tolerance"))
      m_adj_tolerance = config["adjacency_tolerance"];
    if (config.contains("adjacency_threshold"))
      m_adjacency_threshold = config["adjacency_threshold"];
  }
 
}

construct_ta()

TriggerActivity TAMakerMichelElectronAlgorithm::construct_ta ( ) const

private

Definition at line 83 of file TAMakerMichelElectronAlgorithm.cpp.

{
 
  TriggerPrimitive latest_tp_in_window = m_current_window.inputs.back();
 
  TriggerActivity ta;
  ta.time_start = m_current_window.time_start;
  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.
  ta.time_peak = latest_tp_in_window.samples_to_peak * 32 + latest_tp_in_window.time_start;  // FIXME: Replace STP to `time_peak` conversion.
  ta.time_activity = ta.time_peak;
  ta.channel_start = latest_tp_in_window.channel;
  ta.channel_end = latest_tp_in_window.channel;
  ta.channel_peak = latest_tp_in_window.channel;
  ta.adc_integral = m_current_window.adc_integral;
  ta.adc_peak = latest_tp_in_window.adc_peak;
  ta.detid = latest_tp_in_window.detid;
  ta.type = TriggerActivity::Type::kTPC;
  ta.algorithm = TriggerActivity::Algorithm::kMichelElectron;
  ta.inputs = m_current_window.inputs;
 
  return ta;
}

dump_tp()

void TAMakerMichelElectronAlgorithm::dump_tp ( TriggerPrimitive const & input_tp )

private

Definition at line 340 of file TAMakerMichelElectronAlgorithm.cpp.

{
  std::ofstream outfile;
  outfile.open("coldbox_tps.txt", std::ios_base::app);
 
  // Output relevant TP information to file
  outfile << input_tp.time_start << " ";          // Start time of TP
  outfile << input_tp.samples_over_threshold << " "; // in multiples of 25
  outfile << input_tp.samples_to_peak << " ";           //
  outfile << input_tp.channel << " ";             // Offline channel ID
  outfile << input_tp.adc_integral << " ";        // ADC Sum
  outfile << input_tp.adc_peak << " ";            // ADC Peak Value
  outfile << input_tp.detid << " ";               // Det ID - Identifies detector element, APA or PDS part etc...
  outfile.close();
 
  return;
}

dump_window_record()

void TAMakerMichelElectronAlgorithm::dump_window_record ( )

private

Definition at line 313 of file TAMakerMichelElectronAlgorithm.cpp.

{
  // FIX ME: Need to index this outfile in the name by detid or something similar.
  std::ofstream outfile;
  outfile.open("window_record_tam.csv", std::ios_base::app);
 
  for (auto window : m_window_record) {
    outfile << window.time_start << ",";
    outfile << window.inputs.back().time_start << ",";
    outfile << window.inputs.back().time_start - window.time_start << ","; // window_length - from TP start times
    outfile << window.adc_integral << ",";
    outfile << window.n_channels_hit() << ",";       // Number of unique channels with hits
    outfile << window.inputs.size() << ",";          // Number of TPs in Window
    outfile << window.inputs.back().channel << ",";  // Last TP Channel ID
    outfile << window.inputs.front().channel << ","; // First TP Channel ID
    outfile << longest_activity().size() << std::endl;             // New adjacency value for the window
  }
 
  outfile.close();
 
  m_window_record.clear();
 
  return;
}

longest_activity()

std::vector< TriggerPrimitive > TAMakerMichelElectronAlgorithm::longest_activity ( ) const

private

Definition at line 107 of file TAMakerMichelElectronAlgorithm.cpp.

{
  // This function attempts to return a vector of hits that correspond to the longest
  // piece of activity in the current window. The logic follows that from the HMA
  // check_adjacency() function and further details can be found there.
  std::vector<TriggerPrimitive> trackHits;
  std::vector<TriggerPrimitive> finalHits; // The vector of track hits, which we return
 
  uint16_t adj = 1;              // Initialise adjacency, 1 for the first wire.
  uint16_t max = 0;
  unsigned int channel = 0;      // Current channel ID
  unsigned int next_channel = 0; // Next channel ID
  unsigned int next = 0;         // The next position in the hit channels vector
  unsigned int tol_count = 0;    // Tolerance count, should not pass adj_tolerance
 
  // Generate a channelID ordered list of hit channels for this window
  std::vector<TriggerPrimitive> hitList;
  for (auto tp : m_current_window.inputs) {
    hitList.push_back(tp);
  }
  std::sort(hitList.begin(), hitList.end(), [](TriggerPrimitive a, TriggerPrimitive b)
            { return a.channel < b.channel; });
 
  // ADAJACENCY LOGIC ====================================================================
  // =====================================================================================
  // Adjcancency Tolerance = Number of times prepared to skip missed hits before resetting 
  // the adjacency count. This accounts for things like dead channels / missed TPs. The 
  // maximum gap is 4 which comes from tuning on December 2021 coldbox data, and June 2022 
  // coldbox runs.
  for (int i = 0; i < hitList.size(); ++i) {
 
    next = (i + 1) % hitList.size(); // Loops back when outside of channel list range
    channel = hitList.at(i).channel;
    next_channel = hitList.at(next).channel; // Next channel with a hit
 
    if (trackHits.size() == 0 ){ trackHits.push_back(hitList.at(i)); }
 
    // End of vector condition.
    if (next_channel == 0) { next_channel = channel - 1; }
 
    // Skip same channel hits for adjacency counting, but add to the track!
    if (next_channel == channel) { 
      trackHits.push_back(hitList.at(next));
      continue; }
 
    // If next hit is on next channel, increment the adjacency count.
    else if (next_channel == channel + 1){ 
      trackHits.push_back(hitList.at(next));
      ++adj; }
 
    // If next channel is not on the next hit, but the 'second next', increase adjacency 
    // but also tally up with the tolerance counter.
    else if (((next_channel == channel + 2) || (next_channel == channel + 3) ||
              (next_channel == channel + 4) || (next_channel == channel + 5))
             && (tol_count < m_adj_tolerance)) {
      trackHits.push_back(hitList.at(next));
      ++adj;
      for (int i = 0 ; i < next_channel-channel ; ++i){ ++tol_count; }
    }
 
    // If next hit isn't within reach, end the adjacency count and check for a new max.
    // Reset variables for next iteration.
    else {
      if (adj > max) { 
        max = adj;
        finalHits.clear(); // Clear previous track
        for (auto h : trackHits){
          finalHits.push_back(h);
        }  
      }
      adj = 1;
      tol_count = 0;
      trackHits.clear();
    }
  }
 
  return finalHits;
}

process()

void TAMakerMichelElectronAlgorithm::process ( const TriggerPrimitive & input_tp, std::vector< TriggerActivity > & output_ta )

virtual

TP processing function that creates & fills TAs.

Parameters

input_tp[in]	Input TP for the triggering algorithm
output_ta[out]	Output vector of TAs to fill by the algorithm

Implements triggeralgs::TriggerActivityMaker.

Definition at line 20 of file TAMakerMichelElectronAlgorithm.cpp.

{
 
  // The first time operator() is called, reset the window object.
  if (m_current_window.is_empty()) {
    m_current_window.reset(input_tp);
    m_primitive_count++;
    return;
  }
 
  // If the difference between the current TP's start time and the start of the window
  // is less than the specified window size, add the TP to the window.
  if ((input_tp.time_start - m_current_window.time_start) < m_window_length) {
    m_current_window.add(input_tp);
  }
 
  // Check Michel Candidate ========================================================
  // We've filled the window, now require a sufficient length track AND that the track
  // has a potential Bragg P, and then a kink.
  else if (longest_activity().size() > m_adjacency_threshold) {
     
     
     // We have a good length acitivity, now search for Bragg peak and kinks
     std::vector<TriggerPrimitive> trackHits = longest_activity();
     
     if (check_bragg_peak(trackHits)){
       if (check_kinks(trackHits)){
         TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TAM:ME] Emitting a trigger for candidate Michel event.";
         output_ta.push_back(construct_ta());
         m_current_window.reset(input_tp);
       } // Kinks 
     } // Bragg peak
 
  }
 
  // Otherwise, slide the window along using the current TP.
  else {
    m_current_window.move(input_tp, m_window_length);
  }
 
  m_primitive_count++;
  return;
}

Member Data Documentation

debug

bool triggeralgs::TAMakerMichelElectronAlgorithm::debug = false

private

Definition at line 116 of file TAMakerMichelElectronAlgorithm.hpp.

index

int triggeralgs::TAMakerMichelElectronAlgorithm::index = 0

private

Definition at line 120 of file TAMakerMichelElectronAlgorithm.hpp.

m_adj_tolerance

uint16_t triggeralgs::TAMakerMichelElectronAlgorithm::m_adj_tolerance = 3

private

Definition at line 119 of file TAMakerMichelElectronAlgorithm.hpp.

m_adjacency_threshold

uint16_t triggeralgs::TAMakerMichelElectronAlgorithm::m_adjacency_threshold = 15

private

Definition at line 117 of file TAMakerMichelElectronAlgorithm.hpp.

m_current_window

Window triggeralgs::TAMakerMichelElectronAlgorithm::m_current_window

private

Definition at line 111 of file TAMakerMichelElectronAlgorithm.hpp.

m_max_adjacency

int triggeralgs::TAMakerMichelElectronAlgorithm::m_max_adjacency = 0

private

Definition at line 118 of file TAMakerMichelElectronAlgorithm.hpp.

m_primitive_count

uint64_t triggeralgs::TAMakerMichelElectronAlgorithm::m_primitive_count = 0

private

Definition at line 112 of file TAMakerMichelElectronAlgorithm.hpp.

m_window_length

timestamp_t triggeralgs::TAMakerMichelElectronAlgorithm::m_window_length = 50000

private

Definition at line 123 of file TAMakerMichelElectronAlgorithm.hpp.

m_window_record

std::vector<Window> triggeralgs::TAMakerMichelElectronAlgorithm::m_window_record

private

Definition at line 129 of file TAMakerMichelElectronAlgorithm.hpp.

ta_adc

uint16_t triggeralgs::TAMakerMichelElectronAlgorithm::ta_adc = 0

private

Definition at line 121 of file TAMakerMichelElectronAlgorithm.hpp.

ta_channels

uint16_t triggeralgs::TAMakerMichelElectronAlgorithm::ta_channels = 0

private

Definition at line 122 of file TAMakerMichelElectronAlgorithm.hpp.

---

The documentation for this class was generated from the following files:

• /github/workspace/dunedaq/sourcecode/triggeralgs/include/triggeralgs/MichelElectron/TAMakerMichelElectronAlgorithm.hpp
• /github/workspace/dunedaq/sourcecode/triggeralgs/src/TAMakerMichelElectronAlgorithm.cpp