LArPandoraTrackCreation_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"

#include "fhiclcpp/ParameterSet.h"

#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/TrackHitMeta.h"

#include "larpandoracontent/LArObjects/LArPfoObjects.h"

#include <memory>

namespace lar_pandora
{

class LArPandoraTrackCreation : public art::EDProducer
{
public:
    explicit LArPandoraTrackCreation(fhicl::ParameterSet const &pset);

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

    void produce(art::Event &evt) override;

private:
    recob::Track BuildTrack(const int id, const lar_content::LArTrackStateVector &trackStateVector) const;

    std::string     m_pfParticleLabel;              
    unsigned int    m_minTrajectoryPoints;          
    unsigned int    m_slidingFitHalfWindow;         
    bool            m_useAllParticles;              
};

DEFINE_ART_MODULE(LArPandoraTrackCreation)

} // namespace lar_pandora

//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows

#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"

#include "art/Persistency/Common/PtrMaker.h"

#include "canvas/Utilities/InputTag.h"

#include "larcore/Geometry/Geometry.h"

#include "lardata/Utilities/AssociationUtil.h"

#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Vertex.h"

#include "messagefacility/MessageLogger/MessageLogger.h"

#include "larpandoracontent/LArHelpers/LArPfoHelper.h"

#include "larpandora/LArPandoraInterface/LArPandoraHelper.h"

#include <iostream>

namespace lar_pandora
{

LArPandoraTrackCreation::LArPandoraTrackCreation(fhicl::ParameterSet const &pset) :
    m_pfParticleLabel(pset.get<std::string>("PFParticleLabel")),
    m_minTrajectoryPoints(pset.get<unsigned int>("MinTrajectoryPoints", 2)),
    m_slidingFitHalfWindow(pset.get<unsigned int>("SlidingFitHalfWindow", 20)),
    m_useAllParticles(pset.get<bool>("UseAllParticles", false))
{
    produces< std::vector<recob::Track> >();
    produces< art::Assns<recob::PFParticle, recob::Track> >();
    produces< art::Assns<recob::Track, recob::Hit> >();
    produces< art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();

    if (m_minTrajectoryPoints<2) throw cet::exception("LArPandoraTrackCreation") << "MinTrajectoryPoints should not be smaller than 2!";

}

//------------------------------------------------------------------------------------------------------------------------------------------

void LArPandoraTrackCreation::produce(art::Event &evt)
{
    std::unique_ptr< std::vector<recob::Track> > outputTracks( new std::vector<recob::Track> );
    std::unique_ptr< art::Assns<recob::PFParticle, recob::Track> > outputParticlesToTracks( new art::Assns<recob::PFParticle, recob::Track> );
    std::unique_ptr< art::Assns<recob::Track, recob::Hit> > outputTracksToHits( new art::Assns<recob::Track, recob::Hit> );
    std::unique_ptr< art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> > outputTracksToHitsWithMeta( new art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> );

    // 'wirePitchW` is here used only to provide length scale for binning hits and performing sliding/local linear fits.
    // Fits should be robust against the precise choice, provided length scale is comparable to the granularity of the images.
    art::ServiceHandle<geo::Geometry> theGeometry;
    const unsigned int nWirePlanes(theGeometry->MaxPlanes());

    if (nWirePlanes > 3)
        throw cet::exception("LArPandoraTrackCreation") << " LArPandoraTrackCreation::produce --- More than three wire planes present ";

    if ((0 == theGeometry->Ncryostats()) || (0 == theGeometry->NTPC(0)))
        throw cet::exception("LArPandoraTrackCreation") << " LArPandoraTrackCreation::produce --- unable to access first tpc in first cryostat ";

    std::unordered_set<geo::_plane_proj> planeSet;
    for (unsigned int iPlane = 0; iPlane < nWirePlanes; ++iPlane)
        (void) planeSet.insert(theGeometry->TPC(0, 0).Plane(iPlane).View());

    if ((nWirePlanes != planeSet.size()) || !planeSet.count(geo::kU) || !planeSet.count(geo::kV) || (planeSet.count(geo::kW) && planeSet.count(geo::kY)))
        throw cet::exception("LArPandoraTrackCreation") << " LArPandoraTrackCreation::produce --- expect to find u and v views; if there is one further view, it must be w or y ";

    const bool useYPlane((nWirePlanes > 2) && planeSet.count(geo::kY));

    const float wirePitchU(theGeometry->WirePitch(geo::kU));
    const float wirePitchV(theGeometry->WirePitch(geo::kV));
    const float wirePitchW((nWirePlanes < 3) ? 0.5f * (wirePitchU + wirePitchV) : (useYPlane) ? theGeometry->WirePitch(geo::kY) : theGeometry->WirePitch(geo::kW));

    int trackCounter(0);
    const art::PtrMaker<recob::Track> makeTrackPtr(evt, *this);

    // Organise inputs
    PFParticleVector pfParticleVector, extraPfParticleVector;
    PFParticlesToSpacePoints pfParticlesToSpacePoints;
    PFParticlesToClusters pfParticlesToClusters;
    LArPandoraHelper::CollectPFParticles(evt, m_pfParticleLabel, pfParticleVector, pfParticlesToSpacePoints);
    LArPandoraHelper::CollectPFParticles(evt, m_pfParticleLabel, extraPfParticleVector, pfParticlesToClusters);

    VertexVector vertexVector;
    PFParticlesToVertices pfParticlesToVertices;
    LArPandoraHelper::CollectVertices(evt, m_pfParticleLabel, vertexVector, pfParticlesToVertices);

    for (const art::Ptr<recob::PFParticle> pPFParticle : pfParticleVector)
    {
        // Select track-like pfparticles
        if (!m_useAllParticles && !LArPandoraHelper::IsTrack(pPFParticle))
            continue;

        // Obtain associated spacepoints
        PFParticlesToSpacePoints::const_iterator particleToSpacePointIter(pfParticlesToSpacePoints.find(pPFParticle));

        if (pfParticlesToSpacePoints.end() == particleToSpacePointIter)
        {
            mf::LogDebug("LArPandoraTrackCreation") << "No spacepoints associated to particle ";
            continue;
        }

        // Obtain associated clusters
        PFParticlesToClusters::const_iterator particleToClustersIter(pfParticlesToClusters.find(pPFParticle));

        if (pfParticlesToClusters.end() == particleToClustersIter)
        {
            mf::LogDebug("LArPandoraShowerCreation") << "No clusters associated to particle ";
            continue;
        }

        // Obtain associated vertex
        PFParticlesToVertices::const_iterator particleToVertexIter(pfParticlesToVertices.find(pPFParticle));

        if ((pfParticlesToVertices.end() == particleToVertexIter) || (1 != particleToVertexIter->second.size()))
        {
            mf::LogDebug("LArPandoraTrackCreation") << "Unexpected number of vertices for particle ";
            continue;
        }

        // Copy information into expected pandora form
        pandora::CartesianPointVector cartesianPointVector;
        for (const art::Ptr<recob::SpacePoint> spacePoint : particleToSpacePointIter->second)
            cartesianPointVector.emplace_back(pandora::CartesianVector(spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]));

        double vertexXYZ[3] = {0., 0., 0.};
        particleToVertexIter->second.front()->XYZ(vertexXYZ);
        const pandora::CartesianVector vertexPosition(vertexXYZ[0], vertexXYZ[1], vertexXYZ[2]);

        // Call pandora "fast" track fitter
        lar_content::LArTrackStateVector trackStateVector;
        pandora::IntVector indexVector;
        try
        {
            lar_content::LArPfoHelper::GetSlidingFitTrajectory(cartesianPointVector, vertexPosition, m_slidingFitHalfWindow, wirePitchW, trackStateVector, &indexVector);
        }
        catch (const pandora::StatusCodeException &)
        {
            mf::LogDebug("LArPandoraTrackCreation") << "Unable to extract sliding fit trajectory";
            continue;
        }

        if (trackStateVector.size() < m_minTrajectoryPoints)
        {
            mf::LogDebug("LArPandoraTrackCreation") << "Insufficient input trajectory points to build track: " << trackStateVector.size();
            continue;
        }

        HitVector hitsFromSpacePoints, hitsFromClusters, hitsInParticle;
        HitSet hitsInParticleSet;

        LArPandoraHelper::GetAssociatedHits(evt, m_pfParticleLabel, particleToSpacePointIter->second, hitsFromSpacePoints, &indexVector);
        LArPandoraHelper::GetAssociatedHits(evt, m_pfParticleLabel, particleToClustersIter->second, hitsFromClusters);
        //ATTN: hits ordered from space points if available, rest added at the end
        for (unsigned int hitIndex = 0; hitIndex < hitsFromSpacePoints.size(); hitIndex++)
        {
            hitsInParticle.push_back(hitsFromSpacePoints.at(hitIndex));
            (void) hitsInParticleSet.insert(hitsFromSpacePoints.at(hitIndex));
        }

        for (unsigned int hitIndex = 0; hitIndex < hitsFromClusters.size(); hitIndex++)
        {
            if (hitsInParticleSet.count(hitsFromClusters.at(hitIndex)) == 0)
                hitsInParticle.push_back(hitsFromClusters.at(hitIndex));
        }

        // Add invalid points at the end of the vector, so that the number of the trajectory points is the same as the number of hits
        if (trackStateVector.size()>hitsFromSpacePoints.size())
        {
            throw cet::exception("LArPandoraTrackCreation") << "trackStateVector.size() is greater than hitsFromSpacePoints.size()";
        }
        const unsigned int nInvalidPoints = hitsInParticle.size()-trackStateVector.size();
        for (unsigned int i=0;i<nInvalidPoints;++i) {
            trackStateVector.push_back(lar_content::LArTrackState(pandora::CartesianVector(util::kBogusF,util::kBogusF,util::kBogusF),
                                                                  pandora::CartesianVector(util::kBogusF,util::kBogusF,util::kBogusF), nullptr));
        }

        // Output objects
        outputTracks->emplace_back(LArPandoraTrackCreation::BuildTrack(trackCounter++, trackStateVector));
        art::Ptr<recob::Track> pTrack(makeTrackPtr(outputTracks->size() - 1));

        // Output associations, after output objects are in place
        util::CreateAssn(*this, evt, pTrack, pPFParticle, *(outputParticlesToTracks.get()));
        util::CreateAssn(*this, evt, *(outputTracks.get()), hitsInParticle, *(outputTracksToHits.get()));

        //ATTN: metadata added with index from space points if available, null for others
        for (unsigned int hitIndex = 0; hitIndex < hitsInParticle.size(); hitIndex++)
        {
            const art::Ptr<recob::Hit> pHit(hitsInParticle.at(hitIndex));
            const int index((hitIndex < hitsFromSpacePoints.size()) ? hitIndex : std::numeric_limits<int>::max());
            recob::TrackHitMeta metadata(index, -std::numeric_limits<double>::max());
            outputTracksToHitsWithMeta->addSingle(pTrack, pHit, metadata);
        }
    }

    mf::LogDebug("LArPandoraTrackCreation") << "Number of new tracks: " << outputTracks->size() << std::endl;

    evt.put(std::move(outputTracks));
    evt.put(std::move(outputTracksToHits));
    evt.put(std::move(outputTracksToHitsWithMeta));
    evt.put(std::move(outputParticlesToTracks));
}

//------------------------------------------------------------------------------------------------------------------------------------------

recob::Track LArPandoraTrackCreation::BuildTrack(const int id, const lar_content::LArTrackStateVector &trackStateVector) const
{
    if (trackStateVector.empty())
        throw cet::exception("LArPandoraTrackCreation") << "BuildTrack - No input trajectory points provided ";

    recob::tracking::Positions_t xyz;
    recob::tracking::Momenta_t pxpypz;
    recob::TrackTrajectory::Flags_t flags;

    for (const lar_content::LArTrackState &trackState : trackStateVector)
    {
        xyz.emplace_back(recob::tracking::Point_t(trackState.GetPosition().GetX(), trackState.GetPosition().GetY(), trackState.GetPosition().GetZ()));
        pxpypz.emplace_back(recob::tracking::Vector_t(trackState.GetDirection().GetX(), trackState.GetDirection().GetY(), trackState.GetDirection().GetZ()));
        // Set flag NoPoint if point has bogus coordinates, otherwise use clean flag set
        if (std::fabs(trackState.GetPosition().GetX()-util::kBogusF)<std::numeric_limits<float>::epsilon() &&
            std::fabs(trackState.GetPosition().GetY()-util::kBogusF)<std::numeric_limits<float>::epsilon() &&
            std::fabs(trackState.GetPosition().GetZ()-util::kBogusF)<std::numeric_limits<float>::epsilon())
        {
            flags.emplace_back(recob::TrajectoryPointFlags(recob::TrajectoryPointFlags::InvalidHitIndex, recob::TrajectoryPointFlagTraits::NoPoint));
        } else {
            flags.emplace_back(recob::TrajectoryPointFlags());
        }
    }

    // note from gc: eventually we should produce a TrackTrajectory, not a Track with empty covariance matrix and bogus chi2, etc.
    return recob::Track(recob::TrackTrajectory(std::move(xyz), std::move(pxpypz), std::move(flags), false),
                        util::kBogusI, util::kBogusF, util::kBogusI, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), id);
}

} // namespace lar_pandora