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/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 788ca75..376c51d 100644 --- a/src/ACTSSeededCKFTrackingProc.cxx +++ b/src/ACTSSeededCKFTrackingProc.cxx @@ -72,6 +72,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", @@ -682,7 +691,7 @@ void ACTSSeededCKFTrackingProc::processEvent(LCEvent *evt) { // 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..051b48f 100644 --- a/src/Helpers.cxx +++ b/src/Helpers.cxx @@ -6,10 +6,24 @@ #include #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 +315,161 @@ 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 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]; + 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; + + //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.ok()) { + auto end_parameters = resultProp.value().endParameters; + 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 { + 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()); + 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(