DumpTracks_module.cc

Go to the documentation of this file.

// LArSoft includes
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/PFParticle.h"

// support libraries
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Comment.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// art libraries
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Utilities/InputTag.h"

// C//C++ standard libraries
#include <string>
#include <sstream>
#include <iomanip> // std::setw()
#include <algorithm> // std::max(), std::sort(), std::transform()
#include <iterator> // std::back_inserter()
#include <functional> // std::mem_fn()
#include <memory> // std::unique_ptr()


namespace {
  
  template <typename T>
  size_t StringLength(const T& value) {
    std::ostringstream sstr;
    sstr << value;
    return sstr.str().length();
  } // StringLength()
  
  template <typename STREAM, typename CONT>
  void PrintCompactIndexTable(
    STREAM& log, const CONT& Indices, unsigned int IndicesPerLine,
    std::string IndentStr = "")
  {
    unsigned int Padding = StringLength(Indices.back());
    
    typename CONT::const_iterator iIndex = Indices.begin(),
      iend = Indices.end();
    size_t RangeStart = *iIndex, RangeStop = RangeStart;
    std::ostringstream output_line;
    size_t nItemsInLine = 0;
    while (++iIndex != iend) {
      
      if (*iIndex == RangeStop + 1) {
        ++RangeStop;
      }
      else {
        // the new item does not belong to the current range:
        // - print the current range
        if (nItemsInLine) output_line << "  ";
        if (RangeStart == RangeStop) {
          output_line << std::setw(Padding) << RangeStart;
          ++nItemsInLine;
        }
        else {
          char fill = (RangeStart + 1 == RangeStop)? ' ': '-';
          output_line << std::setw(Padding) << RangeStart
            << fill << fill
            << std::setw(Padding) << std::setfill(fill) << RangeStop
            << std::setfill(' ');
          nItemsInLine += 2;
        }
        // - start a new one
        RangeStart = RangeStop = *iIndex;
      } // if ... else
      
      // if we have enough stuff in the buffer, let's print it
      if (nItemsInLine >= IndicesPerLine) {
        nItemsInLine = 0;
        log << IndentStr << output_line.str() << "\n";
        output_line.str("");
      }
      
    } // while
    
    // print what we have accumulated so far of the last line
    log << IndentStr << output_line.str();
    // print the last range (or single element)
    if (nItemsInLine) log << "  ";
    if (RangeStart == RangeStop)
      log << std::setw(Padding) << RangeStart;
    else {
      char fill = (RangeStart + 1 == RangeStop)? ' ': '-';
      log << std::setw(Padding) << RangeStart
        << fill << fill
        << std::setw(Padding) << std::setfill(fill) << RangeStop
        << std::setfill(' ');
    }
    log << std::endl;
  } // PrintCompactIndexTable()
  
  
  
  template <typename STREAM, typename CONT, typename GETINDEX>
  void PrintCompactIndexTable(
    STREAM& log, const CONT& Objects, unsigned int IndicesPerLine,
    GETINDEX IndexExtractor, std::string IndentStr)
  {
    if ((IndicesPerLine == 0) || Objects.empty()) return;
    
    std::vector<size_t> Indices;
    Indices.reserve(Objects.size());
    std::transform(Objects.begin(), Objects.end(), std::back_inserter(Indices),
      IndexExtractor);
    std::sort(Indices.begin(), Indices.end());
    PrintCompactIndexTable(log, Indices, IndicesPerLine, IndentStr);
    
  } // PrintCompactIndexTable()
  
  template <typename STREAM, typename T>
  inline void PrintAssociatedIndexTable(
    STREAM& log, const std::vector<art::Ptr<T>>& Objects,
    unsigned int IndicesPerLine, std::string IndentStr = ""
  ) {
    PrintCompactIndexTable(
      log, Objects, IndicesPerLine, std::mem_fn(&art::Ptr<T>::key),
      IndentStr
      );
  } // PrintAssociatedIndexTable()
  
  
  template <typename STREAM, typename T>
  inline void PrintAssociatedIDTable(
    STREAM& log, const std::vector<art::Ptr<T>>& Objects,
    unsigned int IndicesPerLine, std::string IndentStr = ""
  ) {
    PrintCompactIndexTable(
      log, Objects, IndicesPerLine,
      [](const art::Ptr<T>& ptr) { return ptr->ID(); },
      IndentStr
      );
  } // PrintAssociatedIDTable()
  
} // local namespace


namespace recob {

  class DumpTracks : public art::EDAnalyzer {
      public:

    struct Config {
      using Comment = fhicl::Comment;
      using Name = fhicl::Name;
      
      fhicl::Atom<art::InputTag> TrackModuleLabel{
        Name("TrackModuleLabel"),
        Comment("input tag for the tracks to be dumped")
        };
      fhicl::Atom<std::string> OutputCategory{
        Name("OutputCategory"),
        Comment("name of the category used for message facility output"),
        "DumpTracks"
        };
      fhicl::Atom<unsigned int> WayPoints{
        Name("WayPoints"),
        Comment("number of points along the trajectory printed"),
        10U
        };
      fhicl::Atom<bool> SpacePointAssociations{
        Name("SpacePointAssociations"),
        Comment("prints the number of space points associated to the track"),
        true
        };
      fhicl::Atom<bool> PrintSpacePoints{
        Name("PrintSpacePoints"),
        Comment("prints the index of all space points associated to the track"),
        false
        };
      fhicl::Atom<bool> HitAssociations{
        Name("HitAssociations"),
        Comment("prints the number of hits associated to the track"),
        true
        };
      fhicl::Atom<bool> PrintHits{
        Name("PrintHits"),
        Comment("prints the index of all hits associated to the track"),
        false
        };
      fhicl::Atom<bool> ParticleAssociations{
        Name("ParticleAssociations"),
        Comment("prints the number of PF particles associated to the track"),
        true
        };
      
    }; // Config
    
    using Parameters = art::EDAnalyzer::Table<Config>;
    
    explicit DumpTracks(Parameters const& config); 
    
    void analyze (const art::Event& evt);

      private:

    art::InputTag fTrackModuleLabel; 
    std::string fOutputCategory; 
    unsigned int fPrintWayPoints; 
    
    bool fPrintNHits; 
    bool fPrintNSpacePoints; 
    bool fPrintNParticles; 
    bool fPrintHits; 
    bool fPrintSpacePoints; 
    bool fPrintParticles; 

    void DumpTrack(unsigned int iTrack, recob::Track const& track) const;
    
  }; // class DumpTracks

} // namespace recob


//------------------------------------------------------------------------------
namespace {
  
  template <typename STREAM>
  bool PrintdQdXinView(
    STREAM& log,
    const recob::Track& track, geo::View_t view, std::string ViewName
  ) {
    auto nData = 0;
    try { nData = track.NumberdQdx(view); }
    catch (std::out_of_range) { return false; }
    log << " " << ViewName << ": " << nData;
    return true;
  } // PrintdQdXinView()
  
} // local namespace


namespace recob {

  //-------------------------------------------------
  DumpTracks::DumpTracks(Parameters const& config) 
    : EDAnalyzer        (config) 
    , fTrackModuleLabel (config().TrackModuleLabel())
    , fOutputCategory   (config().OutputCategory())
    , fPrintWayPoints   (config().WayPoints())
    , fPrintNHits       (config().HitAssociations())
    , fPrintNSpacePoints(config().SpacePointAssociations())
    , fPrintNParticles  (config().ParticleAssociations())
    , fPrintHits        (config().PrintHits())
    , fPrintSpacePoints (config().PrintSpacePoints())
    , fPrintParticles   (config().ParticleAssociations())
    {}
  
  //-------------------------------------------------
  void DumpTracks::analyze(const art::Event& evt) {
    
    // fetch the data to be dumped on screen
    auto Tracks
      = evt.getValidHandle<std::vector<recob::Track>>(fTrackModuleLabel);
    
    mf::LogInfo(fOutputCategory)
      << "The event contains " << Tracks->size() << " '"
      << fTrackModuleLabel.encode() << "'tracks";
    
    std::unique_ptr<art::FindManyP<recob::Hit>> pHits(
      fPrintNHits?
      new art::FindManyP<recob::Hit>(Tracks, evt, fTrackModuleLabel):
      nullptr
      );
    if (pHits && !pHits->isValid()) {
      throw art::Exception(art::errors::ProductNotFound)
        << "No hit associated with '" << fTrackModuleLabel.encode()
        << "' tracks.\n";
    }
    
    std::unique_ptr<art::FindManyP<recob::SpacePoint>> pSpacePoints(
      fPrintNSpacePoints?
      new art::FindManyP<recob::SpacePoint>(Tracks, evt, fTrackModuleLabel):
      nullptr
      );
    if (pSpacePoints && !pSpacePoints->isValid()) {
      throw art::Exception(art::errors::ProductNotFound)
        << "No space point associated with '" << fTrackModuleLabel.encode()
        << "' tracks.\n";
    }
    
    std::unique_ptr<art::FindManyP<recob::PFParticle>> pPFParticles(
      fPrintNParticles?
      new art::FindManyP<recob::PFParticle>(Tracks, evt, fTrackModuleLabel):
      nullptr
      );
    if (pPFParticles && !pPFParticles->isValid()) {
      throw art::Exception(art::errors::ProductNotFound)
        << "No particle-flow particle associated with '"
        << fTrackModuleLabel.encode() << "' tracks.\n";
    }
    
    for (unsigned int iTrack = 0; iTrack < Tracks->size(); ++iTrack) {
      const recob::Track& track = Tracks->at(iTrack);
      
      // print track information
      DumpTrack(iTrack, track);
      
      mf::LogVerbatim log(fOutputCategory);
      if (pHits || pSpacePoints || pPFParticles) {
        log << "\n  associated with:";
        if (pHits)
          log << " " << pHits->at(iTrack).size() << " hits;";
        if (pSpacePoints)
          log << " " << pSpacePoints->at(iTrack).size() << " space points;";
        if (pPFParticles)
          log << " " << pPFParticles->at(iTrack).size() << " PF particles;";
      } // if we have any association
      
      if (pHits && fPrintHits) {
        const auto& Hits = pHits->at(iTrack);
        log << "\n  hit indices (" << Hits.size() << "):\n";
        PrintAssociatedIndexTable(log, Hits, 10 /* 10 hits per line */, "    ");
      } // if print individual hits
      
      if (pSpacePoints && fPrintSpacePoints) {
        const auto& SpacePoints = pSpacePoints->at(iTrack);
        log << "\n  space point IDs (" << SpacePoints.size() << "):\n";
        PrintAssociatedIDTable
          (log, SpacePoints, 10 /* 10 hits per line */, "    ");
      } // if print individual space points
      
      if (pPFParticles && fPrintParticles) {
        const auto& PFParticles = pPFParticles->at(iTrack);
        log << "\n  particle indices (" << PFParticles.size() << "):\n";
        // currently a particle has no ID
        PrintAssociatedIndexTable
          (log, PFParticles, 10 /* 10 hits per line */, "    ");
      } // if print individual particles
    } // for tracks
  } // DumpTracks::analyze()
  
  
  //---------------------------------------------------------------------------
  void DumpTracks::DumpTrack
    (unsigned int iTrack, recob::Track const& track) const
  {
    // print a header for the track
    const unsigned int nPoints = track.NumberTrajectoryPoints();
    mf::LogVerbatim log(fOutputCategory);
    log
      << "Track #" << iTrack << " ID: " << track.ID()
        << std::fixed << std::setprecision(3)
        << " theta: " << track.Theta() << " rad, phi: " << track.Phi()
        << " rad, length: " << track.Length() << " cm"
      << "\n  start at: ( " << track.Vertex().X()
        << " ; " << track.Vertex().Y()
        << " ; " << track.Vertex().Z()
        << " ), direction: ( " << track.VertexDirection().X()
        << " ; " << track.VertexDirection().Y()
        << " ; " << track.VertexDirection().Z() << " )"
      << "\n  end at:   ( " << track.End().X()
        << " ; " << track.End().Y()
        << " ; " << track.End().Z()
        << " ), direction: ( " << track.EndDirection().X()
        << " ; " << track.EndDirection().Y()
        << " ; " << track.EndDirection().Z()
        << " )"
      << "\n  with "
        << nPoints << " trajectory points, "
        << track.NumberCovariance() << " covariance matrices";
    unsigned int nViews = 0;
    std::ostringstream sstr;
    if (PrintdQdXinView(sstr, track, geo::kU, "U")) ++nViews;
    if (PrintdQdXinView(sstr, track, geo::kV, "V")) ++nViews;
    if (PrintdQdXinView(sstr, track, geo::kZ, "Z")) ++nViews;
    if (PrintdQdXinView(sstr, track, geo::k3D, "3D")) ++nViews;
    if (nViews)
      log << ", dQ/dx in " << nViews << " views: " << sstr.str();
    else
      log << ", no dQ/dx";
    
    if (fPrintWayPoints > 0) {
      // print up to 10 (actually, 8 or 9) way points
      log << "\n  passes through:";
      unsigned int skip = std::max(nPoints / fPrintWayPoints, 1U);
      unsigned int iPoint = 0;
      while ((iPoint += skip) < nPoints) {
        const TVector3& point = track.LocationAtPoint(iPoint);
        log << "\n    [#" << iPoint << "] ("
          << point.X() << ", " << point.Y() << ", " << point.Z()
          << ")";
      } // while (iPoint)
    } // if print way points
  } // DumpTracks::DumpTrack()
  
  
  //---------------------------------------------------------------------------
  DEFINE_ART_MODULE(DumpTracks)

} // namespace recob