From b7d3b2e561441b5c1c641b3f2ccfe25ce7b38321 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 10 Jul 2024 11:52:29 +0200 Subject: [PATCH 1/5] attempt at adding extrapolation to calo surface --- ACTSTracking/ACTSSeededCKFTrackingProc.hxx | 4 ++ src/ACTSSeededCKFTrackingProc.cxx | 55 ++++++++++++++++++++++ 2 files changed, 59 insertions(+) diff --git a/ACTSTracking/ACTSSeededCKFTrackingProc.hxx b/ACTSTracking/ACTSSeededCKFTrackingProc.hxx index e123f3e..642eaf9 100644 --- a/ACTSTracking/ACTSSeededCKFTrackingProc.hxx +++ b/ACTSTracking/ACTSSeededCKFTrackingProc.hxx @@ -60,6 +60,10 @@ class ACTSSeededCKFTrackingProc : public ACTSProcBase { bool _runCKF = true; bool _propagateBackward = false; + // Extrapolation to calo settings + float _caloFaceR = 1857; //mm + float _caloFaceZ = 2307; //mm + // Seed finding configuration float _seedFinding_rMax = 150; float _seedFinding_deltaRMin = 5; diff --git a/src/ACTSSeededCKFTrackingProc.cxx b/src/ACTSSeededCKFTrackingProc.cxx index 788ca75..291153a 100644 --- a/src/ACTSSeededCKFTrackingProc.cxx +++ b/src/ACTSSeededCKFTrackingProc.cxx @@ -27,6 +27,10 @@ #include #include +#include +#include +#include + using namespace Acts::UnitLiterals; #include "Helpers.hxx" @@ -72,6 +76,15 @@ ACTSSeededCKFTrackingProc::ACTSSeededCKFTrackingProc() "Track error estimate, local position (mm).", _initialTrackError_pos, 10_um); + // Extrapolation to calo surface + registerProcessorParameter("CaloFace_Radius", + "ECAL Inner Radius (mm).", + _caloFaceR, _caloFaceR); + + registerProcessorParameter("CaloFace_Z", + "ECAL half length (mm).", + _caloFaceZ, _caloFaceZ); + // Seeding configurations registerProcessorParameter( "SeedingLayers", @@ -680,6 +693,48 @@ void ACTSSeededCKFTrackingProc::processEvent(LCEvent *evt) { streamlog_out(DEBUG) << "\tnStates " << trackTip.nTrackStates() << std::endl; + //FM: attempt at propagating to calo surface + double radius = _caloFaceR * Acts::UnitConstants::mm; // Radius in mm + double halfLength = _caloFaceZ * Acts::UnitConstants::mm; // Half-height in mm + + // Create the CaloSurface + auto caloCylinder = std::make_shared(radius, halfLength); + auto caloSurface = Acts::Surface::makeShared(Acts::Transform3::Identity(), caloCylinder); + + streamlog_out(DEBUG) << "Created the calo cylindrical surface" << std::endl; + + // define start parameters - swap this out with some smart call + Acts::BoundVector params = trackTip.parameters(); + double d0 = params[Acts::eBoundLoc0]; + double z0 = params[Acts::eBoundLoc1]; + double phi = params[Acts::eBoundPhi]; + double theta = params[Acts::eBoundTheta]; + double qoverp = params[Acts::eBoundQOverP]; + double time = params[Acts::eBoundTime]; + + Acts::Vector3 pos(d0 * cos(phi), d0 * sin(phi), z0); + Acts::BoundMatrix cov = trackTip.covariance(); + + Acts::CurvilinearTrackParameters start(Acts::VectorHelpers::makeVector4(pos, time), phi, theta, qoverp, cov, Acts::ParticleHypothesis::pion()); + streamlog_out(DEBUG) << "Retrieved the start track parameters" << std::endl; + + // Set propagator options + Acts::PropagatorOptions myCaloPropOptions(geometryContext(), magneticFieldContext()); + myCaloPropOptions.pathLimit = 20 * Acts::UnitConstants::m; + std::cout << "init prop option" << std::endl; + + auto resultProp = propagator.propagate(start, *caloSurface, myCaloPropOptions); + if (resultProp.ok()) { + std::cout << "Good!" << std::endl; + //FM: should read back the result of the propagation? + std::cout << "CaloPhi:" << phi << " " << resultProp.value().endParameters->parameters()[Acts::eBoundPhi] << std::endl; + // FM the line below is not correct, but it's what I'd like to do + //trackTip.addTrackState(result); + } + else { + std::cout << "Propagation failed" << std::endl; + } + // Make track object EVENT::Track *track = ACTSTracking::ACTS2Marlin_track( trackTip, magneticField(), magCache); From 8a1bcc81d89a52f9df08c1dfc1967ecb13a0e0c6 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 10 Jul 2024 17:46:42 +0200 Subject: [PATCH 2/5] moved to Helper --- ACTSTracking/Helpers.hxx | 5 ++ src/ACTSSeededCKFTrackingProc.cxx | 48 +--------- src/Helpers.cxx | 141 ++++++++++++++++++++++++++++++ 3 files changed, 147 insertions(+), 47 deletions(-) diff --git a/ACTSTracking/Helpers.hxx b/ACTSTracking/Helpers.hxx index 3cdd4cd..c84943a 100644 --- a/ACTSTracking/Helpers.hxx +++ b/ACTSTracking/Helpers.hxx @@ -83,6 +83,11 @@ EVENT::Track* ACTS2Marlin_track( std::shared_ptr magneticField, Acts::MagneticFieldProvider::Cache& magCache); +EVENT::Track* ACTS2Marlin_track( + const TrackResult& fitter_res, + std::shared_ptr magneticField, + Acts::MagneticFieldProvider::Cache& magCache, double caloFaceR, double caloFaceZ, Acts::GeometryContext geoContext, Acts::MagneticFieldContext magFieldContext, std::shared_ptr trackingGeo); + //! Convert ACTS track state class to Marlin class /** * \param location Location where the track state is defined (ie: `AtIP`) diff --git a/src/ACTSSeededCKFTrackingProc.cxx b/src/ACTSSeededCKFTrackingProc.cxx index 291153a..376c51d 100644 --- a/src/ACTSSeededCKFTrackingProc.cxx +++ b/src/ACTSSeededCKFTrackingProc.cxx @@ -27,10 +27,6 @@ #include #include -#include -#include -#include - using namespace Acts::UnitLiterals; #include "Helpers.hxx" @@ -693,51 +689,9 @@ void ACTSSeededCKFTrackingProc::processEvent(LCEvent *evt) { streamlog_out(DEBUG) << "\tnStates " << trackTip.nTrackStates() << std::endl; - //FM: attempt at propagating to calo surface - double radius = _caloFaceR * Acts::UnitConstants::mm; // Radius in mm - double halfLength = _caloFaceZ * Acts::UnitConstants::mm; // Half-height in mm - - // Create the CaloSurface - auto caloCylinder = std::make_shared(radius, halfLength); - auto caloSurface = Acts::Surface::makeShared(Acts::Transform3::Identity(), caloCylinder); - - streamlog_out(DEBUG) << "Created the calo cylindrical surface" << std::endl; - - // define start parameters - swap this out with some smart call - Acts::BoundVector params = trackTip.parameters(); - double d0 = params[Acts::eBoundLoc0]; - double z0 = params[Acts::eBoundLoc1]; - double phi = params[Acts::eBoundPhi]; - double theta = params[Acts::eBoundTheta]; - double qoverp = params[Acts::eBoundQOverP]; - double time = params[Acts::eBoundTime]; - - Acts::Vector3 pos(d0 * cos(phi), d0 * sin(phi), z0); - Acts::BoundMatrix cov = trackTip.covariance(); - - Acts::CurvilinearTrackParameters start(Acts::VectorHelpers::makeVector4(pos, time), phi, theta, qoverp, cov, Acts::ParticleHypothesis::pion()); - streamlog_out(DEBUG) << "Retrieved the start track parameters" << std::endl; - - // Set propagator options - Acts::PropagatorOptions myCaloPropOptions(geometryContext(), magneticFieldContext()); - myCaloPropOptions.pathLimit = 20 * Acts::UnitConstants::m; - std::cout << "init prop option" << std::endl; - - auto resultProp = propagator.propagate(start, *caloSurface, myCaloPropOptions); - if (resultProp.ok()) { - std::cout << "Good!" << std::endl; - //FM: should read back the result of the propagation? - std::cout << "CaloPhi:" << phi << " " << resultProp.value().endParameters->parameters()[Acts::eBoundPhi] << std::endl; - // FM the line below is not correct, but it's what I'd like to do - //trackTip.addTrackState(result); - } - else { - std::cout << "Propagation failed" << std::endl; - } - // Make track object EVENT::Track *track = ACTSTracking::ACTS2Marlin_track( - trackTip, magneticField(), magCache); + trackTip, magneticField(), magCache, _caloFaceR, _caloFaceZ, geometryContext(), magneticFieldContext(), trackingGeometry()); // Save results trackCollection->addElement(track); diff --git a/src/Helpers.cxx b/src/Helpers.cxx index 2eae839..39091e9 100644 --- a/src/Helpers.cxx +++ b/src/Helpers.cxx @@ -6,10 +6,23 @@ #include #include +#include +#include +#include + #include #include "config.h" +#include +#include +#include + + +using Stepper = Acts::EigenStepper<>; +using Navigator = Acts::Navigator; +using Propagator = Acts::Propagator; + namespace ACTSTracking { std::string findFile(const std::string& inpath) { @@ -301,6 +314,134 @@ EVENT::Track* ACTS2Marlin_track( return track; } +// Conversion with extrapolation to calorimeter face +EVENT::Track* ACTS2Marlin_track( + const TrackResult& fitter_res, + std::shared_ptr magneticField, + Acts::MagneticFieldProvider::Cache& magCache, double caloFaceR, double caloFaceZ, + Acts::GeometryContext geoContext, Acts::MagneticFieldContext magFieldContext, + std::shared_ptr trackingGeo) +{ + IMPL::TrackImpl* track = new IMPL::TrackImpl; + + track->setChi2(fitter_res.chi2()); + track->setNdf(fitter_res.nDoF()); + + const Acts::Vector3 zeroPos(0, 0, 0); + Acts::Result fieldRes = magneticField->getField(zeroPos, magCache); + if (!fieldRes.ok()) { + throw std::runtime_error("Field lookup error: " + fieldRes.error().value()); + } + Acts::Vector3 field = *fieldRes; + + const Acts::BoundVector& params = fitter_res.parameters(); + const Acts::BoundMatrix& covariance = fitter_res.covariance(); + EVENT::TrackState* trackStateAtIP = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtIP, params, covariance, field[2] / Acts::UnitConstants::T); + track->trackStates().push_back(trackStateAtIP); + + EVENT::TrackerHitVec hitsOnTrack; + EVENT::TrackStateVec statesOnTrack; + + for (const auto& trk_state : fitter_res.trackStatesReversed()) + { + if (!trk_state.hasUncalibratedSourceLink()) continue; + + auto sl = trk_state.getUncalibratedSourceLink() + .get(); + EVENT::TrackerHit* curr_hit = sl.lciohit(); + hitsOnTrack.push_back(curr_hit); + + const Acts::Vector3 hitPos(curr_hit->getPosition()[0], + curr_hit->getPosition()[1], + curr_hit->getPosition()[2]); + + Acts::Result fieldRes = + magneticField->getField(hitPos, magCache); + if (!fieldRes.ok()) { + throw std::runtime_error("Field lookup error: " + + fieldRes.error().value()); + } + Acts::Vector3 field = *fieldRes; + + EVENT::TrackState* trackState = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtOther, trk_state.smoothed(), + trk_state.smoothedCovariance(), field[2] / Acts::UnitConstants::T); + statesOnTrack.push_back(trackState); + } + + // Create the CaloSurface + auto caloCylinder = std::make_shared(caloFaceR, caloFaceZ); + auto caloSurface = Acts::Surface::makeShared(Acts::Transform3::Identity(), caloCylinder); + + // define start parameters - swap this out with some smart call + double d0 = params[Acts::eBoundLoc0]; + double z0 = params[Acts::eBoundLoc1]; + double phi = params[Acts::eBoundPhi]; + double theta = params[Acts::eBoundTheta]; + double qoverp = params[Acts::eBoundQOverP]; + double time = params[Acts::eBoundTime]; + + Acts::Vector3 pos(d0 * cos(phi), d0 * sin(phi), z0); + + Acts::CurvilinearTrackParameters start(Acts::VectorHelpers::makeVector4(pos, time), phi, theta, qoverp, covariance, Acts::ParticleHypothesis::pion()); + + // Set propagator options + Acts::PropagatorOptions myCaloPropOptions(geoContext, magFieldContext); + myCaloPropOptions.pathLimit = 20 * Acts::UnitConstants::m; + std::cout << "init prop option" << std::endl; + + //Let's try with our private propagator + // Configurations + Navigator::Config navigatorCfg{trackingGeo}; + navigatorCfg.resolvePassive = false; + navigatorCfg.resolveMaterial = true; + navigatorCfg.resolveSensitive = true; + + // construct all components for the fitter + Stepper stepper(magneticField); + Navigator navigator(navigatorCfg); + Propagator mypropagator(std::move(stepper), std::move(navigator)); + auto resultProp = mypropagator.propagate(start, *caloSurface, myCaloPropOptions); + if (resultProp.value().endParameters.has_value()) { + std::cout << "Good!" << std::endl; + //FM: should read back the result of the propagation? + std::cout << "CaloPhi:" << phi << " " << resultProp.value().endParameters->parameters()[Acts::eBoundPhi] << std::endl; + // FM the line below is not correct, but it's what I'd like to do + //trackTip.addTrackState(result); + } + else { + } + + /* + auto atCalo_params = resultProp.value().endParameters->parameters(); + auto atCalo_covariance = resultProp.value().endParameters->covariance(); + + EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtCalorimeter, atCalo_params, atCalo_covariance, field[2] / Acts::UnitConstants::T); + track->trackStates().push_back(trackStateAtCalo); +*/ + std::reverse(hitsOnTrack.begin(), hitsOnTrack.end()); + std::reverse(statesOnTrack.begin(), statesOnTrack.end()); + + for (EVENT::TrackerHit* hit : hitsOnTrack) { + track->addHit(hit); + } + + if (statesOnTrack.size() > 0) { + dynamic_cast(statesOnTrack.back()) + ->setLocation(EVENT::TrackState::AtLastHit); + dynamic_cast(statesOnTrack.front()) + ->setLocation(EVENT::TrackState::AtFirstHit); + } + + EVENT::TrackStateVec& myTrackStates = track->trackStates(); + myTrackStates.insert(myTrackStates.end(), statesOnTrack.begin(), + statesOnTrack.end()); + + return track; +} + EVENT::TrackState* ACTS2Marlin_trackState( From 85a5b4fb5585817d0401bffcfd00882d8a1e485f Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 10 Jul 2024 18:36:45 +0200 Subject: [PATCH 3/5] added atCalo track state --- src/Helpers.cxx | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/src/Helpers.cxx b/src/Helpers.cxx index 39091e9..1c3edca 100644 --- a/src/Helpers.cxx +++ b/src/Helpers.cxx @@ -402,25 +402,20 @@ EVENT::Track* ACTS2Marlin_track( Stepper stepper(magneticField); Navigator navigator(navigatorCfg); Propagator mypropagator(std::move(stepper), std::move(navigator)); - auto resultProp = mypropagator.propagate(start, *caloSurface, myCaloPropOptions); - if (resultProp.value().endParameters.has_value()) { - std::cout << "Good!" << std::endl; - //FM: should read back the result of the propagation? - std::cout << "CaloPhi:" << phi << " " << resultProp.value().endParameters->parameters()[Acts::eBoundPhi] << std::endl; - // FM the line below is not correct, but it's what I'd like to do - //trackTip.addTrackState(result); + auto end_parameters = mypropagator.propagate(start, *caloSurface, myCaloPropOptions).value().endParameters; + if (end_parameters.has_value()) { + std::cout << "CaloPhi:" << phi << " " << end_parameters->parameters()[Acts::eBoundPhi] << std::endl; + const Acts::BoundMatrix& atCaloCovariance = *(end_parameters->covariance()); + + EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtCalorimeter, end_parameters->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T); + track->trackStates().push_back(trackStateAtCalo); } else { + std::cout << "Failed propagation!" << std::endl; } + - /* - auto atCalo_params = resultProp.value().endParameters->parameters(); - auto atCalo_covariance = resultProp.value().endParameters->covariance(); - - EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( - EVENT::TrackState::AtCalorimeter, atCalo_params, atCalo_covariance, field[2] / Acts::UnitConstants::T); - track->trackStates().push_back(trackStateAtCalo); -*/ std::reverse(hitsOnTrack.begin(), hitsOnTrack.end()); std::reverse(statesOnTrack.begin(), statesOnTrack.end()); From 815bff332344b0ebfa8b346116c2d99e1a57d123 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 10 Jul 2024 19:09:32 +0200 Subject: [PATCH 4/5] added atCalo track state --- src/Helpers.cxx | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/Helpers.cxx b/src/Helpers.cxx index 1c3edca..4afc535 100644 --- a/src/Helpers.cxx +++ b/src/Helpers.cxx @@ -389,7 +389,6 @@ EVENT::Track* ACTS2Marlin_track( // Set propagator options Acts::PropagatorOptions myCaloPropOptions(geoContext, magFieldContext); myCaloPropOptions.pathLimit = 20 * Acts::UnitConstants::m; - std::cout << "init prop option" << std::endl; //Let's try with our private propagator // Configurations @@ -402,9 +401,9 @@ EVENT::Track* ACTS2Marlin_track( Stepper stepper(magneticField); Navigator navigator(navigatorCfg); Propagator mypropagator(std::move(stepper), std::move(navigator)); - auto end_parameters = mypropagator.propagate(start, *caloSurface, myCaloPropOptions).value().endParameters; - if (end_parameters.has_value()) { - std::cout << "CaloPhi:" << phi << " " << end_parameters->parameters()[Acts::eBoundPhi] << std::endl; + auto resultProp = mypropagator.propagate(start, *caloSurface, myCaloPropOptions); + if (resultProp.ok()) { + auto end_parameters = resultProp.value().endParameters; const Acts::BoundMatrix& atCaloCovariance = *(end_parameters->covariance()); EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( @@ -415,7 +414,6 @@ EVENT::Track* ACTS2Marlin_track( std::cout << "Failed propagation!" << std::endl; } - std::reverse(hitsOnTrack.begin(), hitsOnTrack.end()); std::reverse(statesOnTrack.begin(), statesOnTrack.end()); From 80c3b6eff5934135f39edc22a04ce87fbaab7d4e Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 10 Jul 2024 23:13:53 +0200 Subject: [PATCH 5/5] added endcaps to extrapolation --- src/Helpers.cxx | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/src/Helpers.cxx b/src/Helpers.cxx index 4afc535..051b48f 100644 --- a/src/Helpers.cxx +++ b/src/Helpers.cxx @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -373,6 +374,14 @@ EVENT::Track* ACTS2Marlin_track( // Create the CaloSurface auto caloCylinder = std::make_shared(caloFaceR, caloFaceZ); auto caloSurface = Acts::Surface::makeShared(Acts::Transform3::Identity(), caloCylinder); + + // Define the circle dimensions (circles at both ends of the cylinder) + Acts::Translation3 circlePosition1(0, 0, -caloFaceZ); // circle at -z end + Acts::Translation3 circlePosition2(0, 0, caloFaceZ); // circle at +z end + + // Create the circle surfaces + auto circleSurface1 = Acts::Surface::makeShared(Acts::Transform3(circlePosition1), 0. ,caloFaceR); + auto circleSurface2 = Acts::Surface::makeShared(Acts::Transform3(circlePosition2), 0., caloFaceR); // define start parameters - swap this out with some smart call double d0 = params[Acts::eBoundLoc0]; @@ -383,7 +392,6 @@ EVENT::Track* ACTS2Marlin_track( double time = params[Acts::eBoundTime]; Acts::Vector3 pos(d0 * cos(phi), d0 * sin(phi), z0); - Acts::CurvilinearTrackParameters start(Acts::VectorHelpers::makeVector4(pos, time), phi, theta, qoverp, covariance, Acts::ParticleHypothesis::pion()); // Set propagator options @@ -411,7 +419,34 @@ EVENT::Track* ACTS2Marlin_track( track->trackStates().push_back(trackStateAtCalo); } else { - std::cout << "Failed propagation!" << std::endl; + if (theta > M_PI/2){ + auto resultPropP = mypropagator.propagate(start, *circleSurface1, myCaloPropOptions); + if (resultPropP.ok()) { + auto end_parametersP = resultPropP.value().endParameters; + const Acts::BoundMatrix& atCaloCovariance = *(end_parametersP->covariance()); + + EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtCalorimeter, end_parametersP->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T); + track->trackStates().push_back(trackStateAtCalo); + } + else{ + std::cout << "Failed propagation! " << std::endl; + } + } + else{ + auto resultPropM = mypropagator.propagate(start, *circleSurface2, myCaloPropOptions); + if (resultPropM.ok()) { + auto end_parametersM = resultPropM.value().endParameters; + const Acts::BoundMatrix& atCaloCovariance = *(end_parametersM->covariance()); + + EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState( + EVENT::TrackState::AtCalorimeter, end_parametersM->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T); + track->trackStates().push_back(trackStateAtCalo); + } + else{ + std::cout << "Failed propagation!" << std::endl; + } + } } std::reverse(hitsOnTrack.begin(), hitsOnTrack.end());