trigger::TrigMipReco Class Reference

Public Member Functions
	TrigMipReco (const std::string &name, framework::Process &process)
void	configure (framework::config::Parameters &ps) override
	Callback for the EventProcessor to configure itself from the given set of parameters.
void	produce (framework::Event &event) override
	Process the event and put new data products into it.
void	onNewRun (const ldmx::RunHeader &rh) override
	onNewRun is the first function called for each processor after the conditions are fully configured and accessible.
void	onProcessEnd () override
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
Public Member Functions inherited from framework::Producer
	Producer (const std::string &name, Process &process)
	Class constructor.
virtual void	process (Event &event) final
	Processing an event for a Producer is calling produce.
Public Member Functions inherited from framework::EventProcessor
	DECLARE_FACTORY (EventProcessor, EventProcessor *, const std::string &, Process &)
	declare that we have a factory for this class
	EventProcessor (const std::string &name, Process &process)
	Class constructor.
virtual	~EventProcessor ()=default
	Class destructor.
virtual void	beforeNewRun (ldmx::RunHeader &run_header)
	Callback for Producers to add parameters to the run header before conditions are initialized.
virtual void	onFileOpen (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.
virtual void	onFileClose (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.
virtual void	onProcessStart ()
	Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.
template<class T >
const T &	getCondition (const std::string &condition_name)
	Access a conditions object for the current event.
TDirectory *	getHistoDirectory ()
	Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.
void	setStorageHint (framework::StorageControl::Hint hint)
	Mark the current event as having the given storage control hint from this module_.
void	setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
	Mark the current event as having the given storage control hint from this module and the given purpose string.
int	getLogFrequency () const
	Get the current logging frequency from the process.
int	getRunNumber () const
	Get the run number from the process.
std::string	getName () const
	Get the processor name.
void	createHistograms (const std::vector< framework::config::Parameters > &histos)
	Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.

Private Attributes
int	nevents_ {0}
double	processing_time_ {0.}
std::map< std::string, double >	profiling_map_
std::string	hit_coll_name_
std::string	hit_coll_passname_
std::string	pass_coll_name_
bool	calorimeter_type_is_hcal_
int	max_layer_
int	min_track_length_
float	hcal_min_energy_
float	ecal_min_energy_
float	ecal_max_energy_
float	search_radius_
float	isolation_e_cut_
float	hole_fraction_max_

Additional Inherited Members
Protected Member Functions inherited from framework::EventProcessor
void	abortEvent ()
	Abort the event immediately.
Protected Attributes inherited from framework::EventProcessor
HistogramPool	histograms_
	helper object for making and filling histograms
NtupleManager &	ntuple_ {NtupleManager::getInstance()}
	Manager for any ntuples.
logging::logger	the_log_
	The logger for this EventProcessor.

Detailed Description

Definition at line 22 of file TrigMipReco.h.

Constructor & Destructor Documentation

TrigMipReco()

trigger::TrigMipReco::TrigMipReco ( const std::string & name, framework::Process & process )

inline

Definition at line 24 of file TrigMipReco.h.

      : framework::Producer(name, process) {}

Member Function Documentation

configure()

void trigger::TrigMipReco::configure ( framework::config::Parameters & parameters )

overridevirtual

Callback for the EventProcessor to configure itself from the given set of parameters.

The parameters a processor has access to are the member variables of the python class in the sequence that has class_name equal to the EventProcessor class name.

For an example, look at MyProcessor.

Parameters

parameters

Parameters for configuration.

Reimplemented from framework::EventProcessor.

Definition at line 14 of file TrigMipReco.cxx.

                                                         {
  hit_coll_name_ = ps.get<std::string>("hit_coll_name");
  pass_coll_name_ = ps.get<std::string>("pass_coll_name");
  hit_coll_passname_ = ps.get<std::string>("hit_coll_passname");
  calorimeter_type_is_hcal_ = ps.get<bool>("calorimeter_type_is_hcal");
 
  max_layer_ = ps.get<int>("max_layer", 32);
  // mm
  search_radius_ = ps.get<double>("search_radius", 50.0);
  min_track_length_ = ps.get<int>("min_track_length", 5);
  // MeV; Change as needed
  isolation_e_cut_ = ps.get<double>("isolation_e_cut", 180.0);
  hole_fraction_max_ = ps.get<double>("hole_fraction_max", 0.2);
 
  if (calorimeter_type_is_hcal_) {
    // MIP peak is 10-11 MeV
    hcal_min_energy_ = ps.get<double>("hcal_min_energy", 8.0);
  } else {  // ECAL
    // MIP peak around 17 MeV
    ecal_min_energy_ = ps.get<double>("ecal_min_energy", 3.0);
    ecal_max_energy_ = ps.get<double>("ecal_max_energy", 26.0);
  }
}

References framework::config::Parameters::get().

onNewRun()

void trigger::TrigMipReco::onNewRun ( const ldmx::RunHeader & rh )

overridevirtual

onNewRun is the first function called for each processor after the conditions are fully configured and accessible.

This is where you could create single-processors, multi-event calculation objects.

Reimplemented from framework::EventProcessor.

Definition at line 10 of file TrigMipReco.cxx.

                                                  {
  profiling_map_["processing_time_"] = 0;
}

onProcessEnd()

void trigger::TrigMipReco::onProcessEnd ( )

overridevirtual

Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.

Reimplemented from framework::EventProcessor.

Definition at line 308 of file TrigMipReco.cxx.

                               {
  ldmx_log(info) << "AVG Time/Event: " << std::fixed << std::setprecision(3)
                 << processing_time_ / nevents_ << " ms";
}

produce()

void trigger::TrigMipReco::produce ( framework::Event & event )

overridevirtual

Process the event and put new data products into it.

Parameters

event

The Event to process.

Implements framework::Producer.

Definition at line 38 of file TrigMipReco.cxx.

                                               {
  auto start = std::chrono::high_resolution_clock::now();
  nevents_++;
 
  if (!event.exists(hit_coll_name_, hit_coll_passname_)) return;
 
  const auto calo_hits = event.getObject<std::vector<TrigCaloHit>>(
      hit_coll_name_, hit_coll_passname_);
 
  if (calorimeter_type_is_hcal_) {  // HCAL MIP Reconstruction
    std::vector<TrigCaloHit> sorted_hits;
    int even_matrix[24][5] = {};
    int even_start[5] = {99, 99, 99, 99, 99};
    int even_end[5] = {};
    int even_counts[5] = {};
    int odd_matrix[24][5] = {};
    int odd_start[5] = {99, 99, 99, 99, 99};
    int odd_end[5] = {};
    int odd_counts[5] = {};
    // Start in first 5 layers
    constexpr int layer_start = 80;
 
    for (const auto& tp : calo_hits) {
      if (tp.section() > 0 || tp.energy() < hcal_min_energy_ ||
          tp.layer() > 47) {
        continue;
      }
 
      sorted_hits.push_back(tp);
      const int layer_index = tp.layer() / 2;
      const int strip = tp.strip();
      if (tp.layer() % 2) {
        odd_matrix[layer_index][strip] = 1;
        odd_start[strip] = std::min(odd_start[strip], tp.layer());
        odd_end[strip] = std::max(odd_end[strip], tp.layer());
      } else {
        even_matrix[layer_index][strip] = 1;
        even_start[strip] = std::min(even_start[strip], tp.layer());
        even_end[strip] = std::max(even_end[strip], tp.layer());
      }
    }
 
    for (int i = 0; i < 24; i++) {
      for (int j = 0; j < 5; j++) {
        if (even_matrix[i][j]) {
          even_counts[j]++;
        }
        if (odd_matrix[i][j]) {
          odd_counts[j]++;
        }
      }
    }
 
    // straight MIP reco
    std::vector<TrigMip> mips;
    // 5 elements in the even/odd_start matrices
    for (int i = 0; i < 5; i++) {
      if (odd_start[i] < layer_start) {
        TrigMip m;
        m.setStartLayer(odd_start[i]);
        m.setEndLayer(odd_end[i]);
        m.setNHits(odd_counts[i]);
        m.setLength(odd_end[i] - odd_start[i] + 1);
        int holes = m.length() / 2 - m.nHits();
        if (holes < 0) {
          holes = 0;
        }
        m.setNHoles(holes);
        mips.push_back(m);
      }
      if (even_start[i] < layer_start) {
        TrigMip m;
        m.setStartLayer(even_start[i]);
        m.setEndLayer(even_end[i]);
        m.setNHits(even_counts[i]);
        m.setLength(even_end[i] - even_start[i] + 1);
        int holes = m.length() / 2 - m.nHits();
        if (holes < 0) {
          holes = 0;
        }
        m.setNHoles(holes);
        mips.push_back(m);
      }
    }
 
    std::sort(mips.begin(), mips.end());
    event.add(pass_coll_name_, mips);
 
    auto end = std::chrono::high_resolution_clock::now();
    auto time_diff = end - start;
    processing_time_ +=
        std::chrono::duration<double, std::milli>(time_diff).count();
 
    return;
    // ECAL MIP Reconstruction
  } else {
    const float radius_cut_2 = search_radius_ * search_radius_;
    std::map<int, std::vector<TrigCaloHit>> layer_hits;
    std::set<const TrigCaloHit*> used_hits;
    std::vector<std::vector<const TrigCaloHit*>> candidate_tracks;
    std::map<const TrigCaloHit*, size_t> hit_to_best_track;
 
    // Filter for section = 0 and hits < 33 layers
    for (const auto& hit : calo_hits) {
      if (hit.section() > 0 || hit.layer() > max_layer_) continue;
      layer_hits[hit.layer()].push_back(hit);
    }
 
    // Find mip seeds
    for (const auto& [seed_layer, seeds] : layer_hits) {
      for (const auto& seed : seeds) {
        // Skip if hit already used or outside mip energy range
        if (used_hits.count(&seed) || seed.energy() < ecal_min_energy_ ||
            seed.energy() > ecal_max_energy_) {
          continue;
        }
 
        std::vector<const TrigCaloHit*> track{&seed};
        // Most recent hit in track
        const TrigCaloHit* last = &seed;
        int holes = 0;
        float growth_factor = 1.0f;
 
        // Look layer by layer for next hit within dR
        for (int l = seed.layer() + 1; l <= max_layer_; ++l) {
          const TrigCaloHit* best_hit = nullptr;
          // Grow search window if there is a hole
          float best_d_r_2 = radius_cut_2 * growth_factor * growth_factor;
 
          for (const auto& cand : layer_hits[l]) {
            if (used_hits.count(&cand) || cand.energy() < ecal_min_energy_ ||
                cand.energy() > ecal_max_energy_) {
              continue;
            }
 
            // Corrects for layer shift in x-direction, calculated as 4.82 mm
            const float layer_shift_last =
                (last->layer() % 2 == 0) ? 0.0f : 4.82f;
            const float layer_shift_cand =
                (cand.layer() % 2 == 0) ? 0.0f : 4.82f;
            const float dx = (cand.positionX() - layer_shift_cand) -
                             (last->positionX() - layer_shift_last);
            const float dy = cand.positionY() - last->positionY();
            const float d_r_2 = dx * dx + dy * dy;
 
            if (d_r_2 < best_d_r_2) {
              best_d_r_2 = d_r_2;
              // Closest unused hit in next layer
              best_hit = &cand;
            }
          }
 
          if (best_hit) {
            // Builds track from best hits
            track.push_back(best_hit);
            last = best_hit;
            holes = 0;
            // Reset search window
            growth_factor = 1.0f;
          } else {
            holes++;
            // Keep expanding
            growth_factor = static_cast<float>(holes + 1);
          }
        }
 
        bool is_isolated = true;
        if (track.size() >= min_track_length_) {
          // Isolation area energy check
          for (const auto* hit : track) {
            const int layer = hit->layer();
            const float hit_x = hit->positionX();
            const float hit_y = hit->positionY();
            float sum_e = 0.0f;
 
            for (const auto& cand : layer_hits[layer]) {
              // Skips self
              if (&cand == hit) continue;
 
              const float dx = cand.positionX() - hit_x;
              const float dy = cand.positionY() - hit_y;
              const float d_r_2 = dx * dx + dy * dy;
 
              if (d_r_2 < radius_cut_2) {
                sum_e += cand.energy();
              }
            }
 
            if (sum_e >= isolation_e_cut_) {
              is_isolated = false;
              // Adds used hits to vector so they cannot be used again
              used_hits.insert(hit);
              break;
            }
          }
 
          if (!is_isolated) {
            // Reject track if any hit is not isolated
            continue;
          }
 
          const size_t i = candidate_tracks.size();
          candidate_tracks.push_back(track);
 
          for (const auto* hit : track) {
            if (!hit_to_best_track.count(hit) ||
                candidate_tracks[i].size() >
                    candidate_tracks[hit_to_best_track[hit]].size()) {
              hit_to_best_track[hit] = i;
              // Adds used hits to vector so they cannot be used again
              used_hits.insert(hit);
            }
          }
        }
      }
    }
 
    std::set<size_t> valid_track_i_ds;
    for (const auto& [hit, idx] : hit_to_best_track) {
      valid_track_i_ds.insert(idx);
    }
 
    std::vector<TrigMip> mips;
    for (const size_t idx : valid_track_i_ds) {
      const auto& track = candidate_tracks[idx];
      TrigMip mip;
      mip.setStartLayer(track.front()->layer());
      mip.setEndLayer(track.back()->layer());
      mip.setNHits(track.size());
      mip.setLength(track.back()->layer() - track.front()->layer());
      int holes = mip.length() - mip.nHits();
      if (holes < 0) {
        holes = 0;
      }
      mip.setNHoles(holes);
 
      const float hole_fraction =
          static_cast<float>(mip.nHoles()) / mip.length();
      // Remove mip tracks with hole fraction > 0.2
      if (hole_fraction >= hole_fraction_max_) continue;
 
      float total_isolation_e_sum = 0.0f;
      for (const auto* hit : track) {
        const int layer = hit->layer();
        const float hit_x = hit->positionX();
        const float hit_y = hit->positionY();
        for (const auto& cand : layer_hits[layer]) {
          if (&cand == hit) continue;
          const float dx = cand.positionX() - hit_x;
          const float dy = cand.positionY() - hit_y;
          const float d_r_2 = dx * dx + dy * dy;
          if (d_r_2 < radius_cut_2) {
            total_isolation_e_sum += cand.energy();
          }
        }
      }
      mip.setSumEinIsolationRegion(total_isolation_e_sum);
      mips.push_back(mip);
    }
 
    std::sort(mips.begin(), mips.end());
    event.add(pass_coll_name_, mips);
  }
 
  auto end = std::chrono::high_resolution_clock::now();
  auto time_diff = end - start;
  processing_time_ +=
      std::chrono::duration<double, std::milli>(time_diff).count();
}

References framework::Event::exists().

Member Data Documentation

calorimeter_type_is_hcal_

bool trigger::TrigMipReco::calorimeter_type_is_hcal_

private

Definition at line 54 of file TrigMipReco.h.

ecal_max_energy_

float trigger::TrigMipReco::ecal_max_energy_

private

Definition at line 63 of file TrigMipReco.h.

ecal_min_energy_

float trigger::TrigMipReco::ecal_min_energy_

private

Definition at line 61 of file TrigMipReco.h.

hcal_min_energy_

float trigger::TrigMipReco::hcal_min_energy_

private

Definition at line 59 of file TrigMipReco.h.

hit_coll_name_

std::string trigger::TrigMipReco::hit_coll_name_

private

Definition at line 49 of file TrigMipReco.h.

hit_coll_passname_

std::string trigger::TrigMipReco::hit_coll_passname_

private

Definition at line 50 of file TrigMipReco.h.

hole_fraction_max_

float trigger::TrigMipReco::hole_fraction_max_

private

Definition at line 70 of file TrigMipReco.h.

isolation_e_cut_

float trigger::TrigMipReco::isolation_e_cut_

private

Definition at line 67 of file TrigMipReco.h.

max_layer_

int trigger::TrigMipReco::max_layer_

private

Definition at line 56 of file TrigMipReco.h.

min_track_length_

int trigger::TrigMipReco::min_track_length_

private

Definition at line 57 of file TrigMipReco.h.

nevents_

int trigger::TrigMipReco::nevents_ {0}

private

Definition at line 45 of file TrigMipReco.h.

{0};

pass_coll_name_

std::string trigger::TrigMipReco::pass_coll_name_

private

Definition at line 52 of file TrigMipReco.h.

processing_time_

double trigger::TrigMipReco::processing_time_ {0.}

private

Definition at line 46 of file TrigMipReco.h.

{0.};

profiling_map_

std::map<std::string, double> trigger::TrigMipReco::profiling_map_

private

Definition at line 47 of file TrigMipReco.h.

search_radius_

float trigger::TrigMipReco::search_radius_

private

Definition at line 65 of file TrigMipReco.h.

---

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

• Trigger/Algo/include/Trigger/TrigMipReco.h
• Trigger/Algo/src/Trigger/TrigMipReco.cxx