33 #include "range/v3/view.hpp" 36 using ranges::views::filter;
37 using ranges::views::transform;
43 ,
fCrawl{p.get<
bool>(
"Crawl")}
73 std::vector<art::Ptr<recob::Hit>> ChHits;
77 unsigned int skipNoise(0);
79 for (
size_t h = 0; h < ChHits.size(); h++) {
92 unsigned int apa(0), cryo(0);
100 std::pair<double, double> ChanTime(hit->
Channel() * 1., hit->
PeakTime() * 1.);
109 <<
"\nSkipped " << skipNoise <<
" induction noise hits using the BackTrackerService.\n" 110 <<
"This is only to temporarily deal with the excessive amount of noise due to the bad " 113 mf::LogVerbatim(
"RunDisambig") <<
"\n~~~~~~~~~~~ Running Disambiguation ~~~~~~~~~~~\n";
115 std::map<unsigned int, std::vector<art::Ptr<recob::Hit>>>
::iterator APA_it;
117 unsigned int apa = APA_it->first;
132 <<
" Trivial Disambig --> " <<
fnDUSoFar[apa] <<
" / " <<
fnUSoFar[apa] <<
" U, " 141 <<
fnDVSoFar[apa] <<
" / " << fnVSoFar[apa] <<
" V";
150 <<
fnDVSoFar[apa] <<
" / " << fnVSoFar[apa] <<
" V";
154 unsigned int nDisambig(1);
155 while (nDisambig > 0) {
162 <<
fnDVSoFar[apa] <<
" / " << fnVSoFar[apa] <<
" V";
166 for (
size_t i = 0; i <
fAPAToDHits[apa].size(); i++)
178 std::pair<double, double> ChanTime(hit->
Channel() * 1., hit->
PeakTime() * 1.);
182 mf::LogWarning(
"InvalidWireID") <<
"wid is invalid, hit not being made\n";
228 return (AsT <= BsT && BsT <= AeT) || (AsT <= BeT && BeT <= AeT) || (BsT <= AsT && AsT <= BeT) ||
229 (BsT <= AeT && AeT <= BeT);
240 auto const&
hit = *hitPtr;
242 unsigned int peakT =
hit.PeakTime();
245 std::vector<bool> IsReasonableWid(hitwids.size(),
false);
246 unsigned short nPossibleWids(0);
247 for (
size_t w = 0;
w < hitwids.size();
w++) {
251 unsigned int side(wid.
TPC % 2), cryo(wid.
Cryostat);
252 double zminPos(start.Z()), zmaxPos(
end.Z());
256 2 * apa + side - cryo *
geom->
NTPC();
260 auto Min = tpcCenter;
262 auto Max = tpcCenter;
271 if (chan <= ZminChan || ZmaxChan <= chan)
continue;
274 IsReasonableWid[
w] =
true;
282 if (nPossibleWids == 0) {
283 std::vector<double> xyz;
292 <<
"U/V hit inconsistent with Z info; peak time is " << peakT <<
" in APA " << apa
293 <<
" on channel " <<
hit.Channel();
295 else if (nPossibleWids == 1) {
296 for (
size_t d = 0;
d < hitwids.size();
d++)
299 else if (nPossibleWids == 2) {
319 throw cet::exception(
"MakeCloseHits") <<
"Function not meant for non-wrapped channels.\n";
324 int tempchan = Dchan + ext;
325 if (tempchan < (
int)firstChan) tempchan += ChanPerView;
326 if (tempchan > (
int)(firstChan + ChanPerView - 1)) tempchan -= ChanPerView;
333 unsigned int apa(0), cryo(0);
335 unsigned int MakeCount(0);
342 if (!(Dmin <= st && st <= Dmax) && !(Dmin <= et && et <= Dmax))
continue;
346 for (
size_t w = 0;
w < wids.size();
w++) {
347 if (wids[
w].TPC != Dwid.
TPC)
continue;
348 if ((
int)(wids[
w].Wire) - (
int)(Dwid.
Wire) != ext)
continue;
352 std::pair<double, double> ChanTime(closeHit->
Channel() * 1., closeHit->
PeakTime() * 1.);
373 unsigned int nExtended(1);
374 while (nExtended > 0) {
378 for (
size_t h = 0; h < hits.size(); h++) {
379 std::pair<double, double> ChanTime(hits[h]->Channel() * 1., hits[h]->PeakTime() * 1.);
381 double stD = hits[h]->PeakTimePlusRMS(-1.);
382 double etD = hits[h]->PeakTimePlusRMS(+1.);
383 double hitWindow = etD - stD;
388 unsigned int extensions = 0;
389 for (
unsigned int ext = 1; ext <
fNChanJumps + 1; ext++) {
392 double timeExt = hitWindow * ext;
393 N +=
MakeCloseHits((
int)(-ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
394 N +=
MakeCloseHits((
int)(ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
397 nExtended += extensions;
414 double pi = 3.14159265;
417 for (
size_t h = 0; h <
fAPAToHits[apa].size(); h++) {
420 unsigned int plane = 0;
421 if (view ==
geo::kV) { plane = 1; }
424 std::vector<double> ChanTimeCenter(2, 0.);
427 ChanTimeCenter[0] = relchan * wire_pitch;
433 std::vector<std::vector<double>> CloseHitsChanTime;
434 std::vector<double> FurthestCloseChanTime(2, 0.);
435 std::vector<double> ClosestChanTime(2, 0.);
439 for (
size_t c = 0; c <
fAPAToHits[apa].size(); c++) {
441 if (view != closehit->
View())
continue;
443 unsigned int plane = 0;
444 if (view ==
geo::kV) { plane = 1; }
447 std::vector<double> ChanTimeClose(2, 0.);
448 unsigned int relchanclose =
450 ChanTimeClose[0] = relchanclose * wire_pitch;
455 if (ChanTimeClose == ChanTimeCenter)
continue;
457 double ChanDist = ChanTimeClose[0] - ChanTimeCenter[0];
458 if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
460 double distance = std::hypot(ChanDist, ChanTimeClose[1] - ChanTimeCenter[1]);
462 if (distance <=
fCloseHitsRadius) CloseHitsChanTime.push_back(ChanTimeClose);
464 if (distance < minDist) {
465 ClosestChanTime = ChanTimeClose;
471 if (CloseHitsChanTime.size() < 5)
continue;
473 double minRad(2 * pi + 1.), maxRad(0.);
474 bool CloseToNegPi(
false), CloseToPosPi(
false);
475 for (
size_t i = 0; i < CloseHitsChanTime.size(); i++) {
476 std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
477 double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
478 if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
479 double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
480 if (hitrad > maxRad) maxRad = hitrad;
481 if (hitrad < minRad) minRad = hitrad;
482 if (hitrad + fMaxEndPRadRange > pi)
484 else if (hitrad - fMaxEndPRadRange < -pi)
489 if (CloseToPosPi && CloseToNegPi) {
490 for (
size_t i = 0; i < CloseHitsChanTime.size(); i++) {
491 std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
492 double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
493 if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
494 double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
495 if (hitrad > 0) hitrad = pi - hitrad;
496 if (hitrad < 0) hitrad = -pi - hitrad;
497 if (hitrad > maxRad) maxRad = hitrad;
498 if (hitrad < minRad) minRad = hitrad;
502 if (maxRad - minRad < fMaxEndPRadRange)
fAPAToEndPHits[apa].push_back(centhit);
508 <<
" endpoint hits in apa " << apa << std::endl;
512 <<
" at time " << epHit->
PeakTime() << std::endl;
525 mf::LogVerbatim(
"UseEndPts") <<
" APA " << apa <<
" has no endpoints.";
528 std::vector<art::Ptr<recob::Hit>>
const& endPts =
fAPAToEndPHits[apa];
530 std::vector<std::vector<art::Ptr<recob::Hit>>> EndPMatch;
531 unsigned short nZendPts(0);
535 for (
auto const& ZHitPtr : endPts |
filter(on_z_plane)) {
536 auto const& ZHit = *ZHitPtr;
539 unsigned short Umatch(0), Vmatch(0);
543 for (
auto const& hitPtr : endPts |
filter(not_on_z_plane)) {
544 auto const&
hit = *hitPtr;
557 unsigned int tpc(ZHit.WireID().TPC), cryo(ZHit.WireID().Cryostat);
560 if (Umatch == 1 && Vmatch == 1) {
562 std::vector<double> yzEndPt =
565 geo::Point_t const intersect{tpcCenter.X(), yzEndPt[0], yzEndPt[1]};
571 else if (Umatch == 1 && Vmatch != 1) {
573 std::vector<geo::WireIDIntersection> widIntersects;
575 if (widIntersects.size() == 0)
577 else if (widIntersects.size() == 1) {
578 geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
584 for (
size_t i = 0; i < widIntersects.size(); i++) {
589 else if (Umatch == 1 && Vmatch != 1) {
591 std::vector<geo::WireIDIntersection> widIntersects;
593 if (widIntersects.size() == 0)
595 else if (widIntersects.size() == 1) {
596 geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
604 if (nZendPts == 0 && endPts.size() == 2 &&
HitsOverlapInTime(detProp, *endPts[0], *endPts[1])) {
605 std::vector<geo::WireIDIntersection> widIntersects;
607 if (widIntersects.size() == 1) {
608 unsigned int cryo = endPts[0]->WireID().Cryostat;
609 unsigned int tpc = widIntersects[0].TPC;
611 geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
612 unsigned int plane0(0), plane1(0);
613 if (endPts[0]->View() ==
geo::kV) plane0 = 1;
614 if (endPts[1]->View() ==
geo::kV) plane1 = 1;
631 unsigned int nU(0), nV(0);
640 unsigned int nDU(0), nDV(0);
641 for (
size_t h = 0; h <
fAPAToDHits[apa].size(); h++) {
662 unsigned int nDisambiguations(0);
666 auto const& ambighit = *ambighitPtr;
668 std::pair<double, double> ambigChanTime(ambigchan * 1., ambighit.PeakTime());
672 std::vector<unsigned int> widDcounts(ambigwids.size(), 0);
673 std::vector<unsigned int> widAcounts(ambigwids.size(), 0);
683 std::pair<double, double> ChanTime(chan * 1.,
hit.PeakTime());
685 for (
size_t a = 0; a < ambigwids.size(); a++)
693 for (
size_t a = 0; a < ambigwids.size(); a++)
694 for (
size_t w = 0;
w < wids.size();
w++)
695 if (ambigwids[a].TPC == wids[
w].TPC &&
703 unsigned int Dcount(0), Acount(0);
704 for (
size_t d = 0;
d < widDcounts.size();
d++)
705 Dcount += widDcounts[
d];
706 for (
size_t a = 0; a < widAcounts.size(); a++)
707 Acount += widAcounts[a];
708 for (
size_t d = 0;
d < widDcounts.size();
d++) {
709 if (Dcount == widDcounts[
d] && Dcount > 0 && Acount == 0) {
716 return nDisambiguations;
std::map< unsigned int, std::map< std::pair< double, double >, bool > > fHasBeenDisambiged
Convenient way to keep track of disambiguation so far.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::map< raw::ChannelID_t, std::vector< art::Ptr< recob::Hit > > > fChannelToHits
constexpr to_element_t to_element
Encapsulate the construction of a single cyostat .
bool HitsOverlapInTime(detinfo::DetectorPropertiesData const &detProp, recob::Hit const &hitA, recob::Hit const &hitB)
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToZHits
std::map< unsigned int, unsigned int > fnDVSoFar
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
apa::APAGeometryAlg fAPAGeo
Declaration of signal hit object.
std::map< unsigned int, unsigned int > fnVSoFar
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
uint32_t FirstChannelInView(geo::View_t geoview, unsigned int apa, unsigned int cryo) const
void Crawl(unsigned int apa)
Extend what we disambiguation we do have in apa.
double TimeOffsetZ() const
CryostatID_t Cryostat
Index of cryostat.
Planes which measure Z direction.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
geo::View_t View() const
View for the plane of the hit.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
void RunDisambig(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::Handle< std::vector< recob::Hit >> GausHits)
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToEndPHits
\ todo: Channel/APA to hits can be done in a unified way
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
bool APAChannelsIntersect(uint32_t chan1, uint32_t chan2, std::vector< geo::WireIDIntersection > &IntersectVector) const
If the channels intersect, get all intersections.
art::ServiceHandle< geo::Geometry const > geom
void TrivialDisambig(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Make the easiest and safest disambiguations in apa.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Point_t GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
unsigned int ChannelToAPA(uint32_t chan) const
Get number of the APA containing the given channel.
void AssessDisambigSoFar(unsigned int apa)
See how much disambiguation has been done in this apa so far.
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToUVHits
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToHits
DisambigAlg(fhicl::ParameterSet const &pset)
unsigned int ChannelsInView(geo::View_t geoview) const
virtual raw::ChannelID_t PlaneWireToChannel(WireID const &wireID) const =0
Returns the channel ID a wire is connected to.
std::map< unsigned int, double > fVeffSoFar
void MakeDisambigHit(art::Ptr< recob::Hit > const &hit, geo::WireID, unsigned int apa)
Makes a disambiguated hit while keeping track of what has already been disambiguated.
The data type to uniquely identify a TPC.
unsigned int fNChanJumps
Number of channels the crawl can jump over.
void UseEndPts(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Detector simulation of raw signals on wires.
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
float PeakTime() const
Time of the signal peak, in tick units.
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
Contains all timing reference information for the detector.
View_t View(raw::ChannelID_t const channel) const
Returns the view (wire orientation) on the specified TPC channel.
geo::WireID NearestWireIDOnChan(geo::Point_t const &WorldLoc, uint32_t chan, geo::PlaneID const &planeID) const
bool fCrawl
\ todo: Write function that compares hits more detailedly
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
double TimeOffsetU() const
std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
The final list of hits to pass back to be made.
unsigned int CompareViews(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Compare U and V to see if one says something about the other.
std::map< unsigned int, double > fUeffSoFar
geo::WireReadoutGeom const * fWireReadoutGeom
2D representation of charge deposited in the TDC/wire plane
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
std::vector< double > ThreeChanPos(uint32_t u, uint32_t v, uint32_t z) const
Find the center of the 3 intersections, choose best if multiple.
std::map< unsigned int, unsigned int > fnUSoFar
unsigned int FindChanTimeEndPts(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Basic endpoint-hit finder per apa.
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
std::map< unsigned int, std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > > fAPAToDHits
Hold the disambiguations per APA.
unsigned int MakeCloseHits(int ext, geo::WireID wid, double Dmin, double Dmax)
Having disambiguated a time range on a wireID, extend to neighboring channels.
double TimeOffsetV() const
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane .
std::map< std::pair< double, double >, geo::WireID > fChanTimeToWid
If a hit is disambiguated, map its chan and peak time to the chosen wireID.
std::map< unsigned int, unsigned int > fnDUSoFar