dqm::EcalClusterAnalyzer Class Reference

Public Member Functions
	EcalClusterAnalyzer (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	analyze (const framework::Event &event) override
	Process the event and make histograms or summaries.
Public Member Functions inherited from framework::Analyzer
	Analyzer (const std::string &name, Process &process)
	Class constructor.
virtual void	process (Event &event) final
	Processing an event for an Analyzer is calling analyze.
virtual void	beforeNewRun (ldmx::RunHeader &run_header) final
	Don't allow Analyzers to add parameters to the run header.
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	onNewRun (const ldmx::RunHeader &run_header)
	Callback for the EventProcessor to take any necessary action when the run being processed changes.
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.
virtual void	onProcessEnd ()
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
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
bool	use_simulated_electron_number_
	Use the number of simulated electrons instead of the number of determined by the TS track counting.
int	nbr_of_electrons_
	What is the number of electrons in the event?
std::string	ecal_sim_hit_coll_
std::string	ecal_sim_hit_pass_
std::string	rec_hit_coll_name_
std::string	rec_hit_pass_name_
std::string	cluster_coll_name_
std::string	cluster_pass_name_
std::string	ecal_sp_hits_coll_name_
std::string	ecal_sp_hits_pass_name_
double	mixed_hit_cutoff_
bool	inverse_skim_
int	n_ecal_clusters_min_

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 21 of file EcalClusterAnalyzer.h.

Constructor & Destructor Documentation

EcalClusterAnalyzer()

dqm::EcalClusterAnalyzer::EcalClusterAnalyzer ( const std::string & name, framework::Process & process )

inline

Definition at line 23 of file EcalClusterAnalyzer.h.

      : Analyzer(name, process) {}

Member Function Documentation

analyze()

void dqm::EcalClusterAnalyzer::analyze ( const framework::Event & event )

overridevirtual

Process the event and make histograms or summaries.

Parameters

event

The Event to analyze

Implements framework::Analyzer.

Definition at line 39 of file EcalClusterAnalyzer.cxx.

                                                             {
  const auto& ecal_rec_hits{event.getCollection<ldmx::EcalHit>(
      rec_hit_coll_name_, rec_hit_pass_name_)};
  const auto& ecal_sim_hits{event.getCollection<ldmx::SimCalorimeterHit>(
      ecal_sim_hit_coll_, ecal_sim_hit_pass_)};
  const auto& ecal_clusters{event.getCollection<ldmx::EcalCluster>(
      cluster_coll_name_, cluster_pass_name_)};
 
  // Determine the number of recoil electrons in the event
  // By default from the TS track counting
  int nbr_of_electrons{event.getElectronCount()};
  // If configured to use the simulated electron number, use that instead
  if (use_simulated_electron_number_) {
    nbr_of_electrons = nbr_of_electrons_;
  }
 
  std::map<int, int> layer_cluster_count;
  for (const auto& cluster : ecal_clusters) {
    auto layer = cluster.getLayer();
    layer_cluster_count[layer]++;
  }
 
  int total_clusters = 0;
  for (const auto& [layer, count] : layer_cluster_count) {
    total_clusters += count;
  }
 
  int n_ecal_clusters = 0;
  if (layer_cluster_count.size() != 0) {
    n_ecal_clusters = static_cast<int>(std::round(
        static_cast<double>(total_clusters) / layer_cluster_count.size()));
  }
 
  ldmx_log(info) << "Avg number of clusters per layer: " << n_ecal_clusters;
  // Fill histograms with the number of clusters
  histograms_.fill("number_of_clusters", total_clusters);
  histograms_.fill("number_of_clusters_per_layer", n_ecal_clusters);
  histograms_.fill("number_of_clusters_first_layer", layer_cluster_count[0]);
 
  // Fill simplied 3-bin histogram to check the prediction
  if (n_ecal_clusters == nbr_of_electrons) {
    // correct
    histograms_.fill("correctly_predicted_events", 1);
  } else if (n_ecal_clusters < nbr_of_electrons) {
    // undercounting
    histograms_.fill("correctly_predicted_events", 0);
  } else if (n_ecal_clusters > nbr_of_electrons) {
    // overcounting
    histograms_.fill("correctly_predicted_events", 2);
  }
 
  std::unordered_map<int, std::pair<int, std::vector<double>>> hit_info;
  hit_info.reserve(ecal_rec_hits.size());
 
  // Determine the truth information for the recoil electron
  std::vector<std::vector<float>> sp_electron_positions;
  const auto& ecal_sp_hits{event.getCollection<ldmx::SimTrackerHit>(
      ecal_sp_hits_coll_name_, ecal_sp_hits_pass_name_)};
 
  std::vector<ldmx::SimTrackerHit> sorted_sp_hits = ecal_sp_hits;
  std::sort(sorted_sp_hits.begin(), sorted_sp_hits.end(),
            [](const ldmx::SimTrackerHit& a, const ldmx::SimTrackerHit& b) {
              return a.getTrackID() < b.getTrackID();
            });
 
  ldmx_log(trace) << "Number of ECal Scoring Plane Hits: "
                  << sorted_sp_hits.size();
 
  // Collect positions of all recoil electrons on the SP
  // relying on the track ID to identify them
  unsigned int n_filled = 0;
  for (const ldmx::SimTrackerHit& sp_hit : sorted_sp_hits) {
    if (sp_hit.getPdgID() != 11) continue;
    if (sp_hit.getMomentum()[2] <= 0) continue;
    ldmx::SimSpecialID hit_id(sp_hit.getID());
    // Ecal scoring plane is plane 31
    if (hit_id.plane() != 31) continue;
    if (n_filled < nbr_of_electrons) {
      ldmx_log(trace) << "\tSP Hit to be added with Track ID : "
                      << sp_hit.getTrackID() << ", SP Hit Position ("
                      << sp_hit.getPosition()[0] << ", "
                      << sp_hit.getPosition()[1] << ", "
                      << sp_hit.getPosition()[2] << ") mm";
      sp_electron_positions.push_back(sp_hit.getPosition());
      n_filled++;
    }
  }
 
  ldmx_log(info) << "Number of ECal CLUE clusters: " << n_ecal_clusters
                 << ", TS counted electrons: " << nbr_of_electrons
                 << ", SP electrons: " << sp_electron_positions.size();
 
  std::map<int, float> true_energy;
  std::map<int, float> delta_energy;
  double sp_ele_dist{9999.};
  if (nbr_of_electrons == 2 && sp_electron_positions.size() > 1) {
    // Measures sp_ele_distance between two electrons in the ECal scoring plane
    // TODO: generalize for n electrons
    std::vector<float> pos1;
    std::vector<float> pos2;
    pos1 = sp_electron_positions[0];
    pos2 = sp_electron_positions[1];
    sp_ele_dist = std::sqrt((pos1[0] - pos2[0]) * (pos1[0] - pos2[0]) +
                            (pos1[1] - pos2[1]) * (pos1[1] - pos2[1]));
 
    histograms_.fill("sp_distance", sp_ele_dist);
 
  }  // end block about the scoring plane hits
 
  ldmx_log(trace) << "Distance between the two e- in the ECal scoring plane: "
                  << sp_ele_dist << " mm";
 
  double tot_event_energy = 0;
  std::vector<double> tot_origin_edep;
  tot_origin_edep.resize(nbr_of_electrons_ + 1);
  int n_mixed = 0;
 
  // Loop over the rechits and find the matching simhits
  ldmx_log(trace) << "Loop over the rechits and find the matching simhits";
  for (const auto& hit : ecal_rec_hits) {
    auto it = std::find_if(
        ecal_sim_hits.begin(), ecal_sim_hits.end(),
        [&hit](const auto& sim_hit) { return sim_hit.getID() == hit.getID(); });
    if (it != ecal_sim_hits.end()) {
      // if found a simhit matching this rechit
      ldmx_log(trace) << "\tFound simhit matching rechit with ID"
                      << hit.getID();
      int ancestor = 0;
      int prev_ancestor = 0;
      bool tagged = false;
      int tag = 0;
      std::vector<double> edep;
      edep.resize(nbr_of_electrons_ + 1);
      double e_tot = 0;  // keep track of total from all counted ancestors
      ldmx_log(trace) << "\t\tIt has " << it->getNumberOfContribs()
                      << " contribs. ";
      for (int i = 0; i < it->getNumberOfContribs(); i++) {
        // for each contrib in this simhit
        const auto& contrib = it->getContrib(i);
        // get origin electron ID
        ancestor = contrib.origin_id_;
        ldmx_log(trace) << "\t\t\tAncestor ID " << ancestor << " with edep "
                        << contrib.edep_;
        tot_event_energy += contrib.edep_;
        // store energy from this contrib at index = origin electron ID
        if (ancestor <= nbr_of_electrons) {
          edep[ancestor] += contrib.edep_;
          tot_origin_edep[ancestor] += contrib.edep_;
          e_tot += contrib.edep_;
        }
        if (!tagged && i != 0 && prev_ancestor != ancestor) {
          // if origin electron ID does not match previous origin electron ID
          // this hit has contributions from several electrons, ie mixed case
          tag = 0;
          tagged = true;
          ldmx_log(trace) << "\t\t\t\tMixed hit! Ancestor ID changed to "
                          << ancestor;
        }
        prev_ancestor = ancestor;
      }  // over contribs
      // now check if mixed really means mixed, i.e. more than small fraction
      // from a second electron.
      if (tagged) {
        for (int i = 1; i < nbr_of_electrons_ + 1; i++) {
          if (edep[i] / e_tot >
              1 - mixed_hit_cutoff_) {  // one ancestor contributes at least the
                                        // complement to the allowed mixing
                                        // fraction
            tagged = false;
            ancestor =
                distance(edep.begin(), max_element(edep.begin(), edep.end()));
            ldmx_log(trace)
                << "\t\t\t\tUndid mixed hit tagging, now ancestor = "
                << ancestor;
            break;
          }
        }
      }
      if (!tagged) {
        // if not tagged, hit was from a single electron (within acceptable
        // purity)
        tag = ancestor;  // prev_ancestor;
      } else
        n_mixed++;
      histograms_.fill("ancestors", tag);
      hit_info.insert({hit.getID(), std::make_pair(tag, edep)});
    }  // end if simhit found
  }  // end loop on the rechits
 
  // Loop over the clusters
  int clustered_hits = 0;
  ldmx_log(trace) << "Loop over the clusters, N = " << n_ecal_clusters;
  histograms_.fill("tag0frac_vs_SPdist", sp_ele_dist,
                   (float)n_mixed / ecal_rec_hits.size());
  ldmx_log(debug) << "Got " << n_mixed << " mixed hits, a fraction of "
                  << (float)n_mixed / ecal_rec_hits.size();
 
  if (ecal_clusters.size() >= 2) {
    float d_x =
        ecal_clusters[0].getCentroidX() - ecal_clusters[1].getCentroidX();
    float d_y =
        ecal_clusters[0].getCentroidY() - ecal_clusters[1].getCentroidY();
    float d_r = std::sqrt(d_x * d_x + d_y * d_y);
    histograms_.fill("cluster_distance", d_r);
    ldmx_log(trace) << "Gt cluster distance (0,1) = " << d_r;
  }
 
  for (const auto& cl : ecal_clusters) {
    auto layer = cl.getLayer();
    ldmx_log(trace) << "Cluster in layer " << layer
                    << ", energy: " << cl.getEnergy()
                    << ", number of hits: " << cl.getHitIDs().size();
    auto cluster_centroid_x = cl.getCentroidX();
    auto cluster_centroid_y = cl.getCentroidY();
    auto cluster_rms_x = cl.getRMSX();
    auto cluster_rms_y = cl.getRMSY();
 
    // Find the closest sp_electron_positions to the cluster centroid
    double min_distance = 9999.;
    double sp_clue_x_residuals = 9999.;
    double sp_clue_y_residuals = 9999.;
    for (const auto& sp_pos : sp_electron_positions) {
      double distance = std::sqrt(
          (sp_pos[0] - cluster_centroid_x) * (sp_pos[0] - cluster_centroid_x) +
          (sp_pos[1] - cluster_centroid_y) * (sp_pos[1] - cluster_centroid_y));
      if (distance < min_distance) {
        min_distance = distance;
        sp_clue_x_residuals = sp_pos[0] - cluster_centroid_x;
        sp_clue_y_residuals = sp_pos[1] - cluster_centroid_y;
      }
    }  // end loop on the scoring plane electron positions
    // Fill histogram with the distance to the closest scoring plane electron
    ldmx_log(trace) << "\tCluster centroid: (" << cluster_centroid_x << " +/- "
                    << cluster_rms_x << ", " << cluster_centroid_y << " +/- "
                    << cluster_rms_y
                    << " mm; min distance to SP electron: " << min_distance
                    << " mm";
    if (layer == 0) {
      histograms_.fill("sp_clue_distance", min_distance);
      histograms_.fill("sp_clue_x_residual", sp_clue_x_residuals);
      histograms_.fill("sp_clue_y_residual", sp_clue_y_residuals);
    }
    histograms_.fill("sp_clue_distance_vs_layer", layer, min_distance);
 
    // for each cluster
    // total number of hits coming from electron, index = electron ID
    std::vector<double> n_hits_from_electron;
    n_hits_from_electron.resize(nbr_of_electrons + 2);
    // total number of energy coming from electron, index = electron ID
    std::vector<double> energy_from_electron;
    energy_from_electron.resize(nbr_of_electrons + 2);
    double energy_sum = 0.;
    double n_sum = 0.;
    ldmx_log(trace) << "Looping over hits in the cluster";
    const auto& hit_ids = cl.getHitIDs();
    for (const auto& id : hit_ids) {
      // for each hit in cluster, find previously stored info
      auto it = hit_info.find(id);
      if (it != hit_info.end()) {
        auto t = it->second;
        // origin electron ID (or 0 for mixed)
        auto id_electron = t.first;
        // energy vector
        auto energies = t.second;
        // increment number of hits coming from this electron
        n_hits_from_electron[id_electron]++;
        n_sum++;
 
        double hit_energy_sum = 0.;
        for (int i = 1; i < nbr_of_electrons + 1; i++) {
          // loop through energy vector
          if (energies[i] > 0.) {
            energy_sum += energies[i];
            // add energy from electron i in this hit to total energy from
            // electron i in cluster
            energy_from_electron[i] += energies[i];
          }
        }
        // if mixed hit, add the total energy of this hit to mixed hit energy
        // counter
        if (id_electron == 0) energy_from_electron[0] += hit_energy_sum;
        energy_sum += hit_energy_sum;
 
        clustered_hits++;
      }  // end if hit info found
    }  // end loop on the hit IDs in the cluster
 
    if (energy_sum > 0) {
      // get largest energy contribution
      double max_energy_contribution = *max_element(
          energy_from_electron.begin(), energy_from_electron.end());
      std::string to_log;
      for (auto nb : energy_from_electron)
        to_log.append(std::to_string(nb) + " ");
      ldmx_log(debug) << "Energies vector is " << to_log;
 
      // energy purity = largest contribution / all energy
      histograms_.fill("energy_percentage",
                       100. * (max_energy_contribution / energy_sum));
      if (energy_from_electron[0] > 0.)
        histograms_.fill("mixed_hit_energy",
                         100. * (energy_from_electron[0] / energy_sum));
 
      histograms_.fill("total_energy_vs_hits", energy_sum,
                       cl.getHitIDs().size());
      histograms_.fill("total_energy_vs_purity", energy_sum,
                       100. * (max_energy_contribution / energy_sum));
 
      if (nbr_of_electrons == 2) {
        histograms_.fill("sp_ele_distance_vs_purity", sp_ele_dist,
                         100. * (max_energy_contribution / energy_sum));
      }
    }
    if (n_sum > 0) {
      double n_max = *max_element(n_hits_from_electron.begin(),
                                  n_hits_from_electron.end());
      histograms_.fill("same_ancestor", 100. * (n_max / n_sum));
    }
    // find the main contributor
    auto elt = distance(
        energy_from_electron.begin(),
        max_element(energy_from_electron.begin(), energy_from_electron.end()));
    ldmx_log(debug) << "Found that the maximum contributing trackID is " << elt;
    delta_energy[elt] = cl.getEnergy() - true_energy[elt];
    //    delta_energy[2]=e[2]-true_energy[2];
    histograms_.fill("cluster_RMSX", cl.getRMSX());
  }  // end loop on the clusters
  std::string more_log;
  for (auto nb : tot_origin_edep) more_log.append(std::to_string(nb) + " ");
  ldmx_log(debug) << "Edep per ancestor in event is " << more_log;
  ldmx_log(debug) << "Total energy deposited in event: " << tot_event_energy;
 
  histograms_.fill("dE_cl2_vs_cl1", delta_energy[1], delta_energy[2]);
  histograms_.fill("unclustered_hits", (ecal_rec_hits.size() - clustered_hits));
  histograms_.fill("total_rechits_in_event", ecal_rec_hits.size());
  histograms_.fill(
      "unclustered_hits_percentage",
      100. * (ecal_rec_hits.size() - clustered_hits) / ecal_rec_hits.size());
 
  if (inverse_skim_) {
    // inverse operation: drop events with enough clusters
    if (n_ecal_clusters > n_ecal_clusters_min_) {
      setStorageHint(framework::HINT_SHOULD_DROP);
    } else {
      setStorageHint(framework::HINT_SHOULD_KEEP);
    }
  } else {
    // normal operation: keep events with enough clusters
    if (n_ecal_clusters > n_ecal_clusters_min_) {
      setStorageHint(framework::HINT_SHOULD_KEEP);
    } else {
      setStorageHint(framework::HINT_SHOULD_DROP);
    }
  }
}

References framework::HistogramPool::fill(), framework::HINT_SHOULD_DROP, framework::HINT_SHOULD_KEEP, framework::EventProcessor::histograms_, nbr_of_electrons_, ldmx::SimSpecialID::plane(), framework::EventProcessor::setStorageHint(), and use_simulated_electron_number_.

configure()

void dqm::EcalClusterAnalyzer::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 EcalClusterAnalyzer.cxx.

                                                                 {
  use_simulated_electron_number_ =
      ps.get<bool>("use_simulated_electron_number");
  nbr_of_electrons_ = ps.get<int>("nbr_of_electrons");
 
  ecal_sim_hit_coll_ = ps.get<std::string>("ecal_sim_hit_coll");
  ecal_sim_hit_pass_ = ps.get<std::string>("ecal_sim_hit_pass");
 
  rec_hit_coll_name_ = ps.get<std::string>("rec_hit_coll_name");
  rec_hit_pass_name_ = ps.get<std::string>("rec_hit_pass_name");
 
  cluster_coll_name_ = ps.get<std::string>("cluster_coll_name");
  cluster_pass_name_ = ps.get<std::string>("cluster_pass_name");
 
  ecal_sp_hits_coll_name_ =
      ps.getParameter<std::string>("ecal_sp_hits_coll_name");
  ecal_sp_hits_pass_name_ =
      ps.getParameter<std::string>("ecal_sp_hits_pass_name");
  mixed_hit_cutoff_ = ps.getParameter<double>("mixed_hit_cutoff");
 
  inverse_skim_ = ps.get<bool>("inverse_skim");
  n_ecal_clusters_min_ = ps.get<int>("n_ecal_clusters_min");
  return;
}

References framework::config::Parameters::get(), nbr_of_electrons_, and use_simulated_electron_number_.

Member Data Documentation

cluster_coll_name_

std::string dqm::EcalClusterAnalyzer::cluster_coll_name_

private

Definition at line 50 of file EcalClusterAnalyzer.h.

cluster_pass_name_

std::string dqm::EcalClusterAnalyzer::cluster_pass_name_

private

Definition at line 53 of file EcalClusterAnalyzer.h.

ecal_sim_hit_coll_

std::string dqm::EcalClusterAnalyzer::ecal_sim_hit_coll_

private

Definition at line 38 of file EcalClusterAnalyzer.h.

ecal_sim_hit_pass_

std::string dqm::EcalClusterAnalyzer::ecal_sim_hit_pass_

private

Definition at line 41 of file EcalClusterAnalyzer.h.

ecal_sp_hits_coll_name_

std::string dqm::EcalClusterAnalyzer::ecal_sp_hits_coll_name_

private

Definition at line 56 of file EcalClusterAnalyzer.h.

ecal_sp_hits_pass_name_

std::string dqm::EcalClusterAnalyzer::ecal_sp_hits_pass_name_

private

Definition at line 58 of file EcalClusterAnalyzer.h.

inverse_skim_

bool dqm::EcalClusterAnalyzer::inverse_skim_

private

Definition at line 64 of file EcalClusterAnalyzer.h.

mixed_hit_cutoff_

double dqm::EcalClusterAnalyzer::mixed_hit_cutoff_

private

Definition at line 61 of file EcalClusterAnalyzer.h.

n_ecal_clusters_min_

int dqm::EcalClusterAnalyzer::n_ecal_clusters_min_

private

Definition at line 67 of file EcalClusterAnalyzer.h.

nbr_of_electrons_

int dqm::EcalClusterAnalyzer::nbr_of_electrons_

private

What is the number of electrons in the event?

Definition at line 35 of file EcalClusterAnalyzer.h.

Referenced by analyze(), and configure().

rec_hit_coll_name_

std::string dqm::EcalClusterAnalyzer::rec_hit_coll_name_

private

Definition at line 44 of file EcalClusterAnalyzer.h.

rec_hit_pass_name_

std::string dqm::EcalClusterAnalyzer::rec_hit_pass_name_

private

Definition at line 47 of file EcalClusterAnalyzer.h.

use_simulated_electron_number_

bool dqm::EcalClusterAnalyzer::use_simulated_electron_number_

private

Use the number of simulated electrons instead of the number of determined by the TS track counting.

Definition at line 32 of file EcalClusterAnalyzer.h.

Referenced by analyze(), and configure().

---

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

• DQM/include/DQM/EcalClusterAnalyzer.h
• DQM/src/DQM/EcalClusterAnalyzer.cxx