Skip to content

exclude region feature #384

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
miquelnebot merged 2 commits into from
Oct 2, 2023
Merged

exclude region feature #384

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
81291bc
exclude region feature
Sep 15, 2023
daab112
restructure, skip out of fv hits
Sep 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 96 additions & 6 deletions sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
#include "TrajectoryMCSFitter.h"
#include "lardataobj/RecoBase/Track.h"
#include "larcorealg/Geometry/geo_vectors_utils.h"
#include "larcore/Geometry/Geometry.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"

Expand All @@ -15,7 +18,8 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
vector<size_t> breakpoints;
vector<float> segradlengths;
vector<float> cumseglens;
breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
vector<bool> breakpointsgood;
breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens, breakpointsgood);
//
// Fit segment directions, and get 3D angles between them
//
Expand All @@ -28,7 +32,11 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
if (p>0) {
if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
dtheta.push_back(-999.);
} else {
}
if (!breakpointsgood[p] || !breakpointsgood[p-1]) {
dtheta.push_back(-999.);
}
else {
const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
//assert(std::abs(cosval)<=1);
//units are mrad
Expand All @@ -49,15 +57,21 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
}
const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true, momDepConst, pid);
const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid);
//
// std::cout << "fwdResult.p=" << fwdResult.p << " fwdResult.pUnc=" << fwdResult.pUnc << " fwdResult.logL=" << fwdResult.logL << std::endl;
return recob::MCSFitResult(pid,
fwdResult.p,fwdResult.pUnc,fwdResult.logL,
bwdResult.p,bwdResult.pUnc,bwdResult.logL,
segradlengths,dtheta);
}

void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens) const {
//
void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens, vector<bool>& breakpointsgood) const {
// set fiducial volume
std::vector<geo::BoxBoundedGeo> fiducialVolumes;
fiducialVolumes = setFiducialVolumes();
// set excluded volumes
std::vector<geo::BoxBoundedGeo> excludeVolumes;
excludeVolumes = setExcludeVolumes();

const double trajlen = traj.Length();
const int nseg = std::max(minNSegs_,int(trajlen/segLen_));
const double thisSegLen = trajlen/double(nseg);
Expand All @@ -69,6 +83,8 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
auto nextValid=traj.FirstValidPoint();
breakpoints.push_back(nextValid);
auto pos0 = traj.LocationAtPoint(nextValid);
bool pos0good = isInVolume(fiducialVolumes, pos0) && !isInVolume(excludeVolumes, pos0);
breakpointsgood.push_back(pos0good);
nextValid = traj.NextValidPoint(nextValid+1);
int npoints = 0;
while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
Expand All @@ -77,7 +93,9 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
pos0=pos1;
npoints++;
if (thislen>=thisSegLen) {
bool pos1good = isInVolume(fiducialVolumes, pos1) && !isInVolume(excludeVolumes, pos1);
breakpoints.push_back(nextValid);
breakpointsgood.push_back(pos1good);
if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
else segradlengths.push_back(-999.);
cumseglens.push_back(cumseglens.back()+thislen);
Expand All @@ -88,7 +106,10 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
}
//then add last segment
if (thislen>0.) {
auto endpointpos = traj.LocationAtPoint(nextValid);
bool endpointposgood = isInVolume(fiducialVolumes, endpointpos) && !isInVolume(excludeVolumes, endpointpos);
breakpoints.push_back(traj.LastValidPoint()+1);
breakpointsgood.push_back(endpointposgood);
Comment on lines +109 to +112
Copy link
Member

@PetrilloAtWork PetrilloAtWork May 3, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question: isn't nextValid here always an recob::TrackTrajectory::InvalidIndex?
The behaviour of LocationAtPoint() on an InvalidIndex is undefined, so we should make sure we are not triggering it.

segradlengths.push_back(thislen*lar_radl_inv);
cumseglens.push_back(cumseglens.back()+thislen);
}
Expand Down Expand Up @@ -202,7 +223,6 @@ double TrajectoryMCSFitter::mcsLikelihood(double p, double theta0x, std::vector<
double result = 0;
for (int i = beg; i != end; i+=incr ) {
if (dthetaij[i]<0) {
//cout << "skip segment with too few points" << endl;
continue;
}
//
Expand Down Expand Up @@ -302,3 +322,73 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra
}
return current_E;
}

std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setFiducialVolumes() const {
std::vector<geo::BoxBoundedGeo> fiducialVolumes;
std::vector<std::vector<geo::BoxBoundedGeo>> TPCVolumes;
std::vector<geo::BoxBoundedGeo> ActiveVolumes;

const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
for (auto const &cryoID: geometry->Iterate<geo::CryostatID>()) {
std::vector<geo::BoxBoundedGeo> thisTPCVolumes;
for (auto const& tpc: geometry->Iterate<geo::TPCGeo>(cryoID)) {
thisTPCVolumes.push_back(tpc.ActiveBoundingBox());
}
TPCVolumes.push_back(std::move(thisTPCVolumes));
}
for (const std::vector<geo::BoxBoundedGeo> &tpcs: TPCVolumes) {
double xMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
double xMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
double yMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
double yMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
double zMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
double zMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
ActiveVolumes.emplace_back(xMin, xMax, yMin, yMax, zMin, zMax);
}
double fidInsetMinX = fiducialVolumeInsets_[0];
double fidInsetMaxX = fiducialVolumeInsets_[1];
double fidInsetMinY = fiducialVolumeInsets_[2];
double fidInsetMaxY = fiducialVolumeInsets_[3];
double fidInsetMinZ = fiducialVolumeInsets_[4];
double fidInsetMaxZ = fiducialVolumeInsets_[5];

if (fiducialVolumeInsets_.size() != 6) {
std::cout << "Error: fiducialVolumeInsets vector must have length of 6, not fiducializing" << std::endl;
fidInsetMinX = 0.0;
fidInsetMaxX = 0.0;
fidInsetMinY = 0.0;
fidInsetMaxY = 0.0;
fidInsetMinZ = 0.0;
fidInsetMaxZ = 0.0;
}
for (const geo::BoxBoundedGeo &AV: ActiveVolumes) {
fiducialVolumes.emplace_back(AV.MinX() + fidInsetMinX, AV.MaxX() - fidInsetMaxX,
AV.MinY() + fidInsetMinY, AV.MaxY() - fidInsetMaxY,
AV.MinZ() + fidInsetMinZ, AV.MaxZ() - fidInsetMaxZ);
}
return fiducialVolumes;
}

std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setExcludeVolumes() const {
std::vector<geo::BoxBoundedGeo> excludeVolumes;
if (excludeVolumes_.size()%6 != 0) {
std::cout << "Error: excludeVolumes vector must have length multiple of 6, not excluding any regions" << std::endl;
excludeVolumes.emplace_back(-9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0);
return excludeVolumes;
}
for (unsigned int i=0; i<excludeVolumes_.size()/6; i++) {
excludeVolumes.emplace_back(excludeVolumes_[6*i+0], excludeVolumes_[6*i+1], excludeVolumes_[6*i+2], excludeVolumes_[6*i+3], excludeVolumes_[6*i+4], excludeVolumes_[6*i+5]);
}
return excludeVolumes;
}

bool TrajectoryMCSFitter::isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const {
for (const geo::BoxBoundedGeo &volume: volumes) {
if (point.X()>=volume.MinX() && point.X()<=volume.MaxX() &&
point.Y()>=volume.MinY() && point.Y()<=volume.MaxY() &&
point.Z()>=volume.MinZ() && point.Z()<=volume.MaxZ()) {
return true;
}
}
return false;
}
25 changes: 22 additions & 3 deletions sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,15 @@

#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Sequence.h"
#include "fhiclcpp/types/Table.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "lardataobj/RecoBase/MCSFitResult.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardata/RecoObjects/TrackState.h"
#include "larcore/Geometry/Geometry.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcorealg/Geometry/GeometryCore.h"

namespace trkf::sbn {
/**
Expand Down Expand Up @@ -87,10 +91,18 @@ namespace trkf::sbn {
Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."),
3.0
};
fhicl::Sequence<double> fiducialVolumeInsets {
Name("fiducialVolumeInsets"),
Comment("insets from the TPC boundaries to exclude from the fit")
};
fhicl::Sequence<double> excludeVolumes {
Name("excludeVolumes"),
Comment("regions to exclude")
};
};
using Parameters = fhicl::Table<Config>;
//
TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol){
TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vector<double>fiducialVolumeInsets, std::vector<double> excludeVolumes){
pIdHyp_ = pIdHyp;
minNSegs_ = minNSegs;
segLen_ = segLen;
Expand All @@ -101,9 +113,11 @@ namespace trkf::sbn {
pMax_ = pMax;
pStep_ = pStep;
angResol_ = angResol;
fiducialVolumeInsets_ = fiducialVolumeInsets;
excludeVolumes_ = excludeVolumes;
}
explicit TrajectoryMCSFitter(const Parameters & p)
: TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {}
: TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().fiducialVolumeInsets(),p().excludeVolumes()) {}
//
recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, bool momDepConst = true) const { return fitMcs(traj,pIdHyp_,momDepConst); }
recob::MCSFitResult fitMcs(const recob::Track& track, bool momDepConst = true) const { return fitMcs(track,pIdHyp_,momDepConst); }
Expand All @@ -117,7 +131,7 @@ namespace trkf::sbn {
return fitMcs(tt,pid,momDepConst);
}
//
void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens) const;
void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens, std::vector<bool>& breakpointsgood) const;
void linearRegression(const recob::TrackTrajectory& traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t& pcdir) const;
double mcsLikelihood(double p, double theta0x, std::vector<float>& dthetaij, std::vector<float>& seg_nradl, std::vector<float>& cumLen, bool fwd, bool momDepConst, int pid) const;
//
Expand Down Expand Up @@ -146,6 +160,9 @@ namespace trkf::sbn {
double energyLossLandau(const double mass2,const double E2, const double x) const;
//
double GetE(const double initial_E, const double length_travelled, const double mass) const;
std::vector<geo::BoxBoundedGeo> setFiducialVolumes() const;
std::vector<geo::BoxBoundedGeo> setExcludeVolumes() const;
bool isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const;
//
private:
int pIdHyp_;
Expand All @@ -158,6 +175,8 @@ namespace trkf::sbn {
double pMax_;
double pStep_;
double angResol_;
std::vector<double> fiducialVolumeInsets_;
std::vector<double> excludeVolumes_;
};
}

Expand Down
5 changes: 4 additions & 1 deletion sbncode/LArRecoProducer/mcsproducer.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,10 @@
BEGIN_PROLOG
mcs_sbn: {
module_type: MCSFitAllPID
MCS: {}
MCS: {
fiducialVolumeInsets: []
excludeVolumes: []
}
TrackLabel: pandoraTrack
}

Expand Down
13 changes: 13 additions & 0 deletions sbncode/LArRecoProducer/mcsproducer_icarus.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
BEGIN_PROLOG
mcs_sbn: {
module_type: MCSFitAllPID
MCS: {
eLossMode: 2
angResol: 10.0
fiducialVolumeInsets: [10.0, 10.0, 10.0, 10.0, 10.0, 10.0]
excludeVolumes: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0]
}
TrackLabel: pandoraTrack
}

END_PROLOG