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