CosmicVeto_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
// \file    CosmicVeto_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"
// 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 CosmicVeto: public art::EDProducer
  {
  public:
    explicit CosmicVeto(fhicl::ParameterSet const &pset);
    virtual ~CosmicVeto();
    void beginRun(art::Run& run);
    virtual void produce(art::Event& evt);
    
  protected:
    
    std::string  fSlicerLabel;
    std::string  fTrackType;
    bool         fRunFilter;
     
    bool         fDoTrackAngleCut;
    bool         fDoNCellLowCut;
    bool         fDoNCellHighCut;
    bool         fDoContainment;
    bool         fDoOutsideCut;
    bool         fDoSlicePlanesCut;
    bool         fUseBetterCut;
    bool         fDoMinPCut;
    bool         fKeepNuMIInTime;

    unsigned int fNCellLow;
    unsigned int fNCellHigh;
    unsigned int fFwdCut;
    unsigned int fBakCut;
    unsigned int fMinPZ;
    double       fRatio;
    double       fAngVarCutValue;

    fhicl::ParameterSet fPSetNDOS, fPSetND, fPSetFD, fPSetTB;
  private:
    
    cosrej::CosRejFxs*         fCosRejFxs;
    
  };
  
}//End namespace presel

namespace presel
{
  //----------------------------------------------------------------------
  CosmicVeto::CosmicVeto(fhicl::ParameterSet const &pset) :
    fSlicerLabel    (pset.get< std::string  >("SlicerLabel")),
    fTrackType      (pset.get< std::string  >("TrackType")  ),
    fRunFilter      (pset.get< bool         >("RunFilter")),
    fDoTrackAngleCut(0),
    fDoNCellLowCut  (0),
    fDoNCellHighCut (0),
    fDoContainment  (0),
    fDoOutsideCut   (0),
    fDoSlicePlanesCut(0),
    fUseBetterCut   (0),
    fDoMinPCut      (0),
    fKeepNuMIInTime (0),
    fNCellLow       (0),
    fNCellHigh      (0),
    fFwdCut         (0),
    fBakCut         (0),
    fMinPZ          (0),
    fRatio          (0),
    fAngVarCutValue (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> >();
    produces< std::vector<Veto> >();
    produces< art::Assns<Veto, rb::Cluster> >();
  }
  
  //----------------------------------------------------------------------
  CosmicVeto::~CosmicVeto()
  {
  }
  
  //----------------------------------------------------------------------
  void CosmicVeto::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;
        break;
      case novadaq::cnv::kTESTBEAM:
        pset = fPSetTB;
        break;
      default:
        assert(0 && "Unknown detector");
    }
    fDoTrackAngleCut = pset.get< bool         >("DoTrackAngleCut");
    fDoNCellLowCut   = pset.get< bool         >("DoNCellLowCut");
    fDoNCellHighCut  = pset.get< bool         >("DoNCellHighCut");
    fDoContainment   = pset.get< bool         >("DoContainment");
    fDoOutsideCut    = pset.get< bool         >("DoOutsideCut");
    fDoSlicePlanesCut= pset.get< bool         >("DoSlicePlanesCut"); 

    fUseBetterCut    = pset.get< bool         >("UseBetterCut");
    fDoMinPCut       = pset.get< bool         >("DoMinPCut");
    fKeepNuMIInTime  = pset.get< bool         >("KeepNuMIInTime");
    fNCellLow        = pset.get< unsigned int >("NCellLow");
    fNCellHigh       = pset.get< unsigned int >("NCellHigh");
    fFwdCut          = pset.get< unsigned int >("FwdCut");
    fBakCut          = pset.get< unsigned int >("BakCut");
    fMinPZ           = pset.get< unsigned int >("MinPZ");
    fRatio           = pset.get< double       >("Ratio");
    fAngVarCutValue  = pset.get< double       >("AngVarCutValue");
    std::cout << "CosmicVeto settings: angle cut: "; // add one line of config to log file
    if(fUseBetterCut) std::cout << "v2, throughgoing cut: v1, nhit cut: ";
    if(!fUseBetterCut) std::cout << "v1, throughgoing cut: v1, nhit cut: ";
    if(!fDoNCellLowCut) std::cout << "none" << std::endl;
    if(fDoNCellLowCut) std::cout << fNCellLow << std::endl;
    //if(fKeepNuMIInTime) 
    if(fKeepNuMIInTime)
      std::cout<<"Not applying cuts to events in time with NuMI\n";
  }

  //----------------------------------------------------------------------
  void CosmicVeto::produce(art::Event& evt)
  {
    std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
    
    std::unique_ptr< std::vector<presel::Veto> > vetocol(new std::vector<presel::Veto>);
    
    std::unique_ptr< art::Assns<presel::Veto, rb::Cluster> > vetoinfo(new art::Assns<presel::Veto, 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::Track> fmpTrack(slices, evt, fTrackType);
    
    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];
      
      Veto sliceinfo;
      if(slice.IsNoise()) continue;
      
      float fwddist = 100.0, bakdist = 100.0;
      double ratio = 1, diry = 0, angle = 1;
      bool passanglenew=true, passangleold=true, passthru=true;
      
      const double nCells = slice.NCell(); //Number of cells in slice
      sliceinfo.SetNCell(nCells);
      // Selected unless we fail any cuts below
      bool sel = true;

      if(fmpTrack.isValid()){
        std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(sliceIdx);
        if(tracks.size()>0) {
          const art::Ptr<rb::Track> track = tracks[0];

          art::ServiceHandle<geo::LiveGeometry>  livegeom;

          // get all the important track information
          TVector3 start = track->Start();
          TVector3 end = track->Stop();
          TVector3 dirfor = track->Dir().Unit();
          TVector3 dirbak = -dirfor;

          fwddist = livegeom->ProjectedLiveDistanceToEdge(end,dirfor);
          bakdist = livegeom->ProjectedLiveDistanceToEdge(start,dirbak);
          sliceinfo.SetFwdDist(fwddist);
          sliceinfo.SetBakDist(bakdist);
          diry = track->Dir().Y();
          sliceinfo.SetDirY(diry);
          angle = fCosRejFxs->getAngle(track->Dir());
          sliceinfo.SetAngle(angle);
          ratio = (track->NCell()*1.0) / (1.0*nCells);
          sliceinfo.SetRatio(ratio);
          double anglediff = (angle*angle) - (diry*diry);
          double anglevar = fabs(angle) * (diry+1.0);
          sliceinfo.SetAngleVarOld(anglediff);
          sliceinfo.SetAngleVar(anglevar);

          // would this pass the old, flawed, through-going check?
          if(fwddist < fFwdCut &&
             bakdist < fBakCut) 
            passthru = false;
          sliceinfo.SetPassThruOld(passthru);

          // how about the fixed cut we actually use?
          if(!passthru && 
             ratio > 0.8) 
            sliceinfo.SetPassThru(false);

          if(fDoContainment){
            // apply cut if asked for (fixed only)
            if(!passthru && ratio > 0.8) sel = false;
          }

          // would this pass the old angle cut used for first analysis?
          if(fmpTrack.isValid()) {
            if(ratio > 0.88) {
              if(slice.NCell() > 50 && slice.NCell() <= 135 &&
                 anglediff < ((1.1/85)*slice.NCell()-1.65)) passangleold = false;
              if(slice.NCell() > 135 && slice.NCell() <= 250 &&
                 anglediff < ((.6/115)*slice.NCell()-.6)) passangleold = false;
              if(slice.NCell() > 250 && slice.NCell() <= 500 &&
                 anglediff < ((.1/250)*slice.NCell()+0.6)) passangleold = false;
              if(slice.NCell() > 500 && anglediff < 0.8) passangleold = false;
            }
            sliceinfo.SetPassAngleFirstAna(passangleold);

            // how about the latest version of the angle cut?
            if( !fDoOutsideCut){
              if(ratio > 0.75       && 
                 slice.MaxY() > 400 && 
                 anglevar < fAngVarCutValue) 
                passanglenew = false;
            }
            else{
              // For nue cosmic veto, we want the angle cut applied
              // to not only the events that are close to the dete-
              // -ctor top, but close to any detector edge.
              // see doc 14654 
              
              // For Z containment cut, use LiveGeometry otherwise
              // cosmic rejection efficiency drops for fewer numb-
              // -ers of diblocks.
              double zDistToFront = livegeom->DistToFront(slice.MinXYZ());
              double zDistToBack  = livegeom->DistToBack(slice.MaxXYZ());

              bool fromOutside = (slice.MaxY() > 400  ||
                                  slice.MinY() < -700 ||
                                  slice.MaxX() > 650  ||
                                  slice.MinX() < -650 ||
                                  zDistToFront < 100  || 
                                  zDistToBack  < 150 );
              
              if(ratio > 0.75 && 
                 anglevar < fAngVarCutValue &&
                 fromOutside) 
                passanglenew = false;
            }
            sliceinfo.SetPassAngle(passanglenew);
          }

          // do the track angle cut if asked for:

          if (fDoTrackAngleCut){
            if(!fUseBetterCut) { 
              // old, first analysis tuning
              if(!passangleold) sel = false;
            }
            if(fUseBetterCut) { 
              // better tuning, much better cosmic removal, much simpler
              if(!passanglenew) sel = false;
            }
          }
          
        } // end loop over tracks (should only be one)
      } // end if track object is valid

      // do the NCell cuts if asked for:
      if(fDoNCellLowCut){
        if(nCells < fNCellLow) sel = false;
      }
      if(fDoNCellHighCut){
        if(nCells > fNCellHigh) sel = false;
      }
      if( fDoSlicePlanesCut ){
        // see doc 14654 
        if( slice.ExtentPlane(geo::kXorY) > 125 ||
            slice.ExtentPlane(geo::kXorY) < 5)
          sel = false;
      }
      // do ND (simulated half det, bad for full ND)
      if(fDoMinPCut){
        TVector3 minp = slice.MinXYZ();
        if(minp.Z() < fMinPZ && ratio > fRatio) sel = false;
      }

      // record whether or not we should keep this slice
      // before checking the timing cuts. Timing cuts are
      // tracked separately through the PassNuMICut flag.
      sliceinfo.SetKeep(sel);

      if( util::IsInBeamWindow(evt.run(), slice.MeanTNS()) ){
        
        sliceinfo.SetPassNuMICut(true);

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

      vetocol->push_back(sliceinfo);
      util::CreateAssn(*this,evt,*vetocol,slicePtr,*vetoinfo);
      
    } // end for sliceIdx

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

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

namespace presel
{
  DEFINE_ART_MODULE(CosmicVeto)
}