MCMatchAlg.cxx

Go to the documentation of this file.

#include "MCMatchAlg.h"

#include "TString.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "larcore/Geometry/WireReadout.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larreco/MCComp/MCBTAlg.h"
#include "larreco/MCComp/MCBTException.h"

#include <string>

namespace btutil {

  MCMatchAlg::MCMatchAlg()
  {
    _view_to_plane.clear();
  }

  bool MCMatchAlg::BuildMap(detinfo::DetectorClocksData const& clockData,
                            const std::vector<unsigned int>& g4_trackid_v,
                            const std::vector<sim::SimChannel>& simch_v,
                            const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
  {
    fBTAlgo.Reset(g4_trackid_v, simch_v);

    return BuildMap(clockData, cluster_v);
  }

  bool MCMatchAlg::BuildMap(detinfo::DetectorClocksData const& clockData,
                            const std::vector<std::vector<unsigned int>>& g4_trackid_v,
                            const std::vector<sim::SimChannel>& simch_v,
                            const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
  {

    fBTAlgo.Reset(g4_trackid_v, simch_v);

    return BuildMap(clockData, cluster_v);
  }

  bool MCMatchAlg::BuildMap(detinfo::DetectorClocksData const& clockData,
                            const std::vector<std::vector<art::Ptr<recob::Hit>>>& cluster_v)
  {
    size_t num_mcobj = fBTAlgo.NumParts();
    size_t num_cluster = cluster_v.size();

    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();

    //
    // Perform back-tracking
    //
    // (1) Load cluster/hit data product
    // (2) Loop over all clusters and find charge fraction to each MC object
    // (3) For each MC object, find a cluster per plane with the highest charge

    // Loop over clusters & get charge fraction per MC object
    _summed_mcq.clear();
    _cluster_mcq_v.clear();
    _cluster_plane_id.clear();

    _summed_mcq.resize(num_mcobj + 1, std::vector<double>(wireReadoutGeom.Nplanes(), 0));
    _cluster_mcq_v.reserve(num_cluster);
    _cluster_plane_id.reserve(num_cluster);

    for (auto const& hit_v : cluster_v) {

      size_t plane = wireReadoutGeom.Nplanes();

      // Create hit list
      std::vector<WireRange_t> wr_v;
      wr_v.reserve(hit_v.size());

      for (auto const& h : hit_v) {

        WireRange_t wr;
        wr.ch = h->Channel();
        wr.start = h->StartTick();
        wr.end = h->EndTick();
        wr_v.push_back(wr);
        if (plane == wireReadoutGeom.Nplanes()) plane = h->WireID().Plane;
      }

      _cluster_plane_id.push_back(plane);

      auto mcq_v = fBTAlgo.MCQ(clockData, wr_v);

      for (size_t i = 0; i < mcq_v.size(); ++i)

        _summed_mcq[i][plane] += mcq_v[i];

      _cluster_mcq_v.push_back(mcq_v);
    }

    //
    // Find the best matched pair (and its charge) per MC object
    //
    _bmatch_id.clear();
    _bmatch_id.resize(num_mcobj, std::vector<int>(wireReadoutGeom.Nplanes(), -1));

    std::vector<std::vector<double>> bmatch_mcq(num_mcobj,
                                                std::vector<double>(wireReadoutGeom.Nplanes(), 0));

    for (size_t mc_index = 0; mc_index < num_mcobj; ++mc_index) {

      for (size_t cluster_index = 0; cluster_index < num_cluster; ++cluster_index) {

        auto plane = _cluster_plane_id.at(cluster_index);

        double q = _cluster_mcq_v[cluster_index][mc_index];

        if (bmatch_mcq[mc_index][plane] < q) {

          bmatch_mcq[mc_index][plane] = q;
          _bmatch_id[mc_index][plane] = cluster_index;
        }
      }
    }

    return true;
  }

  double MCMatchAlg::ClusterCorrectness(const size_t cluster_index, const size_t mc_index) const
  {

    if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");

    if (cluster_index >= _cluster_mcq_v.size())
      throw MCBTException(Form(
        "Input cluster index (%zu) out of range (%zu)!", cluster_index, _cluster_mcq_v.size()));

    if (mc_index >= _bmatch_id.size())
      throw MCBTException(
        Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));

    auto plane = _cluster_plane_id.at(cluster_index);

    auto best_cluster_index = _bmatch_id.at(mc_index).at(plane);

    if (best_cluster_index < 0) return -1;

    return _cluster_mcq_v.at(cluster_index).at(mc_index) /
           _cluster_mcq_v.at(best_cluster_index).at(mc_index);
  }

  std::pair<size_t, double> MCMatchAlg::ShowerCorrectness(
    const std::vector<unsigned int> cluster_indices) const
  {

    if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");

    if (cluster_indices.size() > _cluster_mcq_v.size())
      throw MCBTException(
        Form("Input cluster indices length (%zu) > # of registered clusters (%zu)!",
             cluster_indices.size(),
             _cluster_mcq_v.size()));

    if (!cluster_indices.size()) throw MCBTException("Input cluster indices empty!");

    // Compute efficiency per MC
    std::vector<double> match_eff(_bmatch_id.size(), 1);

    for (auto const& cluster_index : cluster_indices) {

      for (size_t mc_index = 0; mc_index < _bmatch_id.size(); ++mc_index) {

        double correctness = ClusterCorrectness(cluster_index, mc_index);

        if (correctness >= 0) match_eff.at(mc_index) *= correctness;
      }
    }

    std::pair<size_t, double> result(0, -1);

    // Find the best qratio
    for (size_t mc_index = 0; mc_index < match_eff.size(); ++mc_index) {

      if (match_eff.at(mc_index) < result.second) continue;

      result.second = match_eff.at(mc_index);

      result.first = mc_index;
    }
    return result;
  }

  std::pair<double, double> MCMatchAlg::ClusterEP(const size_t cluster_index,
                                                  const size_t mc_index) const
  {
    if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");

    if (cluster_index >= _cluster_mcq_v.size())
      throw MCBTException(Form(
        "Input cluster index (%zu) out of range (%zu)!", cluster_index, _cluster_mcq_v.size()));

    if (mc_index >= _bmatch_id.size())
      throw MCBTException(
        Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));

    // efficiency = this cluster's mcq for specified mc obj / total mcq by this mcshower in all clusters
    // purity = this cluster's mcq for this mcshower / total mcq from all mc obj in this cluster

    std::pair<double, double> result;

    auto plane = _cluster_plane_id[cluster_index];

    result.first = _cluster_mcq_v[cluster_index][mc_index] / _summed_mcq[mc_index][plane];

    double cluster_mcq_total = 0;
    for (auto const& q : _cluster_mcq_v[cluster_index])
      cluster_mcq_total += q;

    result.second = _cluster_mcq_v[cluster_index][mc_index] / cluster_mcq_total;

    return result;
  }

  const std::vector<int>& MCMatchAlg::BestClusters(const size_t mc_index) const
  {
    if (!_bmatch_id.size()) throw MCBTException("Preparation not done yet!");

    if (mc_index >= _bmatch_id.size())
      throw MCBTException(
        Form("Input MC object index (%zu) out of range (%zu)!", mc_index, _bmatch_id.size()));

    return _bmatch_id[mc_index];
  }

  std::pair<double, double> MCMatchAlg::BestClusterEP(const size_t mc_index,
                                                      const size_t plane_id) const
  {

    auto c_index_v = BestClusters(mc_index);

    if (c_index_v.size() <= plane_id)
      throw MCBTException(Form(
        "Plane ID %zu exceeds # of planes recorded in data (%zu)...", plane_id, c_index_v.size()));

    std::pair<double, double> result(0, 0);

    if (c_index_v[plane_id] < 0) return result;

    return ClusterEP(c_index_v[plane_id], mc_index);
  }

}