RockPresel_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
// \file    RockPresel_module.cc
// \brief   Preselection to filter rock muon events
// \author  Andrew Sutton
///////////////////////////////////////////////////////////////////////////
#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/Vertex.h"
#include "Preselection/PreselObj.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/func/NuMITimeParams.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 RockPresel: public art::EDProducer
  {
  public:
    explicit RockPresel(fhicl::ParameterSet const &pset);
    virtual ~RockPresel();
    void beginRun(art::Run& run);
    virtual void produce(art::Event& evt);
    
  protected:
    
    std::string  fSlicerLabel;
    std::string  fVertexLabel;
    bool         fRunFilter;
  
  private:
    bool isND = false;

  };
  
}//End namespace presel

namespace presel
{
  //----------------------------------------------------------------------
  RockPresel::RockPresel(fhicl::ParameterSet const &pset) :
    fSlicerLabel    (pset.get< std::string  >("SlicerLabel")),
    fVertexLabel    (pset.get< std::string  >("VertexLabel")  ),
    fRunFilter      (pset.get< bool         >("RunFilter"))
  {
    produces< rb::FilterList<rb::Cluster> >();
    produces< std::vector<PreselObj> >();
    produces< art::Assns<PreselObj, rb::Cluster> >();
  }
  
  //----------------------------------------------------------------------
  RockPresel::~RockPresel()
  {
  }
  
  //----------------------------------------------------------------------
  void RockPresel::beginRun(art::Run& run)
  {
    art::ServiceHandle<geo::Geometry>  geom;

    // Only apply the filter to ND events
    if(geom->DetId() == novadaq::cnv::kNEARDET)
      isND = true;
  }

  //----------------------------------------------------------------------
  void RockPresel::produce(art::Event& evt)
  {
    std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
    std::unique_ptr< std::vector<presel::PreselObj> > preselcol(new std::vector<presel::PreselObj>);
    std::unique_ptr< art::Assns<presel::PreselObj, rb::Cluster> > preselinfo(new art::Assns<presel::PreselObj, rb::Cluster>);

    
    // Grab event's slices
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSlicerLabel, slices);
    
    // Object holding tracks associated with slice
    art::FindManyP<rb::Vertex> fmpVertex(slices, evt, fVertexLabel);
    
    // Loop over slices
    for(uint sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      const rb::Cluster& slice = (*slices)[sliceIdx];

      if(slice.IsNoise()) continue;      

      // initialize presel object and set the default selection to true
      presel::PreselObj preselobj;
      bool sel = true;

      // check filter conditions
      // only check ND events
      if(isND && fmpVertex.isValid()){
        std::vector< art::Ptr<rb::Vertex> > vertexes = fmpVertex.at(sliceIdx);

        // assuming only one vertex
        if(vertexes.size()) {
          double vertX = vertexes[0]->GetX();
          double vertY = vertexes[0]->GetY();
          double vertZ = vertexes[0]->GetZ();

          sel = vertX < 179 && vertY < 179 && vertX > -179 && vertY > -179 && vertZ > 21;
        }
      }

      if(!isND && !sel) {
        std::cout << "\n*******************************************************"
                  << "\n   Not ND but sel is false, this shouldn't happen!!!   " 
                  << "\n*******************************************************"
                  << std::endl;
      }
      
      if(!sel && fRunFilter)
        filtcol->Add(slices, sliceIdx);

      preselobj.SetPassPresel(sel);      
      preselcol->push_back(preselobj);

      util::CreateAssn(*this, evt, *preselcol, art::Ptr<rb::Cluster>(slices,sliceIdx), *preselinfo);
      
    } // end of slice loop

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

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

namespace presel
{
  DEFINE_ART_MODULE(RockPresel)
}