Skip to content

Commit

Permalink
Merge pull request #1310 from bonventre/combohit
Browse files Browse the repository at this point in the history 
Allow ComboHitCollections to have a parent at the same level
  • Loading branch information
@kutschke
kutschke committed Aug 7, 2024
2 parents 3661cc3 + f9d98e0 commit 80a19a0
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 60 deletions.
72 changes: 28 additions & 44 deletions CosmicReco/src/LineFinder_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ namespace mu2e{

ProditionsHandle<Tracker> _alignedTracker_h;

int findLine(const ComboHitCollection& shC, art::Event const& event, CosmicTrackSeed &tseed);
int findLine(const ComboHitCollection& shC, std::vector<StrawHitIndex> const& shiv, art::Event const& event, CosmicTrackSeed& tseed);
};


Expand Down Expand Up @@ -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<StrawHitIndex> 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<TimeCluster>(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;

Expand All @@ -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<StrawHitIndex> const& shiv, art::Event const& event, CosmicTrackSeed& tseed){

//mu2e::GeomHandle<mu2e::Tracker> th;
auto tracker = _alignedTracker_h.getPtr(event.id());//th.get();
Expand All @@ -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;i<shC.size();i++){
for (size_t i=0;i<shiv.size();i++){
size_t iloc = shiv[i];
if (found_all)
break;
Straw const& strawi = tracker->getStraw(shC[i].strawId());
for (size_t j=i+1;j<shC.size();j++){
Straw const& strawi = tracker->getStraw(shC[iloc].strawId());
for (size_t j=i+1;j<shiv.size();j++){
size_t jloc = shiv[j];
if (found_all)
break;
Straw const& strawj = tracker->getStraw(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();
Expand All @@ -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;k<shC.size();k++){
Straw const& strawk = tracker->getStraw(shC[k].strawId());
for (size_t k=0;k<shiv.size();k++){
size_t kloc = shiv[k];
Straw const& strawk = tracker->getStraw(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;
Expand All @@ -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;k<shC.size();k++){
Straw const& strawk = tracker->getStraw(shC[k].strawId());
for (size_t k=0;k<shiv.size();k++){
size_t kloc = shiv[k];
Straw const& strawk = tracker->getStraw(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;
Expand Down Expand Up @@ -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<tseed._straw_chits.size(); ich++){
std::vector<StrawHitIndex> 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;
}

Expand Down
21 changes: 16 additions & 5 deletions RecoDataProducts/inc/ComboHit.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<uint16_t> const& indices, CHCIter& iters) const;
Expand All @@ -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<ComboHitCollection> const& ohandle);
void setAsSubset(art::ValidHandle<ComboHitCollection> 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<ComboHitCollection> const& ohandle, StrawIdMask::Level level, bool stopatfirst=false);
void setAsSubset(art::ValidHandle<ComboHitCollection> const& ohandle, StrawIdMask::Level level, bool stopatfirst=false);
#endif
// accessors
auto const& parent() const { return _parent; }
Expand Down
64 changes: 53 additions & 11 deletions RecoDataProducts/src/ComboHit.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<ComboHitCollection> 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<ComboHitCollection> 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<ComboHitCollection> 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<ComboHitCollection> 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; }

Expand All @@ -75,10 +117,10 @@ namespace mu2e {
_parent = CHCPTR(phandle);
}

void ComboHitCollection::fillStrawDigiIndices(size_t chindex, vector<StrawDigiIndex>& shiv) const {
void ComboHitCollection::fillStrawDigiIndices(size_t chindex, vector<StrawDigiIndex>& 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()){
Expand All @@ -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()){
Expand All @@ -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){
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 80a19a0

Please sign in to comment.