CosmicCVNVeto_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
// \file    CosmicCVNVeto_module.cc
// \brief   Preselection to filter cosmic events for reconstruction/PID
//          for everything
// \author  Gavin Davies, Kirk Bays
///////////////////////////////////////////////////////////////////////////
#include "TMath.h"
#include "TVector3.h"

//Novasoft includes
#include "CosRej/CosRejFxs.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "MCCheater/BackTracker.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Track.h"
#include "Preselection/Veto.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/func/NuMITimeParams.h"
#include "TensorFlowProducts/CVNCosmicFiltList.h"
// Framework includes
#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/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "canvas/Persistency/Common/Assns.h"

namespace presel
{
  class CosmicCVNVeto: public art::EDProducer
  {
  public:
    explicit CosmicCVNVeto(fhicl::ParameterSet const &pset);
    virtual ~CosmicCVNVeto();
    void beginRun(art::Run& run);
    virtual void produce(art::Event& evt);
    
  protected:
    
    std::string  fSlicerLabel;
    std::string  fCosmicCVNLabel;
    bool         fRunFilter;
     
    bool         fDoNCellHighCut;
    bool         fDoCosmicCVNCut;
    bool         fKeepNuMIInTime;
    bool         isFD = false;
    unsigned int fNCellHigh;

    fhicl::ParameterSet fPSetNDOS, fPSetND, fPSetFD, fPSetTB;
  private:
    
    
  };
  
}//End namespace presel

namespace presel
{
  //----------------------------------------------------------------------
  CosmicCVNVeto::CosmicCVNVeto(fhicl::ParameterSet const &pset) :
    fSlicerLabel    (pset.get< std::string  >("SlicerLabel")),
    fCosmicCVNLabel (pset.get< std::string  >("CosmicCVNLabel")  ),
    fRunFilter      (pset.get< bool         >("RunFilter")),
    fDoNCellHighCut (0),
    fDoCosmicCVNCut (0),
    fKeepNuMIInTime (0),
    fNCellHigh      (0),
    fPSetNDOS       (pset.get<fhicl::ParameterSet>("ndos")),
    fPSetND         (pset.get<fhicl::ParameterSet>("nd")),
    fPSetFD         (pset.get<fhicl::ParameterSet>("fd")),
    fPSetTB         (pset.get<fhicl::ParameterSet>("tb"))
  
  {
    produces< rb::FilterList<rb::Cluster> >();
  }
  
  //----------------------------------------------------------------------
  CosmicCVNVeto::~CosmicCVNVeto()
  {
  }
  
  //----------------------------------------------------------------------
  void CosmicCVNVeto::beginRun(art::Run& run)
  {
    art::ServiceHandle<geo::Geometry>  geom;
    
    fhicl::ParameterSet pset;
    
    switch(geom->DetId()){
      case novadaq::cnv::kNDOS:
        pset = fPSetNDOS;
        break;
      case novadaq::cnv::kNEARDET:
        pset = fPSetND;
        break;
      case novadaq::cnv::kFARDET:
        pset = fPSetFD;
        isFD = true;
        break;
      case novadaq::cnv::kTESTBEAM:
        pset = fPSetTB;
        break;
      default:
        assert(0 && "Unknown detector");
    }
    fDoNCellHighCut  = pset.get< bool         >("DoNCellHighCut");
    fDoCosmicCVNCut  = pset.get< bool         >("DoCosmicCVNCut");
    fKeepNuMIInTime  = pset.get< bool         >("KeepNuMIInTime");
    fNCellHigh       = pset.get< unsigned int >("NCellHigh");
    if(fKeepNuMIInTime)
      std::cout<<"Not applying cuts to events in time with NuMI\n";
  }

  //----------------------------------------------------------------------
  void CosmicCVNVeto::produce(art::Event& evt)
  {
    std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
   
    if (!isFD){
      evt.put(std::move(filtcol));
      return;
    }
  
    
    //Grab event's slices
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSlicerLabel, slices);

    //Grab the cosmic CVN filter objects
    art::Handle<std::vector<cvntf::CVNCosmicFiltList> > cosmiccvn;
    evt.getByLabel(fCosmicCVNLabel, cosmiccvn);
    std::vector<cvntf::CVNCosmicFiltList> cosmicList = *cosmiccvn;

    
    const int sliceMax = slices->size();
    
    for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
      art::Ptr<rb::Cluster> slicePtr(slices,sliceIdx); // Used to write association
      
      const rb::Cluster& slice = (*slices)[sliceIdx];
      
      if(slice.IsNoise()) continue;
      
      
      const double nCells = slice.NCell(); //Number of cells in slice
      // Selected unless we fail any cuts below
      bool sel = true;

      if(fDoNCellHighCut){
        if(nCells > fNCellHigh) sel = false;
      }
      if( fDoCosmicCVNCut ){
        double cVal = 999.0;
        bool keep = false;
        for (unsigned int iC = 0; iC < cosmicList[0].ListSize(); ++iC){
          cvntf::CVNCosmicFilt cvn = cosmicList[0].GetTimeSlice(iC);
          if ((slice.MeanTNS() > cvn.timeWinMin) &&
             (slice.MeanTNS() < cvn.timeWinMax)){
            if (cvn.cosmicVal < cVal){
              cVal = cvn.cosmicVal;
              keep = cvn.passSel;
            }
          }
        }
        if (!keep) sel = false;
      }

      if( util::IsInBeamWindow(evt.run(), slice.MeanTNS()) ){
        
        if(fKeepNuMIInTime){
          // if this slice is in the numi beam window, 
          // keep it.
          sel = true;     
        }
      }
      
      if(!sel && fRunFilter)
        filtcol->Add(slices, sliceIdx);

    } // end for sliceIdx

    evt.put(std::move(filtcol));

  }
  //............................................................................................
}//end namespace presel

namespace presel
{
  DEFINE_ART_MODULE(CosmicCVNVeto)
}