File indexing completed on 2024-09-27 07:02:23
0001 #include <cmath>
0002 #include <algorithm>
0003 #include <unordered_map>
0004
0005 #include "Acts/ActsVersion.hpp"
0006 #include "Acts/Definitions/Units.hpp"
0007 #include "Acts/Definitions/Common.hpp"
0008 #include "Acts/Geometry/TrackingGeometry.hpp"
0009
0010 #include "Acts/Seeding/BinFinder.hpp"
0011 #include "Acts/Seeding/BinnedSPGroup.hpp"
0012 #include "Acts/Seeding/SeedfinderConfig.hpp"
0013 #include "Acts/Seeding/SpacePointGrid.hpp"
0014 #include "Acts/Seeding/SeedFilter.hpp"
0015 #include "Acts/Seeding/Seedfinder.hpp"
0016 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0017 #include "Acts/Surfaces/PerigeeSurface.hpp"
0018
0019
0020 #include "GaudiAlg/GaudiAlgorithm.h"
0021 #include "GaudiKernel/ToolHandle.h"
0022 #include "GaudiAlg/Transformer.h"
0023 #include "GaudiAlg/GaudiTool.h"
0024 #include "GaudiKernel/RndmGenerators.h"
0025 #include "Gaudi/Property.h"
0026
0027 #include "JugBase/DataHandle.h"
0028 #include "JugBase/IGeoSvc.h"
0029 #include "JugBase/BField/DD4hepBField.h"
0030 #include "JugTrack/IndexSourceLink.hpp"
0031 #include "JugTrack/Measurement.hpp"
0032 #include "JugTrack/Track.hpp"
0033
0034 #include "edm4eic/TrackerHitCollection.h"
0035
0036 #include "Math/Vector3D.h"
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 namespace Jug::Reco {
0052
0053
0054
0055
0056
0057
0058 class TrackParamACTSSeeding : public GaudiAlgorithm {
0059 public:
0060 DataHandle<IndexSourceLinkContainer>
0061 m_inputSourceLinks { "inputSourceLinks",
0062 Gaudi::DataHandle::Reader, this };
0063 DataHandle<MeasurementContainer>
0064 m_inputMeasurements { "inputMeasurements",
0065 Gaudi::DataHandle::Reader, this};
0066 DataHandle<edm4eic::TrackerHitCollection>
0067 m_inputHitCollection { "inputHitCollection",
0068 Gaudi::DataHandle::Reader, this };
0069 DataHandle<TrackParametersContainer>
0070 m_outputInitialTrackParameters {
0071 "outputInitialTrackParameters",
0072 Gaudi::DataHandle::Writer, this };
0073
0074 SmartIF<IGeoSvc> m_geoSvc;
0075 Acts::GeometryContext m_geoContext;
0076 std::shared_ptr<const Jug::BField::DD4hepBField> m_BField =
0077 nullptr;
0078 Acts::MagneticFieldContext m_fieldContext;
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089 class SpacePoint : edm4eic::TrackerHitData {
0090 public:
0091 int32_t _measurementIndex;
0092
0093
0094 SpacePoint(const edm4eic::TrackerHit h,
0095 const int32_t measurementIndex)
0096 : _measurementIndex(measurementIndex)
0097 {
0098 position = h.getPosition();
0099 positionError = h.getPositionError();
0100 }
0101 constexpr float x() const { return position.x; }
0102 constexpr float y() const { return position.y; }
0103 constexpr float z() const { return position.z; }
0104 float r() const { return std::hypot(x(), y()); }
0105 float varianceR() const
0106 {
0107 return (std::pow(x(), 2) * positionError.xx +
0108 std::pow(y(), 2) * positionError.yy) /
0109 (std::pow(x(), 2) + std::pow(y(), 2));
0110 }
0111 constexpr float varianceZ() const { return positionError.zz; }
0112 constexpr uint32_t measurementIndex() const {
0113 return _measurementIndex; }
0114 };
0115
0116 static bool spCompare(SpacePoint r, SpacePoint s)
0117 {
0118 return
0119 std::hypot(r.x(), r.y(), r.z()) <
0120 std::hypot(s.x(), s.y(), s.z());
0121 }
0122
0123
0124 using SeedContainer = std::vector<Acts::Seed<SpacePoint>>;
0125
0126
0127
0128 using ProtoTrack = std::vector<Index>;
0129
0130
0131 using ProtoTrackContainer = std::vector<ProtoTrack>;
0132
0133 using SpacePointContainer = std::vector<SpacePoint>;
0134
0135 struct Config {
0136
0137
0138
0139
0140
0141
0142
0143 std::vector<std::string> inputSpacePoints;
0144
0145 std::string outputSeeds;
0146
0147 std::string outputProtoTracks;
0148
0149 float bFieldInZ = 3 * Acts::UnitConstants::T;
0150 float minPt = 500 * Acts::UnitConstants::MeV;
0151 float rMax = 440 * Acts::UnitConstants::mm;
0152 float zMin = -1500 * Acts::UnitConstants::mm;
0153 float zMax = 1700 * Acts::UnitConstants::mm;
0154 float deltaRMin = 1 * Acts::UnitConstants::mm;
0155 float deltaRMax = 600 * Acts::UnitConstants::mm;
0156 float cotThetaMax = sinh(4.01);
0157
0158 float collisionRegionMin = -250 * Acts::UnitConstants::mm;
0159 float collisionRegionMax = 250 * Acts::UnitConstants::mm;
0160 float maxSeedsPerSpM = 6;
0161 float sigmaScattering = 5;
0162 float radLengthPerSeed = 0.1;
0163 float beamPosX = 0 * Acts::UnitConstants::mm;
0164 float beamPosY = 0 * Acts::UnitConstants::mm;
0165 float impactMax = 3 * Acts::UnitConstants::mm;
0166
0167
0168
0169 double bFieldMin = 0.1 * Acts::UnitConstants::T;
0170
0171
0172 double sigmaLoc0 = 25 * Acts::UnitConstants::um;
0173
0174 double sigmaLoc1 = 100 * Acts::UnitConstants::um;
0175
0176 double sigmaPhi = 0.02 * Acts::UnitConstants::degree;
0177
0178 double sigmaTheta = 0.02 * Acts::UnitConstants::degree;
0179
0180 double sigmaQOverP = 0.1 / Acts::UnitConstants::GeV;
0181
0182 double sigmaT0 = 1400 * Acts::UnitConstants::s;
0183
0184 int numPhiNeighbors = 1;
0185
0186
0187 std::vector<std::pair<int, int> > zBinNeighborsTop;
0188 std::vector<std::pair<int, int> > zBinNeighborsBottom;
0189 } m_cfg;
0190 Acts::SpacePointGridConfig m_gridCfg;
0191 Acts::SeedfinderConfig<SpacePoint> m_finderCfg;
0192
0193
0194 Acts::BoundSymMatrix m_covariance =
0195 Acts::BoundSymMatrix::Zero();
0196
0197 public:
0198 TrackParamACTSSeeding(const std::string &name,
0199 ISvcLocator* svcLoc)
0200 : GaudiAlgorithm(name, svcLoc) {
0201 declareProperty("inputSourceLinks",
0202 m_inputSourceLinks, "");
0203 declareProperty("inputMeasurements",
0204 m_inputMeasurements, "");
0205 declareProperty("inputHitCollection",
0206 m_inputHitCollection, "");
0207 declareProperty("outputInitialTrackParameters",
0208 m_outputInitialTrackParameters, "");
0209 }
0210
0211 StatusCode initialize() override;
0212
0213 void
0214 findSeed(SeedContainer &seeds,
0215 const edm4eic::TrackerHitCollection *hits,
0216 const IndexSourceLinkContainer *sourceLinks,
0217 const MeasurementContainer *measurements,
0218 Acts::Seedfinder<SpacePoint>::State &state);
0219
0220 StatusCode execute() override;
0221 };
0222
0223
0224 StatusCode TrackParamACTSSeeding::initialize()
0225 {
0226 if (GaudiAlgorithm::initialize().isFailure()) {
0227 return StatusCode::FAILURE;
0228 }
0229
0230 m_geoSvc = service("GeoSvc");
0231 if (m_geoSvc == nullptr) {
0232 error() << "Unable to locate Geometry Service. " << endmsg;
0233 return StatusCode::FAILURE;
0234 }
0235
0236 m_BField = std::dynamic_pointer_cast<
0237 const Jug::BField::DD4hepBField>(
0238 m_geoSvc->getFieldProvider());
0239 m_fieldContext = Jug::BField::BFieldVariant(m_BField);
0240
0241 m_gridCfg.bFieldInZ = m_cfg.bFieldInZ;
0242 m_gridCfg.minPt = m_cfg.minPt;
0243 m_gridCfg.rMax = m_cfg.rMax;
0244 m_gridCfg.zMax = m_cfg.zMax;
0245 m_gridCfg.zMin = m_cfg.zMin;
0246 m_gridCfg.cotThetaMax = m_cfg.cotThetaMax;
0247
0248 m_gridCfg.deltaRMax = m_cfg.deltaRMax;
0249
0250
0251 Acts::SeedFilterConfig filterCfg;
0252 filterCfg.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0253 m_finderCfg.seedFilter =
0254 std::make_unique<Acts::SeedFilter<SpacePoint>>(
0255 Acts::SeedFilter<SpacePoint>(filterCfg));
0256
0257 m_finderCfg.rMax = m_cfg.rMax;
0258 m_finderCfg.deltaRMin = m_cfg.deltaRMin;
0259 m_finderCfg.deltaRMax = m_cfg.deltaRMax;
0260 m_finderCfg.collisionRegionMin = m_cfg.collisionRegionMin;
0261 m_finderCfg.collisionRegionMax = m_cfg.collisionRegionMax;
0262 m_finderCfg.zMin = m_cfg.zMin;
0263 m_finderCfg.zMax = m_cfg.zMax;
0264 m_finderCfg.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0265 m_finderCfg.cotThetaMax = m_cfg.cotThetaMax;
0266 m_finderCfg.sigmaScattering = m_cfg.sigmaScattering;
0267 m_finderCfg.radLengthPerSeed = m_cfg.radLengthPerSeed;
0268 m_finderCfg.minPt = m_cfg.minPt;
0269 m_finderCfg.bFieldInZ = m_cfg.bFieldInZ;
0270 m_finderCfg.beamPos =
0271 Acts::Vector2(m_cfg.beamPosX, m_cfg.beamPosY);
0272 m_finderCfg.impactMax = m_cfg.impactMax;
0273
0274
0275
0276 m_covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) =
0277 std::pow(m_cfg.sigmaLoc0, 2);
0278 m_covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) =
0279 std::pow(m_cfg.sigmaLoc1, 2);
0280 m_covariance(Acts::eBoundPhi, Acts::eBoundPhi) =
0281 std::pow(m_cfg.sigmaPhi, 2);
0282 m_covariance(Acts::eBoundTheta, Acts::eBoundTheta) =
0283 std::pow(m_cfg.sigmaTheta, 2);
0284 m_covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) =
0285 std::pow(m_cfg.sigmaQOverP, 2);
0286 m_covariance(Acts::eBoundTime, Acts::eBoundTime) =
0287 std::pow(m_cfg.sigmaT0, 2);
0288
0289 return StatusCode::SUCCESS;
0290 }
0291
0292 void TrackParamACTSSeeding::
0293 findSeed(SeedContainer &seeds,
0294 const edm4eic::TrackerHitCollection *hits,
0295 const IndexSourceLinkContainer *sourceLinks,
0296 const MeasurementContainer *measurements,
0297 Acts::Seedfinder<SpacePoint>::State &state)
0298 {
0299
0300
0301 std::vector<SpacePoint> spacePoint;
0302 std::vector<const SpacePoint *> spacePointPtrs;
0303
0304 Acts::Extent rRangeSPExtent;
0305
0306 std::shared_ptr<const Acts::TrackingGeometry>
0307 trackingGeometry = m_geoSvc->trackingGeometry();
0308
0309 if (msgLevel(MSG::DEBUG)) {
0310 debug() << __FILE__ << ':' << __LINE__ << ": "
0311 << sourceLinks->size() << ' '
0312 << measurements->size() << ' '
0313 << hits->size() << endmsg;
0314 }
0315 auto its = sourceLinks->begin();
0316 auto itm = measurements->begin();
0317 for (; its != sourceLinks->end() &&
0318 itm != measurements->end();
0319 its++, itm++) {
0320 const Acts::Surface *surface = trackingGeometry->findSurface(its->get().geometryId());
0321 if (surface != nullptr) {
0322 Acts::Vector3 v = surface->localToGlobal(m_geoContext, {std::get<Acts::Measurement<Acts::BoundIndices, 2>>(*itm).parameters()[0], std::get<Acts::Measurement<Acts::BoundIndices, 2>>(*itm).parameters()[1]}, {0, 0, 0});
0323 if (msgLevel(MSG::DEBUG)) {
0324 debug() << __FILE__ << ':' << __LINE__ << ": "
0325 << its - sourceLinks->begin() << ' '
0326
0327 << v[0] << ' ' << v[1] << ' ' << v[2]
0328 << endmsg;
0329 }
0330 #ifdef USE_LOCAL_COORD
0331 spacePoint.push_back(
0332 SpacePoint(
0333 edm4eic::TrackerHit(
0334 static_cast<uint64_t>(spacePoint.size()),
0335 {v[0], v[1], v[2]},
0336 {25.0e-6 / 3.0,
0337 25.0e-6 / 3.0, 0.0},
0338 0.0, 10.0, 0.05, 0.0),
0339 static_cast<int32_t>(spacePoint.size())));
0340 spacePointPtrs.push_back(&spacePoint.back());
0341 rRangeSPExtent.check({ spacePoint.back().x(),
0342 spacePoint.back().y(),
0343 spacePoint.back().z() });
0344 #endif
0345 }
0346 }
0347
0348 for(const auto &h : *hits) {
0349 if (msgLevel(MSG::DEBUG)) {
0350 debug() << __FILE__ << ':' << __LINE__ << ": "
0351 << ' ' << h.getPosition().x
0352 << ' ' << h.getPosition().y
0353 << ' ' << h.getPosition().z
0354 << ' ' << h.getPositionError().xx
0355 << ' ' << h.getPositionError().yy
0356 << ' ' << h.getPositionError().zz
0357 << ' ' << h.getTime()
0358 << ' ' << h.getTimeError()
0359 << ' ' << h.getEdep()
0360 << ' ' << h.getEdepError()
0361 << endmsg;
0362 }
0363 #ifndef USE_LOCAL_COORD
0364 spacePoint.push_back(SpacePoint(h, static_cast<int32_t>(spacePoint.size())));
0365 spacePointPtrs.push_back(&spacePoint.back());
0366 rRangeSPExtent.check({ spacePoint.back().x(),
0367 spacePoint.back().y(),
0368 spacePoint.back().z() });
0369 #endif
0370 }
0371 if (msgLevel(MSG::DEBUG)) {
0372 debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0373 }
0374
0375 auto extractGlobalQuantities =
0376 [=](const SpacePoint& sp, float, float, float) ->
0377 std::pair<Acts::Vector3, Acts::Vector2> {
0378 Acts::Vector3 position { sp.x(), sp.y(), sp.z() };
0379 Acts::Vector2 covariance {
0380 sp.varianceR(), sp.varianceZ() };
0381 return std::make_pair(position, covariance);
0382 };
0383 if (msgLevel(MSG::DEBUG)) {
0384 debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0385 }
0386
0387 auto bottomBinFinder =
0388 std::make_shared<Acts::BinFinder<SpacePoint>>(
0389 Acts::BinFinder<SpacePoint>(m_cfg.zBinNeighborsBottom,
0390 m_cfg.numPhiNeighbors));
0391 if (msgLevel(MSG::DEBUG)) {
0392 debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0393 }
0394 auto topBinFinder =
0395 std::make_shared<Acts::BinFinder<SpacePoint>>(
0396 Acts::BinFinder<SpacePoint>(m_cfg.zBinNeighborsTop,
0397 m_cfg.numPhiNeighbors));
0398 if (msgLevel(MSG::DEBUG)) {
0399 debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0400 }
0401
0402 auto grid =
0403 Acts::SpacePointGridCreator::createGrid<SpacePoint>(
0404 m_gridCfg);
0405 if (msgLevel(MSG::DEBUG)) {
0406 debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0407 }
0408
0409 auto spacePointsGrouping =
0410 Acts::BinnedSPGroup<SpacePoint>(
0411 spacePointPtrs.begin(), spacePointPtrs.end(),
0412 extractGlobalQuantities, bottomBinFinder,
0413 topBinFinder, std::move(grid), m_finderCfg);
0414 auto finder = Acts::Seedfinder<SpacePoint>(m_finderCfg);
0415
0416 if (msgLevel(MSG::DEBUG)) {
0417 debug() << __FILE__ << ':' << __LINE__
0418 << ": spacePointsGrouping.size() = "
0419 << spacePointsGrouping.size() << endmsg;
0420 }
0421 #if 0
0422 topBinFinder.get();
0423 bottomBinFinder.get();
0424 #endif
0425
0426 seeds.clear();
0427
0428 auto group = spacePointsGrouping.begin();
0429 auto groupEnd = spacePointsGrouping.end();
0430 #if 1
0431 for (; !(group == groupEnd); ++group) {
0432 finder.createSeedsForGroup(
0433 state, std::back_inserter(seeds),
0434 group.bottom(), group.middle(), group.top(),
0435 rRangeSPExtent);
0436 }
0437 #endif
0438
0439 if (msgLevel(MSG::DEBUG)) {
0440 debug() << "seeds.size() = " << seeds.size() << endmsg;
0441 }
0442 }
0443
0444 StatusCode TrackParamACTSSeeding::execute()
0445 {
0446 const edm4eic::TrackerHitCollection *hits =
0447 m_inputHitCollection.get();
0448 const IndexSourceLinkContainer *sourceLinks =
0449 m_inputSourceLinks.get();
0450 const MeasurementContainer *measurements =
0451 m_inputMeasurements.get();
0452
0453 auto initTrackParameters =
0454 m_outputInitialTrackParameters.createAndPut();
0455
0456 static thread_local SeedContainer seeds;
0457 static thread_local Acts::Seedfinder<SpacePoint>::State state;
0458
0459 findSeed(seeds, hits, sourceLinks, measurements, state);
0460
0461 TrackParametersContainer trackParameters;
0462 ProtoTrackContainer tracks;
0463 trackParameters.reserve(seeds.size());
0464 tracks.reserve(seeds.size());
0465
0466 std::shared_ptr<const Acts::TrackingGeometry>
0467 trackingGeometry = m_geoSvc->trackingGeometry();
0468 std::shared_ptr<const Acts::MagneticFieldProvider>
0469 magneticField = m_geoSvc->getFieldProvider();
0470
0471 if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0472 auto bCache = magneticField->makeCache(m_fieldContext);
0473
0474 std::unordered_map<size_t, bool> spTaken;
0475
0476 for (size_t iseed = 0; iseed < seeds.size(); iseed++) {
0477 const auto &seed = seeds[iseed];
0478
0479 const auto bottomSP = seed.sp().front();
0480 auto hitIdx = bottomSP->measurementIndex();
0481
0482
0483
0484
0485
0486
0487
0488 hitIdx = std::min(hitIdx, static_cast<uint32_t>(
0489 sourceLinks->size() - 1));
0490 const Acts::Surface *surface = nullptr;
0491 for (auto &s : *sourceLinks) {
0492 surface = trackingGeometry->findSurface(s.get().geometryId());
0493 if (surface != nullptr &&
0494 surface->isOnSurface(
0495 m_geoContext,
0496 {bottomSP->x(), bottomSP->y(), bottomSP->z()},
0497 {0, 0, 0})) {
0498 break;
0499 }
0500 }
0501 if (surface == nullptr && msgLevel(MSG::DEBUG)) {
0502 debug() << "hit " << hitIdx
0503 << " is not found in the tracking gemetry"
0504 << endmsg;
0505 continue;
0506 }
0507
0508 if (msgLevel(MSG::DEBUG)) {
0509 debug() << __FILE__ << ':' << __LINE__
0510 << ": iseed = " << iseed << ", "
0511 << surface->type() << ", "
0512 << surface->center(m_geoContext).x() << ", "
0513 << surface->center(m_geoContext).y() << ", "
0514 << surface->center(m_geoContext).z()
0515 << endmsg;
0516 }
0517
0518
0519 auto fieldRes = magneticField->getField(
0520 {bottomSP->x(), bottomSP->y(), bottomSP->z()},
0521 bCache);
0522 if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0523
0524 auto optParams = Acts::estimateTrackParamsFromSeed(
0525 m_geoContext, seed.sp().begin(), seed.sp().end(),
0526 *surface, *fieldRes, m_cfg.bFieldMin);
0527 if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0528 if (not optParams.has_value()) {
0529 debug() << "Estimation of track parameters for seed "
0530 << iseed << " failed." << endmsg;
0531 continue;
0532 }
0533 else if (!(spTaken[seed.sp()[0]->measurementIndex()] ||
0534 spTaken[seed.sp()[1]->measurementIndex()] ||
0535 spTaken[seed.sp()[2]->measurementIndex()])) {
0536 const auto& params = optParams.value();
0537 const double charge =
0538 std::copysign(1, params[Acts::eBoundQOverP]);
0539 initTrackParameters->emplace_back(
0540 surface->getSharedPtr(), params, charge,
0541 m_covariance);
0542
0543 #if 0
0544 spTaken[seed.sp()[0]->measurementIndex()] = true;
0545 spTaken[seed.sp()[1]->measurementIndex()] = true;
0546 spTaken[seed.sp()[2]->measurementIndex()] = true;
0547 #endif
0548 }
0549 if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0550 }
0551
0552 return StatusCode::SUCCESS;
0553 }
0554
0555 DECLARE_COMPONENT(TrackParamACTSSeeding)
0556 }