DumpPFParticles_module.cc

Go to the documentation of this file.

// LArSoft includes
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardata/ArtDataHelper/Dumpers/hexfloat.h"

// art libraries
#include "canvas/Utilities/InputTag.h"
#include "canvas/Persistency/Provenance/EventID.h"
#include "art/Framework/Principal/Provenance.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/Event.h"

// support libraries
#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/OptionalAtom.h"

// C//C++ standard libraries
#include <string>

// ... and more in the implementation part

namespace recob {
  
  class DumpPFParticles: public art::EDAnalyzer {
      public:
    
    struct Config {
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;
      
      fhicl::Atom<art::InputTag> PFModuleLabel {
        Name("PFModuleLabel"),
        Comment("label of producer of the recob::PFParticle to be dumped")
      };
      
      fhicl::Atom<std::string> OutputCategory {
        Name("OutputCategory"),
        Comment("message facility category used for output (for filtering)"),
        "DumpPFParticles"
      };
      
      fhicl::Atom<bool> PrintHexFloats {
        Name("PrintHexFloats"),
        Comment("print all the floating point numbers in base 16"),
        false
      };
      
      fhicl::OptionalAtom<unsigned int> MaxDepth {
        Name("MaxDepth"),
        Comment("at most this number of particle generations will be printed")
      };
      
      fhicl::Atom<bool> MakeParticleGraphs {
        Name("MakeParticleGraphs"),
        Comment("creates a DOT file with particle information for each event"),
        false
      };
      
    }; // struct Config
    
    using Parameters = art::EDAnalyzer::Table<Config>;
    
    explicit DumpPFParticles(Parameters const& config); 
    
    virtual void analyze (const art::Event& evt) override;

      private:

    art::InputTag fInputTag; 
    std::string fOutputCategory; 
    bool fPrintHexFloats; 
    unsigned int fMaxDepth; 
    bool fMakeEventGraphs; 
    
    
    static std::string DotFileName
      (art::EventID const& evtID, art::Provenance const& prodInfo);
    
    void MakePFParticleGraph(
      art::Event const& event,
      art::ValidHandle<std::vector<recob::PFParticle>> const& handle
      ) const;

  }; // class DumpPFParticles
  
} // namespace recob


//==============================================================================
//===  Implementation section
//==============================================================================
// LArSoft includes
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Seed.h"
#include "lardataobj/RecoBase/PCAxis.h"

// art libraries
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindOne.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "art/Framework/Principal/Handle.h"

// support libraries
#include "messagefacility/MessageLogger/MessageLogger.h"

// C//C++ standard libraries
#include <vector>
#include <fstream>
#include <utility> // std::swap()
#include <algorithm> // std::count(), std::find()
#include <limits> // std::numeric_limits<>


namespace {
  template <typename T>
  class int_map: private std::vector<T> {
    using map_type = std::vector<T>;
      public:
    using this_type = int_map<T>;
    
    using typename map_type::value_type;
    using typename map_type::reference;
    using typename map_type::const_reference;
    using typename map_type::iterator;
    using typename map_type::const_iterator;
    
    using typename map_type::size_type;
    
    
    int_map() = default;
    
    int_map(value_type invalid_value): invalid(invalid_value) {}
    
    
    void resize(size_t new_size)
      { get_map().resize(new_size, invalid_value()); }
    
    reference operator[] (size_t pos)
      { if (pos >= size()) resize(pos + 1); return get_map()[pos]; }
    
    const_reference operator[] (size_t pos) const
      { return (pos >= size())? invalid: get_map()[pos]; }
    
    using map_type::size;
    using map_type::empty;
    
    using map_type::begin;
    using map_type::end;
    using map_type::cbegin;
    using map_type::cend;
    using map_type::rbegin;
    using map_type::rend;
    using map_type::crbegin;
    using map_type::crend;
    
    void swap(this_type& data)
      { swap(data.get_map()); std::swap(data.invalid, invalid); }
    
    void swap(std::vector<T>& data) { map_type::swap(data); }
    
    // non-standard calls follow
    size_type n_invalid() const
      { return std::count(begin(), end(), invalid); }
    
    size_type n_valid() const { return size() - n_invalid(); }
    
    bool is_valid(size_type pos) const
      { return (pos < size())? is_valid_value(get_map()[pos]): false; }
    
    bool is_valid_value(value_type const& v) const { return v != invalid; }
    
    value_type const& invalid_value() const { return invalid; }
    
    
      private:
    value_type const invalid; 
    
    map_type& get_map() { return static_cast<map_type&>(*this); }
    map_type const& get_map() const
      { return static_cast<map_type const&>(*this); }
    
  }; // int_map<>
  
  
  //----------------------------------------------------------------------------
  int_map<size_t> CreateMap(std::vector<recob::PFParticle> const& particles) {
    
    int_map<size_t> pmap(std::numeric_limits<size_t>::max());
    pmap.resize(particles.size());
    
    size_t const nParticles = particles.size();
    for (size_t iPart = 0; iPart < nParticles; ++iPart)
      pmap[particles[iPart].Self()] = iPart;
    
    return pmap;
  } // CreateMap()
  
  
  //----------------------------------------------------------------------------
  bool hasDaughter(recob::PFParticle const& particle, size_t partID) {
    auto const& daughters = particle.Daughters();
    return daughters.end()
      != std::find(daughters.begin(), daughters.end(), partID);
  } // hasDaughter()
  
  
  //----------------------------------------------------------------------------
  class ParticleDumper {
      public:
    
    struct PrintOptions_t {
      bool hexFloats = false; 
      unsigned int maxDepth = std::numeric_limits<unsigned int>::max();
      std::string streamName;
    }; // PrintOptions_t
    
    
    ParticleDumper(std::vector<recob::PFParticle> const& particle_list)
      : ParticleDumper(particle_list, {})
      {}
    
    ParticleDumper(
      std::vector<recob::PFParticle> const& particle_list,
      PrintOptions_t print_options
      )
      : particles(particle_list)
      , options(print_options)
      , visited(particles.size(), 0U)
      , particle_map (CreateMap(particles))
      {}
    
    
    void SetVertices(art::FindOne<recob::Vertex> const* vtx_query)
      { vertices = vtx_query; }
    
    void SetTracks(art::FindMany<recob::Track> const* trk_query)
      { tracks = trk_query; }
    
    void SetClusters(art::FindMany<recob::Cluster> const* cls_query)
      { clusters = cls_query; }
    
    void SetSeeds(art::FindMany<recob::Seed> const* seed_query)
      { seeds = seed_query; }
    
    void SetSpacePoints(art::FindMany<recob::SpacePoint> const* sp_query)
      { spacepoints = sp_query; }
    
    void SetPCAxes(art::FindMany<recob::PCAxis> const* pca_query)
      { pcaxes = pca_query; }
    
    
    template <typename Stream>
    void DumpParticle(
      Stream&& out, size_t iPart, std::string indentstr = "",
      unsigned int gen = 0
      ) const;
    
    
    template <typename Stream>
    void DumpParticleWithID(
      Stream&& out, size_t pID, std::string indentstr = "",
      unsigned int gen = 0
      ) const;
    
    
    void DumpAllPrimaries(std::string indentstr = "") const;
    
    
    void DumpAllParticles(std::string indentstr = "") const;
    
    
    template <typename Stream>
    static void DumpPDGID(Stream&& out, int ID);
    
    
      protected:
    std::vector<recob::PFParticle> const& particles; 
    
    PrintOptions_t options; 
    
    art::FindOne<recob::Vertex> const* vertices = nullptr;
    
    art::FindMany<recob::Track> const* tracks = nullptr;
    
    art::FindMany<recob::Cluster> const* clusters = nullptr;
    
    art::FindMany<recob::Seed> const* seeds = nullptr;
    
    art::FindMany<recob::SpacePoint> const* spacepoints = nullptr;
    
    art::FindMany<recob::PCAxis> const* pcaxes = nullptr;
    
    mutable std::vector<unsigned int> visited;
    
    int_map<size_t> const particle_map; 
    
    
    template <typename Stream>
    void DumpPFParticleInfo(
      Stream&& out,
      recob::PFParticle const& part,
      unsigned int iPart,
      std::string indentstr
      ) const;
    
    template <typename Stream>
    void DumpVertex(
      Stream&& out,
      cet::maybe_ref<recob::Vertex const> VertexRef,
      std::string indentstr
      ) const;
    
    template <typename Stream>
    void DumpPCAxisDirection
      (Stream&& out, recob::PCAxis const& axis) const;
    
    template <typename Stream>
    void DumpPCAxes(
      Stream&& out,
      std::vector<recob::PCAxis const*> const& myAxes,
      std::string indentstr
      ) const;
    
    template <typename Stream>
    void DumpTrack(Stream&& out, recob::Track const& track) const;
    
    template <typename Stream>
    void DumpTracks(
      Stream&& out,
      std::vector<recob::Track const*> const& myTracks,
      std::string indentstr
      ) const;

    template <typename Stream>
    void DumpSeed
      (Stream&& out, recob::Seed const& seed, std::string indentstr) const;
    
    template <typename Stream>
    void DumpSeeds(
      Stream&& out,
      std::vector<recob::Seed const*> const& mySeeds,
      std::string indentstr
      ) const;
    
    template <typename Stream>
    void DumpSpacePoint(Stream&& out, recob::SpacePoint const& sp) const;
      
    template <typename Stream>
    void DumpSpacePoints(
      Stream&& out,
      std::vector<recob::SpacePoint const*> const& mySpacePoints,
      std::string indentstr
      ) const;
      
    template <typename Stream>
    void DumpClusterSummary(Stream&& out, recob::Cluster const& track) const;
    
    template <typename Stream>
    void DumpClusters(
      Stream&& out,
      std::vector<recob::Cluster const*> const& myClusters,
      std::string indentstr
      ) const;
    
  }; // class ParticleDumper
  
  
  //----------------------------------------------------------------------------
  class PFParticleGraphMaker {
      public:
    
    PFParticleGraphMaker() = default;
    
    template <typename Stream>
    void MakeGraph
      (Stream&& out, std::vector<recob::PFParticle> const& particles) const;
    
    template <typename Stream>
    void WriteGraphHeader(Stream&& out) const;
    
    template <typename Stream>
    void WriteParticleRelations
      (Stream&& out, std::vector<recob::PFParticle> const& particles) const;
    
    template <typename Stream>
    void WriteParticleNodes
      (Stream&& out, std::vector<recob::PFParticle> const& particles) const;
    
    template <typename Stream>
    void WriteParticleEdges
      (Stream&& out, std::vector<recob::PFParticle> const& particles) const;
    
    template <typename Stream>
    void WriteGraphFooter(Stream&& out) const;
    
  }; // class PFParticleGraphMaker
  
  
  //----------------------------------------------------------------------------
  //--- ParticleDumper implementation
  //--- 
  template <typename Stream>
  void ParticleDumper::DumpParticle(
    Stream&& out, size_t iPart, std::string indentstr /* = "" */,
    unsigned int gen /* = 0 */
    ) const
  {
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    
    recob::PFParticle const& part = particles.at(iPart);
    ++visited[iPart];
    
    if (visited[iPart] > 1) {
      out << indentstr << "particle " << part.Self()
        << " already printed!!!";
      return;
    }
    
    //
    // intro
    //
    DumpPFParticleInfo(std::forward<Stream>(out), part, iPart, indentstr);
    
    //
    // vertex information
    //
    if (vertices)
      DumpVertex(std::forward<Stream>(out), vertices->at(iPart), indentstr);
    
    // daughters tag
    if (part.NumDaughters() > 0)
      out << " with " << part.NumDaughters() << " daughters";
    else
      out << " with no daughter";
    
    //
    // axis
    //
    if (pcaxes)
      DumpPCAxes(std::forward<Stream>(out), pcaxes->at(iPart), indentstr);
    
    //
    // tracks
    //
    if (tracks)
      DumpTracks(std::forward<Stream>(out), tracks->at(iPart), indentstr);
    
    //
    // seeds
    //
    if (seeds)
      DumpSeeds(std::forward<Stream>(out), seeds->at(iPart), indentstr);
    
    //
    // space points
    //
    if (spacepoints) {
      DumpSpacePoints
        (std::forward<Stream>(out), spacepoints->at(iPart), indentstr);
    }
    
    //
    // clusters
    //
    if (clusters) {
      DumpClusters
        (std::forward<Stream>(out), clusters->at(iPart), indentstr);
    }
    
    //
    // daughters
    //
    auto const PartID = part.Self();
    if (part.NumDaughters() > 0) {
      std::vector<size_t> const& daughters = part.Daughters();
      out << "\n" << indentstr
        << "  " << part.NumDaughters() << " particle daughters";
      if (gen > 0) {
        out << ":";
        for (size_t DaughterID: daughters) {
          if (DaughterID == PartID) {
            out << "\n" << indentstr
              << "    oh, just great: this particle is its own daughter!";
          }
          else {
            out << '\n';
            DumpParticleWithID(out, DaughterID, indentstr + "  ", gen - 1);
          }
        }
      } // if descending
      else {
        out << " (further descend suppressed)";
      }
    } // if has daughters
    
    //
    // warnings
    //
    if (visited[iPart] == 2) {
      out << "\n" << indentstr << "  WARNING: particle ID=" << PartID
        << " connected more than once!";
    }
    
    //
    // done
    //
    
  } // ParticleDumper::DumpParticle()
  
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpParticleWithID(
    Stream&& out, size_t pID, std::string indentstr /* = "" */,
    unsigned int gen /* = 0 */
  ) const {
    size_t const pos = particle_map[pID];
    if (particle_map.is_valid_value(pos)) {
      DumpParticle(out, pos, indentstr, gen);
    }
    else {
      out /* << "\n" */ << indentstr << "<ID=" << pID << " not found>";
    }
  } // ParticleDumper::DumpParticleWithID()
    
    
  //----------------------------------------------------------------------------
  void ParticleDumper::DumpAllPrimaries(std::string indentstr /* = "" */) const
  {
    indentstr += "  ";
    size_t const nParticles = particles.size();
    unsigned int nPrimaries = 0;
    for (size_t iPart = 0; iPart < nParticles; ++iPart) {
      if (!particles[iPart].IsPrimary()) continue;
      DumpParticle(
        mf::LogVerbatim(options.streamName),
        iPart, indentstr, options.maxDepth
        );
    } // for
    if (nPrimaries == 0) {
      mf::LogVerbatim(options.streamName)
        << indentstr << "No primary particle found";
    }
  } // ParticleDumper::DumpAllPrimaries()
    
    
  //----------------------------------------------------------------------------
  void ParticleDumper::DumpAllParticles(std::string indentstr /* = "" */) const
  {
    // first print all the primary particles
    DumpAllPrimaries(indentstr);
    // then find out if there are any that are "disconnected"
    unsigned int const nDisconnected
      = std::count(visited.begin(), visited.end(), 0U);
    if (nDisconnected) {
      mf::LogVerbatim(options.streamName) << indentstr
        << nDisconnected << " particles not coming from primary ones:";
      size_t const nParticles = visited.size();
      for (size_t iPart = 0; iPart < nParticles; ++iPart) {
        if (visited[iPart] > 0) continue;
        DumpParticle(
          mf::LogVerbatim(options.streamName), iPart, indentstr + "  ",
          options.maxDepth
          );
      } // for unvisited particles
      mf::LogVerbatim(options.streamName) << indentstr
        << "(end of " << nDisconnected << " particles not from primaries)";
    } // if there are disconnected particles
    // TODO finally, note if there are multiply-connected particles
    
  } // ParticleDumper::DumpAllParticles()
    
    
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpPDGID(Stream&& out, int ID) {
    switch (ID) {
      case -11: out << "e+";  break;
      case  11: out << "e-";  break;
      case -13: out << "mu+"; break;
      case  13: out << "mu-"; break;
      default:  out << "MCID=" << ID; break;
    } // switch
  } // DumpPDGID()
    
    
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpPFParticleInfo(
    Stream&& out,
    recob::PFParticle const& part,
    unsigned int iPart,
    std::string indentstr
    ) const
  {
    auto const PartID = part.Self();
    out << indentstr << "ID=" << PartID;
    if (iPart != PartID) out << " [#" << iPart << "]";
    out << " (type: ";
    DumpPDGID(out, part.PdgCode());
    out << ")";
    if (part.IsPrimary()) out << " (primary)";
    else                  out << " from ID=" << part.Parent();
  } // ParticleDumper::DumpPFParticleInfo()
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpVertex(
    Stream&& out,
    cet::maybe_ref<recob::Vertex const> VertexRef,
    std::string indentstr
    ) const
  {
    if (!VertexRef.isValid()) {
      out << " [no vertex found!]";
      return;
    }
    recob::Vertex const& vertex = VertexRef.ref();
    std::array<double, 3> vtx_pos;
    vertex.XYZ(vtx_pos.data());
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    out << " [decay at ("
      << hexfloat(vtx_pos[0]) << "," << hexfloat(vtx_pos[1])
      << "," << hexfloat(vtx_pos[2])
      << "), ID=" << vertex.ID() << "]";
    
  } // ParticleDumper::DumpVertex()
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpPCAxisDirection
    (Stream&& out, recob::PCAxis const& axis) const
  {
    if (!axis.getSvdOK()) {
      out << "axis is invalid";
      return;
    }
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    out << "axis ID=" << axis.getID() << ", principal: ("
        << hexfloat(axis.getEigenVectors()[0][0])
        << ", " << hexfloat(axis.getEigenVectors()[0][1])
        << ", " << hexfloat(axis.getEigenVectors()[0][2]) << ")";
  } // ParticleDumper::DumpPCAxisDirection()
  
  template <typename Stream>
  void ParticleDumper::DumpPCAxes(
    Stream&& out,
    std::vector<recob::PCAxis const*> const& myAxes,
    std::string indentstr
    ) const
  {
    out << "\n" << indentstr;
    if (myAxes.empty()) {
      out << " [no axis found!]";
      return;
    }
    if (myAxes.size() == 1) {
      DumpPCAxisDirection(std::forward<Stream>(out), *(myAxes.front()));
    }
    else {
      out << "  " << myAxes.size() << " axes present:";
      for (recob::PCAxis const* axis: myAxes) {
        out << "\n" << indentstr << "    ";
        DumpPCAxisDirection(std::forward<Stream>(out), *axis);
      } // for axes
    } // else
  } // ParticleDumper::DumpPCAxes()
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpSeed
    (Stream&& out, recob::Seed const& seed, std::string indentstr) const
  {
    if (!seed.IsValid()) {
      out << "  <invalid>";
      return;
    }
    std::array<double, 3> start, dir;
    seed.GetDirection(dir.data(), nullptr);
    seed.GetPoint(start.data(), nullptr);
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    out << "\n" << indentstr
      << "    (" << hexfloat(start[0]) << "," << hexfloat(start[1])
      << "," << hexfloat(start[2])
      << ")=>(" << hexfloat(dir[0]) << "," << hexfloat(dir[1])
      << "," << hexfloat(dir[2])
      << "), " << hexfloat(seed.GetLength()) << " cm"
      ;
  } // ParticleDumper::DumpSeed()
  
  template <typename Stream>
  void ParticleDumper::DumpSeeds(
    Stream&& out,
    std::vector<recob::Seed const*> const& mySeeds,
    std::string indentstr
    ) const
  {
    if (mySeeds.empty()) return;
    out << "\n" << indentstr << "  "
      << mySeeds.size() << " seeds:";
    for (recob::Seed const* seed: mySeeds) 
      DumpSeed(std::forward<Stream>(out), *seed, indentstr);
  } // ParticleDumper::DumpSeeds()
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpSpacePoint
    (Stream&& out, recob::SpacePoint const& sp) const
  {
    const double* pos = sp.XYZ();
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    out << "  ID=" << sp.ID()
      << " at (" << hexfloat(pos[0]) << "," << hexfloat(pos[1])
      << "," << hexfloat(pos[2]) << ") cm"
      ;
  } // ParticleDumper::DumpSpacePoints()
    
  template <typename Stream>
  void ParticleDumper::DumpSpacePoints(
    Stream&& out,
    std::vector<recob::SpacePoint const*> const& mySpacePoints,
    std::string indentstr
    ) const
  {
    out << "\n" << indentstr << "  ";
    if (mySpacePoints.empty()) {
      out << "no space points";
      return;
    }
    constexpr unsigned int PointsPerLine = 2;
    out << mySpacePoints.size() << " space points:";
    for (size_t iSP = 0; iSP < mySpacePoints.size(); ++iSP) {
      if ((iSP % PointsPerLine) == 0)
        out << "\n" << indentstr << "  ";
      
      DumpSpacePoint(std::forward<Stream>(out), *(mySpacePoints[iSP]));
    } // for points
  } // ParticleDumper::DumpSpacePoints()
    
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpTrack(Stream&& out, recob::Track const& track) const
  {
    lar::OptionalHexFloat hexfloat(options.hexFloats);
    out << " length " << hexfloat(track.Length())
      << "cm from (" << hexfloat(track.Vertex().X())
        << ";" << hexfloat(track.Vertex().Y())
        << ";" << hexfloat(track.Vertex().Z())
      << ") to (" << hexfloat(track.End().X())
        << ";" << hexfloat(track.End().Y())
        << ";" << hexfloat(track.End().Z())
      << ") (ID=" << track.ID() << ")";
  } // ParticleDumper::DumpTrack()
    
  template <typename Stream>
  void ParticleDumper::DumpTracks(
    Stream&& out,
    std::vector<recob::Track const*> const& myTracks,
    std::string indentstr
    ) const
  {
    if (myTracks.empty()) return;
    out << "\n" << indentstr << "  "
      << myTracks.size() << " tracks:";
    for (recob::Track const* track: myTracks) {
      if (myTracks.size() > 1) out << "\n" << indentstr << "   ";
      DumpTrack(std::forward<Stream>(out), *track);
    } // for tracks
  } // ParticleDumper::DumpTracks()
    
  //----------------------------------------------------------------------------
  template <typename Stream>
  void ParticleDumper::DumpClusterSummary
    (Stream&& out, recob::Cluster const& cluster) const
  {
    out << "  " << cluster.NHits() << " hits on " << cluster.Plane()
        << " (ID=" << cluster.ID() << ")";
  } // ParticleDumper::DumpClusterSummary()
  
  template <typename Stream>
  void ParticleDumper::DumpClusters(
    Stream&& out,
    std::vector<recob::Cluster const*> const& myClusters,
    std::string indentstr
    ) const
  {
    if (myClusters.empty()) return;
    out << "\n" << indentstr << "  "
      << myClusters.size() << " clusters:";
    for (recob::Cluster const* cluster: myClusters) {
      if (myClusters.size() > 1) out << "\n" << indentstr << "   ";
      DumpClusterSummary(std::forward<Stream>(out), *cluster);
    }
  } // ParticleDumper::DumpClusters()
  
  
  //----------------------------------------------------------------------------
  //---  PFParticleGraphMaker
  //---
  template <typename Stream>
  void PFParticleGraphMaker::MakeGraph
    (Stream&& out, std::vector<recob::PFParticle> const& particles) const
  {
    
    WriteGraphHeader      (std::forward<Stream>(out));
    WriteParticleRelations(std::forward<Stream>(out), particles);
    WriteGraphFooter      (std::forward<Stream>(out));
    
  } // PFParticleGraphMaker::MakeGraph()
  
  
  //----------------------------------------------------------------------------
  template <typename Stream>
  void PFParticleGraphMaker::WriteGraphHeader(Stream&& out) const {
    out
      << "\ndigraph {"
      << "\n";
  } // PFParticleGraphMaker::WriteGraphHeader()
  
  
  template <typename Stream>
  void PFParticleGraphMaker::WriteParticleNodes
    (Stream&& out, std::vector<recob::PFParticle> const& particles) const
  {
    
    for (auto const& particle: particles) {
      
      std::ostringstream label;
      label << "#" << particle.Self();
      if (particle.PdgCode() != 0) {
        label << ", ";
        ParticleDumper::DumpPDGID(label, particle.PdgCode());
      }
      
      out << "\n  P" << particle.Self()
        << " ["
        << " label = \"" << label.str() << "\"";
      if (particle.IsPrimary()) out << ", style = bold";
      out << " ]";
    } // for
    
  } // PFParticleGraphMaker::WriteParticleNodes()
  
  
  template <typename Stream>
  void PFParticleGraphMaker::WriteParticleEdges
    (Stream&& out, std::vector<recob::PFParticle> const& particles) const
  {
    auto particle_map = CreateMap(particles);
    
    out
      << "\n  "
      << "\n  // relations"
      << "\n  // "
      << "\n  // the arrow is close to the information provider,;"
      << "\n  // and it points from parent to daughter"
      << "\n  // "
      << "\n  // "
      << "\n  "
      ;
    
    for (auto const& particle: particles) {
      auto const partID = particle.Self();
      
      // draw parent line
      if (!particle.IsPrimary()) {
        auto const parentID = particle.Parent();
        size_t const iParent = particle_map[parentID];
        if (!particle_map.is_valid_value(iParent)) {
          // parent is ghost
          out
            << "\nP" << parentID
              << " [ style = dashed, color = red"
              << ", label = \"(#" << parentID << ")\" ] // ghost particle"
            << "\nP" << parentID << " -> P" << partID
              << " [ dir = both, arrowhead = empty, arrowtail = none ]";
        }
        else {
          // parent is a known particle
          recob::PFParticle const& parent = particles[iParent];
          
          // is the relation bidirectional?
          if (hasDaughter(parent, partID)) {
            out
              << "\nP" << parentID << " -> P" << partID
              << " [ dir = both, arrowtail = inv ]";
          } // if bidirectional
          else {
            out
              << "\nP" << parentID << " -> P" << partID
                << " [ dir = forward, color = red ]"
                << " // claimed parent";
          } // if parent disowns
          
        } // if ghost ... else
      } // if not primary
      
      // print daughter relationship only if daughters do not recognise us
      for (auto daughterID: particle.Daughters()) {
        size_t const iDaughter = particle_map[daughterID];
        if (!particle_map.is_valid_value(iDaughter)) {
          // daughter is ghost
          out
            << "\nP" << daughterID
              << " [ style = dashed, color = red"
              << ", label = \"(#" << daughterID << ")\" ] // ghost daughter"
            << "\nP" << partID << " -> P" << daughterID
              << " [ dir = both, arrowhead = none, arrowtail = invempty ]";
        }
        else {
          // daughter is a known particle
          recob::PFParticle const& daughter = particles[iDaughter];
          
          // is the relation bidirectional? (if so, the daughter will draw)
          if (daughter.Parent() != partID) {
            out
              << "\nP" << partID << " -> P" << daughterID
              << " [ dir = both, arrowhead = none, arrowtail = inv, color = orange ]"
              << " // claimed daughter";
          } // if parent disowns
          
        } // if ghost ... else
        
        
      } // for all daughters
      
    } // for all particles
    
  } // PFParticleGraphMaker::WriteParticleNodes()
  
  
  template <typename Stream>
  void PFParticleGraphMaker::WriteParticleRelations
    (Stream&& out, std::vector<recob::PFParticle> const& particles) const
  {
    out << "\n  // " << particles.size() << " particles (nodes)";
    WriteParticleNodes(std::forward<Stream>(out), particles);
    
    out
      << "\n"
      << "\n  // parent/children relations";
    WriteParticleEdges(std::forward<Stream>(out), particles);
    
  } // PFParticleGraphMaker::WriteParticleRelations()
  
  
  template <typename Stream>
  void PFParticleGraphMaker::WriteGraphFooter(Stream&& out) const {
    out
      << "\n"
      << "\n} // digraph"
      << std::endl;
  } // PFParticleGraphMaker::WriteGraphFooter()
    
  
  
  //----------------------------------------------------------------------------
  
} // local namespace



namespace recob {
  
  //----------------------------------------------------------------------------
  DumpPFParticles::DumpPFParticles(Parameters const& config) 
    : EDAnalyzer(config)
    , fInputTag(config().PFModuleLabel())
    , fOutputCategory(config().OutputCategory())
    , fPrintHexFloats(config().PrintHexFloats())
    , fMaxDepth(std::numeric_limits<unsigned int>::max())
    , fMakeEventGraphs(config().MakeParticleGraphs())
    {
      // here we are handling the optional configuration key as it had just a
      // default value
      if (!config().MaxDepth(fMaxDepth))
        fMaxDepth = std::numeric_limits<unsigned int>::max();
    }
  
  
  //----------------------------------------------------------------------------
  std::string DumpPFParticles::DotFileName
    (art::EventID const& evtID, art::Provenance const& prodInfo)
  {
    return prodInfo.processName()
      + '_' + prodInfo.moduleLabel()
      + '_' + prodInfo.productInstanceName()
      + "_Run" + std::to_string(evtID.run())
      + "_Subrun" + std::to_string(evtID.subRun())
      + "_Event" + std::to_string(evtID.event())
      + "_particles.dot";
  } // DumpPFParticles::DotFileName()
  
  
  //----------------------------------------------------------------------------
  void DumpPFParticles::MakePFParticleGraph(
    art::Event const& event,
    art::ValidHandle<std::vector<recob::PFParticle>> const& handle
  ) const {
    art::EventID const eventID = event.id();
    std::string fileName = DotFileName(eventID, *(handle.provenance()));
    std::ofstream outFile(fileName); // overwrite by default
    
    outFile
      <<   "// " << fileName
      << "\n// "
      << "\n// Created for run " << eventID.run()
               << " subrun " << eventID.subRun()
               << " event " << eventID.event()
      << "\n// "
      << "\n// dump of " << handle->size() << " particles"
      << "\n// "
      << std::endl;
    
    PFParticleGraphMaker graphMaker;
    graphMaker.MakeGraph(outFile, *handle);
    
  } // DumpPFParticles::MakePFParticleGraph()
  
  
  //----------------------------------------------------------------------------
  void DumpPFParticles::analyze(const art::Event& evt) {
    
    //
    // collect all the available information
    //
    // fetch the data to be dumped on screen
    art::ValidHandle<std::vector<recob::PFParticle>> PFParticles
      = evt.getValidHandle<std::vector<recob::PFParticle>>(fInputTag);
    
    if (fMakeEventGraphs)
      MakePFParticleGraph(evt, PFParticles);
    
    art::FindOne<recob::Vertex> const ParticleVertices
      (PFParticles, evt, fInputTag);
    art::FindMany<recob::Track> const ParticleTracks
      (PFParticles, evt, fInputTag);
    art::FindMany<recob::Cluster> const ParticleClusters
      (PFParticles, evt, fInputTag);
    art::FindMany<recob::Seed> const ParticleSeeds
      (PFParticles, evt, fInputTag);
    art::FindMany<recob::SpacePoint> const ParticleSpacePoints
      (PFParticles, evt, fInputTag);
    art::FindMany<recob::PCAxis> const ParticlePCAxes
      (PFParticles, evt, fInputTag);
    
    size_t const nParticles = PFParticles->size();
    mf::LogVerbatim(fOutputCategory) << "Event " << evt.id()
      << " contains " << nParticles << " particles from '"
      << fInputTag.encode() << "'";
    
    // prepare the dumper
    ParticleDumper::PrintOptions_t options;
    options.hexFloats = fPrintHexFloats;
    options.maxDepth = fMaxDepth;
    options.streamName = fOutputCategory;
    ParticleDumper dumper(*PFParticles, options);
    if (ParticleVertices.isValid()) dumper.SetVertices(&ParticleVertices);
    else mf::LogPrint("DumpPFParticles") << "WARNING: vertex information not available";
    if (ParticleTracks.isValid()) dumper.SetTracks(&ParticleTracks);
    else mf::LogPrint("DumpPFParticles") << "WARNING: track information not available";
    if (ParticleClusters.isValid()) dumper.SetClusters(&ParticleClusters);
    else mf::LogPrint("DumpPFParticles") << "WARNING: cluster information not available";
    if (ParticleSeeds.isValid()) dumper.SetSeeds(&ParticleSeeds);
    else mf::LogPrint("DumpPFParticles") << "WARNING: seed information not avaialble";
    if (ParticleSpacePoints.isValid())
      dumper.SetSpacePoints(&ParticleSpacePoints);
    else {
      mf::LogPrint("DumpPFParticles")
        << "WARNING: space point information not available";
    }
    if (ParticlePCAxes.isValid())
      dumper.SetPCAxes(&ParticlePCAxes);
    else {
      mf::LogPrint("DumpPFParticles")
        << "WARNING: principal component axis not available";
    }
    dumper.DumpAllParticles("  ");
    
    mf::LogVerbatim(fOutputCategory) << "\n"; // two empty lines
    
  } // DumpPFParticles::analyze()

  DEFINE_ART_MODULE(DumpPFParticles)

} // namespace recob