CosmicEff_module.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////////////////
/// \brief Calculate efficiency for cosmic rays to identify bad channels
/// \author Christopher Backhouse - bckhouse@caltech.edu
/////////////////////////////////////////////////////////////////////////////////////////

#include <vector>

// ROOT includes
#include "TH1D.h"
#include "TH2D.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes
#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/OfflineChan.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/HitMap.h"

namespace chaninfo
{
  /// \brief Calculate efficiency for cosmic rays to identify bad channels
  class CosmicEff: public art::EDAnalyzer
  {
  public:
    explicit CosmicEff(const fhicl::ParameterSet& pset);
    virtual ~CosmicEff();

    virtual void analyze(const art::Event& evt);
    virtual void reconfigure(const fhicl::ParameterSet& pset);

    virtual void beginRun(const art::Run& run);
    virtual void endJob();

  protected:
    std::string fTrackLabel;

    TH2F *fHits, *fMiss;
    TH2F *fHitsXYLoose, *fMissXYLoose;
    TH2F *fHitsXYTight, *fMissXYTight;
    TH2F *fHitsZLoose, *fMissZLoose;
    TH2F *fHitsZTight, *fMissZTight;
  };

  //......................................................................
  //  HitEfficiency::HitEfficiency(fhicl::ParameterSet const & p)
  //    : EDAnalyzer(p)
  CosmicEff::CosmicEff(const fhicl::ParameterSet& pset)
    : EDAnalyzer(pset), fHits(0) // This one indicates all the others are uninitialized too
  {
    reconfigure(pset);
  }

  //......................................................................
  CosmicEff::~CosmicEff()
  {
  }

  //......................................................................
  void CosmicEff::reconfigure(const fhicl::ParameterSet& pset)
  {
    fTrackLabel = pset.get<std::string>("TrackLabel");
  }

  //......................................................................
  void CosmicEff::beginRun(const art::Run& /*run*/)
  {
    if(fHits) return; // Only make the plots once

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

    // Make histograms big enough to hold all detector cells
    const unsigned int maxCell = std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());
    fHits = tfs->make<TH2F>("hits", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    fMiss = tfs->make<TH2F>("miss", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);

    fHitsXYLoose = tfs->make<TH2F>("hitsXYLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    fMissXYLoose = tfs->make<TH2F>("missXYLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);

    fHitsXYTight = tfs->make<TH2F>("hitsXYTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    fMissXYTight = tfs->make<TH2F>("missXYTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);

    fHitsZLoose = tfs->make<TH2F>("hitsZLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    fMissZLoose = tfs->make<TH2F>("missZLoose", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);

    fHitsZTight = tfs->make<TH2F>("hitsZTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    fMissZTight = tfs->make<TH2F>("missZTight", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
  }

  //......................................................................
  void CosmicEff::analyze(const art::Event& evt)
  {
    art::ServiceHandle<geo::Geometry> geom;

    art::Handle< std::vector<rb::Track> > trks;
    evt.getByLabel(fTrackLabel, trks);
    if(trks->empty()) return;

    const unsigned int trkMax = trks->size();

    for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
      art::Ptr<rb::Track> track(trks, trkIdx);

      rb::HitMap hmap(track.get());

      // Which cells does this track's trajectory pass through?
      std::vector<geo::OfflineChan> xhits, yhits;
      geom->CountCellsOnLineFast(track->Start(), track->Stop(), xhits, yhits);

      // If there's enough hits in this view, then count how many tracks were
      // hit when they should be and how many were missed.
      if(track->NXCell() > 2){
        for(unsigned int n = 0; n < xhits.size(); ++n){
          const unsigned int plane = xhits[n].Plane();
          const unsigned int cell = xhits[n].Cell();
          if(hmap.CellExists(plane, cell)){
            fHits->Fill(plane, cell);
          }
          else{
            fMiss->Fill(plane, cell);
          }
        }
      }
      if(track->NYCell() > 2){
        for(unsigned int n = 0; n < yhits.size(); ++n){
          const unsigned int plane = yhits[n].Plane();
          const unsigned int cell = yhits[n].Cell();
          if(hmap.CellExists(plane, cell)){
            fHits->Fill(plane, cell);
          }
          else{
            fMiss->Fill(plane, cell);
          }
        }
      }

      // Min/Max cell number of hits on this plane
      std::map<int, int> minsXY, maxsXY;
      // Channels that have a hit below/above them
      std::set<geo::OfflineChan> belows, aboves;
      // Min/max plane number by view and cell index
      std::map<int, int> minsZ[2], maxsZ[2];
      // Channels that have a hit upstream/downstream of them
      std::set<geo::OfflineChan> lefts, rights;

      for(unsigned int chitIdx = 0; chitIdx < track->NCell(); ++chitIdx){
        const art::Ptr<rb::CellHit>& chit = track->Cell(chitIdx);
        const int plane = chit->Plane();
        const int cell = chit->Cell();
        const geo::View_t view = chit->View();

        // Have to be more careful on min because value defaults to zero
        if(minsXY.find(plane) == minsXY.end() || cell < minsXY[plane])
          minsXY[plane] = cell;
        if(cell > maxsXY[plane]) maxsXY[plane] = cell;

        if(minsZ[view].find(cell) == minsZ[view].end() || plane < minsZ[view][cell])
          minsZ[view][cell] = plane;
        if(plane > maxsZ[view][cell]) maxsZ[view][cell] = plane;

        belows.insert(geo::OfflineChan(plane, cell+1));
        aboves.insert(geo::OfflineChan(plane, cell-1));

        const int nextplane = geom->NextPlaneInView(plane, +1);
        lefts.insert(geo::OfflineChan(nextplane, cell));
        const int prevplane = geom->NextPlaneInView(plane, -1);
        rights.insert(geo::OfflineChan(prevplane, cell));
      } // end for chitIdx

      for(std::map<int, int>::iterator it = minsXY.begin(); it != minsXY.end(); ++it){
        const int plane = it->first;
        const int from = it->second;
        const int to = maxsXY[plane];

        // Expect a hit if there is a hit above and below somewhere on this plane
        for(int cell = from+1; cell <= to-1; ++cell){
          if(hmap.CellExists(plane, cell))
            fHitsXYLoose->Fill(plane, cell);
          else
            fMissXYLoose->Fill(plane, cell);
        }
      } // end for it

      for(int view = 0; view < 2; ++view){
        for(std::map<int, int>::iterator it = minsZ[view].begin(); it != minsZ[view].end(); ++it){
          const int cell = it->first;
          const int from = it->second;
          const int to = maxsZ[view][cell];

          if(to == from) continue;

          // Except a hit if there's one on an earlier and later plane
          int plane = geom->NextPlaneInView(from, +1);
          while(plane != to){
            if(hmap.CellExists(plane, cell))
              fHitsZLoose->Fill(plane, cell);
            else
              fMissZLoose->Fill(plane, cell);

            plane = geom->NextPlaneInView(plane, +1);
          }
        } // end for it
      } // end for view

      for(std::set<geo::OfflineChan>::iterator it = belows.begin(); it != belows.end(); ++it){
        const int plane = it->Plane();
        const int cell = it->Cell();
        if(!aboves.count(geo::OfflineChan(plane, cell))) continue;

        // Expect a hit if both our neighbours are hit
        if(hmap.CellExists(plane, cell))
          fHitsXYTight->Fill(plane, cell);
        else
          fMissXYTight->Fill(plane, cell);
      }

      for(std::set<geo::OfflineChan>::iterator it = lefts.begin(); it != lefts.end(); ++it){
        const int plane = it->Plane();
        const int cell = it->Cell();
        if(!rights.count(geo::OfflineChan(plane, cell))) continue;

        // Expect a hit if both our neighbouring planes are hit here
        if(hmap.CellExists(plane, cell))
          fHitsZTight->Fill(plane, cell);
        else
          fMissZTight->Fill(plane, cell);
      }
    } // end for trkIdx
  }

  //......................................................................
  void CosmicEff::endJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<chaninfo::BadChanList> bcl;

    // Store what the occupancy cuts think about the channels
    const unsigned int maxCell = std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());
    TH2* hKnown = tfs->make<TH2F>("known", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    //    TH2* hOccs = tfs->make<TH2F>("occs", ";Plane;Cell", geom->NPlanes(), 0, geom->NPlanes(), maxCell, 0, maxCell);
    for(unsigned int p = 0; p < geom->NPlanes(); ++p){
      for(unsigned int c = 0; c < geom->Plane(p)->Ncells(); ++c){
        if(bcl->IsBad(p, c)){
          hKnown->Fill(p, c);
        }
        //        hOccs->Fill(p, c, bcl->GetOccupancy(p, c));
      }
    }

    // Make sums, ratios, and projections of what we already have:

    TH2F* tot = tfs->make<TH2F>(*fHits);
    tot->SetName("tot");
    tot->Add(fMiss);

    TH2F* totXYLoose = tfs->make<TH2F>(*fHitsXYLoose);
    totXYLoose->SetName("totXYLoose");
    totXYLoose->Add(fMissXYLoose);

    TH2F* totXYTight = tfs->make<TH2F>(*fHitsXYTight);
    totXYTight->SetName("totXYTight");
    totXYTight->Add(fMissXYTight);

    TH2F* totZLoose = tfs->make<TH2F>(*fHitsZLoose);
    totZLoose->SetName("totZLoose");
    totZLoose->Add(fMissZLoose);

    TH2F* totZTight = tfs->make<TH2F>(*fHitsZTight);
    totZTight->SetName("totZTight");
    totZTight->Add(fMissZTight);

    TH2F* frac = tfs->make<TH2F>(*fMiss);
    frac->SetName("frac");
    frac->Divide(tot);

    TH2F* fracXYLoose = tfs->make<TH2F>(*fMissXYLoose);
    fracXYLoose->SetName("fracXYLoose");
    fracXYLoose->Divide(totXYLoose);

    TH2F* fracXYTight = tfs->make<TH2F>(*fMissXYTight);
    fracXYTight->SetName("fracXYTight");
    fracXYTight->Divide(totXYTight);

    TH2F* fracZLoose = tfs->make<TH2F>(*fMissZLoose);
    fracZLoose->SetName("fracZLoose");
    fracZLoose->Divide(totZLoose);

    TH2F* fracZTight = tfs->make<TH2F>(*fMissZTight);
    fracZTight->SetName("fracZTight");
    fracZTight->Divide(totZTight);

    TH1F* proj = tfs->make<TH1F>("proj", "", 101, 0, 1.01);
    TH1F* projXYLoose = tfs->make<TH1F>("projXYLoose", "", 101, 0, 1.01);
    TH1F* projXYTight = tfs->make<TH1F>("projXYTight", "", 101, 0, 1.01);
    TH1F* projZLoose = tfs->make<TH1F>("projZLoose", "", 101, 0, 1.01);
    TH1F* projZTight = tfs->make<TH1F>("projZTight", "", 101, 0, 1.01);

    // Projections don't include channels already known to be bad.

    const int X = frac->GetNbinsX()+1;
    const int Y = frac->GetNbinsY()+1;
    for(int x = 1; x < X; ++x){
      for(int y = 1; y < Y; ++y){
        const bool isBad = hKnown->GetBinContent(x, y);

        if(isBad){
          frac->SetBinContent(x, y, 0);
          fracXYLoose->SetBinContent(x, y, 0);
          fracXYTight->SetBinContent(x, y, 0);
          fracZLoose->SetBinContent(x, y, 0);
          fracZTight->SetBinContent(x, y, 0);
          continue;
        }

        const double fr = frac->GetBinContent(x, y);
        if(tot->GetBinContent(x, y)) proj->Fill(fr);

        const double frLoXY = fracXYLoose->GetBinContent(x, y);
        if(totXYLoose->GetBinContent(x, y)) projXYLoose->Fill(frLoXY);

        const double frTiXY = fracXYTight->GetBinContent(x, y);
        if(totXYTight->GetBinContent(x, y)) projXYTight->Fill(frTiXY);

        const double frLoZ = fracZLoose->GetBinContent(x, y);
        if(totZLoose->GetBinContent(x, y)) projZLoose->Fill(frLoZ);

        const double frTiZ = fracZTight->GetBinContent(x, y);
        if(totZTight->GetBinContent(x, y)) projZTight->Fill(frTiZ);
      } // end for Y
    } // end for X
  }

  DEFINE_ART_MODULE(CosmicEff)

} // end namespace chaninfo
///////////////////////////////////////////////////////////////////////////