95 std::vector<recob::Slice> slcCol;
97 std::unique_ptr<art::Assns<recob::Slice, recob::Hit>>
100 std::unique_ptr<art::Assns<recob::Slice, recob::SpacePoint>>
107 std::vector< art::Ptr<recob::Hit> >
hits;
111 std::vector< art::Ptr<recob::SpacePoint> > sps;
115 if(!spFromHit.isValid()) {
116 std::cout<<
"spFromHit is invalid\n";
120 static bool first =
true;
122 bool success =
false;
123 for(
auto &
hit : hits) {
124 auto &sps = spFromHit.at(
hit.key());
125 if(sps.empty())
continue;
126 success = (sps[0].id() == spsHandle.id());
129 if(!success)
throw cet::exception(
"DBCluster3D")<<
"HitModuleLabel, SpacePointModuleLabel and SPHitAssnLabel are inconsistent\n";
134 if(!hitFromSp.isValid()) {
135 std::cout<<
"hitFromSp is invalid\n";
148 if (maxid>=0) nslc = maxid + 1;
151 std::vector<std::vector<art::Ptr<recob::Hit>>> slcHits(nslc);
153 std::map<geo::PlaneID, std::vector<std::pair<art::Ptr<recob::Hit>,
unsigned int>>> hitmap;
154 for (
auto &
hit : hits){
155 auto &sps = spFromHit.at(
hit.key());
165 for (
auto &
hit : hits){
167 for (
size_t i = 0; i<slcHits.size(); ++i){
168 if (std::find(slcHits[i].
begin(), slcHits[i].
end(),
hit) != slcHits[i].end()){
176 double x0 =
hit->WireID().Wire;
177 double y0 =
hit->PeakTime() * UnitsPerTick;
178 double mindis = DBL_MAX;
179 unsigned slcIndex = UINT_MAX;
181 double dx = hit2.first->WireID().Wire - x0;
182 double dy = hit2.first->PeakTime() * UnitsPerTick - y0;
183 double dis = dx*dx + dy*dy;
186 slcIndex = hit2.second;
189 if (slcIndex != UINT_MAX && mindis <
fMinHitDis){
190 slcHits[slcIndex].push_back(
hit);
196 std::vector<std::vector<art::Ptr<recob::SpacePoint>>> sps_in_slc(nslc);
204 for(
size_t isl = 0; isl < nslc; ++isl) {
205 double sum = sps_in_slc[isl].size();
207 double center[3] = {0.};
209 float charge = slcHits[isl].size();
210 for(
auto& spt : sps_in_slc[isl]) {
211 for(
unsigned short xyz = 0; xyz < 3; ++xyz) center[xyz] += spt->XYZ()[xyz];
213 for(
unsigned short xyz = 0; xyz < 3; ++xyz) center[xyz] /= sum;
215 double sumx = 0, sumy = 0., sumz = 0., sumx2 = 0, sumy2 = 0, sumz2 = 0;
216 double sumxy = 0, sumxz = 0;
217 for(
auto& spt : sps_in_slc[isl]) {
218 double xx = spt->XYZ()[0] - center[0];
219 double yy = spt->XYZ()[1] - center[1];
220 double zz = spt->XYZ()[2] - center[2];
230 double delta = sum * sumx2 - sumx * sumx;
231 if(delta <= 0)
continue;
233 double dydx = (sumxy * sum - sumx * sumy) / delta;
234 double dzdx = (sumxz * sum - sumx * sumz) / delta;
237 double norm = std::sqrt(1 + dydx * dydx + dzdx * dzdx);
238 direction[0] = 1 /
norm;
239 direction[1] = dydx /
norm;
240 direction[2] = dzdx /
norm;
242 unsigned int imax = 0, imin = 0;
243 float minAlong = 1E6, maxAlong = -1E6;
245 for(
unsigned int ipt = 0; ipt < sps_in_slc[isl].size(); ++ipt) {
246 auto& spt = sps_in_slc[isl][ipt];
247 for(
unsigned short xyz = 0; xyz < 3; ++xyz) tmpVec[xyz] = spt->XYZ()[xyz] - center[xyz];
249 for(
unsigned short xyz = 0; xyz < 3; ++xyz) dotp += tmpVec[xyz] * direction[xyz];
250 if(dotp < minAlong) { minAlong = dotp; imin = ipt; }
251 if(dotp > maxAlong) { maxAlong = dotp; imax = ipt; }
255 float aspectRatio = 0;
256 double arg = sum * sumy2 - sumy * sumy;
258 aspectRatio += std::abs(sum * sumxy - sumx * sumy) / sqrt(delta * arg);
261 arg = sum * sumz2 - sumz * sumz;
263 aspectRatio += std::abs(sum * sumxz - sumx * sumz) / sqrt(delta * arg);
266 if(cnt > 1) aspectRatio /= cnt;
268 Point_t ctr(center[0], center[1], center[2]);
269 Vector_t dir(direction[0], direction[1], direction[2]);
270 auto pos0 = sps_in_slc[isl][imin]->XYZ();
271 Point_t ep0(pos0[0], pos0[1], pos0[2]);
272 auto pos1 = sps_in_slc[isl][imax]->XYZ();
273 Point_t ep1(pos1[0], pos1[1], pos1[2]);
274 slcCol.emplace_back(
id, ctr,
dir, ep0, ep1, aspectRatio, charge);
279 std::unique_ptr<std::vector<recob::Slice> > scol(
new std::vector<recob::Slice>(std::move(slcCol)));
280 evt.put(std::move(scol));
281 evt.put(std::move(slc_hit_assn));
282 evt.put(std::move(slc_sps_assn));
geo::GeometryCore const * fGeom
The data type to uniquely identify a Plane.
const art::InputTag fSpacePointModuleLabel
const art::InputTag fHitModuleLabel
const art::InputTag fSPHitAssnLabel
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
recob::tracking::Point_t Point_t
void init(const std::vector< art::Ptr< recob::SpacePoint >> &sps, art::FindManyP< recob::Hit > &hitFromSp)
std::vector< evd::details::RawDigitInfo_t >::const_iterator begin(RawDigitCacheDataClass const &cache)
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
Detector simulation of raw signals on wires.
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
std::vector< point_t > points
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
recob::tracking::Vector_t Vector_t
cet::coded_exception< error, detail::translate > exception
1.8.11 
