MVAReader.h

Go to the documentation of this file.

// \version
//
// \brief Wrappers for accessing MVA results and associated data products
//
// \author robert.sulej@cern.ch
//
#ifndef ANAB_MVAREADER_H
#define ANAB_MVAREADER_H

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Utilities/InputTag.h"

#include "lardata/ArtDataHelper/MVAWrapperBase.h"

namespace anab {

  template <class T, size_t N>
  class FVectorReader : public FVectorWrapperBase {
  public:
    static std::unique_ptr<FVectorReader> create(const art::Event& evt, const art::InputTag& tag)
    {
      bool success;
      std::unique_ptr<FVectorReader> ptr(new FVectorReader(evt, tag, success));
      if (success) { return ptr; }
      else {
        return nullptr;
      }
    }

    FVectorReader(const art::Event& evt, const art::InputTag& tag);

    T const& item(size_t key) const { return (*fDataHandle)[key]; }
    std::vector<T> const& items() const { return *fDataHandle; }

    std::vector<FeatureVector<N>> const& vectors() const { return *fVectors; }

    //const std::array<float, N> & getVector(size_t key) const { return (*fVectors)[key].data(); }

    std::array<float, N> getVector(size_t key) const
    {
      std::array<float, N> vout;
      for (size_t i = 0; i < N; ++i)
        vout[i] = (*fVectors)[key][i];
      return vout;
    }

    std::array<float, N> getVector(art::Ptr<T> const& item) const { return getVector(item.key()); }

    size_t size() const { return fVectors->size(); }

    size_t length() const { return N; }

    const std::string& dataTag() const { return fDescription->dataTag(); }

    const art::Handle<std::vector<T>>& dataHandle() const { return fDataHandle; }

    const std::string& columnName(size_t index) const { return fDescription->outputName(index); }

    int getIndex(const std::string& name) const { return fDescription->getIndex(name); }

    friend std::ostream& operator<<(std::ostream& o, FVectorReader const& a)
    {
      o << "FVectorReader:" << std::endl << *(a.fDescription) << std::endl;
      return o;
    }

  protected:
    FVectorReader(const art::Event& evt, const art::InputTag& tag, bool& success);

  private:
    FVecDescription<N> const* fDescription;
    std::vector<FeatureVector<N>> const* fVectors;
    art::Handle<std::vector<T>> fDataHandle;
  };

  template <class T, size_t N>
  class MVAReader : public FVectorReader<T, N>, public MVAWrapperBase {
  public:
    static std::unique_ptr<MVAReader> create(const art::Event& evt, const art::InputTag& tag)
    {
      bool success;
      std::unique_ptr<MVAReader> ptr(new MVAReader(evt, tag, success));
      if (success) { return ptr; }
      else {
        return nullptr;
      }
    }

    MVAReader(const art::Event& evt, const art::InputTag& tag) : FVectorReader<T, N>(evt, tag) {}

    std::vector<FeatureVector<N>> const& outputs() const { return FVectorReader<T, N>::vectors(); }

    std::array<float, N> getOutput(size_t key) const { return FVectorReader<T, N>::getVector(key); }

    std::array<float, N> getOutput(art::Ptr<T> const& item) const
    {
      return FVectorReader<T, N>::getVector(item.key());
    }

    std::array<float, N> getOutput(std::vector<art::Ptr<T>> const& items) const
    {
      return pAccumulate(items, FVectorReader<T, N>::vectors());
    }

    std::array<float, N> getOutput(std::vector<art::Ptr<T>> const& items,
                                   std::vector<float> const& weights) const
    {
      return pAccumulate(items, weights, FVectorReader<T, N>::vectors());
    }

    std::array<float, N> getOutput(std::vector<art::Ptr<T>> const& items,
                                   std::function<float(T const&)> fweight) const
    {
      return pAccumulate(items, fweight, FVectorReader<T, N>::vectors());
    }

    const std::string& outputName(size_t index) const
    {
      return FVectorReader<T, N>::columnName(index);
    }

  private:
    MVAReader(const art::Event& evt, const art::InputTag& tag, bool& success)
      : FVectorReader<T, N>(evt, tag, success)
    {}
  };

} // namespace anab

//----------------------------------------------------------------------------
// FVectorReader functions.
//
template <class T, size_t N>
anab::FVectorReader<T, N>::FVectorReader(const art::Event& evt, const art::InputTag& tag)
  : fDescription(0)
{
  if (!N) { throw cet::exception("FVectorReader") << "Vector size should be > 0." << std::endl; }

  auto descriptionHandle = evt.getValidHandle<std::vector<anab::FVecDescription<N>>>(tag);

  // search for FVecDescription<N> produced for the type T, with the instance name from the tag
  std::string outputInstanceName = tag.instance() + getProductName(typeid(T));
  for (auto const& dscr : *descriptionHandle) {
    if (dscr.outputInstance() == outputInstanceName) {
      fDescription = &dscr;
      break;
    }
  }
  if (!fDescription) {
    throw cet::exception("FVectorReader")
      << "Vectors description not found for " << outputInstanceName << std::endl;
  }

  fVectors = &*(evt.getValidHandle<std::vector<FeatureVector<N>>>(
    art::InputTag(tag.label(), fDescription->outputInstance(), tag.process())));

  if (!evt.getByLabel(fDescription->dataTag(), fDataHandle)) {
    throw cet::exception("FVectorReader")
      << "Associated data product handle failed: " << *(fDataHandle.whyFailed()) << std::endl;
  }

  if (fVectors->size() != fDataHandle->size()) {
    throw cet::exception("FVectorReader")
      << "Feature vectors and data products sizes inconsistent: " << fVectors->size()
      << "!=" << fDataHandle->size() << std::endl;
  }
}
//----------------------------------------------------------------------------

template <class T, size_t N>
anab::FVectorReader<T, N>::FVectorReader(const art::Event& evt,
                                         const art::InputTag& tag,
                                         bool& success)
  : fDescription(0)
{
  success = false; // until all is done correctly

  if (!N) {
    std::cout << "FVectorReader: Vector size should be > 0." << std::endl;
    return;
  }

  art::Handle<std::vector<anab::FVecDescription<N>>> descriptionHandle;
  if (!evt.getByLabel(tag, descriptionHandle)) { return; }

  // search for FVecDescription<N> produced for the type T, with the instance name from the tag
  std::string outputInstanceName = tag.instance() + getProductName(typeid(T));
  for (auto const& dscr : *descriptionHandle) {
    if (dscr.outputInstance() == outputInstanceName) {
      fDescription = &dscr;
      break;
    }
  }
  if (!fDescription) {
    std::cout << "FVectorReader: Vectors description not found for " << outputInstanceName
              << std::endl;
    return;
  }

  fVectors = &*(evt.getValidHandle<std::vector<FeatureVector<N>>>(
    art::InputTag(tag.label(), fDescription->outputInstance(), tag.process())));

  if (!evt.getByLabel(fDescription->dataTag(), fDataHandle)) {
    std::cout << "FVectorReader: Associated data product handle failed: "
              << *(fDataHandle.whyFailed()) << std::endl;
    return;
  }

  if (fVectors->size() != fDataHandle->size()) {
    std::cout << "FVectorReader: Feature vectors and data products sizes inconsistent: "
              << fVectors->size() << "!=" << fDataHandle->size() << std::endl;
    return;
  }

  success = true; // ok, all data found in the event
}
//----------------------------------------------------------------------------

#endif //ANAB_MVAREADER