Monopole_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       Monopole
// Module Type: analyzer
// File:        Monopole_module.cc
//
// Generated at Tue Jun  9 10:22:10 2015 by Martin Frank using artmod
// from cetpkgsupport v1_08_05.
////////////////////////////////////////////////////////////////////////

#include <art/Framework/Core/EDAnalyzer.h>
// #include <canvas/Persistency/Common/FindManyP.h>
#include "art/Framework/Core/FileBlock.h"
#include <art/Framework/Core/FindManyP.h>
#include <art/Framework/Core/ModuleMacros.h>
#include <art/Framework/Principal/Event.h>
#include <art/Framework/Principal/Handle.h>
#include <art/Framework/Principal/Run.h>
#include <art/Framework/Principal/SubRun.h>
#include <art/Framework/Services/Optional/TFileService.h>
#include <art/Persistency/Provenance/ProcessHistoryRegistry.h>
// #include <canvas/Utilities/InputTag.h>
#include <art/Utilities/InputTag.h>
#include <fhiclcpp/ParameterSet.h>
#include <fhiclcpp/ParameterSetRegistry.h>
#include <messagefacility/MessageLogger/MessageLogger.h>

#include <DDTBaseDataProducts/TriggerDecision.h>

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "Monopole/Constants.h"
#include "Monopole/Track3D.h"
#include "RawData/RawTrigger.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/HoughResult.h"
#include "Simulation/Particle.h"
// #include "nusimdata/SimulationBase/MCParticle.h"
// #include "nusimdata/SimulationBase/MCTruth.h"
#include "SimulationBase/MCParticle.h"
#include "SimulationBase/MCTruth.h"

#include <TH1.h>
#include <TTree.h>

#include <algorithm>
#include <bitset>
#include <map>
#include <vector>

namespace mono
{
  class Monopole;
}

class mono::Monopole : public art::EDAnalyzer
{
public:
  explicit Monopole(fhicl::ParameterSet const & p);

  Monopole(Monopole const &) = delete;
  Monopole(Monopole &&) = delete;
  Monopole & operator = (Monopole const &) = delete;
  Monopole & operator = (Monopole &&) = delete;

  void beginJob() override;
  void analyze(art::Event const & e) override;
  void respondToOpenInputFile(art::FileBlock const& fb) override;
  
private:

  int         extent_cell(std::vector<sim::FLSHit> const& hits) const;
  int         extent_plane(std::vector<sim::FLSHit> const& hits) const;
  TVector3    get_position(sim::FLSHit const& hit) const;
  unsigned    n_mc_hits(mono::Track3D const& track) const;
  unsigned    n_slices(std::vector<rb::Cluster> const& slices) const;
  double      phi(sim::Particle const* particle) const;
  void        split_by_view
              (std::map<std::string, std::vector<sim::FLSHit> > & hitmap) const;
  double      theta(sim::Particle const* particle) const;
  geo::View_t view(unsigned plane) const;
 
  
  void tset(std::string const& branch_name, double const& value);

  std::string current_file_name_;
  
  double branch_default_value_;
  unsigned n_tracks_to_record_;
  TTree *tree_;
  std::map<std::string, TH1*> h_;
  std::map<std::string, double> t_;
  art::ServiceHandle<art::TFileService> tfs_;
  art::ServiceHandle<cheat::BackTracker> bt_;
  art::ServiceHandle<geo::Geometry> geo_;
};



void mono::Monopole::tset(std::string const& branch_name, double const& value)
{
  auto it = t_.find(branch_name);
  if (it == t_.end())
    throw art::Exception(art::errors::UnimplementedFeature)
      << __FILE__ << ":" << __LINE__ << "\n"
      << "Branch " << branch_name << " does not exist!"
      << std::endl;
  
  it->second = value;
}



mono::Monopole::Monopole(fhicl::ParameterSet const & p) :
  EDAnalyzer(p),
  current_file_name_("not set"),
  branch_default_value_(-9999),
  n_tracks_to_record_(3)
{
}



void mono::Monopole::respondToOpenInputFile(art::FileBlock const& fb)
{
  current_file_name_ = fb.fileName();

  std::cout << "mono::Monopole::respondToOpenInputFile: "
            << "The current file name is "
            << current_file_name_
            << std::endl;
}



void mono::Monopole::beginJob()
{
  tree_ = tfs_->make<TTree>("Event", "Monopole Event Tree");

  std::vector<std::string> branch_names =
    { "run_number", "subrun_number", "event_number", "event_length", "n_hits",
      "n_clusters", "cluster_n_hits", "cluster_n_hits_x", "cluster_n_hits_y",
      "reco_n_tracks", "is_mc", "n_slices_monopole", "n_slices_higheremoved",
      "n_hough_peaks_x_max", "n_hough_peaks_y_max",
      "slice_sum_adc_max",
      "gen_dxdz", "gen_dydz", "gen_beta",
      "mc_dxdz", "mc_dydz", "mc_dxdt", "mc_dydt", "mc_dzdt",
      "mc_velocity", "mc_beta",
      "mc_length", "mc_phi", "mc_theta",
      "mc_center_average_edep", "mc_center_n_hits",
      "mc_near_end_average_edep", "mc_near_end_n_hits",
      "mc_far_end_average_edep", "mc_far_end_n_hits",
      "mc_monopole_hit_detector",
      "mc_dx", "mc_dy", "mc_dz",
      "mc_dplane_all", "mc_dplane_x", "mc_dplane_y",
      "mc_dcell_all", "mc_dcell_x", "mc_dcell_y",
      "mc_n_hits_all", "mc_n_hits_x", "mc_n_hits_y",
      "ddt_trigger_decisions",
      "n_hits_ge_100_ADC",
      "n_hits_ge_150_ADC",
      "mc_n_hits_all_without_noise",
      "mc_n_hits_ge_100_ADC",
      "mc_n_hits_ge_150_ADC",
      "mc_n_hits_in_slice4d",
      "mc_n_hits_in_slice4d_ge_100_ADC",
      "mc_n_hits_on_cosmic_track",
      "mc_n_hits_on_cosmic_track_ge_100_ADC",
      "mc_adc_total",
      "mc_adc_mean",
      "mc_adc_variance",
      "mc_adc_min",
      "mc_adc_max",
      "mc_dEdx",
      "mc_dEdx_Birks",
      "mc_n_cerenkov",
      "mc_path_length_sum"
    };

  std::vector<std::string> track_names =
    { "dxdz", "dydz", "dxdt", "dydt", "dzdt", "velocity", "beta",
      "dzdt_x", "dzdt_y", "r2_xt", "r2_yt", "r2_zt_x", "r2_zt_y",
      "n_hits", "n_hits_x", "n_hits_y", "n_mc_hits", "sum_adc",
      "dt", "dx", "dy", "dz", "length", "length_xz", "length_yz",
      "tstart", "tend",
      "max_time_gap_xz", "max_time_gap_yz",
      "most_missing_planes_xz", "most_missing_planes_yz",
      "dcell_x", "dcell_y", "dplane_x", "dplane_y"
    };

  for (auto const& tname : track_names)
    for (unsigned n_track = 1; n_track <= n_tracks_to_record_; ++n_track)
      branch_names.
        push_back("reco_track_" + std::to_string(n_track) + "_" + tname);
  
  for (auto const& bname : branch_names)
    tree_->Branch(bname.c_str(), &t_[bname], (bname + "/D").c_str());

  tree_->Branch("input_file_name", &current_file_name_);


  h_["adc_mc"] =
    tfs_->make<TH1F>("adc_mc", "ADC Distribution - MC", 4096, -0.5, 4095.5);
  h_["adc_data"] =
    tfs_->make<TH1F>("adc_data", "ADC Distribution - Data", 4096, -0.5, 4095.5);


  // 64,000 TDC = 1000 us
  h_["tdc_mc"] =
    tfs_->make<TH1F>("tdc_mc", "TDC Distribution - MC", 1000, 0, 64000);

  h_["dEdx"] =
    tfs_->make<TH1F>("dEdx", "dE/dx Distribution;dE/dx [MeV/cm]", 200, 0, 40);
}



void mono::Monopole::analyze(art::Event const & e)
{
  //
  // Extract Objects from the Event
  //
  // Raw Trigger Header
  art::Handle<std::vector<rawdata::RawTrigger> > raw_triggers;
  e.getByLabel("daq", raw_triggers);
  if (raw_triggers->size() != 1)
    throw art::Exception(art::errors::DataCorruption)
      << __FILE__ << ":" << __LINE__ << "\n"
      << "RawTrigger vector has incorrect size: " << raw_triggers->size()
      << std::endl;
  unsigned long long
    event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
  uint32_t event_length(raw_triggers->front().fTriggerRange_TriggerLength);

  // Offline Cell Hits
  art::Handle<std::vector<rb::CellHit> > offline_hits;
  e.getByLabel("calhit", offline_hits);
  std::vector<art::Ptr<rb::CellHit> > offline_hitptrs;
  for (unsigned i = 0; i < offline_hits->size(); ++i)
    offline_hitptrs.emplace_back(offline_hits, i);

  /*
    EXCLUDE EXTRA INFO PARSING FOR NOW

  uint32_t extra_info(raw_triggers->front().fTriggerHeader_SourceSubID);
  std::bitset<novaddt::smt::Extra_Info::BITSET_SIZE> extra_bits(extra_info);

  std::cout << "MF: extra_info = " << extra_info << std::endl;
  std::cout << "MF: extra_bits = " << extra_bits << std::endl;
  if (extra_bits.test(novaddt::smt::Extra_Info::INVALID_OLD_DEFAULT_BIT))
    std::cout << "MF: extra_info is invalid (i.e. set to old default)"
              << std::endl;
  else
    std::cout << "MF: (Check, Timeout, Invalid) = (" << std::boolalpha
              << extra_bits.test(novaddt::smt::Extra_Info::TIMEOUT_BIT)
              << ", " << extra_bits.test(novaddt::smt::Extra_Info::CHECK_BIT)
              << ", " << extra_bits.
      test(novaddt::smt::Extra_Info::INVALID_OLD_DEFAULT_BIT)
              << ")" << std::endl;
  */

  // DDT Decisions
  art::Handle<std::vector<novaddt::TriggerDecision> > ddt_decisions;
  e.getByLabel("slowmonopoletrigger", ddt_decisions);

  // NOvA Slicer (Slicer4D)
  art::Handle<std::vector<rb::Cluster> > nova_slices;
  e.getByLabel("slicer", nova_slices);
  
  // Cosmic Tracks
  art::Handle<std::vector<rb::Track> > cosmic_tracks;
  e.getByLabel("cosmictrack", cosmic_tracks);
  
  // Monopole Cluster
  art::Handle<std::vector<rb::Cluster> > mono_clusters;
  e.getByLabel("monocluster", mono_clusters);

  // Monopole Slicer
  art::Handle<std::vector<rb::Cluster> > mono_slices;
  e.getByLabel("monoslicer", mono_slices);

  // High Energy Removed Slices
  art::Handle<std::vector<rb::Cluster> > higheremoved_slices;
  e.getByLabel("higheremoval", higheremoved_slices);

  // Reconstructed Tracks
  art::Handle<std::vector<mono::Track3D> > reco_tracks;
  e.getByLabel("monotrack", reco_tracks);

  // Generator Information
  art::Handle<std::vector<simb::MCTruth> > gen_info;
  e.getByLabel("generator", gen_info);


  //
  // Print Hits
  //
  // Monopole Track Hits
  // std::cout << "MF: N_tracks = " << reco_tracks->size() << std::endl;
  // if (reco_tracks->size() == 1)
  //   std::cout << reco_tracks->front() << std::endl;

  // // All Hits
  // art::ServiceHandle<geo::Geometry> geometry;
  // for (auto const& hit : offline_hitptrs)
  // {
  //   std::string view;
  //   if (hit->View() == geo::View_t::kX)
  //     view = "XZ";
  //   else if (hit->View() == geo::View_t::kY)
  //     view = "YZ";
  //   else
  //     view = "BAD";

  //   double xyz[3];
  //   geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
  //   TVector3 hit_vector(xyz);
  //   double v = geometry->CellTpos(hit->Plane(), hit->Cell());
  //   double z = hit_vector.z();
  //   double t = hit->TNS();
    
  //   if (z >   1050 && z <  1400 &&
  //    t > -16000 && t < -9000)
  //     std::cout << "AllHits: View " << view
  //            << " Hit: t, v, z  = " << t << ", " << v << ", " << z << "\n";
  // }


  
  //
  // Fill Histograms (and add some hit-based quantities together)
  //
  int n_hits_ge_100_ADC = 0;
  int n_hits_ge_150_ADC = 0;

  double mc_n_hits_all_without_noise = 0.0;
  int mc_n_hits_ge_100_ADC = 0;
  int mc_n_hits_ge_150_ADC = 0;

  double mc_adc_total = 0.0;
  double mc_adc_mean = 0.0;
  double mc_adc_variance = 0.0;
  double mc_adc_max = -9999;
  double mc_adc_min = -9999;
  if (!offline_hits.failedToGet())
  {
    std::vector<int32_t> mc_tdc_values;
    std::vector<double> mc_adc_values;
    for (auto const& hit : offline_hitptrs)
    {
      if (hit->ADC() >= 100)
        ++n_hits_ge_100_ADC;

      if (hit->ADC() >= 150)
        ++n_hits_ge_150_ADC;

      if (hit->IsMC())
      {
        if (!bt_->IsNoise(hit))
        {
          ++mc_n_hits_all_without_noise;
          
          h_.at("adc_mc")->Fill(hit->ADC());
          mc_adc_values.push_back(hit->ADC());
          mc_tdc_values.push_back(hit->TDC());

          if (hit->ADC() >= 100)
            ++mc_n_hits_ge_100_ADC;

          if (hit->ADC() >= 150)
            ++mc_n_hits_ge_150_ADC;
        }
      } else {
        h_.at("adc_data")->Fill(hit->ADC());
      } 
    }

    if (!mc_tdc_values.empty())
    {
      auto tdc_min =
        *std::min_element(mc_tdc_values.cbegin(), mc_tdc_values.cend());
      for (auto const& tdc : mc_tdc_values)
        h_.at("tdc_mc")->Fill(tdc - tdc_min);
    }

    if (mc_n_hits_all_without_noise > 1)
    {
      for (auto const& adc : mc_adc_values)
        mc_adc_total += adc;

      mc_adc_mean = mc_adc_total / mc_n_hits_all_without_noise;
      
      for (auto const& adc : mc_adc_values)
        mc_adc_variance += std::pow(adc - mc_adc_mean, 2);

      mc_adc_variance = std::sqrt(mc_adc_variance /
                                  (mc_n_hits_all_without_noise - 1));

      mc_adc_min =
        *std::min_element(mc_adc_values.cbegin(), mc_adc_values.cend());
      mc_adc_max =
        *std::max_element(mc_adc_values.cbegin(), mc_adc_values.cend());
    }
  }

  // Slicer4D Hit Counts
  int mc_n_hits_in_slice4d = 0;
  int mc_n_hits_in_slice4d_ge_100_ADC = 0;
  if (!nova_slices.failedToGet())
  {
    for (auto const& nova_slice : *nova_slices)
    {
      if (!nova_slice.IsNoise())
      {
        for (auto const& nova_slice_hit : nova_slice.AllCells())
        {
          if (nova_slice_hit->IsMC() &&
              !bt_->IsNoise(nova_slice_hit))
          {
            ++mc_n_hits_in_slice4d;

            if (nova_slice_hit->ADC() >= 100)
              ++mc_n_hits_in_slice4d_ge_100_ADC;
          }
        }
      }
    }
  }

  // Cosmic Track Hit Counts
  int mc_n_hits_on_cosmic_track = 0;
  int mc_n_hits_on_cosmic_track_ge_100_ADC = 0;
  if (!cosmic_tracks.failedToGet())
  {
    for (auto const& cosmic_track : *cosmic_tracks)
    {
      for (auto const& cosmic_track_hit : cosmic_track.AllCells())
      {
        if (cosmic_track_hit->IsMC() &&
            !bt_->IsNoise(cosmic_track_hit))
        {
          ++mc_n_hits_on_cosmic_track;

          if (cosmic_track_hit->ADC() >= 100)
            ++mc_n_hits_on_cosmic_track_ge_100_ADC;
        }
      }
    }
  }

  

  //
  // Fill Tree
  //
  for (auto & branch : t_)
    branch.second = branch_default_value_;

  // Event Information
  tset("run_number", e.run());
  tset("subrun_number", e.subRun());
  tset("event_number", e.event());
  tset("event_length", event_length);

  // DDT Information
  if (!ddt_decisions.failedToGet())
    tset("ddt_trigger_decisions", ddt_decisions->size());

  // Hit-Level Information (from histogram filling above)
  tset("n_hits_ge_100_ADC", n_hits_ge_100_ADC);
  tset("n_hits_ge_150_ADC", n_hits_ge_150_ADC);

  tset("mc_n_hits_all_without_noise", mc_n_hits_all_without_noise);
  tset("mc_n_hits_ge_100_ADC", mc_n_hits_ge_100_ADC);
  tset("mc_n_hits_ge_150_ADC", mc_n_hits_ge_150_ADC);

  tset("mc_adc_total", mc_adc_total);
  tset("mc_adc_mean", mc_adc_mean);
  tset("mc_adc_variance", mc_adc_variance);
  tset("mc_adc_min", mc_adc_min);
  tset("mc_adc_max", mc_adc_max);

  // NOvA Slice (Slicer4D) Information
  tset("mc_n_hits_in_slice4d", mc_n_hits_in_slice4d);
  tset("mc_n_hits_in_slice4d_ge_100_ADC", mc_n_hits_in_slice4d_ge_100_ADC);

  // Cosmic Track Information
  tset("mc_n_hits_on_cosmic_track", mc_n_hits_on_cosmic_track);
  tset("mc_n_hits_on_cosmic_track_ge_100_ADC",
       mc_n_hits_on_cosmic_track_ge_100_ADC);
  
  // Reconstruction Information
  if (!offline_hits.failedToGet())
    tset("n_hits", offline_hits->size());

  if (!mono_clusters.failedToGet())
  {
    tset("n_clusters", mono_clusters->size());

    if (mono_clusters->size() == 1)
    {
      auto cluster = mono_clusters->front();
      tset("cluster_n_hits", cluster.NCell());
      tset("cluster_n_hits_x", cluster.NXCell());
      tset("cluster_n_hits_y", cluster.NYCell());
    }
  }

  if (!mono_slices.failedToGet())
  {
    tset("n_slices_monopole", n_slices(*mono_slices));

    auto max_it = std::max_element
      (mono_slices->begin(), mono_slices->end(),
       [](auto const lhs, auto const rhs)
       { return lhs.TotalADC() < rhs.TotalADC(); });
    tset("slice_sum_adc_max", max_it->TotalADC());
  }

  if (!higheremoved_slices.failedToGet())
  {
    tset("n_slices_higheremoved", n_slices(*higheremoved_slices));

    // Retrieve Hough information from the event:
    art::FindManyP<rb::HoughResult>
    find_hough_results(higheremoved_slices, e, "hough");
    
    // Record the maximum number of Hough peaks:
    std::map<geo::View_t, unsigned> n_peaks_max =
      { {geo::View_t::kX, 0},
        {geo::View_t::kY, 0} };

    for (size_t n_slice = 0; n_slice != higheremoved_slices->size(); ++n_slice)
    {
      if (!higheremoved_slices->at(n_slice).IsNoise())
      {
        std::vector<art::Ptr<rb::HoughResult> > hough_results;
        find_hough_results.get(n_slice, hough_results);

        for (auto const& hough_result : hough_results)
        {
          auto view = hough_result->fView;
          unsigned n_peak = hough_result->fPeak.size();
          n_peaks_max.at(view) = std::max(n_peak, n_peaks_max.at(view));
        }
      }
    }
    tset("n_hough_peaks_x_max", n_peaks_max.at(geo::View_t::kX));
    tset("n_hough_peaks_y_max", n_peaks_max.at(geo::View_t::kY));
  }
  
  if (!reco_tracks.failedToGet())
  {
    tset("reco_n_tracks", reco_tracks->size());

    unsigned n_tracks =
      std::min(n_tracks_to_record_, unsigned(reco_tracks->size()));
    for (unsigned n_track = 1; n_track <= n_tracks; ++n_track)
    {
      auto track = reco_tracks->at(n_track - 1);

      std::string name = "reco_track_" + std::to_string(n_track) + "_";
      tset(name + "dxdz", track.dxdz());
      tset(name + "dydz", track.dydz());

      tset(name + "dxdt", track.velocity().x());
      tset(name + "dydt", track.velocity().y());
      tset(name + "dzdt", track.velocity().z());
      tset(name + "velocity", track.velocity().Mag());
      tset(name + "beta", track.beta());
    
      tset(name + "dzdt_x", track.dzdt(geo::View_t::kX));
      tset(name + "dzdt_y", track.dzdt(geo::View_t::kY));

      // squared linear correlation factors:
      tset(name + "r2_xt", track.r2_xt());
      tset(name + "r2_yt", track.r2_yt());
      tset(name + "r2_zt_x", track.r2_zt(geo::View_t::kX));
      tset(name + "r2_zt_y", track.r2_zt(geo::View_t::kY));
      
      tset(name + "n_hits", track.cluster().NCell());
      tset(name + "n_hits_x", track.cluster().NXCell());
      tset(name + "n_hits_y", track.cluster().NYCell());
      tset(name + "n_mc_hits", n_mc_hits(track));
      
      tset(name + "sum_adc", track.cluster().TotalADC());
      
      tset(name + "tstart", track.cluster().MinTNS());
      tset(name + "tend", track.cluster().MaxTNS());

      tset(name + "dcell_x", track.cluster().ExtentCell(geo::View_t::kX));
      tset(name + "dcell_y", track.cluster().ExtentCell(geo::View_t::kY));

      tset(name + "dplane_x", track.cluster().ExtentPlane(geo::View_t::kX));
      tset(name + "dplane_y", track.cluster().ExtentPlane(geo::View_t::kY));
      
      TVector3 ext = track.cluster().ExtentXYZ();
      tset(name + "dx", ext.x());
      tset(name + "dy", ext.y());
      tset(name + "dz", ext.z());
      tset(name + "length", ext.Mag());
      tset(name + "length_xz", util::pythag(ext.x(), ext.z()));
      tset(name + "length_yz", util::pythag(ext.y(), ext.z()));
      tset(name + "dt", track.cluster().ExtentTNS());

      tset(name + "max_time_gap_xz", track.max_time_gap_xz());
      tset(name + "max_time_gap_yz", track.max_time_gap_yz());

      tset(name + "most_missing_planes_xz",
           track.cluster().MostMissingPlanes(geo::View_t::kX));
      tset(name + "most_missing_planes_yz",
           track.cluster().MostMissingPlanes(geo::View_t::kY));
    }
  }

  if (bt_->HaveTruthInfo())
  {
    tset("is_mc", true);
    tset("mc_monopole_hit_detector", false);
    
    auto c = Constants::SPEED_OF_LIGHT;
    
    // Generator Information
    /*
      While working with the generator information, I found that its
      information is incorrect (e.g. m = -2e9 GeV/c^2 for monopole).
      The information extracted from the particle navigator, however,
      is correct.

      TO DO:
      1. since the generator mass is incorrect, implement it manually:
      const double MONOPOLE_MASS =  1e16; // GeV / c^2

      2. use the hard-coded mass to calculate all of the generator quantities

      3. add all of the MC quantities at the generator level

      4. add the mass, px, py, pz, and p to MC and Gen level information
         (makes it easy to calculate all quantities with the tree)
    */
    if (gen_info->size() == 1)
    {
      simb::MCParticle const& particle = gen_info->front().GetParticle(0);
      unsigned const trajectory = 0;

      auto m = particle.Mass();
      auto p = particle.Momentum(trajectory);

      tset("gen_dxdz", p.Px() / p.Pz());
      tset("gen_dydz", p.Py() / p.Pz());
      tset("gen_beta", p.P()  / m);
    } else {
      throw art::Exception(art::errors::LogicError)
        << __FILE__ << ":" << __LINE__ << "\n"
        << "The generator object contains " << gen_info->size()
        << " elements, there should only be 1." << std::endl;
    }  
    
    // MC Information
    sim::ParticleNavigator const& pn = bt_->ParticleNavigator();
    std::cout << "MF: Number of MC Particles in Navigator = "
              << pn.size() << std::endl;
    if (!pn.empty())
    {
      tset("mc_monopole_hit_detector", true);

      sim::Particle const* particle = pn.Primary(0);
      unsigned const trajectory = 0;

      auto m = particle->Mass();
      auto p = particle->Momentum(trajectory);
      
      tset("mc_phi", phi(particle));
      tset("mc_theta", theta(particle));

      tset("mc_dxdz", p.Px() / p.Pz());
      tset("mc_dydz", p.Py() / p.Pz());

      tset("mc_dxdt",     c * p.Px() / m);
      tset("mc_dydt",     c * p.Py() / m);
      tset("mc_dzdt",     c * p.Pz() / m);
      tset("mc_velocity", c * p.P()  / m);
      tset("mc_beta",         p.P()  / m);

      std::map<std::string, std::vector<sim::FLSHit> > fls_hits;
      fls_hits["all"] = bt_->ParticleToFLSHit(particle->TrackId());
      split_by_view(fls_hits);

      // dE/dx
      double Edep_sum        = 0.0;
      double Edep_Birks_sum  = 0.0;
      double path_length_sum = 0.0;
      double n_cerenkov_sum  = 0.0;
      for (auto hit : fls_hits["all"])
      {
        double path_length = hit.GetTotalPathLength();
        if (path_length > 5.0)
        {
          Edep_sum        += hit.GetEdep();
          Edep_Birks_sum  += hit.GetEdepBirks();
          path_length_sum += path_length;

          h_.at("dEdx")->Fill(hit.GetEdep() / path_length * 1e3); // MeV/cm
        }
      }

      tset("mc_path_length_sum", path_length_sum); // cm
      if (path_length_sum > 0)
      {
        tset("mc_dEdx",       Edep_sum       / path_length_sum); // GeV/cm
        tset("mc_dEdx_Birks", Edep_Birks_sum / path_length_sum); // GeV/cm
        tset("mc_n_cerenkov", n_cerenkov_sum / path_length_sum);
      }

      for (auto const& hits : fls_hits)
      {
        std::string suffix = hits.first;

        tset("mc_n_hits_" + suffix, hits.second.size());
        tset("mc_dplane_" + suffix, extent_plane(hits.second));
        tset("mc_dcell_" + suffix, extent_cell(hits.second));
      }      
      
      std::sort(fls_hits.at("all").begin(), fls_hits.at("all").end(),
                [](sim::FLSHit const& lhs, sim::FLSHit const& rhs)
                { return lhs.GetEntryT() < rhs.GetEntryT(); });

      auto entry_hit = fls_hits.at("all").front();
      auto exit_hit = fls_hits.at("all").back();
      TVector3 entry_position = get_position(entry_hit);
      TVector3 exit_position = get_position(exit_hit);
      TVector3 track = exit_position - entry_position;
      tset("mc_dx", track.X());
      tset("mc_dy", track.Y());
      tset("mc_dz", track.Z());
      tset("mc_length", track.Mag());
    }

    // Energy Deposition
    std::map<std::string, std::vector<double> > edeps;
    for (auto const& hit : offline_hitptrs)
    { 
      auto fls_hits = bt_->HitToFLSHit(hit);
      for (auto const& fls_hit : fls_hits)
      {
        double edep =
          static_cast<double>(hit->ADC()) / fls_hit.GetTotalPathLength();

        double z = (fls_hit.GetEntryZ() + fls_hit.GetExitZ()) / 2.0;
        if (-25.0 < z && z < 25.0)
          edeps["center"].push_back(edep);
        else if (-650.0 < z && z < -600.0)
          edeps["far_end"].push_back(edep);
        else if (600.0 < z && z < 650)
          edeps["near_end"].push_back(edep);
      }
    }

    for (auto const& position : edeps)
    {
      double sum = 0.0;
      for (auto const& value : position.second)
        sum += value;

      double average_edep = sum / static_cast<double>(position.second.size());
      tset("mc_" + position.first + "_average_edep", average_edep);
      tset("mc_" + position.first + "_n_hits", position.second.size());
    }
  } else {
    tset("is_mc", false);
  }
  
  tree_->Fill();
}



int mono::Monopole::extent_cell
(std::vector<sim::FLSHit> const& hits) const
{
  if (hits.empty())
    return branch_default_value_;

  auto it = std::minmax_element
    (hits.cbegin(), hits.cend(),
     [](auto const lhs, auto const rhs)
     { return lhs.GetCellID() < rhs.GetCellID(); });

  // we add 1 to be consistent with rb::Cluster::ExtentCell()
  return std::abs(it.first->GetCellID() - it.second->GetCellID()) + 1;
}
 
 

int mono::Monopole::extent_plane
(std::vector<sim::FLSHit> const& hits) const
{
  if (hits.empty())
    return branch_default_value_;

  auto it = std::minmax_element
    (hits.cbegin(), hits.cend(),
     [](auto const lhs, auto const rhs)
     { return lhs.GetPlaneID() < rhs.GetPlaneID(); });

  // we add 1 to be consistent with rb::Cluster::ExtentPlane()
  return std::abs(it.first->GetPlaneID() - it.second->GetPlaneID()) + 1;
}



TVector3 mono::Monopole::get_position(sim::FLSHit const& hit) const
{
  TVector3 result;
  geo_->Plane(hit.GetPlaneID())->Cell(hit.GetCellID())->
    GetCenter(result, hit.GetZAverage());
  
  return result;
}



unsigned mono::Monopole::n_mc_hits(mono::Track3D const& track) const
{
  unsigned result = 0;

  if (!bt_->HaveTruthInfo())
    return result;

  for (auto const& hit : track.cluster().AllCells())
    if (hit->IsMC())
      ++result;

  return result;
}



unsigned mono::Monopole::n_slices(std::vector<rb::Cluster> const& slices) const
{
  unsigned result = std::count_if
    (slices.begin(), slices.end(),
     [](rb::Cluster const& slice){ return !slice.IsNoise(); });
  
  return result;
}



double mono::Monopole::phi(sim::Particle const* particle) const
{
  double Px(particle->Px()), Py(particle->Py());

  return std::atan2(Py, Px);
}



void mono::Monopole::split_by_view
(std::map<std::string, std::vector<sim::FLSHit> > & hitmap) const
{
  hitmap["x"].clear();
  hitmap["y"].clear();
  
  for (auto const & hit : hitmap.at("all"))
  {
    if (view(hit.GetPlaneID()) == geo::View_t::kX)
      hitmap.at("x").push_back(hit);
    else
      hitmap.at("y").push_back(hit);
  }
}



double mono::Monopole::theta(sim::Particle const* particle) const
{
  double Px(particle->Px()), Py(particle->Py()), Pz(particle->Pz());
  double Pt = std::sqrt(Px * Px + Py * Py);

  return std::atan2(Pt, Pz);
}



geo::View_t mono::Monopole::view(unsigned plane) const
{
  // x planes are odd, y planes are even
  if (plane % 2 == 1)
    return geo::View_t::kX;

  return geo::View_t::kY;
}



DEFINE_ART_MODULE(mono::Monopole)