Skip to content

Commit

Permalink
set noise hit mcParticleIndex to -1
Browse files Browse the repository at this point in the history
  • Loading branch information
XiaocongAi committed Aug 25, 2024
1 parent c164938 commit 3b25620
Show file tree
Hide file tree
Showing 9 changed files with 54 additions and 51 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ struct SimpleOutlierFinder {
if (not state.hasCalibrated() or not state.hasPredicted()) {
return false;
}
//const size_t measdim = state.calibratedSize();
// const size_t measdim = state.calibratedSize();

// auto geoId = state.referenceSurface().geometryId();

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,6 @@ ProcessCode TruthSeedSelector::execute(const AlgorithmContext& ctx) const {
const auto& hits = makeRange(particleHitsMap.equal_range(p.particleId()));
// number of recorded hits
size_t nHits = hits.size();
std::cout << "nHits = " << nHits << std::endl;
return within(rho, 0., m_cfg.rhoMax) and
within(p.position().z(), m_cfg.zMin, m_cfg.zMax) and
within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) and
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ class RootSTCFMeasurementReader_BKG : public IReader {
/// Whether the events are ordered or not
bool orderedEvents = true;
/// Whether ignore noise hits when reading
bool ignoreNoiseHits = false;
bool ignoreNoiseHits = false;
};

/// Constructor
Expand Down
8 changes: 4 additions & 4 deletions Examples/Io/Root/src/RootSTCFMeasurementReader.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -258,9 +258,9 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader::read(
// Reading ITK hits
for (size_t i = 0; i < ITKpositionX->GetSize(); ++i) {
int parentID = (*ITKparentID)[i];
if (parentID != 0){
if (parentID != 0) {
continue;
}
}
nITKHits++;

Acts::Vector3 pos((*ITKpositionX)[i], (*ITKpositionY)[i],
Expand Down Expand Up @@ -331,9 +331,9 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader::read(
// Reading MDC hits
for (size_t i = 0; i < MDCcellID->GetSize(); ++i) {
int parentID = (*MDCparentID)[i];
if (parentID != 0){
if (parentID != 0) {
continue;
}
}
nMDCHits++;

Acts::Vector3 pos((*MDCpositionX)[i], (*MDCpositionY)[i],
Expand Down
67 changes: 32 additions & 35 deletions Examples/Io/Root/src/RootSTCFMeasurementReader_BKG.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
std::map<int, int> particleITKAllHitIdx;
std::map<int, int> particleITKSigHitIdx;
std::map<int, int> particleHitIdx;

/////////////////////////////////////////////////////////////////////////////////////////////////
// I. read sim hits
/////////////////////////////////////////////////////////////////////////////////////////////////
Expand All @@ -269,20 +269,23 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
/////////////////////////////////////////////////////////////////////////////////////////////////
// 1) Reading ITK hits
/////////////////////////////////////////////////////////////////////////////////////////////////
//for (size_t i = 0; i < ITKpositionX->GetSize(); ++i) {
// for (size_t i = 0; i < ITKpositionX->GetSize(); ++i) {
for (size_t i = 0; i < ITKtype->GetSize(); ++i) {
int parentID = (*ITKparentID)[i];
if (parentID != 0){
if (parentID != 0) {
continue;
}
}
int type = (*ITKtype)[i];
// std::cout<<"ITK type is :"<<type<<std::endl;
// std::cout<<"ITK type is :"<<type<<std::endl;
int particleId = (*ITKparticleId)[i];
//bool isNoiseHit = ((*ITKmomentumX)[i] ==0 and (*ITKmomentumY)[i]==0 and (*ITKmomentumZ)[i] == 0);
bool isNoiseHit = (type == -1);
if(m_cfg.ignoreNoiseHits and isNoiseHit) {
continue;
}
bool isNoiseHit = (type == -1);
// Reset the noise particle index to -1
if (isNoiseHit) {
particleId = -1;
}
if (m_cfg.ignoreNoiseHits and isNoiseHit) {
continue;
}
nITKHits++;
Acts::Vector3 pos((*ITKpositionX)[i], (*ITKpositionY)[i],
(*ITKpositionZ)[i]);
Expand Down Expand Up @@ -342,7 +345,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(

// The last parameter could be incorrect!!!
ActsFatras::Hit hit(geoId, particleId, pos4, mom4, mom4 + delta4,
particleHitIdx[particleId]);
particleHitIdx[particleId], isNoiseHit);
unordered_hits.push_back(std::move(hit));

particleHitIdx[particleId]++;
Expand All @@ -358,17 +361,20 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
/////////////////////////////////////////////////////////////////////////////////////////////////
for (size_t i = 0; i < MDCcellID->GetSize(); ++i) {
int parentID = (*MDCparentID)[i];
if (parentID != 0){
if (parentID != 0) {
continue;
}
}
int type = (*MDCtype)[i];
// std::cout<<"MDC type is :"<<type<<std::endl;
int particleId = (*MDCparticleId)[i];
//bool isNoiseHit = ((*MDCmomentumX)[i]==0 and (*MDCmomentumY)[i]==0 and (*MDCmomentumZ)[i]==0);
bool isNoiseHit = (type == -1);
if(m_cfg.ignoreNoiseHits and isNoiseHit) {
continue;
}
bool isNoiseHit = (type == -1);
// Reset the noise particle index to -1
if (isNoiseHit) {
particleId = -1;
}
if (m_cfg.ignoreNoiseHits and isNoiseHit) {
continue;
}
nMDCHits++;

Acts::Vector3 pos((*MDCpositionX)[i], (*MDCpositionY)[i],
Expand Down Expand Up @@ -405,7 +411,6 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
// localPos.y() " << lPosition[1] << std::endl;
}


ActsFatras::Hit::Vector4 pos4{
posUpdated.x() * Acts::UnitConstants::mm,
posUpdated.y() * Acts::UnitConstants::mm,
Expand All @@ -431,7 +436,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(

// The last parameter is not simHitIdx?
ActsFatras::Hit hit(geoId, particleId, pos4, mom4, mom4 + delta4,
particleHitIdx[particleId]);
particleHitIdx[particleId], isNoiseHit);
unordered_hits.push_back(std::move(hit));
particleHitIdx[particleId]++;
}
Expand Down Expand Up @@ -472,14 +477,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
auto pos = simHit.position();
auto dir = simHit.unitDirection();
auto particleId = simHit.particleId();

// This is used to decide whether this hit is a noise (noise hit have
// zero momentum)
auto momBefore = simHit.momentum4Before();
bool isNoise = false;
if (momBefore[0] == 0 and momBefore[1] == 0 and momBefore[2] == 0) {
isNoise = true;
}
bool isNoiseHit = simHit.isNoise();

Index measurementIdx = measurements.size();
sourceLinkStorage.emplace_back(moduleGeoId, measurementIdx);
Expand All @@ -495,7 +493,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
cov(1, 1) = 0.4 * 0.4;

// Don't smear noise
if (not isNoise) {
if (not isNoiseHit) {
Acts::ActsVector<2> par{m_ITKRadius[moduleGeoId.layer() / 2 - 1] *
Acts::VectorHelpers::phi(pos) +
0.1 * stdNormal(rng),
Expand All @@ -517,7 +515,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
std::array<Acts::BoundIndices, 1> indices = {Acts::eBoundLoc0};

// Gives the drift distance a sign if it's not noise
if (not isNoise) {
if (not isNoiseHit) {
auto lpResult =
surfacePtr->globalToLocal(context.geoContext, pos, dir);
if (not lpResult.ok()) {
Expand All @@ -541,7 +539,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(

Acts::ActsSymMatrix<1> cov = Acts::ActsSymMatrix<1>::Identity();
cov(0, 0) = sigma * sigma;
if (not isNoise) {
if (not isNoiseHit) {
Acts::ActsVector<1> par{driftDistance + sigma * stdNormal(rng)};
measurements.emplace_back(
Acts::Measurement<Acts::BoundIndices, 1>(sourceLink, indices,
Expand Down Expand Up @@ -569,9 +567,8 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
}
std::cout << "nITKHits = " << nITKHits << ", nMDCHits = " << nMDCHits
<< std::endl;
} // Finish reading sim hits and transform them into measurements


} // Finish reading sim hits and transform them into measurements

/////////////////////////////////////////////////////////////////////////////////////////////////
// II. now read particles
/////////////////////////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -648,7 +645,7 @@ ActsExamples::ProcessCode ActsExamples::RootSTCFMeasurementReader_BKG::read(
}
std::cout << "nParticles = " << nParticles << std::endl;
}

// Write the collections to the EventStore
context.eventStore.add(m_cfg.outputSourceLinks, std::move(sourceLinks));
context.eventStore.add(m_cfg.outputSourceLinks + "__storage",
Expand Down
3 changes: 1 addition & 2 deletions Examples/Run/Reconstruction/STCFRecCKFTracks_BKG.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -54,8 +54,7 @@ void addRecCKFOptions(ActsExamples::Options::Description& desc) {
using boost::program_options::value;

auto opt = desc.add_options();
opt("reco-ignore-noise", bool_switch(),
"Don't read noise hits");
opt("reco-ignore-noise", bool_switch(), "Don't read noise hits");
opt("ckf-truth-smeared-seeds", bool_switch(),
"Use track parameters smeared from truth particles for steering CKF");
opt("ckf-truth-estimated-seeds", bool_switch(),
Expand Down
8 changes: 7 additions & 1 deletion Fatras/include/ActsFatras/EventData/Hit.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,11 @@ class Hit {
/// position on the given surface.
Hit(Acts::GeometryIdentifier geometryId, Barcode particleId,
const Vector4& pos4, const Vector4& before4, const Vector4& after4,
int32_t index_ = -1)
int32_t index_ = -1, bool isNoise_ = false)
: m_geometryId(geometryId),
m_particleId(particleId),
m_index(index_),
m_isNoise(isNoise_),
m_pos4(pos4),
m_before4(before4),
m_after4(after4) {}
Expand All @@ -67,6 +68,9 @@ class Hit {
/// @retval negative if the hit index is undefined.
constexpr int32_t index() const { return m_index; }

/// Added by xiaocong. Whether this is a noise
bool isNoise() const { return m_isNoise; }

/// Space-time position four-vector.
const Vector4& fourPosition() const { return m_pos4; }
/// Three-position, i.e. spatial coordinates without the time.
Expand Down Expand Up @@ -113,6 +117,8 @@ class Hit {
Vector4 m_before4 = Vector4::Zero();
/// Global particle energy-momentum four-vector after the hit.
Vector4 m_after4 = Vector4::Zero();
/// Whether the hit is a noise
bool m_isNoise = false;
};

} // namespace ActsFatras
10 changes: 6 additions & 4 deletions Fatras/include/ActsFatras/EventData/Particle.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,9 @@ class Particle {
Particle &setCharge(Scalar charge) { return m_charge = charge, *this; }

///**************set the particle ITKHits
Particle &setITKHits(int ITKHits) { return m_ITKHits = ITKHits, *this; }///***
Particle &setITKHits(int ITKHits) {
return m_ITKHits = ITKHits, *this;
} ///***

/// Change the energy by the given amount.
///
Expand All @@ -131,8 +133,8 @@ class Particle {
}
return *this;
}
///***********get particle ITKHits
constexpr int particleITKHits() const { return m_ITKHits; }/////***
///***********get particle ITKHits
constexpr int particleITKHits() const { return m_ITKHits; } /////***
/////**************
/// Particle identifier within an event.
constexpr Barcode particleId() const { return m_particleId; }
Expand Down Expand Up @@ -224,7 +226,7 @@ class Particle {
Scalar m_pathInX0 = Scalar(0);
Scalar m_pathInL0 = Scalar(0);
///*******particle ITKHits
int m_ITKHits;///************
int m_ITKHits; ///************
};

std::ostream &operator<<(std::ostream &os, const Particle &particle);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ class CombineSelectors {
///
/// @tparam Ts the types of the selected inputs
template <typename... Ts>
bool operator()(const Ts &... things) const {
bool operator()(const Ts &...things) const {
static_assert(
(true && ... && std::is_same_v<bool, decltype(Selectors()(things...))>),
"Not all selectors conform to the expected interface (bool)(const "
Expand All @@ -52,7 +52,7 @@ class CombineSelectors {
std::tuple<Selectors...> m_selectors;

template <std::size_t... Is, typename... Ts>
bool impl(std::index_sequence<Is...>, const Ts &... things) const {
bool impl(std::index_sequence<Is...>, const Ts &...things) const {
Combine combine;
// compute status for all selectors
bool status[] = {std::get<Is>(m_selectors)(things...)...};
Expand Down

0 comments on commit 3b25620

Please sign in to comment.