1015 std::vector<int> skip;
1018 lariov::ChannelStatusProvider
const* channelStatus =
1019 lar::providerFrom<lariov::ChannelStatusService>();
1021 CLHEP::RandFlat flat(engine);
1023 std::vector<art::Ptr<recob::Hit>>
hit;
1029 hit.push_back((*ii));
1031 if (hit.size() == 0) {
1041 if (hit.size() == 0) {
1046 auto const num_planes = geom->
Nplanes(tpcid);
1047 std::vector<double> wire_pitch(num_planes, 0.);
1049 auto const p = plane.ID().Plane;
1050 wire_pitch[p] = plane.WirePitch();
1054 std::vector<double> xyScale(num_planes, 0.);
1057 double driftVelFactor = 0.001 * detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
1059 for (
size_t p = 0; p < xyScale.size(); ++p)
1060 xyScale[p] = driftVelFactor *
sampling_rate(clockData) / wire_pitch[p];
1064 const int dy = detProp.ReadOutWindowSize();
1066 skip.resize(hit.size());
1067 std::vector<int> listofxmax;
1068 std::vector<int> listofymax;
1069 std::vector<int> hitTemp;
1070 std::vector<int> sequenceHolder;
1071 std::vector<int> currentHits;
1072 std::vector<int> lastHits;
1074 float indcolscaling = 0.;
1076 float centerofmassx = 0;
1077 float centerofmassy = 0;
1079 float intercept = 0.;
1086 int accDx(0), accDy(0);
1095 unsigned int count = hit.size();
1096 std::vector<unsigned int> accumPoints;
1097 accumPoints.resize(hit.size());
1099 unsigned int nLinesFound = 0;
1101 for (; count > 0; count--) {
1104 unsigned int randInd = (
unsigned int)(flat.fire() * hit.size());
1106 MF_LOG_DEBUG(
"HoughBaseAlg") <<
"randInd=" << randInd <<
" and size is " << hit.size();
1109 if (skip.at(randInd) == 1)
continue;
1112 if (accumPoints.at(randInd))
continue;
1113 accumPoints.at(randInd) = 1;
1116 for (
unsigned int i = 0; i < listofxmax.size(); ++i) {
1118 if (yClearStart < 0) yClearStart = 0;
1121 if (yClearEnd >= accDy) yClearEnd = accDy - 1;
1124 if (xClearStart < 0) xClearStart = 0;
1127 if (xClearEnd >= accDx) xClearEnd = accDx - 1;
1129 for (
y = yClearStart;
y <= yClearEnd; ++
y) {
1130 for (x = xClearStart; x <= xClearEnd; ++
x) {
1138 unsigned int wireMax = hit.at(randInd)->WireID().Wire;
1141 std::array<int, 3> max = c.AddPointReturnMax(wireMax, (
int)(hit.at(randInd)->PeakTime()));
1142 maxCell = max.at(0);
1151 denom = centerofmassx = centerofmassy = 0;
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);
1158 centerofmassx += j * cell;
1159 centerofmassy += i * cell;
1162 centerofmassx /= denom;
1163 centerofmassy /= denom;
1166 centerofmassx = centerofmassy = 0;
1169 listofxmax.push_back(xMax);
1170 listofymax.push_back(yMax);
1173 c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
1174 MF_LOG_DEBUG(
"HoughBaseAlg") <<
"Transform(I) found maximum at (d=" << rho <<
" a=" << theta
1176 " from absolute maximum " 1177 << c.GetCell(yMax, xMax) <<
" at (d=" << yMax <<
", a=" << xMax
1179 slope = -1. / tan(theta);
1180 intercept = (rho / sin(theta));
1191 if (!std::isinf(slope) && !std::isnan(slope)) {
1192 sequenceHolder.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) -
1197 (std::sqrt(pow(xyScale[hit.at(i)->WireID().Plane] * slope, 2) + 1)));
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());
1206 if (hitTemp.size() < 2)
continue;
1207 currentHits.clear();
1210 currentHits.push_back(0);
1211 for (
size_t i = 0; i + 1 < sequenceHolder.size(); ++i) {
1213 while ((channelStatus->IsBad(sequenceHolder.at(i) + j)) ==
true)
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();
1222 currentHits.clear();
1225 if (currentHits.size() > lastHits.size()) lastHits = currentHits;
1233 double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
1236 int lastHitsChannel = 0;
1237 int nHitsPerChannel = 1;
1239 MF_LOG_DEBUG(
"HoughBaseAlg") <<
"filling the pCorner arrays around here..." 1240 <<
"\n but first, lastHits size is " << lastHits.size()
1241 <<
" and lastHitsChannel=" << lastHitsChannel;
1245 unsigned int lastChannel = hit.at(hitTemp.at(lastHits.at(0)))->Channel();
1247 for (
size_t i = 0; i < lastHits.size() - 1; ++i) {
1248 bool newChannel =
false;
1250 if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() != lastChannel) {
1253 if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() == lastChannel) nHitsPerChannel++;
1259 if (slope > 0 || (!newChannel && nHitsPerChannel <= 1)) {
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)) >
1270 else if (slope < 0 && newChannel && nHitsPerChannel > 1) {
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)) >
1280 lastChannel = hit.at(hitTemp.at(lastHits.at(i)))->Channel();
1281 lastHitsChannel = i + 1;
1282 nHitsPerChannel = 0;
1288 clusterHits.
clear();
1289 if (lastHits.size() < 5)
continue;
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;
1303 clusHitsOut.push_back(clusterHits);
1304 slopevec.push_back(slope);
1305 totalQvec.emplace_back(integralQ.
Sum(),
1315 if (nLinesFound > (
unsigned int)
fMaxLines)
break;
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;
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) {
1336 if (maxCell > 0) pix = (int)((1500 * c.GetCell(
y, x)) / maxCell);
1337 outPix[PicIndex++] = pix;
1348 return clusHitsOut.size();
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
float fRhoResolutionFactor
Factor determining the resolution in rho.
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.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
float fMissedHitsDistance
Distance between hits in a hough line before a hit is considered missed.
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.
int fSaveAccumulator
Save bitmap image of accumulator for debugging?
void push_back(Ptr< U > const &p)
Signal from induction planes.
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...
float fMaxDistance
Max distance that a hit can be from a line to be considered part of that line.
float fMaxSlope
Max slope a line can have.
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Description of geometry of one entire detector.
int fRhoZeroOutRange
Range in rho over which to zero out area around previously found lines in the accumulator.
Detector simulation of raw signals on wires.
constexpr TPCID const & asTPCID() const
Conversion to TPCID (for convenience of notation).
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
int fThetaZeroOutRange
Range in theta over which to zero out area around previously found lines in the accumulator.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Collects statistics on a single quantity (weighted)
float fMissedHitsToLineSize
Ratio of missed hits to line size for a line to be considered a fake.
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)