ShowerQuality_module.cc

Go to the documentation of this file.

// Class:       ShowerQuality
// Module Type: analyzer
// File:        ShowerQuality_module.cc
//
// Generated at Tue Mar 31 07:22:17 2015 by Kazuhiro Terao using artmod
// from cetpkgsupport v1_08_05.

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

#include "ShowerRecoException.h"
#include "larcore/Geometry/Geometry.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardataobj/MCBase/MCShower.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Shower.h"
#include "larreco/MCComp/MCMatchAlg.h"

#include "TH1D.h"
#include "TH2D.h"
#include "TTree.h"

#include <map>

class ShowerQuality : public art::EDAnalyzer {
public:
  explicit ShowerQuality(fhicl::ParameterSet const& p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  ShowerQuality(ShowerQuality const&) = delete;
  ShowerQuality(ShowerQuality&&) = delete;
  ShowerQuality& operator=(ShowerQuality const&) = delete;
  ShowerQuality& operator=(ShowerQuality&&) = delete;

private:
  // Required functions.
  void analyze(art::Event const& e) override;

  void beginJob() override;

  void SetShowerProducer(const std::string name) { fShowerProducer = name; }
  void SetMCShowerProducer(const std::string name) { fMCShowerProducer = name; }
  void SetSimChannelProducer(const std::string name) { fSimChannelProducer = name; }

  void SetMaxEnergyCut(const double energy) { _mc_energy_max = energy; }

  void SetMinEnergyCut(const double energy) { _mc_energy_min = energy; }

  ::btutil::MCMatchAlg fBTAlg;

  double _mc_energy_min;

  double _mc_energy_max;

  std::string fShowerProducer;

  std::string fMCShowerProducer;

  std::string fSimChannelProducer;

  TH1D* hMatchCorrectness;

  TH1D* hVtxDX; 
  TH1D* hVtxDY; 
  TH1D* hVtxDZ; 
  TH1D* hVtxDR; 

  TH1D* hDCosX;       
  TH1D* hDCosY;       
  TH1D* hDCosZ;       
  TH1D* h3DAngleDiff; 

  TH2D* hEnergyCorr; 

  TH1D* hEnergyAssym; 
  TH1D* hEnergyDiff;  

  TH1D* hMatchedClusterEff; 
  TH1D* hMatchedClusterPur; 

  std::map<int, TH1D*> mDEDX;

  TH1D* hBestPlane;

  struct TreeParams_t {

    double reco_x, reco_y, reco_z;
    double reco_dcosx, reco_dcosy, reco_dcosz;
    double reco_energy;
    double reco_dedx;
    double reco_dedx_U;
    double reco_dedx_V;
    double reco_dedx_Y;
    int best_plane_id;

    double mc_x, mc_y, mc_z;
    double mc_dcosx, mc_dcosy, mc_dcosz;
    double mc_energy;
    int mc_pdgid;

    double mc_reco_anglediff;
    double mc_reco_dist;

    double mc_containment;
    double match_correctness;
    double cluster_eff;
    double cluster_pur;

  } fTreeParams;

  TTree* fTree;

  void InitializeAnaTree();

  // Declare member data here.
};

ShowerQuality::ShowerQuality(fhicl::ParameterSet const& p) : EDAnalyzer(p) // ,
// More initializers here.
{

  //fShowerProducer   = "";
  //fMCShowerProducer = "";
  //fSimChannelProducer = "";
  SetShowerProducer(p.get<std::string>("ShowerProducer"));
  SetMCShowerProducer(p.get<std::string>("MCShowerProducer"));
  SetSimChannelProducer(p.get<std::string>("SimChannelProducer"));
  SetMinEnergyCut(p.get<double>("MCShowerEnergyMin"));
  SetMaxEnergyCut(p.get<double>("MCShowerEnergyMax"));

  hMatchCorrectness = nullptr;

  hVtxDX = nullptr;
  hVtxDY = nullptr;
  hVtxDZ = nullptr;
  hVtxDR = nullptr;

  hDCosX = nullptr;
  hDCosY = nullptr;
  hDCosZ = nullptr;
  h3DAngleDiff = nullptr;

  hEnergyCorr = nullptr;
  hEnergyAssym = nullptr;
  hEnergyDiff = nullptr;

  hMatchedClusterPur = nullptr;
  hMatchedClusterEff = nullptr;

  mDEDX.clear();
  hBestPlane = nullptr;

  fTree = nullptr;
}

void ShowerQuality::beginJob()
{

  if (fShowerProducer.empty() || fMCShowerProducer.empty() || fSimChannelProducer.empty()) {
    std::string msg;
    msg += "\033[93m[ERROR]\033[00m <<";
    msg += __FUNCTION__;
    msg += ">> Producer's name not set!";
    std::cout << msg.c_str() << std::endl;
    throw ::showerreco::ShowerRecoException(msg.c_str());
  }

  art::ServiceHandle<geo::Geometry const> geo;
  art::ServiceHandle<art::TFileService const> tfs;

  if (fTree) delete fTree;
  fTree = tfs->make<TTree>("fShowerQualityTree", "");

  //
  // Matching correctness histogram initialization
  //
  if (hMatchCorrectness) delete hMatchCorrectness;
  hMatchCorrectness =
    tfs->make<TH1D>("hMatchCorrectness",
                    "Shower 2D Cluster Matching Correctness; Correctness; Showers",
                    101,
                    -0.005,
                    1.005);

  //
  // 3D Vtx (start point) MC/Reco comparison histogram initialization
  //
  if (hVtxDX) delete hVtxDX;
  if (hVtxDY) delete hVtxDY;
  if (hVtxDZ) delete hVtxDZ;
  if (hVtxDR) delete hVtxDR;

  hVtxDX = tfs->make<TH1D>(
    "hVtxDX", "Reco - MC Start X [cm] Displacement; #DeltaX [cm]; Showers", 200, -100, 100);

  hVtxDY = tfs->make<TH1D>(
    "hVtxDY", "Reco - MC Start Y [cm] Displacement; #DeltaY [cm]; Showers", 200, -100, 100);

  hVtxDZ = tfs->make<TH1D>(
    "hVtxDZ", "Reco - MC Start Z [cm] Displacement; #DeltaZ [cm]; Showers", 200, -100, 100);

  hVtxDR = tfs->make<TH1D>(
    "hVtxDR", "Reco - MC Start 3D Vtx Displacement; #DeltaR [cm]; Showers", 200, -100, 100);

  //
  // 3D Angular MC/Reco comparison histogram initialization
  //
  if (hDCosX) delete hDCosX;
  if (hDCosY) delete hDCosY;
  if (hDCosZ) delete hDCosZ;
  if (h3DAngleDiff) delete h3DAngleDiff;

  hDCosX = tfs->make<TH1D>(
    "hDCosX", "Direction Unit Vector Reco - MC #DeltaX; #DeltaCosX; Showers", 100, -2, 2);

  hDCosY = tfs->make<TH1D>(
    "hDCosY", "Direction Unit Vector Reco - MC #DeltaY; #DeltaCosY; Showers", 100, -2, 2);

  hDCosZ = tfs->make<TH1D>(
    "hDCosZ", "Direction Unit Vector Reco - MC #DeltaZ; #DeltaCosZ; Showers", 100, -2, 2);

  h3DAngleDiff =
    tfs->make<TH1D>("h3DAngleDiff",
                    "3D Opening Angle Between Reco & MC; Opening Angle [degrees]; Showers",
                    181,
                    -0.5,
                    180.5);

  //
  // Energy MC/Reco comparison histogram initialization
  //
  if (hEnergyCorr) delete hEnergyCorr;
  if (hEnergyAssym) delete hEnergyAssym;
  if (hEnergyDiff) delete hEnergyDiff;

  hEnergyCorr =
    tfs->make<TH2D>("hEnergyCorr",
                    "Reco (x) vs. MC (y) Energy Comparison; Reco Energy [MeV]; MC Energy [MeV]",
                    200,
                    0,
                    1000,
                    200,
                    0,
                    1000);

  hEnergyAssym = tfs->make<TH1D>("hEnergyAssym",
                                 "MC - Reco Energy Fractional Difference; Assymetry; Showers",
                                 201,
                                 -1.005,
                                 1.005);

  hEnergyDiff = tfs->make<TH1D>(
    "hEnergyDiff", "MC - Reco Energy Difference; Energy Difference [MeV]; Showers", 200, 0, 1000);

  //
  // Shower cluster purity & efficiency histograms initialization
  //
  if (hMatchedClusterEff) delete hMatchedClusterEff;
  if (hMatchedClusterPur) delete hMatchedClusterPur;

  hMatchedClusterEff =
    tfs->make<TH1D>("hMatchedClusterEff_PlaneCombo",
                    "Matched Shower Cluster's Charge Efficiency; Efficiency; Clusters",
                    101,
                    -0.005,
                    1.005);

  hMatchedClusterPur = tfs->make<TH1D>("hMatchedClusterPur_PlaneCombo",
                                       "Matched Shower Cluster's Charge Purity; Purity; Clusters",
                                       101,
                                       -0.005,
                                       1.005);

  //
  // Best plane ID histogram initialization
  //
  hBestPlane = tfs->make<TH1D>("hBestPlane",
                               "Best Plane (for energy & dE/dx estimate); Plane ID; Showers",
                               geo->Nplanes(),
                               -0.5,
                               geo->Nplanes() - 0.5);

  InitializeAnaTree();
}

void ShowerQuality::analyze(art::Event const& e)
{
  // Retrieve mcshower data product
  auto mcsHandle = e.getValidHandle<std::vector<sim::MCShower>>(fMCShowerProducer);
  auto resHandle = e.getValidHandle<std::vector<recob::Shower>>(fShowerProducer);
  auto schHandle = e.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelProducer);
  const std::vector<sim::MCShower>& ev_mcs(*mcsHandle);
  const std::vector<recob::Shower>& ev_shower(*resHandle);
  const std::vector<sim::SimChannel>& ev_simch(*schHandle);

  if (ev_shower.empty()) return;

  // Get the whole clusters + associated clusters
  art::Handle<std::vector<recob::Cluster>> clsHandle;
  art::FindManyP<recob::Cluster> cluster_m(resHandle, e, fShowerProducer);
  const std::vector<recob::Cluster>& ev_cluster = cluster_m.at(0).front().parentAs<std::vector>();

  // Make clusters in terms of hit vector to feed into BT algorithm
  art::FindManyP<recob::Hit> hit_m(clsHandle, e, clsHandle.provenance()->moduleLabel());
  std::vector<std::vector<art::Ptr<recob::Hit>>> ev_cluster_hit;
  ev_cluster_hit.reserve(clsHandle->size());
  std::map<art::Ptr<recob::Cluster>, size_t> cluster_ptr_map;
  for (size_t i = 0; i < ev_cluster.size(); ++i) {
    const art::Ptr<recob::Cluster> cluster_ptr(clsHandle, i);
    cluster_ptr_map[cluster_ptr] = ev_cluster_hit.size();
    ev_cluster_hit.push_back(hit_m.at(i));
  }

  // Create ass_cluster_v index vector
  std::vector<std::vector<unsigned int>> ass_cluster_v;
  ass_cluster_v.reserve(ev_shower.size());
  for (size_t shower_index = 0; shower_index < ev_shower.size(); ++shower_index) {
    ass_cluster_v.push_back(std::vector<unsigned int>());
    for (auto const& p : cluster_m.at(shower_index))
      ass_cluster_v.back().push_back(cluster_ptr_map[p]);
  }

  // Create G4 track ID vector for which we are interested in
  std::vector<std::vector<unsigned int>> g4_trackid_v;
  std::vector<unsigned int> mc_index_v;
  g4_trackid_v.reserve(ev_mcs.size());
  for (size_t mc_index = 0; mc_index < ev_mcs.size(); ++mc_index) {
    auto const& mcs = ev_mcs[mc_index];
    double energy = mcs.DetProfile().E();
    std::vector<unsigned int> id_v;
    id_v.reserve(mcs.DaughterTrackID().size());
    if (_mc_energy_min < energy && energy < _mc_energy_max) {
      for (auto const& id : mcs.DaughterTrackID()) {
        if (id == mcs.TrackID()) continue;
        id_v.push_back(id);
      }
      id_v.push_back(mcs.TrackID());
      g4_trackid_v.push_back(id_v);
      mc_index_v.push_back(mc_index);
    }
  }

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
  if (!fBTAlg.BuildMap(clockData, g4_trackid_v, ev_simch, ev_cluster_hit)) {
    std::cerr << "\033[93m[ERROR]\033[00m <<ShowerQuality::analyze>> Failed to "
                 "build back-tracking map for MC..."
              << std::endl;
    return;
  }

  // Find the best-representative reco-ed Shower given an MCShower
  std::vector<std::vector<double>> shower_mcq_vv(ev_shower.size(),
                                                 std::vector<double>(mc_index_v.size(), 0));

  for (size_t shower_index = 0; shower_index < ass_cluster_v.size(); ++shower_index) {

    auto const& ass_cluster = ass_cluster_v[shower_index];

    std::vector<::btutil::WireRange_t> w_v;

    for (auto const& cluster_index : ass_cluster) {

      auto const& ass_hit = ev_cluster_hit[cluster_index];

      w_v.reserve(ass_hit.size() + w_v.size());

      for (auto const& hit_ptr : ass_hit) {
        w_v.emplace_back(hit_ptr->Channel(), hit_ptr->StartTick(), hit_ptr->EndTick());
      }
    }

    auto mcq_v = fBTAlg.BTAlg().MCQ(clockData, w_v);

    auto& shower_mcq_v = shower_mcq_vv[shower_index];

    for (size_t mcs_index = 0; mcs_index < (mcq_v.size() - 1); ++mcs_index) {

      shower_mcq_v[mcs_index] = mcq_v[mcs_index];
    }
  }

  // Loop over MCShower and inspect corresponding shower quality
  for (size_t mcs_index = 0; mcs_index < mc_index_v.size(); ++mcs_index) {

    auto const& mc_shower = ev_mcs[mc_index_v[mcs_index]];

    // Search for the best representative shower
    size_t best_shower_index = shower_mcq_vv.size();
    double max_mcq = 0;
    for (size_t shower_index = 0; shower_index < shower_mcq_vv.size(); ++shower_index) {

      if (shower_mcq_vv[shower_index][mcs_index] > max_mcq) best_shower_index = shower_index;
    }

    if (best_shower_index == shower_mcq_vv.size()) {
      std::string msg;
      std::cerr << "\033[93m[ERROR]\033[00m <<ShowerQuality::analyze>> "
                << "Failed to find a corresponding shower for MCShower " << mc_index_v[mcs_index]
                << std::endl;
      continue;
    }

    auto const& reco_shower = ev_shower[best_shower_index];

    auto res = fBTAlg.ShowerCorrectness(ass_cluster_v[best_shower_index]);

    fTreeParams.match_correctness = res.second;

    if (fTreeParams.match_correctness < 0) {
      std::cerr << "\033[93m[ERROR]\033[00m <<ShowerQuality::analyze>> "
                << "Failed to find a corresponding MCShower for shower " << best_shower_index
                << std::endl;
      continue;
    }

    // MC Info
    fTreeParams.mc_x = mc_shower.DetProfile().X();
    fTreeParams.mc_y = mc_shower.DetProfile().Y();
    fTreeParams.mc_z = mc_shower.DetProfile().Z();

    fTreeParams.mc_energy = mc_shower.DetProfile().E();
    fTreeParams.mc_pdgid = mc_shower.PdgCode();
    fTreeParams.mc_containment = mc_shower.DetProfile().E() / mc_shower.Start().E();

    //fTreeParams.mc_dcosx = mc_shower.DetProfile().Px() / fTreeParams.mc_energy;
    //fTreeParams.mc_dcosy = mc_shower.DetProfile().Py() / fTreeParams.mc_energy;
    //fTreeParams.mc_dcosz = mc_shower.DetProfile().Pz() / fTreeParams.mc_energy;
    fTreeParams.mc_dcosx = mc_shower.Start().Px() / mc_shower.Start().E();
    fTreeParams.mc_dcosy = mc_shower.Start().Py() / mc_shower.Start().E();
    fTreeParams.mc_dcosz = mc_shower.Start().Pz() / mc_shower.Start().E();

    // Reco vtx
    fTreeParams.reco_x = reco_shower.ShowerStart()[0];
    fTreeParams.reco_y = reco_shower.ShowerStart()[1];
    fTreeParams.reco_z = reco_shower.ShowerStart()[2];

    // Reco angle
    fTreeParams.reco_dcosx = reco_shower.Direction()[0];
    fTreeParams.reco_dcosy = reco_shower.Direction()[1];
    fTreeParams.reco_dcosz = reco_shower.Direction()[2];

    // Reco - MC angle diff
    fTreeParams.mc_reco_anglediff = acos(fTreeParams.reco_dcosx * fTreeParams.mc_dcosx +
                                         fTreeParams.reco_dcosy * fTreeParams.mc_dcosy +
                                         fTreeParams.reco_dcosz * fTreeParams.mc_dcosz) /
                                    3.14159265359 * 180.;
    // Reco - MC vtx distance
    fTreeParams.mc_reco_dist = sqrt(pow(fTreeParams.reco_x - fTreeParams.mc_x, 2) +
                                    pow(fTreeParams.reco_y - fTreeParams.mc_y, 2) +
                                    pow(fTreeParams.reco_z - fTreeParams.mc_z, 2));

    // Reco cluster efficiency & purity
    fTreeParams.cluster_eff = 1.;
    fTreeParams.cluster_pur = 1.;
    for (auto const& cluster_index : ass_cluster_v[best_shower_index]) {
      auto ep = fBTAlg.ClusterEP(cluster_index, mcs_index);
      if (ep.first == 0 && ep.second == 0) continue;
      fTreeParams.cluster_eff *= ep.first;
      fTreeParams.cluster_pur *= ep.second;
    }

    // Reco energy & dedx info
    fTreeParams.best_plane_id = reco_shower.best_plane();

    /*
      int best_plane_index = -1;

      for(size_t i=0; i < ass_cluster_v[best_shower_index].size(); ++i) {

        size_t cluster_index = ass_cluster_v[best_shower_index][i];
        //std::cout<<best_plane_index<<" : "<<ev_cluster->at(cluster_index).View()<<std::endl;
        if( ev_cluster->at(cluster_index).View() == reco_shower.best_plane() ) {
          best_plane_index = i;
          break;
        }
      }

      if(best_plane_index < 0) {
      throw ::showerreco::ShowerRecoException(Form("Failed to identify the best plane for shower %zu",
      best_shower_index)
      );
      }
    */

    fTreeParams.reco_energy = reco_shower.Energy().at(reco_shower.best_plane());
    //fTreeParams.reco_dedx_U   = reco_shower.dEdx().at(0);
    //fTreeParams.reco_dedx_V   = reco_shower.dEdx().at(1);
    //fTreeParams.reco_dedx_Y   = reco_shower.dEdx().at(2);
    fTreeParams.reco_dedx = reco_shower.dEdx().at(reco_shower.best_plane());

    //
    // Fill histograms
    //
    hMatchCorrectness->Fill(fTreeParams.match_correctness);

    hVtxDX->Fill(fTreeParams.reco_x - fTreeParams.mc_x);
    hVtxDY->Fill(fTreeParams.reco_y - fTreeParams.mc_y);
    hVtxDZ->Fill(fTreeParams.reco_z - fTreeParams.mc_z);
    hVtxDR->Fill(fTreeParams.mc_reco_dist);

    // Angular info
    hDCosX->Fill(fTreeParams.reco_dcosx - fTreeParams.mc_dcosx);
    hDCosY->Fill(fTreeParams.reco_dcosy - fTreeParams.mc_dcosy);
    hDCosZ->Fill(fTreeParams.reco_dcosz - fTreeParams.mc_dcosz);
    h3DAngleDiff->Fill(fTreeParams.mc_reco_anglediff);

    hEnergyCorr->Fill(fTreeParams.reco_energy, fTreeParams.mc_energy);

    hEnergyAssym->Fill((fTreeParams.reco_energy - fTreeParams.mc_energy) /
                       (fTreeParams.reco_energy + fTreeParams.mc_energy) * 2.);

    hEnergyDiff->Fill(fTreeParams.mc_energy - fTreeParams.reco_energy);

    if (mDEDX.find(fTreeParams.mc_pdgid) == mDEDX.end())

      mDEDX.insert(std::make_pair(
        fTreeParams.mc_pdgid,
        new TH1D(Form("hdEdx_PDG_%d", fTreeParams.mc_pdgid),
                 Form("Reco dE/dx for PDG = %d; dE/dx [MeV/cm]; Showers", fTreeParams.mc_pdgid),
                 100,
                 0,
                 50)));

    hMatchedClusterEff->Fill(fTreeParams.cluster_eff);
    hMatchedClusterPur->Fill(fTreeParams.cluster_pur);

    mDEDX[fTreeParams.mc_pdgid]->Fill(fTreeParams.reco_dedx);

    hBestPlane->Fill(fTreeParams.best_plane_id);

    // Fill Tree
    fTree->Fill();
  }
}

void ShowerQuality::InitializeAnaTree()
{

  fTree->Branch("reco_x", &fTreeParams.reco_x, "reco_x/D");
  fTree->Branch("reco_y", &fTreeParams.reco_y, "reco_y/D");
  fTree->Branch("reco_z", &fTreeParams.reco_z, "reco_z/D");
  fTree->Branch("reco_dcosx", &fTreeParams.reco_dcosx, "reco_dcosx/D");
  fTree->Branch("reco_dcosy", &fTreeParams.reco_dcosy, "reco_dcosy/D");
  fTree->Branch("reco_dcosz", &fTreeParams.reco_dcosz, "reco_dcosz/D");
  fTree->Branch("reco_energy", &fTreeParams.reco_energy, "reco_energy/D");

  fTree->Branch("best_plane_id", &fTreeParams.best_plane_id, "best_plane_id/i");

  fTree->Branch("mc_x", &fTreeParams.mc_x, "mc_x/D");
  fTree->Branch("mc_y", &fTreeParams.mc_y, "mc_y/D");
  fTree->Branch("mc_z", &fTreeParams.mc_z, "mc_z/D");
  fTree->Branch("mc_dcosx", &fTreeParams.mc_dcosx, "mc_dcosx/D");
  fTree->Branch("mc_dcosy", &fTreeParams.mc_dcosy, "mc_dcosy/D");
  fTree->Branch("mc_dcosz", &fTreeParams.mc_dcosz, "mc_dcosz/D");
  fTree->Branch("mc_energy", &fTreeParams.mc_energy, "mc_energy/D");

  fTree->Branch("reco_dedx", &fTreeParams.reco_dedx, "reco_dedx_/D");
  fTree->Branch("reco_dedx_U", &fTreeParams.reco_dedx_U, "reco_dedx_U/D");
  fTree->Branch("reco_dedx_V", &fTreeParams.reco_dedx_V, "reco_dedx_V/D");
  fTree->Branch("reco_dedx_Y", &fTreeParams.reco_dedx_Y, "reco_dedx_Y/D");
  fTree->Branch("mc_pdgid", &fTreeParams.mc_pdgid, "mc_pdgid/i");

  fTree->Branch("mc_reco_anglediff", &fTreeParams.mc_reco_anglediff, "mc_reco_anglediff/D");
  fTree->Branch("mc_reco_dist", &fTreeParams.mc_reco_dist, "mc_reco_dist/D");

  fTree->Branch("mc_containment", &fTreeParams.mc_containment, "mc_containment/D");

  fTree->Branch("match_correctness", &fTreeParams.match_correctness, "match_correctness/D");
  fTree->Branch("cluster_eff", &fTreeParams.cluster_eff, "cluster_eff/D");
  fTree->Branch("cluster_pur", &fTreeParams.cluster_pur, "cluster_pur/D");
}

DEFINE_ART_MODULE(ShowerQuality)