TCShowerAnalysis_module.cc

Go to the documentation of this file.

// -------------------------------------------------
// tcshower analysis tree
//
// Author: Rory Fitzpatrick (roryfitz@umich.edu)
// Created: 7/16/18
// -------------------------------------------------

// Framework includes
#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 "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Shower.h"
#include "lardataobj/Simulation/SimChannel.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "larsim/MCCheater/BackTrackerService.h"
#include "larsim/MCCheater/ParticleInventoryService.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "TTree.h"

constexpr int kMaxShowers = 1000; // maximum number of showers

namespace shower {

  class TCShowerAnalysis : public art::EDAnalyzer {
  public:
    explicit TCShowerAnalysis(fhicl::ParameterSet const& pset);

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

    void reset();

    void truthMatcher(detinfo::DetectorClocksData const& clockData,
                      std::vector<art::Ptr<recob::Hit>> all_hits,
                      std::vector<art::Ptr<recob::Hit>> shower_hits,
                      const simb::MCParticle*& MCparticle,
                      double& Efrac,
                      double& Ecomplet);

    TTree* fTree;
    int run;
    int subrun;
    int event;
    int nuPDG_truth;
    int ccnc_truth;
    int mode_truth;
    int nshws;                   // number of showers
    int shwid[kMaxShowers];      // recob::Shower::ID()
    float shwdcosx[kMaxShowers]; // shower direction cosin
    float shwdcosy[kMaxShowers];
    float shwdcosz[kMaxShowers];
    float shwstartx[kMaxShowers]; // shower start position (cm)
    float shwstarty[kMaxShowers];
    float shwstartz[kMaxShowers];
    double shwdedx[kMaxShowers][2]; // shower dE/dx of the initial track
                                    // measured on the 3 plane (MeV/cm)
    int shwbestplane[kMaxShowers];  // recommended plane for energy and dE/dx
                                    // information

    int highestHitsPDG;
    double highestHitsFrac;

    std::string fHitModuleLabel;
    std::string fShowerModuleLabel;
    std::string fGenieGenModuleLabel;
    std::string fDigitModuleLabel;

    calo::CalorimetryAlg fCalorimetryAlg;

  }; // class TCShowerAnalysis

} // shower

// -------------------------------------------------

shower::TCShowerAnalysis::TCShowerAnalysis(fhicl::ParameterSet const& pset)
  : EDAnalyzer(pset)
  , fHitModuleLabel(pset.get<std::string>("HitModuleLabel", "trajcluster"))
  , fShowerModuleLabel(pset.get<std::string>("ShowerModuleLabel", "tcshower"))
  , fGenieGenModuleLabel(pset.get<std::string>("GenieGenModuleLabel", "generator"))
  , fDigitModuleLabel(pset.get<std::string>("DigitModuleLabel", "largeant"))
  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
{} // TCShowerTemplateMaker

// -------------------------------------------------

void shower::TCShowerAnalysis::beginJob()
{

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

  fTree = tfs->make<TTree>("tcshowerana", "tcshowerana");
  fTree->Branch("run", &run, "run/I");
  fTree->Branch("subrun", &subrun, "subrun/I");
  fTree->Branch("event", &event, "event/I");

  fTree->Branch("nuPDG_truth", &nuPDG_truth, "nuPDG_truth/I");
  fTree->Branch("ccnc_truth", &ccnc_truth, "ccnc_truth/I");
  fTree->Branch("mode_truth", &mode_truth, "mode_truth/I");

  fTree->Branch("nshws", &nshws, "nshws/I");
  fTree->Branch("shwid", shwid, "shwid[nshws]/I");
  fTree->Branch("shwdcosx", shwdcosx, "shwdcosx[nshws]/F");
  fTree->Branch("shwdcosy", shwdcosy, "shwdcosy[nshws]/F");
  fTree->Branch("shwdcosz", shwdcosz, "shwdcosz[nshws]/F");
  fTree->Branch("shwstartx", shwstartx, "shwstartx[nshws]/F");
  fTree->Branch("shwstarty", shwstarty, "shwstarty[nshws]/F");
  fTree->Branch("shwstartz", shwstartz, "shwstartz[nshws]/F");
  fTree->Branch("shwdedx", shwdedx, "shwdedx[nshws][2]/D");
  fTree->Branch("shwbestplane", shwbestplane, "shwbestplane[nshws]/I");

  fTree->Branch("highestHitsPDG", &highestHitsPDG, "highestHitsPDG/I");
  fTree->Branch("highestHitsFrac", &highestHitsFrac, "highestHitsFrac/D");

} // beginJob

// -------------------------------------------------

void shower::TCShowerAnalysis::analyze(const art::Event& evt)
{

  run = evt.run();
  subrun = evt.subRun();
  event = evt.id().event();

  art::Handle<std::vector<recob::Hit>> hitListHandle;
  std::vector<art::Ptr<recob::Hit>> hitlist;
  if (evt.getByLabel(fHitModuleLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle);

  art::Handle<std::vector<sim::SimChannel>> scListHandle;
  std::vector<art::Ptr<sim::SimChannel>> simchanlist;
  if (evt.getByLabel(fDigitModuleLabel, scListHandle))
    art::fill_ptr_vector(simchanlist, scListHandle);

  art::Handle<std::vector<recob::Shower>> showerListHandle;
  std::vector<art::Ptr<recob::Shower>> showerlist;
  if (evt.getByLabel(fShowerModuleLabel, showerListHandle))
    art::fill_ptr_vector(showerlist, showerListHandle);

  art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
  std::vector<art::Ptr<simb::MCTruth>> mclist;
  if (evt.getByLabel(fGenieGenModuleLabel, mctruthListHandle))
    art::fill_ptr_vector(mclist, mctruthListHandle);

  art::FindManyP<recob::Hit> shwfm(showerListHandle, evt, fShowerModuleLabel);

  //shower information
  if (showerListHandle.isValid()) {
    nshws = showerlist.size();

    for (int i = 0; i < std::min(int(showerlist.size()), kMaxShowers); ++i) {
      shwid[i] = showerlist[i]->ID();
      shwdcosx[i] = showerlist[i]->Direction().X();
      shwdcosy[i] = showerlist[i]->Direction().Y();
      shwdcosz[i] = showerlist[i]->Direction().Z();
      shwstartx[i] = showerlist[i]->ShowerStart().X();
      shwstarty[i] = showerlist[i]->ShowerStart().Y();
      shwstartz[i] = showerlist[i]->ShowerStart().Z();
      for (size_t j = 0; j < (showerlist[i]->dEdx()).size(); ++j) {
        shwdedx[i][j] = showerlist[i]->dEdx()[j];
      }
      shwbestplane[i] = showerlist[i]->best_plane();
    }

  } // shower info

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);

  if (mclist.size()) {
    art::Ptr<simb::MCTruth> mctruth = mclist[0];
    if (mctruth->NeutrinoSet()) {

      nuPDG_truth = mctruth->GetNeutrino().Nu().PdgCode();
      ccnc_truth = mctruth->GetNeutrino().CCNC();
      mode_truth = mctruth->GetNeutrino().Mode();

      if (showerlist.size()) { // only looks at the first shower since this is
                               // for tcshower
        std::vector<art::Ptr<recob::Hit>> showerhits = shwfm.at(0);
        // get shower truth info
        const simb::MCParticle* particle;
        double tmpEfrac = 0.0;
        double tmpEcomplet = 0;
        truthMatcher(clockData, hitlist, showerhits, particle, tmpEfrac, tmpEcomplet);
        if (particle) {
          std::cout << "shower pdg: " << particle->PdgCode() << " efrac " << tmpEfrac << std::endl;

          highestHitsPDG = particle->PdgCode();
          highestHitsFrac = tmpEfrac;
        }
      }
    }
  }

  fTree->Fill();

} // analyze

// -------------------------------------------------

void shower::TCShowerAnalysis::reset()
{

  run = -99999;
  subrun = -99999;
  event = -99999;

  nuPDG_truth = -99999;
  ccnc_truth = -99999;
  mode_truth = -99999;

  nshws = 0;
  for (int i = 0; i < kMaxShowers; ++i) {
    shwid[i] = -99999;
    shwdcosx[i] = -99999;
    shwdcosy[i] = -99999;
    shwdcosz[i] = -99999;
    shwstartx[i] = -99999;
    shwstarty[i] = -99999;
    shwstartz[i] = -99999;
    for (int j = 0; j < 2; ++j) {
      shwdedx[i][j] = -99999;
    }
    shwbestplane[i] = -99999;
  }

  highestHitsPDG = -99999;
  highestHitsFrac = -99999;

} // reset

// -------------------------------------------------

void shower::TCShowerAnalysis::truthMatcher(detinfo::DetectorClocksData const& clockData,
                                            std::vector<art::Ptr<recob::Hit>> all_hits,
                                            std::vector<art::Ptr<recob::Hit>> shower_hits,
                                            const simb::MCParticle*& MCparticle,
                                            double& Efrac,
                                            double& Ecomplet)
{

  MCparticle = 0;
  Efrac = 1.0;
  Ecomplet = 0;

  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
  std::map<int, double> trkID_E;
  for (size_t j = 0; j < shower_hits.size(); ++j) {
    art::Ptr<recob::Hit> hit = shower_hits[j];
    // For know let's use collection plane to look at the shower reconstruction
    if (hit->View() != 1) continue;
    std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
    for (size_t k = 0; k < TrackIDs.size(); k++) {
      if (trkID_E.find(std::abs(TrackIDs[k].trackID)) == trkID_E.end())
        trkID_E[std::abs(TrackIDs[k].trackID)] = 0;
      trkID_E[std::abs(TrackIDs[k].trackID)] += TrackIDs[k].energy;
    }
  }
  double max_E = -999.0;
  double total_E = 0.0;
  int TrackID = -999;
  double partial_E = 0.0;

  if (!trkID_E.size()) return; // Ghost shower???
  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end(); ++ii) {
    total_E += ii->second;
    if ((ii->second) > max_E) {
      partial_E = ii->second;
      max_E = ii->second;
      TrackID = ii->first;
    }
  }
  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);

  Efrac = partial_E / total_E;

  // completeness
  double totenergy = 0;
  for (size_t k = 0; k < all_hits.size(); ++k) {
    art::Ptr<recob::Hit> hit = all_hits[k];
    std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
    for (size_t l = 0; l < TrackIDs.size(); ++l) {
      if (std::abs(TrackIDs[l].trackID) == TrackID) { totenergy += TrackIDs[l].energy; }
    }
  }
  Ecomplet = partial_E / totenergy;

} // truthMatcher

DEFINE_ART_MODULE(shower::TCShowerAnalysis)