Filter_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    Filter.cxx
// \brief   Filter events based on simple expressions
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "NovaDAQConventions/DAQConventions.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "Simulation/FLSHitList.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "Utilities/func/MathUtil.h"
#include "MCCheater/BackTracker.h"
// Framework includes
#include "art/Framework/Core/EDFilter.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 "TFormula.h"

namespace filter
{
  enum kPars{
    kNEvts,  ///< Number of true neutrino interactions
    kPdg,    ///< PDG code of the neutrino
    kCcNc,   ///< Value of Nu().CCNC()
    kIact,   ///< Value of Nu().InteractionType()
    kEnu,    ///< True neutrino energy (GeV)
    kX,      ///< True vertex x position (cm)
    kY,      ///< True vertex y position (cm)
    kZ,      ///< True vertex z position (cm)
    kEdge,   ///< Distance to detector edge (cm)
    kInsideFiducialVolume,  ///< Is event inside fiducial volume? (fiducial volume is defined in Geometry.fcl)
    kNhit,   ///< Total number of hits in the event
    kTotAdc, ///< Total ADC in the event
    kNSlice, ///< Total number of slices in the event (including noise slice)
    kSliceNXHit, ///< Number of hits in X view of the biggest physics slice
    kSliceNYHit, ///< Number of hits in X view of the biggest physics slice
    kVertexInFiducial, ///< Is track vertex (start position) inside the fiducial volume defined in filter
    kNFLShit,
    kNumPars
  };

  // This list must be kept in sync with the enum above. This sucks, but
  // apparently C++ doesn't have designated initialization (C does).
  const TString kParNames[kNumPars] = {"nevts", "pdg", "ccnc", "iact", "enu",
                                       "x", "y", "z", "edge", "is_inside_fiducial_volume",
                                       "nhit", "totadc",
                                       "nslice", "slicenxhit", "slicenyhit","vertex_in_fiducial","nFLShit"};

  /// Filter events based on simple expressions
  class Filter: public art::EDFilter
  {
  public:
    explicit Filter(const fhicl::ParameterSet& pset);
    ~Filter();

    bool filter(art::Event& evt);

    void reconfigure(const fhicl::ParameterSet& pset);

    void beginJob();

  protected:
    void FillTruthVariables(const art::Event& evt, double* vars) const;
    void FillDigitVariables(const art::Event& evt, double* vars) const;
    void FillG4GenVariables(const art::Event& evt, double* vars) const;
    void FillSliceVariables(const art::Event& evt, double* vars) const;
    void FillTrackVariables(const art::Event& evt, double* vars) const;

    bool isVertexInFiducial(const rb::Track& trk) const;

    TFormula* fFormula;

    std::string fGenLabel;     ///< Default "generator"
    std::string fG4GenLabel;   ///< Default "geantgen"
    std::string fCellHitLabel; ///< Default "daq"
    std::string fSlicerLabel;  ///< Default "slicer"
    std::string fTrackLabel;   ///< Default "kalmantrack"
  };

  //......................................................................
  Filter::Filter(const fhicl::ParameterSet& pset)
    : fFormula(0)
  {
    reconfigure(pset);
  }

  //......................................................................
  Filter::~Filter()
  {
    delete fFormula;
  }

  bool compareBySecondLength(const std::pair<int, std::string>& a,
                             const std::pair<int, std::string>& b)
  {
    return a.second.size() < b.second.size();
  }

  //......................................................................
  void Filter::reconfigure(const fhicl::ParameterSet& pset)
  {
    fGenLabel = pset.get<std::string>("GeneratorModuleLabel", "generator");
    fG4GenLabel = pset.get<std::string>("G4GenModuleLabel", "geantgen");
    fCellHitLabel = pset.get<std::string>("CalHitModuleLabel", "calhit");
    fSlicerLabel = pset.get<std::string>("SlicerModuleLabel", "slicer");
    fTrackLabel = pset.get<std::string>("TrackModuleLabel", "kalmantrack");

    TString expr = pset.get<std::string>("Expression");

    // Empty string evaluates true
    if(expr == "") expr = "1";

    // Need to sort the parameters from longest to shortest string, doing the
    // replacements in any other order causes everything to go wrong...
    std::vector<std::pair<int, std::string> > sortedPars;
    for(int n = 0; n < kNumPars; ++n){
      sortedPars.emplace_back(n, std::string(kParNames[n]));
    }
    std::sort(sortedPars.rbegin(), sortedPars.rend(), compareBySecondLength);
    
    // Rewrite the expression with numeric indices in place of variable names
    for(int n = 0; n < kNumPars; ++n){
      expr.ReplaceAll(sortedPars[n].second,
                      TString::Format("[%d]", sortedPars[n].first));
    }

    delete fFormula;
    fFormula = new TFormula("form", expr);
  }

  //......................................................................
  void Filter::beginJob()
  {
  }

  //......................................................................
  bool Filter::filter(art::Event& evt)
  {
    double vars[kNumPars] = {0};

    FillTruthVariables(evt, vars);
    FillG4GenVariables(evt, vars);
    FillDigitVariables(evt, vars);
    FillSliceVariables(evt, vars);
    FillTrackVariables(evt, vars);

    return fFormula->EvalPar(0, vars);
  }

  //......................................................................
  void Filter::FillTruthVariables(const art::Event& evt, double* vars) const
  {
    art::ServiceHandle<cheat::BackTracker> bt;
    if( !bt->HaveTruthInfo() ) return;

    art::Handle<std::vector<simb::MCTruth> > mcTruths;
    evt.getByLabel(fGenLabel, mcTruths);
    if(mcTruths.failedToGet()) return;    

    vars[kNEvts] = mcTruths->size();
    // Abort if no truth info (data?)
    if(mcTruths->empty()) return;
    // If there are two events, take the first one...
    simb::MCTruth mcTruth = (*mcTruths)[0];
    // Abort if not a neutrino-induced event (cosmics?)
    if(!mcTruth.NeutrinoSet()) return;

    const simb::MCNeutrino& nu = mcTruth.GetNeutrino();
    vars[kCcNc] = nu.CCNC();
    vars[kIact] = nu.InteractionType();

    const simb::MCParticle& nupart = nu.Nu();
    vars[kPdg] = nupart.PdgCode();
    vars[kEnu] = nupart.E();

    const simb::MCParticle& lept = nu.Lepton();
    const double x = lept.Vx();
    const double y = lept.Vy();
    const double z = lept.Vz();
    vars[kX] = x; vars[kY] = y; vars[kZ] = z;

    art::ServiceHandle<geo::Geometry> geom;
    const double dx = geom->DetHalfWidth()-fabs(x);                           // If outside, then dx is negative
    const double dy = geom->DetHalfHeight()-fabs(y);                          // If outside, then dy is negative
    const double dz = 0.5*geom->DetLength() - fabs(z + 0.5*geom->DetLength());// If outside, then dz is negative
    const double kEdge_mod = std::min(std::min(fabs(dx), fabs(dy)), fabs(dz));// minimum of the absolute values
    vars[kEdge] = (dx < 0.0 || dy <0.0 || dz<0.0) ? -kEdge_mod : kEdge_mod;   // If outside, then kEdge is negative

    /// Fill whether the point is inside fiducial volume
    vars[kInsideFiducialVolume] = geom->isInsideFiducialVolume(x, y, z) ? 1.0 : 0.0;
  }

 //......................................................................
  void Filter::FillDigitVariables(const art::Event& evt, double* vars) const
  {
    double totadc = 0;
    // get the RawDigit list
    art::Handle< std::vector<rb::CellHit> > digitcol;
    evt.getByLabel(fCellHitLabel, digitcol);
    if(digitcol.failedToGet()) return;    

    const unsigned int nhit = digitcol->size();
    for(unsigned int i = 0; i < nhit; ++i){
      art::Ptr<rb::CellHit> ch(digitcol, i);
      totadc += ch->ADC();
    }

    vars[kNhit] = nhit;
    vars[kTotAdc] = totadc;
  }

  //......................................................................

  void Filter::FillG4GenVariables(const art::Event& evt, double* vars) const
  {
    // get the FLSHit list
    art::Handle< std::vector<sim::FLSHitList> > FLScol;
    evt.getByLabel(fG4GenLabel, FLScol);
    if (FLScol.failedToGet()) return;
    
    const unsigned int nhit = FLScol->size();
    vars[kNFLShit] = nhit;
  }

  //......................................................................
  void Filter::FillSliceVariables(const art::Event& evt, double* vars) const
  {
    art::Handle<std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSlicerLabel, slices);
    if(slices.failedToGet()) return;

    const unsigned int sliceMax = slices->size();
    vars[kNSlice] = sliceMax;

    // Record the properties of the biggest non-noise slice
    unsigned int most = 0;
    for(unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;
      if(slice.NCell() > most){
        most = slice.NCell();
        vars[kSliceNXHit] = slice.NXCell();
        vars[kSliceNYHit] = slice.NYCell();
      }
    } // end for sliceIdx
  }

  //......................................................................
  void Filter::FillTrackVariables(const art::Event& evt, double* vars) const
  {
    art::Handle<std::vector<rb::Track> > trk;
    evt.getByLabel(fTrackLabel, trk);
    if(trk.failedToGet()) return;
    
    unsigned int trkMax = trk->size();
    // Record the properties of tracks inside a defined fiducial volume
    for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
      const rb::Track& track = (*trk)[trkIdx];
      //if(track.ExtentPlane() > 10) return;
      const TVector3 dir(track.Dir());
      const double coscosmic = TMath::Abs(-dir.Y()/util::pythag(dir.X(),dir.Y(),dir.Z()));
      if(coscosmic < 0.7) return;
      //if(track.TotalLength() > 70) return;
      const bool check = isVertexInFiducial(track);
      vars[kVertexInFiducial] = check;
    } // end for trkIdx
  }
  
  bool Filter::isVertexInFiducial(const rb::Track& trk) const
  {
    // Start and end don't mean much, pick whichever is highest in y.
    
    const TVector3 v1 = trk.Start();
    const TVector3 v2 = trk.Stop();
    // assume down-going in order to reject cosmics
    const TVector3 start = (v1.Y() > v2.Y()) ? v1 : v2;
    //const TVector3 start = v1;

    const double x = start.X();
    const double y = start.Y();
    const double z = start.Z();

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


    switch(geom->DetId()){
    case novadaq::cnv::kNDOS:
    case novadaq::cnv::kNEARDET:
      if(fabs(x) > 115) return false; // Fiducial X
      if(fabs(y) > 150) return false; // Fiducial Y
      if(z < 50 || z > 770) return false; // Fiducial Z

      if(v1.X()==v2.X()||v1.Y()==v2.Y()) return false; // dx or dy = 0 means bad fit
      break;
    case novadaq::cnv::kFARDET:
      if(fabs(x) > 730 || fabs(y) > 730) return false;
      if(z < 50 || z > 6300) return false;

      if(v1.X()==v2.X()||v1.Y()==v2.Y()) return false;
      break;
    default:
      assert(0 && "Unknown detector");
    }

    return true;
  }

  DEFINE_ART_MODULE(Filter)

} // namespace
////////////////////////////////////////////////////////////////////////