ITS3: fix study its3-tracking-study-workflow

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h

@@ -80,6 +80,7 @@ class ITSTrackingInterface
   TimeFrameN* mTimeFrame = nullptr;
 
  protected:
+  virtual void overrideParameters(std::vector<TrackingParameters>& t, std::vector<VertexingParameters>& v) {}
   virtual void requestTopologyDictionary(framework::ProcessingContext& pc);
   virtual void loadROF(gsl::span<const itsmft::ROFRecord>& trackROFspan,
                        gsl::span<const itsmft::CompClusterExt> clusters,

Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx

@@ -50,6 +50,7 @@ void ITSTrackingInterface::initialise()
   }
   auto trackParams = TrackingMode::getTrackingParameters(mMode);
   auto vertParams = TrackingMode::getVertexingParameters(mMode);
+  overrideParameters(trackParams, vertParams);
   LOGP(info, "Initializing tracker in {} phase reconstruction with {} passes for tracking and {}/{} for vertexing", TrackingMode::toString(mMode), trackParams.size(), o2::its::VertexerParamConfig::Instance().nIterations, vertParams.size());
   mTracker->setParameters(trackParams);
   mVertexer->setParameters(vertParams);

Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx

@@ -61,8 +61,11 @@ float Vertexer<NLayers>::clustersToVertices(LogFunc logger)
       mMemoryPool->setMaxMemory(mVertParams[iteration].MaxMemory);
       unsigned int nTracklets01{0}, nTracklets12{0};
       logger(fmt::format("=== ITS {} Seeding vertexer iteration {} summary:", mTraits->getName(), iteration));
-      trkPars.PhiBins = mTraits->getVertexingParameters()[0].PhiBins;
-      trkPars.ZBins = mTraits->getVertexingParameters()[0].ZBins;
+      const auto& currentVtxPars = mTraits->getVertexingParameters()[iteration];
+      trkPars.PhiBins = currentVtxPars.PhiBins;
+      trkPars.ZBins = currentVtxPars.ZBins;
+      trkPars.LayerZ = currentVtxPars.LayerZ;
+      trkPars.LayerRadii = currentVtxPars.LayerRadii;
       trkPars.PassFlags = mVertParams[iteration].PassFlags;
       trkPars.PassFlags.set(IterationStep::FirstPass, IterationStep::RebuildClusterLUT);
       auto timeInitIteration = evaluateTask(&Vertexer::initialiseVertexer, StateNames[mCurStep = Init], iteration, evalLog, trkPars);

Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h

@@ -27,6 +27,7 @@ class ITS3TrackingInterface final : public its::ITSTrackingInterface
   void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) final;
 
  protected:
+  void overrideParameters(std::vector<o2::its::TrackingParameters>& t, std::vector<o2::its::VertexingParameters>& v) final;
   void requestTopologyDictionary(framework::ProcessingContext& pc) final;
   void loadROF(gsl::span<const itsmft::ROFRecord>& trackROFspan,
                gsl::span<const itsmft::CompClusterExt> clusters,

Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx

@@ -72,19 +72,20 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf,
 
   // check for missing/empty/unset rofs
   // the code requires consistent monotonically increasing input without gaps
-  const auto& timing = tf->getROFOverlapTableView().getLayer(layer);
+  const auto& timing = tf->getROFOverlapTableView().getLayer(layer >= 0 ? layer : 0);
   if (timing.mNROFsTF != rofs.size()) {
     LOGP(fatal, "Received inconsistent number of rofs on layer:{} expected:{} received:{}", layer, timing.mNROFsTF, rofs.size());
   }
 
-  its::bounded_vector<uint8_t> clusterSizeVec(clusters.size(), tf->getMemoryPool().get());
+  its::bounded_vector<uint8_t> clusterSizeVec(tf->getMemoryPool().get());
+  clusterSizeVec.reserve(clusters.size());
 
   for (size_t iRof{0}; iRof < rofs.size(); ++iRof) {
     const auto& rof = rofs[iRof];
     for (int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) {
       const auto& c = clusters[clusterId];
       const auto sensorID = c.getSensorID();
-      const auto layer = geom->getLayer(sensorID);
+      const auto lay = geom->getLayer(sensorID);
 
       float sigmaY2{0}, sigmaZ2{0}, sigmaYZ{0};
       uint8_t clusterSize{0};
@@ -110,30 +111,37 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf,
       }
       math_utils::detail::bringToPMPi(alpha); // alpha is defined on -Pi,Pi
 
-      tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), x, alpha,
+      tf->addTrackingFrameInfoToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), x, alpha,
                                       std::array<float, 2>{y, trkXYZ.z()},
                                       std::array<float, 3>{sigmaY2, sigmaYZ, sigmaZ2});
 
       /// Rotate to the global frame
-      tf->addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), tf->getUnsortedClusters()[layer].size());
-      tf->addClusterExternalIndexToLayer(layer, clusterId);
+      tf->addClusterToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), tf->getUnsortedClusters()[lay].size());
+      tf->addClusterExternalIndexToLayer(lay, clusterId);
+    }
+    // effectively calculating an exclusive sum
+    if (layer >= 0) {
+      tf->mROFramesClusters[layer][iRof + 1] = tf->mUnsortedClusters[layer].size();
+    } else {
+      for (unsigned int iL{0}; iL < tf->mUnsortedClusters.size(); ++iL) {
+        tf->mROFramesClusters[iL][iRof + 1] = tf->mUnsortedClusters[iL].size();
+      }
     }
-    tf->mROFramesClusters[layer][iRof + 1] = (int)tf->getUnsortedClusters()[layer].size();
   }
 
   tf->setClusterSize(layer, clusterSizeVec);
 
-  if (layer == 1) {
-    for (auto& v : tf->mNTrackletsPerCluster) {
-      v.resize(tf->getUnsortedClusters()[1].size());
-    }
-    for (auto& v : tf->mNTrackletsPerClusterSum) {
-      v.resize(tf->getUnsortedClusters()[1].size() + 1);
+  if (layer == 1 || layer == -1) {
+    for (auto i = 0; i < tf->mNTrackletsPerCluster.size(); ++i) {
+      tf->mNTrackletsPerCluster[i].resize(tf->mUnsortedClusters[1].size());
+      tf->mNTrackletsPerClusterSum[i].resize(tf->mUnsortedClusters[1].size() + 1);
     }
   }
 
   if (mcLabels != nullptr) {
-    tf->mClusterLabels[layer] = mcLabels;
+    tf->mClusterLabels[layer >= 0 ? layer : 0] = mcLabels;
+  } else {
+    tf->mClusterLabels[layer >= 0 ? layer : 0] = nullptr;
   }
   return 0;
 }

Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx

@@ -9,18 +9,31 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+#include "ITS3Base/SpecsV2.h"
 #include "ITS3Reconstruction/TrackingInterface.h"
 #include "ITS3Reconstruction/IOUtils.h"
 #include "ITSBase/GeometryTGeo.h"
-#include "ITStracking/FastMultEstConfig.h"
-#include "ITStracking/TrackingConfigParam.h"
 #include "DataFormatsITSMFT/DPLAlpideParam.h"
 #include "DetectorsBase/GRPGeomHelper.h"
-#include "Framework/DeviceSpec.h"
 
 namespace o2::its3
 {
 
+void ITS3TrackingInterface::overrideParameters(std::vector<o2::its::TrackingParameters>& t, std::vector<o2::its::VertexingParameters>& v)
+{
+  // only override IB radii
+  for (auto& tt : t) {
+    tt.LayerRadii[0] = constants::radii[0];
+    tt.LayerRadii[1] = constants::radii[1];
+    tt.LayerRadii[2] = constants::radii[2];
+  }
+  for (auto& vv : v) {
+    vv.LayerRadii[0] = constants::radii[0];
+    vv.LayerRadii[1] = constants::radii[1];
+    vv.LayerRadii[2] = constants::radii[2];
+  }
+}
+
 void ITS3TrackingInterface::requestTopologyDictionary(framework::ProcessingContext& pc)
 {
   pc.inputs().get<o2::its3::TopologyDictionary*>("itscldict"); // just to trigger the finaliseCCDB

Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx

@@ -96,6 +96,9 @@ void Digitizer::process(const std::vector<itsmft::Hit>* hits, int evID, int srcI
               return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID();
             });
   for (int i : hitIdx | std::views::filter([&](int idx) {
+                 if (layer < 0) {
+                   return true;
+                 }
                  return mGeometry->getLayer((*hits)[idx].GetDetectorID()) == layer;
                })) {
     processHit((*hits)[i], mROFrameMax, evID, srcID, layer);
@@ -170,7 +173,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast, int layer)
     auto& extra = *(mExtraBuff.front().get());
     for (size_t iChip{0}; iChip < mChips.size(); ++iChip) {
       auto& chip = mChips[iChip];
-      if (chip.isDisabled() || mGeometry->getLayer(chip.getChipIndex()) != layer) {
+      if (chip.isDisabled() || (layer >= 0 && mGeometry->getLayer(chip.getChipIndex()) != layer)) {
         continue;
       }
       chip.addNoise(mROFrameMin, mROFrameMin, &mParams);
@@ -216,7 +219,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast, int layer)
   }
 }
 
-void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID, int layer)
+void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID, int lay)
 {
   // convert single hit to digits
   auto chipID = hit.GetDetectorID();
@@ -248,18 +251,16 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
   }
   float tTot = mParams.getSignalShape().getMaxDuration();
   // frame of the hit signal start wrt event ROFrame
-  int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv(layer));
+  int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv(lay));
   // frame of the hit signal end  wrt event ROFrame: in the triggered mode we read just 1 frame
-  uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv(layer) : roFrameRel;
+  uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv(lay) : roFrameRel;
   int nFrames = roFrameRelMax + 1 - roFrameRel;
   uint32_t roFrameMax = mNewROFrame + roFrameRelMax;
-  if (roFrameMax > maxFr) {
-    maxFr = roFrameMax; // if signal extends beyond current maxFrame, increase the latter
-  }
+  maxFr = std::max(roFrameMax, maxFr); // if signal extends beyond current maxFrame, increase the latter
 
   // here we start stepping in the depth of the sensor to generate charge diffision
-  int detID{hit.GetDetectorID()};
-  const auto& matrix = mGeometry->getMatrixL2G(detID);
+  const int layer = mGeometry->getLayer(chipID);
+  const auto& matrix = mGeometry->getMatrixL2G(chipID);
   int nSteps = chip.isIB() ? mParams.getIBNSimSteps() : mParams.getNSimSteps();
   float nStepsInv = chip.isIB() ? mParams.getIBNSimStepsInv() : mParams.getNSimStepsInv();
   math_utils::Vector3D<float> xyzLocS, xyzLocE;
@@ -422,7 +423,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
         continue;
       }
       uint16_t colIS = icol + colS;
-      registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl, layer);
+      registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl, lay);
     }
   }
 }

Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h

@@ -36,7 +36,7 @@ struct ITS3TrackingStudyParam : o2::conf::ConfigurableParamHelper<ITS3TrackingSt
   int minTPCCls{110};
 
   // propagator
-  o2::base::PropagatorImpl<float>::MatCorrType CorrType = o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrLUT;
+  o2::base::PropagatorImpl<float>::MatCorrType CorrType = o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrTGeo;
 
   /// studies
   bool doDCA = false;

Detectors/Upgrades/ITS3/study/macros/PlotPulls.C

@@ -9,11 +9,14 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-/// \file PlotDCA.C
-/// \brief Simple macro to plot ITS3 impact parameter resolution
+/// \file PlotPulls.C
+/// \brief Simple macro to plot ITS3 pulls
 
 #if !defined(__CLING__) || defined(__ROOTCLING__)
+#include <array>
+#include <cmath>
 #include <memory>
+#include <vector>
 
 #include <TROOT.h>
 #include <TCanvas.h>
@@ -27,7 +30,6 @@
 
 #include "ReconstructionDataFormats/GlobalTrackID.h"
 #include "ReconstructionDataFormats/Track.h"
-#include "ReconstructionDataFormats/DCA.h"
 #include "SimulationDataFormat/MCTrack.h"
 #endif
 
@@ -49,7 +51,13 @@ void PlotPulls(const char* fName = "its3TrackStudy.root")
 {
   TH1::SetDefaultSumw2();
   std::unique_ptr<TFile> inFile(TFile::Open(fName));
+  if (!inFile || inFile->IsZombie()) {
+    return;
+  }
   auto tree = inFile->Get<TTree>("pull");
+  if (!tree) {
+    return;
+  }
 
   uint8_t src; // track type
   tree->SetBranchAddress("src", &src);
@@ -133,13 +141,13 @@ void PlotPulls(const char* fName = "its3TrackStudy.root")
   std::vector<TH1D*> projs;
   const char* fitOpt{"QWMERSB"};
   for (int i{0}; i < o2::track::kNParams; ++i) {
-    for (auto iPt{0}; iPt < nPtBins; ++iPt) {
+    for (auto iPt{0}; iPt < nPtBins - 1; ++iPt) {
       auto hProj = pulls[i]->ProjectionY(Form("%s_%d", pulls[i]->GetName(), iPt), iPt + 1, iPt + 1);
       hProj->SetName(Form("p%s_pt%d", pNames[i], iPt));
       hProj->SetTitle(Form("Pull %s #it{p}_{T}#in[%.2f, %.2f)", pNames[i], ptLimits[iPt], ptLimits[iPt + 1]));
       projs.push_back(hProj);
       if (hProj->GetEntries() < 100) {
-        return;
+        continue;
       }
       fGaus->SetParameter(1, 0);
       fGaus->SetParameter(2, 1);
@@ -153,14 +161,18 @@ void PlotPulls(const char* fName = "its3TrackStudy.root")
     }
   }
 
-  hMahDist2->Scale(1. / hMahDist2->Integral("width"));
   TF1* fchi2Fit = new TF1("fchi2_fit", chi2_pdf, 0.1, 6, 3);
   fchi2Fit->SetParNames("A", "k", "s");
   fchi2Fit->SetParameter(0, 1);
   fchi2Fit->SetParameter(1, 5);
   fchi2Fit->SetParameter(2, 1);
-  auto fitres = hMahDist2->Fit(fchi2Fit, "RMQS");
-  fitres->Print();
+  if (hMahDist2->Integral("width") > 0.) {
+    hMahDist2->Scale(1. / hMahDist2->Integral("width"));
+    auto fitres = hMahDist2->Fit(fchi2Fit, "RMQS");
+    if (fitres.Get()) {
+      fitres->Print();
+    }
+  }
 
   TFile outFile("plotPulls.root", "RECREATE");
   for (int i{0}; i < o2::track::kNParams; ++i) {