diff --git a/TrkHitReco/src/CombineStrawHits_module.cc b/TrkHitReco/src/CombineStrawHits_module.cc index 22814b28da..88d456510f 100644 --- a/TrkHitReco/src/CombineStrawHits_module.cc +++ b/TrkHitReco/src/CombineStrawHits_module.cc @@ -21,6 +21,8 @@ #include "TMath.h" #include +#include +#include namespace mu2e { @@ -61,7 +63,10 @@ namespace mu2e { void produce( art::Event& e); private: - void combine(EventWindowMarker const& ewm, ComboHitCollection const& chcOrig, ComboHitCollection& chcol); + // origIndex maps a position in chcOrig to the corresponding index in chcol's parent + // collection; it is the identity unless the input was re-sorted first + void combine(EventWindowMarker const& ewm, ComboHitCollection const& chcOrig, ComboHitCollection& chcol, + std::vector const& origIndex); void combineHits(const ComboHitCollection& chcOrig, ComboHit& combohit); int _debug; @@ -143,20 +148,29 @@ namespace mu2e { // select hits based on flag panels[ch.strawId().uniquePanel()].push_back(ish); } + // chcolNew's parent is chcH (the unsorted input), so record where each sorted hit came + // from: without this the ComboHits would store positions in chcolsort while resolving + // them against chcOrig, silently corrupting per-hit provenance + std::vector origIndex; + origIndex.reserve(nsh); for (uint16_t ipanel=0;ipanel origIndex(chcOrig.size()); + std::iota(origIndex.begin(), origIndex.end(), 0); + combine(ewm, chcOrig, *chcolNew, origIndex); } event.put(std::move(chcolNew)); } - void CombineStrawHits::combine(EventWindowMarker const& ewm, ComboHitCollection const& chcOrig, ComboHitCollection& chcol) + void CombineStrawHits::combine(EventWindowMarker const& ewm, ComboHitCollection const& chcOrig, ComboHitCollection& chcol, + std::vector const& origIndex) { float minT = _minT; @@ -175,7 +189,11 @@ namespace mu2e { if ( _testflag && hit1.flag().hasAnyProperty(StrawHitFlag::dead)) continue; if ( _testflag && (!hit1.flag().hasAllProperties(_shsel) || hit1.flag().hasAnyProperty(_shmask))) continue; ComboHit combohit; - combohit.init(hit1,ich); + combohit.init(hit1,origIndex[ich]); + // TimeDivision is a quality tag: StrawHitRecoUtils sets it only when the longitudinal + // position was actually measured, and leaves it off when StrawResponse had to clamp the + // position to the straw end. It must therefore be propagated from the constituents. + bool tdiv = hit1.flag().hasAllProperties(StrawHitFlag::tdiv); int panel1 = hit1.strawId().uniquePanel(); for (size_t jch=ich+1;jch _maxwdchi) continue; - bool ok = combohit.addIndex(jch); + bool ok = combohit.addIndex(origIndex[jch]); if (!ok){ std::cout << "CombineStrawHits past limit" << std::endl; } else { isUsed[jch]= true; + tdiv &= hit2.flag().hasAllProperties(StrawHitFlag::tdiv); } } - // clear the flag bits; they are reset later - const static StrawHitFlag initialFlag("TimeDivision"); - combohit._flag = initialFlag; + // clear the flag bits; they are reset below + combohit._flag = StrawHitFlag(); + if(tdiv) combohit._flag.merge(StrawHitFlag::tdiv); int nch = combohit.nCombo(); if(nch < _minN || nch > _maxN){ - if(_filter)break; + if(_filter)continue; } else combohit._flag.merge(StrawHitFlag::nhitsel); // actually combine the hits if necessar, and make the cuts @@ -214,19 +233,19 @@ namespace mu2e { auto time = _useTOT ? combohit.correctedTime() : combohit.time(); if (time < minT || time > maxT ){ - if(_filter)break; + if(_filter)continue; } else combohit._flag.merge(StrawHitFlag::timesel); auto energy = combohit.energyDep(); if( energy > _maxE || energy < _minE ) { - if(_filter)break; + if(_filter)continue; } else combohit._flag.merge(StrawHitFlag::energysel); auto r2 = combohit.pos().Perp2(); if( r2 < _minR2 || r2 > _maxR2 ) { - if(_filter)break; + if(_filter)continue; } else combohit._flag.merge(StrawHitFlag::radsel); combohit._mask = _mask; @@ -249,7 +268,7 @@ namespace mu2e { { size_t index = combohit.index(ich); if (_debug > 3)std::cout << index << ", "; - if (index > chcOrig.size()) + if (index >= chcOrig.size()) throw cet::exception("RECO")<<"mu2e::CombineStrawHits: inconsistent index "< phirange(-M_PI,M_PI); + std::uniform_real_distribution phirange(-M_PI,M_PI); std::normal_distribution urand{0.0, ures}; std::normal_distribution vrand{0.0, vres}; diff --git a/TrkHitReco/src/StrawHitReco_module.cc b/TrkHitReco/src/StrawHitReco_module.cc index 1071e86aa4..21220f4d15 100644 --- a/TrkHitReco/src/StrawHitReco_module.cc +++ b/TrkHitReco/src/StrawHitReco_module.cc @@ -5,6 +5,7 @@ // Merged with flag and position creation B. Echenard, CalTech // // framework +#include "cetlib_except/exception.h" #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Handle.h" #include "Offline/GeometryService/inc/GeomHandle.hh" @@ -147,6 +148,9 @@ namespace mu2e { produces(); produces(); if (_writesh) produces(); + // flagCrossTalk reads the StrawHitCollection, which is only created when it is written out + if (_flagXT && !_writesh) + throw cet::exception("RECO")<<"mu2e::StrawHitReco: FlagCrossTalk requires WriteStrawHitCollection"< 0) std::cout << "In StrawHitReco constructor " << std::endl; } @@ -213,10 +217,11 @@ namespace mu2e { _shrUtils.createComboHit(ewm, isd, chCol, shCol, caloClusters, pbtOffset, digi.strawId(), digi.TDC(), digi.TOT(), pmp, trackerStatus, srep, tt); - //flag straw and electronic cross-talk - if(_flagXT){ - _shrUtils.flagCrossTalk(shCol, chCol); - } + } + //flag straw and electronic cross-talk: this scans the whole collection, so it must run once + //per event, not once per digi + if(_flagXT){ + _shrUtils.flagCrossTalk(shCol, chCol); } if(_writesh)event.put(std::move(shCol)); intInfo->setNTrackerHits(chCol->size());