Preselection_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
// \file    Preselection_module.cc
// \brief   Preselection to skip events that fail cuts for optimisation
// \author  Christopher Backhouse - bckhouse@caltech.edu
///////////////////////////////////////////////////////////////////////////

#include "Calibrator/Calibrator.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/GeometryBase.h"

#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/PID.h"

#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"

#include "Utilities/AssociationUtil.h"

#include "art/Framework/Core/EDProducer.h"
#include "canvas/Persistency/Common/FindManyP.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/Framework/Services/Optional/TFileService.h"

namespace lem
{
  class Preselection: public art::EDProducer
  {
  public:
    explicit Preselection(const fhicl::ParameterSet& pset);
    virtual ~Preselection();

    virtual void produce(art::Event& evt);

  protected:
    std::string fSlicerLabel;

    bool   fTruthMode;
    bool   fSkipTrackCut;
    bool   fSkipEnergyCut;
    double fSliceNCellLow;
    double fSliceNCellHigh;
    double fSliceExtentPlane;
    double fTrackLengthCut;
    double fRecoELow;
    double fRecoEHigh;
    bool   fInFidCut;
  };

  //----------------------------------------------------------------------
  Preselection::Preselection(const fhicl::ParameterSet& pset):
    fSlicerLabel(pset.get<std::string>("SlicerLabel")),
    fTruthMode(pset.get<bool>("TruthMode")),
    fSkipTrackCut(pset.get<bool>("SkipTrackCut")),
    fSkipEnergyCut(pset.get<bool>("SkipEnergyCut")),
    fSliceNCellLow(pset.get<double>("SliceNCellLow")),
    fSliceNCellHigh(pset.get<double>("SliceNCellHigh")),
    fSliceExtentPlane(pset.get<double>("SliceExtentPlane")),
    fTrackLengthCut(pset.get<double>("TrackLengthCut")),
    fRecoELow(pset.get<double>("RecoELow")),
    fRecoEHigh(pset.get<double>("RecoEHigh")),
    fInFidCut(pset.get<bool>("InFidCut"))
  {

    produces<rb::FilterList<rb::Cluster> >();
    produces<std::vector<rb::PID> >();
    produces<art::Assns<rb::PID, rb::Cluster> >();

  }

  //----------------------------------------------------------------------
  Preselection::~Preselection()
  {
  }

  //----------------------------------------------------------------------
  void Preselection::produce(art::Event& evt)
  {
    std::unique_ptr<rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);

    std::unique_ptr<std::vector<rb::PID> > pidcol(new std::vector<rb::PID>);
    std::unique_ptr<art::Assns<rb::PID, rb::Cluster> > assns(new art::Assns<rb::PID, rb::Cluster>);

    art::ServiceHandle<calib::Calibrator> cal;
    art::ServiceHandle<geo::Geometry> geom;

    // TODO: break truth selection out into its own module
    /*
    if(fTruthMode){
      art::Handle<std::vector<simb::MCTruth> > truths;
      evt.getByLabel("generator", truths);
      if(truths->empty()) return false;
      if(truths->size() > 1) return true;

      const simb::MCTruth& truth = (*truths)[0];

      const simb::MCNeutrino& nu = truth.GetNeutrino();
      const double trueE = nu.Nu().E();
      // These cuts are empirically observed to remove almost nothing that
      // would survive the preselection below. Approximately 0.1% lost to the
      // low cuts and 0.5% to the high cuts.
      if(nu.CCNC() == simb::kCC){
        return (trueE > 1 && trueE < 3.5);
      }
      else{
        const double hadY = nu.Y();
        return (trueE*hadY > .75 && trueE*hadY < 4.5);
      }

      assert(0 && "Not reached");
    }
    */

    art::Handle<std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSlicerLabel, slices);

    const int sliceMax = slices->size();
    for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
      if(rb::IsFiltered(evt, slices, sliceIdx)) continue;

      const rb::Cluster& slice = (*slices)[sliceIdx];

      if(slice.IsNoise()) continue;

      double longestTrack = 0;
      if(!fSkipTrackCut){
        // This is maybe inefficient, but for now FindManyP throws an exception
        // on finding no assns. We don't want that.
        art::FindManyP<rb::Track> fmt(slices, evt, std::string("discretetrack"));
        // Get all the tracks associated with this slice
        std::vector<art::Ptr<rb::Track> > tracks = fmt.at(sliceIdx);

        for(unsigned int n = 0; n < tracks.size(); ++n){
          const double L = tracks[n]->TotalLength();
          if(L > longestTrack) longestTrack = L;
        }
      }

      const double nCells = slice.NCell();
      const int lengthPlanes = slice.ExtentPlane();

      const TVector3 mean = slice.MeanXYZ();
      const bool inFid = geom->isInsideFiducialVolume(mean.X(), mean.Y(), mean.Z());

      const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
      const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;

      double totalRecoGeV = 0;
      for(unsigned int hitIdx = 0; hitIdx < slice.NCell(); ++hitIdx){
        // This is slow because we have to keep recalculating the mean
        // positions. The manual version below is much better.
        // const rb::RecoHit rhit = slice.RecoHit(hitIdx);

        const art::Ptr<rb::CellHit>& chit = slice.Cell(hitIdx);
        const rb::RecoHit rhit = cal->MakeRecoHit(*chit,
                                                  chit->View() == geo::kX ? wx : wy);
        if(!rhit.IsCalibrated()){
          std::cout << "Not calibrated?!" << std::endl;
          continue;
        }
        totalRecoGeV += rhit.GeV();
      }

      bool sel = true;
      if(nCells < fSliceNCellLow || nCells > fSliceNCellHigh) sel = false;
      if(longestTrack > fTrackLengthCut) sel = false;
      if(lengthPlanes > fSliceExtentPlane) sel = false;
      if(!fSkipEnergyCut) {
        if(2*totalRecoGeV < fRecoELow || 2*totalRecoGeV > fRecoEHigh) sel = false;
      }
      if(!inFid && fInFidCut == true) sel = false;

      if(!sel) filtcol->Add(slices, sliceIdx);

      // We're a nue preselection
      pidcol->push_back(rb::PID(12, sel));

      util::CreateAssn(*this, evt, *pidcol, art::Ptr<rb::Cluster>(slices, sliceIdx), *assns);
    } // end for sliceIdx

    evt.put(std::move(pidcol));
    evt.put(std::move(assns));
    evt.put(std::move(filtcol));
  }

  DEFINE_ART_MODULE(Preselection)

} // end namespace