ShowerFinder_module.cc

Go to the documentation of this file.

// ShowerFinder class
//
// \author roxanne.guenette@yale.edu
//

#include <string>

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"

#include <math.h>

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft Includes
#include "larcore/Geometry/Geometry.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/EndPoint2D.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Shower.h"

// ROOT
#include "TMath.h"

namespace shwf {

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

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

    std::string fVertexModuleLabel;         
    std::string fClusterModuleLabel;        
    std::string fHoughLineModuleLabel;      
    std::string fVertexStrengthModuleLabel; 
    double fRcone;                          
    double fLcone;                          
  };                                        // class ShowerFinder

}

namespace shwf {

  //-------------------------------------------------
  ShowerFinder::ShowerFinder(fhicl::ParameterSet const& pset) : EDProducer{pset}
  {
    fVertexModuleLabel = pset.get<std::string>("VertexModuleLabel");
    fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
    fHoughLineModuleLabel = pset.get<std::string>("HoughLineModuleLabel");
    fVertexStrengthModuleLabel = pset.get<std::string>("VertexStrengthModuleLabel");
    fRcone = pset.get<double>("Rcone");
    fLcone = pset.get<double>("Lcone");

    produces<std::vector<recob::Shower>>();
    produces<art::Assns<recob::Shower, recob::Cluster>>();
    produces<art::Assns<recob::Shower, recob::Hit>>();
  }

  //
  //-------------------------------------------------
  void ShowerFinder::produce(art::Event& evt)
  {

    std::unique_ptr<std::vector<recob::Shower>> showercol(new std::vector<recob::Shower>);
    std::unique_ptr<art::Assns<recob::Shower, recob::Cluster>> cassn(
      new art::Assns<recob::Shower, recob::Cluster>);
    std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> hassn(
      new art::Assns<recob::Shower, recob::Hit>);

    // Read in the vertex List object(s).
    art::Handle<std::vector<recob::EndPoint2D>> vertexListHandle;
    evt.getByLabel(fVertexModuleLabel, vertexListHandle);

    // Read in the hough line List object(s).
    art::Handle<std::vector<recob::Cluster>> houghListHandle;
    evt.getByLabel(fHoughLineModuleLabel, houghListHandle);

    // Read in the cluster List object(s).
    art::Handle<std::vector<recob::Cluster>> clusterListHandle;
    evt.getByLabel(fClusterModuleLabel, clusterListHandle);

    art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterModuleLabel);

    // Read in the vertex Strength List object(s).
    art::Handle<std::vector<recob::EndPoint2D>> vertexStrengthListHandle;
    evt.getByLabel(fVertexStrengthModuleLabel, vertexStrengthListHandle);

    std::vector<size_t> protoShowers; //vector of indices of clusters associated to a cone

    std::vector<art::Ptr<recob::Hit>> clusterhits; //hits in the cluster

    art::ServiceHandle<geo::Geometry const> geom;

    //This vector will contain all strong and strongest vertices
    art::PtrVector<recob::EndPoint2D> vertSel;

    //This loop is going over all the vertices in the event
    //and is interested in ONLY strong and strongest vertices.
    mf::LogInfo("ShowerFinder") << "Vertex STRENGTH list size = "
                                << vertexStrengthListHandle->size()
                                << " AND vertices:" << vertexListHandle->size()
                                << "\nCLUSTER list size = " << clusterListHandle->size()
                                << " AND Hough: :" << houghListHandle->size();

    for (size_t iv = 0; iv < vertexListHandle->size(); ++iv) {
      art::Ptr<recob::EndPoint2D> vertex(vertexListHandle, iv);
      //std::cout << "Vertex " << iv << " :  str = " << vertex->ID() << std::endl;
      //if(vertex->Strength() == 4 || vertex->Strength() == 3){
      if (vertex->ID() == 1 || vertex->Strength() == 3) { //only use Strongest and strong
        vertSel.push_back(vertex);
      }
      else
        continue;
    }

    //Definition of the geometry of the cone (which is basically a triangle)
    double scan_angle = 0; //angle of the scan steps
    double xa_cone = 0;    // x coordinate of the cone's apex (wire number)
    double ya_cone = 0;    // y coordinate of the cone's apex (drift time)
    double x1_cone = 0;    // x coordinate of the cone's top right point (wire number)
    double y1_cone = 0;    // y coordinate of the cone's top right point (drift time)
    double x2_cone = 0;    // x coordinate of the cone's top left point (wire number)
    double y2_cone = 0;    // y coordinate of the cone's top left point  (drift time)

    //The length of the side of the cone
    double fScone = std::sqrt((fRcone * fRcone) + (fLcone * fLcone));

    // Opening angle of the cone (defined from input parameters)
    double cone_angle = (TMath::ATan(fRcone / fLcone)) / 2.0;
    mf::LogInfo("ShowerFinder") << "Cone Opening Angle: " << (180.0 * cone_angle) / TMath::Pi();
    double compl_angle = 0;

    unsigned int n_scan = 1 + (int)(TMath::Pi() / (2.0 * cone_angle));
    mf::LogInfo("ShowerFinder") << "N scan: " << n_scan;

    double x_hit = 0; //x coordinate of hit
    double y_hit = 0; //y coordinate of hit

    int hits_cluster_counter = 0; //count the number of hits in a cluster that is inside a cone
    //int hits_cluster_Total = 0; //The total number of hits in a cluster

    // For EVERY vertex, the algorithm is going to scan the plane to find clusters
    // contained in the scanning cones

    for (size_t ivert = 0; ivert < vertSel.size(); ++ivert) {

      mf::LogInfo("ShowerFinder") << "Number of STRONG vertices = " << vertSel.size();

      //get the coordinates of the vertex for the summit of the cone
      xa_cone = vertSel[ivert]->WireID().Wire; //for update to EndPoint2D ... WK 4/22/13
      ya_cone = vertSel[ivert]->DriftTime();

      mf::LogInfo("ShowerFinder") << "Vertex at: (" << xa_cone << ", " << ya_cone << ")";

      //Beginning of the scan!
      for (unsigned int iscan = 0; iscan < n_scan; ++iscan) {

        mf::LogInfo("ShowerFinder") << ">>>> Start SCAN: " << iscan;

        //define the scan anlge
        scan_angle = (TMath::Pi() / 2.0) - (iscan * (2.0 * cone_angle));

        mf::LogInfo("ShowerFinder") << "Scan Angle: " << (180. * scan_angle) / TMath::Pi();

        //get the complementary angle for geometry puurposes
        compl_angle = scan_angle - cone_angle;

        //Calculate the coordinates of the top right corner of the cone
        x1_cone = xa_cone + fScone * (std::cos(compl_angle));
        y1_cone = ya_cone + fScone * (std::sin(compl_angle));

        //Calculate the coordinates of the top left corner of the cone
        x2_cone = xa_cone + fScone * (std::cos(scan_angle + cone_angle));
        y2_cone = ya_cone + fScone * (std::sin(scan_angle + cone_angle));

        //Looking if a cluster is in this cone (loop over all hits of all clusters)
        protoShowers.clear();
        for (size_t iclust = 0; iclust < clusterListHandle->size(); ++iclust) {

          //  art::Ptr<recob::Cluster> clust(clusterListHandle, iclust);
          recob::Cluster const& clust = clusterListHandle->at(iclust);

          //Get the hits vector from the cluster
          clusterhits = fmh.at(iclust);
          if (clusterhits.size() == 0) continue;

          //Loop over ALL hits in the cluster. Looking if the cluster's
          // hit is comprised in the cone
          for (size_t ihits = 0; ihits < clusterhits.size(); ++ihits) {

            x_hit = clusterhits[ihits]->WireID().Wire; //for update to EndPoint2D ... WK 4/22/13
            y_hit = clusterhits[ihits]->PeakTime();

            // Check in hits is INSIDE cone

            //define the 2 line equations:
            if (y_hit <= ((y2_cone - ya_cone) / (x2_cone - xa_cone)) * x_hit + ya_cone &&
                y_hit >= ((y1_cone - ya_cone) / (x1_cone - xa_cone)) * x_hit + ya_cone) {
              hits_cluster_counter++;
            }

          } //end hits loop

          //If there is more than 50% if the cluster INSIDE the cone, this is a protoshower
          if (clusterhits.size() == 0) continue;
          if (((double)hits_cluster_counter / (double)clusterhits.size()) >= 0.5) {
            mf::LogInfo("ShowerFinder") << "GOT A SHOWER!!!  in scan " << iscan
                                        << "  cluster: " << iclust << " : " << clust.ID();


            protoShowers.push_back(iclust);
          }
          clusterhits.clear();
          hits_cluster_counter = 0;

        } //end cluster loop

        if (protoShowers.empty()) continue;

        showercol->push_back(recob::Shower());

        // associate the shower with its clusters
        util::CreateAssn(
          evt, *cassn, showercol->size() - 1, protoShowers.begin(), protoShowers.end());

        // get the hits associated with each cluster and associate those with the shower
        for (size_t p = 0; p < protoShowers.size(); ++p) {
          const size_t psIndex = protoShowers[p];
          std::vector<art::Ptr<recob::Hit>> hits = fmh.at(psIndex);
          util::CreateAssn(evt, *showercol, hits, *hassn);
        }

      } //end scan loop
    }   //end vertices loop

    //NEED TO SEPARATE THE SHOWERS FROM THE DIFFERENT VERTEX!!!!

    mf::LogInfo("ShowerFinder") << "--->  Recorded shower =  " << showercol->size();
    //if the shower is stand alone ok, else, erase the next one
    //shower.SetID(is);
    //shower.SetVertexCoord(xa_cone, ya_cone);

    vertSel.clear();

    evt.put(std::move(showercol));
    evt.put(std::move(cassn));
    evt.put(std::move(hassn));

    return;
  } // end of produce

} // end of namespace

namespace shwf {

  DEFINE_ART_MODULE(ShowerFinder)

} //end of shower namespace