diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
index 392eb9edc43dd2f5bb4aead0087a9abe7e06143c..ff163e105d51b02267ce35e3fdfe43737c5a184c 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.cxx
@@ -20,6 +20,7 @@
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "StoreGate/ReadHandle.h"
+#include <list>
namespace Tracker
{
@@ -52,7 +53,7 @@ StatusCode SegmentFitAlg::initialize() {
}
ATH_CHECK(m_clusterContainerKey.initialize());
-// ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_trackCollection.initialize());
// Get the SCT ID helper
ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
@@ -67,7 +68,7 @@ StatusCode SegmentFitAlg::initialize() {
}
// Process one station
-StatusCode SegmentFitAlg::reconstructStation(ClusterContainerHandle container, int index) const
+StatusCode SegmentFitAlg::reconstructStation(ClusterContainerHandle container, int index, std::unique_ptr<TrackCollection>& tracks) const
{
int station = m_detMgr->numerology().stationId(index);
ATH_MSG_VERBOSE("Processing station " << station);
@@ -119,14 +120,16 @@ StatusCode SegmentFitAlg::reconstructStation(ClusterContainerHandle container, i
// try adding clusters to the good seeds
FitMap allFits = enumerateFits(goodSeeds, badSeeds);
- ATH_MSG_VERBOSE("Found " << allFits.size() << " good Fits after enumeration from " << goodSeeds.size() << " seeds.");
- for (auto f : allFits)
+ ATH_MSG_VERBOSE("Found " << allFits.size() << " fits after enumeration from " << goodSeeds.size() << " seeds.");
+ FitMap goodFits = checkFitQuality(allFits);
+ ATH_MSG_VERBOSE(goodFits.size() << " out of " << allFits.size() << " fits pass quality cut.");
+ for (auto f : goodFits)
{
auto fit = f.second;
if (fit->candidates.size() > 4)
{
ATH_MSG_VERBOSE("Fit Parameters: (" << fit->fitParams(0) << ", " << fit->fitParams(1) << ", " << fit->fitParams(2) << ", " << fit->fitParams(3) << "); chi2 = " << fit->fitChi2 << "; DoF = " << fit->candidates.size() - 4 );
- for (size_t i = fit->clusterMask.find_first(); i != boost::dynamic_bitset<>::npos; i = fit->clusterMask.find_next(i))
+ for (size_t i = fit->clusterMask.find_first(); i != ClusterSet::npos; i = fit->clusterMask.find_next(i))
{
Amg::Vector3D seedPos = l[i]->cluster.globalPosition();
ATH_MSG_VERBOSE("Fit cluster [" << l[i]->layer << "|" << l[i]->eta << "|" << l[i]->phi <<
@@ -134,7 +137,20 @@ StatusCode SegmentFitAlg::reconstructStation(ClusterContainerHandle container, i
}
}
}
-
+ FitMap selectedFits = selectFits(goodFits);
+ ATH_MSG_VERBOSE("Selected " << selectedFits.size() << " fits for output.");
+ for (auto f : selectedFits)
+ {
+ auto fit = f.second;
+ ATH_MSG_VERBOSE("Fit Parameters: (" << fit->fitParams(0) << ", " << fit->fitParams(1) << ", " << fit->fitParams(2) << ", " << fit->fitParams(3) << "); chi2 = " << fit->fitChi2 << "; DoF = " << fit->candidates.size() - 4 );
+ for (size_t i = fit->clusterMask.find_first(); i != ClusterSet::npos; i = fit->clusterMask.find_next(i))
+ {
+ Amg::Vector3D seedPos = l[i]->cluster.globalPosition();
+ ATH_MSG_VERBOSE("Fit cluster [" << l[i]->layer << "|" << l[i]->eta << "|" << l[i]->phi <<
+ "|" << l[i]->side << "] at (" << seedPos.x() << ", " << seedPos.y() << ", " << seedPos.z() << ")");
+ }
+ ATH_CHECK(AddTrack(tracks, fit));
+ }
return StatusCode::SUCCESS;
}
@@ -279,7 +295,7 @@ SegmentFitAlg::enumerateSeeds(const ClusterList& clusters) const
}
auto i_compat = clusters[i]->compatible;
std::vector<size_t> i_vec {i};
- for (size_t j = i_compat.find_next(i); j != boost::dynamic_bitset<>::npos; j = i_compat.find_next(j))
+ for (size_t j = i_compat.find_next(i); j != ClusterSet::npos; j = i_compat.find_next(j))
{
size_t ij_views[2] {i_views[0], i_views[1]};
if (clusters[j]->view > 0)
@@ -291,7 +307,7 @@ SegmentFitAlg::enumerateSeeds(const ClusterList& clusters) const
ij_views[0]++;
}
auto ij_compat = clusters[j]->compatible & i_compat;
- for (size_t k = ij_compat.find_next(j); k != boost::dynamic_bitset<>::npos; k = ij_compat.find_next(k))
+ for (size_t k = ij_compat.find_next(j); k != ClusterSet::npos; k = ij_compat.find_next(k))
{
size_t ijk_views[2] {ij_views[0], ij_views[1]};
if (clusters[k]->view > 0)
@@ -304,7 +320,7 @@ SegmentFitAlg::enumerateSeeds(const ClusterList& clusters) const
}
if (ijk_views[0] == 0 || ijk_views[1] == 0) continue; // Need two hits in each view to fit
auto ijk_compat = clusters[k]->compatible & ij_compat;
- for (size_t l = ijk_compat.find_next(k); l != boost::dynamic_bitset<>::npos; l = ijk_compat.find_next(l))
+ for (size_t l = ijk_compat.find_next(k); l != ClusterSet::npos; l = ijk_compat.find_next(l))
{
size_t ijkl_views[2] {ijk_views[0], ijk_views[1]};
if (clusters[l]->view > 0)
@@ -355,7 +371,7 @@ SegmentFitAlg::findGoodSeeds(FitList& allSeeds, FitMap& goodSeeds, MaskSet& badS
void
SegmentFitAlg::fitClusters(fitInfo& fit) const
{
-
+ m_numberOfFits++;
Eigen::Matrix< double, 5, 5 > s;
s << fit.sums[1] , fit.sums[2] , fit.sums[4] , fit.sums[5] , fit.sums[10] ,
fit.sums[2] , fit.sums[3] , fit.sums[5] , fit.sums[6] , fit.sums[11] ,
@@ -391,6 +407,7 @@ SegmentFitAlg::checkFit(const Tracker::SegmentFitAlg::fitInfo& seed) const
ATH_MSG_VERBOSE("Rejected fit without covariance matrix.");
return false;
}
+ if (m_tanThetaCut > 0.0 && pow(seed.fitParams(2),2.0) + pow(seed.fitParams(3), 2.0) > m_tanThetaCut * m_tanThetaCut) return false;
for (auto i : seed.candidates)
{
double z = seed.clusters[i]->cluster.globalPosition().z();
@@ -420,11 +437,11 @@ SegmentFitAlg::enumerateFits(const FitMap& seeds, const MaskSet& veto) const
for (auto p : newSeeds)
{
auto seed = p.second;
- for (size_t i = seed->compatibleMask.find_first(); i != boost::dynamic_bitset<>::npos; i = seed->compatibleMask.find_next(i))
+ for (size_t i = seed->compatibleMask.find_first(); i != ClusterSet::npos; i = seed->compatibleMask.find_next(i))
{
if (abs(seed->fitParams(1) + seed->fitParams(3) * (seed->clusters[i]->z - fitInfo::zCenter) - seed->clusters[i]->yCenter) > m_yResidualCut)
continue;
- boost::dynamic_bitset<> mask = seed->clusterMask;
+ ClusterSet mask { seed->clusterMask };
mask.set(i);
if (allFits.find(mask) == allFits.end() && !vetoFit(allFailures, mask))
{
@@ -450,13 +467,66 @@ SegmentFitAlg::enumerateFits(const FitMap& seeds, const MaskSet& veto) const
return allFits;
}
-bool
-SegmentFitAlg::vetoFit(const MaskSet& badSeeds, const boost::dynamic_bitset<>& hypothesis) const
+SegmentFitAlg::FitMap
+SegmentFitAlg::checkFitQuality(const FitMap& fits) const
{
- for (auto bad : badSeeds)
+ FitMap goodFits { };
+ for (auto f : fits)
{
- if (bad == hypothesis) return true;
+ auto candidate = f.second;
+ size_t nDoF = candidate->clusterMask.count() - 4;
+ if (nDoF == 0 || candidate->fitChi2/nDoF <= m_reducedChi2Cut) goodFits[candidate->clusterMask] = candidate;
}
+ m_numberOfGoodFits += goodFits.size();
+ return goodFits;
+}
+
+SegmentFitAlg::FitMap
+SegmentFitAlg::selectFits(const FitMap& fits) const
+{
+ // Select the input fit with the most clusters (using lowest chi2 as tie-breaker)
+ // Discard any other input fits using those clusters
+ // Repeat until no input fits remain
+ FitMap selectedFits { };
+ std::list<std::shared_ptr<fitInfo> > info {};
+ for (auto fit : fits)
+ {
+ info.push_back(fit.second);
+ }
+ info.sort([](const std::shared_ptr<fitInfo>& left, const std::shared_ptr<fitInfo>& right)
+ {
+ if (left->clusterMask.count() > right->clusterMask.count()) return true;
+ if (left->clusterMask.count() < right->clusterMask.count()) return false;
+ if (left->fitChi2 < right->fitChi2) return true;
+ return false;
+ });
+
+ ATH_MSG_VERBOSE("Sorted input fits:");
+ for (auto it = info.begin(); it != info.end(); ++it)
+ {
+ auto content = *it;
+ ATH_MSG_VERBOSE("clusters: " << content->clusterMask.count() << " / chi2: " << content->fitChi2) ;
+ }
+ ClusterSet usedClusters { clusterInfo::nClusters };
+
+ while (info.size() > 0)
+ {
+ auto selected = info.front();
+ ATH_MSG_VERBOSE("Selected nclust = " << selected->clusterMask.count() << " with chi2 = " << selected->fitChi2);
+ selectedFits[selected->clusterMask] = selected;
+ usedClusters |= selected->clusterMask;
+
+ info.remove_if([&](std::shared_ptr<fitInfo> p) {return (p->clusterMask & ~usedClusters) != p->clusterMask;});
+ }
+
+ m_numberOfSelectedFits += selectedFits.size();
+ return selectedFits;
+}
+
+bool
+SegmentFitAlg::vetoFit(const MaskSet& badSeeds, const ClusterSet& hypothesis) const
+{
+ if (badSeeds.count(hypothesis) > 0) return true;
return false;
}
@@ -468,8 +538,8 @@ SegmentFitAlg::execute(const EventContext& ctx) const
++m_numberOfEvents;
// setFilterPassed(false, ctx);
-// SG::WriteHandle trackContainer{m_trackCollection, ctx};
-// std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+ SG::WriteHandle trackContainer{m_trackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
ClusterContainerHandle clusterContainer{m_clusterContainerKey, ctx};
ATH_CHECK(clusterContainer.isValid());
@@ -477,7 +547,7 @@ SegmentFitAlg::execute(const EventContext& ctx) const
for (int stationIdx = 0; stationIdx < m_detMgr->numerology().numStations(); stationIdx++)
{
- ATH_CHECK(reconstructStation(clusterContainer, stationIdx));
+ ATH_CHECK(reconstructStation(clusterContainer, stationIdx, outputTracks));
}
// Loop over stations
@@ -530,7 +600,7 @@ SegmentFitAlg::execute(const EventContext& ctx) const
// ATH_MSG_DEBUG("No stations with clusters found.");
// }
-// ATH_CHECK(trackContainer.record(std::move(outputTracks)));
+ ATH_CHECK(trackContainer.record(std::move(outputTracks)));
// // Clean up
// for (auto entry : clusterMap)
@@ -569,77 +639,79 @@ StatusCode SegmentFitAlg::finalize()
ATH_MSG_INFO( m_numberOfSeeds << " seeds fit.");
ATH_MSG_INFO( m_numberOfGoodSeeds << " good seed fits.");
ATH_MSG_INFO( m_numberOfFits << " fits performed." );
+ ATH_MSG_INFO( m_numberOfGoodFits << " fits passed quality cuts." );
+ ATH_MSG_INFO( m_numberOfSelectedFits << " fits selected for output." );
return StatusCode::SUCCESS;
}
-// // Method to create and store Trk::Track from fit results
-// StatusCode
-// SegmentFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
-// const Eigen::Matrix< double, 4, 1 >& fitResult,
-// const Eigen::Matrix< double, 4, 4 >& fitCovariance,
-// std::vector<const clusterInfo*>& fitClusters,
-// double chi2, int ndof) const
-// {
-// Trk::TrackInfo i { Trk::TrackInfo::TrackFitter::Unknown, Trk::ParticleHypothesis::muon };
-// Trk::FitQuality* q = new Trk::FitQuality {chi2, ndof};
-// DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
+// Method to create and store Trk::Track from fit results
+StatusCode
+SegmentFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
+ const std::shared_ptr<fitInfo>& theFit) const
+{
+ Trk::TrackInfo i { Trk::TrackInfo::TrackFitter::Unknown, Trk::ParticleHypothesis::muon };
+ Trk::FitQuality* q = new Trk::FitQuality {theFit->fitChi2, static_cast<double>(theFit->clusterMask.count() - 4)};
+ DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
-// int station = m_idHelper->station(fitClusters[0]->cluster->detectorElement()->identify());
+ // translate parameters to nominal fit point
+ s->push_back( GetState(theFit->fitParams, theFit->fitCovariance, nullptr) );
-// // translate parameters to nominal fit point
-// s->push_back( GetState(fitResult, fitCovariance, nullptr, station) );
+ for (int clusterIndex : theFit->candidates)
+ {
+ s->push_back( GetState( theFit->fitParams, theFit->fitCovariance, &theFit->clusters[clusterIndex]->cluster) );
+ }
-// for (const clusterInfo* cInfo : fitClusters)
-// {
-// s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
-// }
+ // for (const clusterInfo* cInfo : fitClusters)
+ // {
+ // s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
+ // }
-// // Create and store track
-// std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink(s);
-// tracks->push_back(new Trk::Track(i, std::move(*sink) , q));
-// return StatusCode::SUCCESS;
-// }
+ // Create and store track
+ std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink(s);
+ tracks->push_back(new Trk::Track(i, std::move(*sink) , q));
+ return StatusCode::SUCCESS;
+}
-// Trk::TrackStateOnSurface*
-// SegmentFitAlg::GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
-// const Eigen::Matrix< double, 4, 4 >& fitCovariance,
-// const FaserSCT_Cluster* fitCluster, int station ) const
-// {
-// // position of fit point:
-// // int station = m_idHelper->station(fitCluster->detectorElement()->identify());
-// double zFit = m_zCenter[station];
-// if (fitCluster != nullptr) zFit = fitCluster->globalPosition()[2];
-// Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - m_zCenter[station]), fitResult[1] + fitResult[3] * (zFit - m_zCenter[station]), zFit };
-// double phi = atan2( fitResult[3], fitResult[2] );
-// double theta = atan( sqrt(fitResult[2]*fitResult[2] + fitResult[3]*fitResult[3]) );
-// double qoverp = 1.0/100000.0;
-
-// Eigen::Matrix< double, 4, 4 > jacobian = Eigen::Matrix< double, 4, 4 >::Zero();
-// jacobian << 1, 0, zFit - m_zCenter[station], 0,
-// 0, 1, 0, zFit - m_zCenter[station],
-// 0, 0, fitResult[3]/(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]),
-// fitResult[2]/(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]),
-// 0, 0, fitResult[2]/(sqrt(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])*(1+fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])),
-// fitResult[3]/(sqrt(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])*(1+fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]));
-// Eigen::Matrix< double, 4, 4 > covPar {jacobian * fitCovariance * jacobian.transpose()};
-// std::unique_ptr<AmgSymMatrix(5)> covPar5 { new AmgSymMatrix(5){ AmgSymMatrix(5)::Zero() } };
-// covPar5->block<4,4>(0,0) = covPar;
-// covPar5->block<1,1>(4,4) = Eigen::Matrix< double, 1, 1 > { (50000.0)*qoverp*qoverp };
-
-// FaserSCT_ClusterOnTrack* rot = nullptr;
-// if (fitCluster != nullptr)
-// {
-// rot = new FaserSCT_ClusterOnTrack{ fitCluster,
-// Trk::LocalParameters { Trk::DefinedParameter { fitCluster->localPosition()[0], Trk::loc1 },
-// Trk::DefinedParameter { fitCluster->localPosition()[1], Trk::loc2 } },
-// fitCluster->localCovariance(),
-// m_idHelper->wafer_hash(fitCluster->detectorElement()->identify())};
-// }
-// std::unique_ptr<Trk::TrackParameters> p { new Trk::CurvilinearParameters { pos, phi, theta, qoverp, *(covPar5.release()) } };
-// return new Trk::TrackStateOnSurface { rot, p.release() };
-// }
+Trk::TrackStateOnSurface*
+SegmentFitAlg::GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
+ const Eigen::Matrix< double, 4, 4 >& fitCovariance,
+ const FaserSCT_Cluster* fitCluster) const
+{
+ // position of fit point:
+ // int station = m_idHelper->station(fitCluster->detectorElement()->identify());
+ double zFit = fitInfo::zCenter;
+ if (fitCluster != nullptr) zFit = fitCluster->globalPosition()[2];
+ Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - fitInfo::zCenter), fitResult[1] + fitResult[3] * (zFit - fitInfo::zCenter), zFit };
+ double phi = atan2( fitResult[3], fitResult[2] );
+ double theta = atan( sqrt(fitResult[2]*fitResult[2] + fitResult[3]*fitResult[3]) );
+ double qoverp = 1.0/100000.0;
+
+ Eigen::Matrix< double, 4, 4 > jacobian = Eigen::Matrix< double, 4, 4 >::Zero();
+ jacobian << 1, 0, zFit - fitInfo::zCenter, 0,
+ 0, 1, 0, zFit - fitInfo::zCenter,
+ 0, 0, fitResult[3]/(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]),
+ fitResult[2]/(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]),
+ 0, 0, fitResult[2]/(sqrt(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])*(1+fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])),
+ fitResult[3]/(sqrt(fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3])*(1+fitResult[2]*fitResult[2]+fitResult[3]*fitResult[3]));
+ Eigen::Matrix< double, 4, 4 > covPar {jacobian * fitCovariance * jacobian.transpose()};
+ std::unique_ptr<AmgSymMatrix(5)> covPar5 { new AmgSymMatrix(5){ AmgSymMatrix(5)::Zero() } };
+ covPar5->block<4,4>(0,0) = covPar;
+ covPar5->block<1,1>(4,4) = Eigen::Matrix< double, 1, 1 > { (50000.0)*qoverp*qoverp };
+
+ FaserSCT_ClusterOnTrack* rot = nullptr;
+ if (fitCluster != nullptr)
+ {
+ rot = new FaserSCT_ClusterOnTrack{ fitCluster,
+ Trk::LocalParameters { Trk::DefinedParameter { fitCluster->localPosition()[0], Trk::loc1 },
+ Trk::DefinedParameter { fitCluster->localPosition()[1], Trk::loc2 } },
+ fitCluster->localCovariance(),
+ m_idHelper->wafer_hash(fitCluster->detectorElement()->identify())};
+ }
+ std::unique_ptr<Trk::TrackParameters> p { new Trk::CurvilinearParameters { pos, phi, theta, qoverp, *(covPar5.release()) } };
+ return new Trk::TrackStateOnSurface { rot, p.release() };
+}
// void
// SegmentFitAlg::Residuals(std::vector<const clusterInfo*>& fitClusters) const
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
index 29f184a33d35f008d89e2c12849e1616f64cbcb5..eb41946cd0a9094a24919eaf6824705094afa234 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/src/SegmentFitAlg.h
@@ -47,6 +47,7 @@ namespace TrackerDD
namespace Tracker
{
+ typedef boost::dynamic_bitset<> ClusterSet;
/**
* @class SCT_Clusterization
@@ -72,6 +73,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
virtual bool isClonable() const override { return true; };
//@}
+
private:
/** @name Disallow default instantiation, copy, assignment */
//@{
@@ -186,7 +188,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
static size_t nClusters;
uint index;
const FaserSCT_Cluster& cluster;
- boost::dynamic_bitset<> compatible;
+ ClusterSet compatible;
double sinAlpha;
double cosAlpha;
double u;
@@ -223,6 +225,21 @@ class SegmentFitAlg : public AthReentrantAlgorithm
addCluster(index);
}
}
+ // fitInfo(const fitInfo& f)
+ // : clusters { f.clusters }
+ // , candidates { f.candidates }
+ // , clusterMask { f.clusterMask }
+ // , compatibleMask { f.compatibleMask }
+ // , fitParams { f.fitParams }
+ // , fitCovariance { f.fitCovariance }
+ // , fitChi2 { f.fitChi2 }
+ // , hasCovariance { f.hasCovariance }
+ // {
+ // for (size_t i = 0; i < 15; i++)
+ // {
+ // sums[i] = f.sums[i];
+ // }
+ // }
void addCluster(size_t n)
{
candidates.push_back(n);
@@ -255,7 +272,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
fitCovariance.setZero();
fitChi2 = 0;
hasCovariance = false;
- }
+ }
}
HepGeom::Point3D<double> extrapolateFit(double z) const
{
@@ -263,8 +280,8 @@ class SegmentFitAlg : public AthReentrantAlgorithm
}
const ClusterList& clusters;
std::vector<size_t> candidates;
- boost::dynamic_bitset<> clusterMask;
- boost::dynamic_bitset<> compatibleMask;
+ ClusterSet clusterMask;
+ ClusterSet compatibleMask;
double sums[15];
static double zCenter;
Eigen::Matrix<double, 4, 1> fitParams;
@@ -272,13 +289,14 @@ class SegmentFitAlg : public AthReentrantAlgorithm
double fitChi2;
bool hasCovariance;
};
+
typedef std::vector<std::shared_ptr<fitInfo> > FitList;
- typedef std::unordered_map<boost::dynamic_bitset<>, std::shared_ptr<fitInfo> > FitMap;
- typedef std::unordered_set<boost::dynamic_bitset<> > MaskSet;
+ typedef std::unordered_map<ClusterSet, std::shared_ptr<fitInfo> > FitMap;
+ typedef std::unordered_set<ClusterSet> MaskSet;
typedef SG::ReadHandle<FaserSCT_ClusterContainer> ClusterContainerHandle;
- StatusCode reconstructStation(ClusterContainerHandle clusters, int index) const;
+ StatusCode reconstructStation(ClusterContainerHandle clusters, int index, std::unique_ptr<TrackCollection>& tracks) const;
size_t stationOccupancy(ClusterContainerHandle clusters, int station) const;
double findCenter(int station) const;
ClusterList tabulateClusters(ClusterContainerHandle clusters, int station) const;
@@ -288,7 +306,9 @@ class SegmentFitAlg : public AthReentrantAlgorithm
void fitClusters(fitInfo& fit) const;
bool checkFit(const fitInfo& fit) const;
FitMap enumerateFits(const FitMap& seeds, const MaskSet& veto) const;
- bool vetoFit(const MaskSet& badSeeds, const boost::dynamic_bitset<>& hypothesis) const;
+ FitMap checkFitQuality(const FitMap& fits) const;
+ FitMap selectFits(const FitMap& fits) const;
+ bool vetoFit(const MaskSet& badSeeds, const ClusterSet& hypothesis) const;
int getStation(const Tracker::FaserSCT_ClusterCollection* collection) const;
// struct layerCombo
@@ -369,16 +389,12 @@ class SegmentFitAlg : public AthReentrantAlgorithm
// }
// };
- // StatusCode AddTrack(std::unique_ptr<TrackCollection>& tracks,
- // const Eigen::Matrix< double, 4, 1 >& fitResult,
- // const Eigen::Matrix< double, 4, 4>& fitCovariance,
- // std::vector<const clusterInfo*>& fitClusters,
- // double chi2, int ndof) const;
+ StatusCode AddTrack(std::unique_ptr<TrackCollection>& tracks,
+ const std::shared_ptr<fitInfo>& theFit) const;
- // Trk::TrackStateOnSurface* GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
- // const Eigen::Matrix< double, 4, 4 >& fitCovariance,
- // const FaserSCT_Cluster* fitCluster,
- // int station ) const;
+ Trk::TrackStateOnSurface* GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
+ const Eigen::Matrix< double, 4, 4 >& fitCovariance,
+ const FaserSCT_Cluster* fitCluster ) const;
// void
// Residuals(std::vector<const clusterInfo*>& fitClusters) const;
@@ -387,7 +403,7 @@ class SegmentFitAlg : public AthReentrantAlgorithm
const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
SG::ReadHandleKey<FaserSCT_ClusterContainer> m_clusterContainerKey { this, "ClustersName", "SCT_ClusterContainer", "FaserSCT cluster container" };
- // SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "ClusterFit", "Output track collection name" };
+ SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "SegmentFit", "Output track collection name" };
UnsignedIntegerProperty m_minClustersPerStation { this, "MinClusters", 4, "Minimum number of clusters allowed in a station" };
UnsignedIntegerProperty m_maxClustersPerStation { this, "MaxClusters", 36, "Maximum number of clusters allowed in a station" };
@@ -396,6 +412,8 @@ class SegmentFitAlg : public AthReentrantAlgorithm
IntegerProperty m_edgeStrips { this, "EdgeStrips", 32, "Number of strips from either edge of sensor considered for overlap" };
DoubleProperty m_waferTolerance { this, "WaferTolerance", 3 * Gaudi::Units::mm, "Wafer boundary tolerance for extrapolated fit." };
DoubleProperty m_yResidualCut { this, "ResidualCut", 3 * Gaudi::Units::mm, "Cluster y tolerance for compatibility with extrapolated fit." };
+ DoubleProperty m_tanThetaCut { this, "TanThetaCut", 0.0, "Maximum accepted tangent of track angle from beam axis; 0 means no cut."};
+ DoubleProperty m_reducedChi2Cut { this, "ReducedChi2Cut", 10.0, "Maximum accepted chi^2 per degree of freedom for final fits; 0 means no cut." };
mutable std::atomic<int> m_numberOfEvents{0};
mutable std::atomic<int> m_numberExcessOccupancy{0};
@@ -410,6 +428,8 @@ class SegmentFitAlg : public AthReentrantAlgorithm
mutable std::atomic<int> m_numberOfSeeds{0};
mutable std::atomic<int> m_numberOfGoodSeeds{0};
mutable std::atomic<int> m_numberOfFits{0};
+ mutable std::atomic<int> m_numberOfGoodFits{0};
+ mutable std::atomic<int> m_numberOfSelectedFits{0};
};
inline int SegmentFitAlg::getStation(const Tracker::FaserSCT_ClusterCollection* collection) const
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/test/TestBeamSegmentFitDbg.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/test/TestBeamSegmentFitDbg.py
index 450b6313f2bd6bebf115f5bda613e28cec15a153..c4fcd2a24f74dee7e00cf1d1317882f088bbbf20 100755
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/test/TestBeamSegmentFitDbg.py
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/test/TestBeamSegmentFitDbg.py
@@ -14,7 +14,7 @@ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
#from Digitization.DigitizationParametersConfig import writeDigitizationMetadata
-from ScintRecAlgs.ScintRecAlgsConfig import WaveformReconstructionCfg
+from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg
from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
@@ -26,10 +26,10 @@ Configurable.configurableRun3Behavior = True
# Configure
ConfigFlags.Input.Files = [
- 'tb.raw',
+ 'tbMu.raw',
]
-ConfigFlags.Output.ESDFileName = "tb.ESD.pool.root"
-ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
+ConfigFlags.Output.ESDFileName = "tbMu.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00" # Always needed; must match FaserVersion
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
@@ -38,7 +38,7 @@ ConfigFlags.Common.isOnline = False
ConfigFlags.GeoModel.Align.Dynamic = False
ConfigFlags.Beam.NumberOfCollisions = 0.
-ConfigFlags.Detector.GeometryFaserSCT = True
+# ConfigFlags.Detector.GeometryFaserSCT = True
ConfigFlags.lock()
@@ -53,34 +53,64 @@ from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnv
acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
acc.merge(WaveformReconstructionCfg(ConfigFlags))
acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_EDGEMODE_RDOs"))
-acc.merge(SegmentFitAlgCfg(ConfigFlags))
+acc.merge(SegmentFitAlgCfg(ConfigFlags, TanThetaCut = 0.0, MaxClusters = 25))
acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
acc.getEventAlgo("Tracker::SegmentFitAlg").OutputLevel = VERBOSE
-from AthenaConfiguration.ComponentFactory import CompFactory
-decoderTool = CompFactory.ScintWaveformDecoderTool(name = "ScintWaveformDecoderTool",
- CaloChannels = [0, 1, 2, 3, 4, 5],
- PreshowerChannels = [6, 7],
- TriggerChannels = [8, 9],
- VetoChannels=[])
-acc.addPublicTool(decoderTool)
+# from AthenaConfiguration.ComponentFactory import CompFactory
+# decoderTool = CompFactory.ScintWaveformDecoderTool(name = "ScintWaveformDecoderTool",
+# CaloChannels = [0, 1, 2, 3, 4, 5],
+# PreshowerChannels = [6, 7],
+# TriggerChannels = [8, 9],
+# VetoChannels=[])
+# acc.addPublicTool(decoderTool)
-# explicitly save RDO information
+# # explicitly save RDO information
+# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+# itemList = [ "xAOD::EventInfo#*",
+# "xAOD::EventAuxInfo#*",
+# "FaserSCT_RDO_Container#*",
+# "xAOD::FaserTriggerData#*",
+# "xAOD::FaserTriggerDataAux#*",
+# "ScintWaveformContainer#*",
+# "TrackCollection#*",
+# "xAOD::WaveformHitContainer#*",
+# "xAOD::WaveformHitAuxContainer#*",
+# "xAOD::WaveformClock#*",
+# "xAOD::WaveformClockAuxInfo#*",
+# ]
+# acc.merge(OutputStreamCfg(ConfigFlags, "ESD", itemList))
+# acc.getEventAlgo("OutputStreamESD").AcceptAlgs = ["Tracker::SegmentFitAlg"]
+
+#
+# Configure output
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
-itemList = [ "xAOD::EventInfo#*",
- "xAOD::EventAuxInfo#*",
- "FaserSCT_RDO_Container#*",
- "xAOD::FaserTriggerData#*",
- "xAOD::FaserTriggerDataAux#*",
- "ScintWaveformContainer#*",
- "TrackCollection#*",
- "xAOD::WaveformHitContainer#*",
- "xAOD::WaveformHitAuxContainer#*",
- "xAOD::WaveformClock#*",
- "xAOD::WaveformClockAuxInfo#*",
- ]
-acc.merge(OutputStreamCfg(ConfigFlags, "ESD", itemList))
-acc.getEventAlgo("OutputStreamESD").AcceptAlgs = ["Tracker::SegmentFitAlg"]
+itemList = [ "xAOD::EventInfo#*"
+ , "xAOD::EventAuxInfo#*"
+ , "xAOD::FaserTriggerData#*"
+ , "xAOD::FaserTriggerDataAux#*"
+ , "FaserSCT_RDO_Container#*"
+ , "Tracker::FaserSCT_ClusterContainer#*"
+ #, "Tracker::SCT_SpacePointContainer#*"
+ #, "Tracker::SCT_SpacePointOverlapCollection#*"
+ , "TrackCollection#*"
+]
+acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
+
+# Waveform reconstruction
+from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
+acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
+
+
+
+# Hack to avoid problem with our use of MC databases when isMC = False
+replicaSvc = acc.getService("DBReplicaSvc")
+replicaSvc.COOLSQLiteVetoPattern = ""
+replicaSvc.UseCOOLSQLite = True
+replicaSvc.UseCOOLFrontier = False
+replicaSvc.UseGeomSQLite = True
+
+
# Timing
#acc.merge(MergeRecoTimingObjCfg(ConfigFlags))