LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
cluster::HoughBaseAlg Class Reference

#include "HoughBaseAlg.h"

Classes

struct  ChargeInfo_t
 Data structure collecting charge information to be filled in cluster. More...
 

Public Member Functions

 HoughBaseAlg (fhicl::ParameterSet const &pset)
 
virtual ~HoughBaseAlg ()=default
 
size_t FastTransform (const std::vector< art::Ptr< recob::Cluster >> &clusIn, std::vector< recob::Cluster > &ccol, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, CLHEP::HepRandomEngine &engine, art::Event const &evt, std::string const &label)
 
size_t Transform (const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, CLHEP::HepRandomEngine &engine, std::vector< unsigned int > *fpointId_to_clusterId, unsigned int clusterId, unsigned int *nClusters, std::vector< protoTrack > *protoTracks)
 
size_t FastTransform (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &clusIn, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, CLHEP::HepRandomEngine &engine)
 
size_t FastTransform (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &clusIn, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, CLHEP::HepRandomEngine &engine, std::vector< double > &slope, std::vector< ChargeInfo_t > &totalQ)
 
size_t Transform (std::vector< art::Ptr< recob::Hit >> const &hits)
 
size_t Transform (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, double &slope, double &intercept)
 

Private Member Functions

void HLSSaveBMPFile (char const *, unsigned char *, int, int)
 

Private Attributes

int fMaxLines
 Max number of lines that can be found. More...
 
int fMinHits
 
int fSaveAccumulator
 Save bitmap image of accumulator for debugging? More...
 
int fNumAngleCells
 that fall into the "correct" bin will be small and consistent with noise. More...
 
float fMaxDistance
 Max distance that a hit can be from a line to be considered part of that line. More...
 
float fMaxSlope
 Max slope a line can have. More...
 
int fRhoZeroOutRange
 Range in rho over which to zero out area around previously found lines in the accumulator. More...
 
int fThetaZeroOutRange
 Range in theta over which to zero out area around previously found lines in the accumulator. More...
 
float fRhoResolutionFactor
 Factor determining the resolution in rho. More...
 
int fPerCluster
 
int fMissedHits
 
float fMissedHitsDistance
 Distance between hits in a hough line before a hit is considered missed. More...
 
float fMissedHitsToLineSize
 Ratio of missed hits to line size for a line to be considered a fake. More...
 

Friends

class HoughTransformClus
 

Detailed Description

Definition at line 455 of file HoughBaseAlg.h.

Constructor & Destructor Documentation

cluster::HoughBaseAlg::HoughBaseAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 251 of file HoughBaseAlg.cxx.

References util::abs(), cluster::HoughTransform::AddPointReturnMax(), geo::PlaneID::asPlaneID(), geo::TPCID::asTPCID(), detinfo::DetectorPropertiesData::DriftVelocity(), e, detinfo::DetectorPropertiesData::Efield(), fhicl::ParameterSet::get(), cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), hits(), cluster::HoughTransform::Init(), protoTrack::Init(), geo::GeometryCore::Iterate(), geo::kInduction, MF_LOG_DEBUG, geo::GeometryCore::Nplanes(), detinfo::DetectorPropertiesData::NumberTimeSamples(), geo::GeometryCore::Plane(), detinfo::sampling_rate(), cluster::HoughTransform::SetCell(), geo::GeometryCore::SignalType(), sqr(), cluster::HoughTransform::SubtractPoint(), detinfo::DetectorPropertiesData::Temperature(), Transform(), lar::dump::vector(), w, x, and y.

252 {
253  fMaxLines = pset.get<int>("MaxLines");
254  fMinHits = pset.get<int>("MinHits");
255  fSaveAccumulator = pset.get<int>("SaveAccumulator");
256  fNumAngleCells = pset.get<int>("NumAngleCells");
257  fMaxDistance = pset.get<float>("MaxDistance");
258  fMaxSlope = pset.get<float>("MaxSlope");
259  fRhoZeroOutRange = pset.get<int>("RhoZeroOutRange");
260  fThetaZeroOutRange = pset.get<int>("ThetaZeroOutRange");
261  fRhoResolutionFactor = pset.get<float>("RhoResolutionFactor");
262  fPerCluster = pset.get<int>("HitsPerCluster");
263  fMissedHits = pset.get<int>("MissedHits");
264  fMissedHitsDistance = pset.get<float>("MissedHitsDistance");
265  fMissedHitsToLineSize = pset.get<float>("MissedHitsToLineSize");
266 }
float fRhoResolutionFactor
Factor determining the resolution in rho.
Definition: HoughBaseAlg.h:534
int fMaxLines
Max number of lines that can be found.
Definition: HoughBaseAlg.h:519
float fMissedHitsDistance
Distance between hits in a hough line before a hit is considered missed.
Definition: HoughBaseAlg.h:542
int fNumAngleCells
that fall into the "correct" bin will be small and consistent with noise.
Definition: HoughBaseAlg.h:523
int fSaveAccumulator
Save bitmap image of accumulator for debugging?
Definition: HoughBaseAlg.h:522
float fMaxDistance
Max distance that a hit can be from a line to be considered part of that line.
Definition: HoughBaseAlg.h:528
float fMaxSlope
Max slope a line can have.
Definition: HoughBaseAlg.h:529
int fRhoZeroOutRange
Range in rho over which to zero out area around previously found lines in the accumulator.
Definition: HoughBaseAlg.h:531
int fThetaZeroOutRange
Range in theta over which to zero out area around previously found lines in the accumulator.
Definition: HoughBaseAlg.h:533
float fMissedHitsToLineSize
Ratio of missed hits to line size for a line to be considered a fake.
Definition: HoughBaseAlg.h:544
virtual cluster::HoughBaseAlg::~HoughBaseAlg ( )
virtualdefault

Member Function Documentation

size_t cluster::HoughBaseAlg::FastTransform ( const std::vector< art::Ptr< recob::Cluster >> &  clusIn,
std::vector< recob::Cluster > &  ccol,
std::vector< art::PtrVector< recob::Hit >> &  clusHitsOut,
CLHEP::HepRandomEngine &  engine,
art::Event const &  evt,
std::string const &  label 
)

Definition at line 866 of file HoughBaseAlg.cxx.

References art::PtrVector< T >::end(), hits(), cluster::ClusterParamsImportWrapper< Algo >::ImportHits(), MF_LOG_DEBUG, recob::Hit::PeakTime(), geo::WireID::planeID(), recob::Cluster::Sentry, recob::Hit::SigmaPeakTime(), sw, recob::Hit::View(), geo::GeometryCore::Views(), geo::WireID::Wire, recob::Hit::WireID(), and xx.

Referenced by cluster::HoughLineFinder::produce().

872 {
873  std::vector<int> skip;
874 
875  art::FindManyP<recob::Hit> fmh(clusIn, evt, label);
876 
877  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
878  HoughTransform c;
879 
880  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
881  auto const detProp =
883  util::GeometryUtilities const gser{*geom, clockData, detProp};
884 
885  // prepare the algorithm to compute the cluster characteristics;
886  // we use the "standard" one here; configuration would happen here,
887  // but we are using the default configuration for that algorithm
888  ClusterParamsImportWrapper<StandardClusterParamsAlg> ClusterParamAlgo;
889 
890  std::vector<art::Ptr<recob::Hit>> hit;
891 
892  for (auto view : geom->Views()) {
893 
894  MF_LOG_DEBUG("HoughBaseAlg") << "Analyzing view " << view;
895 
896  art::PtrVector<recob::Cluster>::const_iterator clusterIter = clusIn.begin();
897  int clusterID = 0; //the unique ID of the cluster
898 
899  size_t cinctr = 0;
900  while (clusterIter != clusIn.end()) {
901 
902  MF_LOG_DEBUG("HoughBaseAlg")
903  << "Analyzing cinctr=" << cinctr << " which is in view " << (*clusterIter)->View();
904 
905  hit.clear();
906  if (fPerCluster) {
907  if ((*clusterIter)->View() == view) hit = fmh.at(cinctr);
908  }
909  else {
910  while (clusterIter != clusIn.end()) {
911  if ((*clusterIter)->View() == view) {
912 
913  hit = fmh.at(cinctr);
914  } // end if cluster is in correct view
915  clusterIter++;
916  ++cinctr;
917  } //end loop over clusters
918  } //end if not fPerCluster
919  if (hit.size() == 0) {
920  if (fPerCluster) {
921  clusterIter++;
922  ++cinctr;
923  }
924  continue;
925  }
926 
927  MF_LOG_DEBUG("HoughBaseAlg") << "We have all the hits..." << hit.size();
928 
929  std::vector<double> slopevec;
930  std::vector<ChargeInfo_t> totalQvec;
931  std::vector<art::PtrVector<recob::Hit>> planeClusHitsOut;
932  this->FastTransform(clockData, detProp, hit, planeClusHitsOut, engine, slopevec, totalQvec);
933 
934  MF_LOG_DEBUG("HoughBaseAlg") << "Made it through FastTransform" << planeClusHitsOut.size();
935 
936  for (size_t xx = 0; xx < planeClusHitsOut.size(); ++xx) {
937  auto const& hits = planeClusHitsOut.at(xx);
938  recob::Hit const& FirstHit = *hits.front();
939  recob::Hit const& LastHit = *hits.back();
940  const unsigned int sw = FirstHit.WireID().Wire;
941  const unsigned int ew = LastHit.WireID().Wire;
942 
943  // feed the algorithm with all the cluster hits
944  ClusterParamAlgo.ImportHits(gser, hits);
945 
946  // create the recob::Cluster directly in the vector;
947  // NOTE usually we would use cluster::ClusterCreator to save some typing
948  // and some mistakes. In this case, we don't want to pull in the
949  // dependency on ClusterFinder, where ClusterCreator currently lives
950  ccol.emplace_back(float(sw), // start_wire
951  0., // sigma_start_wire
952  FirstHit.PeakTime(), // start_tick
953  FirstHit.SigmaPeakTime(), // sigma_start_tick
954  ClusterParamAlgo.StartCharge(gser).value(),
955  ClusterParamAlgo.StartAngle().value(),
956  ClusterParamAlgo.StartOpeningAngle().value(),
957  float(ew), // end_wire
958  0., // sigma_end_wire,
959  LastHit.PeakTime(), // end_tick
960  LastHit.SigmaPeakTime(), // sigma_end_tick
961  ClusterParamAlgo.EndCharge(gser).value(),
962  ClusterParamAlgo.EndAngle().value(),
963  ClusterParamAlgo.EndOpeningAngle().value(),
964  ClusterParamAlgo.Integral().value(),
965  ClusterParamAlgo.IntegralStdDev().value(),
966  ClusterParamAlgo.SummedADC().value(),
967  ClusterParamAlgo.SummedADCStdDev().value(),
968  ClusterParamAlgo.NHits(),
969  ClusterParamAlgo.MultipleHitDensity(),
970  ClusterParamAlgo.Width(gser),
971  clusterID,
972  FirstHit.View(),
973  FirstHit.WireID().planeID(),
975 
976  ++clusterID;
977  clusHitsOut.push_back(planeClusHitsOut.at(xx));
978  }
979 
980  hit.clear();
981  if (clusterIter != clusIn.end()) {
982  clusterIter++;
983  ++cinctr;
984  }
985  // listofxmax.clear();
986  // listofymax.clear();
987  } // end loop over clusters
988 
989  } // end loop over views
990 
991  return ccol.size();
992 }
Double_t xx
Definition: macro.C:12
size_t FastTransform(const std::vector< art::Ptr< recob::Cluster >> &clusIn, std::vector< recob::Cluster > &ccol, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, CLHEP::HepRandomEngine &engine, art::Event const &evt, std::string const &label)
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:276
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
void hits()
Definition: readHits.C:15
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
Detector simulation of raw signals on wires.
Float_t sw
Definition: plot.C:20
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
#define MF_LOG_DEBUG(id)
std::set< View_t > const & Views() const
Returns a list of possible views in the detector.
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:224
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
constexpr PlaneID const & planeID() const
Definition: geo_types.h:620
TCEvent evt
Definition: DataStructs.cxx:8
size_t cluster::HoughBaseAlg::FastTransform ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> const &  clusIn,
std::vector< art::PtrVector< recob::Hit >> &  clusHitsOut,
CLHEP::HepRandomEngine &  engine 
)

Definition at line 995 of file HoughBaseAlg.cxx.

1000 {
1001  std::vector<double> slopevec;
1002  std::vector<ChargeInfo_t> totalQvec;
1003  return FastTransform(clockData, detProp, clusIn, clusHitsOut, engine, slopevec, totalQvec);
1004 }
size_t FastTransform(const std::vector< art::Ptr< recob::Cluster >> &clusIn, std::vector< recob::Cluster > &ccol, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, CLHEP::HepRandomEngine &engine, art::Event const &evt, std::string const &label)
size_t cluster::HoughBaseAlg::FastTransform ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> const &  clusIn,
std::vector< art::PtrVector< recob::Hit >> &  clusHitsOut,
CLHEP::HepRandomEngine &  engine,
std::vector< double > &  slope,
std::vector< ChargeInfo_t > &  totalQ 
)
Todo:
explain where the 0.001 comes from
Todo:
: the collection plane's characteristic hit width's are,
Todo:
: on average, about 5 time samples wider than the induction plane's.
Todo:
: this is hard-coded for now.
Todo:
should it really be wire_pitch[0] in the if statements below, or the pitch for the plane of the hit?

Definition at line 1007 of file HoughBaseAlg.cxx.

References util::abs(), lar::util::StatCollector< T, W >::add(), cluster::HoughTransform::AddPointReturnMax(), geo::PlaneID::asPlaneID(), geo::TPCID::asTPCID(), art::PtrVector< T >::back(), art::PtrVector< T >::clear(), detinfo::DetectorPropertiesData::DriftVelocity(), detinfo::DetectorPropertiesData::Efield(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), cluster::HoughTransform::Init(), geo::GeometryCore::Iterate(), geo::kInduction, MF_LOG_DEBUG, geo::GeometryCore::Nplanes(), geo::GeometryCore::Plane(), art::PtrVector< T >::push_back(), detinfo::DetectorPropertiesData::ReadOutWindowSize(), lar::util::StatCollector< T, W >::RMS(), detinfo::sampling_rate(), cluster::HoughTransform::SetCell(), geo::GeometryCore::SignalType(), lar::util::StatCollector< T, W >::Sum(), detinfo::DetectorPropertiesData::Temperature(), lar::dump::vector(), x, and y.

1014 {
1015  std::vector<int> skip;
1016 
1017  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
1018  lariov::ChannelStatusProvider const* channelStatus =
1019  lar::providerFrom<lariov::ChannelStatusService>();
1020 
1021  CLHEP::RandFlat flat(engine);
1022 
1023  std::vector<art::Ptr<recob::Hit>> hit;
1024 
1025  size_t cinctr = 0;
1026  hit.clear();
1027  for (std::vector<art::Ptr<recob::Hit>>::const_iterator ii = clusIn.begin(); ii != clusIn.end();
1028  ii++)
1029  hit.push_back((*ii)); // this is new
1030 
1031  if (hit.size() == 0) {
1032  if (fPerCluster) { ++cinctr; }
1033  return -1;
1034  }
1035 
1036  geo::WireID const& wireid = hit.at(0)->WireID();
1037  auto const& tpcid = wireid.asPlaneID().asTPCID();
1038 
1039  geo::SigType_t sigt = geom->SignalType(wireid);
1040 
1041  if (hit.size() == 0) {
1042  if (fPerCluster) { ++cinctr; }
1043  return -1; // continue;
1044  }
1045 
1046  auto const num_planes = geom->Nplanes(tpcid);
1047  std::vector<double> wire_pitch(num_planes, 0.);
1048  for (auto const& plane : geom->Iterate<geo::PlaneGeo>(tpcid)) {
1049  auto const p = plane.ID().Plane;
1050  wire_pitch[p] = plane.WirePitch();
1051  }
1052 
1053  //factor to make x and y scale the same units
1054  std::vector<double> xyScale(num_planes, 0.);
1055 
1057  double driftVelFactor = 0.001 * detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
1058 
1059  for (size_t p = 0; p < xyScale.size(); ++p)
1060  xyScale[p] = driftVelFactor * sampling_rate(clockData) / wire_pitch[p];
1061 
1062  int x = 0, y = 0;
1063  int dx = geom->Plane(wireid.asPlaneID()).Nwires(); // number of wires
1064  const int dy = detProp.ReadOutWindowSize(); // number of time samples.
1065  skip.clear();
1066  skip.resize(hit.size());
1067  std::vector<int> listofxmax;
1068  std::vector<int> listofymax;
1069  std::vector<int> hitTemp; //indecies ofcandidate hits
1070  std::vector<int> sequenceHolder; //channels of hits in list
1071  std::vector<int> currentHits; //working vector of hits
1072  std::vector<int> lastHits; //best list of hits
1073  art::PtrVector<recob::Hit> clusterHits;
1074  float indcolscaling = 0.; //a parameter to account for the different
1075  //characteristic hit width of induction and collection plane
1076  float centerofmassx = 0;
1077  float centerofmassy = 0;
1078  float denom = 0;
1079  float intercept = 0.;
1080  float slope = 0.;
1081  //this array keeps track of the hits that have already been associated with a line.
1082  int xMax = 0;
1083  int yMax = 0;
1084  float rho;
1085  float theta;
1086  int accDx(0), accDy(0);
1087 
1088  HoughTransform c;
1089  //Init specifies the size of the two-dimensional accumulator
1090  //(based on the arguments, number of wires and number of time samples).
1091  //adds all of the hits (that have not yet been associated with a line) to the accumulator
1092  c.Init(dx, dy, fRhoResolutionFactor, fNumAngleCells);
1093 
1094  // count is how many points are left to randomly insert
1095  unsigned int count = hit.size();
1096  std::vector<unsigned int> accumPoints;
1097  accumPoints.resize(hit.size());
1098  int nAccum = 0;
1099  unsigned int nLinesFound = 0;
1100 
1101  for (; count > 0; count--) {
1102 
1103  // The random hit we are examining
1104  unsigned int randInd = (unsigned int)(flat.fire() * hit.size());
1105 
1106  MF_LOG_DEBUG("HoughBaseAlg") << "randInd=" << randInd << " and size is " << hit.size();
1107 
1108  // Skip if it's already in a line
1109  if (skip.at(randInd) == 1) continue;
1110 
1111  // If we have already accumulated the point, skip it
1112  if (accumPoints.at(randInd)) continue;
1113  accumPoints.at(randInd) = 1;
1114 
1115  // zeroes out the neighborhood of all previous lines
1116  for (unsigned int i = 0; i < listofxmax.size(); ++i) {
1117  int yClearStart = listofymax.at(i) - fRhoZeroOutRange;
1118  if (yClearStart < 0) yClearStart = 0;
1119 
1120  int yClearEnd = listofymax.at(i) + fRhoZeroOutRange;
1121  if (yClearEnd >= accDy) yClearEnd = accDy - 1;
1122 
1123  int xClearStart = listofxmax.at(i) - fThetaZeroOutRange;
1124  if (xClearStart < 0) xClearStart = 0;
1125 
1126  int xClearEnd = listofxmax.at(i) + fThetaZeroOutRange;
1127  if (xClearEnd >= accDx) xClearEnd = accDx - 1;
1128 
1129  for (y = yClearStart; y <= yClearEnd; ++y) {
1130  for (x = xClearStart; x <= xClearEnd; ++x) {
1131  c.SetCell(y, x, 0);
1132  }
1133  }
1134  } // end loop over size of listxmax
1135 
1136  //find the weightiest cell in the accumulator.
1137  int maxCell = 0;
1138  unsigned int wireMax = hit.at(randInd)->WireID().Wire;
1139 
1140  // Add the randomly selected point to the accumulator
1141  std::array<int, 3> max = c.AddPointReturnMax(wireMax, (int)(hit.at(randInd)->PeakTime()));
1142  maxCell = max.at(0);
1143  xMax = max.at(1);
1144  yMax = max.at(2);
1145  nAccum++;
1146 
1147  // Continue if the biggest maximum for the randomly selected point is smaller than fMinHits
1148  if (maxCell < fMinHits) continue;
1149 
1150  // Find the center of mass of the 3x3 cell system (with maxCell at the center).
1151  denom = centerofmassx = centerofmassy = 0;
1152 
1153  if (xMax > 0 && yMax > 0 && xMax + 1 < accDx && yMax + 1 < accDy) {
1154  for (int i = -1; i < 2; ++i) {
1155  for (int j = -1; j < 2; ++j) {
1156  int cell = c.GetCell(yMax + i, xMax + j);
1157  denom += cell;
1158  centerofmassx += j * cell;
1159  centerofmassy += i * cell;
1160  }
1161  }
1162  centerofmassx /= denom;
1163  centerofmassy /= denom;
1164  }
1165  else
1166  centerofmassx = centerofmassy = 0;
1167 
1168  //fill the list of cells that have already been found
1169  listofxmax.push_back(xMax);
1170  listofymax.push_back(yMax);
1171 
1172  // Find the lines equation
1173  c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
1174  MF_LOG_DEBUG("HoughBaseAlg") << "Transform(I) found maximum at (d=" << rho << " a=" << theta
1175  << ")"
1176  " from absolute maximum "
1177  << c.GetCell(yMax, xMax) << " at (d=" << yMax << ", a=" << xMax
1178  << ")";
1179  slope = -1. / tan(theta);
1180  intercept = (rho / sin(theta));
1181  double distance;
1185  if (sigt == geo::kInduction)
1186  indcolscaling = 5.;
1187  else
1188  indcolscaling = 0.;
1189 
1190  // note this doesn't work for wrapped wire planes!
1191  if (!std::isinf(slope) && !std::isnan(slope)) {
1192  sequenceHolder.clear();
1193  hitTemp.clear();
1194  for (size_t i = 0; i < hit.size(); ++i) {
1195  distance = (TMath::Abs(hit.at(i)->PeakTime() - slope * (double)(hit.at(i)->WireID().Wire) -
1196  intercept) /
1197  (std::sqrt(pow(xyScale[hit.at(i)->WireID().Plane] * slope, 2) + 1)));
1198 
1199  if (distance < fMaxDistance + hit.at(i)->RMS() + indcolscaling && skip.at(i) != 1) {
1200  hitTemp.push_back(i);
1201  sequenceHolder.push_back(hit.at(i)->Channel());
1202  }
1203 
1204  } // end loop over hits
1205 
1206  if (hitTemp.size() < 2) continue;
1207  currentHits.clear();
1208  lastHits.clear();
1209  int j;
1210  currentHits.push_back(0);
1211  for (size_t i = 0; i + 1 < sequenceHolder.size(); ++i) {
1212  j = 1;
1213  while ((channelStatus->IsBad(sequenceHolder.at(i) + j)) == true)
1214  j++;
1215  if (sequenceHolder.at(i + 1) - sequenceHolder.at(i) <= j + fMissedHits)
1216  currentHits.push_back(i + 1);
1217  else if (currentHits.size() > lastHits.size()) {
1218  lastHits = currentHits;
1219  currentHits.clear();
1220  }
1221  else
1222  currentHits.clear();
1223  }
1224 
1225  if (currentHits.size() > lastHits.size()) lastHits = currentHits;
1226 
1227  // Check if lastHits has hits with big gaps in it
1228  // lastHits[i] is ordered in increasing channel and then increasing peak time,
1229  // as a consequence, if the line has a negative slope and there are multiple hits in the line for a channel,
1230  // we have to go back to the first hit (in terms of lastHits[i]) of that channel to find the distance
1231  // between hits
1232  //std::cout << "New line" << std::endl;
1233  double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
1234  tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
1235  int missedHits = 0;
1236  int lastHitsChannel = 0; //lastHits.at(0);
1237  int nHitsPerChannel = 1;
1238 
1239  MF_LOG_DEBUG("HoughBaseAlg") << "filling the pCorner arrays around here..."
1240  << "\n but first, lastHits size is " << lastHits.size()
1241  << " and lastHitsChannel=" << lastHitsChannel;
1242 
1243  double pCorner0[2];
1244  double pCorner1[2];
1245  unsigned int lastChannel = hit.at(hitTemp.at(lastHits.at(0)))->Channel();
1246 
1247  for (size_t i = 0; i < lastHits.size() - 1; ++i) {
1248  bool newChannel = false;
1249  if (slope < 0) {
1250  if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() != lastChannel) {
1251  newChannel = true;
1252  }
1253  if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() == lastChannel) nHitsPerChannel++;
1254  }
1255 
1258 
1259  if (slope > 0 || (!newChannel && nHitsPerChannel <= 1)) {
1260 
1261  pCorner0[0] = (hit.at(hitTemp.at(lastHits.at(i)))->Channel()) * wire_pitch[0];
1262  pCorner0[1] = hit.at(hitTemp.at(lastHits.at(i)))->PeakTime() * tickToDist;
1263  pCorner1[0] = (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel()) * wire_pitch[0];
1264  pCorner1[1] = hit.at(hitTemp.at(lastHits.at(i + 1)))->PeakTime() * tickToDist;
1265  if (std::sqrt(pow(pCorner0[0] - pCorner1[0], 2) + pow(pCorner0[1] - pCorner1[1], 2)) >
1267  missedHits++;
1268  }
1269 
1270  else if (slope < 0 && newChannel && nHitsPerChannel > 1) {
1271 
1272  pCorner0[0] =
1273  (hit.at(hitTemp.at(lastHits.at(lastHitsChannel)))->Channel()) * wire_pitch[0];
1274  pCorner0[1] = hit.at(hitTemp.at(lastHits.at(lastHitsChannel)))->PeakTime() * tickToDist;
1275  pCorner1[0] = (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel()) * wire_pitch[0];
1276  pCorner1[1] = hit.at(hitTemp.at(lastHits.at(i + 1)))->PeakTime() * tickToDist;
1277  if (std::sqrt(pow(pCorner0[0] - pCorner1[0], 2) + pow(pCorner0[1] - pCorner1[1], 2)) >
1279  missedHits++;
1280  lastChannel = hit.at(hitTemp.at(lastHits.at(i)))->Channel();
1281  lastHitsChannel = i + 1;
1282  nHitsPerChannel = 0;
1283  }
1284  }
1285 
1286  if ((double)missedHits / ((double)lastHits.size() - 1) > fMissedHitsToLineSize) continue;
1287 
1288  clusterHits.clear();
1289  if (lastHits.size() < 5) continue;
1290 
1291  // reduce rounding errors by using double (RMS is very sensitive to them)
1292  lar::util::StatCollector<double> integralQ, summedQ;
1293 
1294  for (size_t i = 0; i < lastHits.size(); ++i) {
1295  clusterHits.push_back(hit.at(hitTemp.at(lastHits.at(i))));
1296  integralQ.add(clusterHits.back()->Integral());
1297  summedQ.add(clusterHits.back()->SummedADC());
1298  skip.at(hitTemp.at(lastHits.at(i))) = 1;
1299  }
1300  //protection against very steep uncorrelated hits
1301  if (std::abs(slope) > fMaxSlope) continue;
1302 
1303  clusHitsOut.push_back(clusterHits);
1304  slopevec.push_back(slope);
1305  totalQvec.emplace_back(integralQ.Sum(),
1306  integralQ.RMS(), // TODO biased value; should unbias?
1307  summedQ.Sum(),
1308  summedQ.RMS() // TODO biased value; should unbias?
1309  );
1310 
1311  } // end if !std::isnan
1312 
1313  nLinesFound++;
1314 
1315  if (nLinesFound > (unsigned int)fMaxLines) break;
1316 
1317  } // end loop over hits
1318 
1319  // saves a bitmap image of the accumulator (useful for debugging),
1320  // with scaling based on the maximum cell value
1321  if (fSaveAccumulator) {
1322  //finds the maximum cell in the accumulator for image scaling
1323  int cell, pix = 0, maxCell = 0;
1324  for (y = 0; y < accDy; ++y) {
1325  for (x = 0; x < accDx; ++x) {
1326  cell = c.GetCell(y, x);
1327  if (cell > maxCell) maxCell = cell;
1328  }
1329  }
1330 
1331  std::unique_ptr<unsigned char[]> outPix(new unsigned char[accDx * accDy]);
1332  unsigned int PicIndex = 0;
1333  for (y = 0; y < accDy; ++y) {
1334  for (x = 0; x < accDx; ++x) {
1335  //scales the pixel weights based on the maximum cell value
1336  if (maxCell > 0) pix = (int)((1500 * c.GetCell(y, x)) / maxCell);
1337  outPix[PicIndex++] = pix;
1338  }
1339  }
1340 
1341  HLSSaveBMPFile("houghaccum.bmp", outPix.get(), accDx, accDy);
1342  } // end if saving accumulator
1343 
1344  hit.clear();
1345  lastHits.clear();
1346  listofxmax.clear();
1347  listofymax.clear();
1348  return clusHitsOut.size();
1349 }
Float_t x
Definition: compare.C:6
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
float fRhoResolutionFactor
Factor determining the resolution in rho.
Definition: HoughBaseAlg.h:534
Float_t y
Definition: compare.C:6
constexpr auto abs(T v)
Returns the absolute value of the argument.
Weight_t RMS() const
Returns the root mean square.
int fMaxLines
Max number of lines that can be found.
Definition: HoughBaseAlg.h:519
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
float fMissedHitsDistance
Distance between hits in a hough line before a hit is considered missed.
Definition: HoughBaseAlg.h:542
reference back()
Definition: PtrVector.h:387
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
int fNumAngleCells
that fall into the "correct" bin will be small and consistent with noise.
Definition: HoughBaseAlg.h:523
int fSaveAccumulator
Save bitmap image of accumulator for debugging?
Definition: HoughBaseAlg.h:522
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
Signal from induction planes.
Definition: geo_types.h:151
Weight_t Sum() const
Returns the weighted sum of the values.
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
float fMaxDistance
Max distance that a hit can be from a line to be considered part of that line.
Definition: HoughBaseAlg.h:528
float fMaxSlope
Max slope a line can have.
Definition: HoughBaseAlg.h:529
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
int fRhoZeroOutRange
Range in rho over which to zero out area around previously found lines in the accumulator.
Definition: HoughBaseAlg.h:531
Detector simulation of raw signals on wires.
constexpr TPCID const & asTPCID() const
Conversion to TPCID (for convenience of notation).
Definition: geo_types.h:438
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
#define MF_LOG_DEBUG(id)
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
int fThetaZeroOutRange
Range in theta over which to zero out area around previously found lines in the accumulator.
Definition: HoughBaseAlg.h:533
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
void clear()
Definition: PtrVector.h:533
Collects statistics on a single quantity (weighted)
float fMissedHitsToLineSize
Ratio of missed hits to line size for a line to be considered a fake.
Definition: HoughBaseAlg.h:544
void add(Data_t value, Weight_t weight=Weight_t(1.0))
Adds one entry with specified value and weight.
void HLSSaveBMPFile(char const *, unsigned char *, int, int)
void cluster::HoughBaseAlg::HLSSaveBMPFile ( char const *  ,
unsigned char *  ,
int  ,
int   
)
private

Definition at line 832 of file HoughBaseAlg.cxx.

References util::size().

833 {
834  std::ofstream bmpFile(fileName, std::ios::binary);
835  bmpFile.write("B", 1);
836  bmpFile.write("M", 1);
837  int bitsOffset = 54 + 256 * 4;
838  int size = bitsOffset + dx * dy; //header plus 256 entry LUT plus pixels
839  bmpFile.write((const char*)&size, 4);
840  int reserved = 0;
841  bmpFile.write((const char*)&reserved, 4);
842  bmpFile.write((const char*)&bitsOffset, 4);
843  int bmiSize = 40;
844  bmpFile.write((const char*)&bmiSize, 4);
845  bmpFile.write((const char*)&dx, 4);
846  bmpFile.write((const char*)&dy, 4);
847  short planes = 1;
848  bmpFile.write((const char*)&planes, 2);
849  short bitCount = 8;
850  bmpFile.write((const char*)&bitCount, 2);
851  int i, temp = 0;
852  for (i = 0; i < 6; i++)
853  bmpFile.write((const char*)&temp, 4); // zero out optional color info
854  // write a linear LUT
855  char lutEntry[4]; // blue,green,red
856  lutEntry[3] = 0; // reserved part
857  for (i = 0; i < 256; i++) {
858  lutEntry[0] = lutEntry[1] = lutEntry[2] = i;
859  bmpFile.write(lutEntry, sizeof lutEntry);
860  }
861  // write the actual pixels
862  bmpFile.write((const char*)pix, dx * dy);
863 }
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
size_t cluster::HoughBaseAlg::Transform ( const detinfo::DetectorClocksData clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> const &  hits,
CLHEP::HepRandomEngine &  engine,
std::vector< unsigned int > *  fpointId_to_clusterId,
unsigned int  clusterId,
unsigned int *  nClusters,
std::vector< protoTrack > *  protoTracks 
)

Referenced by HoughBaseAlg().

size_t cluster::HoughBaseAlg::Transform ( std::vector< art::Ptr< recob::Hit >> const &  hits)
size_t cluster::HoughBaseAlg::Transform ( detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> const &  hits,
double &  slope,
double &  intercept 
)
Todo:
could eventually refine this method to throw out hits that are
Todo:
far from the hough line and refine the fit

Definition at line 1352 of file HoughBaseAlg.cxx.

References cluster::HoughTransform::AddPointReturnMax(), cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), cluster::HoughTransform::GetMax(), hits(), cluster::HoughTransform::Init(), MF_LOG_DEBUG, geo::GeometryCore::Nwires(), and detinfo::DetectorPropertiesData::ReadOutWindowSize().

1356 {
1357  HoughTransform c;
1358 
1360 
1361  int dx = geom->Nwires(geo::PlaneID{0, 0, 0}); // number of wires
1362  const int dy = detProp.ReadOutWindowSize(); // number of time samples.
1363 
1364  c.Init(dx, dy, fRhoResolutionFactor, fNumAngleCells);
1365 
1366  for (unsigned int i = 0; i < hits.size(); ++i) {
1367  c.AddPointReturnMax(hits[i]->WireID().Wire, (int)(hits[i]->PeakTime()));
1368  } // end loop over hits
1369 
1370  //gets the actual two-dimensional size of the accumulator
1371  int accDx = 0;
1372  int accDy = 0;
1373  c.GetAccumSize(accDy, accDx);
1374 
1375  //find the weightiest cell in the accumulator.
1376  int xMax = 0;
1377  int yMax = 0;
1378  c.GetMax(yMax, xMax);
1379 
1380  //find the center of mass of the 3x3 cell system (with maxCell at the center).
1381  float centerofmassx = 0.;
1382  float centerofmassy = 0.;
1383  float denom = 0.;
1384 
1385  if (xMax > 0 && yMax > 0 && xMax + 1 < accDx && yMax + 1 < accDy) {
1386  for (int i = -1; i < 2; ++i) {
1387  for (int j = -1; j < 2; ++j) {
1388  denom += c.GetCell(yMax + i, xMax + j);
1389  centerofmassx += j * c.GetCell(yMax + i, xMax + j);
1390  centerofmassy += i * c.GetCell(yMax + i, xMax + j);
1391  }
1392  }
1393  centerofmassx /= denom;
1394  centerofmassy /= denom;
1395  }
1396  else
1397  centerofmassx = centerofmassy = 0;
1398 
1399  float rho = 0.;
1400  float theta = 0.;
1401  c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
1402  MF_LOG_DEBUG("HoughBaseAlg") << "Transform(III) found maximum at (d=" << rho << " a=" << theta
1403  << ")"
1404  " from absolute maximum "
1405  << c.GetCell(yMax, xMax) << " at (d=" << yMax << ", a=" << xMax
1406  << ")";
1407  slope = -1. / tan(theta);
1408  intercept = rho / sin(theta);
1409 
1412 
1413  return hits.size();
1414 }
float fRhoResolutionFactor
Factor determining the resolution in rho.
Definition: HoughBaseAlg.h:534
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
int fNumAngleCells
that fall into the "correct" bin will be small and consistent with noise.
Definition: HoughBaseAlg.h:523
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
#define MF_LOG_DEBUG(id)
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.

Friends And Related Function Documentation

friend class HoughTransformClus
friend

Definition at line 514 of file HoughBaseAlg.h.

Member Data Documentation

float cluster::HoughBaseAlg::fMaxDistance
private

Max distance that a hit can be from a line to be considered part of that line.

Definition at line 528 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMaxLines
private

Max number of lines that can be found.

Definition at line 519 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMaxSlope
private

Max slope a line can have.

Definition at line 529 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMinHits
private

Min number of hits in the accumulator to consider (number of hits required to be considered a line).

Definition at line 520 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMissedHits
private

Number of wires that are allowed to be missed before a line is broken up into segments

Definition at line 539 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMissedHitsDistance
private

Distance between hits in a hough line before a hit is considered missed.

Definition at line 542 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMissedHitsToLineSize
private

Ratio of missed hits to line size for a line to be considered a fake.

Definition at line 544 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fNumAngleCells
private

that fall into the "correct" bin will be small and consistent with noise.

Number of angle cells in the accumulator (a measure of the angular resolution of the line finder). If this number is too large than the number of votes

Definition at line 523 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fPerCluster
private

Tells the original Hough algorithm to look at clusters individually, or all hits at once

Definition at line 536 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fRhoResolutionFactor
private

Factor determining the resolution in rho.

Definition at line 534 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fRhoZeroOutRange
private

Range in rho over which to zero out area around previously found lines in the accumulator.

Definition at line 531 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fSaveAccumulator
private

Save bitmap image of accumulator for debugging?

Definition at line 522 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fThetaZeroOutRange
private

Range in theta over which to zero out area around previously found lines in the accumulator.

Definition at line 533 of file HoughBaseAlg.h.


The documentation for this class was generated from the following files: