ThroughgoingSelection_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    ThroughgoingSelection.cxx
// \brief   Selects throughgoing cosmics
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDFilter.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

// NOvASoft includes
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "MEFinder/MEClusters.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"

namespace calib
{
  /// Selects throughgoing cosmics
  class ThroughgoingSelection: public art::EDFilter
  {
  public:
    explicit ThroughgoingSelection(const fhicl::ParameterSet& pset);
    ~ThroughgoingSelection();

    virtual bool filter(art::Event& evt);

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

    void beginJob();
    void endJob();

  protected:
    bool IsUncontained(const rb::Track& trk) const;
    bool IsTrulyUncontained(const rb::Track& trk) const;

    int fTable[2][2];

  };

  //......................................................................
  ThroughgoingSelection::ThroughgoingSelection(const fhicl::ParameterSet& pset)
  {
    produces< std::vector<rb::Track> >();

    reconfigure(pset);

    for(int i = 0; i < 2; ++i) for(int j = 0; j < 2; ++j) fTable[i][j] = 0;
  }

  //......................................................................
  ThroughgoingSelection::~ThroughgoingSelection()
  {
  }

  //......................................................................
  void ThroughgoingSelection::reconfigure(const fhicl::ParameterSet& pset)
  {
  }

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

  //......................................................................
  bool ThroughgoingSelection::filter(art::Event& evt)
  {
    std::unique_ptr<std::vector<rb::Track> > outTracks(new std::vector<rb::Track>);

    art::Handle<std::vector<rb::Track> > tracks;
    evt.getByLabel("cosmictrack", tracks);

    art::FindManyP<me::TrkME> fmp(tracks, evt, "mefinder");

    for(size_t trIdx = 0; trIdx < tracks->size(); ++trIdx){
      // If there's any Michel associated, don't take this track
      if(!fmp.at(trIdx).empty()) continue;

      const rb::Track& track = (*tracks)[trIdx];

      const bool thru = IsUncontained(track);

      if(thru) outTracks->push_back(track);

      if(!evt.isRealData()){
        const bool trueThru = IsTrulyUncontained(track);
        ++fTable[thru][trueThru];
      }
    } // end for trIdx

    const bool ret = !outTracks->empty();
    evt.put(std::move(outTracks));
    return ret;
  }

  //......................................................................
  void ThroughgoingSelection::endJob()
  {
    if(fTable[0][0]+fTable[1][0] == 0 || fTable[0][1]+fTable[1][1] == 0) return;
    std::cout << "Throughgoing efficiency table. x-axis selection, y-axis truth:" << std::endl;
    std::cout << fTable[0][0] << "\t" << fTable[1][0] << "\t" << .1*(1000*fTable[0][0]/(fTable[0][0]+fTable[1][0])) << "%" << std::endl;
    std::cout << fTable[0][1] << "\t" << fTable[1][1] << "\t" << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl;
    std::cout << .1*(1000*fTable[0][0]/(fTable[0][0]+fTable[0][1])) << "%\t" << .1*(1000*fTable[1][1]/(fTable[1][0]+fTable[1][1])) << "%" << std::endl;
    std::cout << "Purity = " << .1*(1000*fTable[1][1]/(fTable[1][0]+fTable[1][1])) << "%" << std::endl;
    std::cout << "Efficiency = " << .1*(1000*fTable[1][1]/(fTable[0][1]+fTable[1][1])) << "%" << std::endl;
  }

  //......................................................................
  bool ThroughgoingSelection::IsUncontained(const rb::Track& trk) const
  {
    // Start and end don't mean much, pick whichever is lowest in y.
    const TVector3 v1 = trk.Start();
    const TVector3 v2 = trk.Stop();
    const TVector3 stop = (v1.Y() < v2.Y()) ? v1 : v2;

    if(v1.X() == v2.X() ||
       v1.Y() == v2.Y()) return false; // dx or dy = 0 means bad fit

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

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

    if(fabs(x) > geom->DetHalfWidth()-20) return true;
    if(fabs(y) > geom->DetHalfHeight()-20) return true;
    if(z < 20 || z > geom->DetLength()-20) return true;

    return false;
  }

  //......................................................................
  bool ThroughgoingSelection::
  IsTrulyUncontained(const rb::Track& trk) const
  {
    art::ServiceHandle<cheat::BackTracker> bt;

    std::vector<const sim::Particle*> parts = bt->HitsToParticle(trk.AllCells());
    if(parts.empty()) return false;

    const sim::Particle* part = parts.front();
    if(abs(part->PdgCode()) != 13) return false;

    const TVector3 stop = part->EndPosition().Vect();

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

    if(fabs(stop.X()) > geom->DetHalfWidth()) return true;
    if(fabs(stop.Y()) > geom->DetHalfHeight()) return true;
    if(stop.Z() < 0) return true;
    if(stop.Z() > geom->DetLength()) return true;

    return false;
  }

  DEFINE_ART_MODULE(ThroughgoingSelection)

} // end namespace
////////////////////////////////////////////////////////////////////////