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 "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h" 
#include "art/Framework/Services/Optional/TFileService.h" 
#include "art/Framework/Services/Optional/TFileDirectory.h" 
#include "canvas/Persistency/Common/FindManyP.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// larsoft
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "lardata/ArtDataHelper/TrackUtils.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "larcorealg/Geometry/CryostatGeo.h"
#include "larcorealg/Geometry/TPCGeo.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/Shower.h"
#include "larreco/RecoAlg/ProjectionMatchingAlg.h"
#include "larreco/RecoAlg/PMAlg/PmaTrack3D.h"
#include "larreco/RecoAlg/PMAlg/Utilities.h"
#include "larreco/RecoAlg/ShowerEnergyAlg.h"

// C++
#include <iostream>
#include <map>
#include <iterator>

// ROOT
#include "TVector2.h"
#include "TMath.h"

//temp
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TLine.h"
#include "TCanvas.h"
#include "TString.h"
#include "TF1.h"
#include "larsim/MCCheater/BackTrackerService.h"
#include "TH1I.h"
#include "TText.h"
#include "TFile.h"
#include "TPrincipal.h"
#include "TProfile.h"

namespace shower {
  class EMShowerAlg;
  class HitPosition;
}

class shower::EMShowerAlg {
public:

  EMShowerAlg(fhicl::ParameterSet const& pset);

  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);

  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);

  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);

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

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

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

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

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

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

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

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

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

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

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

  int fDebug;

private:

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

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

  TVector3 Construct3DPoint(art::Ptr<recob::Hit> const& hit1, art::Ptr<recob::Hit> const& hit2);

  double FinddEdx(std::vector<art::Ptr<recob::Hit> > const& trackHits, std::unique_ptr<recob::Track> const& track);

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

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

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

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

  TVector2 HitCoordinates(art::Ptr<recob::Hit> const& hit);

  TVector2 HitPosition(art::Ptr<recob::Hit> const& hit);

  TVector2 HitPosition(TVector2 const& pos, geo::PlaneID planeID);

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

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

  TVector2 Project3DPointOntoPlane(TVector3 const& point, int plane, int cryostat = 0);

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

  TVector2 ShowerCentre(const std::vector<art::Ptr<recob::Hit> >& showerHits);

  TVector2 ShowerDirection(const std::vector<art::Ptr<recob::Hit> >& showerHits);

  double ShowerHitRMS(const std::vector<art::Ptr<recob::Hit> >& showerHits);

  double ShowerHitRMSGradient(const std::vector<art::Ptr<recob::Hit> >& showerHits, TVector2 trueStart = TVector2(0,0));

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


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

  // Services used by this class
  art::ServiceHandle<geo::Geometry> fGeom;
  detinfo::DetectorProperties const* fDetProp;
  art::ServiceHandle<art::TFileService> tfs;

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

  std::string fDetector;


  // tmp
  int FindTrueParticle(const std::vector<art::Ptr<recob::Hit> >& showerHits);
  int FindParticleID(const art::Ptr<recob::Hit>& hit);
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  TH1I* hTrueDirection;
  TProfile* hNumHitsInSegment, *hNumSegments;
  void MakePicture();
  bool fMakeGradientPlot, fMakeRMSGradientPlot;
  int fNumShowerSegments;

};

class shower::HitPosition {
 public:

  TVector2 WireTick;
  TVector2 Cm;

  // default constructor
  HitPosition();
  // contructors
  HitPosition(TVector2 wiretick, TVector2 cm): HitPosition()
  {
    WireTick = wiretick;
    Cm = cm;
  }
  HitPosition(TVector2 wiretick, geo::PlaneID planeID): HitPosition() {
    WireTick = wiretick;
    Cm = ConvertWireTickToCm(wiretick, planeID);
  }

  TVector2 ConvertWireTickToCm(TVector2 wiretick, geo::PlaneID planeID) {
    return TVector2(wiretick.X() * fGeom->WirePitch(planeID),
                    fDetProp->ConvertTicksToX(wiretick.Y(), planeID));
  }

 private:

  geo::GeometryCore const* fGeom = nullptr;
  detinfo::DetectorProperties const* fDetProp = nullptr;

};

#endif