MRCCAna_module.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
// MRCCAna_module.cc
// Mike Wallbank, University of Cincinnati <wallbank@fnal.gov>, July 2020
//
// Class: MRCCAna
// Type: analyzer
// Description: MRCC validation module to analyze the hits which were
//              removed by the MuonRemoval process, in particular with
//              regards to NC interactions.
/////////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "fhiclcpp/ParameterSet.h"

// nova
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/FitSum.h"
#include "RecoBase/FilterList.h"

// stl
#include <string>
#include <vector>
#include <map>

// ROOT
#include "TH1I.h"
#include "TH1F.h"
#include "TProfile.h"

// --------------------------------------------------------------------------
namespace murem {
  class MRCCNeutrino;
  class MRCCParticle;
  class MRCCAna;
}

// --------------------------------------------------------------------------
class murem::MRCCNeutrino {
public:

  MRCCNeutrino(art::Ptr<simb::MCTruth> neutrino_int,
               const std::vector<const sim::Particle*>& particles);

  void AddAssociatedCluster(unsigned int clusterID);
  void AddAssociatedMRCCCluster(unsigned int clusterID);
  void AddParticleHit(unsigned int trackID, art::Ptr<rb::CellHit> hit);
  void AddParticleMRCCHit(unsigned int trackID, art::Ptr<rb::CellHit> hit);

  float Energy() const;
  bool CC() const;
  float ClusterCharge() const;
  float ClusterMRCCCharge() const;
  float ClusterHadronicCharge() const;
  float ClusterMRCCHadronicCharge() const;
  float ClusterLeptonicCharge() const;
  float ClusterMRCCLeptonicCharge() const;

  unsigned int NumAssociatedClusters() const;
  unsigned int NumAssociatedMRCCClusters() const;

  art::Ptr<simb::MCTruth> Truth() const;
  std::map<unsigned int, MRCCParticle*> Particles() const;
  int Muon() const;

private:

  art::Ptr<simb::MCTruth> fTruth;

  std::map<unsigned int, MRCCParticle*> fParticles;
  unsigned int* fMuon;

  std::vector<unsigned int> fAssociatedClusters;
  std::vector<unsigned int> fAssociatedMRCCClusters;

};

// --------------------------------------------------------------------------
class murem::MRCCParticle {
public:

  MRCCParticle(const simb::MCParticle* particle);

  void AddHit(art::Ptr<rb::CellHit> hit);
  void AddMRCCHit(art::Ptr<rb::CellHit> hit);

  float Energy() const;
  float DepositedCharge() const;
  float DepositedMRCCCharge() const;

  const simb::MCParticle* TrueParticle() const;

private:

  float fTotalCharge;
  float fTotalMRCCCharge;

  const simb::MCParticle* fTrueParticle;
  std::vector<art::Ptr<rb::CellHit> > fHits;
  std::vector<art::Ptr<rb::CellHit> > fMRCCHits;

};

// --------------------------------------------------------------------------
murem::MRCCNeutrino::MRCCNeutrino(art::Ptr<simb::MCTruth> neutrino_int,
                                  const std::vector<const sim::Particle*>& particles) {
  fTruth = neutrino_int;
  fMuon = nullptr;
  for (std::vector<const sim::Particle*>::const_iterator particleIt = particles.begin();
       particleIt != particles.end(); ++particleIt) {
    fParticles[(*particleIt)->TrackId()] = new MRCCParticle(*particleIt);
    if (!neutrino_int->GetNeutrino().CCNC() and !fMuon and (*particleIt)->PdgCode() == 13)
      fMuon = new unsigned int((*particleIt)->TrackId());
  }
}

// --------------------------------------------------------------------------
void murem::MRCCNeutrino::AddAssociatedCluster(unsigned int clusterID) {
  fAssociatedClusters.push_back(clusterID);
  return;
}

// --------------------------------------------------------------------------
void murem::MRCCNeutrino::AddAssociatedMRCCCluster(unsigned int clusterID) {
  fAssociatedMRCCClusters.push_back(clusterID);
  return;
}

// --------------------------------------------------------------------------
void murem::MRCCNeutrino::AddParticleHit(unsigned int trackID, art::Ptr<rb::CellHit> hit) {
  fParticles[trackID]->AddHit(hit);
  return;
}

// --------------------------------------------------------------------------
void murem::MRCCNeutrino::AddParticleMRCCHit(unsigned int trackID, art::Ptr<rb::CellHit> hit) {
  fParticles[trackID]->AddMRCCHit(hit);
  return;
}

// --------------------------------------------------------------------------
unsigned int murem::MRCCNeutrino::NumAssociatedClusters() const {
  return fAssociatedClusters.size();
}

// --------------------------------------------------------------------------
unsigned int murem::MRCCNeutrino::NumAssociatedMRCCClusters() const {
  return fAssociatedMRCCClusters.size();
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::Energy() const {
  return fTruth->GetNeutrino().Nu().E();
}

// --------------------------------------------------------------------------
bool murem::MRCCNeutrino::CC() const {
  return !fTruth->GetNeutrino().CCNC();
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterCharge() const {
  float charge = 0.;
  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
       particleIt != fParticles.end(); ++particleIt)
    charge += particleIt->second->DepositedCharge();
  return charge;
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterMRCCCharge() const {
  float charge = 0.;
  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
       particleIt != fParticles.end(); ++particleIt)
    charge += particleIt->second->DepositedMRCCCharge();
  return charge;
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterHadronicCharge() const {
  float charge = 0.;
  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
       particleIt != fParticles.end(); ++particleIt) {
    if (fMuon and particleIt->first == *fMuon)
      continue;
    charge += particleIt->second->DepositedCharge();
  }
  return charge;
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterMRCCHadronicCharge() const {
  float charge = 0.;
  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
       particleIt != fParticles.end(); ++particleIt) {
    if (fMuon and particleIt->first == *fMuon)
      continue;
    charge += particleIt->second->DepositedMRCCCharge();
  }
  return charge;
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterLeptonicCharge() const {
  float charge = 0.;
  if (fMuon)
    charge = fParticles.at(*fMuon)->DepositedCharge();
  return charge;
}

// --------------------------------------------------------------------------
float murem::MRCCNeutrino::ClusterMRCCLeptonicCharge() const {
  float charge = 0.;
  if (fMuon)
    charge = fParticles.at(*fMuon)->DepositedMRCCCharge();
  return charge;
}

// --------------------------------------------------------------------------
art::Ptr<simb::MCTruth> murem::MRCCNeutrino::Truth() const {
  return fTruth;
}

// --------------------------------------------------------------------------
std::map<unsigned int, murem::MRCCParticle*> murem::MRCCNeutrino::Particles() const {
  return fParticles;
}

// --------------------------------------------------------------------------
int murem::MRCCNeutrino::Muon() const {
  if (fMuon) return *fMuon;
  return -1;
}

// --------------------------------------------------------------------------
murem::MRCCParticle::MRCCParticle(const simb::MCParticle* particle) :
  fTotalCharge(0.),
  fTotalMRCCCharge(0.),
  fTrueParticle(particle) {
}

// --------------------------------------------------------------------------
void murem::MRCCParticle::AddHit(art::Ptr<rb::CellHit> hit) {
  fHits.push_back(hit);
  fTotalCharge += hit->ADC();
  return;
}

// --------------------------------------------------------------------------
void murem::MRCCParticle::AddMRCCHit(art::Ptr<rb::CellHit> hit) {
  fMRCCHits.push_back(hit);
  fTotalMRCCCharge += hit->ADC();
  return;
}

// --------------------------------------------------------------------------
float murem::MRCCParticle::Energy() const {
  return fTrueParticle->E();
}

// --------------------------------------------------------------------------
float murem::MRCCParticle::DepositedCharge() const {
  return fTotalCharge;
}

// --------------------------------------------------------------------------
float murem::MRCCParticle::DepositedMRCCCharge() const {
  return fTotalMRCCCharge;
}

// --------------------------------------------------------------------------
const simb::MCParticle* murem::MRCCParticle::TrueParticle() const {
  return fTrueParticle;
}

// --------------------------------------------------------------------------
class murem::MRCCAna : public art::EDAnalyzer {
public:

  explicit MRCCAna(const fhicl::ParameterSet& pset);
  virtual ~MRCCAna();

  void analyze(const art::Event& evt);
  void beginJob();

private:

  // methods
  bool SliceSelection(const art::Event& evt, const art::Ptr<rb::Cluster>& slice);

  // labels
  std::string fTruthModuleLabel;
  std::string fHitsModuleLabel;
  std::string fClusterModuleLabel;
  std::string fElasticModuleLabel;
  std::string fFuzzy3DModuleLabel;
  std::string fBPFTrackModuleLabel;
  std::string fMRCCHitsModuleLabel;
  std::string fMRCCClusterModuleLabel;

  // services
  art::ServiceHandle<cheat::BackTracker> bt;
  art::ServiceHandle<art::TFileService> tfs;

  // data products
  TH1I* hNumNuClusters;
  TH1I* hNumNuClustersMRCC;
  TProfile* hSliceHadronicFrac;
  TProfile* hSliceHadronicFracMRCC;
  TProfile* hAnyMRCCNeutrinoEnergy;
  TProfile* hAnyMRCCMuonEnergy;
  TProfile* hSliceEFracRemoved;
  TProfile* hSliceHadEFracRemoved;
  TProfile* hSliceLepEFracRemoved;
  TH1F* hSliceDeltaE;
  TH1F* hSliceDeltaLepE;
  TH1F* hSliceDeltaHadE;
  TProfile* hLeptonEFracRemoved;
  TProfile* hHadronEFracRemoved;
  TProfile* hMRCCSliceEFracRemoved;
  TProfile* hMRCCSliceHadEFracRemoved;
  TProfile* hMRCCSliceLepEFracRemoved;
  TProfile* hMRCCLeptonEFracRemoved;
  TProfile* hMRCCHadronEFracRemoved;
  TH1D* hHadYMRCC;
  TH1D* hHadYNoMRCC;
  TProfile* hSliceEFracRemovedQE;
  TProfile* hSliceEFracRemovedDIS;
  TProfile* hSliceEFracRemovedRES;

};

DEFINE_ART_MODULE(murem::MRCCAna)

// --------------------------------------------------------------------------
murem::MRCCAna::MRCCAna(const fhicl::ParameterSet& pset) :
  EDAnalyzer(pset),
  fTruthModuleLabel(pset.get<std::string>("TruthModuleLabel")),
  fHitsModuleLabel(pset.get<std::string>("HitsModuleLabel")),
  fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel")),
  fElasticModuleLabel(pset.get<std::string>("ElasticModuleLabel")),
  fFuzzy3DModuleLabel(pset.get<std::string>("Fuzzy3DModuleLabel")),
  fBPFTrackModuleLabel(pset.get<std::string>("BPFTrackModuleLabel")),
  fMRCCHitsModuleLabel(pset.get<std::string>("MRCCHitsModuleLabel")),
  fMRCCClusterModuleLabel(pset.get<std::string>("MRCCClusterModuleLabel")) {
}

// --------------------------------------------------------------------------
murem::MRCCAna::~MRCCAna() {
}

// --------------------------------------------------------------------------
void murem::MRCCAna::beginJob() {

  hNumNuClusters
    = tfs->make<TH1I>("NumNuClusters", ";Number of Clusters Associated with True Neutrino;", 5,0,5);
  hNumNuClusters->GetXaxis()->CenterLabels();
  hNumNuClustersMRCC
    = tfs->make<TH1I>("NumNuClustersMRCC", ";Number of Clusters Associated with True Neutrino;", 5,0,5);
  hNumNuClustersMRCC->GetXaxis()->CenterLabels();
  hSliceHadronicFrac
    = tfs->make<TProfile>("SliceHadronicFrac", ";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
  hSliceHadronicFrac->SetMinimum(0);
  hSliceHadronicFrac->SetMaximum(1);
  hSliceHadronicFracMRCC
    = tfs->make<TProfile>("SliceHadronicFracMRCC", ";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
  hSliceHadronicFracMRCC->SetMinimum(0);
  hSliceHadronicFracMRCC->SetMaximum(1);
  hAnyMRCCNeutrinoEnergy
    = tfs->make<TProfile>("AnyMRCCNeutrinoEnergy", ";Neutrino Energy (GeV);Fraction of Neutrino Interactions with Charge Removed;", 100,0,40);
  hAnyMRCCNeutrinoEnergy->SetMinimum(0);
  hAnyMRCCNeutrinoEnergy->SetMaximum(1);
  hAnyMRCCMuonEnergy
    = tfs->make<TProfile>("AnyMRCCMuonEnergy", ";Muon Energy (GeV);Fraction of CC Muons with Charge Removed;", 100,0,40);
  hAnyMRCCMuonEnergy->SetMinimum(0);
  hAnyMRCCMuonEnergy->SetMaximum(1);
  hSliceEFracRemoved
    = tfs->make<TProfile>("SliceEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
  hSliceEFracRemoved->SetMinimum(0);
  hSliceEFracRemoved->SetMaximum(1);
  hSliceHadEFracRemoved
    = tfs->make<TProfile>("SliceHadEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
  hSliceHadEFracRemoved->SetMinimum(0);
  hSliceHadEFracRemoved->SetMaximum(1);
  hSliceLepEFracRemoved
    = tfs->make<TProfile>("SliceLepEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
  hSliceLepEFracRemoved->SetMinimum(0);
  hSliceLepEFracRemoved->SetMaximum(1);
  hSliceDeltaE
    = tfs->make<TH1F>("SliceDeltaE", ";Fraction of Slice Charge Removed (GeV);", 100,0,1);
  hSliceDeltaLepE
    = tfs->make<TH1F>("SliceDeltaLepE", ";Fraction of Slice Leptonic Charge Removed (GeV);", 100,0,1);
  hSliceDeltaHadE
    = tfs->make<TH1F>("SliceDeltaHadE", ";Fraction of Slice Hadronic Charge Removed (GeV);", 100,0,1);
  hLeptonEFracRemoved
    = tfs->make<TProfile>("LeptonEFracRemoved", ";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
  hLeptonEFracRemoved->SetMinimum(0);
  hLeptonEFracRemoved->SetMaximum(1);
  hHadronEFracRemoved
    = tfs->make<TProfile>("HadronEFracRemoved", ";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
  hHadronEFracRemoved->SetMinimum(0);
  hHadronEFracRemoved->SetMaximum(1);
  hMRCCSliceEFracRemoved
    = tfs->make<TProfile>("MRCCSliceEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
  hMRCCSliceEFracRemoved->SetMinimum(0);
  hMRCCSliceEFracRemoved->SetMaximum(1);
  hMRCCSliceHadEFracRemoved
    = tfs->make<TProfile>("MRCCSliceHadEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
  hMRCCSliceHadEFracRemoved->SetMinimum(0);
  hMRCCSliceHadEFracRemoved->SetMaximum(1);
  hMRCCSliceLepEFracRemoved
    = tfs->make<TProfile>("MRCCSliceLepEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
  hMRCCSliceLepEFracRemoved->SetMinimum(0);
  hMRCCSliceLepEFracRemoved->SetMaximum(1);
  hMRCCLeptonEFracRemoved
    = tfs->make<TProfile>("MRCCLeptonEFracRemoved", ";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
  hMRCCLeptonEFracRemoved->SetMinimum(0);
  hMRCCLeptonEFracRemoved->SetMaximum(1);
  hMRCCHadronEFracRemoved
    = tfs->make<TProfile>("MRCCHadronEFracRemoved", ";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
  hMRCCHadronEFracRemoved->SetMinimum(0);
  hMRCCHadronEFracRemoved->SetMaximum(1);
  hHadYMRCC = tfs->make<TH1D>("HadYMRCC", ";Hadronic Y;", 100, 0, 1);
  hHadYNoMRCC = tfs->make<TH1D>("HadYNoMRCC", ";Hadronic Y;", 100, 0, 1);
  hSliceEFracRemovedQE
    = tfs->make<TProfile>("SliceEFracRemovedQE", ";Fraction of Slice Charge Removed (QE);", 100,0,40);
  hSliceEFracRemovedDIS
    = tfs->make<TProfile>("SliceEFracRemovedDIS", ";Fraction of Slice Charge Removed (DIS);", 100,0,40);
  hSliceEFracRemovedRES
    = tfs->make<TProfile>("SliceEFracRemovedRES", ";Fraction of Slice Charge Removed (RES);", 100,0,40);

  return;

}

// --------------------------------------------------------------------------
void murem::MRCCAna::analyze(const art::Event& evt) {

  // Get the truth information
  art::Handle<std::vector<simb::MCTruth> > truthHandle;
  std::vector<art::Ptr<simb::MCTruth> > truths;
  if (evt.getByLabel(fTruthModuleLabel, truthHandle))
    art::fill_ptr_vector(truths, truthHandle);

  // Get all the hits before muon removal
  art::Handle<std::vector<rb::CellHit> > hitsHandle;
  std::vector<art::Ptr<rb::CellHit> > hits;
  if (evt.getByLabel(fHitsModuleLabel, hitsHandle))
    art::fill_ptr_vector(hits, hitsHandle);

  // Get the clusters before muon removal
  art::Handle<std::vector<rb::Cluster> > clusterHandle;
  std::vector<art::Ptr<rb::Cluster> > clusters;
  if (evt.getByLabel(fClusterModuleLabel, clusterHandle))
    art::fill_ptr_vector(clusters, clusterHandle);

  // Get the hits following muon removal
  art::Handle<std::vector<rb::CellHit> > mrccHitsHandle;
  std::vector<art::Ptr<rb::CellHit> > mrccHits;
  if (evt.getByLabel(fMRCCHitsModuleLabel, mrccHitsHandle))
    art::fill_ptr_vector(mrccHits, mrccHitsHandle);

  // Get the clusters following muon removal
  art::Handle<std::vector<rb::Cluster> > mrccClusterHandle;
  std::vector<art::Ptr<rb::Cluster> > mrccClusters;
  if (evt.getByLabel(fMRCCClusterModuleLabel, mrccClusterHandle))
    art::fill_ptr_vector(mrccClusters, mrccClusterHandle);

  // Container to hold true neutrino information
  std::map<unsigned int, MRCCNeutrino*> true_neutrinos;

  // Grab and save the neutrino interaction truth information
  for (std::vector<art::Ptr<simb::MCTruth> >::const_iterator truthIt = truths.begin(); truthIt != truths.end(); ++truthIt) {
    std::vector<const sim::Particle*> sim_particles = bt->MCTruthToParticles(*truthIt);
    true_neutrinos[truthIt->key()] = new MRCCNeutrino(*truthIt, sim_particles);
  }

  // Add all hits to the relevant true particle for each neutrino interaction, before muon removal
  for (std::vector<art::Ptr<rb::Cluster> >::const_iterator clusterIt = clusters.begin();
       clusterIt != clusters.end(); ++clusterIt) {
    if ((*clusterIt)->IsNoise() or !SliceSelection(evt, *clusterIt))
      continue;
    std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, hits);
    if (nus.size() != 1)
      continue;
    true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedCluster(clusterIt->key());
    art::PtrVector<rb::CellHit> cluster_hits = (*clusterIt)->AllCells();
    for (art::PtrVector<rb::CellHit>::const_iterator hitIt = cluster_hits.begin();
         hitIt != cluster_hits.end(); ++hitIt) {
      const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
      if (true_particle == nullptr)
        continue;
      true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleHit(true_particle->TrackId(), *hitIt);
    }
  }

  // Add all hits to the relevant true particle for each neutrino interaction, after muon removal
  for (std::vector<art::Ptr<rb::Cluster> >::const_iterator clusterIt = mrccClusters.begin();
       clusterIt != mrccClusters.end(); ++clusterIt) {
    if ((*clusterIt)->IsNoise())// or !SliceSelection(evt, *clusterIt))
      continue;
    std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, mrccHits);
    if (nus.size() != 1)
      continue;
    true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedMRCCCluster(clusterIt->key());
    art::PtrVector<rb::CellHit> cluster_hits = (*clusterIt)->AllCells();
    for (art::PtrVector<rb::CellHit>::const_iterator hitIt = cluster_hits.begin();
         hitIt != cluster_hits.end(); ++hitIt) {
      const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
      if (true_particle == nullptr)
        continue;
      true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleMRCCHit(true_particle->TrackId(), *hitIt);
    }
  }

  // Plot things!
  for (std::map<unsigned int, MRCCNeutrino*>::const_iterator neutrinoIt = true_neutrinos.begin();
       neutrinoIt != true_neutrinos.end(); ++neutrinoIt) {

    const MRCCNeutrino* true_nu = neutrinoIt->second;

    // Neutrino-level plots
    hNumNuClusters->Fill(true_nu->NumAssociatedClusters());
    hNumNuClustersMRCC->Fill(true_nu->NumAssociatedMRCCClusters());

    // For the slice plots, only consider unique neutrino slices which weren't affected by the muon removal
    if (true_nu->NumAssociatedClusters() != 1 or true_nu->NumAssociatedClusters() != true_nu->NumAssociatedMRCCClusters())
      continue;

    // Slice-level plots
    hSliceHadronicFrac->Fill(true_nu->Energy(), true_nu->ClusterHadronicCharge()/true_nu->ClusterCharge());
    hSliceHadronicFracMRCC->Fill(true_nu->Energy(), true_nu->ClusterMRCCHadronicCharge()/true_nu->ClusterMRCCCharge());

    float charge_before = true_nu->ClusterCharge(), charge_after = true_nu->ClusterMRCCCharge();
    bool mrcc = charge_before != charge_after;

    if (true_nu->CC())
      hAnyMRCCNeutrinoEnergy->Fill(true_nu->Energy(), mrcc);

    if (charge_before) {
      hSliceEFracRemoved->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
      hSliceDeltaE->Fill((charge_before-charge_after)/charge_before);
      if (true_nu->Truth()->GetNeutrino().InteractionType() == 1001)
        hSliceEFracRemovedQE->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
      if (true_nu->Truth()->GetNeutrino().InteractionType() >= 1003 and true_nu->Truth()->GetNeutrino().InteractionType() <= 1090)
        hSliceEFracRemovedRES->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
      if (true_nu->Truth()->GetNeutrino().InteractionType() == 1091)
        hSliceEFracRemovedDIS->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
      if (mrcc)
        hMRCCSliceEFracRemoved->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
    }

    float had_charge_before = true_nu->ClusterHadronicCharge(), had_charge_after = true_nu->ClusterMRCCHadronicCharge();
    if (had_charge_before) {
      hSliceHadEFracRemoved->Fill(true_nu->Energy(), (had_charge_before-had_charge_after)/had_charge_before);
      hSliceDeltaHadE->Fill((had_charge_before-had_charge_after)/had_charge_before);
      if (mrcc)
        hMRCCSliceHadEFracRemoved->Fill(true_nu->Energy(), (had_charge_before-had_charge_after)/had_charge_before);
    }

    float lep_charge_before = true_nu->ClusterLeptonicCharge(), lep_charge_after = true_nu->ClusterMRCCLeptonicCharge();
    if (lep_charge_before) {
      hSliceLepEFracRemoved->Fill(true_nu->Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
      hSliceDeltaLepE->Fill((lep_charge_before-lep_charge_after)/lep_charge_before);
      if (mrcc)
        hMRCCSliceLepEFracRemoved->Fill(true_nu->Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
    }

    if (mrcc)
      hHadYMRCC->Fill(true_nu->Truth()->GetNeutrino().Y());
    else
      hHadYNoMRCC->Fill(true_nu->Truth()->GetNeutrino().Y());

    // Particle-level plots
    std::map<unsigned int, MRCCParticle*> particles = true_nu->Particles();
    int muon = true_nu->Muon();
    for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = particles.begin();
         particleIt != particles.end(); ++particleIt) {
      float p_charge_before = particleIt->second->DepositedCharge();
      float p_charge_after = particleIt->second->DepositedMRCCCharge();
      if (!p_charge_before)
        continue;
      if ((int)particleIt->first == muon) {
        hLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
        hAnyMRCCMuonEnergy->Fill(particleIt->second->Energy(), mrcc);
        if (mrcc)
          hMRCCLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
      } else {
        hHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
        if (mrcc)
          hMRCCHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
      }
    }

  }

  return;

}

// --------------------------------------------------------------------------
bool murem::MRCCAna::SliceSelection(const art::Event& evt, const art::Ptr<rb::Cluster>& slice) {

  // Copied from MuonRemove_module

  // rock filtered?
  if (rb::IsFiltered(evt, slice, "rockpresel"))
    return false;

  // find vertex
  art::FindManyP<rb::Vertex> fmElastic({slice}, evt, fElasticModuleLabel);
  if (!fmElastic.isValid())
    return false;

  std::vector<art::Ptr<rb::Vertex> > elastics = fmElastic.at(0);
  if (elastics.empty())
    return false;

  // find prongs
  art::FindManyP<rb::Prong> fmFuzzyProng3D(elastics, evt, fFuzzy3DModuleLabel);
  std::vector<art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
  if (prongs.empty())
    return false;

  // find energy of parent muon
  double muonEmax = -1;
  art::FindManyP<rb::Track> fmt({slice}, evt, fBPFTrackModuleLabel);
  std::vector<art::Ptr<rb::Track> > tracks = fmt.at(0);
  for (unsigned int trackIt = 0; trackIt < tracks.size(); ++trackIt) {
    art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[trackIt]}, evt, fBPFTrackModuleLabel);
    if (fmBPFFitSums.isValid()) {
      std::vector<art::Ptr<rb::FitSum> > bpfFitSums = fmBPFFitSums.at(0);
      if (abs(bpfFitSums[0]->PDG()) == 13) {
        double muonE = bpfFitSums[0]->FourMom().E();
        if (muonE > muonEmax)
          muonEmax = muonE;
      }
    }
  }
  if (muonEmax <= 0)
    return false;

  return true;

}