89 auto scol = std::make_unique<std::vector<recob::Slice>>();
90 auto slc_hit_assn = std::make_unique<art::Assns<recob::Slice, recob::Hit>>();
91 auto slc_sps_assn = std::make_unique<art::Assns<recob::Slice, recob::SpacePoint>>();
97 std::vector<art::Ptr<recob::Hit>>
hits;
101 std::vector<art::Ptr<recob::SpacePoint>> sps;
105 if (!spFromHit.isValid()) {
106 std::cout <<
"spFromHit is invalid\n";
110 static bool first =
true;
112 bool success =
false;
113 bool foundsps =
false;
114 for (
auto&
hit : hits) {
115 auto& sps = spFromHit.at(
hit.key());
116 if (sps.empty())
continue;
117 success = (sps[0].id() == spsHandle.id());
121 if ((!success) && foundsps)
123 <<
"HitModuleLabel, SpacePointModuleLabel and SPHitAssnLabel are inconsistent\n";
128 if (!hitFromSp.isValid()) {
129 std::cout <<
"hitFromSp is invalid\n";
141 if (maxid >= 0) nslc = maxid + 1;
144 std::vector<std::vector<art::Ptr<recob::Hit>>> slcHits(nslc);
146 std::map<geo::PlaneID, std::vector<std::pair<art::Ptr<recob::Hit>,
unsigned int>>> hitmap;
147 for (
auto&
hit : hits) {
148 auto& sps = spFromHit.at(
hit.key());
159 for (
auto&
hit : hits) {
161 for (
size_t i = 0; i < slcHits.size(); ++i) {
162 if (std::find(slcHits[i].
begin(), slcHits[i].
end(),
hit) != slcHits[i].end()) {
170 double x0 =
hit->WireID().Wire;
171 double y0 =
hit->PeakTime() * UnitsPerTick;
172 double mindis = DBL_MAX;
173 unsigned slcIndex = UINT_MAX;
175 double dx = hit2.first->WireID().Wire - x0;
176 double dy = hit2.first->PeakTime() * UnitsPerTick - y0;
177 double dis = dx * dx + dy * dy;
180 slcIndex = hit2.second;
183 if (slcIndex != UINT_MAX && mindis <
fMinHitDis) { slcHits[slcIndex].push_back(
hit); }
188 std::vector<std::vector<art::Ptr<recob::SpacePoint>>> sps_in_slc(nslc);
196 for (
size_t isl = 0; isl < nslc; ++isl) {
197 double sum = sps_in_slc[isl].size();
199 double center[3] = {0.};
201 float charge = slcHits[isl].size();
202 for (
auto& spt : sps_in_slc[isl]) {
203 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
204 center[xyz] += spt->XYZ()[xyz];
206 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
209 double sumx = 0, sumy = 0., sumz = 0., sumx2 = 0, sumy2 = 0, sumz2 = 0;
210 double sumxy = 0, sumxz = 0;
211 for (
auto& spt : sps_in_slc[isl]) {
212 double xx = spt->XYZ()[0] - center[0];
213 double yy = spt->XYZ()[1] - center[1];
214 double zz = spt->XYZ()[2] - center[2];
224 double delta = sum * sumx2 - sumx * sumx;
225 if (delta <= 0)
continue;
227 double dydx = (sumxy * sum - sumx * sumy) / delta;
228 double dzdx = (sumxz * sum - sumx * sumz) / delta;
231 double norm = std::sqrt(1 + dydx * dydx + dzdx * dzdx);
232 direction[0] = 1 /
norm;
233 direction[1] = dydx /
norm;
234 direction[2] = dzdx /
norm;
236 unsigned int imax = 0, imin = 0;
237 float minAlong = 1E6, maxAlong = -1E6;
239 for (
unsigned int ipt = 0; ipt < sps_in_slc[isl].size(); ++ipt) {
240 auto& spt = sps_in_slc[isl][ipt];
241 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
242 tmpVec[xyz] = spt->XYZ()[xyz] - center[xyz];
244 for (
unsigned short xyz = 0; xyz < 3; ++xyz)
245 dotp += tmpVec[xyz] * direction[xyz];
246 if (dotp < minAlong) {
250 if (dotp > maxAlong) {
257 float aspectRatio = 0;
258 double arg = sum * sumy2 - sumy * sumy;
260 aspectRatio +=
std::abs(sum * sumxy - sumx * sumy) / sqrt(delta * arg);
263 arg = sum * sumz2 - sumz * sumz;
265 aspectRatio +=
std::abs(sum * sumxz - sumx * sumz) / sqrt(delta * arg);
268 if (cnt > 1) aspectRatio /= cnt;
270 Point_t ctr(center[0], center[1], center[2]);
271 Vector_t dir(direction[0], direction[1], direction[2]);
272 auto pos0 = sps_in_slc[isl][imin]->XYZ();
273 Point_t ep0(pos0[0], pos0[1], pos0[2]);
274 auto pos1 = sps_in_slc[isl][imax]->XYZ();
275 Point_t ep1(pos1[0], pos1[1], pos1[2]);
276 scol->emplace_back(
id, ctr,
dir, ep0, ep1, aspectRatio, charge);
281 evt.put(std::move(scol));
282 evt.put(std::move(slc_hit_assn));
283 evt.put(std::move(slc_sps_assn));
geo::GeometryCore const * fGeom
The data type to uniquely identify a Plane.
constexpr auto abs(T v)
Returns the absolute value of the argument.
const art::InputTag fSpacePointModuleLabel
const art::InputTag fHitModuleLabel
const art::InputTag fSPHitAssnLabel
recob::tracking::Point_t Point_t
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
void init(const std::vector< art::Ptr< recob::SpacePoint >> &sps, art::FindManyP< recob::Hit > &hitFromSp)
Detector simulation of raw signals on wires.
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
std::vector< point_t > points
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
recob::tracking::Vector_t Vector_t
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception
1.8.11 
