EMShowerAlg.h

Go to the documentation of this file.

// Implementation of the EMShower algorithm
//
// Forms EM showers from clusters and associated tracks.
// Also provides methods for finding the vertex and further
// properties of the shower.
//
// Mike Wallbank (m.wallbank@sheffield.ac.uk), September 2015

#ifndef EMShowerAlg_hxx
#define EMShowerAlg_hxx

// framework
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/fwd.h"

// larsoft
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardataalg/DetectorInfo/DetectorProperties.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Shower.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "larreco/RecoAlg/ProjectionMatchingAlg.h"
#include "larreco/RecoAlg/ShowerEnergyAlg.h"

// C++
#include <map>
#include <memory>
#include <string>
#include <vector>

// ROOT
#include "RtypesCore.h"
#include "TVector2.h"
#include "TVector3.h"

namespace detinfo {
  class DetectorClocksData;
  class DetectorPropertiesData;
}

namespace shower {
  class EMShowerAlg;
  struct HitPosition;
}

class shower::EMShowerAlg {
public:
  EMShowerAlg(fhicl::ParameterSet const& pset, int debug = 0);

  void AssociateClustersAndTracks(std::vector<art::Ptr<recob::Cluster>> const& clusters,
                                  art::FindManyP<recob::Hit> const& fmh,
                                  art::FindManyP<recob::Track> const& fmt,
                                  std::map<int, std::vector<int>>& clusterToTracks,
                                  std::map<int, std::vector<int>>& trackToClusters) const;

  void AssociateClustersAndTracks(std::vector<art::Ptr<recob::Cluster>> const& clusters,
                                  art::FindManyP<recob::Hit> const& fmh,
                                  art::FindManyP<recob::Track> const& fmt,
                                  std::vector<int> const& clustersToIgnore,
                                  std::map<int, std::vector<int>>& clusterToTracks,
                                  std::map<int, std::vector<int>>& trackToClusters) const;

  std::vector<int> CheckShowerPlanes(std::vector<std::vector<int>> const& initialShowers,
                                     std::vector<art::Ptr<recob::Cluster>> const& clusters,
                                     art::FindManyP<recob::Hit> const& fmh) const;

  std::unique_ptr<recob::Track> ConstructTrack(
    detinfo::DetectorPropertiesData const& detProp,
    std::vector<art::Ptr<recob::Hit>> const& track1,
    std::vector<art::Ptr<recob::Hit>> const& track2,
    std::map<int, TVector2> const& showerCentreMap) const;

  std::unique_ptr<recob::Track> ConstructTrack(
    detinfo::DetectorPropertiesData const& detProp,
    std::vector<art::Ptr<recob::Hit>> const& track1,
    std::vector<art::Ptr<recob::Hit>> const& track2) const;

  void FindInitialTrack(detinfo::DetectorPropertiesData const& detProp,
                        std::map<int, std::vector<art::Ptr<recob::Hit>>> const& hits,
                        std::unique_ptr<recob::Track>& initialTrack,
                        std::map<int, std::vector<art::Ptr<recob::Hit>>>& initialTrackHits) const;

  std::vector<std::vector<int>> FindShowers(
    std::map<int, std::vector<int>> const& trackToClusters) const;

  recob::Shower MakeShower(
    detinfo::DetectorClocksData const& clockData,
    detinfo::DetectorPropertiesData const& detProp,
    art::PtrVector<recob::Hit> const& hits,
    std::unique_ptr<recob::Track> const& initialTrack,
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialTrackHits) const;

  recob::Shower MakeShower(detinfo::DetectorClocksData const& clockData,
                           detinfo::DetectorPropertiesData const& detProp,
                           art::PtrVector<recob::Hit> const& hits,
                           art::Ptr<recob::Vertex> const& vertex,
                           int& iok) const;

  std::vector<recob::SpacePoint> MakeSpacePoints(
    detinfo::DetectorPropertiesData const& detProp,
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& hits,
    std::vector<std::vector<art::Ptr<recob::Hit>>>& hitAssns) const;

  std::map<int, std::vector<art::Ptr<recob::Hit>>> OrderShowerHits(
    detinfo::DetectorPropertiesData const& detProp,
    art::PtrVector<recob::Hit> const& shower,
    int plane) const;

  void FindInitialTrackHits(std::vector<art::Ptr<recob::Hit>> const& showerHits,
                            art::Ptr<recob::Vertex> const& vertex,
                            std::vector<art::Ptr<recob::Hit>>& trackHits) const;

  Int_t WeightedFit(const Int_t n,
                    const Double_t* x,
                    const Double_t* y,
                    const Double_t* w,
                    Double_t* parm) const;

  bool isCleanShower(std::vector<art::Ptr<recob::Hit>> const& hits) const;

  int fDebug;

private:
  void CheckIsolatedHits_(std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const;

  bool CheckShowerHits_(
    detinfo::DetectorPropertiesData const& detProp,
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const;

  geo::Point_t Construct3DPoint_(detinfo::DetectorPropertiesData const& detProp,
                                 art::Ptr<recob::Hit> const& hit1,
                                 art::Ptr<recob::Hit> const& hit2) const;

  double FinddEdx_(detinfo::DetectorClocksData const& clockData,
                   detinfo::DetectorPropertiesData const& detProp,
                   std::vector<art::Ptr<recob::Hit>> const& trackHits,
                   std::unique_ptr<recob::Track> const& track) const;

  std::vector<art::Ptr<recob::Hit>> FindOrderOfHits_(detinfo::DetectorPropertiesData const& detProp,
                                                     std::vector<art::Ptr<recob::Hit>> const& hits,
                                                     bool perpendicular = false) const;

  std::map<int, std::vector<art::Ptr<recob::Hit>>> FindShowerStart_(
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& orderedShowerMap) const;

  std::map<int, std::map<int, bool>> GetPlanePermutations_(
    const detinfo::DetectorPropertiesData& detProp,
    const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const;

  double GlobalWire_(const geo::WireID& wireID) const;

  TVector2 HitCoordinates_(recob::Hit const& hit) const;

  TVector2 HitPosition_(detinfo::DetectorPropertiesData const& detProp,
                        recob::Hit const& hit) const;

  TVector2 HitPosition_(detinfo::DetectorPropertiesData const& detProp,
                        TVector2 const& pos,
                        geo::PlaneID planeID) const;

  std::unique_ptr<recob::Track> MakeInitialTrack_(
    detinfo::DetectorPropertiesData const& detProp,
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialHitsMap,
    std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const;

  void OrderShowerHits_(detinfo::DetectorPropertiesData const& detProp,
                        std::vector<art::Ptr<recob::Hit>> const& shower,
                        std::vector<art::Ptr<recob::Hit>>& orderedShower,
                        art::Ptr<recob::Vertex> const& vertex) const;

  TVector2 Project3DPointOntoPlane_(detinfo::DetectorPropertiesData const& detProp,
                                    geo::Point_t const& point,
                                    int plane,
                                    int cryostat = 0) const;

  std::map<double, int> RelativeWireWidth_(
    const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const;

  TVector2 ShowerCenter_(detinfo::DetectorPropertiesData const& detProp,
                         std::vector<art::Ptr<recob::Hit>> const& showerHits) const;

  TVector2 ShowerDirection_(detinfo::DetectorPropertiesData const& detProp,
                            const std::vector<art::Ptr<recob::Hit>>& showerHits) const;

  double ShowerHitRMS_(detinfo::DetectorPropertiesData const& detProp,
                       std::vector<art::Ptr<recob::Hit>> const& showerHits) const;

  double ShowerHitRMSGradient_(detinfo::DetectorPropertiesData const& detProp,
                               const std::vector<art::Ptr<recob::Hit>>& showerHits) const;

  int WorstPlane_(const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const;

  // Parameters
  double const fMinTrackLength;
  double const fdEdxTrackLength;
  double const fSpacePointSize;

  // Parameters to fit wire vs time
  unsigned int const fNfitpass;
  std::vector<unsigned int> const fNfithits;
  std::vector<double> const fToler;

  // Services used by this class
  art::ServiceHandle<geo::Geometry const> fGeom;
  geo::WireReadoutGeom const* fWireReadoutGeom{
    &art::ServiceHandle<geo::WireReadout const>()->Get()};

  // Algs used by this class
  shower::ShowerEnergyAlg const fShowerEnergyAlg;
  calo::CalorimetryAlg const fCalorimetryAlg;
  pma::ProjectionMatchingAlg const fProjectionMatchingAlg;

  std::string const fDetector;

  // tmp
  bool const fMakeGradientPlot;
  bool const fMakeRMSGradientPlot;
  int const fNumShowerSegments;
};

struct shower::HitPosition {
  TVector2 WireTick;
  TVector2 Cm;
};

#endif