DBCluster3D_module.cc

Go to the documentation of this file.

// Class:       DBCluster3D
// Plugin Type: producer (art v2_11_02)
// File:        DBCluster3D_module.cc
//
// Generated at Sat Jun 16 07:43:32 2018 by Tingjun Yang using cetskelgen
// from cetlib version v3_03_01.

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Slice.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "larreco/RecoAlg/DBScan3DAlg.h"
#include "larcore/Geometry/Geometry.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"

#include <memory>

namespace cluster {
  class DBCluster3D;
}


class cluster::DBCluster3D : public art::EDProducer {
  
  using Point_t = recob::tracking::Point_t;
  using Vector_t  = recob::tracking::Vector_t;
  
public:
  explicit DBCluster3D(fhicl::ParameterSet const & p);
  // The compiler-generated destructor is fine for non-base
  // classes without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  DBCluster3D(DBCluster3D const &) = delete;
  DBCluster3D(DBCluster3D &&) = delete;
  DBCluster3D & operator = (DBCluster3D const &) = delete;
  DBCluster3D & operator = (DBCluster3D &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;

private:

  const art::InputTag fHitModuleLabel;
  const art::InputTag fSpacePointModuleLabel;
  const art::InputTag fSPHitAssnLabel;

  DBScan3DAlg fDBScan;

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

  double tickToDist;
  double fMinHitDis;

};


cluster::DBCluster3D::DBCluster3D(fhicl::ParameterSet const & p)
  : fHitModuleLabel(p.get< art::InputTag >("HitModuleLabel"))
  , fSpacePointModuleLabel(p.get< art::InputTag >("SpacePointModuleLabel"))
  , fSPHitAssnLabel(p.get< art::InputTag >("SPHitAssnLabel"))
  , fDBScan(p.get< fhicl::ParameterSet >("DBScan3DAlg"))
  , fMinHitDis(p.get< double >("MinHitDis"))
{
  produces< std::vector<recob::Slice> >();
  produces< art::Assns<recob::Slice, recob::Hit> >();
  produces< art::Assns<recob::Slice, recob::SpacePoint> >();

  fGeom = &*(art::ServiceHandle<geo::Geometry>());
  fDetProp = lar::providerFrom<detinfo::DetectorPropertiesService>();

  tickToDist = fDetProp->DriftVelocity(fDetProp->Efield(),fDetProp->Temperature());
  tickToDist *= 1.e-3 * fDetProp->SamplingRate(); // 1e-3 is conversion of 1/us to 1/ns     
  fMinHitDis *= fMinHitDis;
}

void cluster::DBCluster3D::produce(art::Event & evt)
{

  std::vector<recob::Slice> slcCol;

  std::unique_ptr<art::Assns<recob::Slice, recob::Hit>> 
    slc_hit_assn(new art::Assns<recob::Slice, recob::Hit>);

  std::unique_ptr<art::Assns<recob::Slice, recob::SpacePoint>> 
    slc_sps_assn(new art::Assns<recob::Slice, recob::SpacePoint>);

  auto hitsHandle = evt.getValidHandle< std::vector<recob::Hit> >(fHitModuleLabel);
  auto spsHandle = evt.getValidHandle< std::vector<recob::SpacePoint> >(fSpacePointModuleLabel);

  // all hits in the collection
  std::vector< art::Ptr<recob::Hit> > hits;
  art::fill_ptr_vector(hits, hitsHandle);

  // all space points in the collection
  std::vector< art::Ptr<recob::SpacePoint> > sps;
  art::fill_ptr_vector(sps, spsHandle);

  art::FindManyP< recob::SpacePoint > spFromHit(hitsHandle, evt, fSPHitAssnLabel);
  if(!spFromHit.isValid()) {
    std::cout<<"spFromHit is invalid\n";
    return;
  }
  // Find the first Hit - SpacePoint assn and check consistency on the first event
  static bool first = true;
  if(first) {
    bool success = false;
    for(auto &hit : hits) {
      auto &sps = spFromHit.at(hit.key());
      if(sps.empty()) continue;
      success = (sps[0].id() == spsHandle.id());
      break;
    } // hit
    if(!success) throw cet::exception("DBCluster3D")<<"HitModuleLabel, SpacePointModuleLabel and SPHitAssnLabel are inconsistent\n";
    first = false;
  } // first

  art::FindManyP< recob::Hit > hitFromSp(spsHandle, evt, fSPHitAssnLabel);
  if(!hitFromSp.isValid()) {
    std::cout<<"hitFromSp is invalid\n";
    return;
  }
  fDBScan.init(sps, hitFromSp);
  fDBScan.dbscan();

  //Find number of slices
  int maxid = INT_MIN;
  for (size_t i = 0; i<fDBScan.points.size(); ++i){
//    std::cout<<"Space point index "<<i<<" "<<fDBScan.points[i].sp.key()<<" "<<fDBScan.points[i].cluster_id<<std::endl;
    if (fDBScan.points[i].cluster_id > maxid) maxid = fDBScan.points[i].cluster_id;
  }
  size_t nslc = 0;
  if (maxid>=0) nslc = maxid + 1;

  //Save hits associated with each slice
  std::vector<std::vector<art::Ptr<recob::Hit>>> slcHits(nslc);
  //Save hits on each PlaneID with pfparticle index
  std::map<geo::PlaneID, std::vector<std::pair<art::Ptr<recob::Hit>, unsigned int>>> hitmap;
  for (auto &hit : hits){
    auto &sps = spFromHit.at(hit.key());
    if (sps.size()){//found associated space point
      if (fDBScan.points[sps[0].key()].cluster_id>=0){
        slcHits[fDBScan.points[sps[0].key()].cluster_id].push_back(hit);
        hitmap[geo::PlaneID(hit->WireID())].push_back(std::make_pair(hit, fDBScan.points[sps[0].key()].cluster_id));
      }
    } // sps.size()
  } // hit

  //Save hits not associated with any spacepoints
  for (auto &hit : hits){
    bool found = false;
    for (size_t i = 0; i<slcHits.size(); ++i){
      if (std::find(slcHits[i].begin(), slcHits[i].end(), hit) != slcHits[i].end()){
        found = true;
        break;
      }
    }
    if (!found){
      double wirePitch = fGeom->WirePitch(hit->WireID());
      double UnitsPerTick = tickToDist / wirePitch;
      double x0 = hit->WireID().Wire;
      double y0 = hit->PeakTime() * UnitsPerTick;
      double mindis = DBL_MAX;
      unsigned slcIndex = UINT_MAX;
      for (auto &hit2 : hitmap[geo::PlaneID(hit->WireID())]) {
        double dx = hit2.first->WireID().Wire - x0;
        double dy = hit2.first->PeakTime() * UnitsPerTick - y0;
        double dis = dx*dx + dy*dy;
        if (dis<mindis){
          mindis = dis;
          slcIndex = hit2.second;
        }
      }
      if (slcIndex != UINT_MAX && mindis < fMinHitDis){
        slcHits[slcIndex].push_back(hit);
      }
    }
  }

  //Save spacepoints for each slice
  std::vector<std::vector<art::Ptr<recob::SpacePoint>>> sps_in_slc(nslc);
  for (size_t i = 0; i<fDBScan.points.size(); ++i){
    if (fDBScan.points[i].cluster_id>=0){
      sps_in_slc[fDBScan.points[i].cluster_id].push_back(sps[i]);
    }
  } // i

  // calculate the properties of the slice
  for(size_t isl = 0; isl < nslc; ++isl) {
    double sum = sps_in_slc[isl].size();
    // find the center
    double center[3] = {0.};
    // TODO: calculate the charge using recob::PointCharge instead of just counting hits
    float charge = slcHits[isl].size();
    for(auto& spt : sps_in_slc[isl]) {
      for(unsigned short xyz = 0; xyz < 3; ++xyz) center[xyz] += spt->XYZ()[xyz];
    } // spt
    for(unsigned short xyz = 0; xyz < 3; ++xyz) center[xyz] /= sum;
    // do a linear fit
    double sumx = 0, sumy = 0., sumz = 0., sumx2 = 0, sumy2 = 0, sumz2 = 0;
    double sumxy = 0, sumxz = 0;
    for(auto& spt : sps_in_slc[isl]) {
      double xx = spt->XYZ()[0] - center[0];
      double yy = spt->XYZ()[1] - center[1];
      double zz = spt->XYZ()[2] - center[2];
      sumx += xx;
      sumy += yy;
      sumz += zz;
      sumx2 += xx * xx;
      sumy2 += yy * yy;
      sumz2 += zz * zz;
      sumxy += xx * yy;
      sumxz += xx * zz;
    } // spt
    double delta = sum * sumx2 - sumx * sumx;
    if(delta <= 0) continue;
    // calculate the slopes
    double dydx = (sumxy * sum  - sumx * sumy) / delta;
    double dzdx = (sumxz * sum  - sumx * sumz) / delta;
    // and convert to direction vector
    double direction[3];
    double norm = std::sqrt(1 + dydx * dydx + dzdx * dzdx);
    direction[0] = 1 / norm;
    direction[1] = dydx / norm;
    direction[2] = dzdx / norm;
    // Find the points that are farthest away from the center along the slice axis
    unsigned int imax = 0, imin = 0;
    float minAlong = 1E6, maxAlong = -1E6;
    double tmpVec[3];
    for(unsigned int ipt = 0; ipt < sps_in_slc[isl].size(); ++ipt) {
      auto& spt = sps_in_slc[isl][ipt];
      for(unsigned short xyz = 0; xyz < 3; ++xyz) tmpVec[xyz] = spt->XYZ()[xyz] - center[xyz];
      double dotp = 0;
      for(unsigned short xyz = 0; xyz < 3; ++xyz) dotp += tmpVec[xyz] * direction[xyz];
      if(dotp < minAlong) { minAlong = dotp; imin = ipt; }
      if(dotp > maxAlong) { maxAlong = dotp; imax = ipt; }
    } // spt
    // Find the aspect ratio
    float cnt = 0.;
    float aspectRatio = 0;
    double arg = sum * sumy2 - sumy * sumy;
    if(arg > 0) {
      aspectRatio += std::abs(sum * sumxy - sumx * sumy) / sqrt(delta * arg);
      ++cnt;
    }
    arg = sum * sumz2 - sumz * sumz;
    if(arg > 0) {
      aspectRatio += std::abs(sum * sumxz - sumx * sumz) / sqrt(delta * arg);
      ++cnt;
    }
    if(cnt > 1) aspectRatio /= cnt;
    int id = isl + 1;
    Point_t ctr(center[0], center[1], center[2]);
    Vector_t dir(direction[0], direction[1], direction[2]);
    auto pos0 = sps_in_slc[isl][imin]->XYZ();
    Point_t ep0(pos0[0], pos0[1], pos0[2]);
    auto pos1 = sps_in_slc[isl][imax]->XYZ();
    Point_t ep1(pos1[0], pos1[1], pos1[2]);
    slcCol.emplace_back(id, ctr, dir, ep0, ep1, aspectRatio, charge);
    util::CreateAssn(*this, evt, slcCol, slcHits[isl], *slc_hit_assn);
    util::CreateAssn(*this, evt, slcCol, sps_in_slc[isl], *slc_sps_assn);
  } // isl
  
  std::unique_ptr<std::vector<recob::Slice> > scol(new std::vector<recob::Slice>(std::move(slcCol)));
  evt.put(std::move(scol));
  evt.put(std::move(slc_hit_assn));
  evt.put(std::move(slc_sps_assn));
}

DEFINE_ART_MODULE(cluster::DBCluster3D)