Add extra info with charge and timing to unbinned residuals

Author: shahor02

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h

@@ -128,8 +128,9 @@ class ResidualAggregatorDevice : public o2::framework::Task
     updateTimeDependentParams(pc);
     std::chrono::duration<double, std::milli> ccdbUpdateTime = std::chrono::high_resolution_clock::now() - runStartTime;
 
-    // we always require the unbinned residuals and the associated track references
+    // we always require the unbinned residuals and the associated detector info and track references
     auto residualsData = pc.inputs().get<gsl::span<o2::tpc::UnbinnedResid>>("unbinnedRes");
+    auto residualsDataDet = pc.inputs().get<gsl::span<o2::tpc::DetInfoResid>>("detinfoRes");
     auto trackRefs = pc.inputs().get<gsl::span<o2::tpc::TrackDataCompact>>("trackRefs");
 
     // track data input is optional
@@ -151,7 +152,7 @@ class ResidualAggregatorDevice : public o2::framework::Task
 
     o2::base::TFIDInfoHelper::fillTFIDInfo(pc, mAggregator->getCurrentTFInfo());
     LOG(detail) << "Processing TF " << mAggregator->getCurrentTFInfo().tfCounter << " with " << trkData->size() << " tracks and " << residualsData.size() << " unbinned residuals associated to them";
-    mAggregator->process(residualsData, trackRefs, trkDataPtr, lumi);
+    mAggregator->process(residualsData, residualsDataDet, trackRefs, trkDataPtr, lumi);
     std::chrono::duration<double, std::milli> runDuration = std::chrono::high_resolution_clock::now() - runStartTime;
     LOGP(debug, "Duration for run method: {} ms. From this taken for time dependent param update: {} ms",
          std::chrono::duration_cast<std::chrono::milliseconds>(runDuration).count(),
@@ -222,6 +223,7 @@ DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, b
   auto& inputs = dataRequest->inputs;
   o2::tpc::VDriftHelper::requestCCDBInputs(inputs);
   inputs.emplace_back("unbinnedRes", "GLO", "UNBINNEDRES");
+  inputs.emplace_back("detinfoRes", "GLO", "DETINFORES");
   inputs.emplace_back("trackRefs", "GLO", "TRKREFS");
   if (trackInput) {
     inputs.emplace_back("trkData", "GLO", "TRKDATA");

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCUnbinnedResidualReaderSpec.h

@@ -43,6 +43,7 @@ class TPCUnbinnedResidualReader : public o2::framework::Task
   std::string mInFileName;
   std::string mInTreeName;
   std::vector<UnbinnedResid> mUnbinnedResid, *mUnbinnedResidPtr = &mUnbinnedResid;
+  std::vector<DetInfoResid> mDetInfoUnbRes, *mDetInfoUnbResPtr = &mDetInfoUnbRes;
   std::vector<TrackData> mTrackData, *mTrackDataPtr = &mTrackData;
   std::vector<TrackDataCompact> mTrackDataCompact, *mTrackDataCompactPtr = &mTrackDataCompact;
 };

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx

@@ -143,6 +143,7 @@ void TPCInterpolationDPL::run(ProcessingContext& pc)
     }
   }
   pc.outputs().snapshot(Output{"GLO", "UNBINNEDRES", 0}, mInterpolation.getClusterResiduals());
+  pc.outputs().snapshot(Output{"GLO", "DETINFORES", 0}, mInterpolation.getClusterResidualsDetInfo());
   pc.outputs().snapshot(Output{"GLO", "TRKREFS", 0}, mInterpolation.getTrackDataCompact());
   if (mSendTrackData) {
     pc.outputs().snapshot(Output{"GLO", "TRKDATA", 0}, mInterpolation.getReferenceTracks());
@@ -188,6 +189,7 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas
     }
   }
   outputs.emplace_back("GLO", "UNBINNEDRES", 0, Lifetime::Timeframe);
+  outputs.emplace_back("GLO", "DETINFORES", 0, Lifetime::Timeframe);
   outputs.emplace_back("GLO", "TRKREFS", 0, Lifetime::Timeframe);
   if (sendTrackData) {
     outputs.emplace_back("GLO", "TRKDATA", 0, Lifetime::Timeframe);

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCUnbinnedResidualReaderSpec.cxx

@@ -44,6 +44,11 @@ void TPCUnbinnedResidualReader::connectTree()
   assert(mTreeIn);
   mTreeIn->SetBranchAddress("residuals", &mUnbinnedResidPtr);
   mTreeIn->SetBranchAddress("trackRefs", &mTrackDataCompactPtr);
+  if (mTreeIn->GetBranch("detInfo")) {
+    mTreeIn->SetBranchAddress("detInfo", &mDetInfoUnbResPtr);
+  } else {
+    LOGP(warn, "No detInfo branch found in the unbinned residuals tree, empty vector will be sent");
+  }
   if (mTrackInput) {
     mTreeIn->SetBranchAddress("tracks", &mTrackDataPtr);
   }
@@ -58,6 +63,7 @@ void TPCUnbinnedResidualReader::run(ProcessingContext& pc)
   LOG(info) << "Pushing " << mUnbinnedResid.size() << " unbinned residuals at entry " << currEntry;
   pc.outputs().snapshot(Output{"GLO", "UNBINNEDRES", 0}, mUnbinnedResid);
   pc.outputs().snapshot(Output{"GLO", "TRKREFS", 0}, mTrackDataCompact);
+  pc.outputs().snapshot(Output{"GLO", "DETINFORES", 0}, mDetInfoUnbRes);
   if (mTrackInput) {
     LOG(info) << "Pushing " << mTrackData.size() << " reference tracks for these residuals";
     pc.outputs().snapshot(Output{"GLO", "TRKDATA", 0}, mTrackData);
@@ -73,6 +79,7 @@ DataProcessorSpec getUnbinnedTPCResidualsReaderSpec(bool trkInput)
 {
   std::vector<OutputSpec> outputs;
   outputs.emplace_back("GLO", "UNBINNEDRES", 0, Lifetime::Timeframe);
+  outputs.emplace_back("GLO", "DETINFORES", 0, Lifetime::Timeframe);
   outputs.emplace_back("GLO", "TRKREFS", 0, Lifetime::Timeframe);
   if (trkInput) {
     outputs.emplace_back("GLO", "TRKDATA", 0, Lifetime::Timeframe);

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h

@@ -49,7 +49,7 @@ struct ResidualsContainer {
   void fillStatisticsBranches();
   uint64_t getNEntries() const { return nResidualsTotal; }
 
-  void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span<const UnbinnedResid> resid, const gsl::span<const o2::tpc::TrackDataCompact> trkRefsIn, const gsl::span<const o2::tpc::TrackData>* trkDataIn, const o2::ctp::LumiInfo* lumiInput);
+  void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span<const UnbinnedResid> resid, const gsl::span<const DetInfoResid> detInfoRes, const gsl::span<const o2::tpc::TrackDataCompact> trkRefsIn, const gsl::span<const o2::tpc::TrackData>* trkDataIn, const o2::ctp::LumiInfo* lumiInput);
   void merge(ResidualsContainer* prev);
   void print();
   void writeToFile(bool closeFileAfterwards);
@@ -64,6 +64,7 @@ struct ResidualsContainer {
   std::vector<uint32_t> sumUnbinnedResid, *sumUnbinnedResidPtr{&sumUnbinnedResid}; ///< sum of unbinned residuals for each TF
   std::vector<o2::ctp::LumiInfo> lumi, *lumiPtr{&lumi};                      ///< luminosity information from CTP per TF
   std::vector<UnbinnedResid> unbinnedRes, *unbinnedResPtr{&unbinnedRes};     ///< unbinned residuals which are sent to the aggregator
+  std::vector<DetInfoResid> detInfoUnbRes, *detInfoUnbResPtr{&detInfoUnbRes}; ///< detector info associated to unbinned residuals which are sent to the aggregator
   std::vector<TrackData> trkData, *trkDataPtr{&trkData};                     ///< track data and cluster ranges
   std::vector<TrackDataCompact> trackInfo, *trackInfoPtr{&trackInfo};        ///< allows to obtain track type for each unbinned residual downstream
   o2::ctp::LumiInfo lumiTF;                                                  ///< for each processed TF we store the lumi information in the tree of unbinned residuals

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h

@@ -101,6 +101,50 @@ struct UnbinnedResid {
   ClassDefNV(UnbinnedResid, 2);
 };
 
+struct DetInfoResid { // detector info associated with residual
+  uint32_t word = 0;  // container interpreted in a different way depending on the detector type
+  //
+  // TPC view: qTot and qMax of the cluster
+  uint16_t qTotTPC() const { return static_cast<uint16_t>(word & 0xFFFFu); }
+  uint16_t qMaxTPC() const { return static_cast<uint16_t>((word >> 16) & 0xFFFFu); }
+  void setTPC(uint16_t qTot, uint16_t qMax) { word = (static_cast<uint32_t>(qMax) << 16) | static_cast<uint32_t>(qTot); }
+  //
+  // TRD view: q0, q1, q2 + calibrated slope (truncated to in +-3.5 range)
+  static constexpr uint32_t TRDQ0NB = 7, TRDQ1NB = 7, TRDQ2NB = 6, TRDSlpNB = 12;
+  static constexpr uint32_t TRDQ0Msk = (1 << TRDQ0NB) - 1, TRDQ1Msk = (1 << TRDQ1NB) - 1, TRDQ2Msk = ((1 << TRDQ2NB) - 1), TRDSlpMsk = (1 << TRDSlpNB) - 1;
+  static constexpr float TRDMaxSlope = 3.5, TRDSlope2Int = ((1 << TRDSlpNB) - 1) / (2 * TRDMaxSlope), TRDInt2Slope = 1.f / TRDSlope2Int;
+  uint16_t q0TRD() const { return static_cast<uint16_t>(word & TRDQ0Msk); }
+  uint16_t q1TRD() const { return static_cast<uint16_t>((word >> TRDQ0NB) & TRDQ1Msk); }
+  uint16_t q2TRD() const { return static_cast<uint16_t>((word >> (TRDQ0NB + TRDQ1NB)) & TRDQ2Msk); }
+  float slopeTRD() const { return ((word >> (TRDQ0NB + TRDQ1NB + TRDQ2NB)) & TRDSlpMsk) * TRDInt2Slope - TRDMaxSlope; }
+  void setTRD(uint8_t q0, uint8_t q1, uint8_t q2, float slope)
+  {
+    float rslope = (slope + TRDMaxSlope) * TRDSlope2Int;
+    if (rslope < 0.f) {
+      rslope = 0;
+    } else if (rslope > TRDSlpMsk) {
+      rslope = TRDSlpMsk;
+    }
+    uint32_t slpI = std::round(rslope);
+    word = (static_cast<uint32_t>(slpI << (TRDQ0NB + TRDQ1NB + TRDQ2NB)) |
+            static_cast<uint32_t>((q2 & TRDQ2Msk) << (TRDQ0NB + TRDQ1NB)) |
+            static_cast<uint32_t>((q1 & TRDQ1Msk) << TRDQ0NB) |
+            static_cast<uint32_t>(q0 & TRDQ0Msk));
+  }
+  //
+  // TOF view (time difference in \mus wrt seeding ITS-TPC track)
+  float timeTOF() const { return std::bit_cast<float>(word); }
+  void setTOF(float t) { word = std::bit_cast<uint32_t>(t); }
+  //
+  // No info for ITS is stored
+  //
+  // PV view (time difference in \mus wrt contributing ITS-TPC track)
+  float timePV() const { return std::bit_cast<float>(word); }
+  void setPV(float t) { word = std::bit_cast<uint32_t>(t); }
+
+  ClassDefNV(DetInfoResid, 1);
+};
+
 /// Structure for the information required to associate each residual with a given track type (ITS-TPC-TRD-TOF, etc)
 struct TrackDataCompact {
   TrackDataCompact() = default;
@@ -296,10 +340,13 @@ class TrackInterpolation
 
   void setExtDetResid(bool v) { mExtDetResid = v; }
 
-  int processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TrackParCov& trkWork, std::array<float, 2>* trkltTRDYZ = nullptr, std::array<float, 3>* trkltTRDCov = nullptr, TrackData* trkData = nullptr);
+  int processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TrackParCov& trkWork, std::array<float, 2>* trkltTRDYZ = nullptr,
+                      std::array<float, 3>* trkltTRDCov = nullptr, TrackData* trkData = nullptr,
+                      o2::trd::Tracklet64* trk64 = nullptr, o2::trd::CalibratedTracklet* trkCalib = nullptr);
 
   // --------------------------------- output ---------------------------------------------
   std::vector<UnbinnedResid>& getClusterResiduals() { return mClRes; }
+  std::vector<DetInfoResid>& getClusterResidualsDetInfo() { return mDetInfoRes; }
   std::vector<TrackDataCompact>& getTrackDataCompact() { return mTrackDataCompact; }
   std::vector<TrackDataExtended>& getTrackDataExtended() { return mTrackDataExtended; }
   std::vector<TrackData>& getReferenceTracks() { return mTrackData; }
@@ -348,6 +395,7 @@ class TrackInterpolation
   std::vector<TrackDataCompact> mTrackDataCompact{};   ///< required to connect each residual to a global track
   std::vector<TrackDataExtended> mTrackDataExtended{}; ///< full tracking information for debugging
   std::vector<UnbinnedResid> mClRes{};                 ///< residuals for each available TPC cluster of all tracks
+  std::vector<DetInfoResid> mDetInfoRes{};             ///< packed detector info associated with each residual
   std::vector<TrackData> mTrackDataUnfiltered{};       ///< same as mTrackData, but for all tracks before outlier filtering
   std::vector<TPCClusterResiduals> mClResUnfiltered{}; ///< same as mClRes, but for all residuals before outlier filtering
 

Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx

@@ -124,6 +124,7 @@ void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string
     treeOutResidualsUnbinned->Branch("trackInfo", &trackInfoPtr);
     treeOutResidualsUnbinned->Branch("CTPLumi", &lumiTF);
     treeOutResidualsUnbinned->Branch("timeMS", &timeMS);
+    treeOutResidualsUnbinned->Branch("detInfo", &detInfoUnbResPtr);
   }
   if (writeTrackData) {
     treeOutTrackData = std::make_unique<TTree>("trackData", "Track information incl cluster range ref");
@@ -170,7 +171,7 @@ void ResidualsContainer::fillStatisticsBranches()
   }
 }
 
-void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span<const UnbinnedResid> resid, const gsl::span<const o2::tpc::TrackDataCompact> trkRefsIn, const gsl::span<const o2::tpc::TrackData>* trkDataIn, const o2::ctp::LumiInfo* lumiInput)
+void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span<const UnbinnedResid> resid, const gsl::span<const DetInfoResid> detInfoRes, const gsl::span<const o2::tpc::TrackDataCompact> trkRefsIn, const gsl::span<const o2::tpc::TrackData>* trkDataIn, const o2::ctp::LumiInfo* lumiInput)
 {
   // receives large vector of unbinned residuals and fills the sector-wise vectors
   // with binned residuals and statistics
@@ -185,13 +186,14 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp
     firstSeenTF = ti.tfCounter;
   }
   for (const auto& residIn : resid) {
-    ++nUnbinnedResidualsInTF;
     bool counterIncremented = false;
     if (writeUnbinnedResiduals) {
       unbinnedRes.push_back(residIn);
+      detInfoUnbRes.push_back(detInfoRes.size() ? detInfoRes[nUnbinnedResidualsInTF] : DetInfoResid{});
       ++nResidualsTotal;
       counterIncremented = true;
     }
+    ++nUnbinnedResidualsInTF;
     if (!writeBinnedResid) {
       continue;
     }
@@ -247,6 +249,7 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp
     timeMS = orbitReset + ti.tfCounter * o2::constants::lhc::LHCOrbitMUS * 1.e-3;
     treeOutResidualsUnbinned->Fill();
     unbinnedRes.clear();
+    detInfoUnbRes.clear();
     trackInfo.clear();
   }
   tfOrbits.push_back(ti.firstTForbit);
@@ -338,6 +341,9 @@ void ResidualsContainer::merge(ResidualsContainer* prev)
   if (writeUnbinnedResiduals) {
     prev->treeOutResidualsUnbinned->SetBranchAddress("res", &unbinnedResPtr);
     prev->treeOutResidualsUnbinned->SetBranchAddress("trackInfo", &trackInfoPtr);
+    prev->treeOutResidualsUnbinned->SetBranchAddress("CTPLumi", &lumiTF);
+    prev->treeOutResidualsUnbinned->SetBranchAddress("timeMS", &timeMS);
+    prev->treeOutResidualsUnbinned->SetBranchAddress("detInfo", &detInfoUnbResPtr);
     for (int i = 0; i < treeOutResidualsUnbinned->GetEntries(); ++i) {
       treeOutResidualsUnbinned->GetEntry(i);
       prev->treeOutResidualsUnbinned->Fill();

Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h

@@ -29,7 +29,9 @@
 #pragma link C++ class o2::tpc::TrackResiduals::VoxRes + ;
 #pragma link C++ class o2::tpc::TrackResiduals::VoxStats + ;
 #pragma link C++ class o2::tpc::UnbinnedResid + ;
+#pragma link C++ class o2::tpc::DetInfoResid + ;
 #pragma link C++ class std::vector < o2::tpc::UnbinnedResid> + ;
+#pragma link C++ class std::vector < o2::tpc::DetInfoResid> + ;
 #pragma link C++ class std::vector < o2::tpc::TrackResiduals::LocalResid> + ;
 #pragma link C++ class std::vector < o2::tpc::TrackResiduals::VoxStats> + ;
 #pragma link C++ class o2::tpc::ResidualAggregator + ;