DirOfGamma.cxx

Go to the documentation of this file.

#include "DirOfGamma.h"

#include "larreco/RecoAlg/PMAlg/Utilities.h"

#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom<>()
#include "larcore/Geometry/WireReadout.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "lardataalg/DetectorInfo/DetectorPropertiesData.h"

#include "messagefacility/MessageLogger/MessageLogger.h"

#include "TMath.h"

ems::Hit2D::Hit2D(detinfo::DetectorPropertiesData const& detProp, art::Ptr<recob::Hit> src)
  : fHit(src)
{
  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
  auto const& wireID = src->WireID();

  auto const wireCenter = wireReadoutGeom.Wire(wireID).GetCenter();
  double const x = detProp.ConvertTicksToX(src->PeakTime(), wireID);

  double const globalWire = wireReadoutGeom.Plane({wireID.Cryostat, wireID.TPC % 2, wireID.Plane})
                              .WireCoordinate(wireCenter);
  fPoint.Set(globalWire, x);
  fCharge = src->ROISummedADC();
}

ems::Bin2D::Bin2D(const TVector2& center) : fCenter2D(center), fTotCharge(0.0), fSize(0) {}

void ems::Bin2D::Add(Hit2D* hit)
{
  fHits2D.push_back(hit);
  fTotCharge += hit->GetCharge();
  fSize = fHits2D.size();
  SortLess();
}

void ems::Bin2D::Sort()
{
  return std::sort(fHits2D.begin(), fHits2D.end(), bDistCentMore2D(fCenter2D));
}

void ems::Bin2D::SortLess()
{
  return std::sort(fHits2D.begin(), fHits2D.end(), bDistCentLess2D(fCenter2D));
}

std::vector<art::Ptr<recob::Hit>> ems::Bin2D::GetIniHits(const double radius,
                                                         const unsigned int nhits) const
{

  std::vector<art::Ptr<recob::Hit>> vec;
  for (unsigned int i = 0; i < fHits2D.size(); i++) {
    if (pma::Dist2(fHits2D[i]->GetPointCm(), fCenter2D) < radius * radius) {
      vec.push_back(fHits2D[i]->GetHitPtr());
      if (vec.size() == nhits) break;
    }
  }

  return vec;
}

ems::EndPoint::EndPoint(const Hit2D& center, const std::vector<Hit2D*>& hits, unsigned int nbins)
  : fCenter2D(center), fPoints2D(hits), fNbins(nbins)
{

  for (unsigned int i = 0; i < fNbins; i++) {
    fBins.push_back(Bin2D(center.GetPointCm()));
  }

  FillBins();
  ComputeMaxCharge();
  ComputeMeanCharge();

  fPlane = center.GetHitPtr()->WireID().Plane;
  fTpc = center.GetHitPtr()->WireID().TPC;
  fCryo = center.GetHitPtr()->WireID().Cryostat;
}

void ems::EndPoint::FillBins()
{
  TVector2 vstart(0, 1);

  unsigned int saveid = 0;
  bool exist = false;
  for (unsigned int i = 0; i < fPoints2D.size(); i++) {
    if (fPoints2D[i]->GetHitPtr().key() != fCenter2D.GetHitPtr().key()) {
      TVector2 pos(fPoints2D[i]->GetPointCm());
      TVector2 centre(fCenter2D.GetPointCm());
      TVector2 vecp = pos - centre;
      float sinsign = (vstart.X() * vecp.Y()) - (vstart.Y() * vecp.X());
      float cosine = (vstart * vecp) / vecp.Mod();
      float theta = 180.0F * (std::acos(cosine)) / TMath::Pi();

      unsigned int id = 0;
      double bin = double(360.0) / double(fNbins);

      if (sinsign >= 0)
        id = int(theta / bin);
      else if (sinsign < 0)
        id = int(theta / bin) + (fNbins / 2);
      if (id > (fNbins - 1)) id = (fNbins - 1);

      fBins[id].Add(fPoints2D[i]);
      fBins[(id + 1) % fNbins].Add(fPoints2D[i]);
    }
    else {
      saveid = i;
      exist = true;
    }
  }

  if (exist)
    for (unsigned int id = 0; id < fNbins; id++)
      fBins[id].Add(fPoints2D[saveid]);
}

void ems::EndPoint::ComputeMaxCharge()
{
  fMaxCharge = 0.0;
  unsigned int saveid = 0;
  for (unsigned int i = 0; i < fNbins; i++)
    if (fBins[i].Size() && (fMaxCharge < fBins[i].GetTotCharge())) {
      fMaxCharge = fBins[i].GetTotCharge();
      saveid = i;
    }

  fMaxChargeIdBin = saveid;
}

void ems::EndPoint::ComputeMeanCharge()
{
  fMeanCharge = 0.0;
  if (fNbins == 0) return;

  unsigned int iprev, inext;

  if (fMaxChargeIdBin > 0)
    iprev = fMaxChargeIdBin - 1;
  else
    iprev = fNbins - 1;

  inext = (fMaxChargeIdBin + 1) % fNbins;

  double sumcharge = 0.0;
  for (unsigned int i = 0; i < fNbins; i++)
    if ((i != fMaxChargeIdBin) && (i != iprev) && (i != inext))
      sumcharge += fBins[i].GetTotCharge();

  fMeanCharge = sumcharge / double(fNbins);
}

double ems::EndPoint::GetAsymmetry() const
{
  if ((fMaxCharge + fMeanCharge) == 0) return 0.0;
  return ((fMaxCharge - fMeanCharge) / (fMaxCharge + fMeanCharge));
}

ems::DirOfGamma::DirOfGamma(detinfo::DetectorPropertiesData const& detProp,
                            const std::vector<art::Ptr<recob::Hit>>& src,
                            unsigned int nbins,
                            unsigned int idcl)
  : fNbins(nbins), fIdCl(idcl), fCandidateID(0), fIsCandidateIDset(false)
{
  fHits = src;

  for (unsigned int i = 0; i < src.size(); i++) {
    Hit2D* hit = new Hit2D(detProp, src[i]);
    fPoints2D.push_back(hit);
  }

  ComputeBaryCenter();

  for (unsigned int i = 0; i < fNbins; i++)
    fBins.push_back(Bin2D(fBaryCenter));

  FillBins();
  ComputeMaxDist();
  if (FindCandidates()) {
    ComputeMaxCharge();
    FindInitialPart();
  }
}

void ems::DirOfGamma::ComputeBaryCenter()
{
  double nomx = 0.0;
  double nomy = 0.0;
  double denom = 0.0;
  for (unsigned int i = 0; i < fPoints2D.size(); i++) {
    nomx += fPoints2D[i]->GetPointCm().X() * fPoints2D[i]->GetCharge();
    nomy += fPoints2D[i]->GetPointCm().Y() * fPoints2D[i]->GetCharge();
    denom += fPoints2D[i]->GetCharge();
  }

  double bx = nomx / denom;
  double by = nomy / denom;
  fBaryCenter.Set(bx, by);
}

void ems::DirOfGamma::FillBins()
{
  TVector2 vstart(0, 1);

  for (unsigned int i = 0; i < fPoints2D.size(); i++) {
    TVector2 pos(fPoints2D[i]->GetPointCm().X(), fPoints2D[i]->GetPointCm().Y());
    TVector2 vecp = pos - fBaryCenter;
    float sinsign = (vstart.X() * vecp.Y()) - (vstart.Y() * vecp.X());
    float cosine = (vstart * vecp) / (vstart.Mod() * vecp.Mod());
    float theta = 180.0F * (std::acos(cosine)) / TMath::Pi();

    unsigned int id = 0;
    double bin = double(360.0) / double(fNbins);

    if (sinsign >= 0)
      id = int(theta / bin);
    else if (sinsign < 0)
      id = int(theta / bin) + (fNbins / 2);
    if (id > (fNbins - 1)) id = (fNbins - 1);

    fBins[id].Add(fPoints2D[i]);
  }

  for (unsigned int id = 0; id < fBins.size(); id++)
    fBins[id].Sort();
}

void ems::DirOfGamma::ComputeMaxDist()
{
  double maxdist2 = 0.0;

  for (unsigned int id = 0; id < fBins.size(); id++) {

    if (!fBins[id].Size()) continue;

    Hit2D* candidate = fBins[id].GetHits2D().front();
    if (candidate) {
      double dist2 = pma::Dist2(candidate->GetPointCm(), fBaryCenter);
      if (dist2 > maxdist2) { maxdist2 = dist2; }
    }
  }

  fNormDist = std::sqrt(maxdist2);
}

bool ems::DirOfGamma::FindCandidates()
{
  float rad = 0.5F * fNormDist;
  unsigned int nbins = fNbins * 4;
  for (unsigned int id = 0; id < fNbins; id++) {

    if (!fBins[id].Size()) continue;

    std::vector<Hit2D*> points;
    Hit2D* candidate2D = fBins[id].GetHits2D().front();

    for (unsigned int i = 0; i < fPoints2D.size(); i++) {
      double distnorm = std::sqrt(pma::Dist2(candidate2D->GetPointCm(), fBaryCenter)) / fNormDist;
      double dist2 = pma::Dist2(candidate2D->GetPointCm(), fPoints2D[i]->GetPointCm());

      if ((distnorm > 0.5) && (dist2 < rad * rad)) points.push_back(fPoints2D[i]);
    }

    if (fBins[id].Size() > 1) {
      EndPoint ep(*candidate2D, points, nbins);
      fCandidates.push_back(ep);
    }
  }
  if (fCandidates.size())
    return true;
  else
    return false;
}

void ems::DirOfGamma::ComputeMaxCharge()
{
  fNormCharge = 0.0;
  for (unsigned int i = 0; i < fCandidates.size(); i++) {
    if (fCandidates[i].GetMaxCharge() > fNormCharge) {
      fNormCharge = fCandidates[i].GetMaxCharge();
    }
  }
}

void ems::DirOfGamma::FindInitialPart()
{
  double max_asymmetry = 0.0;
  unsigned int saveid = 0;
  bool found = false;

  double maxdist2 = 0.0;
  double maxcharge = 0.0;
  unsigned int idmaxdist = 0;
  unsigned int idmaxcharge = 0;

  for (unsigned int i = 0; i < fCandidates.size(); i++) {
    double dist2 = pma::Dist2(fCandidates[i].GetPosition(), fBaryCenter);
    double charge = fCandidates[i].GetMaxCharge();
    if (dist2 > maxdist2) {
      maxdist2 = dist2;
      idmaxdist = i;
    }
    if (charge > maxcharge) {
      maxcharge = charge;
      idmaxcharge = i;
    }
  }

  maxdist2 = 0.0;
  unsigned int idmaxdistb = 0;
  for (size_t i = 0; i < fCandidates.size(); i++) {
    if ((i == idmaxdist) || (i == idmaxcharge)) continue;

    double dist2 = pma::Dist2(fCandidates[i].GetPosition(), fCandidates[idmaxdist].GetPosition());
    if (dist2 > maxdist2) {
      maxdist2 = dist2;
      idmaxdistb = i;
    }
  }

  if (fCandidates.size() > 2) {
    for (unsigned int i = 0; i < fCandidates.size(); i++) {
      double asymmetry = fCandidates[i].GetAsymmetry();

      if ((i == idmaxdist) || (i == idmaxcharge) || (i == idmaxdistb)) {
        if (asymmetry > max_asymmetry) {
          max_asymmetry = asymmetry;
          saveid = i;
          found = true;
        }
      }
    }
  }
  else {
    for (unsigned int i = 0; i < fCandidates.size(); i++) {
      double asymmetry = fCandidates[i].GetAsymmetry();

      if ((i == idmaxdist) || (i == idmaxdistb)) {
        if (asymmetry > max_asymmetry) {
          max_asymmetry = asymmetry;
          saveid = i;
          found = true;
        }
      }
    }
  }

  if (!found)
    mf::LogError("DirOfGamma") << fCandidates.size() << "DirOfGamma - Find Initial Part problem.";

  fStartHit = fCandidates[saveid].GetHit();
  fStartPoint = fCandidates[saveid].GetPosition();
  fIniHits = fCandidates[saveid].MaxChargeBin().GetIniHits();
  fCandidateID = saveid;
}