ParticleDecayId_module.cc

Go to the documentation of this file.

// Class:       ParticleDecayId
// Module Type: producer
// File:        ParticleDecayId_module.cc
// Authors:     dorota.stefan@cern.ch pplonski86@gmail.com robert.sulej@cern.ch
//
// THIS IS STIIL DEVELOPMENT NOW - CODE MAKING VERTICES MAY STILL CHANGE STRATEGY
//

#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/Wire.h"
#include "larrecodnn/ImagePatternAlgs/Tensorflow/PointIdAlg/PointIdAlg.h"

#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/Principal/Run.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Table.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include <TVector3.h>

#include <iostream>
#include <map>
#include <memory>
#include <utility>
#include <vector>

namespace nnet {

  class ParticleDecayId : public art::EDProducer {
  public:
    struct Config {
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;

      fhicl::Table<nnet::PointIdAlg::Config> PointIdAlg{Name("PointIdAlg")};

      fhicl::Atom<art::InputTag> WireLabel{
        Name("WireLabel"),
        Comment("tag of deconvoluted ADC on wires (recob::Wire)")};

      fhicl::Atom<art::InputTag> TrackModuleLabel{
        Name("TrackModuleLabel"),
        Comment("tag of tracks where decays points are to be found")};

      fhicl::Atom<double> RoiThreshold{
        Name("RoiThreshold"),
        Comment("search for decay points where the net output > ROI threshold")};

      fhicl::Atom<double> PointThreshold{
        Name("PointThreshold"),
        Comment("tag decay point if it is detected in at least two planes with net outputs product "
                "> POINT threshold")};

      fhicl::Atom<int> SkipView{Name("SkipView"),
                                Comment("use all views to find decays if -1, or skip the view with "
                                        "provided index and use only the two other views")};
    };
    using Parameters = art::EDProducer::Table<Config>;
    explicit ParticleDecayId(Parameters const& p);

    ParticleDecayId(ParticleDecayId const&) = delete;
    ParticleDecayId(ParticleDecayId&&) = delete;
    ParticleDecayId& operator=(ParticleDecayId const&) = delete;
    ParticleDecayId& operator=(ParticleDecayId&&) = delete;

  private:
    void produce(art::Event& e) override;

    bool DetectDecay(detinfo::DetectorClocksData const& clockData,
                     detinfo::DetectorPropertiesData const& detProp,
                     const std::vector<recob::Wire>& wires,
                     const std::vector<art::Ptr<recob::Hit>>& hits,
                     std::map<size_t, TVector3>& spoints,
                     std::vector<std::pair<TVector3, double>>& result);

    PointIdAlg fPointIdAlg;

    art::InputTag fWireProducerLabel;
    art::InputTag fTrackModuleLabel;

    double fRoiThreshold, fPointThreshold;

    int fSkipView;
  };
  // ------------------------------------------------------

  ParticleDecayId::ParticleDecayId(ParticleDecayId::Parameters const& config)
    : EDProducer{config}
    , fPointIdAlg(config().PointIdAlg())
    , fWireProducerLabel(config().WireLabel())
    , fTrackModuleLabel(config().TrackModuleLabel())
    , fRoiThreshold(config().RoiThreshold())
    , fPointThreshold(config().PointThreshold())
    , fSkipView(config().SkipView())
  {
    produces<std::vector<recob::Vertex>>();
    produces<art::Assns<recob::Vertex, recob::Track>>();
  }
  // ------------------------------------------------------

  void ParticleDecayId::produce(art::Event& evt)
  {
    std::cout << std::endl << "event " << evt.id().event() << std::endl;

    auto vtxs = std::make_unique<std::vector<recob::Vertex>>();
    auto vtx2trk = std::make_unique<art::Assns<recob::Vertex, recob::Track>>();

    auto wireHandle = evt.getValidHandle<std::vector<recob::Wire>>(fWireProducerLabel);
    auto trkListHandle = evt.getValidHandle<std::vector<recob::Track>>(fTrackModuleLabel);
    auto spListHandle = evt.getValidHandle<std::vector<recob::SpacePoint>>(fTrackModuleLabel);

    art::FindManyP<recob::Hit> hitsFromTracks(trkListHandle, evt, fTrackModuleLabel);
    art::FindManyP<recob::SpacePoint> spFromTracks(trkListHandle, evt, fTrackModuleLabel);
    art::FindManyP<recob::Hit> hitsFromSPoints(spListHandle, evt, fTrackModuleLabel);

    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);

    std::vector<std::pair<TVector3, double>> decays;
    for (size_t i = 0; i < hitsFromTracks.size(); ++i) {
      auto hits = hitsFromTracks.at(i);
      auto spoints = spFromTracks.at(i);
      if (hits.empty()) continue;

      std::map<size_t, TVector3> trkSpacePoints;
      for (const auto& p : spoints) {
        auto sp_hits = hitsFromSPoints.at(p.key());
        for (const auto& h : sp_hits) {
          trkSpacePoints[h.key()] = TVector3(p->XYZ()[0], p->XYZ()[1], p->XYZ()[2]);
        }
      }

      DetectDecay(clockData, detProp, *wireHandle, hits, trkSpacePoints, decays);
    }

    double xyz[3];
    for (const auto& p3d : decays) {

      xyz[0] = p3d.first.X();
      xyz[1] = p3d.first.Y();
      xyz[2] = p3d.first.Z();
      std::cout << "   detected: [" << xyz[0] << ", " << xyz[1] << ", " << xyz[2]
                << "] p:" << p3d.second << std::endl;

      size_t vidx = vtxs->size();
      vtxs->push_back(recob::Vertex(xyz, vidx));

      // to do: assn to eg. appropriate track
      // selected among set of connected tracks
      //art::Ptr<recob::Track> tptr(trkListHandle, i);
      //art::Ptr<recob::Vertex> vptr(vid, vidx, evt.productGetter(vid));
      //vtx2trk->addSingle(vptr, tptr);
    }

    evt.put(std::move(vtxs));
    evt.put(std::move(vtx2trk));
  }
  // ------------------------------------------------------

  bool ParticleDecayId::DetectDecay(detinfo::DetectorClocksData const& clockData,
                                    detinfo::DetectorPropertiesData const& detProp,
                                    const std::vector<recob::Wire>& wires,
                                    const std::vector<art::Ptr<recob::Hit>>& hits,
                                    std::map<size_t, TVector3>& spoints,
                                    std::vector<std::pair<TVector3, double>>& result)
  {
    const size_t nviews = 3;

    std::vector<art::Ptr<recob::Hit>> wire_drift[nviews];
    for (size_t i = 0; i < hits.size(); ++i) // split hits between views
    {
      wire_drift[hits[i]->View()].push_back(hits[i]);
    }

    std::vector<float> outputs[nviews];
    for (size_t v = 0; v < nviews;
         ++v) // calculate nn outputs for each view (hopefully not changing cryo/tpc many times)
    {
      outputs[v].resize(wire_drift[v].size(), 0);
      for (size_t i = 0; i < wire_drift[v].size(); ++i) {
        int tpc = wire_drift[v][i]->WireID().TPC;
        int cryo = wire_drift[v][i]->WireID().Cryostat;

        fPointIdAlg.setWireDriftData(clockData, detProp, wires, v, tpc, cryo);

        outputs[v][i] = fPointIdAlg.predictIdVector(wire_drift[v][i]->WireID().Wire,
                                                    wire_drift[v][i]->PeakTime())[0]; // p(decay)
      }
    }

    std::vector<std::pair<size_t, float>> candidates2d[nviews];
    std::vector<std::pair<TVector3, float>> candidates3d[nviews];
    for (size_t v = 0; v < nviews; ++v) {
      size_t idx = 0;
      while (idx < outputs[v].size()) {
        if (outputs[v][idx] > fRoiThreshold) {
          size_t ci = idx;
          float max = outputs[v][idx];
          ++idx;

          while ((idx < outputs[v].size()) && (outputs[v][idx] > fRoiThreshold)) {
            if (outputs[v][idx] > max) {
              max = outputs[v][idx];
              ci = idx;
            }
            ++idx;
          }
          candidates2d[v].emplace_back(ci, max);
          candidates3d[v].emplace_back(spoints[wire_drift[v][ci].key()], max);
        }
        else
          ++idx;
      }
    }

    double min_dist =
      2.0; // [cm], threshold for today to distinguish between two different candidates,
           // if belo threshold, then use 3D point corresponding to higher cnn output

    // need coincidence of high cnn out in two views, then look if there is another close candidate
    // and again select by cnn output value, would like to have few strong candidates
    bool found = false;
    for (size_t v = 0; v < nviews - 1; ++v) {
      for (size_t i = 0; i < candidates3d[v].size(); ++i) {
        TVector3 c0(candidates3d[v][i].first);
        float p0 = candidates3d[v][i].second;

        for (size_t u = v + 1; u < nviews; ++u) {
          for (size_t j = 0; j < candidates3d[v].size(); ++j) {
            TVector3 c1(candidates3d[v][j].first);
            float p1 = candidates3d[v][j].second;

            if ((c0 - c1).Mag() < min_dist) {
              TVector3 c(c0);
              if (p1 > p0) { c = c1; }
              double p = p0 * p1;

              if (p > fPointThreshold) {
                double d, dmin = min_dist;
                size_t kmin = 0;
                for (size_t k = 0; k < result.size(); ++k) {
                  d = (result[k].first - c).Mag();
                  if (d < dmin) {
                    dmin = d;
                    kmin = k;
                  }
                }
                if (dmin < min_dist) {
                  if (result[kmin].second < p) // replace previously found point
                  {
                    result[kmin].first = c;
                    result[kmin].second = p;
                    found = true;
                  }
                }
                else // nothing close in the list, add new point
                {
                  result.emplace_back(c, p);
                  found = true;
                }
              } // if (p > fPointThreshold)
            }   // coincidence: points from views u and v are close
          }     // loop over points in view u
        }       // loop over views u
      }         // loop over points in view v
    }           // loop over views v
    return found;
  }

  DEFINE_ART_MODULE(ParticleDecayId)

}