ShowerTrackSpacePointDirection_tool.cc

Go to the documentation of this file.

//############################################################################
//### Name:        ShowerTrackSpacePointDirection                          ###
//### Author:      Dominic Barker                                          ###
//### Date:        13.05.19                                                ###
//### Description: Tool for finding the shower direction using the         ###
//###              the average direction of theinitial track spacepoints.  ###
//############################################################################

//Framework Includes
#include "art/Utilities/ToolMacros.h"

//LArSoft Includes
#include "larpandora/LArPandoraEventBuilding/LArPandoraShower/Tools/IShowerTool.h"

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

namespace ShowerRecoTools {

  class ShowerTrackSpacePointDirection : IShowerTool {

  public:
    ShowerTrackSpacePointDirection(const fhicl::ParameterSet& pset);

    //Calculate the direction using the initial track spacepoints
    int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
                         art::Event& Event,
                         reco::shower::ShowerElementHolder& ShowerEleHolder) override;

  private:
    //fcl
    int fVerbose;
    bool fUsePandoraVertex; //Direction from point defined as
    //(Position of SP - Vertex) rather than
    //(Position of SP - Track Start Point).

    std::string fInitialTrackSpacePointsInputLabel;
    std::string fShowerStartPositionInputLabel;
    std::string fInitialTrackInputLabel;
    std::string fShowerDirectionOutputLabel;
  };

  ShowerTrackSpacePointDirection::ShowerTrackSpacePointDirection(const fhicl::ParameterSet& pset)
    : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
    , fVerbose(pset.get<int>("Verbose"))
    , fUsePandoraVertex(pset.get<bool>("UsePandoraVertex"))
    , fInitialTrackSpacePointsInputLabel(pset.get<std::string>("InitialTrackSpacePointsInputLabel"))
    , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
    , fInitialTrackInputLabel(pset.get<std::string>("InitialTrackInputLabel"))
    , fShowerDirectionOutputLabel(pset.get<std::string>("ShowerDirectionOutputLabel"))
  {}

  int ShowerTrackSpacePointDirection::CalculateElement(
    const art::Ptr<recob::PFParticle>& pfparticle,
    art::Event& Event,
    reco::shower::ShowerElementHolder& ShowerEleHolder)
  {

    //Check the Track Hits has been defined
    if (!ShowerEleHolder.CheckElement(fInitialTrackSpacePointsInputLabel)) {
      if (fVerbose)
        mf::LogError("ShowerTrackSpacePointDirection")
          << "Initial track spacepoints not set" << std::endl;
      return 0;
    }

    //Check the start position is set.
    if (fUsePandoraVertex && !ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
      if (fVerbose)
        mf::LogError("ShowerTrackSpacePointDirection")
          << "Start position not set, returning " << std::endl;
      return 0;
    }

    //Get the start poistion
    geo::Point_t StartPosition = {-999, -999, -999};
    if (fUsePandoraVertex) {
      ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, StartPosition);
    }
    else {
      //Check the Tracks has been defined
      if (!ShowerEleHolder.CheckElement(fInitialTrackInputLabel)) {
        if (fVerbose)
          mf::LogError("ShowerTrackSpacePointDirection") << "Initial track not set" << std::endl;
        return 0;
      }
      recob::Track InitialTrack;
      ShowerEleHolder.GetElement(fInitialTrackInputLabel, InitialTrack);
      geo::Point_t Start_point = InitialTrack.Start();
      StartPosition = {Start_point.X(), Start_point.Y(), Start_point.Z()};
    }

    //Get the initial track hits.
    std::vector<art::Ptr<recob::SpacePoint>> intitaltrack_sp;
    ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, intitaltrack_sp);

    //Calculate the mean direction and the the standard deviation
    float sumX = 0, sumX2 = 0;
    float sumY = 0, sumY2 = 0;
    float sumZ = 0, sumZ2 = 0;

    //Get the spacepoints associated to the track hit
    for (auto const& sp : intitaltrack_sp) {

      //Get the direction relative to the start positon
      auto const pos = sp->position() - StartPosition;
      if (pos.R() == 0) { continue; }

      sumX = pos.X();
      sumX2 += pos.X() * pos.X();
      sumY = pos.Y();
      sumY2 += pos.Y() * pos.Y();
      sumZ = pos.Z();
      sumZ2 += pos.Z() * pos.Z();
    }

    float NumSps = intitaltrack_sp.size();
    auto const Mean = geo::Vector_t{sumX / NumSps, sumY / NumSps, sumZ / NumSps}.Unit();

    float RMSX = 999;
    float RMSY = 999;
    float RMSZ = 999;
    if (sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))) > 0) {
      RMSX = std::sqrt(sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))));
    }
    if (sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))) > 0) {
      RMSY = std::sqrt(sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))));
    }
    if (sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))) > 0) {
      RMSZ = std::sqrt(sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))));
    }

    //Loop over the spacepoints and remove ones the relative direction is not within one sigma.
    geo::Vector_t Direction_Mean = {0, 0, 0};
    int N = 0;
    for (auto const sp : intitaltrack_sp) {
      auto const Direction = sp->position() - StartPosition;
      if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
          (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
          (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
        if (Direction.R() == 0) { continue; }
        ++N;
        Direction_Mean += Direction;
      }
    }

    if (N > 0) {
      //Take the mean value
      Direction_Mean = Direction_Mean.Unit();
      ShowerEleHolder.SetElement(Direction_Mean, fShowerDirectionOutputLabel);
    }
    else {
      if (fVerbose)
        mf::LogError("ShowerTrackSpacePointDirection")
          << "None of the points are within 1 sigma" << std::endl;
      return 1;
    }
    return 0;
  }
}

DEFINE_ART_CLASS_TOOL(ShowerRecoTools::ShowerTrackSpacePointDirection)