Extend time gain calibration

- allow for using TPC tracks from global track matches
- use DCA from global track in case global track match is required
- improve method for reading back a previous calibration

Author: wiechula

Detectors/TPC/calibration/src/CalibdEdx.cxx

@@ -754,17 +754,28 @@ void CalibdEdx::dumpToFile(const char* outFile)
 
 CalibdEdx CalibdEdx::readFromFile(const char* inFile)
 {
-  TFile f(inFile, "READ");
-  auto* obj = (CalibdEdx*)f.Get("calib");
+  std::unique_ptr<TFile> f(TFile::Open(inFile));
+  if (!f || f->IsZombie()) {
+    LOGP(error, "Could not open file: {}", inFile);
+    CalibdEdx calTmp;
+    return calTmp;
+  }
+
+  auto obj = f->Get<CalibdEdx>("calib");
   if (!obj) {
+    LOGP(error, "Could not read CalibdEdx object from file: {}", inFile);
     CalibdEdx calTmp;
     return calTmp;
   }
+
+  THnF* hTmp = f->Get<THnF>("histogram_data");
+
   CalibdEdx cal(*obj);
-  THnF* hTmp = (THnF*)f.Get("histogram_data");
+  delete obj;
+
   if (!hTmp) {
-    CalibdEdx calTmp;
-    return calTmp;
+    LOGP(warning, "Could not read histogram from file: {}. Returning empty histogram", inFile);
+    return cal;
   }
   cal.setFromRootHist(hTmp);
   return cal;

Detectors/TPC/workflow/CMakeLists.txt

@@ -200,7 +200,7 @@ o2_add_executable(idc-test-ft
 o2_add_executable(miptrack-filter
                   COMPONENT_NAME tpc
                   SOURCES src/tpc-miptrack-filter.cxx
-                  PUBLIC_LINK_LIBRARIES O2::TPCWorkflow)
+                  PUBLIC_LINK_LIBRARIES O2::TPCWorkflow O2::GlobalTrackingWorkflow)
 
 o2_add_executable(track-and-cluster-filter
                   COMPONENT_NAME tpc

Detectors/TPC/workflow/include/TPCWorkflow/MIPTrackFilterSpec.h

@@ -17,14 +17,16 @@
 #define O2_TPC_MIPTRACKFILTERSPEC_H_
 
 #include "Framework/DataProcessorSpec.h"
+#include "ReconstructionDataFormats/GlobalTrackID.h"
+using GID = o2::dataformats::GlobalTrackID;
 
 using namespace o2::framework;
 
 namespace o2::tpc
 {
 
 /// create a processor spec
-o2::framework::DataProcessorSpec getMIPTrackFilterSpec();
+o2::framework::DataProcessorSpec getMIPTrackFilterSpec(GID::mask_t srcTracks = GID::getSourcesMask("TPC"));
 
 } // namespace o2::tpc
 

Detectors/TPC/workflow/src/MIPTrackFilterSpec.cxx

@@ -16,15 +16,14 @@
 #include "TPCWorkflow/MIPTrackFilterSpec.h"
 
 #include <algorithm>
-#include <iterator>
 #include <vector>
 #include <memory>
 #include <random>
 
 // o2 includes
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/TrackCuts.h"
-#include "DetectorsCalibration/Utils.h"
+#include "Framework/CCDBParamSpec.h"
 #include "Framework/Logger.h"
 #include "DetectorsBase/GRPGeomHelper.h"
 #include "Framework/Task.h"
@@ -33,16 +32,23 @@
 #include "Framework/ConfigParamRegistry.h"
 #include "TPCWorkflow/ProcessingHelpers.h"
 #include "Headers/DataHeader.h"
+#include "DataFormatsGlobalTracking/RecoContainer.h"
+#include "ReconstructionDataFormats/PrimaryVertex.h"
+#include "DataFormatsCalibration/MeanVertexObject.h"
+#include "ReconstructionDataFormats/VtxTrackRef.h"
 
 using namespace o2::framework;
+using DataRequest = o2::globaltracking::DataRequest;
+using GID = o2::dataformats::GlobalTrackID;
 
 namespace o2::tpc
 {
 
 class MIPTrackFilterDevice : public Task
 {
  public:
-  MIPTrackFilterDevice(std::shared_ptr<o2::base::GRPGeomRequest> gr) : mGRPGeomRequest(gr) {}
+  MIPTrackFilterDevice(std::shared_ptr<o2::base::GRPGeomRequest> gr, std::shared_ptr<DataRequest> dr, GID::mask_t trackSourcesMask)
+    : mGRPGeomRequest(gr), mDataRequest(dr), mTrackSourcesMask(trackSourcesMask) {}
 
   void init(framework::InitContext& ic) final;
   void run(ProcessingContext& pc) final;
@@ -53,16 +59,20 @@ class MIPTrackFilterDevice : public Task
   void sendOutput(DataAllocator& output);
 
   std::shared_ptr<o2::base::GRPGeomRequest> mGRPGeomRequest;
-  TrackCuts mCuts{};                  ///< Tracks cuts object
-  std::vector<TrackTPC> mMIPTracks;   ///< Filtered MIP tracks
-  unsigned int mProcessEveryNthTF{1}; ///< process every Nth TF only
-  int mMaxTracksPerTF{-1};            ///< max number of MIP tracks processed per TF
-  uint32_t mTFCounter{0};             ///< counter to keep track of the TFs
-  int mProcessNFirstTFs{0};           ///< number of first TFs which are not sampled
-  float mDCACut{-1};                  ///< DCA cut
-  bool mSendDummy{false};             ///< send empty data in case TF is skipped
-
-  bool acceptDCA(const TrackTPC& track);
+  std::shared_ptr<DataRequest> mDataRequest;
+  GID::mask_t mTrackSourcesMask;
+  TrackCuts mCuts{};                      ///< Tracks cuts object
+  std::vector<TrackTPC> mMIPTracks;       ///< Filtered MIP tracks
+  o2::dataformats::MeanVertexObject mVtx; ///< Mean vertex object
+  unsigned int mProcessEveryNthTF{1};     ///< process every Nth TF only
+  int mMaxTracksPerTF{-1};                ///< max number of MIP tracks processed per TF
+  uint32_t mTFCounter{0};                 ///< counter to keep track of the TFs
+  int mProcessNFirstTFs{0};               ///< number of first TFs which are not sampled
+  float mDCACut{-1};                      ///< DCA cut
+  float mDCAZCut{-1};                     ///< DCA z cut
+  bool mSendDummy{false};                 ///< send empty data in case TF is skipped
+
+  bool acceptDCA(o2::track::TrackPar propTrack, o2::math_utils::Point3D<float> refPoint, bool useDCAz = false);
 };
 
 void MIPTrackFilterDevice::init(framework::InitContext& ic)
@@ -100,13 +110,16 @@ void MIPTrackFilterDevice::init(framework::InitContext& ic)
   mCuts.setCutLooper(cutLoopers);
 
   mDCACut = ic.options().get<float>("dca-cut");
+  mDCAZCut = ic.options().get<float>("dca-z-cut");
 
   o2::base::GRPGeomHelper::instance().setRequest(mGRPGeomRequest);
 }
 
 void MIPTrackFilterDevice::run(ProcessingContext& pc)
 {
   o2::base::GRPGeomHelper::instance().checkUpdates(pc);
+  pc.inputs().get<o2::dataformats::MeanVertexObject*>("meanvtx");
+
   const auto currentTF = processing_helpers::getCurrentTF(pc);
   if ((mTFCounter++ % mProcessEveryNthTF) && (currentTF >= mProcessNFirstTFs)) {
     LOGP(info, "Skipping TF {}", currentTF);
@@ -117,19 +130,61 @@ void MIPTrackFilterDevice::run(ProcessingContext& pc)
     return;
   }
 
-  const auto tracks = pc.inputs().get<gsl::span<TrackTPC>>("tracks");
-  const auto nTracks = tracks.size();
 
-  if ((mMaxTracksPerTF != -1) && (nTracks > mMaxTracksPerTF)) {
-    // indices to good tracks
-    std::vector<size_t> indices;
-    indices.reserve(nTracks);
+  o2::globaltracking::RecoContainer recoData;
+  recoData.collectData(pc, *mDataRequest);
+  const auto tracksTPC = recoData.getTPCTracks();
+  const auto nTracks = tracksTPC.size();
+
+  // indices to good tracks
+  std::vector<size_t> indices;
+  indices.reserve(nTracks);
+
+  const auto useGlobalTracks = mTrackSourcesMask[GID::ITSTPC];
+  o2::math_utils::Point3D<float> vertex = mVtx.getXYZ();
+
+  if (useGlobalTracks) {
+    auto trackIndex = recoData.getPrimaryVertexMatchedTracks();                      // Global ID's for associated tracks
+    auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs();                      // references from vertex to these track IDs
+    std::vector<GID::Source> selSrc{GID::ITSTPC, GID::ITSTPCTRD, GID::ITSTPCTRDTOF}; // for Instance
+    // LOGP(info, "Number of vertex tracks: {}", vtxRefs.size());
+    const auto nv = (vtxRefs.size()>0) ? vtxRefs.size() - 1 : 0;                     // note: the last entry groups the tracks which were not related to any vertex, to skip them, use vtxRefs.size()-1
+
+    for (int iv = 0; iv < nv; iv++) {
+      const auto& vtref = vtxRefs[iv];
+      //LOGP(info, "Processing vertex {} with {} tracks", iv, vtref.getEntries());
+      vertex = recoData.getPrimaryVertex(iv).getXYZ();
+      //LOGP(info, "Vertex position: x={} y={} z={}", vertex.x(), vertex.y(), vertex.z());
+
+      for (auto src : selSrc) {
+        int idMin = vtxRefs[iv].getFirstEntryOfSource(src), idMax = idMin + vtxRefs[iv].getEntriesOfSource(src);
+        //LOGP(info, "Source {}: idMin={} idMax={}", GID::getSourceName(src), idMin, idMax);
+
+        for (int i = idMin; i < idMax; i++) {
+          auto vid = trackIndex[i];
+          const auto& track = recoData.getTrackParam(vid); // this is a COPY of the track param which we will modify during DCA calculation
+          auto gidTPC = recoData.getTPCContributorGID(vid);
+          if (gidTPC.isSourceSet()) {
+            const auto idxTPC = gidTPC.getIndex();
+            if (mCuts.goodTrack(tracksTPC[idxTPC]) && acceptDCA(tracksTPC[idxTPC], vertex, true)) {
+              indices.emplace_back(idxTPC);
+            }
+          }
+        }
+      }
+    }
+
+  } else {
     for (size_t i = 0; i < nTracks; ++i) {
-      if (mCuts.goodTrack(tracks[i]) && acceptDCA(tracks[i])) {
+      if (mCuts.goodTrack(tracksTPC[i]) && acceptDCA(tracksTPC[i], vertex)) {
         indices.emplace_back(i);
       }
     }
+  }
+
+  size_t nTracksSel = indices.size();
 
+  if ((mMaxTracksPerTF != -1) && (nTracksSel > mMaxTracksPerTF)) {
     // in case no good tracks have been found
     if (indices.empty()) {
       mMIPTracks.clear();
@@ -144,15 +199,14 @@ void MIPTrackFilterDevice::run(ProcessingContext& pc)
     std::shuffle(indices.begin(), indices.end(), rng);
 
     // copy good tracks
-    const int loopEnd = (mMaxTracksPerTF > indices.size()) ? indices.size() : mMaxTracksPerTF;
-    for (int i = 0; i < loopEnd; ++i) {
-      mMIPTracks.emplace_back(tracks[indices[i]]);
-    }
-  } else {
-    std::copy_if(tracks.begin(), tracks.end(), std::back_inserter(mMIPTracks), [this](const auto& track) { return mCuts.goodTrack(track) && acceptDCA(track); });
+    nTracksSel = (mMaxTracksPerTF > indices.size()) ? indices.size() : mMaxTracksPerTF;
+  }
+
+  for (int i = 0; i < nTracksSel; ++i) {
+    mMIPTracks.emplace_back(tracksTPC[indices[i]]);
   }
 
-  LOGP(info, "Filtered {} MIP tracks out of {} total tpc tracks", mMIPTracks.size(), tracks.size());
+  LOGP(info, "Filtered {} / {} MIP tracks out of {} total tpc tracks, using {}", mMIPTracks.size(), indices.size(), tracksTPC.size(), useGlobalTracks ? "global tracks" : "TPC only tracks");
   sendOutput(pc.outputs());
   mMIPTracks.clear();
 }
@@ -162,6 +216,11 @@ void MIPTrackFilterDevice::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
   if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
     return;
   }
+  if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) {
+    LOG(info) << "Setting new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString();
+    mVtx = *(const o2::dataformats::MeanVertexObject*)obj;
+    return;
+  }
 }
 
 void MIPTrackFilterDevice::sendOutput(DataAllocator& output) { output.snapshot(Output{header::gDataOriginTPC, "MIPS", 0}, mMIPTracks); }
@@ -171,44 +230,46 @@ void MIPTrackFilterDevice::endOfStream(EndOfStreamContext& eos)
   LOG(info) << "Finalizig MIP Tracks filter";
 }
 
-bool MIPTrackFilterDevice::acceptDCA(const TrackTPC& track)
+bool MIPTrackFilterDevice::acceptDCA(o2::track::TrackPar propTrack, o2::math_utils::Point3D<float> refPoint, bool useDCAz)
 {
   if (mDCACut < 0) {
     return true;
   }
 
   auto propagator = o2::base::Propagator::Instance();
   std::array<float, 2> dca;
-  const o2::math_utils::Point3D<float> refPoint{0, 0, 0};
-  o2::track::TrackPar propTrack(track);
   const auto ok = propagator->propagateToDCABxByBz(refPoint, propTrack, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca);
   const auto dcar = std::abs(dca[0]);
 
-  return ok && (dcar < mDCACut);
+  return ok && (dcar < mDCACut) && (!useDCAz || (std::abs(dca[1]) < mDCAZCut));
 }
 
-DataProcessorSpec getMIPTrackFilterSpec()
+DataProcessorSpec getMIPTrackFilterSpec(GID::mask_t srcTracks)
 {
   std::vector<OutputSpec> outputs;
   outputs.emplace_back(header::gDataOriginTPC, "MIPS", 0, Lifetime::Sporadic);
 
-  std::vector<InputSpec> inputs;
-  inputs.emplace_back("tracks", "TPC", "TRACKS");
+  const auto useMC = false;
+  auto dataRequest = std::make_shared<DataRequest>();
+  dataRequest->requestTracks(srcTracks, useMC);
+  dataRequest->requestPrimaryVertices(useMC);
 
   auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false,                             // orbitResetTime
                                                               true,                              // GRPECS=true
                                                               false,                             // GRPLHCIF
                                                               true,                              // GRPMagField
                                                               true,                              // askMatLUT
                                                               o2::base::GRPGeomRequest::Aligned, // geometry
-                                                              inputs,
+                                                              dataRequest->inputs,
                                                               true);
 
+  dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, o2::framework::ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1));
+
   return DataProcessorSpec{
     "tpc-miptrack-filter",
-    inputs,
+    dataRequest->inputs,
     outputs,
-    adaptFromTask<MIPTrackFilterDevice>(ggRequest),
+    adaptFromTask<MIPTrackFilterDevice>(ggRequest, dataRequest, srcTracks),
     Options{
       {"min-momentum", VariantType::Double, 0.35, {"minimum momentum cut"}},
       {"max-momentum", VariantType::Double, 0.55, {"maximum momentum cut"}},
@@ -220,7 +281,8 @@ DataProcessorSpec getMIPTrackFilterSpec()
       {"process-first-n-TFs", VariantType::Int, 1, {"Number of first TFs which are not sampled"}},
       {"send-dummy-data", VariantType::Bool, false, {"Send empty data in case TF is skipped"}},
       {"dont-cut-loopers", VariantType::Bool, false, {"Do not cut loopers by comparing zout-zin"}},
-      {"dca-cut", VariantType::Float, 3.f, {"DCA cut in cm, < 0 to disable"}},
+      {"dca-cut", VariantType::Float, 3.f, {"DCA cut in xy (cm), < 0 to disable cut in xy and z"}},
+      {"dca-z-cut", VariantType::Float, 5.f, {"DCA cut in z (cm)"}},
     }};
 }
 

Detectors/TPC/workflow/src/tpc-miptrack-filter.cxx

@@ -13,6 +13,10 @@
 #include "TPCWorkflow/MIPTrackFilterSpec.h"
 #include "Framework/ConfigParamSpec.h"
 #include "DataFormatsTPC/TrackTPC.h"
+#include "GlobalTrackingWorkflowHelpers/InputHelper.h"
+#include "ReconstructionDataFormats/GlobalTrackID.h"
+
+using GID = o2::dataformats::GlobalTrackID;
 
 template <typename T>
 using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition<T>;
@@ -21,6 +25,8 @@ void customize(std::vector<ConfigParamSpec>& workflowOptions)
 {
   std::vector<ConfigParamSpec> options{
     {"enable-writer", VariantType::Bool, false, {"selection string input specs"}},
+    {"use-global-tracks", VariantType::Bool, false, {"use global matched tracks instead of TPC only"}},
+    {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}},
   };
 
   std::swap(workflowOptions, options);
@@ -33,8 +39,17 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config)
 {
   using namespace o2::tpc;
 
+  const auto useGlobal = config.options().get<bool>("use-global-tracks");
   WorkflowSpec workflow;
-  workflow.emplace_back(getMIPTrackFilterSpec());
+
+  const auto useMC = false;
+  auto srcTracks = GID::getSourcesMask("TPC");
+  const auto srcCls = GID::getSourcesMask("");
+  if (useGlobal) {
+    srcTracks = GID::getSourcesMask("ITS,TPC,ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF");
+  }
+
+  workflow.emplace_back(getMIPTrackFilterSpec(srcTracks));
 
   if (config.options().get<bool>("enable-writer")) {
     const char* processName = "tpc-mips-writer";