CosmicTrack_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \brief   Track finder for cosmic rays
/// \author  brebel@fnal.gov, gsdavies@iastate.edu
/// \date
///////////////////////////////////////////////////////////////////////
#include <cmath>
#include <string>
#include <vector>

// ROOT includes
#include "TVector3.h"

// NOvA includes
#include "TrackFit/CosmicTrackAlg.h"
#include "TrackFit/WindowTrackingAlg.h"
#include "Calibrator/Calibrator.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Track.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"

// Framework includes
#include "art/Framework/Core/EDProducer.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"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"


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

    void produce(art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& p);
    
  private:
    std::string             fClusterInput;     ///< Input folder from cluster reco
    int                     fTrackAlgEnum;     ///< which algorithm to use 0: CosmicTrackAlg
                                               ///< 1: WindowTrackingAlg
    trk::TrackAlg          *fTrkAlg;           ///< Algorithm to carry out track fitting

  };
  
  //......................................................................
  CosmicTrack::CosmicTrack(fhicl::ParameterSet const& pset)
    : fTrkAlg(nullptr)
  {
    reconfigure(pset);

    produces< std::vector<rb::Track> >();
    // Associate Tracks with the Slice they came from
    produces< art::Assns<rb::Track, rb::Cluster> >();

    switch(fTrackAlgEnum){
      case 0:
        fTrkAlg = new trk::CosmicTrackAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
        break;
      case 1:
        fTrkAlg = new trk::WindowTrackingAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
        break;
      default:
        fTrkAlg = new trk::WindowTrackingAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
        break;
    }

  }
  
  //......................................................................
  CosmicTrack::~CosmicTrack()
  {
    if (fTrkAlg) delete fTrkAlg;
  }
  
  //......................................................................
  void CosmicTrack::reconfigure(const fhicl::ParameterSet& pset)
  {
    fClusterInput = pset.get< std::string >("ClusterInput"             );
    fTrackAlgEnum = pset.get< int         >("TrackingAlgEnumeration", 1); 
  }

  //......................................................................
  ///
  /// Reconstruction method for this module. Provides read-write access
  /// to the event data.
  ///
  /// \param evt - Read-write access to the event data
  ///
  void CosmicTrack::produce(art::Event& evt) 
  {
    // Declare a container for Track objects to be stored in the art::event 
    std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
    // Associate Tracks with the Slice they came from
    std::unique_ptr< art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);

    // define vector holding the clusters
    art::Handle< std::vector<rb::Cluster> > hSlices;
    evt.getByLabel(fClusterInput, hSlices);

    std::vector< art::Ptr<rb::Cluster> > slices;
    art::fill_ptr_vector(slices, hSlices);

    // now we have time slices of cellhits in clusters from Slicer- fit each
    // one and make the tracks
    const int sliceMax = slices.size();

    for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
      
      // ignore clusters that are tagged as noise
      if(slices[sliceIdx]->IsNoise() ) continue;
      
      LOG_DEBUG("CosmicTrack") << "On slice " << sliceIdx << " of " << sliceMax;

      // Use the cosmic track algorithm to build a track from the slice
      std::vector<rb::Track> tracks = fTrkAlg->MakeTrack(slices[sliceIdx].get());
      for(auto track : tracks){
        if ( track.NCell() == 0 ) continue; // failed to make a track
        trackcol->push_back(track);
          //for(auto const tp : track.Trajectory()) tp.Print();
        util::CreateAssn(*this, evt, *trackcol, slices[sliceIdx], *assns);
      }
    }// end loop over slices

    evt.put(std::move(trackcol));
    evt.put(std::move(assns));
  } // end produce

} // end trk namespace

namespace trk{

  DEFINE_ART_MODULE(CosmicTrack)

}