[WIP] Adds: O2 DCA histograms in tracks task

Author: Bhawani Singh

Detectors/TPC/qc/include/TPCQC/Tracks.h

@@ -21,6 +21,8 @@
 #include <string_view>
 #include <unordered_map>
 #include <memory>
+#include "CCDB/CcdbApi.h"
+#include "CCDB/BasicCCDBManager.h"
 
 // root includes
 #include "TH1.h"
@@ -77,6 +79,13 @@ class Tracks
     mCutMinnCls = nClusterCut;
     mCutMindEdxTot = dEdxTot;
   }
+
+  // To set the runNumber
+  void setRunNumber(int RunNumber)
+  {
+    runNumber = RunNumber;
+  }
+
   // Just for backward compatibility with crrent QC, temporary, will be removed in the next PR
   /// get 1D histograms
   std::vector<TH1F>& getHistograms1D() { return mHist1D; }
@@ -100,6 +109,7 @@ class Tracks
   float mCutAbsEta = 1.f;      // Eta cut
   int mCutMinnCls = 60;        // minimum N clusters
   float mCutMindEdxTot = 20.f; // dEdxTot min value
+  int runNumber = -999;        // runNumber
   std::unordered_map<std::string_view, std::unique_ptr<TH1>> mMapHist;
   std::vector<TH1F> mHist1D{};      ///< Initialize vector of 1D histograms
   std::vector<TH2F> mHist2D{};      ///< Initialize vector of 2D histograms

Detectors/TPC/qc/src/Tracks.cxx

@@ -21,13 +21,18 @@
 // o2 includes
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/dEdxInfo.h"
+#include "GPUCommonArray.h"
+#include "DetectorsBase/Propagator.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "TGeoManager.h"
 #include "TPCQC/Tracks.h"
 #include "TPCQC/Helpers.h"
 
 ClassImp(o2::tpc::qc::Tracks);
 
 using namespace o2::tpc::qc;
-
+// DCA histograms
+const std::vector<std::string_view> types{"A_Pos", "A_Neg", "C_Pos", "C_Neg"};
 //______________________________________________________________________________
 void Tracks::initializeHistograms()
 {
@@ -78,6 +83,12 @@ void Tracks::initializeHistograms()
   mMapHist["hPhiAsideRatio"] = std::make_unique<TH1F>("hPhiAsideRatio", "Azimuthal angle, A side, ratio neg./pos. ;phi", 360, 0., 2 * M_PI);
   mMapHist["hPhiCsideRatio"] = std::make_unique<TH1F>("hPhiCsideRatio", "Azimuthal angle, C side, ratio neg./pos. ;phi", 360, 0., 2 * M_PI);
   mMapHist["hPtRatio"] = std::make_unique<TH1F>("hPtRatio", "Transverse momentum, ratio neg./pos. ;p_T", logPtBinning.size() - 1, logPtBinning.data());
+
+  // DCA Histograms
+  for (const auto type : types) {
+  mMapHist[fmt::format("hDCAr_{}", type).data()] = std::make_unique<TH2F>(fmt::format("hDCAr_{}", type).data(), fmt::format("DCAr {};phi;DCAr (cm)", type).data(), 360, 0, o2::math_utils::twoPid(), 100, -3., 3.);
+  mMapHist[fmt::format("hDCAz_{}", type).data()] = std::make_unique<TH2F>(fmt::format("hDCAz_{}", type).data(), fmt::format("DCAz {};phi;DCAr (cm)", type).data(), 360, 0, o2::math_utils::twoPid(), 100, -3., 3.);
+  }
 }
 //______________________________________________________________________________
 void Tracks::resetHistograms()
@@ -117,6 +128,45 @@ bool Tracks::processTrack(const o2::tpc::TrackTPC& track)
     mMapHist["hNClustersAfterCuts"]->Fill(nCls);
     mMapHist["hEta"]->Fill(eta);
 
+  auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
+  ccdbmgr.setURL("https://alice-ccdb.cern.ch");
+  auto runDuration = ccdbmgr.getRunDuration(runNumber);
+  auto tRun = runDuration.first + (runDuration.second - runDuration.first) / 2; // time stamp for the middle of the run duration
+  ccdbmgr.setTimestamp(tRun);
+
+  // CTP orbit reset time
+  auto orbitResetTimeNS = ccdbmgr.get<std::vector<int64_t>>("CTP/Calib/OrbitReset");
+  int64_t orbitResetTimeMS = (*orbitResetTimeNS)[0] * 1e-3;
+  LOGP(info, "Orbit reset time in MS is {}", orbitResetTimeMS);
+
+    auto geoAligned = ccdbmgr.get<TGeoManager>("GLO/Config/GeometryAligned");
+    auto magField = ccdbmgr.get<o2::parameters::GRPMagField>("GLO/Config/GRPMagField");
+    const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdbmgr.get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
+    auto propagator = o2::base::Propagator::Instance();
+    propagator->setMatLUT(matLut);
+      //---| propagate to 0,0,0 |---
+      o2::gpu::gpustd::array<float, 2> dca;
+      const o2::math_utils::Point3D<float> refPoint{0, 0, 0};
+      //auto propTrack = TrackTPC(track);
+      o2::track::TrackPar propTrack(track); // Should be cheaper than the one above
+      bool useThisTrack = true;
+      if (!propagator->propagateToDCABxByBz(refPoint, propTrack, 0.5, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) {
+        useThisTrack = false;
+      }
+     ///fine grained propagation most probably not needed
+       if (!propagator->propagateToDCABxByBz(refPoint, propTrack, 0.005, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) {
+        useThisTrack = false;
+      }
+
+      if (useThisTrack)
+      {
+      // ---| fill histos |---
+      const int type = (track.getQ2Pt() < 0) + 2 * track.hasCSideClustersOnly();
+      const auto phi = o2::math_utils::to02PiGen(track.getPhi());
+      mMapHist[fmt::format("hDCAr_{}", types[type]).data()]->Fill(phi, dca[0]);
+      mMapHist[fmt::format("hDCAz_{}", types[type]).data()]->Fill(phi, dca[1]);
+      }
+
     if (hasASideOnly == 1) {
       mMapHist["hPhiAside"]->Fill(phi);
     } else if (hasCSideOnly == 1) {