diff --git a/CosmicReco/src/LineFinder_module.cc b/CosmicReco/src/LineFinder_module.cc index c596ed079..83b06f9c3 100644 --- a/CosmicReco/src/LineFinder_module.cc +++ b/CosmicReco/src/LineFinder_module.cc @@ -77,7 +77,7 @@ namespace mu2e{ ProditionsHandle _alignedTracker_h; - int findLine(const ComboHitCollection& shC, art::Event const& event, CosmicTrackSeed &tseed); + int findLine(const ComboHitCollection& shC, std::vector const& shiv, art::Event const& event, CosmicTrackSeed& tseed); }; @@ -112,24 +112,18 @@ void LineFinder::produce(art::Event& event ) { for (size_t index=0;index< tccol.size();++index) { const auto& tclust = tccol[index]; + std::vector shiv; + // get collection at straw level + auto chcp = chcol.fillStrawHitIndices(tclust.hits(),shiv,StrawIdMask::uniquestraw); + auto chc = *chcp; - // convert (presumably) panel hits in time clusterto straw hits CosmicTrackSeed tseed; - ComboHitCollection tchits; - tchits.reserve(tclust.hits().size()*2); // guess 2 straw hits/panel hit - - int nhits = 0; - ComboHitCollection::CHCIter chids; - chcol.fillComboHits( tclust.hits(), chids); - for (auto const& it : chids){ - tchits.push_back(*it); - nhits += it->nStrawHits(); - } + tseed._straw_chits.setAsSubset(chH,StrawIdMask::uniquestraw); tseed._timeCluster = art::Ptr(tcH,index); tseed._track.converged = true; - int seedSize = findLine(tchits, event, tseed); + int seedSize = findLine(chc, shiv, event, tseed); if (_diag > 1) std::cout << "LineFinder: seedSize = " << seedSize << std::endl; @@ -148,7 +142,7 @@ void LineFinder::produce(art::Event& event ) { event.put(std::move(seed_col)); } -int LineFinder::findLine(const ComboHitCollection& shC, art::Event const& event, CosmicTrackSeed& tseed){ +int LineFinder::findLine(const ComboHitCollection& shC, std::vector const& shiv, art::Event const& event, CosmicTrackSeed& tseed){ //mu2e::GeomHandle th; auto tracker = _alignedTracker_h.getPtr(event.id());//th.get(); @@ -160,18 +154,20 @@ int LineFinder::findLine(const ComboHitCollection& shC, art::Event const& event, bool found_all = false; // lets get the best pairwise vector - for (size_t i=0;igetStraw(shC[i].strawId()); - for (size_t j=i+1;jgetStraw(shC[iloc].strawId()); + for (size_t j=i+1;jgetStraw(shC[j].strawId()); + Straw const& strawj = tracker->getStraw(shC[jloc].strawId()); for (int is=-1*_Nsteps;is<_Nsteps+1;is++){ - CLHEP::Hep3Vector ipos = shC[i].posCLHEP() + strawi.getDirection()*shC[i].wireRes()*_stepSize*is; + CLHEP::Hep3Vector ipos = shC[iloc].posCLHEP() + strawi.getDirection()*shC[iloc].wireRes()*_stepSize*is; for (int js=-1*_Nsteps;js<_Nsteps+1;js++){ - CLHEP::Hep3Vector jpos = shC[j].posCLHEP() + strawj.getDirection()*shC[j].wireRes()*_stepSize*js; + CLHEP::Hep3Vector jpos = shC[jloc].posCLHEP() + strawj.getDirection()*shC[jloc].wireRes()*_stepSize*js; CLHEP::Hep3Vector newdir = (jpos-ipos).unit(); CLHEP::Hep3Vector icross = (jpos-ipos).cross(strawi.getDirection()).unit(); @@ -184,17 +180,18 @@ int LineFinder::findLine(const ComboHitCollection& shC, art::Event const& event, // now loop over all hits and see how many are in this track int count = 0; double ll = 0; - for (size_t k=0;kgetStraw(shC[k].strawId()); + for (size_t k=0;kgetStraw(shC[kloc].strawId()); TwoLinePCA pca( strawk.getMidPoint(), strawk.getDirection(), ipos, newdir); double dist = (pca.point1()-strawk.getMidPoint()).mag(); if (pca.dca() < _maxDOCA && dist < strawk.halfLength()){ count += 1; - ll += pow(dist-shC[k].wireDist(),2)/shC[k].wireVar(); + ll += pow(dist-shC[kloc].wireDist(),2)/shC[kloc].wireVar(); } } - if (count == (int) shC.size()) + if (count == (int) shiv.size()) found_all = true; if (count > bestcount || (count == bestcount && ll < bestll)){ bestcount = count; @@ -219,17 +216,20 @@ int LineFinder::findLine(const ComboHitCollection& shC, art::Event const& event, double avg_t0 = 0; int good_hits = 0; - for (size_t k=0;kgetStraw(shC[k].strawId()); + for (size_t k=0;kgetStraw(shC[kloc].strawId()); TwoLinePCA pca( strawk.getMidPoint(), strawk.getDirection(), seedInt, seedDir); double dist = (pca.point1()-strawk.getMidPoint()).dot(strawk.getDirection()); if (pca.dca() < _maxDOCA && fabs(dist) < strawk.halfLength()){ double traj_time = ((pca.point2() - seedInt).dot(seedDir))/299.9; - double hit_t0 = shC[k].time() - shC[k].driftTime() - shC[k].propTime() - traj_time; + double hit_t0 = shC[kloc].time() - shC[kloc].driftTime() - shC[kloc].propTime() - traj_time; avg_t0 += hit_t0; good_hits++; - tseed._straw_chits.push_back(shC[k]); + ComboHit combohit; + combohit.init(shC[kloc],kloc); + tseed._straw_chits.push_back(std::move(combohit)); } } avg_t0 /= good_hits; @@ -258,22 +258,6 @@ int LineFinder::findLine(const ComboHitCollection& shC, art::Event const& event, tseed._track.SetFitEquation(XYZTrack); tseed._track.SetMinuitEquation(XYZTrack); - // For compatibility FIXME - for(size_t ich= 0; ich shitids; - tseed._straw_chits.fillStrawHitIndices(ich, shitids); - - /* - for(auto const& ids : shitids){ - size_t istraw = (ids); - TrkStrawHitSeed tshs; - tshs._index = istraw; - tshs._t0 = tseed._t0; - tseed._trkstrawhits.push_back(tshs); - } - */ - } - return good_hits; } diff --git a/RecoDataProducts/inc/ComboHit.hh b/RecoDataProducts/inc/ComboHit.hh index 7c9597d11..e4d009a87 100644 --- a/RecoDataProducts/inc/ComboHit.hh +++ b/RecoDataProducts/inc/ComboHit.hh @@ -150,15 +150,16 @@ namespace mu2e { // This function is called recursively, so the the vector must be empty on the top-most call #ifndef __ROOTCLING__ // find the parent at a given level - CHCPTR parent(StrawIdMask::Level level) const; - void fillStrawDigiIndices( size_t chindex, SHIV& shids) const; + // if parent and grandparent are the same level, will select the grandparent unless stopatfirst set + CHCPTR parent(StrawIdMask::Level level, bool stopatfirst=false) const; + void fillStrawDigiIndices( size_t chindex, SHIV& shids, bool stopatfirst=false) const; // Fill indices to the specified level. Return value is the collection to whic // the indices apply. first, given all my hits - ComboHitCollection const* fillStrawHitIndices( SHIV& shiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw) const; + ComboHitCollection const* fillStrawHitIndices( SHIV& shiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw, bool stopatfirst=false) const; // given a specific hit (index) in myself - ComboHitCollection const* fillStrawHitIndices( size_t chindex, SHIV& shiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw) const; + ComboHitCollection const* fillStrawHitIndices( size_t chindex, SHIV& shiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw, bool stopatfirst=false) const; // given a vector of indices - ComboHitCollection const* fillStrawHitIndices(SHIV const& inshiv, SHIV& outshiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw) const; + ComboHitCollection const* fillStrawHitIndices(SHIV const& inshiv, SHIV& outshiv, StrawIdMask::Level clevel=StrawIdMask::uniquestraw, bool stopatfirst=false) const; // the following are deprecated in favor of the more-efficient and self-checking functions above // translate a collection of ComboHits into the lowest-level (straw) combo hits. This function is recursive void fillComboHits( std::vector const& indices, CHCIter& iters) const; @@ -171,6 +172,16 @@ namespace mu2e { void setParent(CHCPTR const& parent); // or set to be the same as another collection void setSameParent(ComboHitCollection const& other); + // if argument has its own parent, set this parent to match. + // otherwise set this parent to be collection in argument + void setAsSubset(CHCPTR const& other); + void setAsSubset(art::Handle const& ohandle); + void setAsSubset(art::ValidHandle const& ohandle); + // optionally specify what level to make as parent + // if parent and grandparent are the same level, will select the grandparent unless stopatfirst set + void setAsSubset(CHCPTR const& optr, StrawIdMask::Level level, bool stopatfirst=false); + void setAsSubset(art::Handle const& ohandle, StrawIdMask::Level level, bool stopatfirst=false); + void setAsSubset(art::ValidHandle const& ohandle, StrawIdMask::Level level, bool stopatfirst=false); #endif // accessors auto const& parent() const { return _parent; } diff --git a/RecoDataProducts/src/ComboHit.cc b/RecoDataProducts/src/ComboHit.cc index 6d85a1322..b3cba2f5c 100644 --- a/RecoDataProducts/src/ComboHit.cc +++ b/RecoDataProducts/src/ComboHit.cc @@ -49,17 +49,59 @@ namespace mu2e { #ifndef __ROOTCLING__ - ComboHitCollection::CHCPTR ComboHitCollection::parent(StrawIdMask::Level level) const { + ComboHitCollection::CHCPTR ComboHitCollection::parent(StrawIdMask::Level level, bool stopatfirst) const { auto retval = _parent; - if(_parent.refCore().isNull() || level == this->level()){ + if(_parent.refCore().isNull()){ throw cet::exception("RECO")<<"mu2e::ComboHitCollection: no such parent" << std::endl; } else { // recursive call - if(retval->level() != level) retval = _parent->parent(level); + if (retval->level() != level || (!stopatfirst && retval->parent().refCore().isNonnull() && retval->parent()->level() == level)){ + retval = _parent->parent(level); + } } return retval; } + void ComboHitCollection::setAsSubset(art::Handle const& ohandle){ + if (ohandle.isValid()){ + auto optr = CHCPTR(ohandle); + setAsSubset(optr); + } else { + throw cet::exception("RECO")<<"mu2e::ComboHitCollection: invalid handle" << std::endl; + } + } + void ComboHitCollection::setAsSubset(art::ValidHandle const& ohandle){ + auto optr = CHCPTR(ohandle); + setAsSubset(optr); + } + void ComboHitCollection::setAsSubset(CHCPTR const& other){ + _parent = other->parent(); + if (_parent.refCore().isNull()){ + _parent = other; + } + } + + void ComboHitCollection::setAsSubset(art::ValidHandle const& ohandle, StrawIdMask::Level level, bool stopatfirst){ + if (ohandle.isValid()){ + auto optr = CHCPTR(ohandle); + setAsSubset(optr,level,stopatfirst); + } else { + throw cet::exception("RECO")<<"mu2e::ComboHitCollection: invalid handle" << std::endl; + } + } + void ComboHitCollection::setAsSubset(art::Handle const& ohandle, StrawIdMask::Level level, bool stopatfirst){ + auto optr = CHCPTR(ohandle); + setAsSubset(optr,level,stopatfirst); + } + void ComboHitCollection::setAsSubset(CHCPTR const& optr, StrawIdMask::Level level, bool stopatfirst){ + if (optr->level() == level && (stopatfirst || optr->parent().refCore().isNull() || optr->parent()->level() != level)){ + _parent = optr; + }else{ + _parent = optr->parent(level,stopatfirst); + } + } + + void ComboHitCollection::setSameParent(ComboHitCollection const& other) { _parent = other.parent(); } void ComboHitCollection::setParent(CHCPTR const& parent) { _parent = parent; } @@ -75,10 +117,10 @@ namespace mu2e { _parent = CHCPTR(phandle); } - void ComboHitCollection::fillStrawDigiIndices(size_t chindex, vector& shiv) const { + void ComboHitCollection::fillStrawDigiIndices(size_t chindex, vector& shiv, bool stopatfirst) const { // if this is a straw-level ComboHit, get the digi index ComboHit const& ch = this->at(chindex); - if(level() == StrawIdMask::uniquestraw) { + if(level() == StrawIdMask::uniquestraw && (stopatfirst || _parent.refCore().isNull() || _parent->level() != StrawIdMask::uniquestraw)){ shiv.push_back(ch.index(0)); // if this collection references a sub-collection, go down recursively } else if(_parent.refCore().isNonnull()){ @@ -90,9 +132,9 @@ namespace mu2e { } } - ComboHitCollection const* ComboHitCollection::fillStrawHitIndices( size_t chindex, SHIV& shiv, StrawIdMask::Level clevel) const { + ComboHitCollection const* ComboHitCollection::fillStrawHitIndices( size_t chindex, SHIV& shiv, StrawIdMask::Level clevel, bool stopatfirst) const { ComboHitCollection const* retval = this; - if(level() == clevel){ + if(level() == clevel && (stopatfirst || _parent.refCore().isNull() || _parent->level() != clevel)){ shiv.push_back(chindex); // if this collection references other collections, go down recursively } else if(_parent.refCore().isNonnull()){ @@ -117,9 +159,9 @@ namespace mu2e { return retval; } - ComboHitCollection const* ComboHitCollection::fillStrawHitIndices(SHIV const& inshiv, SHIV& outshiv, StrawIdMask::Level clevel) const { + ComboHitCollection const* ComboHitCollection::fillStrawHitIndices(SHIV const& inshiv, SHIV& outshiv, StrawIdMask::Level clevel, bool stopatfirst) const { ComboHitCollection const* retval = this; - if(level() == clevel){ + if(level() == clevel && (stopatfirst || _parent.refCore().isNull() || _parent->level() != clevel)){ outshiv = inshiv; } else if(_parent.refCore().isNonnull()){ if(_parent->level() == clevel){ @@ -149,10 +191,10 @@ namespace mu2e { return retval; } - ComboHitCollection const* ComboHitCollection::fillStrawHitIndices( SHIV& shiv, StrawIdMask::Level clevel) const { + ComboHitCollection const* ComboHitCollection::fillStrawHitIndices( SHIV& shiv, StrawIdMask::Level clevel, bool stopatfirst) const { shiv.clear(); const ComboHitCollection* retval = this; - if(level() == clevel){ + if(level() == clevel && (stopatfirst || _parent.refCore().isNull() || _parent->level() != clevel)){ shiv.reserve(this->size()); for(size_t iind = 0;iind < this->size(); ++iind)shiv.push_back(iind); // if this collection references other collections, go down recursively