[Common] Use helix propagation to compute fwdtrack DCA in fwdtrackextension (#14282)

Co-authored-by: Maurice Coquet <mcoquet@alicecerno2.cern.ch>

Author: mcoquet642

Common/TableProducer/CMakeLists.txt

@@ -98,8 +98,8 @@ o2physics_add_dpl_workflow(fwdtrackextension
                     SOURCES fwdtrackextension.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                                           O2::ReconstructionDataFormats
-                                          O2::DetectorsBase
-                                          O2::DetectorsCommonDataFormats
+					  O2::GlobalTracking
+					  O2::MCHTracking
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(track-to-collision-associator

Common/TableProducer/fwdtrackextension.cxx

@@ -10,21 +10,30 @@
 // or submit itself to any jurisdiction.
 
 //
-// Task performing forward track DCA computation
+// \file fwdtrackextension.cxx
+// \brief Task performing forward track DCA computation.
+// \author Maurice Coquet, maurice.louis.coquet@cern.ch
 //
 
+#include "Common/Core/fwdtrackUtilities.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
+#include <CCDB/BasicCCDBManager.h>
+#include <DataFormatsParameters/GRPMagField.h>
+#include <DetectorsBase/GeometryManager.h>
+#include <DetectorsBase/Propagator.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
 #include <Framework/DataTypes.h>
 #include <Framework/runDataProcessing.h>
+#include <GlobalTracking/MatchGlobalFwd.h>
 #include <ReconstructionDataFormats/TrackFwd.h>
 
 #include <Math/MatrixRepresentationsStatic.h>
 #include <Math/SMatrix.h>
 
+#include <string>
 #include <vector>
 
 using namespace o2;
@@ -34,30 +43,62 @@ using namespace o2::framework::expressions;
 using SMatrix55 = ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<double, 5>>;
 using SMatrix5 = ROOT::Math::SVector<double, 5>;
 
+using MuonsWithCov = soa::Join<aod::FwdTracks, aod::FwdTracksCov>;
+
 struct FwdTrackExtension {
-  Produces<aod::FwdTracksDCA> extendedTrackQuantities;
+  Produces<aod::FwdTracksDCA> fwdDCA;
+  Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
+  Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
+  Configurable<std::string> configCcdbUrl{"configCcdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+  Configurable<bool> refitGlobalMuon{"refitGlobalMuon", true, "Recompute parameters of global muons"};
+
+  Service<o2::ccdb::BasicCCDBManager> fCCDB;
+  o2::parameters::GRPMagField* grpmag = nullptr; // for run 3, we access GRPMagField from GLO/Config/GRPMagField
+  int fCurrentRun;                               // needed to detect if the run changed and trigger update of magnetic field
 
-  void process(aod::FwdTracks const& tracks, aod::Collisions const&)
+  void init(o2::framework::InitContext&)
   {
-    for (auto& track : tracks) {
-      float dcaX = -999;
-      float dcaY = -999;
-      if (track.has_collision()) {
-        if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack || track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalForwardTrack || track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) {
+    // Load geometry
+    fCCDB->setURL(configCcdbUrl);
+    fCCDB->setCaching(true);
+    fCCDB->setLocalObjectValidityChecking();
 
-          auto const& collision = track.collision();
-          double chi2 = track.chi2();
-          SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt());
-          std::vector<double> v1;
-          SMatrix55 tcovs(v1.begin(), v1.end());
-          o2::track::TrackParCovFwd pars1{track.z(), tpars, tcovs, chi2};
-          pars1.propagateToZlinear(collision.posZ());
+    if (!o2::base::GeometryManager::isGeometryLoaded()) {
+      LOGF(info, "Load geometry from CCDB");
+      fCCDB->get<TGeoManager>(geoPath);
+    }
+  }
 
-          dcaX = (pars1.getX() - collision.posX());
-          dcaY = (pars1.getY() - collision.posY());
-        }
+  void process(aod::Collisions::iterator const& collision, o2::aod::BCsWithTimestamps const& /*...*/, MuonsWithCov const& tracks, aod::MFTTracks const& /*...*/)
+  {
+    auto bc = collision.template bc_as<o2::aod::BCsWithTimestamps>();
+    if (fCurrentRun != bc.runNumber()) {
+      grpmag = fCCDB->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, bc.timestamp());
+      if (grpmag != nullptr) {
+        LOGF(info, "Init field from GRP");
+        o2::base::Propagator::initFieldFromGRP(grpmag);
+      }
+      LOGF(info, "Set field for muons");
+      o2::mch::TrackExtrap::setField();
+      fCurrentRun = bc.runNumber();
+    }
+    const float zField = grpmag->getNominalL3Field();
+    for (const auto& track : tracks) {
+      const auto trackType = track.trackType();
+      o2::dataformats::GlobalFwdTrack fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwd(track, track);
+      if (refitGlobalMuon && (trackType == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack || trackType == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalForwardTrack)) {
+        auto muontrack = track.template matchMCHTrack_as<MuonsWithCov>();
+        auto mfttrack = track.template matchMFTTrack_as<aod::MFTTracks>();
+        o2::dataformats::GlobalFwdTrack propmuon = o2::aod::fwdtrackutils::propagateMuon(muontrack, muontrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, 0.f, zField);
+        SMatrix5 tpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt());
+        SMatrix55 tcovs{};
+        o2::track::TrackParCovFwd mft{mfttrack.z(), tpars, tcovs, mfttrack.chi2()};
+        fwdtrack = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuon, mft);
       }
-      extendedTrackQuantities(dcaX, dcaY);
+      const auto proptrack = o2::aod::fwdtrackutils::propagateTrackParCovFwd(fwdtrack, trackType, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, 0.f, zField);
+      const float dcaX = (proptrack.getX() - collision.posX());
+      const float dcaY = (proptrack.getY() - collision.posY());
+      fwdDCA(dcaX, dcaY);
     }
   }
 };